yade-users team mailing list archive

[Hossein <question700314@xxxxxxxxxxxxxxxxxxxxx>]

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

    Status: Answered => Open

Hossein is still having a problem:
 Thanks Karol 
your comment was really practical 

I have another problem that when I wanna use Twophase flow engine with adhesion (gamma parameters which I must use) this engine work not very well. So, I have to have two different engine, one of them (Cundall Contact Law)  must be used by the Twophase and another one  (Hertz Mindlin )one must be used terms of loading. The way that came on my mind is that I should turn down an engine while another one must be replaced. I mention my code in below and also mark a line that I think the new engine 
(Hertz Mindlin Include Adhesion  ) must be add. 

Actually, I do need to use Hertz Mindlin (Include Adhesion with defining
surface energy (gamma)) in Stage Two in below cod.


###################################
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=70e7,poisson=.2,frictionAngle=radians(compFricDegree),density=2600,label='spheres'))
O.materials.append(FrictMat(young=70e7,poisson=.2,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.4)

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,gamma=7)],
		[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=[];kz=[]


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

betan=0.04
betas=0.04
gamma=0.7
Mindlin_Law.includeMoment=True
Mindlin_Law.includeAdhesion=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]))
  kz.append((-triax.stress(triax.wall_top_id)[1])/(-triax.stress(triax.wall_right_id)[0]))
  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]), 'kz=',((-triax.stress(triax.wall_right_id)[0])/((-triax.stress(3)[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.