How to make a simulation of triaxial-compress at 2 dimensions ?

Dear All !

I have already success to make a 2-dimensions box with Gen , and got a good result when I make a unaxial-compression simulation. But when I tried to make a simulation of traxial-compression , I met some problem as the stress waved obviously with the increase of strain . In the other hand , I also wondered whether the force or the displacement in Z-direction the system ingored when the fuction force2dComputations() was called .

Here is my script of triaxila-compress simulation :

---------------------------------------------------------------------------------

#rot_compress.py: A uniaxial compression simulation using ESyS-Particle
# Author: D. Weatherley
# Date: 27 December 2008
# Organisation: ESSCC, University of Queensland
# (C) All rights reserved, 2008.
#
#
#import the appropriate ESyS-Particle modules:
from esys.lsm import *
from esys.lsm.util import *
from esys.lsm.geometry import *
from WallLoader import WallLoaderRunnable
from ServoWallLoader import ServoWallLoaderRunnable
from esys.lsm.util import Vec3, BoundingBox

#instantiate a simulation object:
sim = LsmMpi (numWorkerProcesses = 1, mpiDimList = [1,1,1])

modelname="model_2d_100.geo"
xdim=100
ydim=100
zdim=0
minR = 0.2
maxR = 1.0
timestep = 2000000
ramptime = 500
velocity = 0.0002
kn = 4.288e3
c_press = 4
area = xdim*zdim

#initialise the neighbour search algorithm:
sim.initNeighbourSearch (
particleType = "RotSphere",
gridSpacing = 2.5*maxR,
verletDist = 0.2*minR
)
#set the number of timesteps and timestep increment:
sim.setNumTimeSteps (timestep)
sim.setTimeStepSize (0.0005)
sim.force2dComputations (True)
#specify the spatial domain for the simulation
domain = BoundingBox(Vec3(-20,-20,-20), Vec3(200,200,200))
sim.setSpatialDomain (domain)
#create a prism of spherical particles:
sim.readGeometry(modelname)

#create rotational elastic-brittle bonds between clusters:
pp_bonds = sim.createInteractionGroup (
BrittleBeamPrms(
name="pp_bonds",
youngsModulus=kn,
poissonsRatio=0.25,
cohesion=6,
tanAngle=4.0,
tag=0
)
)

#initialise frictional interactions for unbonded particles:
sim.createInteractionGroup (
RotElasticPrms(
name="contaction",
normalK=kn,
scaling = True
)
)

sim.createInteractionGroup (
FrictionPrms(
name="friction",
youngsModulus=kn,
poissonsRatio=0.25,
dynamicMu=0.4,
staticMu=0.6
)
)

#create an exclusion between bonded and frictional interactions:
sim.createExclusion (
interactionName1 = "pp_bonds",
interactionName2 = "friction"
)

sim.createExclusion (
interactionName1 = "pp_bonds",
interactionName2 = "contact"
)

############################## wall #########################################

#create a wall at the bottom of the model:
sim.createWall (
name = "bottom_wall",
posn = Vec3(0.0000, 0.0000, 0.0000),
normal = Vec3(0.0000, 1.0000, 0.0000)
)

#create a wall at the top of the model:
sim.createWall (
name = "top_wall",
posn = Vec3(0.0000, ydim, 0.0000),
normal = Vec3(0.0000, -1.0000, 0.0000)
)

#create 4 wall at the round of the model:
sim.createWall (
   name = "round_wall1",
   posn = Vec3(0.0000, 0.0000, 0.0000),
   normal = Vec3(1.0000, 0.0000, 0.0000)
)

sim.createWall (
   name = "round_wall3",
   posn = Vec3(xdim, 0.00, zdim),
   normal = Vec3(-1.0000, 0.0000, 0.0000)
)

#specify elastic repulsion from the bottom wall:
sim.createInteractionGroup (
NRotElasticWallPrms (
name = "bottom_wall_repel",
wallName = "bottom_wall",
normalK = 10*kn
)
)

#specify elastic repulsion from the top wall:
sim.createInteractionGroup (
NRotElasticWallPrms (
name = "top_wall_repel",
wallName = "top_wall",
normalK = 10*kn
)
)

sim.createInteractionGroup (
NRotElasticWallPrms (
name = "round_wall1_repel",
wallName = "round_wall1",
normalK = 10*kn
)
)

sim.createInteractionGroup (
NRotElasticWallPrms (
name = "round_wall3_repel",
wallName = "round_wall3",
normalK = 10*kn
)
)

#add translational viscous damping:
sim.createInteractionGroup (
LinDampingPrms(
name="damping1",
viscosity=0.02,
maxIterations=50
)
)

#add rotational viscous damping:
sim.createInteractionGroup (
RotDampingPrms(
name="damping2",
viscosity=0.02,
maxIterations=50
)
)

############################ load #################################

##add a wall loader to move the top wall:
wall_loader1 = WallLoaderRunnable(
LsmMpi = sim,
wallName = "top_wall",
vPlate = Vec3 (0.0, -1*velocity, 0.0),
startTime = 500,
rampTime = ramptime
)
sim.addPreTimeStepRunnable (wall_loader1)

#add a wall loader to move the bottom wall:
wall_loader2 = WallLoaderRunnable(
LsmMpi = sim,
wallName = "bottom_wall",
vPlate = Vec3 (0.0, velocity, 0.0),
startTime = 500,
rampTime = ramptime
)
sim.addPreTimeStepRunnable (wall_loader2)

servo_loader1 = ServoWallLoaderRunnable(
   LsmMpi = sim,
   interactionName = "round_wall1_repel",
   force = Vec3 (c_press*ydim, 0.0, 0.0),# 5MPa #
   startTime = 0,
   rampTime = 500
)
sim.addPreTimeStepRunnable (servo_loader1)

servo_loader3 = ServoWallLoaderRunnable(
   LsmMpi = sim,
   interactionName = "round_wall3_repel",
   force = Vec3 (-1*c_press*ydim, 0.0, 0.0),# 5MPa #
   startTime = 0,
   rampTime = 500
)
sim.addPreTimeStepRunnable (servo_loader3)

########################### record ######################################

#create a FieldSaver to store number of bonds:
sim.createFieldSaver (
InteractionScalarFieldSaverPrms(
interactionName="pp_bonds",
fieldName="count",
fileName="c1-nbonds.dat",
fileFormat="SUM",
beginTimeStep=0,
endTimeStep=timestep,
timeStepIncr=100
)
)
#create a FieldSaver to store the total kinetic energy of the particles:
sim.createFieldSaver (
ParticleScalarFieldSaverPrms(
fieldName="e_kin",
fileName="c1-ekin.dat",
fileFormat="SUM",
beginTimeStep=0,
endTimeStep=timestep,
timeStepIncr=100
)
)
#create a FieldSaver to store potential energy stored in bonds:
sim.createFieldSaver (
InteractionScalarFieldSaverPrms(
interactionName="pp_bonds",
fieldName="potential_energy",
fileName="c1-epot.dat",
fileFormat="SUM",
beginTimeStep=0,
endTimeStep=timestep,
timeStepIncr=100
)
)

#create a FieldSaver to wall positions:
posn_saver = WallVectorFieldSaverPrms(
wallName=["bottom_wall", "top_wall"],
fieldName="Position",
fileName="c1-out_Position.dat",
fileFormat="RAW_SERIES",
beginTimeStep=0,
endTimeStep=timestep,
timeStepIncr=100
)
sim.createFieldSaver(posn_saver)

posn_saver2 = WallVectorFieldSaverPrms(
wallName=["round_wall1", "round_wall3"],
fieldName="Position",
fileName="c1-out_Position2.dat",
fileFormat="RAW_SERIES",
beginTimeStep=0,
endTimeStep=timestep,
timeStepIncr=100
)
sim.createFieldSaver(posn_saver2)

#create a FieldSaver to wall forces:
force_saver = WallVectorFieldSaverPrms(
wallName=["bottom_wall", "top_wall"],
fieldName="Force",
fileName="c1-out_Force.dat",
fileFormat="RAW_SERIES",
beginTimeStep=0,
endTimeStep=timestep,
timeStepIncr=100
)
sim.createFieldSaver(force_saver)

#create a FieldSaver to round_wall forces:
force_saver2 = WallVectorFieldSaverPrms(
   wallName=["round_wall1", "round_wall3"],
# wallName=["round_wall1", "round_wall3"],
   fieldName="Force",
   fileName="c1-out_Force2.dat",
   fileFormat="RAW_SERIES",
   beginTimeStep=0,
   endTimeStep=timestep,
   timeStepIncr=100
)
sim.createFieldSaver(force_saver2)

sim.createCheckPointer (
CheckPointPrms (
fileNamePrefix = "c1-snapshot",
beginTimeStep = 0,
endTimeStep = timestep,
timeStepIncr = 1000
)
)

#execute the simulation:
sim.run()

-------------------------------------------------------------------------

With kind regards,

Zhongwen