yade-users team mailing list archive

[Leonard <question698487@xxxxxxxxxxxxxxxxxxxxx>]

New question #698487 on Yade:
https://answers.launchpad.net/yade/+question/698487

Hi,
I'd like to ask is it possible that change model boundary from rigid wall to periodic boundary? My attempt is as follow:

I use script1 (from Bruno [1]) to generate a sample with rigid wall boundary, note that I set periodic=True in makeCloud. And then I save the sample as sample.yade.gz.

Then I import the sample.yade.gz into script2 (from periodic-triax example[2]), and I remove the rigid walls and give it a periodic boundary.

However, I failed in switching the boundary. The sample exploded immediately after one O.step() forward. I am not sure whether it is impossible to switch the boundary or I have not found out the correct way. Thereby, I'd like to ask for you advice for that: 
-Is it possible that change model boundary from rigid wall to periodic boundary after I generate a sample under the condition of rigid wall?
-The reason why I want to switch the boundary is that using the example code[1] (which is rigid wall) is easy to make a sample with target confining presure and porosity by using radii expansion methods (by setting internalCompaction=True), but it seems that there is no radii expansion method for periodic engine. Which makes it uneasy for using only periodic engine to generate target sample. Do you have any idea about how to integrate [1] and [2], specifically, making sample by using [1] and carrying out triaxial loading by using periodic engine as [2]?

Below are the MWE:
################script1 begin##########Around 15 second for running ######
from yade import pack, plot
nRead=readParamsFromTable(
 num_spheres=1000,
 compFricDegree = 30,
 key='_triax_base_',
 unknownOk=True
)
from yade.params import table

num_spheres=table.num_spheres
key=table.key
targetPorosity = 0.43
compFricDegree = table.compFricDegree
finalFricDegree = 30
rate=-0.02
damp=0.2
stabilityThreshold=0.01
young=5e6
mn,mx=Vector3(0,0,0),Vector3(0.7,0.7,1.4)

O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

sp=pack.SpherePack()
sp.makeCloud(mn,mx,-1,0.3333,num_spheres,periodic=True,porosity=0.95,seed=1)
O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])

Gl1_Sphere.quality=3

### DEFINING ENGINES ###
triax=TriaxialStressController(
 maxMultiplier=1.+2e4/young,
 finalMaxMultiplier=1.+2e3/young,
 thickness = 0,
 stressMask = 7,
 internalCompaction=True,## please change True to False here to look at what I say
)

newton=NewtonIntegrator(damping=damp)

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),
 triax,
 newton
]

Gl1_Sphere.stripes=0
if nRead==0: yade.qt.Controller(), yade.qt.View()

### APPLYING CONFINING PRESSURE ###
triax.goal1=triax.goal2=triax.goal3=-10000
while 1:
 O.run(1000, True)
 unb=unbalancedForce()
 print 'unbalanced force:',unb,' mean stress: ',triax.meanStress
 if unb<stabilityThreshold and abs(-10000-triax.meanStress)/10000<0.001:
  break

print "### Isotropic state saved ###"
### REACHING A SPECIFIED POROSITY PRECISELY ###
import sys
while triax.porosity>targetPorosity:
 # we decrease friction value and apply it to all the bodies and contacts
 compFricDegree = 0.95*compFricDegree
 setContactFriction(radians(compFricDegree))
 print "\r Friction: ",compFricDegree," porosity:",triax.porosity,
 sys.stdout.flush()
 O.run(500,1)
print "### Compacted state saved ###"
O.save('sample.yade.gz')
################script1 done##########

################script2 begin########
from __future__ import print_function
sigmaIso = -10e3
from yade import pack, qt, plot

O.load('sample.yade.gz')

for i in range(6):
    O.bodies.erase(i)

O.periodic = True
O.cell.setBox(0.7,0.7,1.4)
sp = pack.SpherePack()
sp.fromSimulation()
sp.toSimulation()

O.engines = [
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb()]),
        InteractionLoop([Ig2_Sphere_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_FrictPhys()], [Law2_ScGeom_FrictPhys_CundallStrack()]),
        PeriTriaxController(
                label='triax',
                # specify target values and whether they are strains or stresses
                goal=(sigmaIso, sigmaIso, sigmaIso),
                stressMask=7,
                # type of servo-control
                dynCell=True,
                maxStrainRate=(10, 10, 10),
                # wait until the unbalanced force goes below this value
                maxUnbalanced=.1,
                relStressTol=1e-3,
                # call this function when goal is reached and the packing is stable
                doneHook='compactionFinished()'
        ),
        NewtonIntegrator(damping=.2),
        # PyRunner(command='addPlotData()', iterPeriod=100),
]
O.dt = .5 * PWaveTimeStep()
O.trackEnergy = True
def compactionFinished():
	print('stage1 finished')
	# set the current cell configuration to be the reference one
	O.cell.trsf = Matrix3.Identity
	# change control type: keep constant confinement in x,z, 20% compression in y
	triax.goal = (sigmaIso,-0.4, sigmaIso)
	triax.stressMask = 5
	# allow faster deformation along x,z to better maintain stresses
	triax.maxStrainRate = (1., 0.5, 1.)
	# next time, call triaxFinished instead of compactionFinished
	triax.doneHook = 'triaxFinished()'
	# do not wait for stabilization before calling triaxFinished
	triax.maxUnbalanced = 10

def triaxFinished():
	print('Finished')
	O.pause()
################script2 done##########

Thanks!
Leonard

[1]https://gitlab.com/yade-dev/trunk/-/blob/mpi/examples/triax-tutorial/script-session1.py
[2]https://gitlab.com/yade-dev/trunk/blob/master/doc/sphinx/tutorial/06-periodic-triaxial-test.py


-- 
You received this question notification because your team yade-users is
an answer contact for Yade.