applying constant stress

Asked by ytang on 2020-10-28

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

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
2020-10-28
Last reply:
2020-10-28

This question was reopened

Jan Stránský (honzik) said : #1

Hello,
servo-control is displacement control, but trying to fit the prescribed force

You a few options here, e.g.:
- play with parameters of the servo-control (normalVEL, 100, 0.95, normalSTRESS, 10, 0.001, ...) to fit better your criteria
- use some more advanced servo-control algorithm (e.g. based on longer history)
- assign mass and apply force to topplate body

cheers
Jan

Can you help with this problem?

Provide an answer of your own, or ask ytang for more information if necessary.

To post a message you must log in.