Once the flow engine is used, the particles run out of the wall

Asked by Zhicheng Gao

There was no problem when I just simulated the triaxial test, but when I activated the FlowEngine, the particles moved out of the wall and the error Failed to triangulate body with id=104 ....
Can someone help with this problem? Thank you for your attention to this matter!
This is my code:

#-*- coding: utf-8 -*-
##______________ First section, generate sample_________

from __future__ import print_function
from yade import pack, qt, plot
from math import *
from builtins import range

nRead=readParamsFromTable(
        ## model parameters
        num_spheres=100,
        targetPorosity= .425,
        confiningPressure=-100000,
        ## material parameters
        compFricDegree=15,#contact friction during the confining phase
        finalFricDegree=30,#contact friction during the deviatoric loading
        young=2e8,
        poisson=.2,
        density=2600,
        alphaKr=7.5,
        alphaKtw=0,
 competaRoll=.22,
        finaletaRoll=.22,
        etaTwist=0,
        normalCohesion=0,
        shearCohesion=0,
        ## fluid parameters
        fluidDensity=1000,
        dynamicViscosity=.001,
        ## control parameters
        damp=0,
        stabilityThreshold=.001,
        ## output specifications
        filename='suffusion',
        unknowOk=True
)

from yade.params.table import *

mn,mx=Vector3(0,0,0),Vector3(0.008,0.008,0.008)

# create materials for spheres
#shear strength is the sum of friction and adhesion, so the momentRotationLaw=True
O.materials.append(CohFrictMat(alphaKr=alphaKr,alphaKtw=alphaKtw,density=density,etaRoll=competaRoll,etaTwist=etaTwist,frictionAngle=radians(compFricDegree),momentRotationLaw=True,normalCohesion=normalCohesion,poisson=poisson,shearCohesion=shearCohesion,young=young,label='spheres'))
O.materials.append(FrictMat(young=10*young,poisson=poisson,frictionAngle=0,density=0,label='walls'))
walls=aabbWalls([mn,mx],thickness=0,material='walls')
wallIds=O.bodies.append(walls)

# generate particles packing
sp=pack.SpherePack()
sp.makeCloud(mn,mx,psdSizes=[0.00008,0.000125,0.0001592,0.0002003,0.0003153,0.000399,0.000502,0.0005743],psdCumm=[0.0,0.00628,0.0565,0.198,0.721,0.915,0.991,1.0],num=num_spheres,seed=1)
sp.toSimulation(material='spheres')

#check the number of bodies
number_generate_body=len(O.bodies)
for i in range(number_generate_body):
    if not O.bodies[i]:
        number_generate_body=number_generate_body-1

O.engines=[
        ForceResetter(),
        InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
        InteractionLoop(
            [Ig2_Sphere_Sphere_ScGeom6D(),Ig2_Box_Sphere_ScGeom()],
            [Ip2_FrictMat_FrictMat_FrictPhys(),Ip2_CohFrictMat_CohFrictMat_CohFrictPhys(label='contact',setCohesionNow=False,setCohesionOnNewContacts=False)],
            [Law2_ScGeom_FrictPhys_CundallStrack(),Law2_ScGeom6D_CohFrictPhys_CohesionMoment(useIncrementalForm=True,always_use_moment_law=True)],
 ),
        FlowEngine(dead=1,label="flow"),#introduced as a dead engine for the moment, see 2nd section
        GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8),
        TriaxialStressController(label='triax',
            # specify target values and whether they are strains or stresses
            goal1=confiningPressure,goal2=confiningPressure,goal3=confiningPressure, stressMask=7,
            # type of servo-control, external walls compaction
            internalCompaction=False,
            thickness=0,
            ),
        NewtonIntegrator(damping=0)
]
qt.View()
import sys
while True:
    O.run(1000,True)
    unb=unbalancedForce()
    print('unbalanced force:',unb,' mean stress: ',triax.meanStress)
    if unb<stabilityThreshold and abs((confiningPressure-triax.meanStress)/confiningPressure)<0.001:
        break
setContactFriction(radians(finalFricDegree))
height=O.bodies[3].state.pos[1]-O.bodies[2].state.pos[1]
print(aabbDim(),height)

#B. Activate flow engine and set boundary conditions in order to get permeability
flow.dead=0
flow.defTolerance=0.3
flow.meshUpdateInterval=200
flow.useSolver=3
flow.permeabilityFactor=1
flow.viscosity=10
flow.bndCondIsPressure=[0,0,1,1,0,0]
flow.bndCondValue=[0,0,0.5,0,0,0]
flow.boundaryUseMaxMin=[0,0,0,0,0,0]
O.dt=0.1e-3
O.dynDt=False

O.run(1,1)
Qin = flow.getBoundaryFlux(2)
Qout = flow.getBoundaryFlux(3)

print("Qin=",Qin," Qout=",Qout)

#check the number of bodies
number_generate_body=len(O.bodies)
for i in range(number_generate_body):
    if not O.bodies[i]:
        number_generate_body=number_generate_body-1
print(number_generate_body)

Question information

Language:
English Edit question
Status:
Solved
For:
Yade Edit question
Assignee:
No assignee Edit question
Solved by:
Robert Caulk
Solved:
Last query:
Last reply:
Revision history for this message
Robert Caulk (rcaulk) said :
#1

Probably a time step instability. So you could test it by decreasing the timestep. Hard to know without more info (e.g. what did you change from the oedometer.py?)

Cheers,

Robert

Revision history for this message
Zhicheng Gao (zhichenggao) said :
#2

Dear Robert Caulk,
I just changed the particle size distribution and the size of the box

Revision history for this message
Zhicheng Gao (zhichenggao) said :
#3

I use the engine "GlobalStiffnessTimeStepper(active=1,timeStepUpdateInterval=100,timestepSafetyCoefficient=0.8)" to calculate the time step automatically. Generally speaking, this time step should be very small.

Revision history for this message
Zhicheng Gao (zhichenggao) said :
#4

@Robert Caulk

Dear Robert Caulk,

Thank you for your kind answer!
Can I use the engine "GlobalStiffnessTimeStepper" to automatically calculate time steps? What adjustments should be made after modifying the particle size distribution?

Best regards,
Zhicheng Gao

Revision history for this message
Best Robert Caulk (rcaulk) said :
#5

GlobalStiffnessTimeStepper does not consider the stability of the implicit flow solver. So no, you cannot use that :-)

>adjustments should be made after modifying the particle size distribution?

Decrease the time step.

Revision history for this message
Zhicheng Gao (zhichenggao) said :
#6

Thanks Robert Caulk, that solved my question.