New question #702694 on Yade:
https://answers.launchpad.net/yade/+question/702694
Hi,
I am modelling the two-phase flow in spherical particles compacted into a
cylinder. I want the liquid to enter from top of the compacts. It is unsatured
initially (dry), I have tracking the saturation and pressure (the cells
(pores)) during flow there is no changes, it is the same over a longer period.
I think there something wrong with how my code are set-up and need to add more
commands, could you please help and guide me on what is wrong with my code
Best regards
MWE
import os
import numpy as np
import pandas as pd
from yade import pack, ymport
from yade import utils, plot, timing
from csv import writer
import statistics
O=Omega()
Diameter = 7.9e-5
D=Diameter
r1 = Diameter/2
rho_PH101 = 1561
k1 = 10000
kp = 140000
kc = k1 * 0.1
ks = k1 * 0.1
Chi1 = 0.34
PhiF1=0.999
fr_PH101=0.41
save=1
young=1.e6
finalFricDegree = 30 # contact friction during the deviatoric loading
errors=0
toleranceWarning =1.e-11
toleranceCritical=1.e-6
O.materials.append(LudingMat(frictionAngle=fr_PH101, density=rho_PH101, k1=k1,
kp=kp, ks=ks, kc=kc, PhiF=PhiF1, G0 = 0.0,label='mat1'))
O.materials.append(FrictMat(young=young,poisson=0.5,frictionAngle=0,density=0,label='walls'))
mn,mx=Vector3(-0.0008,-0.0008,-0.0015),Vector3(0.0008,0.0008,-0.00004)
global walls
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)
Cyl_height=0.006
Tab_rad=7*Diameter
cross_area=math.pi*(Tab_rad**2)
Comp_press=1e8
Comp_force=Comp_press*cross_area
sp=pack.SpherePack()
sp.makeCloud((-5.0*Diameter,-5.0*Diameter,-15.0*Diameter),(5.0*Diameter,5.0*Diameter,40.0*Diameter),
rMean=Diameter/2.0, num=2500)
sp.toSimulation(material='mat1')
######################################################################
die=O.bodies.append(yade.geom.facetCylinder((0,0,0),radius=Tab_rad,height=30*Diameter,segmentsNumber=20,wallMask=6,material='mat1'))
O.engines = [
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.05),
Bo1_Wall_Aabb(),
Bo1_Facet_Aabb()
]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.05),
Ig2_Facet_Sphere_ScGeom(),
Ig2_Wall_Sphere_ScGeom()],
[Ip2_LudingMat_LudingMat_LudingPhys()],
[Law2_ScGeom_LudingPhys_Basic()]
),
NewtonIntegrator(damping=0.1, gravity=[0, 0, -9.81]),
PyRunner(command='checkForce()', realPeriod=1, label="fCheck"),
TwoPhaseFlowEngine(dead=1,label="flow"),
]
def checkForce():
if O.iter < 2200000:
return
timing.reset()
highSphere = 0.0
for b in O.bodies:
if highSphere < b.state.pos[2] and isinstance(b.shape, Sphere):
highSphere = b.state.pos[2]
else:
pass
O.bodies.append(wall(highSphere+0.5*Diameter, axis=2, sense=-1,
material='mat1'))
global plate
plate = O.bodies[-1] # the last particles is the plate
plate.state.vel = (0, 0, -2)
fCheck.command = 'unloadPlate()'
def unloadPlate():
if abs(O.forces.f(plate.id)[2]) > Comp_force:
plate.state.vel *= -1
fCheck.command = 'stopUnloading()'
def stopUnloading():
if abs(O.forces.f(plate.id)[2]) == 0:
O.bodies.erase(plate.id)
for b in O.bodies:
if isinstance(b.shape, Sphere):
#b.state.blockedDOFs = 'xyXY'
r=b.shape.radius
oldm=b.state.mass
oldI=b.state.inertia
m=oldm*3./4./r
b.state.mass=m
b.state.inertia[0] = 15./16./r*oldI[0] #inertia with respect
to x and y axes are not used and the computation here is wrong
b.state.inertia[1] = 15./16./r*oldI[1] #inertia with respect
to x and y axes are not used and the computation here is wrong
b.state.inertia[2] = 15./16./r*oldI[2] #only inertia with
respect to z axis is usefull
for b in O.bodies:
if isinstance(b.shape, Sphere):
continue
elif b.id<max(wallIds)+1:
continue
else:
O.bodies.erase(b.id)
fCheck.command = 'Updates_flow()'
def Updates_flow():
O.engines=O.engines[0:3]+O.engines[4:]
O.engines=O.engines+[NewtonIntegrator(damping=0.1, gravity=[0, 0, 0]),]
save_filename='Flow_test.xml'
O.save(save_filename)
#O.pause()
fCheck.command = 'Lq_flow()'
def Lq_flow():
O.dt=1e-7
fCheck.dead = True
PCPRESSURE = 5000
flow.dead=0
flow.entryPressureMethod = 2 #2
flow.entryMethodCorrection = 2.0 #2
flow.maximumRatioPoreThroatoverPoreBody = 0.9#0.30 #0.9
flow.surfaceTension = 0.072 #/ 0.72 #0.042 * cos(51.0 * 3.14 / 180.0) *
flow.truncationValue = 1e-6
flow.truncationPrecision = 1e-6
flow.useSolver = 3
flow.viscosity = 0.001 #* (permeability / (1.7e-11)) * scale
flow.deltaTimeTruncation = 1e-10
flow.initialPC = 1000.0 #PCPRESSURE #3100.0 * scale
flow.bndCondIsPressure=[0,0,0,0,1,1]
flow.bndCondValue=[0,0,0,0,-5000,0]
flow.bndCondIsWaterReservoir = (False,False,False,False,True,True)
flow.isPhaseTrapped=True
flow.drainageFirst=False#Unsaturated initially, first imbition
flow.isDrainageActivated=False
flow.isImbibitionActivated=True
flow.isCellLabelActivated=True
flow.isActivated=True
flow.initialization()
flow.savePoreNetwork()
O.engines = O.engines+[PyRunner(command='Vtk_flow()', iterPeriod=10000)]
def Vtk_flow():
flow.savePhaseVtk()
print(flow.getSaturation(isSideBoundaryIncluded=False))
press_now=0
n_pore=flow.nCells()
for i in range(n_pore):
press_now+=flow.getCellPressure(i)
print(press_now)
--
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
