yade-users team mailing list archive

[ytang <question693681@xxxxxxxxxxxxxxxxxxxxx>]

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

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

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