Shape anisotropy under isotropic confinement
Hello. I am trying to create samples with shape anisotropy under isotropic confinement. Currently, my code is encountering some issues (explained in detail in the "current problem" section below). I am wondering if there is a problem with my approach. Thank you! Below is the code that reproduces the issue:
# -*- encoding=utf-8 -*-
#states:
#1、Making sample according to the specified radius and quantity within a specific Z-axis area,currently top of the box.
#2、Replace particles with clump immediately, then change and lock clump's orientation if needed
#3、Apply gravity to consolidate
#4、Unlock its orientation if needed
#5、Repeat steps 1-4 until the desired quantity is reached
#6、Delete particles above Z=10, undo gravity, and delete clumps above Z=10 again
#7、Delete old walls and create new walls around current particles
#8、Use TriaxialStressController to apply less pressure to complete sample preparation
#current problem:
#1、only when num_spheres is high enough (for example, num_spheres=2000,rMean=0.45) that the simulation can be run, otherwise the particles will fly out of the model in stage 8
#2、When unlock_after_drop=True, it will need more particles (for example, num_spheres=3000,rMean=0.4), but no matter what the previous value is, the test will appear particles with a radius of nan in stage 8.
import time,sys
from yade import pack,export
import numpy as np
readParamsFromTable(
key='default'
)
from yade.params.table import *
#Physical parameters
num_spheres=3000
rMean=0.4
young=8e8
poisson=0.3
density=2810
radvar=0.33
msgoal=-10e3#Stress during sample preparation
msdegree=0#Friction angle during sample preparation
#Control parameters
makingsample_divid_time=10#Number of times of sample preparation
ifclump=True#Whether to introduce Clump, if not, the following two parameters will be invalid.
ifchangeori=True#Whether the clump is directional. Otherwise, the following parameter will be invalid.
iflockori=True#Whether need to lock the direction after change its orientation
unlock_after_drop=True#Whether unlock the direction after the clump settle down
mn,mx=Vector3(0,0,0),Vector3(10,10,30)
stabilityThreshold=0.01
targetPorosity=0.6
old_color=(0.5,0.5,0.5)
new_clumps_color = (0.5, 0.7, 0.7)
new_spheres_color = (0.2, 0.7, 0.2)
print(time.strftime("%Y-%m-%d %H:%M:%S", time.localtime()))
datadir='./data/makingsampledata-{}/'.format(key)
vtkdir='./data/makingsamplevtk-{}/'.format(key)
O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=radians(msdegree),density=density,label='spheres'))
sp=SpherePack()
O.materials.append(FrictMat(young=young,poisson=poisson,frictionAngle=0,density=0,label='frictionless'))
walls=aabbWalls(extrema=[mn,mx],thickness=1e-10,material='frictionless',color=old_color)
wallIds=O.bodies.append(walls)
vtkrecorder=VTKRecorder(fileName=vtkdir+'3d-vtk-', recorders=['all'], iterPeriod=1000,label='vtkrecorder',ascii=True)
runningchecker=PyRunner(iterPeriod=1000,command='runningcheck(expected_time)',label='runningchecker')
import os
def makingdir(dir):
if os.path.exists(dir):
#Delete all files under the directory
print('Delete all files under the directory:{}'.format(dir))
for root, dirs, files in os.walk(dir, topdown=False):
for name in files:
os.remove(os.path.join(root, name))
for name in dirs:
os.rmdir(os.path.join(root, name))
os.makedirs(dir,exist_ok=True)
makingdir(datadir)
makingdir(vtkdir)
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_MindlinPhys(
# krot=mskr,
# ktwist=ktw,
# eta=eta,
)
],
[Law2_ScGeom_MindlinPhys_Mindlin(
#includeMoment=True,
)]
),
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
#triax,
PyRunner(iterPeriod=1000,command='history()',label='triaxrecorder'),
runningchecker,
vtkrecorder,
NewtonIntegrator(damping=.3,gravity=Vector3(0,0,-9.8*10),label='NewtonG')
]
O.trackEnergy = True
from yade import plot
## a function saving variables
def history():
plot.addData(
e11=-triax.strain[0]*100, e22=-triax.strain[1]*100, e33=-triax.strain[2]*100,
ev=100*(triax.strain[0]+triax.strain[1]+triax.strain[2]),
s11=-triax.stress(0)[0]/1000,
s22=-triax.stress(2)[1]/1000,
s33=-triax.stress(4)[2]/1000,
diffstress=(-triax.stress(2)[1]-(-triax.stress(0)[0]-triax.stress(4)[2])/2)/1000,
Etot=O.energy.total(),
i=O.iter,
**O.energy
)
plot.plots={
'i':('e11','e22','e33',),
'i ':('s11','s22','s33'),
' i':(O.energy.keys,None,'Etot'),
'e22':('s11','s22','s33',None,'ev'),
}
plot.plot()
def runningcheck(expected_time):
for b in O.bodies:
if b==None or not isinstance(b.shape,Sphere):
continue
else:
if isnan(b.shape.radius):
raise Exception('The radius of sphere {} is nan!'.format(b.id))
break
mn_box=Vector3(O.bodies[wallIds[0]].state.pos[0],O.bodies[wallIds[2]].state.pos[1],O.bodies[wallIds[4]].state.pos[2])
mx_box=Vector3(O.bodies[wallIds[1]].state.pos[0],O.bodies[wallIds[3]].state.pos[1],O.bodies[wallIds[5]].state.pos[2])
for i in O.bodies:
if not isinstance(i.shape,Sphere):
continue
location=i.state.pos
for j in range(3):
if location[j]>mx_box[j] or location[j]<mn_box[j]:
print('The current position of the ball {}: {}'.format(i.id,i.state.pos))
raise Exception('ball {} ran outside the model!'.format(i.id))
relRadList3 = [
1,1,1,1,
1,1,1,1,
1,1,1,1,
1,1,1,1,]
relPosList3 = [
[-1.5, 1.5,0],[-0.5, 1.5,0],[ 0.5, 1.5,0],[ 1.5, 1.5,0],
[-1.5, 0.5,0],[-0.5, 0.5,0],[ 0.5, 0.5,0],[ 1.5, 0.5,0],
[-1.5,-0.5,0],[-0.5,-0.5,0],[ 0.5,-0.5,0],[ 1.5,-0.5,0],
[-1.5,-1.5,0],[-0.5,-1.5,0],[ 0.5,-1.5,0],[ 1.5,-1.5,0],
]
templates = []
templates.append(clumpTemplate(relRadii=relRadList3, relPositions=relPosList3))
def mainlogic():
global wallIds
#state1:
triaxrecorder.dead=True
runningchecker.dead=True
makingsample_count=1
mn_2=Vector3(mn[0],mn[1],(mx[2]-mn[2])*(1-1/makingsample_divid_time))
num_spheres_2=int(num_spheres/makingsample_divid_time)
num_actual=sp.makeCloud(
rMean=rMean,
minCorner=mn_2,
maxCorner=mx,
rRelFuzz=radvar,
num=num_spheres_2,
)
if num_actual!=num_spheres_2:
print('Warning! The number of particles may be insufficient!')
while makingsample_count<=makingsample_divid_time:
print('Sample prereration time(s):{}'.format(makingsample_count))
sp.toSimulation(color=old_color,material='spheres')
if ifclump:
#state2:
print('Perform particle replacement')
clumpreplaceresult=O.bodies.replaceByClumps(templates,[1],discretization=20)
if ifchangeori:
for i in clumpreplaceresult:
current_vector=O.bodies[i[0]].state.ori*Vector3(0,0,1)
projected_vector=Vector3(0,0,1)
qua=Quaternion.Identity
Quaternion.setFromTwoVectors(qua,current_vector,projected_vector)
O.bodies[i[0]].state.ori=qua*O.bodies[i[0]].state.ori
if iflockori:
for i in clumpreplaceresult:
if O.bodies[i[0]]!=None:#It may become empty, and you need to bypass the empty one.
O.bodies[i[0]].state.blockedDOFs='XY'
for j in i[1]:
O.bodies[j].state.blockedDOFs='XY'
#Delete the clumps that were in contact with the wall after being replaced. These clumps were moved outside the wall after previous oriented.
count=0
mn_box=Vector3(O.bodies[wallIds[0]].state.pos[0],O.bodies[wallIds[2]].state.pos[1],O.bodies[wallIds[4]].state.pos[2])
mx_box=Vector3(O.bodies[wallIds[1]].state.pos[0],O.bodies[wallIds[3]].state.pos[1],O.bodies[wallIds[5]].state.pos[2])
for i in clumpreplaceresult:
iferaseclump=False
for member_id in O.bodies[i[0]].shape.members.keys():
for j in range(3):
if (O.bodies[member_id].state.pos[j]>mx_box[j]) or (O.bodies[member_id].state.pos[j]<mn_box[j]):#位置在墙外
count+=1
iferaseclump=True
break
if iferaseclump:
break
if iferaseclump:
O.bodies.erase(i[0],eraseClumpMembers=1)
print('Remove {} clumps that were replaced outside the wall'.format(count))
global expected_time
expected_time=1e5
#After running for 100 steps, delete the clusters with members outside the model.
#These clumps are the product of excessive overlap in the previous step of orientation or are knocked out by clump and should be deleted.
O.run(100,True)
calm()
count=0
mn_box=Vector3(O.bodies[wallIds[0]].state.pos[0],O.bodies[wallIds[2]].state.pos[1],O.bodies[wallIds[4]].state.pos[2])
mx_box=Vector3(O.bodies[wallIds[1]].state.pos[0],O.bodies[wallIds[3]].state.pos[1],O.bodies[wallIds[5]].state.pos[2])
for i in O.bodies:
if isinstance(i.shape,Clump):
iferaseclump=False
for member_id in i.shape.members.keys():
for j in range(3):
if (O.bodies[member_id].state.pos[j]>mx_box[j]) or (O.bodies[member_id].state.pos[j]<mn_box[j]):#位置在墙外
count+=1
iferaseclump=True
break
if iferaseclump:
break
if iferaseclump:
O.bodies.erase(i.id,eraseClumpMembers=1)
print('Remove {} clumps that were moved outside the wall'.format(count))
#state3:
#Let the clumps drop down a little to make room for the next cycle
O.run(3000,True)
makingsample_count+=1
#state4:
#Unlock the orientation of the clumps
#Note: This step will cause Error if uncommented
if unlock_after_drop:
for i in clumpreplaceresult:
if O.bodies[i[0]]!=None:
O.bodies[i[0]].state.blockedDOFs=''
for j in i[1]:
O.bodies[j].state.blockedDOFs=''
while 1:
O.run(1000, True)
unb=unbalancedForce()
print('Waiting for stablize. Current unb:{:.3},calculate speed:{:d},iter:{:d}'.format(unb,int(O.speed),O.iter))
calm()
if unb<0.3:#等待颗粒沉降完成
O.save('{}0.1MakingSphere.yade.gz'.format(datadir))
print("### The Clumps have landed ###")
break
#stage6
count=0
for i in O.bodies:
if isinstance(i.shape,Clump):
iferaseclump=False
for member_id in i.shape.members.keys():
if (O.bodies[member_id].state.pos[2]+O.bodies[member_id].shape.radius>10):#位置+半径比Z轴方向10高
count+=1#Count the number of clumps that z-location is greater than 10
iferaseclump=True
break
for j in range(3):#0,1,2
if (O.bodies[member_id].state.pos[j]>mx_box[j]) or (O.bodies[member_id].state.pos[j]<mn_box[j]):#Location is outside the wall
iferaseclump=True
break
if iferaseclump:
break
if iferaseclump:
O.bodies.erase(i.id,eraseClumpMembers=1)
if count<=10:
raise Exception('Too few particles({}), sample preparation failed!'.format(count))
print('Delete {} particles with z greater than 10'.format(count))
calm()
NewtonG.gravity=Vector3(0,0,0)#撤销重力
while 1:
O.run(1000,True)
unb=unbalancedForce()
calm()
print('release gravaty, current unb:{:.3}, calcuate speed:{:d},iter:{:d}'.format(unb,int(O.speed),O.iter))
if unb<0.05 or isnan(unb):
print('End of stress relief phase')
break
#At this time, there will be new particles higher than z=10, which need to be deleted again.
count=0
for i in O.bodies:
if isinstance(i.shape,Clump):
count+=1
iferaseclump=False
for member_id in i.shape.members.keys():
if O.bodies[member_id].state.pos[2]+O.bodies[member_id].shape.radius>10:
iferaseclump=True
break
if iferaseclump:
O.bodies.erase(i.id,eraseClumpMembers=1)
#stage7
#Move the top wall to the bottom.
#Because I cannot directly modify the position of the top wall, deleting top first and then creating a new one will also result in an error.
#Therefore, remove all walls first and then add new ones back.
for wallid in wallIds:
O.bodies.erase(wallid)
print('The original wall has been deleted')
wall2=aabbWalls(thickness=1e-10,material='frictionless',color=old_color)
wallIds=O.bodies.append(wall2)
print('The new wall has been added')
global triax
triax=TriaxialStressController(
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],
internalCompaction=True,
stressMask=7,
goal1=msgoal,
goal2=msgoal,
goal3=msgoal,
maxMultiplier=1. + 1e-4, # spheres growing factor (fast growth)
finalMaxMultiplier=1. + 1e-5, # spheres growing factor (slow growth)
max_vel=10,
label="triaxcompresser",
)
O.engines=O.engines[:4]+[triax]+O.engines[4:]
triaxrecorder.dead=False
print('Real number of particles:{}'.format(len(O.bodies)))
#stage8 compress the sample
expected_time=3e5
runningchecker.dead=False
while 1:
O.run(1000, True)
unb=unbalancedForce()
print('Compressing! Current unb:{:.3},Mean stress:{:.3},porosity:{:.4},calculate speed:{:d},iter:{:d}'.format(unb,triax.meanStress,triax.porosity,int(O.speed),int(O.iter)))
calm()
if unb<stabilityThreshold and abs(msgoal-triax.meanStress)/(-msgoal)<0.001:
O.save('{}0.2MakingClump.yade.gz'.format(datadir))
print("### Compressing is done ###")
print('Current porosity:{:.3}',triax.porosity)
break
runningcheck(expected_time)
runningchecker.dead=False
if targetPorosity-triax.porosity>0.01:
print('The porosity is much smaller! Maybe the sample production failed!')
for i in O.bodies:
if isinstance(i.shape,Sphere) or isinstance(i.shape,Clump):
i.state.blockedDOFs='XYZ'
print('Successfully Block the DOFs of all particles.')
calm()
O.save('{}0.3finalsample.yade.gz'.format(datadir))
print('Real number of particles:{}'.format(calculate_real_body_count()))
plot.reset()
def calculate_real_body_count():
count=0
for b in O.bodies:
if b==None:
continue
if b.isStandalone or isinstance(b.shape,Clump):
count+=1
return count
mainlogic()