yade-users team mailing list archive

[steve <question696047@xxxxxxxxxxxxxxxxxxxxx>]

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

steve posted a new comment:
Hi Jan,
Sorry that I thought the simulation was 2D so I used the 2D porosity. But I forgot that stress is also obtained in 3D. 
1. After correction, the peak values do become closer to each other, but there is still some gap as well as the shape.
2. the error in volume calculation is that I didn't use aabbExtrema() to calculate the volume because if I add cement particle before the roller(moving line78,79 to line 69) , an error called 'Assertion `!swap' failed' occurred.
3. I tried the TesslationWrap() but the value is way much smaller, maybe it's due to the reason that the function also considers the rollers and supports?
line 1-59 is just making 2D packing, line 60-156 is 4-point bending model
Many thanks,
Steve

#################code attached
## compact 2D packing
from yade import pack,plot
nRead=readParamsFromTable(
	num_spheres=2000,
	compFricDegree = 30, 
	key='_triax_base_',
	unknownOk=True
)
from yade.params import table
num_spheres=table.num_spheres
key=table.key
targetPorosity = 0.55 
compFricDegree = table.compFricDegree
finalFricDegree = 30
rate=-0.02 
damp=0.2
stabilityThreshold=0.01 
young=5e6 
mn,mx=Vector3(0,-.1,0),Vector3(0.35,.1,.1)
O.switchScene();O.resetThisScene()
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([(0,-.1,0),(0.35,.1,.1)],thickness=0,material='walls')
wallIds=O.bodies.append(walls)
sp=pack.SpherePack()
sp.makeCloud((0,0,0),(0.35,0,.1),-1,0.3,num_spheres,False, 0.95,seed=1)
O.bodies.append([sphere(center,rad,material='spheres') for center,rad in sp])
triax=TriaxialStressController(
	maxMultiplier=1.+2e4/young, 
	finalMaxMultiplier=1.+2e3/young, 
	thickness = 0,
	stressMask = 7,
	internalCompaction=True,
)
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,
	TriaxialStateRecorder(iterPeriod=100,file='WallStresses'+table.key),
	newton
]
triax.goal1=triax.goal3=-500
triax.goal2=0
while 1:
  O.run(1000, True)
  unb=unbalancedForce()
  print ('unbalanced force:',unb,' mean stress: ',triax.meanStress)
  if unb<stabilityThreshold and abs(-1000/3-triax.meanStress)/1000/3<0.001:
    break
sp=SpherePack(); sp.fromSimulation()
O.switchScene()
## constuct 4-point bending model
import numpy as np
from yade import pack,plot
## define parameters
strainRate=-.02
sphereRadius=.0005
young=2.8e9
sigmaT=1e6
frictionAngle=atan(0.6)
relDuctility=17
omegaThreshold=0.999
epsCrackOnset=.004
damping=.6
damLaw = 0
relFuzz=0.3
factor=1.5
isoPrestress=0
poisson= 0.1
pred=pack.inAlignedBox((-0.4,-1,-1),(0.4,1,1)) 
idConcrete=O.materials.append(CpmMat(young=young,poisson=poisson,frictionAngle=frictionAngle,density=3120,sigmaT=sigmaT,epsCrackOnset=epsCrackOnset,isoPrestress=isoPrestress,relDuctility = relDuctility,damLaw = damLaw))
#

idRoller=O.materials.append(FrictMat(young=9e9,frictionAngle=atan(.5),poisson=.3,density=3e3,label='walls'))
## create top roller clumps
TL=O.bodies.append([sphere(center=(.125-0.175,0,.115-0.05+0.0002),radius=.015,material=idRoller)])
TR=O.bodies.append([sphere(center=(.225-0.175,0,.115-0.05+0.0002),radius=.015,material=idRoller)])
## creat bottom roller
O.bodies.append([sphere(center=(.025-0.175,0,-.015-0.05),radius=.015,material=idRoller,fixed=True)])
O.bodies.append([sphere(center=(.325-0.175,0,-.015-0.05),radius=.015,material=idRoller,fixed=True)])
#
spp=filterSpherePack(pred,sp,material=idConcrete,returnSpherePack=True)
cement=spp.toSimulation(material=idConcrete)
#### calculate the porosity for volume
volume=0
for i in cement:
  volume+=4/3*pi*pow(O.bodies[i].shape.radius,3)
porosity=1-volume/(0.3*0.0045*0.1)
print(porosity)
T=[O.bodies[i] for i in cement if (O.bodies[i].state.pos[2]<(-.05+2*O.bodies[i].shape.radius) and O.bodies[i].state.pos[0]>(-2*O.bodies[i].shape.radius) and O.bodies[i].state.pos[0]<(2*O.bodies[i].shape.radius))]
initial_dis= T[0].state.pos[0]-T[1].state.pos[0]
####grab the particle at the bottom fibre
p=[]
for i in cement:
  if (O.bodies[i].state.pos[2]<(-.05+2*O.bodies[i].shape.radius) and O.bodies[i].state.pos[0]>-.05 and O.bodies[i].state.pos[0]<.05):
    p.append(i)


O.engines=[
  ForceResetter(),
  InsertionSortCollider([
    Bo1_Sphere_Aabb(aabbEnlargeFactor=factor,label='bo1s'),
  ]),
  InteractionLoop(
    [
    Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=factor,label='ig2ss'),
    ],
    [
    Ip2_CpmMat_CpmMat_CpmPhys(),
    Ip2_FrictMat_CpmMat_FrictPhys(),
    Ip2_FrictMat_FrictMat_FrictPhys()
    ],
    [
    Law2_ScGeom_CpmPhys_Cpm(),
    Law2_ScGeom_FrictPhys_CundallStrack()
    ]
  ),
  GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.3),
  NewtonIntegrator(damping=0.5),
  CpmStateUpdater(realPeriod=.5),
  PyRunner(command='initiateStrain()',iterPeriod=1,label='move'),
  PyRunner(command='finish()',realPeriod=5),
  PyRunner(command='addPlotData()',iterPeriod=100),
]	
O.bodies[TL[0]].state.blockedDOFs='xyzXYZ'
O.bodies[TR[0]].state.blockedDOFs='xyzXYZ'
O.step() 
bo1s.aabbEnlargeFactor=1 
ig2ss.interactionDetectionFactor=1
def initiateStrain():
   O.bodies[TL[0]].state.vel[2]=strainRate
   O.bodies[TR[0]].state.vel[2]=strainRate
def finish():
  if (O.bodies[TL[0]].state.refPos[2]-O.bodies[TL[0]].state.pos[2])*60/13>0.03:    
     O.pause()
     pass
def addPlotData():
  f_up=utils.sumForces(ids=TL+TR,direction=(0,0,-1))
  ### directly obtain stress on the bottom fibre
  TW=TesselationWrapper()
  TW.setState()
  TW.computeVolumes()
  sigma=bodyStressTensors()
  tensile_stress=sum(sigma[i][0,0] for i in p) 
  tensile_stress2= sum(sigma[i][0,0]*4/3*pi*O.bodies[i].shape.radius**3/TW.volume(i) for i in p)
  plot.addData(s=-f_up*0.3/(.01*.01),e2=60/13*(O.bodies[TL[0]].state.refPos[2]-O.bodies[TL[0]].state.pos[2]),s2=tensile_stress*(1-porosity)/len(p),s3=tensile_stress2/len(p)) ## e2 is the tensile strain derived from the load point displacement
plot.plots={'e2':'s','e2 ':'s2','e2   ':'s3'} ## s is from formula fl/d^2 which neglect the width, s2 is directly obtained from bottom fibre and scaled by factor, s3 is using TesslationWrap?
plot.plot()
O.run()

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