yade-users team mailing list archive

[Fu zuoguang <question226352@xxxxxxxxxxxxxxxxxxxxx>]

Question #226352 on Yade changed:
https://answers.launchpad.net/yade/+question/226352

    Status: Answered => Open

Fu zuoguang is still having a problem:
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()

-- 
You received this question notification because you are a member of
yade-users, which is an answer contact for Yade.