2D biaxial compression task completion

Hello,

(1) you could solve it yourself :-)

for b in O.bodies:
> if isinstance(b.shape,Sphere):
> pos = b.state.pos
> if pos[0] <3.5 and pos[1] < 3.5: b.shape.color = (1,0,0) #
> area 1
> elif pos[0] >= 3.5 and pos[1] <3.5: b.shape.color = (0,1,0) #
> area 2
> elif pos[0] >= 3.5 and pos[1] >= 3.5: b.shape.color = (0,0,1) #
> area 3
> else: b.shape.pos = (1,1,0) #
> area 4
>

b.shape.pos = (1,1,0) does not assign color ;-) sorry if I wrote it in one
of previous emails..

(2) see https://answers.launchpad.net/yade/+question/219884

cheers
Jan

Dear Jan Stránský:
     Thanks for helping me last time and I can getting a correct color assignation for particles in the cell, it also bring me the right mothed of attributes-determination for particles. And now I have other deep questions for asking your help as follows:
(1). I wanna employ a strict pattern for particles-generation rather than make a randomly distribution and a simple example can be taken here for describing what the problem is like that:
There are 3 types of particles in my specimen and the radii of which are a,b and c. Each type of particles has a equal proportion, which is,of course, 33.3%. The total number of particles in my specimen is d. So, I want to obtain a perfect particles generation with the parameters a,b,c and d.
how can I do for this purpose?
(2). In the process of initial state determination, I wanna use the porosity as simulation stop condition and once the porosity of my specimen is equal to which I specified before the simulation starting, the simulation can be stopped immediately. how can I write this orders in my script.
(3). I use the classes "ThreeDTriaxialEngine" for applying loads in the second step of biaxial compression(named triax02) and I copy this from the examples of YADE without understanding how it works. I am not sure whether it is correct for me to apply such loads as that:
1. In the horizontal direction, I want to apply constant forces 2N/m;
2. In the Vertical direction, I want to apply constant velocities 0.05m/s.

SEEKING for your help for correction this script:

### fundamental details of application ###
# unicode: UTF-8
Filename='2D-simulation'
from yade import pack,os
################################# preprocessing for simulation ##########################################
### prescribing variables and functions for simulation controller ###
# material defination
spheremat = O.materials.append(ViscElMat(kn=4e5,ks=4e5,cn=.0,cs=.0,density=1500,frictionAngle=25.565))
wallmat = O.materials.append(ViscElMat(kn=4e5,ks=4e5,cn=.0,cs=.0,density=2600,frictionAngle=25.565))
# walls defination
mn,mx=Vector3(0,0,0),Vector3(0.07,0.07,0.1)
wallIds=O.bodies.append(utils.aabbWalls([mn,mx],thickness=.0001,material=wallmat))
# ThreeDTriaxialEngine defination for initial-state determination(the first calculation step)
triax01=ThreeDTriaxialEngine(
 wall_bottom_id=wallIds[2],wall_top_id=wallIds[3],
 wall_left_id=wallIds[0],wall_right_id=wallIds[1],
 wall_back_id=wallIds[4],wall_front_id=wallIds[5],
 wall_front_activated = False,wall_back_activated = False,
 internalCompaction=False,
 stressControl_1 = True, stressControl_2 = True,stressControl_3 = True,
 computeStressStrainInterval =10,
 sigma_iso = 1.25e4,
 sigma1 = 1.25e4,
 sigma2 = 1.25e4,
 sigma3 = 1.25e4,
 strainRate1 = 0.03,strainRate2 = 0.03,
)
# ThreeDTriaxialEngine defination for triaxial compression(the second calculation step)
triax02=ThreeDTriaxialEngine(
 wall_bottom_id=wallIds[2],wall_top_id=wallIds[3],
 wall_left_id=wallIds[0],wall_right_id=wallIds[1],
 wall_back_id=wallIds[4],wall_front_id=wallIds[5],
 wall_front_activated = False,wall_back_activated = False,
 internalCompaction=False,
 stressControl_1 = True, stressControl_2 = False,stressControl_3 = True,
 computeStressStrainInterval =10,
 sigma_iso = 3,
 sigma1 = 3,
# sigma2 = 3,
# sigma3 = 1.25e4,
 strainRate1 = -0.008,strainRate2 = 0.05,
)
# Simulation stop controller defination
'''#def checkUnbalanced():
    unb=unbalancedForce()
    meanS=(triax01.stress(triax01.wall_right_id)[0]+triax01.stress(triax01.wall_top_id)[1])/2
    q=unb
    r=abs(meanS-triax01.sigma_iso)/triax01.sigma_iso
    if q<0.01 and r<1e-4:
       O.pause()
       O.save('first-step.xml'.format(Filename))'''
################################# control flow for simulation ##########################################
# particles generation
O.periodic=1
O.cell.setBox(0.07,0.07,0.07)
sp=pack.SpherePack()
sp.makeCloud((0,0,.05),(0.07,0.07,.05),rMean=0.0008,rRelFuzz=0.3,num=840,periodic=True)
sp.toSimulation(material=spheremat)
# determining colors for particles in different aeras of the cell
for b in O.bodies:
    if isinstance(b.shape,Sphere):
         pos = b.state.pos
         if pos[0] <0.035 and pos[1] < 0.035: b.shape.color = (1,0,0) # area 1
         elif pos[0] >= 0.035 and pos[1] <0.035: b.shape.color = (0,1,0) # area 2
         elif pos[0] >= 0.035 and pos[1] >= 0.035: b.shape.color = (0,0,1) # area 3
         elif pos[0] < 0.035 and pos[1] >= 0.035: b.shape.color = (0,10,1) # area 4
         else: b.shape.pos = (1,1,0)
O.periodic=0
# blockedDOFs
for b in O.bodies:
 if isinstance(b.shape,Sphere):
   b.state.blockedDOFs='zXY'
# Simulation assembly for the first step
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Wall_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
  [Ip2_ViscElMat_ViscElMat_ViscElPhys()],
  [Law2_ScGeom_ViscElPhys_Basic()]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax01,
 NewtonIntegrator(damping=.1),
# PyRunner(command='checkUnbalanced()',iterPeriod=200)
]
# first step of simulation startting with a correct inheriting for the next step
O.dt = 2e-4
O.run(27000,True);
O.save('first-step.xml')
O.wait()
f = file("/home/fzg/fu/result-first.dat",'w')
f.write("varibles='X-refcordinate','Y-refcordinate','Z-refcordinate'\n\n")
f.write('%-16s %-16s %-16s\n'%('X-refcordinate','Y-refcordinate','Z-refcordinate'))
for b in O.bodies:
    if isinstance(b.shape,Sphere):
       refpos= b.state.refPos
       poscu = b.state.pos
       displ = b.state.displ()
   # pos0 = poscu[0]-0.035
    # pos1 = poscu[1]-0.035
     # pos2 = poscu[2]-0.035
       rad = b.shape.radius
       strainrate = b.state.vel
 # vel0 = strainrate[0]*pos0
  # vel1 = strainrate[1]*pos1
       f.write('%-16g %-16g %-16g\n'%(poscu[0],poscu[1],poscu[2]))
f.write('################################ ends ##########################################')
f.close()
# loading inheriting
O.load('first-step.xml')
# Simulation assembly for the second step
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Wall_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
  [Ip2_ViscElMat_ViscElMat_ViscElPhys()],
  [Law2_ScGeom_ViscElPhys_Basic()]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax02,
 NewtonIntegrator(damping=.1),
]
# second step of simulation startting
O.dt = 5e-4
O.run(10000,True);
# whole task over
O.save('final-step.xml'.format(Filename));
O.wait()
################################## postprocessing for simulation ######################################################
################################## postprocessing for simulation ######################################################
f = file("/home/fzg/fu/result.dat",'w')
f.write('## This is the result data of 2D simulation\n\n')
f.write('## There are 12 types of varibles in this data as follows:\n\n')
f.write("varibles='X-refcordinate','Y-refcordinate','Z-refcordinate','X-cordinate','Y-cordinate','Z-cordinate','Radius','StrainRate-xx','StrainRate-yy','Velocity-xx','Velocity-yy','Ids'\n\n")
f.write('%-16s %-16s %-16s %-16s %-16s %-16s %-16s %-16s %-16s %-16s %-16s %-16s\n'%('X-refcordinate','Y-refcordinate','Z-refcordinate','X-cordinate','Y-cordinate','Z-cordinate','Radius','Velocity-xx','Velocity-yy','resi-Velocity-xx','resi-Velocity-yy','Ids'))
for b in O.bodies:
    if isinstance(b.shape,Sphere):
       refpos= b.state.refPos
       poscu = b.state.pos
       displ = b.state.displ()
   # pos0 = poscu[0]-0.035
    # pos1 = poscu[1]-0.035
     # pos2 = poscu[2]-0.035
       rad = b.shape.radius
       strainrate = b.state.vel
       resivel0 = strainrate[0]*(1.035-poscu[0])
       resivel1 = strainrate[1]*(1.035-poscu[1])
       f.write('%-16g %-16g %-16g %-16g %-16g %-16g %-16g %-16g %-16g %-16g %-16g %-16g\n'%(refpos[0],refpos[1],refpos[2],poscu[0],poscu[1],poscu[2],rad,strainrate[0],strainrate[1],resivel0,resivel0,b.id))
f.write('################################ ends ##########################################')
f.close()
def rename():
    global Filename
    os.rename("/home/fzg/fu/result.dat","/home/fzg/fu/{0}.plt".format(Filename))
rename()

Hello,

(1). I wanna employ a strict pattern for particles-generation rather than
> make a randomly distribution and a simple example can be taken here for
> describing what the problem is like that:
> There are 3 types of particles in my specimen and the radii of which are
> a,b and c. Each type of particles has a equal proportion, which is,of
> course, 33.3%. The total number of particles in my specimen is d. So, I
> want to obtain a perfect particles generation with the parameters a,b,c and
> d.
> how can I do for this purpose?
>

see [1] and parameters psdSizes and psdCumm. I've never used it myself, but
I think that you can define some psd close to only three discrete particle
sizes. Probably the number of particles would not be exactly what you
defined, but it should be close. Alternatively, you can write your own
myMakeCloud function, it should not be a big deal.

> (2). In the process of initial state determination, I wanna use the
> porosity as simulation stop condition and once the porosity of my specimen
> is equal to which I specified before the simulation starting, the
> simulation can be stopped immediately. how can I write this orders in my
> script.
>

def checkPorosity():
  porosity = someCommandToDeterminePorosity() # maybe something like 1 -
sum(4/3.*pi*b.shape.radius**3 for b in O.bodies)/O.cell.volume
  if porosity > definedPorosity: O.pause()
O.engines = [
  ...
  PyRunner(iterPeriod=10,command="checkPorosity()")
]

> (3). I use the classes "ThreeDTriaxialEngine" for applying loads in the
> second step of biaxial compression(named triax02) and I copy this from the
> examples of YADE without understanding how it works. I am not sure whether
> it is correct for me to apply such loads as that:
> 1. In the horizontal direction, I want to apply constant forces 2N/m;
> 2. In the Vertical direction, I want to apply constant velocities 0.05m/s.
>

personally I have never used ThreeDTriaxialEngine, so perhaps somebody else
can help you

cheers
Jan

[1]
https://yade-dem.org/doc/yade.pack.html#yade._packSpheres.SpherePack.makeCloud

(1) You could do this:

Yade [8]: sp2.makeCloud((0,0,0),(1,1,1),rMean=0.1, rRelFuzz=0,num=10,periodic=1)
Yade [9]: sp2.makeCloud((0,0,0),(1,1,1),rMean=0.07, rRelFuzz=0,num=10,periodic=1)
Yade [10]: sp2.makeCloud((0,0,0),(1,1,1),rMean=0.05, rRelFuzz=0,num=10,periodic=1)
Yade [11]: sp.toSimulation()

(2) -> is shown in triax tutorial script no 1 (#1, L148)

(3) Triaxial engines keeps stress constant, not force. If you want constant force, you will have to update the target stress periodicaly in the script, using the target force and the updated surface:
triax.sigma1 = force/triax.height #or something like that

>I copy this from the examples of YADE without understanding how it works

Maybe try understanding the tutorial scripts in details first, check also the documentation of ThreeDTriax, then back to your own script. It would probably be easier that way, and it would also help to narrow down what is not clear.
Note that ThreeDTriax is deprecated. TriaxialStressController is relatively similar and easier to use.

Bruno

#1. https://github.com/yade/trunk/blob/master/examples/triax-tutorial/script-session1.py

Dear Jan Stránský:
     Thanks for helping me last time and I have some new problems to ask your help, which are about getting output of the deviator stress in biaxial compression simulation. I wanna get the reaults of deviator stress using these script :

### fundamental details of application ###
# unicode: UTF-8
Filename='2D-simulation'
from yade import pack,os
################################# preprocessing for simulation ##########################################
### prescribing variables and functions for simulation controller ###
# material defination
spheremat = O.materials.append(ViscElMat(kn=4e5,ks=4e5,cn=.0,cs=.0,density=1500,frictionAngle=5.8))
wallmat = O.materials.append(ViscElMat(kn=4e5,ks=4e5,cn=.0,cs=.0,density=1500,frictionAngle=5.8))
# walls defination
mn,mx=Vector3(0,0,0),Vector3(0.07,0.07,0.1)
wallIds=O.bodies.append(utils.aabbWalls([mn,mx],thickness=.0001,material=wallmat))
# ThreeDTriaxialEngine defination for initial-state determination(the first calculation step)
triax01=ThreeDTriaxialEngine(
 wall_bottom_id=wallIds[2],wall_top_id=wallIds[3],
 wall_left_id=wallIds[0],wall_right_id=wallIds[1],
 wall_back_id=wallIds[4],wall_front_id=wallIds[5],
 wall_front_activated = False,wall_back_activated = False,
 internalCompaction=False,
 stressControl_1 = True, stressControl_2 = True,stressControl_3 = True,
 computeStressStrainInterval =10,
 sigma_iso = 1.25e4,
 sigma1 = 1.25e4,
 sigma2 = 1.25e4,
 sigma3 = 1.25e4,
 strainRate1 = 0.03,strainRate2 = 0.03,
)
# ThreeDTriaxialEngine defination for triaxial compression(the second calculation step)
triax02=ThreeDTriaxialEngine(
 wall_bottom_id=wallIds[2],wall_top_id=wallIds[3],
 wall_left_id=wallIds[0],wall_right_id=wallIds[1],
 wall_back_id=wallIds[4],wall_front_id=wallIds[5],
 wall_front_activated = False,wall_back_activated = False,
 internalCompaction=False,
 stressControl_1 = True, stressControl_2 = False,stressControl_3 = True,
 computeStressStrainInterval =10,
 sigma_iso = 1.25e4,
 sigma1 = 1.25e4,
# sigma2 = 3,
# sigma3 = 1.25e4,
 strainRate1 = -0.072,strainRate2 = 0.214,
)
# Simulation stop controller defination
'''#def checkUnbalanced():
    unb=unbalancedForce()
    meanS=(triax01.stress(triax01.wall_right_id)[0]+triax01.stress(triax01.wall_top_id)[1])/2
    q=unb
    r=abs(meanS-triax01.sigma_iso)/triax01.sigma_iso
    if q<0.01 and r<1e-4:
       O.pause()
       O.save('first-step.xml'.format(Filename))'''
################################# control flow for simulation ##########################################
# particles generation
O.periodic=1
O.cell.setBox(0.07,0.07,0.07)
sp=pack.SpherePack()
sp.makeCloud((0,0,.05),(0.07,0.07,.05),rMean=0.0008,rRelFuzz=0.3,num=840,periodic=True)
sp.toSimulation(material=spheremat)
# determining colors for particles in different aeras of the cell
for b in O.bodies:
    if isinstance(b.shape,Sphere):
         pos = b.state.pos
         if pos[0] <0.035 and pos[1] < 0.035: b.shape.color = (1,0,0) # area 1
         elif pos[0] >= 0.035 and pos[1] <0.035: b.shape.color = (0,1,0) # area 2
         elif pos[0] >= 0.035 and pos[1] >= 0.035: b.shape.color = (0,0,1) # area 3
         elif pos[0] < 0.035 and pos[1] >= 0.035: b.shape.color = (0,10,1) # area 4
         else: b.shape.pos = (1,1,0)
O.periodic=0
# blockedDOFs
for b in O.bodies:
 if isinstance(b.shape,Sphere):
   b.state.blockedDOFs='zXY'
# Simulation assembly for the first step
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Wall_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
  [Ip2_ViscElMat_ViscElMat_ViscElPhys()],
  [Law2_ScGeom_ViscElPhys_Basic()]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax01,
 NewtonIntegrator(damping=.1),
# PyRunner(command='checkUnbalanced()',iterPeriod=200)
]
# first step of simulation startting with a correct inheriting for the next step
O.dt = 2e-4
O.run(21300,True);
for b in O.bodies:
      b.state.vel=(0,0,0)
O.save('first-step.xml')
O.wait()
f = file("/home/fzg/fu/result-first.dat",'w')
f.write("varibles='X-refcordinate','Y-refcordinate','Z-refcordinate'\n\n")
f.write('%-16s %-16s %-16s\n'%('X-refcordinate','Y-refcordinate','Z-refcordinate'))
for b in O.bodies:
    if isinstance(b.shape,Sphere):
       refpos= b.state.refPos
       poscu = b.state.pos
       displ = b.state.displ()
   # pos0 = poscu[0]-0.035
    # pos1 = poscu[1]-0.035
     # pos2 = poscu[2]-0.035
       rad = b.shape.radius
       strainrate = b.state.vel
 # vel0 = strainrate[0]*pos0
  # vel1 = strainrate[1]*pos1
       f.write('%-16g %-16g %-16g\n'%(poscu[0],poscu[1],poscu[2]))
f.write('################################ ends ##########################################')
f.close()
f=file("/home/fzg/fu/macro.dat",'a')
f.write('%-16s %-16s %-16s %-16s %-16s %-16s %-16s\n'%('iters','time','area','axial-strain','volumetric-strain','porosity','DeviatorS')),
f.close()
# loading inheriting
O.load('first-step.xml')
# Simulation assembly for the second step
def macroresults():
 wall_left = O.bodies[0].state.pos[0]
 wall_right = O.bodies[1].state.pos[0]
 wall_bottom = O.bodies[2].state.pos[1]
 wall_top = O.bodies[3].state.pos[1]
 wall_back = O.bodies[4].state.pos[2]
 wall_front = O.bodies[5].state.pos[2]
 x_dimension = wall_right - wall_left
 y_dimension = wall_top - wall_bottom
 #z_dimension = wall_front - wall_back
 area = x_dimension * y_dimension
        volumetricstrain=(area-0.0021919)/0.0021919
 #area = (x_dimension-thick) * (y_dimension-thick) ## 'thick' is the thickness of the walls
 #vol = x_dimension * y_dimension * z_dimension
 porosity = 1-sum(pi*b.shape.radius**2 for b in O.bodies if isinstance(b.shape,Sphere))/area
        t=(O.iter-21300)*1e-4
        axialstrain =-(y_dimension-0.0468177)/0.0468177
        DeviatorS= (triax02.stress(triax02.wall_top_id)[1]-triax02.stress(triax02.wall_right_id)[0])
        f=file("/home/fzg/fu/macro.dat",'a')
        f.write('%-16g %-16g %-16g %-16g %-16g %-16g %-16g\n'%(O.iter,t,area,axialstrain,volumetricstrain,porosity,DeviatorS)),
        f.close()
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb(),Bo1_Wall_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
  [Ip2_ViscElMat_ViscElMat_ViscElPhys()],
  [Law2_ScGeom_ViscElPhys_Basic()]
 ),
 GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
 triax02,
 NewtonIntegrator(damping=.1),
        PyRunner(iterPeriod=500,command='macroresults()')
]
# second step of simulation startting
O.dt = 0.5e-4
macroresults()
O.run(16000,True);
# whole task over
O.save('final-step.xml'.format(Filename));
O.wait()
macroresults()
##################################################
f = file("/home/fzg/fu/result.dat",'w')
f.write('## This is the result data of 2D simulation\n\n')
f.write('## There are 12 types of varibles in this data as follows:\n\n')
f.write("varibles='X-refcordinate','Y-refcordinate','Z-refcordinate','X-cordinate','Y-cordinate','Z-cordinate','Radius','StrainRate-xx','StrainRate-yy','Velocity-xx','Velocity-yy','Ids'\n\n")
f.write('%-16s %-16s %-16s %-16s %-16s %-16s %-16s %-16s %-16s %-16s %-16s %-16s\n'%('X-refcordinate','Y-refcordinate','Z-refcordinate','X-cordinate','Y-cordinate','Z-cordinate','Radius','Velocity-xx','Velocity-yy','resi-Velocity-xx','resi-Velocity-yy','Ids'))
for b in O.bodies:
    if isinstance(b.shape,Sphere):
       refpos= b.state.refPos
       poscu = b.state.pos
       displ = b.state.displ()
   # pos0 = poscu[0]-0.035
    # pos1 = poscu[1]-0.035
     # pos2 = poscu[2]-0.035
       rad = b.shape.radius
       strainrate = b.state.vel
       resivel0 = strainrate[0]*(1.035-poscu[0])
       resivel1 = strainrate[1]*(1.035-poscu[1])
       f.write('%-16g %-16g %-16g %-16g %-16g %-16g %-16g %-16g %-16g %-16g %-16g %-16g\n'%(refpos[0],refpos[1],refpos[2],poscu[0],poscu[1],poscu[2],rad,strainrate[0],strainrate[1],resivel0,resivel0,b.id))
f.write('################################ ends ##########################################')
f.close()
def rename():
    global Filename
    os.rename("/home/fzg/fu/result.dat","/home/fzg/fu/{0}.plt".format(Filename))
rename()

But the output text bring me a surprise that almost all the numbers of deviator stress are 0, I can not determine where is wrong and seek your help!

Thanks Jan Stránský, that solved my question.