yade-users team mailing list archive

[Ziyu Wang <question703024@xxxxxxxxxxxxxxxxxxxxx>]

Question #703024 on Yade changed:
https://answers.launchpad.net/yade/+question/703024

    Status: Answered => Open

Ziyu Wang is still having a problem:
Hi Robert,

In fact, I think there may be a misunderstanding between us, or I
misunderstood your meaning(Sorry for that..)

My description in 1# is to add the Thermalengine and related parameters
on the basis of the fluid-solid coupling code I wrote myself. The
parameter values refer to the examples, but do not involve the
modification of the example script.

I'll put the fluid solid coupling code that was successfully run before
below. I don't know if it's what you need..

Best wishes.

######
from yade import pack, ymport, plot, utils, export, timing
import numpy as np
import sys
Sy=4e6
damp=0.4

key='_triax_base_'
young=66.2e9
name='JCFPM_triax'
poisson=0.522
name='cylinder'
preStress=-90e6
strainRate = -10
OUT=str(Sy)+'_JCFPM_triax'
r1=0.005
r2=0.008

nw=24
nh=15
rParticle=0.000731723
bcCoeff = 5
width = 0.025
height = 0.05
A=0.25*3.14*width*width

allx,ally,allz,r=np.genfromtxt('packing-cylinder.spheres', unpack=True)
mnx=min(allx)*1.01 #same as aabbExtrema, get the xyz of the lowest corner, multiply by factor to reduce it to let the wall surround all the spheres.
mny=min(ally)*1.01
mnz=min(allz)*0.99
mxx=max(allx)*1.01 #same as aabbExtrema, get the xyz of the top corner, multiply by factor to increase it to let the wall surround all the spheres.
mxy=max(ally)*1.01
mxz=max(allz)*1.01

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='wall1',young=young,poisson=poisson,frictionAngle=radians(25)))
O.materials.append(JCFpmMat(type=1,frictionAngle=0,density=0,label='wall2'))

mn,mx=Vector3(mnx,mny,mnz),Vector3(mxx,mxy,mxz)
walls=aabbWalls([mn,mx],thickness=0,material='wall2')
wallIds=O.bodies.append(walls)

spheres=O.bodies.append(ymport.text('packing-'+name+'.spheres',scale=1,shift=Vector3(0,0,0),material='sphere',color=(0.25,0.25,0.25)))

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='wall1')
		f2 = facet((v1,v3,v4),color=(0,0,1),material='wall1')
		facets.extend((f1,f2))

O.bodies.append(facets)
mass = O.bodies[7].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 bodyByPos(x,y,z):
 cBody = O.bodies[1]
 cDist = Vector3(100,100,100)
 for b in O.bodies:
  if isinstance(b.shape, Sphere):
   dist = b.state.pos - Vector3(x,y,z)
   if np.linalg.norm(dist) < np.linalg.norm(cDist):
    cDist = dist
    cBody = b
 return cBody

bodyOfInterest = bodyByPos(0.0125,0.0125,0.025)

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 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.wait()

O.dt=.003*utils.PWaveTimeStep()

newton=NewtonIntegrator(damping=damp)
O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.5,label='is2aabb'),Bo1_Facet_Aabb(),Bo1_Box_Aabb()]),
	InteractionLoop(
		[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.5,label='ss2sc'),Ig2_Facet_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
		[Ip2_JCFpmMat_JCFpmMat_JCFpmPhys()],
		[Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,Key=OUT+'_Crack',label='interactionLaw')]
	),
	PyRunner(iterPeriod=1,command="addForces()"),
	FlowEngine(dead=1,label="flow"),
	#ThermalEngine(dead=1,label='thermal'),
	newton,
	PyRunner(command='plotAddData()',iterPeriod=1000),
	PyRunner(iterPeriod=1000,command='stopIfDamaged()'),
]

def plotAddData():
	elat_avg=lateral()
	Qin = flow.getBoundaryFlux(4)
	perm = abs(Qin) * flow.viscosity * height / (A * Sy)
	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,
		permeability = perm
	)
	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)
flow.debug=False
flow.permeabilityMap = False
flow.defTolerance=-1
flow.meshUpdateInterval=1000
flow.fluidBulkModulus=2.2e9
flow.useSolver=4
flow.permeabilityFactor=1
flow.viscosity= 0.001
flow.decoupleForces =  False
flow.bndCondIsPressure=[1,1,1,1,1,1]
flow.bndCondValue=[0,0,0,0,0,Sy]
flow.dead=0
flow.emulateAction()

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

By the way,the packing file is generated by the method in #2.Thanks
again!

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