New question #700448 on Yade: https://answers.launchpad.net/yade/+question/700448
Hi Yade friends I do need run below code with two different engine including (Cundall Contact Law, Hertz Mindlin Law). At the stage of one cundall Contact must be used and then in Stage two Hertz Mindlin Contact law must be used only. Actually, I do not know how to turn off an engine at the end of stage one and then turn another one with Hertz contact. ################################### from mpl_toolkits.axes_grid1 import host_subplot 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 =8.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 young1=10e7 psdSizes,psdCumm=[0.00509,0.00704,0.0101,0.0127,0.0166,0.0222,0.0338],[.322,.376,.452,.533,.632,.726,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.0001,0.0001,0.0001) sp.makeCloud(minCorner=mn,maxCorner=mx,psdSizes=psdSizesMeter,psdCumm=psdCumm,distributeMass=True,num=num_spheres,seed=seed) sp.psd(bins=4000,mass=True) key='_Youngmadules_'+str(young1) O.materials.append(FrictMat(young=40e7,poisson=.3,frictionAngle=radians(compFricDegree),density=2600,label='spheres')) O.materials.append(FrictMat(young=40e7,poisson=.3,frictionAngle=0,density=0,label='frictionless')) 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=False, goal1=confiningS, goal2=confiningS, goal3=confiningS, max_vel=10, label="triax" ) newton=NewtonIntegrator(damping=0) O.engines=[ ForceResetter(), InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]), InteractionLoop( [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_FrictMat_FrictMat_MindlinPhys(betan=0.04,betas=0.04)], [Law2_ScGeom_FrictPhys_CundallStrack(label = 'Cundall_Law'),Law2_ScGeom_MindlinPhys_Mindlin(includeAdhesion=True,includeMoment=False,label = 'Mindlin_Law'), ] ), GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8), VTKRecorder(iterPeriod=1000,fileName="./export",recorders=['spheres','intr','color'],dead=1,label='VTK'), triax, TwoPhaseFlowEngine(label="unsat"), newton ] always_use_Cundall_Law=True always_use_Mindlin_Law=False while 1: O.run(1000,True) unb=unbalancedForce() print ('unb=', unb,'meanstress:',abs(triax.goal1-triax.meanStress)/abs(triax.goal1),'e=', triax.porosity/(1-triax.porosity), ) if unb<0.5 and abs(triax.goal1-triax.meanStress)/abs(triax.goal1)<0.5 and (triax.porosity)/(1-triax.porosity)<0.979999999999999999 : break print ("Particle.Distribution.is.stable") ################################## ################################## finalFricDegree = 9.0 triax.internalCompaction=False setContactFriction(radians(finalFricDegree)) while 1: O.run(1000,True) unb=unbalancedForce() if unb<0.5 and abs(triax.goal1-triax.meanStress)/abs(triax.goal1)<0.5: break ei0=-triax.strain[0];ei1=-triax.strain[1];ei2=-triax.strain[2] si0=-triax.stress(0)[0];si1=-triax.stress(2)[1];si2=-triax.stress(4)[2] from yade import plot O.engines=O.engines+[PyRunner(iterPeriod=20,command='history()',label='recorder')] ## a function saving variables def history(): plot.addData(e11=-triax.strain[0]-ei0, e22=-triax.strain[1]-ei1, e33=-triax.strain[2]-ei2, s11=-triax.stress(0)[0]-si0, s22=-triax.stress(2)[1]-si1, s33=-triax.stress(4)[2]-si2, ee=(((triax.porosity)/(1-triax.porosity))*1e-2), syy=((-triax.stress(3)[1])*1e-3), i=O.iter ) #plot.addData( i=O.iter,ee=(triax.porosity)/(1-triax.porosity),s22=-triax.stress(3)[1]*1e-3) #plot.addData(e22=-triax.strain[1],t=O.time,s22=-triax.stress(2)[1],p=flow.MeasurePorePressure((0.5,0.5,0.5))) ##make nice animations: #O.engines=O.engines+[PyRunner(iterPeriod=200,command='flow.saveVtk()')] from yade import plot plot.plots={'syy':('ee')} plot.plot() ################################### #Stage one ( Two phase flow ) ################################### print("Get saturation") triax.stressMask=7 triax.goal1=triax.goal3=0 #confiningS goalTop=triax.stress(3)[1] triax.goal2=goalTop triax.wall_bottom_activated=0 void=[];sy=[];Suction=[];Void2=[];valumstrain=[] meanDiameter=(O.bodies[-1].shape.radius + O.bodies[6].shape.radius) / 2. unsat.bndCondIsPressure=[0,0,1,1,0,0] unsat.bndCondValue=[0,0,-1e8,0,0,0] unsat.isPhaseTrapped=True unsat.initialization() unsat.surfaceTension = 0.0728 #file=open('Results'+key+'.txt',"w") for pg in arange(1.0e-5,2,0.7): unsat.bndCondValue=[0,0,(-1.0)*pg*unsat.surfaceTension/meanDiameter,0,0,0] unsat.invasion() unsat.computeCapillaryForce() for b in O.bodies: O.forces.setPermF(b.id, unsat.fluidForce(b.id)) while 1: O.run(1000,True) unb=unbalancedForce() if unb<0.5: break print ('PC:',((-unsat.bndCondValue[2])*1e-3), 'Sw:',(unsat.getSaturation(False)),'e:',(triax.porosity)/(1-triax.porosity), 'e33=', (triax.strain[2]), 'e22=', (-triax.strain[1]), 's22=', (-triax.stress(3)[1]*1e-3)) print ('----------LOADING--------------') ############################################################### ###Stage Two ( Loading ) ############################################################### O.engines[2].lawDispatcher.functors[1].always_use_Cundall_Law=False O.engines[2].lawDispatcher.functors[1].always_use_Mindlin_Law=True unsat.bndCondIsPressure=[0,0,0,1,0,0] triax.stressMask=2 triax.goal1=triax.goal3=0 goalTop=triax.stress(3)[1] triax.goal2=goalTop triax.internalCompaction=False triax.wall_bottom_activated=False loadingMatrix=[-5e3,-2e5,-8.3e5,-4e5] for i in arange(0,len(loadingMatrix),1): triax.goal2=loadingMatrix[i] O.run(1000,True) void.append((triax.porosity)/(1-triax.porosity)) sy.append((-triax.stress(3)[1])*1e-3) valumstrain.append((-triax.strain[1])) print ('e=',((triax.porosity)/(1-triax.porosity)), 's22=',(-triax.stress(3)[1]*1e-3), 'Load:',(loadingMatrix[i]), 'Sw:',(unsat.getSaturation(False)), 'e22=',(-triax.strain[1]), 'PC=',((-unsat.bndCondValue[2])*1e-3)) plt.figure() plt.subplot(221) plt.grid(True) plt.plot(sy,void,linestyle='-',marker="^", markersize=8,color='r',linewidth=2) plt.legend(('Exprimental data','simulation data'),loc='upper right', shadow=True) plt.show() -- 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