yade-users team mailing list archive

Thread
Date
[Question #695192]: wall.append(O.bodies.append(pfacet) in triaxial test

To: yade-users@xxxxxxxxxxxxxxxxxxx
From: Hanying Zhang <question695192@xxxxxxxxxxxxxxxxxxxxx>
Date: Tue, 26 Jan 2021 14:31:10 -0000
Reply-to: question695192@xxxxxxxxxxxxxxxxxxxxx
Sender: bounces@xxxxxxxxxxxxx
New question #695192 on Yade:
https://answers.launchpad.net/yade/+question/695192

Hi,

I am doing a triaxial test simulating fiber-reinforced sand in which cylinder and Pfacet are used as fiber and wall(triax.wall). After setting strain to be -0.9 I find that Pfacet doesn't move along coordinate axis as regular triaxial wall. Pfacet can deform and try to fold each other. I need pfacet to move along coordinate axis. can anyone help me with that? I am using ubuntu18.04 and my yade version is 2018.02b.

###################################################
from yade.gridpfacet import *
from yade import pack, plot
#from random import random as r
import numpy as np
from numpy import *
import math

#### Parameter ####
mi,ma = (-100e-3,-100e-3,-200e-3),(100e-3,100e-3,200e-3)
young=4.0e6
poisson=3
density=1e3
color=[0.,1.,1.]
stabilityThreshold=0.01
#=============================meterials========================================
O.materials.append(FrictMat(young=4.0e6,poisson=.3,frictionAngle=30,density=2630e0,label='sphereMat'))#for sphere

O.materials.append(FrictMat(young=3e9,poisson=.15,frictionAngle=20,density=910e+0,label='extcylMat'))#for sphere-cylinder
O.materials.append(CohFrictMat(young=3e9,poisson=.15,density=910e0,frictionAngle=20,normalCohesion=1e40,shearCohesion=1e40,momentRotationLaw=True,label='intcylMat'))#for cylinder-cylinder

O.materials.append( CohFrictMat( young=3e6,poisson=0.15,density=910e0,frictionAngle=20,normalCohesion=3e100,shearCohesion=3e100,momentRotationLaw=True,label='gridNodeMat' ) )#for gridNodes
#O.materials.append(CohFrictMat(young=3e9,poisson=.15,density=910e6,frictionAngle=20,normalCohesion=1e40,shearCohesion=1e40,momentRotationLaw=True,label='gridNodeMat'))#for gridNodes
O.materials.append(FrictMat(young=4e6,poisson=0.3,density=1000,frictionAngle=20,label='pFacetMat')) #for pfacet

O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=2630,label='walls'))

triax=TriaxialStressController(
	
	maxMultiplier=1.+2e4/young, # spheres growing factor (fast growth)
	finalMaxMultiplier=1.+2e3/young, # spheres growing factor (slow growth)
	thickness = 0,
	stressMask = 0,
	max_vel = 3.7e-7,
	wall_front_activated=True,
	wall_back_activated=True,
	wall_top_activated=True,
	wall_bottom_activated=True,
	wall_left_activated=True,
	wall_right_activated=True,

	internalCompaction=False, # If true the confining pressure is generated by growing particles
)

#==============================Engines=========================================
O.engines=[
	ForceResetter(),
	InsertionSortCollider([
		Bo1_Sphere_Aabb(),
		Bo1_GridConnection_Aabb(),
		Bo1_PFacet_Aabb(),
		Bo1_Box_Aabb(),
		Bo1_Wall_Aabb()
	]),
	InteractionLoop([
		Ig2_Sphere_Sphere_ScGeom(),
		Ig2_Box_Sphere_ScGeom(),
		Ig2_GridNode_GridNode_GridNodeGeom6D(),
		Ig2_Sphere_GridConnection_ScGridCoGeom(),
		Ig2_GridConnection_GridConnection_GridCoGridCoGeom(),
		Ig2_GridConnection_PFacet_ScGeom(),
		Ig2_PFacet_PFacet_ScGeom(),
		Ig2_Sphere_PFacet_ScGridCoGeom(),
		Ig2_Wall_PFacet_ScGeom(),
		Ig2_Wall_Sphere_ScGeom()
	],
	[
		Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(setCohesionNow=True,setCohesionOnNewContacts=False),	# internal cylinder physics
		Ip2_FrictMat_FrictMat_FrictPhys(),	# physics for external interactions, i.e., cylinder-cylinder, sphere-sphere, cylinder-sphere
		Ip2_FrictMat_CpmMat_FrictPhys(),
	],
	[
		Law2_ScGeom_FrictPhys_CundallStrack(),	# contact law for sphere-sphere
		Law2_ScGridCoGeom_FrictPhys_CundallStrack(),	# contact law for cylinder-sphere
		Law2_ScGeom6D_CohFrictPhys_CohesionMoment(),	# contact law for "internal" cylinder forces
		Law2_GridCoGridCoGeom_FrictPhys_CundallStrack(),	# contact law for cylinder-cylinder interaction
		#Law2_GridCoGridCoGeom_CohFrictPhys_CundallStrack(),
		Law2_ScGeom_CpmPhys_Cpm(),
	]
	),
	GlobalStiffnessTimeStepper(timestepSafetyCoefficient=0.1,label='ts'), 
	triax,
	NewtonIntegrator(gravity=(0,0,0),damping=0.5,label='newton'),
	]

triax.goal1=triax.goal2=triax.goal3=-0.9

RNode=0.00002
nodes=[]
wall=[]
nodes.append(O.bodies.append( gridNode([-100e-3, 100e-3,-200e-3],RNode,wire=False,fixed=True,material='gridNodeMat',color=color) ))
nodes.append(O.bodies.append( gridNode([ 100e-3, 100e-3,-200e-3],RNode,wire=False,fixed=True,material='gridNodeMat',color=color) ))
nodes.append(O.bodies.append( gridNode([ 100e-3,-100e-3,-200e-3],RNode,wire=False,fixed=True,material='gridNodeMat',color=color) ))
nodes.append(O.bodies.append( gridNode([-100e-3,-100e-3,-200e-3],RNode,wire=False,fixed=True,material='gridNodeMat',color=color) ))

nodes.append(O.bodies.append( gridNode([-100e-3, 100e-3, 200e-3],RNode,wire=False,fixed=True,material='gridNodeMat',color=color) ))
nodes.append(O.bodies.append( gridNode([ 100e-3, 100e-3, 200e-3],RNode,wire=False,fixed=True,material='gridNodeMat',color=color) ))
nodes.append(O.bodies.append( gridNode([ 100e-3,-100e-3, 200e-3],RNode,wire=False,fixed=True,material='gridNodeMat',color=color) ))
nodes.append(O.bodies.append( gridNode([-100e-3,-100e-3, 200e-3],RNode,wire=False,fixed=True,material='gridNodeMat',color=color) ))

O.bodies.append(gridConnection(0,1,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(0,2,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(0,3,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(0,4,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(0,5,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(0,7,RNode,color=color,material='gridNodeMat'))

O.bodies.append(gridConnection(6,1,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(6,2,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(6,3,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(6,4,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(6,5,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(6,7,RNode,color=color,material='gridNodeMat'))

O.bodies.append(gridConnection(2,1,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(2,3,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(2,5,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(2,7,RNode,color=color,material='gridNodeMat'))

O.bodies.append(gridConnection(4,1,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(4,3,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(4,7,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(4,5,RNode,color=color,material='gridNodeMat'))

O.bodies.append(gridConnection(3,7,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(5,1,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(5,7,RNode,color=color,material='gridNodeMat'))
O.bodies.append(gridConnection(3,1,RNode,color=color,material='gridNodeMat'))

wall.append(O.bodies.append(pfacet(nodes[0],nodes[1],nodes[2],wire=True,material='pFacetMat',color=color)))
wall.append(O.bodies.append(pfacet(nodes[0],nodes[2],nodes[3],wire=True,material='pFacetMat',color=color)))

wall.append(O.bodies.append(pfacet(nodes[0],nodes[4],nodes[5],wire=True,material='pFacetMat',color=color)))
wall.append(O.bodies.append(pfacet(nodes[5],nodes[1],nodes[0],wire=True,material='pFacetMat',color=color)))

wall.append(O.bodies.append(pfacet(nodes[0],nodes[4],nodes[7],wire=True,material='pFacetMat',color=color)))
wall.append(O.bodies.append(pfacet(nodes[0],nodes[7],nodes[3],wire=True,material='pFacetMat',color=color)))

wall.append(O.bodies.append(pfacet(nodes[6],nodes[5],nodes[4],wire=True,material='pFacetMat',color=color)))
wall.append(O.bodies.append(pfacet(nodes[6],nodes[4],nodes[7],wire=True,material='pFacetMat',color=color)))

wall.append(O.bodies.append(pfacet(nodes[6],nodes[5],nodes[1],wire=True,material='pFacetMat',color=color)))
wall.append(O.bodies.append(pfacet(nodes[6],nodes[1],nodes[2],wire=True,material='pFacetMat',color=color)))

wall.append(O.bodies.append(pfacet(nodes[6],nodes[7],nodes[3],wire=True,material='pFacetMat',color=color)))
wall.append(O.bodies.append(pfacet(nodes[6],nodes[3],nodes[2],wire=True,material='pFacetMat',color=color)))


#=======================================sphere=================================================
sp = yade.pack.SpherePack()
sp.makeCloud(mi,ma,porosity=0.66/1.66,psdSizes=[0.5e-3,0.75e-3,0.8e-3,0.9e-3,1.2e-3,2e-3],psdCumm=[0.,0.2,0.4,0.6,0.8,1.0],num=133)
spheres=sp.toSimulation(color=(0,0.3,0.7),material='sphereMat')

#=======================================fiber=================================================
#below codes just to generate cylinders that don't overlap with each other.
#I tried to simplify it somehow error occurred in terminal.
#I think below codes have nothing to do with my problem. Thanks!
length=12e-3
diameter=0.3e-3
Xw=0.25e-2 # fiber content
fiber_vol=(80e-3*(40e-3*0.5)**2*np.pi)*Xw
Vfiber=length*np.pi*(0.5*diameter)**2
numFiber=fiber_vol/Vfiber
print "the numer of fibres =",int(numFiber)
Ne=int(length/(0.9e-3))
print 'Ne=',Ne
nodesIds=[]
cylIds=[]
numFibre=0
target_num=int(numFiber)
target_num=int(6)
####fiberfilewrite
fiberwr=open('fibers.txt','w')
cylNodes=[]
time=0
HalleluYah=0
for n in range(1000):
    if numFibre<target_num:
		random.seed(n*5)
		z0=random.uniform(-200e-3, 200e-3 )
		random.seed( )
		sita0=random.uniform(-np.pi,np.pi)
		random.seed( n+7)
		l_y=random.uniform(-100e-3,100e-3)
		x0=l_y*sin(sita0)
		y0=l_y*cos(sita0)
		
		random.seed(n+8)
		sita=random.uniform(-0.5*np.pi, 0.5*np.pi)
		random.seed()
		phi=random.uniform(0, 2*np.pi)  
		
		hz=length*sin(sita)
		hx=length*cos(sita)*cos(phi)
		hy=length*cos(sita)*sin(phi)
		
		x1=x0+hx
		y1=y0+hy
		z1=z0+hz

		if 0<abs(x1)<100e-3 and 0<abs(y1)<100e-3 and 0<abs(z1)<200e-3 and 0<abs(x0)<100e-3 and 0<abs(y0)<100e-3 and 0<abs(z0)<200e-3 :
			if len(cylNodes)<=0:
			    fiberwr.write(str(numFibre)+'\t'+str(x0)+'\t'+str(y0)+'\t'+str(z0)+'\t'+str(x1)+'\t'+str(y1)+'\t'+str(z1)+'\n')
			    cylNodes.append([x0,y0,z0,x1,y1,z1])
			    numFibre+=1
			    vertices=[]
			    for i in range(0, Ne+1):
				px=float(i)*hx/float(Ne)+x0; py=float(i)*hy/float(Ne)+y0; pz=float(i)*hz/float(Ne)+z0;
				vertices.append([px,py,pz])        
			    cylinderConnection(vertices,0.15e-3,nodesIds,cylIds,color=[0.5,0.5,0],fixed=False,highlight=False,intMaterial='intcylMat',extMaterial='extcylMat')
			else :
				valid=True
				for c in cylNodes:
					if valid:
						A1=np.array([x0,y0,z0])
						B1=np.array([x1,y1,z1])
						r1=0.15e-3
						A2=np.array([c[0],c[1],c[2]])	  
						B2=np.array([c[3],c[4],c[5]])
						r2=0.15e-3
						l1=B1-A1
						l2=B2-A2
						l3=A2-A1
						
						delta1=(np.dot(l1,l2))/(np.linalg.norm(l1)*np.linalg.norm(l2))
						
						if abs(abs(delta1)-1)<0.0000000000000001:
	    						p=np.linalg.norm(np.cross(l2,l3))/np.linalg.norm(l2)
	    						if p<r1+r2 :
								print("PARALLEL INTERSECT!")
								valid=False
	    						else :
								print("PARALLEL DO NOT INTERSECT!")
								valid=True
						else :
							l4=np.cross(l1,l2)
							ll1=np.cross(l1,l4)
	    						ll2=np.cross(l2,l4)

					    		A11=np.matrix([ll1,[l2[1],-l2[0],0.0],[0.0,l2[2],-l2[1]]])
					    		x11=np.dot(ll1,A1)
					    		x12=l2[1]*A2[0]-l2[0]*A2[1]
					    		x13=l2[2]*A2[1]-l2[1]*A2[2]
					    		b1=np.matrix([[x11],[x12],[x13]])
					    		C1=np.linalg.solve(A11,b1)
						
					    		A21=np.matrix([ll2,[l1[1],-l1[0],0.0],[0.0,l1[2],-l1[1]]])
					    		x21=np.dot(ll2,A2)
					    		x22=l1[1]*A1[0]-l1[0]*A1[1]
					    		x23=l1[2]*A1[1]-l1[1]*A1[2]
					    		b2=np.matrix([[x21],[x22],[x23]])
					    		C2=np.linalg.solve(A21,b2)

					    		t1=(C1[0]-A1[0])/(B1[0]-A1[0])
					    		t2=(C2[0]-A2[0])/(B2[0]-A2[0])
					    		
					    		if t1<0 :
								E=A1
					    		elif t1>=0 and t1<=1 :
								E=C1
						    	else :
								E=B1

					    		if t2<0 :
								F=A2
					    		elif t2>=0 and t2<=1 :
								F=C2
					    		else :
								F=B2
					    		D=np.linalg.norm(E-F)
					    		
					    		if D<(r1 + r2) :
								print("INTERSECT!")
								valid=False
								
					    		else :
								valid=True
								
				if valid:
					print(n,valid)
					fiberwr.write(str(numFibre)+'\t'+str(x0)+'\t'+str(y0)+'\t'+str(z0)+'\t'+str(x1)+'\t'+str(y1)+'\t'+str(z1)+'\n')
			    		cylNodes.append([x0,y0,z0,x1,y1,z1])
			    		numFibre+=1
			    		vertices=[]
			   		for i in range(0, Ne+1):
						px=float(i)*hx/float(Ne)+x0; py=float(i)*hy/float(Ne)+y0; pz=float(i)*hz/float(Ne)+z0;
						vertices.append([px,py,pz])        
			   		cylinderConnection(vertices,0.15e-3,nodesIds,cylIds,color=[0.5,0.5,0],fixed=False,highlight=False,intMaterial='intcylMat',extMaterial='extcylMat')
			   		#print(n,"A1",x0,y0,z0, "B1",x1,y1,z1)

fiberwr.close()

Thanks!


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