applying force on facet

Asked by ytang

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

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
Last query:
Last reply:
Revision history for this message
Klaus Thoeni (klaus.thoeni) said :
#1

Hi,

facets in Yade are massless and non-dynamic. You cannot apply forces, you can only apply velicities/displacements. If you need dynamic facets consider using pFacets.

Klaus

Revision history for this message
Best Jan Stránský (honzik) said :
#2

@ytang:

> the following code doesn't show any errors.

it does show this error:
IOError: [Errno 2] No such file or directory: 'cluster-adddd.txt'

> why the facets wall doesn't move?

The order of engines does matter.
In your case, the forces are reset before being applied (you set them by addForcestop, but then they are reset by ForceResetter).
Put addForcestop PyRunner **before** NewtonIntegrator (as it is in [2])

@Klaus
> facets in Yade are massless and non-dynamic. You cannot apply forces
this is the default state. By setting mass (and inertia) and unblocking desired DOFs, you can easily have dynamic "plain" facets.

But yes, pfacets may be a better choice for various reasons

cheers
Jan

Revision history for this message
ytang (ytang116) said :
#3

Hi Jan, Klaus,

Sorry for the late reply.

For Jan,

the file is in Dropbox.

You are right, I need to put the addForce function before the NewtonIntegrator.
Besides, I think it's very hard to control the facet, it's very to split it into pieces. So, I think I need to change to Pfacet.

best,
Yong

Revision history for this message
ytang (ytang116) said :
#4

Thanks Jan Stránský, that solved my question.

Revision history for this message
Jan Stránský (honzik) said :
#5

> the file is in Dropbox.

yes, it is clear from the link address :-) Nevertheless, it could not be publicly accessed (it is possible now), probably you will delete is some time etc. that's why there is "Please, no external links!" section in [4]

cheers
Jan

[4] https://www.yade-dem.org/wiki/Howtoask