Hi All,
I'm trying to do a CPT simulation[1]. I can do this successfully.
Now, I'm trying to apply constant stress at the top surface. If I did not create the penetrator, the stress will reach a constant value and keep its value. However, when I simulate the penetration process, the pressure did not keep a constant value. it changes back and forth. I apply the stress like the method in this post[2].
#####################################
Here is the code:
#####################################
from yade import pack, plot
normalSTRESS=2.e5
normalVEL=normalSTRESS/1e6
young=4e8
finalcompFricDegree=19.5
from yade import ymport
from yade import pack,export,geom
import itertools
from numpy import *
import numpy as np
from yade import pack, plot, export, utils
import math
young=4e8
damp = 0.8
finalcompFricDegree=19.5
finalFricDegree=19.5
############################################# sphere particles material properties ########################
#sphereMat=O.materials.append(CohFrictMat(young=young,poisson=0.3,frictionAngle=radians(finalcompFricDegree),isCohesive=False,alphaKr=0.2,alphaKtw=0,etaRoll=0.5,momentRotationLaw=True,density=2648))
O.bodies.append(ymport.text("servo.txt"))
#O.materials.append(FrictMat(young=4e9,poisson=0.3,frictionAngle=0,density=0,label='walls'))
#box = O.bodies.append(geom.facetBox((0.5,0.5,0.5),0.5,material='walls'))
O.bodies.append(geom.facetBox((.5,.5,.5),(.5,.5,.5),wallMask=31))
facetMat=O.materials.append(FrictMat(young=4e8,poisson=0.3,frictionAngle=radians(finalcompFricDegree),density=0))
x0=0.5;y0=0.5;z0=1
cylinderIDS= O.bodies.append(geom.facetCylinder((x0,y0,1),
radius=0.025,height=0.15,orientation=Quaternion((1, 0, 0), 0),wallMask=5,segmentsNumber=10, angleRange=None,material=facetMat))
coneIDS= O.bodies.append(geom.facetCone((x0,y0,0.9141745),
radiusTop=0.025,radiusBottom=0.0,height=0.02165,orientation=Quaternion((1, 0, 0), 0),wallMask=6,segmentsNumber=10, angleRange=None,material=facetMat))
O.bodies.append(wall((0.5,0.5,0.38),axis=2,sense=-1))
topplate = O.bodies[-1]
############################ cone and shaft geometry parameter ############################
###########################################################################################
############################## set the global engine for the simulation####################
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
),
GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.8),
NewtonIntegrator(gravity=(0,0,-9.81),damping=0.0),
TranslationEngine(velocity=0.0, translationAxis=(0.,0.,-1.), ids=[topplate.id],label='zTranslation'),
PyRunner(iterPeriod=100,initRun=True,command='servoController()',label='servor'),
PyRunner(iterPeriod=100,initRun=True,command='addPlotData()',label='plotdata'),
HelixEngine(angularVelocity=(50*(2*pi/60)),linearVelocity=-0.04,rotateAroundZero=True,zeroPoint=(x0,y0,z0),rotationAxis=(0,0,1),ids=cylinderIDS),
HelixEngine(angularVelocity=(50*(2*pi/60)),linearVelocity=-0.04,rotateAroundZero=True,zeroPoint=(x0,y0,z0),rotationAxis=(0,0,1),ids=coneIDS),
#VTKRecorder(fileName='3d-vtk-',recorders=['all','bstresses'],iterPeriod=30000),
#############################################################################
############################################################################
]
O.dt=.5*PWaveTimeStep()
def servoController():
global Fn, sigma
Fn = 0
sigma = 0
if O.iter<5000: return
if unbalancedForce()>.1: return
#Fn = sum(O.forces.f(b)[2] for b in topplate)
Fn = abs((O.forces.f(topplate.id))[2])
sigma = Fn/2
if zTranslation.velocity < normalVEL:
zTranslation.velocity+= normalVEL/100
if sigma > (0.95*normalSTRESS):
zTranslation.velocity = 10*normalVEL*((normalSTRESS-sigma)/normalSTRESS)
if abs(((normalSTRESS-sigma)/normalSTRESS)) < 0.001:
zTranslation.velocity = 0
def addPlotData():
Fz = sum(O.forces.f(c)[2] for c in coneIDS)
Fp = abs((O.forces.f(topplate.id))[2] )
zposition = O.bodies[-1].state.pos[2]
#stress = Fz/2
plot.addData(z=zposition,Fz = Fz ,F = Fp,i=O.iter)
plot.plots={'i':('F',None,('Fz','r--'))}
plot.plot()
#O.run(10000,True)
###################
The sample is here:
https://www.dropbox.com/sh/o4hu5tule2j8hdd/AACbLNjk1ytHWuUEZjJjDifTa?dl=0
I also attached two figures under this folder. You can see the trend of stress.
##############################################################
best,
Yong
references: [1] https://en.wikipedia.org/wiki/Cone_penetration_test
[2]https://answers.launchpad.net/yade/+question/693282
