yade-users team mailing list archive

[ytang <question695587@xxxxxxxxxxxxxxxxxxxxx>]

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

Hi All,

I'm trying to simulate the CPT [1] in a deep location. In order to simulate the confining pressure state, I use the method mentioned here [2].

The basic idea is as follows:

(1) import the particles;
(2) create the top, bottom, and cylindrical walls by facets with mass
(3) give the facets with some forces.
In the end, I just want to test the unbalanced force.  The result shows that the unbalanced force doesn't change at all even though I set the gravity force. Besides, it seems the facets don't move to reach the target pressure. (you can see these two situations here
[3] )
######################################################################
the following code doesn't show any errors.
My questions are:

(1) why the unbalanced force doesn't change? [3]
(2) why the facets wall doesn't move? [3]
##########################################################################
##########################################################################
from yade import ymport ,pack,export,geom,plot,export,utils
import itertools
from numpy import *
import numpy as np
import math
young=4e8
finalcompFricDegree=19.5
stabilityThreshold=0.02
############################################# 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("cluster-adddd.txt",material = sphereMat)) #### the particles are here[3]
############################################################################################################
#################################
frictMat = O.materials.append(FrictMat(young=4e8,poisson=0.3,frictionAngle=19.5))
#######################################################################################################
############################################### cylindrical facets ####################################
width  = 0.4 ; height = 0.5 ; preStress = -3e5
#########################facets division
nw = 25 ; nh = 15
############################################### cylindrical facets ####################
facets = []
rCyl2 = .5*width / cos(pi/float(nw))
for r in range(nw):
	for h in range(nh):
		v1 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw))+0.2, rCyl2*sin(2*pi*(r+0)/float(nw))+0.2, height*(h+0)/float(nh) )
		v2 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw))+0.2, rCyl2*sin(2*pi*(r+1)/float(nw))+0.2, height*(h+0)/float(nh) )
		v3 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw))+0.2, rCyl2*sin(2*pi*(r+1)/float(nw))+0.2, height*(h+1)/float(nh) )
		v4 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw))+0.2, rCyl2*sin(2*pi*(r+0)/float(nw))+0.2, height*(h+1)/float(nh) )
		f1 = facet((v1,v2,v3),color=(0,0,1),material=frictMat)
		f2 = facet((v1,v3,v4),color=(0,0,1),material=frictMat)
		facets.extend((f1,f2))
O.bodies.append(facets)
mass = O.bodies[0].state.mass
for f in facets:
	f.state.mass = mass
	f.state.blockedDOFs = 'XYZz'
############################################# apply prestress to facets
def addForces():
	for f in facets:
		n = f.shape.normal
		a = f.shape.area
		O.forces.addF(f.id,preStress*a*n)
##############################################################################################################
################################### top  facets ##############################################################
widthtop  = 0.5 ; heighttop = 0.45 ; preStresstop = -3e5
########################################################  facets division
nxtop = 25 ; nytop = 25
################################################
facetstop = []
for nx in range(nxtop):
	for ny in range(nytop):
		v11 = Vector3(  (nx+0)*widthtop/nxtop-0.05, (ny+0)*widthtop/nytop-0.05, heighttop )
		v12 = Vector3(  (nx+1)*widthtop/nxtop-0.05, (ny+0)*widthtop/nytop-0.05, heighttop )
		v13 = Vector3(  (nx+1)*widthtop/nxtop-0.05, (ny+1)*widthtop/nytop-0.05, heighttop )
		v14 = Vector3(  (nx+0)*widthtop/nxtop-0.05, (ny+1)*widthtop/nytop-0.05, heighttop )
		f11 = facet((v11,v12,v13),color=(0,0,1),material=frictMat)
		f12 = facet((v11,v13,v14),color=(0,0,1),material=frictMat)
		facetstop.extend((f11,f12))
O.bodies.append(facetstop)
masstop = O.bodies[0].state.mass
for f in facetstop:
	f.state.mass = masstop
	f.state.blockedDOFs = 'XYZxy'
############################################# apply prestress to facets
def addForcestop():
	for f in facetstop:
		n = f.shape.normal
		a = f.shape.area
		O.forces.addF(f.id,preStresstop*a*n)
########################################################################################################
################################## bottom facets########################################################
widthbottom  = 0.5 ;  heightbottom = 0.0 ;  preStressbottom = 3e5
#########################facets division
nxbottom = 25 ;  nybottom = 25
facetsbottom = []
for nx in range(nxbottom):
	for ny in range(nybottom):
		v11 = Vector3(  (nx+0)*widthbottom/nxbottom-0.05, (ny+0)*widthbottom/nybottom-0.05, heightbottom )
		v12 = Vector3(  (nx+1)*widthbottom/nxbottom-0.05, (ny+0)*widthbottom/nybottom-0.05, heightbottom )
		v13 = Vector3(  (nx+1)*widthbottom/nxbottom-0.05, (ny+1)*widthbottom/nybottom-0.05, heightbottom )
		v14 = Vector3(  (nx+0)*widthbottom/nxbottom-0.05, (ny+1)*widthbottom/nybottom-0.05, heightbottom )
		f11 = facet((v11,v12,v13),color=(0,0,1),material=frictMat)
		f12 = facet((v11,v13,v14),color=(0,0,1),material=frictMat)
		facetsbottom.extend((f11,f12))
O.bodies.append(facetsbottom)
massbottom = O.bodies[0].state.mass
for f in facetsbottom:
	f.state.mass = massbottom
	f.state.blockedDOFs = 'XYZxy'
############################################# apply prestress to facets
def addForcesbottom():
	for f in facetsbottom:
		n = f.shape.normal
		a = f.shape.area
		O.forces.addF(f.id,preStressbottom*a*n)
########################################################################################################
O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
	InteractionLoop(
			[Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Facet_Sphere_ScGeom()],
			[Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys()],
					[Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=True,label='cohesiveLaw')]
	),
	GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=10,timestepSafetyCoefficient=0.8),
	NewtonIntegrator(gravity=(0,0,-9.81),damping=0.0),
	VTKRecorder(fileName='3d-vtk-',recorders=['all','bstresses'],iterPeriod=20000),
	############################################################################
	############################################################################
	PyRunner(iterPeriod=1,command="addForces()"),
	PyRunner(iterPeriod=1,command="addForcestop()"),
	PyRunner(iterPeriod=1,command="addForcesbottom()"),		
	PyRunner(iterPeriod=1000,command="checkunbalanced()"),				
]
def checkunbalanced():
	unb=unbalancedForce()
	print 'unbalanced force: ',unb
	if unb<stabilityThreshold:
		O.pause()	
export.text("cluster-final.txt")
O.save('make_sample_cluster_final-confining.yade.gz' )
O.saveTmp()
####################################################################################
####################################################################################

references:
[1] https://en.wikipedia.org/wiki/Cone_penetration_test
[2] https://gitlab.com/yade-dev/trunk/-/blob/master/examples/concrete/triax.py
[3] https://www.dropbox.com/sh/r8j8cpviz742y21/AAAMJ9opuVFWkiPLyioeetJ_a?dl=0

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