yade-users team mailing list archive

[Mithushan Soundaranathan <question702694@xxxxxxxxxxxxxxxxxxxxx>]

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.