yade-users team mailing list archive

[Hossein <question703707@xxxxxxxxxxxxxxxxxxxxx>]

New question #703707 on Yade:
https://answers.launchpad.net/yade/+question/703707

Hi 
I'm working on the saturated triaxial test. I should have positive water pressure in the simulation  and need to have its calculation. I assume 2PFV does not appropriate for this aim.
How can I use water pressure in the triaxial test without use of 2PFV?
I mention my dry-triaxial test code in below.
####################################
import matplotlib; matplotlib.rc('axes',grid=True)
import matplotlib.pyplot as plt
from yade import pack
import pylab
from numpy import *
import numpy as np
from math import *

utils.readParamsFromTable(seed=1,num_spheres=10000,compFricDegree = 5.0)
from yade.params import table

seed=table.seed
num_spheres=table.num_spheres# number of spheres
compFricDegree = table.compFricDegree # initial contact friction during the confining phase (will be decreased during the REFD compaction process)
confiningS=-1e5
rate=1


finalFricDegree1=45
## creat a packing with a specific particle side distribution (PSD)
psdSizes,psdCumm=[0.075,0.173,0.192,0.211,0.232,0.252,0.279,0.303],[0.065,0.25,0.57,0.61,0.74,0.8,0.9,1]
psdSizesArray=np.array(psdSizes)
psdSizesMeter=psdSizesArray*0.001 #Convert the size of particles to meter
sp=pack.SpherePack()
mn,mx=Vector3(0,0,0),Vector3(0.004,0.006,0.004)
sp.makeCloud(minCorner=mn,maxCorner=mx,psdSizes=psdSizesMeter,psdCumm=psdCumm,distributeMass=True,num=num_spheres,seed=seed)

key='_finalFricDegree_'+str(finalFricDegree1)

## create material #0, which will be used as default
O.materials.append(FrictMat(young=40e7,poisson=.2,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=40e7,poisson=.2,frictionAngle=0,density=0,label='frictionless'))

## create walls around the packing
walls=aabbWalls((mn,mx),thickness=0,material='frictionless')
wallIds=O.bodies.append(walls)

O.bodies.append([utils.sphere(center,rad,material='spheres') for center,rad in sp])

triax=TriaxialStressController(
	internalCompaction=True,
	goal1=confiningS,
	goal2=confiningS,
	goal3=confiningS,
	max_vel=10,
	label="triax"
)

newton=NewtonIntegrator(damping=0.7)

O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
		[Ip2_FrictMat_FrictMat_FrictPhys()],
		[Law2_ScGeom_FrictPhys_CundallStrack()]
	),
	GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
	triax,
	newton
]

while 1:
  O.run(1000,True)
  unb=unbalancedForce()
  if unb<0.01 and abs(triax.goal1-triax.meanStress)/abs(triax.goal1)<0.001:
    break

print ("Particle Distribution is Stable")

print ('Porosity:',triax.porosity)

#############################
##   REACH NEW EQU. STATE ###
#############################
finalFricDegree = finalFricDegree1 # contact friction during the deviatoric loading

#We move to deviatoric loading, let us turn internal compaction off to keep particles sizes constant
triax.internalCompaction=False
# Change contact friction (remember that decreasing it would generate instantaneous instabilities)
setContactFriction(radians(finalFricDegree))


while 1:
  O.run(1000,True)
  unb=unbalancedForce()
  if unb<0.001 and abs(triax.goal1-triax.meanStress)/abs(triax.goal1)<0.001:
    break  

triax.depth0=triax.depth
triax.height0=triax.height
triax.width0=triax.width

print ("Time for the deviator loads")

print ('Porosity:',triax.porosity)

##########################
##     Start test      ###
##########################

AxialstrainMatrix1=[]
DeviatoricStressMatrix1=[]
VolumetricStrainMatrix1=[]
PrincipalStressRatioMatrix1=[]
MeanStressMatrix1=[]

triax.stressMask=5
triax.goal2=-rate
triax.goal1=confiningS
triax.goal3=confiningS

file=open('Results'+key+'.txt',"w")
while 1:
	O.run(1000,True)
	print ('eyy=',triax.strain[1],'Sxx=', triax.stress(0))
	
	if triax.strain[1] <-0.05:
		break
##################################################33    

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