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. _______________________________________________ Mailing list: https://launchpad.net/~yade-users Post to : yade-users@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp