Question #700611 on Yade changed:
https://answers.launchpad.net/yade/+question/700611
Status: Answered => Open
Luis Barbosa is still having a problem:
For this I have to provide txt files for the creation of the packs.
Sorry for the external links:
====================================================================================
https://drive.google.com/file/d/16y5Ia4kUUm6oap8QlnSvQf9uCuMzJ7TN/view?usp=sharing
https://drive.google.com/file/d/1_rUoO-6b2j6e3LSRsMlh8C4vKRu8txoN/view?usp=sharing
https://drive.google.com/file/d/1oR9bHoru-
Us6wvHwn7UpkZwt9IX6GRvE/view?usp=sharing
https://drive.google.com/file/d/1cLyeXTjNuavGRxuFv27HyPxvtn-2EW5i/view?usp=sharing
=====================================================================================
from yade import pack
from yade import bodiesHandling
from yade import export
from yade import utils
from yade import ymport
import math
############################################
### DEFINING VARIABLES AND MATERIALS ###
############################################
# The following 5 lines will be used later for batch execution
nRead=readParamsFromTable(
num_spheres=3000,# number of spheres
compFricDegree = 1, # contact friction during the confining phase
key='_triax_base_', # put you simulation's name here
unknownOk=True
)
from yade.params import table
num_spheres=table.num_spheres# number of spheres
key=table.key
targetPorosity = 0.50 #the porosity we want for the packing
compFricDegree = table.compFricDegree # initial contact friction during the
confining phase (will be decreased during the REFD compaction process)
finalFricDegree = 30 # contact friction during the deviatoric loading
rate=0 # loading rate (strain rate)
damp=0.8 # damping coefficient
stabilityThreshold=0.01 # we test unbalancedForce against this value in
different loops (see below)
young=80e5 # contact stiffness200e4
young2=80e5
youngcoat=50e5
bondstr=1000#2e7
bondstr2=1000
bondstrcoat=10
## create materials for spheres and plates
mat=O.materials.append(JCFpmMat(type=1,young=young,poisson=0.3,frictionAngle=radians(compFricDegree),density=2000,tensileStrength=bondstr,cohesion=bondstr,jointNormalStiffness=0,jointShearStiffness=0,jointCohesion=bondstr,jointFrictionAngle=radians(0),jointDilationAngle=0.0,label='spheres'))
O.materials.append(JCFpmMat(type=1,young=20e7,poisson=0.3,frictionAngle=radians(0),density=2600,tensileStrength=0,cohesion=0,jointNormalStiffness=0,jointShearStiffness=0,jointCohesion=0,jointFrictionAngle=radians(0),jointDilationAngle=0.0,label='walls'))
O.materials.append(JCFpmMat(type=1,young=youngcoat,poisson=0.3,frictionAngle=radians(1),density=1500,tensileStrength=bondstrcoat,cohesion=bondstrcoat,jointNormalStiffness=0,jointShearStiffness=0,jointCohesion=bondstrcoat,jointFrictionAngle=radians(0),jointDilationAngle=0.0,label='spherescoat'))
## create walls around the packing
mn,mx=Vector3(0,0,0),Vector3(0.0015,0.0015,0.0015)
mnbox,mxbox=Vector3(0,0,0),Vector3(0.002,0.00195,0.002)
walls=aabbWalls([mnbox,mxbox],thickness=0,material='walls')
wallIds=O.bodies.append(walls)
O.bodies.append(ymport.textExt("matrix_vtest6.txt", format='x_y_z_r',
shift=Vector3(0,0.00025,0), scale=1.0,material='spheres',color=(0,1,1)))
O.bodies.append(ymport.textExt("coat_vtest5e5.txt", format='x_y_z_r',
shift=Vector3(0,0,0), scale=1.0,material='spherescoat',color=(0,1,1)))
################Particle substitution by large
aggregate######################################################################
bodid=[]
for b in O.bodies:
if b and isinstance(b.shape,Sphere):
# print (b.shape.radius)
if b.state.pos[1]>0.00175:
bodid.append(b.id)
i=0
for p in bodid:
O.bodies.erase(bodid[i])
i=i+1
bodid=[]
a=[]
for b in O.bodies:# in sp:
if b and isinstance(b.shape,Sphere):
# print (b.shape.radius)
if b.shape.radius==0.0002:
bodid.append(b.id)
a.append(b.state.pos)
i=0
for p in bodid:
t=a[i]
f1=O.bodies.append(ymport.textExt("agg2e4_10e6.txt", format='x_y_z_r',
shift=t-Vector3(0,0,0.0002), scale=1.0,material='spheres',color=(0,1,1)))
O.bodies.erase(bodid[i])
i=i+1
bodiddd=[]
aaa=[]
for bbb in O.bodies:# in sp:
if bbb and isinstance(bbb.shape,Sphere):
# print (b.shape.radius)
if bbb.shape.radius==0.0005:
bodiddd.append(bbb.id)
aaa.append(bbb.state.pos)
iii=0
for ppp in bodiddd:
ttt=aaa[iii]
f3=O.bodies.append(ymport.textExt("agg5e5_18e6.txt", format='x_y_z_r',
shift=ttt-Vector3(0,0,0.0005), scale=1.0,material='spheres',color=(0,1,1)))
O.bodies.erase(bodiddd[iii])
iii=iii+1
##############################################################################################################################
#====================================================================================================================================================================
############################
### DEFINING ENGINES ###
############################
triax=TriaxialStressController(
## TriaxialStressController will be used to control stress and strain. It
controls particles size and plates positions.
## this control of boundary conditions was used for instance in
http://dx.doi.org/10.1016/j.ijengsci.2008.07.002
maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)
finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)
thickness = 0,
stressMask = 7,
internalCompaction=False, # If true the confining pressure is generated by
growing particles
wall_front_activated=True,
wall_back_activated=True,
wall_top_activated=True,
wall_bottom_activated=True,
wall_left_activated=True,
wall_right_activated=True,
goal1=-100,
goal2=-100,
goal3=-100,
)
newton=NewtonIntegrator(damping=damp)
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
[Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=-1,label='Physicspheres')],#,xSectionWeibullShapeParameter=1.5,
xSectionWeibullScaleParameter=1
[Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(smoothJoint=False)]
),
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.5),
triax,
newton
]
while 1:
O.run(500,True)
unb=unbalancedForce()
triax.goal1=triax.goal2=triax.goal3=triax.meanStress*1.1
print ('unbalanced force:',unb,' mean stress:
',triax.meanStress,triax.porosity)
if triax.porosity<targetPorosity:
break
print ("### Compacted state saved ###")
print(triax.stress(3)[1])
setContactFriction(radians(finalFricDegree))
#=======================================
triax.stressMask = 2
triax.wall_bottom_activated=0
#now goal2 is the target strain rate
triax.goal1=rate
triax.goal3=rate
triax.goal2=triax.stress(3)[1]
from yade import plot
O.run(10,True)
#strain is logarithmic strain or true strain which is ls=(ln1+e) where e=dl/L
(strain)
ei0=-triax.strain[0];ei1=-triax.strain[1];ei2=-triax.strain[2]
si0=-triax.stress(0)[0];si1=-triax.stress(2)[1];si2=-triax.stress(4)[2]
## a function saving variables
def history():
plot.addData(e11=-triax.strain[0]-ei0, e22=-triax.strain[1]-ei1,
e33=-triax.strain[2]-ei2,
ev=-triax.strain[0]-triax.strain[1]-triax.strain[2],
s11=-triax.stress(triax.wall_right_id)[0]-si0,
s22=-triax.stress(triax.wall_top_id)[1]-si1,
s33=-triax.stress(triax.wall_front_id)[2]-si2,
e=math.exp(-triax.strain[1]-ei1)-1,
pc=-unsat.bndCondValue[2],
sw=unsat.getSaturation(isSideBoundaryIncluded=False),
z1=O.bodies[3].state.pos[1],
i=O.iter)
if 1:
## include a periodic engine calling that function in the simulation loop
O.engines=O.engines[0:5]+[PyRunner(iterPeriod=20,command='history()',label='recorder')]+O.engines[5:7]
plot.plots={'pc':('sw',None,'e22')}
plot.plot()
#######################################################
## Drainage Test ###
#######################################################
##Instantiate a two-phase engine
unsat=TwoPhaseFlowEngine()
#meanDiameter=(O.bodies[-1].shape.radius + O.bodies[6].shape.radius) / 2.
##set boundary conditions, the drainage is controlled by decreasing W-phase
pressure and keeping NW-phase pressure constant
unsat.bndCondIsPressure=[0,0,1,1,0,0]
unsat.bndCondIsWaterReservoir=[0,0,1,0,0,0]
unsat.bndCondValue=[0,0,-1e8,0,0,0]
unsat.isPhaseTrapped=True #the W-phase can be disconnected from its reservoir
unsat.surfaceTension = 0.0728
unsat.initialization()
#start invasion, the data of normalized pc-sw-strain will be written into
pcSwStrain.txt
f5=open('SwPcTriax710coated3.txt',"w")
ts=O.dt
pgstep= 40
print (pgstep)
pgmax= 9000
for pg in arange(1.0e-8,pgmax,pgstep):
unsat.bndCondValue=[0,0,(-1.0)*pg,0,0,0]
unsat.invasion()
unsat.computeCapillaryForce()
unsat.meshUpdateInterval=500
unsat.defTolerance=-1
unsat.updateTriangulation=True
print(unsat.getSaturation(isSideBoundaryIncluded=False),pg,-triax.strain[1])
for b in O.bodies:
O.forces.setPermF(b.id, unsat.fluidForce(b.id))
while 1:
O.run(100,True)
unb=unbalancedForce()
if unb<0.1:
break
f5.write(str(pg)+" "+str(unsat.getSaturation(isSideBoundaryIncluded=False))+"
"+str(triax.strain[1])+"\n")
f5.close()
--
You received this question notification because your team yade-users is
an answer contact for Yade.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
