← Back to team overview

yade-users team mailing list archive

[Question #697952]: Once the flow engine is used, the particles run out of the wall

 

New question #697952 on Yade:
https://answers.launchpad.net/yade/+question/697952

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)

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.