yade-users team mailing list archive
-
yade-users team
-
Mailing list archive
-
Message #28224
[Question #702983]: TPF cylinder compacts
New question #702983 on Yade:
https://answers.launchpad.net/yade/+question/702983
Hi,
I am modelling two phase flow in a sphere particles compacted into a cylinder.
When I run my code, I get following error/messages:
1. Error calculating permeability
surfneg est 0
PassCompK = 0
meanK = -nan
globalK = -nan
maxKdivKmean*globalK = -nan
minKdivKmean*globalK = -nan
meanTubesRadius = -nan
meanDistance = -nan
PassKcorrect = 0
POS = 0 NEG = 0 pass = 0
PassSTDEV = 0
STATISTIC K
Kmoy = -nan Standard Deviation = -nan
kmin = 0 kmax = -nan
PassKopt = 0
0grains - vTotal = 0 Vgrains = 0 vPoral1 = 0
Vtotalissimo = 0 VSolidTot = 0 vPoral2 = 0 sSolidTot = 0
2.Entry saturation error!1.04611 637 9.02426e-06 4.43353e-06
Simulation is terminated because of an error in entry saturation!
3.Water pressure at: -1000 and air pressure at: 0 InitialPC: 1000
4.Error, saturation from Pc(S) curve is not correct: 1.01431 207 log:-3.92699 -3.8716 pw=-20068.8 19673.4
5.Non existing capillary pressure!Non existing capillary pressure!
How can I solve these errors?
I think might the error on the permeability calculation is causing the rest of the errors. Also the flow.getBoundaryFlux() is equal to zero.
I was wondering how the permeability is calcaluted in Yade, how is the function computePermeability() defined in flow.actionTPF [1]?
Best regards,
Mithu
1. https://gitlab.com/yade-dev/trunk/-/blob/master/pkg/pfv/TwoPhaseFlowEngine.cpp#L2843
MWE
import os
import numpy as np
import pandas as pd
from yade import pack, ymport, export
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=0
wall_pos=0
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.0006,-0.0006,-0.0014),Vector3(0.0006,0.0006,-0.00048)
#global walls
#walls=aabbWalls([mn,mx],thickness=0,material='walls')
#wallIds=O.bodies.append(walls)
if wall_pos==1:
allx,ally,allz,r=np.genfromtxt('Flow_test.txt', unpack=True) #access the packed cyl file (txt) which should be in the same folder of this code file
mnx=min(allx)*0.99 #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)*0.99
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
print ('mn xyz, mx xyz', mnx,mny,mnz,mxx,mxy,mxz)
mn,mx=Vector3(mnx,mny,mnz),Vector3(mxx,mxy,mxz)
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'))
#die=O.bodies.append(yade.geom.facetBox((0,0,0),(Tab_rad,Tab_rad,Cyl_height/5),wallMask=31,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
if wall_pos==0:
for b in O.bodies:
if isinstance(b.shape, Sphere):
continue
else:
O.bodies.erase(b.id)
if wall_pos==1:
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.2, gravity=[0, 0, 0]),]
save_filename='Flow_test.xml'
O.save(save_filename)
if wall_pos==0:
export.text('Flow_test.txt')
O.pause()
if save==1:
O.run()
#load and activate flow
if save==0:
read_filename='Flow_test.xml'
O.load(read_filename)
fCheck.dead = True # (!!!)
for b in O.bodies:
if isinstance(b.shape, Sphere):
b.dynamic = True
PCPRESSURE = 5000
flow.dead=0
PcMax = 50000.0
nrSteps = 80
dPc = PcMax / float(nrSteps)
flow.isPhaseTrapped=True #the W-phase can be disconnected from its reservoir
flow.drainageFirst=False#Unsaturated initially, first imbition
flow.isDrainageActivated=False
flow.initialWetting = False
flow.isImbibitionActivated=True #Start imbition
flow.isCellLabelActivated=True
flow.isInvadeBoundary=True
flow.isActivated=True
flow.initialPC = 1000.0 #PCPRESSURE #3100.0 * scale]
flow.bndCondIsPressure=[0,0,0,0,1,1]
flow.bndCondIsWaterReservoir = [0,0,0,0,1,0]
flow.boundaryUseMaxMin=[0,0,0,0,0,0]
#flow.initialization()
#for i in range(0,nrSteps):
flow.bndCondValue=[0,0,0,0,5000,0]
#print((PcMax-(i+1.0)*dPc))
flow.actionTPF()
#flow.savePhaseVtk()
#history()
flow.entryPressureMethod = 2 #1,2,3
flow.entryMethodCorrection = 2#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.permeabilityFactor=1
flow.viscosity = 0.001 #* (permeability / (1.7e-11)) * scale
flow.deltaTimeTruncation = 1e-7
#flow.bndCondValue=[0,0,0,0,5000,0]
flow.defTolerance=0.3
flow.meshUpdateInterval=1
#flow.debugTPF()=True
#flow.waterBoundaryPressure=2000
#flow.emulateAction()
flow.savePoreNetwork()
O.dt=1e-7
O.engines = O.engines+[PyRunner(command='Vtk_flow()', iterPeriod=50000)]
#data analysing during simulation
def Vtk_flow():
flow.savePhaseVtk()
print(flow.getSaturation(isSideBoundaryIncluded=True))
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)
for i in range(n_pore):
if flow.getCellIsTrapNW(i)=='True':
print(i)
for i in range(n_pore):
if flow.getCellIsTrapW(i)=='True':
print(i)
boundary_flux=0
for i in range(n_pore):
boundary_flux+=flow.getBoundaryFlux(i)
print(boundary_flux)
--
You received this question notification because your team yade-users is
an answer contact for Yade.
