soil structure interaction

Hi, I’m now doing the research about ‘soil structure interaction’, but I met some difficulties that I can’t overcome.
My model is 2D.
1. I generate four boundaries:

id_Mat=O.materials.append(FrictMat(young=100e9,poisson=0.3,density=1000,frictionAngle = radians(10)))
Mat=O.materials[id_Mat]
id_box1 = O.bodies.append(box((0.5,12.5,0.5),(0.5,12.5,0.5),fixed=True,material=Mat))
id_box2 = O.bodies.append(box((101.5,12.5,0.5),(0.5,12.5,0.5),fixed=True,material=Mat))
id_box3 = O.bodies.append(box((51,20.5,0.5),(50,0.5,0.5),fixed=True,material=Mat))
id_box4 = O.bodies.append(box((51,-0.5,0.5),(70,0.5,0.5),fixed=True,material=Mat))

2. Generate particles

sp = pack.SpherePack()
sp.makeCloud(Vector3(1,0,0.5),Vector3(101,20,0.5),num=N_particles,rMean=AvgRadius,rRelFuzz=.33,distributeMass=False)
sp.toSimulation(material = SphereMat)

3. Radius expansion method. ConfPress1=50kPa. I don’t know if it’s proper to use TriaxialStressController here.

triax1=TriaxialStressController(
    maxMultiplier=1.+1.5e5/Young, # spheres growing factor (fast growth)
    finalMaxMultiplier=1.+4e3/Young,
    internalCompaction = True, # If true the confining pressure is generated by growing particles
    stressMask = 7,
    computeStressStrainInterval = 10,
    goal1 = ConfPress1,
goal2 = ConfPress1,

The first phase of work is described above. It’s the ‘generation.py’.

Next, I will write the second file called ‘conpression.py’. The outcome (generation.yade.bz2) of first file ‘generation.py’ will be loaded and I want to apply 100kPa normal stress on the upper boundary and cancel it’s ‘fixed’, but I don’t know how to do it (1. How to apply 100kPa on the upper boundary? 2. How to modify the code of previous file to cancel the ‘fixed’ ?)

When the above two stages are finished, the third file (shear.py) will work. I want the bottom boundary to move to the left and I don’t know how to realize it.

The code of ‘generation.py’:
from __future__ import print_function
from builtins import range
from yade import pack,export,qt,plot,os,timing
import matplotlib; matplotlib.rc('axes',grid=True)
import pylab

#Material constants
Density = 3000
FrictionAngle = 1.5
PoissonRatio = 0.5
Young = 300e6
Damp = 0.5
AvgRadius = 0.1
N_particles = 2000

#Confining variables
ConfPress1 = -50000

#define material for all bodies:
id_Mat=O.materials.append(FrictMat(young=100e9,poisson=0.3,density=1000,frictionAngle = radians(10)))
Mat=O.materials[id_Mat]

id_box1 = O.bodies.append(box((0.5,12.5,0.5),(0.5,12.5,0.5),fixed=True,material=Mat))
id_box2 = O.bodies.append(box((101.5,12.5,0.5),(0.5,12.5,0.5),fixed=True,material=Mat))
id_box3 = O.bodies.append(box((51,20.5,0.5),(50,0.5,0.5),fixed=True,material=Mat))
id_box4 = O.bodies.append(box((51,-0.5,0.5),(70,0.5,0.5),fixed=True,material=Mat))

SphereMat = O.materials.append(FrictMat(young = Young, poisson = PoissonRatio, frictionAngle = radians(FrictionAngle), density = Density))
sp = pack.SpherePack()
sp.makeCloud(Vector3(1,0,0.5),Vector3(101,20,0.5),num=N_particles,rMean=AvgRadius,rRelFuzz=.33,distributeMass=False)
sp.toSimulation(material = SphereMat)

triax1=TriaxialStressController(
    maxMultiplier=1.+1.5e5/Young, # spheres growing factor (fast growth)
    finalMaxMultiplier=1.+4e3/Young,
    internalCompaction = True, # If true the confining pressure is generated by growing particles
    stressMask = 7,
    computeStressStrainInterval = 10,

    goal1 = ConfPress1,
    goal2 = ConfPress1,
)

newton=NewtonIntegrator(damping=Damp)

###engine
O.engines=[
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
    InteractionLoop(
     [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
     [Ip2_FrictMat_FrictMat_FrictPhys()],
     [Law2_ScGeom_FrictPhys_CundallStrack()]
    ),
    GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8, defaultDt=4*utils.PWaveTimeStep()),
    triax1,
    newton,
    PyRunner(realPeriod=10,command='checkUnbalanced()',label='check'),
    PyRunner(command='addPlotData()',iterPeriod=2000,label='record'),
    TriaxialStateRecorder(iterPeriod=2000,file='WallStresses')
]

# Simulation stop conditions defination
def checkUnbalanced():
    unb=unbalancedForce()
    mStress = (triax1.stress(triax1.wall_right_id)[0]+triax1.stress(triax1.wall_top_id)[1])/2.0
    s1 = triax1.stress(triax1.wall_right_id)[0]
    s2 = triax1.stress(triax1.wall_top_id)[1]
    if unb<0.01 and abs(ConfPress1-mStress)/abs(ConfPress1)<0.01:
       O.pause()
       O.save('generation.yade.bz2')

# collect history of data
def addPlotData():
    unb = unbalancedForce()
    mStress = -(triax1.stress(triax1.wall_right_id)[0]+triax1.stress(triax1.wall_top_id)[1])/2.
    area = (O.bodies[1].state.pos[0]-O.bodies[0].state.pos[0])*(O.bodies[3].state.pos[1]-O.bodies[2].state.pos[1])
    Porosity = 1-sum(pi*b.shape.radius**2 for b in O.bodies if isinstance(b.shape,Sphere))/2000

    plot.addData(e11=-triax1.strain[0], e22=-triax1.strain[1], e33=-triax1.strain[2],
     ev=-triax1.strain[0]-triax1.strain[1]-triax1.strain[2],
   s11=-triax1.stress(triax1.wall_right_id)[0],
   s22=-triax1.stress(triax1.wall_top_id)[1],
   s33=-triax1.stress(triax1.wall_front_id)[2],
   ub=unbalancedForce(),
   dstress=(-triax1.stress(triax1.wall_top_id)[1])-(-triax1.stress(triax1.wall_right_id)[0]),
   disx=triax1.width,
   disy=triax1.height,
   disz=triax1.depth,
   i=O.iter,
   porosity=Porosity,
   ke=utils.kineticEnergy(),
   totalEnergy=O.energy.total()
   )

    plot.saveDataTxt('generation.txt.bz2')

Is there any problem in the code?

Above all, my questions are:
1. How to apply 100kPa on the upper boundary?
2. How to modify the code of previous file to cancel the ‘fixed’ ?
3. How to move the bottom boundary to the left?
4. Is anyone willing to share the similar examples’ code?
5. Is there some examples of direct shear (not simple shear)?

Thank you!