New question #703024 on Yade: https://answers.launchpad.net/yade/+question/703024
Hi,everyone! I have completed a relatively good fluid-solid coupling simulation before. Recently, I wanted to add temperature effect, so I added thermalengine to the original script. However, there was a problem as described in the title: I opened the 3D view and found that the particles in the middle of the cylinder had disappeared, so I could not get the desired results. I do not know what happened,thank you for any possible help or suggestion.. ------------------code------------------------- from yade import pack, ymport, plot, utils, export, timing import numpy as np import sys #readParamsFromTable(Sy=4e6) #from yade.params import table #global Sy Sy=4e6 damp=0.4 key='_triax_base_' young=66.2e9 name='JCFPM_triax' poisson=0.522 name='cylinder' preStress=-60e6 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) bottom,top=Vector3(0,0,0),Vector3(0,0,0.05) 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(color=(0.9,0.8,0.6),material='sphere') 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=.01*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=-1 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() flow.bndCondIsTemperature=[1,1,1,1,1,1] flow.thermalEngine=True flow.thermalBndCondValue=[100,100,100,100,100,100] flow.tZero=25 flow.tZero=25 thermal.dead=1 thermal.conduction=True thermal.fluidConduction=True thermal.debug=0 thermal.thermoMech=True thermal.solidThermoMech = True thermal.fluidThermoMech = True thermal.advection=True thermal.useKernMethod=1 thermal.bndCondIsTemperature=[1,1,1,1,1,1] thermal.thermalBndCondValue=[25,25,25,25,25,25] thermal.fluidK = 0.650 thermal.fluidBeta = 2e-5 # 0.0002 thermal.particleT0 = 25 thermal.particleK = 2.0 thermal.particleCp = 710 thermal.particleAlpha = 3.0e-5 thermal.particleDensity = 2640 thermal.tsSafetyFactor = 0.8 #0.01 thermal.uniformReynolds =10 thermal.minimumThermalCondDist=0 thermal.dead=0 plot.plots = { 'e':('s',), 'elateral':('s'),} plot.plot() O.run() ######################## Best wishes! -- 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