A problem about the Stress BC

[zhifengmuyun]

Hi everyone,
I want to use Neumann BC to simulate a simple uniaxial tensile process and here are the codes:
from underworld import UWGeodynamics as GEO
import numpy as np

u = GEO.UnitRegistry

inX = -200.0e3 * u.meter
maxX = 200.0e3 * u.meter
minY = - 80.0e3 * u.meter
maxY = 40.0e3 * u.meter
minZ = -200.0e3 * u.meter
maxZ = 200.0e3 * u.meter

density = 3300 * u.kilogram / u.metre3
gravity = 9.81 * u.meter / u.second2

eta = 1e21 * u.pascal * u.second
mu = 9.5e7 * u.pascal
dt_e = 1e5 * u.year

alpha = eta / mu # Maxwell relaxation time

eta_eff = ( eta * dt_e ) / (alpha + dt_e)

ref_length = 1e5 * u.meter
ref_eta = eta
ref_rho = density
ref_force = density * gravity
ref_time = dt_e / 10

KL = ref_length
Kt = ref_time
KM = ref_force * KL2 * Kt2

GEO.scaling_coefficients["[length]"] = KL
GEO.scaling_coefficients["[time]"] = Kt
GEO.scaling_coefficients["[mass]"]= KM

Lambda = 100.0 * u.kilometer # displacement wavelength
t_relax = 4 * np.pi * eta / (Lambda * density * gravity)
tMax = (t_relax * 5).to_base_units()

Model = GEO.Model(elementRes=(32, 16),
minCoord=(minX, minY),
maxCoord=(maxX, maxY),
periodic=[False, False],
gravity=[0.,-gravity])

Model.outputDir="ViscoElastic_StressBC"

GEO.rcParams["default.outputs"] = [
"velocityField",
"pressureField",
"materialField",
"viscosityField",
"strainRateField",
"projStressTensor",
"projDensityField"

]

air = Model.add_material(name="Air")
material1 = Model.add_material(name="Material1")
material2 = Model.add_material(name="Material2")

import underworld.function as fn

coords = fn.input()

x_total = maxX - minX
x_third = x_total / 3.0

conditions = [

(coords[1] > 0, air.index),

((coords[1] < 0), material1.index),




(True, air.index)

]

Model.materialField.data[:] = fn.branching.conditional(conditions).evaluate(Model.swarm)

air.density = 0.1 * u.kilogram / u.meter**3
air.viscosity = 1e18 * u.pascal * u.second

material1.density = 3500u.kilogram / u.metre**3 # 例如：地壳密度
material1.viscosity = 5e20 u.pascal * u.second

material1.elasticity = GEO.Elasticity(shear_modulus=mu, observation_time=dt_e)

material1.plasticity = GEO.VonMises(cohesion=1e7 * (u.kilogram * u.meter**-1 * u.second**-2))

Model.minViscosity = 1e18 * u.pascal * u.second
Model.maxViscosity = 1e25 * u.pascal * u.second

Model.set_velocityBCs(bottom=[0,0],top=[None,None],left=[None,0],right=[None,0])
Model.set_stressBCs(bottom=[None,None],top=[None,None],left=[15u.pascal,0],right = [-15u.megapascal,0])

import numpy as np
coords = np.ndarray((1, 2))
coords[0, 1] = GEO.nd(minY + 0.15 * (maxY - minY))
coords
Model.add_passive_tracers(name="Single", vertices=coords)

Model.solver.set_inner_method("lu")

tTracer = [GEO.nd(Model.time)]
previousStress_xy = [0.]
totalStress_xy = [0.]

def postSolveHook():
global tTracer, previousStress_xy, totalStress_xy
tTracer.append(GEO.nd(Model.time))
previousStress_xy.append(
Model._previousStressField[2].evaluate(Model.swarm)[0])
totalStress_xy.append(
Model._stressFn[2].evaluate(Model.swarm)[0][0])

Model.post_solve_functions["Stress Record"] = postSolveHook

checkpoint_interval = tMax/5000

Model.run_for(duration=tMax/500,
dt=1e-3*t_relax,
checkpoint_interval=checkpoint_interval)
I think the result will show that The object is compressed, becoming shorter and thicker. However, the figure is this:

I checked the swarms and found the shape resembles the density field. Maybe the stress BC works something wrong. Could you please give me some suggestions?