# applying force on facet

Asked by ytang on 2021-02-17

Hi All,

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

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
 )
######################################################################
the following code doesn't show any errors.
My questions are:

(1) why the unbalanced force doesn't change? 
(2) why the facets wall doesn't move? 
##########################################################################
##########################################################################
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
############################################################################################################
#################################
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.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.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.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()
####################################################################################
####################################################################################

## Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
2021-02-19
Last query:
2021-02-19
Last reply:
2021-02-18
 Klaus Thoeni (klaus.thoeni) said on 2021-02-18: #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 Jan Stránský (honzik) said on 2021-02-18: #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 )

@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

 ytang (ytang116) said on 2021-02-19: #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

 ytang (ytang116) said on 2021-02-19: #4

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

 Jan Stránský (honzik) said on 2021-02-20: #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 

cheers
Jan

To post a message you must log in.