yade-users team mailing list archive

[Ziyu Wang <question702097@xxxxxxxxxxxxxxxxxxxxx>]

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

Hi,
I am referring to concrete/triax.py [1] and encountered problems in simulating rock triaxial compression and measuring transverse strain. I used the method mentioned in [2] to measure transverse strain. The details are described as follows:

In uniaxial compression, this method is used to measure the transverse strain and good results are obtained. However, when the example in [1] is applied, an error may occur due to the large transverse stress I applied (prestress=-30e6 or bigger).As a result, in the whole process of axial compression, the transverse strain will first decrease to a negative value (I guess it is caused by the pressure of the large transverse stress), and then it will start to increase with the compression process, and this situation will become more obvious with the increase of the transverse stress.

I don't know if my guess is correct. If so, how can I offset the influence of transverse stress on the measurement of transverse strain?

Thanks for help!

---------------------code-----------------------
from yade import pack, ymport, plot, utils, export, timing
import numpy as np
import sys

damp=0.4
young=66.2e9
name='JCFPM_triax'
compFricDegree=30
poisson=0.522
name='cylinder'
preStress=-30e6
strainRate = -0.1
OUT=str(preStress)+'_JCFPM_triax'
r1=0.005
r2=0.008
nw=24
nh=15
rParticle=0.000731723
bcCoeff = 5
width = 0.025
height = 0.05

O.materials.append(JCFpmMat(type=1,density=2640,young=young,poisson=poisson,tensileStrength=8e6,cohesion=40e6,frictionAngle=radians(25),label='sphere'))
O.materials.append(JCFpmMat(type=1,label='wall',young=young,poisson=poisson,frictionAngle=radians(25)))

bottom,top=Vector3(0,0,0),Vector3(0,0,0.05)
radius=0.0125
sp=pack.SpherePack()
pred=pack.inCylinder(bottom,top,radius)
sp=pack.randomDensePack(pred,radius=0.0005,rRelFuzz=0,returnSpherePack=True,memoizeDb='/tmp/triax.sqlite')
spheres=sp.toSimulation(material='sphere',color=(0.9,0.8,0.6))

bot = [O.bodies[s] for s in spheres if O.bodies[s].state.pos[2]<rParticle*bcCoeff]
top = [O.bodies[s] for s in spheres if O.bodies[s].state.pos[2]>height-rParticle*bcCoeff]
vel = strainRate*(height-rParticle*2*bcCoeff)
for s in bot:
	s.shape.color = (1,0,0)
	s.state.blockedDOFs = 'xyzXYZ'
	s.state.vel = (0,0,-vel)
for s in top:
	s.shape.color = Vector3(0,1,0)
	s.state.blockedDOFs = 'xyzXYZ'
	s.state.vel = (0,0,vel)

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)), rCyl2*sin(2*pi*(r+0)/float(nw)), height*(h+0)/float(nh) )
		v2 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)), rCyl2*sin(2*pi*(r+1)/float(nw)), height*(h+0)/float(nh) )
		v3 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)), rCyl2*sin(2*pi*(r+1)/float(nw)), height*(h+1)/float(nh) )
		v4 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw)), rCyl2*sin(2*pi*(r+0)/float(nw)), height*(h+1)/float(nh) )
		f1 = facet((v1,v2,v3),color=(0,0,1),material='wall')
		f2 = facet((v1,v3,v4),color=(0,0,1),material='wall')
		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'



def lateral():
 elatTot=0.0
 nTot=0
 for b in O.bodies:
 	x=b.state.refPos[0]
 	y=b.state.refPos[1]
 	d=sqrt(pow(x,2)+pow(y,2))
 	if d > r1 and d < r2:
 		b.shape.color=(0.6,0.5,0.15)
 		xnew=b.state.pos[0]
 		ynew=b.state.pos[1]
 		dnew=sqrt(pow(xnew,2)+pow(ynew,2))
 		elat=(dnew-d)/d
 		elatTot+=elat
 		nTot+=1
 elat_avg=elatTot/nTot
 return elat_avg

def addForces():
	for f in facets:
		n = f.shape.normal
		a = f.shape.area
		O.forces.addF(f.id,preStress*a*n)

def stopIfDamaged(maxEps=0.001):
	extremum = max(abs(s) for s in plot.data['s'])
	s = abs(plot.data['s'][-1])
	e = abs(plot.data['e'][-1])
	if O.iter < 1000 or s > .5*extremum and e < maxEps:
		return
	if abs(s)/abs(extremum) < 0.5 :
		print('Simulation finished')
		presentcohesive_count = 0
		for i in O.interactions:
        		if hasattr(i.phys, 'isCohesive'):
            			if i.phys.isCohesive == True:
                			presentcohesive_count+=1
		print('the number of cohesive bond now is:',presentcohesive_count)
		print('Max stress and strain:',extremum,e)
		O.pause()

O.dt=.5*utils.PWaveTimeStep()
newton=NewtonIntegrator(damping=damp)

O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.5,label='is2aabb'),Bo1_Facet_Aabb()]),
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.5,label='ss2sc'),Ig2_Facet_Sphere_ScGeom()],
		[Ip2_JCFpmMat_JCFpmMat_JCFpmPhys()],
		[Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,Key=OUT+'_Crack',label='interactionLaw')]
	),
	PyRunner(iterPeriod=1,command="addForces()"),
	newton,
	PyRunner(command='plotAddData()',iterPeriod=10),
	PyRunner(iterPeriod=50,command='stopIfDamaged()'),
	VTKRecorder(iterPeriod=2000,initRun=True,fileName='triax/JFFPM-',recorders=['spheres','jcfpm','cracks'],Key=OUT+'_Crack',label='vtk',dead=0),

]

def plotAddData():
	elat_avg=lateral()
	f1 = sum(O.forces.f(b.id)[2] for b in top)
	f2 = sum(O.forces.f(b.id)[2] for b in bot)
	f = .5*(f2-f1)
	s = f/(pi*.25*width*width) 
	e = (top[0].state.displ()[2] - bot[0].state.displ()[2]) / (height-rParticle*2*bcCoeff)
	plot.addData(
		elateral = elat_avg,
		i = O.iter,
		s = -s,
		e = -e,
		tc = interactionLaw.nbTensCracks,
		sc = interactionLaw.nbShearCracks,
	)
	plot.saveDataTxt(OUT)

O.step()
ss2sc.interactionDetectionFactor=-1
is2aabb.aabbEnlargeFactor=-1
cohesiveCount = 0
for i in O.interactions:
	if hasattr(i.phys, 'isCohesive'):
		if i.phys.isCohesive == True:
			cohesiveCount+=1
print('the origin total number of cohesive bond is:',cohesiveCount)

plot.plots = { 'e':('s',), 'elateral':('s'),}
plot.plot()
O.run()
----------------------------------------------------------------------------


[1]https://gitlab.com/yade-dev/trunk/-/blob/master/examples/concrete/triax.py
[2]https://answers.launchpad.net/yade/+question/685862

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