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.
_______________________________________________
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
