Question #707278 on Yade changed: https://answers.launchpad.net/yade/+question/707278
Status: Answered => Open Ruidong LI is still having a problem: Thanks for your reply, Jan! Your answer solved my problem. I rewrote the code for generating the flexible membrane. Moreover, I used your suggested servo-control mode for the top, bottom, and cylindrical walls. I tried to measure whether the confining stress was reached or not. The result indicated that the confining stress in the top and bottom direction reached the target value while that in the cylindrical direction ly reached the half of target value. I don't know where the problem is. The assignment of cylindrical stress or measure? The updated MWE is presented below: from __future__ import print_function from yade import pack, ymport, qt import gts, os.path, locale, plot, random import numpy as np from yade.gridpfacet import * import math locale.setlocale( locale.LC_ALL, 'en_US.UTF-8' ) #gts is locale-dependend. If, for example, german locale is used, gts.read()-function does not import floats normally ############################################ ### 1 ### ### DEFINING VARIABLES AND MATERIALS ### ############################################ # The following lines are used parameter definitions readParamsFromTable( # material parameters young_w = 5e9, # Young's modulus of plates and walls young_g = 1.8e6, # Young's modulus of grains [1] den_w = 7874, # density of plates and walls den_g = 980, # density of grains poi_w = 0.3, # possion's ratio of plates and walls poi_g = 0.25, # possion's ratio of grains friAngle_w = 0.5,#radians(15), # friction angle of plates and walls friAngle_g = 0.5,#radians(29), # friction angle of grains # Parameters of cylindrical walls x_cyl = 0.0547, # x-coordinate of center of cylindrical walls y_cyl = 0.0535, # y-coordinate of center of cylindrical walls z_cyl = 0, # z-coordinate of center of cylindrical walls r_cyl = 0.0358, # radius of the cylinder h_cyl = 0.14, # height of the cylinder ) from yade.params.table import * # create materials for spheres and walls wallMat = O.materials.append(FrictMat(young = young_w, poisson = poi_w, frictionAngle = friAngle_w, density = den_w, label = 'walls')) sphereMat = O.materials.append(FrictMat(young = young_g, poisson = poi_g, frictionAngle = friAngle_g, density = den_g, label = 'spheres')) ############################################### ### 2 ### ### IMPORTING GRAINS AND CREATING CLUMPS ### ############################################### # spheres pred = pack.inCylinder((x_cyl, y_cyl, z_cyl), (x_cyl, y_cyl, h_cyl), r_cyl) sp = SpherePack() sp = pack.randomDensePack(pred, spheresInCell=2000, radius=0.005, returnSpherePack=True) spheres = sp.toSimulation(color=(0, 1, 1), material=sphereMat) # assign unique color for each sphere currentSphereId = O.bodies[0].id color = randomColor() for b in O.bodies: if b.id != currentSphereId:#change color each time clumpId changes currentSphereId = b.id color = randomColor() if isinstance(b.shape,Sphere):#colorize spheres b.shape.color = color # create top and bottom plates h_cyl = max([b.state.pos[2] + b.shape.radius for b in O.bodies if isinstance(b.shape, Sphere)]) h0_cyl = min([b.state.pos[2] - b.shape.radius for b in O.bodies if isinstance(b.shape, Sphere)]) idTplate = O.bodies.append(yade.geom.facetCylinder((x_cyl,y_cyl,h_cyl),radius=r_cyl*1.2,material=wallMat,height=0,segmentsNumber=40,color=(1,1,1),wire=True)) idBplate = O.bodies.append(yade.geom.facetCylinder((x_cyl,y_cyl,h0_cyl),radius=r_cyl*1.2,material=wallMat,height=0,segmentsNumber=40,color=(1,1,1))) ######################### ### 3 ### ### DEFINE ENGINES ### ######################### #target confining stress confiningStress = 5e4 global wAz wAz = pi * (1.2*r_cyl) * (1.2*r_cyl) axialforce = confiningStress * wAz # Define engine O.engines = [ ForceResetter(), InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_PFacet_Aabb(), Bo1_Facet_Aabb()], sortThenCollide=True), InteractionLoop( [ Ig2_GridNode_GridNode_GridNodeGeom6D(), Ig2_GridConnection_GridConnection_GridCoGridCoGeom(), Ig2_Sphere_Sphere_ScGeom6D(), Ig2_Sphere_PFacet_ScGridCoGeom(), Ig2_Facet_Sphere_ScGeom(), ], [ Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=True), Ip2_FrictMat_FrictMat_FrictPhys() ], [ Law2_ScGeom6D_CohFrictPhys_CohesionMoment(), Law2_ScGeom_FrictPhys_CundallStrack(), Law2_ScGridCoGeom_FrictPhys_CundallStrack(), Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(), ] ), NewtonIntegrator(gravity=(0,0,0),damping=0.7,label='newton'), ServoPIDController(axis=[0, 0, 1], maxVelocity=10, iterPeriod=100, ids=idTplate, target=axialforce, kP=1.0, kI=1.0, kD=1.0), ServoPIDController(axis=[0, 0, -1], maxVelocity=10, iterPeriod=100, ids=idBplate, target=axialforce, kP=1.0, kI=1.0, kD=1.0), PyRunner(command='addPlotData()', iterPeriod=100, label='graph'), ] ################################ ### 3 ### ### CREATE FLEXIBLE MEMBRANE ### ################################ ## Generate flexible membrane O.materials.append(CohFrictMat(young=1e7,poisson=1,density=2650,frictionAngle=radians(30),normalCohesion=3e7,shearCohesion=3e7,momentRotationLaw=True,label='gridNodeMat')) O.materials.append( FrictMat( young=1e7,poisson=0.1,density=2650,frictionAngle=radians(30),label='pFacetMat' )) global height pfacets=[] width=2*r_cyl #diameter of cylinder height=h_cyl #height of cylinder nw=40 # number of facets along perimeter nh=25 # number of facets along height rNode=width/100 #radius of grid node color1=[255./255.,102./255.,0./255.] color2=[0,0,0] color3=[248/255,248/255,168/255] color4=[156/255,160/255,98/255] rCyl2 = 0.5*width / cos(pi/float(nw)) vCenter = Vector3(x_cyl, y_cyl, 0) # create gridnodes nodesIds=[] for r in range(nw): for h in range(nh+1): v = Vector3(rCyl2*cos(2*pi*r/float(nw)),rCyl2*sin(2*pi*r/float(nw)), height*h/float(nh))+vCenter V = (O.bodies.append(gridNode(v,rNode,wire=False,fixed=False, material='gridNodeMat',color=color2)) ) # append node ids to seperate matrix for later use nodesIds.append(V) # create grid connection nodeIDmin = min(nodesIds) cylIdsV=[] cylIdsH=[] cylIdsI=[] for r in range(nw): for h in range(nh): if r == nw-1: #determine grid node V1 = nodeIDmin+(r+0)*(nh+1)+h+0 V2 = nodeIDmin+(r+0)*(nh+1)+h+1 V3 = nodeIDmin+(0+0)*(nh+1)+h+1 V4 = nodeIDmin+(0+0)*(nh+1)+h+0 else: #determine grid node V1 = nodeIDmin+(r+0)*(nh+1)+h+0 V2 = nodeIDmin+(r+0)*(nh+1)+h+1 V3 = nodeIDmin+(r+1)*(nh+1)+h+1 V4 = nodeIDmin+(r+1)*(nh+1)+h+0 #create grid connection cylIdsV.append(O.bodies.append(gridConnection(V1,V2,rNode,material='gridNodeMat',color=color3))) cylIdsH.append(O.bodies.append(gridConnection(V1,V4,rNode,material='gridNodeMat',color=color3))) cylIdsI.append(O.bodies.append(gridConnection(V1,V3,rNode,material='gridNodeMat',color=color3))) if h == nh-1: cylIdsH.append(O.bodies.append(gridConnection(V2,V3,rNode,material='gridNodeMat',color=color3))) # create PFacet for r in range(nw): for h in range(nh): if r == nw-1: #determine grid node V1 = nodeIDmin+(r+0)*(nh+1)+h+0 V2 = nodeIDmin+(r+0)*(nh+1)+h+1 V3 = nodeIDmin+(0+0)*(nh+1)+h+1 V4 = nodeIDmin+(0+0)*(nh+1)+h+0 else: #determine grid node V1 = nodeIDmin+(r+0)*(nh+1)+h+0 V2 = nodeIDmin+(r+0)*(nh+1)+h+1 V3 = nodeIDmin+(r+1)*(nh+1)+h+1 V4 = nodeIDmin+(r+1)*(nh+1)+h+0 #create Pfacets pfacets.append(O.bodies.append(pfacet(V1,V3,V2,wire=True,material='pFacetMat',color=color3))) pfacets.append(O.bodies.append(pfacet(V1,V4,V3,wire=True,material='pFacetMat',color=color4))) wAc = pi * 2 * r_cyl * h_cyl nAPFacet = wAc/(nw*nh) def node2servo(node): global x_cyl global y_cyl global nAPFacet global height global confiningStress x, y, z = node.state.pos if z == 0 or z == height: area = nAPFacet * 0.5 else: area = nAPFacet f = confiningStress * area axis = Vector3(x-x_cyl, y-y_cyl, 0) return ServoPIDController(ids=[node.id],axis=axis,iterPeriod=100,target=f,kP=1.0, kI=1.0, kD=1.0) nodes = [b for b in O.bodies if isinstance(b.shape,GridNode)] servos = [node2servo(node) for node in nodes] O.engines = O.engines + servos def addPlotData(): # axial stress global wAz global wAc global x_cyl global y_cyl fwt = Vector3(0, 0, 0) for i in idTplate: fwt += O.forces.f(i) fwb = Vector3(0, 0, 0) for j in idBplate: fwb += O.forces.f(j) # cylindrical stress fcyl = 0 for n in nodes: fnode = O.forces.f(n.id) x,y,z = n.state.pos dirV = Vector3(x-x_cyl,y-y_cyl,0).normalized() fnodeMagn = fnode.dot(dirV) fcyl += fnodeMagn plot.addData(z=O.iter, swt=fwt[2]/wAz, swb=fwb[2]/wAz, swc = fcyl/wAc) qt.View() r = qt.Renderer() r.bgColor = 1, 1, 1 plot.plots = {'z': ('swt', 'swb', 'swc')} plot.plot() O.run() -- 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