Yade-OpenFoam-coupling delete the particles outside the fluid cells

Hello,

Yes it is possible, although what kind of boundary conditions are used in the fluid problem?
One way to do this would be to add a PyRunner function.

In the function :

* check for particles that has lost the cohesion, exclude these ids.
* reset the fluidCoupling.setNumParticles, with the number of cohesive bonds,
* put the ids of the particles you want in fluidCoupling.setIdList

The PyRunner should be set before the fluidCoupling engine, let me know if this works.

Also, could you share a simple script (and the name of the openfoam solver), so I can test this?

Hi Deepak, thank you for your message. I'm not too sure how to put the pyrunner before the fluidcoupling engine, I was also wondering how to reset the aabbEnlargeFactor after the first iteration. When I was working with just Yade, I reset this parameter after O.step and then use the gui to run the rest of the simulation.

This is my yade script

from __future__ import print_function
import sys
from yadeimport import *
from yade.utils import *
from yade import ymport

initMPI() #Initialize the mpi environment, always required.
fluidCoupling = yade.FoamCoupling(); #Initialize the engine
fluidCoupling.getRank(); #part of Initialization.

#example of spheres in shear flow : two-way point force coupling
class simulation():

  def __init__(self):

    O.periodic = True;

    #proppant properties
    FrictAng_p = 0.9
    Density_p = 2650
    Young_p = 100e6
    TensileStr_p=3000
    Cohesion_p=3000

    proppant = JCFpmMat(type=2,young=Young_p,frictionAngle=FrictAng_p,density=Density_p,tensileStrength=TensileStr_p,cohesion=Cohesion_p,label='proppant')

    O.materials.append(proppant)

    proppant_assembly = O.bodies.append(ymport.textExt('new_prop.txt','x_y_z_r',shift=Vector3(0,0,0),scale=1,material=proppant,color=(1.00,0.67,0.50)))

    for b in proppant_assembly:
      if O.bodies[b].state.pos[0] < 0:
        O.bodies.erase(b)

    sphereIDs = [b.id for b in O.bodies if type(b.shape)==Sphere and b.material.label=='proppant']

    #coupling engine settings

    fluidCoupling.setNumParticles(len(sphereIDs))
    fluidCoupling.setIdList(sphereIDs)
    fluidCoupling.isGaussianInterp=False; #use pimpleFoamYade for gaussianInterp

    # Integrator
    newton=NewtonIntegrator(damping=0.0, gravity = (0.0 ,0.0, 0.0))
     # add small damping in case of stability issues.. ~ 0.1 max, also note : If gravity is needed, set it in constant/g dir.

    O.engines=[
  ForceResetter(),
  InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.2, label="bols"),Bo1_Facet_Aabb()], allowBiggerThanPeriod=True),
  InteractionLoop(
    [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.2, label="ig2s"),Ig2_Facet_Sphere_ScGeom()],
    [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1)],
        [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,label='lawFunctor')]
  ),
  GlobalStiffnessTimeStepper(timestepSafetyCoefficient=0.7, label = "ts"),
        fluidCoupling, #to be called after timestepper
        PyRunner(command='sim.printMessage()', iterPeriod= 1000, label='outputMessage'),
  newton,
        #PyRunner(command='sim.deletePar()',iterPeriod=50, label='checkPar'),
        VTKRecorder(fileName='yadep/3d-vtk-',recorders=['spheres','colors'],iterPeriod=1000)
    ]

    if O.iter>1:
      bo1s.aabbEnlargeFactor = 1
      ig2s.interactionDetectionFactor = 1
      print('-------------------reset aabbEnlargeFactor------------------')

  def printMessage(self):

     print("********************************YADE-ITER = " + str(O.iter) +" **********************************")
     if O.iter == 4000:
         maxVel = 0.05
         for b in O.bodies:
             if type(b.shape)==Sphere:
                 bodyVel = abs(b.state.vel.norm())
                 if bodyVel > maxVel:
                     raise ValueError("Body velocity exceeds imposed shear velocity by ", abs(bodyVel-maxVel))

  def deletePar(self):
    print("******YADE-ITER = " + str(O.iter) +" **********")
    for b in proppant_assembly:
      if O.bodies[b].state.pos[0]>0.02:
        print('delete id '+str(b))
        O.bodies.erase(b)
    fluidCoupling.setNumParticles(len(sphereIDs))
    fluidCoupling.setIdList(sphereIDs)

  def irun(self,num):
      O.run(num,1)

if __name__=="__main__":
    sim = simulation()
    sim.irun(5000)
    # print("body id = ", O.bodies[34].id)
    fluidCoupling.killMPI()

import builtins
builtins.sim=sim

proppant.txt
0.00643783930395 0.0153571235196 0.00996462421737 0.000296398704259
0.0063978094014 0.0151208968385 0.00923970556825 0.000295898392322
0.00806594949045 0.0152558278563 0.00630458365066 0.00029536865027
0.00773116324663 0.0151764265584 0.0076636891438 0.000295015488902
0.00896207660218 0.0153221567143 0.0091664883318 0.000294397456509

blockMeshDict:

FoamFile
{
    version 2.0;
    format ascii;
    class dictionary;
    object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

convertToMeters 1;

vertices
(
    (-0.002 0.0132 -0.002)
    (0.022 0.0132 -0.002)
    (0.022 0.0174 -0.002)
    (-0.002 0.0174 -0.002)
    (-0.002 0.0132 0.022)
    (0.022 0.0132 0.022)
    (0.022 0.0174 0.022)
    (-0.002 0.0174 0.022)
);

blocks
(
    hex (0 1 2 3 4 5 6 7) (12 3 12) simpleGrading (1 1 1)
);

edges
(
);

boundary
(
    top
    {
        type wall;
        faces
        (
            (3 7 6 2)
        );
    }

    bottom
    {
        type wall;
        faces
        (
            (1 5 4 0)
        );
    }

    inlet
    {
        type patch;
        faces
        (
            (0 4 7 3)
        );
    }

    outlet
    {
        type patch;
        faces
        (
            (2 6 5 1)
        );
    }

    front_back
    {
        type empty;
        faces
        (
            (0 3 2 1)
            (4 5 6 7)
        );
    }

);

mergePatchPairs
(
);

0/p
FoamFile
{
    version 2.0;
    format ascii;
    class volScalarField;
    object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions [0 2 -2 0 0 0 0];

internalField uniform 0;

boundaryField
{

    inlet
    {
        type fixedValue;
        value uniform 64;
    }

    outlet
    {
        type fixedValue;
        value uniform 0;
    }

    top
    {
        type zeroGradient;
    }

    bottom
    {
        type zeroGradient;
    }

    front_back
    {
        type empty;
    }

}

fvSolutions
solvers
{
    p
    {

        solver PCG;
        preconditioner DIC;
        tolerance 1e-06;
        relTol 0;
    }

    pFinal
    {
        $p;
        relTol 0;
    }

    U
    {

        solver PBiCG;
        preconditioner DILU;
        tolerance 1e-05;
        relTol 0;
    }
}

PISO
{
    nCorrectors 2;
    nNonOrthogonalCorrectors 0;
    pRefCell 0;
    pRefValue 0;
}

FoamFile
{
    version 2.0;
    format ascii;
    class volVectorField;
    object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions [0 1 -1 0 0 0 0];

internalField uniform (0 0 0);

boundaryField
{
    inlet
    {
        type zeroGradient;
    }

    outlet
    {
        type zeroGradient;
    }

    top
    {
        type fixedValue;
        value uniform (0 0 0);
    }

    bottom
    {
        type fixedValue;
        value uniform (0 0 0);
    }

    front_back
    {
        type empty;
    }

}

Forgot to add that the solver was icoFoam

Hello,

Modifications in the source code have been made to insert/delete bodies. You can clone this branch[1] and verify it.
I would suggest to make the following changes in your Yade script :

In function deletePar(self):

  def deletePar(self):
    print("******YADE-ITER = " + str(O.iter) +" **********")
    ids = fluidCoupling.getIdList() # get the ids of bodies in existing coupling
    for b in ids:
      if O.bodies[b].state.pos[0]>0.02:
        print('delete id '+str(b))
        O.bodies.erase(b)
        fluidCoupling.eraseId(b); # erase a specific id from the coupling

This function is called in O.engines as :

    O.engines=[ O.engines=[
  ForceResetter(),
  InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.2, label="bols"),Bo1_Facet_Aabb()], allowBiggerThanPeriod=True),
  InteractionLoop(
    [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.2, label="ig2s"),Ig2_Facet_Sphere_ScGeom()],
    [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1)],
        [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,label='lawFunctor')]
  ),
  GlobalStiffnessTimeStepper(timestepSafetyCoefficient=0.7, label = "ts"),
  PyRunner(command='sim.deletePar()',iterPeriod=10, label='checkPar'),
        fluidCoupling, #to be called after timestepper
        PyRunner(command='sim.printMessage()', iterPeriod= 1000, label='outputMessage'),
  newton,

        VTKRecorder(fileName='yadep/3d-vtk-',recorders=['spheres','colors'],iterPeriod=1000)
    ]

To use these options (at present/immediately) you will have to clone the Yade repository and switch to branch : FoamCouplingOptions[1] and compile also pull the latest updates from Yade-OpenFOAM-coupling repository[2], (git pull origin, followed by Allwmake)

[1] https://gitlab.com/yade-dev/trunk/tree/FoamCouplingOptions
[2] https://github.com/dpkn31/Yade-OpenFOAM-coupling

Hi Deepak,

Thank you so much for your help so far. I've modified my script based on your comment to check if a body does not have any cohesive bonds in its interactions, if true then delete this sphere from fluidcoupling IDs. In the first openFoam time step, some spheres were deleted as expected, however, after the deletions and in the same foam timestep, it seems the program hangs after 500 yade iterations and in the fluidCoupling step (I have put the particle deletion pyrunner before calling fluidCoupling). I have copied my python script below. I am sorry it's probably difficult to run this script as it requires two other txt files that stored the sphere locations. I was wondering if it was due to some obvious errors in my script that I haven't noticed?

from __future__ import print_function
import sys
from yadeimport import *
from yade.utils import *
from yade import ymport

initMPI() #Initialize the mpi environment, always required.
fluidCoupling = yade.FoamCoupling(); #Initialize the engine
fluidCoupling.getRank(); #part of Initialization.

#example of spheres in shear flow : two-way point force coupling
class simulation():

  def __init__(self):

    O.periodic = True;

    #proppant properties
    FrictAng_p = 0.9
    Density_p = 2650
    Young_p = 100e6
    TensileStr_p=3000
    Cohesion_p=3000

    Young = 57e8 #apparent modulus
    FrictAng = 0.5
    Density = 2650
    Poisson = 0.28
    Cohesion = 38e6 # pa
    TensileStr = 38e6 # pa

    rock = JCFpmMat(type=1,young=Young,frictionAngle=FrictAng,density=Density,poisson=Poisson,tensileStrength=TensileStr,cohesion=Cohesion,label='rock')
    proppant = JCFpmMat(type=2,young=Young_p,frictionAngle=FrictAng_p,density=Density_p,tensileStrength=TensileStr_p,cohesion=Cohesion_p,label='proppant')
    O.materials.append(JCFpmMat(young=Young_p,frictionAngle=0,density=0,label='wallmat'))

    O.materials.append(proppant)
    O.materials.append(rock)

    rock_assembly = O.bodies.append(ymport.textExt('new_rock.txt','x_y_z_r',shift=Vector3(0,0,0),scale=1,material=rock))

    for b in O.bodies:
      if b.state.pos[1]>0.017 or b.state.pos[1]<0.0135:
        O.bodies.erase(b.id)
      #if (b.state.pos[1]>0.0135 and b.state.pos[1]<0.015) :
       # b.groupMask=2
      else:
        b.dynamic = False
        b.state.vel=(0, 0, 0)

    proppant_assembly = O.bodies.append(ymport.textExt('new_prop.txt','x_y_z_r',shift=Vector3(0,0,0),scale=1,material=proppant,color=(1.00,0.67,0.50)))
    print ('length of bodies proppant '+str(len(O.bodies)))

    for b in proppant_assembly:
      O.bodies[b].groupMask=2
      if O.bodies[b].state.pos[0] < 0:
        print("found it ")
        O.bodies.erase(b)

    sphereIDs = [b.id for b in O.bodies if (type(b.shape)==Sphere and b.material.label=='proppant')]

    #coupling engine settings
    fluidCoupling.setNumParticles(len(sphereIDs))
    fluidCoupling.setIdList(sphereIDs)
    fluidCoupling.isGaussianInterp=False; #use pimpleFoamYade for gaussianInterp

    # Integrator
    newton=NewtonIntegrator(damping=0.0, gravity = (0.0 ,0.0, 0.0))
     # add small damping in case of stability issues.. ~ 0.1 max, also note : If gravity is needed, set it in constant/g dir.

    O.engines=[
  ForceResetter(),
  InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.2, label="bols"),Bo1_Facet_Aabb()], allowBiggerThanPeriod=True),
  InteractionLoop(
    [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.2, label="ig2s"),Ig2_Facet_Sphere_ScGeom()],
    [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1)],
        [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,label='lawFunctor')]
  ),
  GlobalStiffnessTimeStepper(timestepSafetyCoefficient=0.7, label = "ts"),
        PyRunner(command='sim.deletePar()',iterPeriod=30, label='checkPar'),
        fluidCoupling, #to be called after timestepper
        PyRunner(command='sim.printMessage()', iterPeriod= 1000, label='outputMessage'),
  newton,
        #PyRunner(command='sim.deletePar()',iterPeriod=50, label='checkPar'),
        VTKRecorder(fileName='yadep/3d-vtk-',mask = 2,recorders=['spheres','colors'],iterPeriod=1000)
    ]

  def printMessage(self):

     print("********************************YADE-ITER = " + str(O.iter) +" **********************************")
     if O.iter == 4000:
         maxVel = 0.05
         for b in O.bodies:
             if type(b.shape)==Sphere:
                 bodyVel = abs(b.state.vel.norm())
                 if bodyVel > maxVel:
                     raise ValueError("Body velocity exceeds imposed shear velocity by ", abs(bodyVel-maxVel))

  def deletePar(self):
    ids = fluidCoupling.getIdList()
    for b in ids:
      temp=0
      for i in O.bodies[b].intrs():
        if i.phys.isCohesive==True:
          temp+=1
      if temp==0:
        fluidCoupling.eraseId(b)
        print('deleted '+str(b))
    print("******YADE-ITER = " + str(O.iter) +" **********")

    #if O.bodies[b].state.nbInitBonds==0 or (O.bodies[b].state.nbInitBonds==O.bodies[b].state.nbBrokenBonds):

  def steprun(self):
    O.step()
    bols.aabbEnlargeFactor = 1
    ig2s.interactionDetectionFactor = 1
    print('-------------------reset aabbEnlargeFactor------------------')

  def irun(self,num):
      O.run(num,1)

if __name__=="__main__":
    sim = simulation()
    sim.steprun()
    sim.irun(5000) #run 5000 iteration and wait
    # print("body id = ", O.bodies[34].id)
    fluidCoupling.killMPI()

import builtins
builtins.sim=sim

Hello,

There could be several reasons for this :

Have you checked for the stability of the system ? Do particles 'fly' out
of the domain within a few timesteps ? Have you made sure the courant
number/time step in the fluid side is an acceptable value (<0.5) ?
If all these are okay, then I suspect it is with the deletePar function as
it may not be in sync with coupling. One way to ensure this is to add an
if condition to check whether coupling occurs in the present time step..

def deletePar(self) :
     if (O.iter%fluidCoupling.dataExchangeInterval==0) :
               ...
try it and let me know..

On Tue, Sep 3, 2019 at 10:33 AM Anqi H <email address hidden>
wrote:

> Question #683336 on Yade changed:
> https://answers.launchpad.net/yade/+question/683336
>
> Status: Answered => Open
>
> Anqi H is still having a problem:
> Hi Deepak,
>
> Thank you so much for your help so far. I've modified my script based on
> your comment to check if a body does not have any cohesive bonds in its
> interactions, if true then delete this sphere from fluidcoupling IDs. In
> the first openFoam time step, some spheres were deleted as expected,
> however, after the deletions and in the same foam timestep, it seems the
> program hangs after 500 yade iterations and in the fluidCoupling step (I
> have put the particle deletion pyrunner before calling fluidCoupling). I
> have copied my python script below. I am sorry it's probably difficult
> to run this script as it requires two other txt files that stored the
> sphere locations. I was wondering if it was due to some obvious errors
> in my script that I haven't noticed?
>
>
> from __future__ import print_function
> import sys
> from yadeimport import *
> from yade.utils import *
> from yade import ymport
>
> initMPI() #Initialize the mpi environment,
> always required.
> fluidCoupling = yade.FoamCoupling(); #Initialize the engine
> fluidCoupling.getRank(); #part of Initialization.
>
>
> #example of spheres in shear flow : two-way point force coupling
> class simulation():
>
> def __init__(self):
>
> O.periodic = True;
>
> #proppant properties
> FrictAng_p = 0.9
> Density_p = 2650
> Young_p = 100e6
> TensileStr_p=3000
> Cohesion_p=3000
>
> Young = 57e8 #apparent modulus
> FrictAng = 0.5
> Density = 2650
> Poisson = 0.28
> Cohesion = 38e6 # pa
> TensileStr = 38e6 # pa
>
> rock =
> JCFpmMat(type=1,young=Young,frictionAngle=FrictAng,density=Density,poisson=Poisson,tensileStrength=TensileStr,cohesion=Cohesion,label='rock')
> proppant =
> JCFpmMat(type=2,young=Young_p,frictionAngle=FrictAng_p,density=Density_p,tensileStrength=TensileStr_p,cohesion=Cohesion_p,label='proppant')
>
> O.materials.append(JCFpmMat(young=Young_p,frictionAngle=0,density=0,label='wallmat'))
>
> O.materials.append(proppant)
> O.materials.append(rock)
>
> rock_assembly =
>
> O.bodies.append(ymport.textExt('new_rock.txt','x_y_z_r',shift=Vector3(0,0,0),scale=1,material=rock))
>
> for b in O.bodies:
> if b.state.pos[1]>0.017 or b.state.pos[1]<0.0135:
> O.bodies.erase(b.id)
> #if (b.state.pos[1]>0.0135 and b.state.pos[1]<0.015) :
> # b.groupMask=2
> else:
> b.dynamic = False
> b.state.vel=(0, 0, 0)
>
> proppant_assembly =
> O.bodies.append(ymport.textExt('new_prop.txt','x_y_z_r',shift=Vector3(0,0,0),scale=1,material=proppant,color=(1.00,0.67,0.50)))
> print ('length of bodies proppant '+str(len(O.bodies)))
>
> for b in proppant_assembly:
> O.bodies[b].groupMask=2
> if O.bodies[b].state.pos[0] < 0:
> print("found it ")
> O.bodies.erase(b)
>
>
> sphereIDs = [b.id for b in O.bodies if (type(b.shape)==Sphere and
> b.material.label=='proppant')]
>
> #coupling engine settings
> fluidCoupling.setNumParticles(len(sphereIDs))
> fluidCoupling.setIdList(sphereIDs)
> fluidCoupling.isGaussianInterp=False; #use pimpleFoamYade for
> gaussianInterp
>
> # Integrator
> newton=NewtonIntegrator(damping=0.0, gravity = (0.0 ,0.0, 0.0))
> # add small damping in case of stability issues.. ~ 0.1 max, also
> note : If gravity is needed, set it in constant/g dir.
>
>
> O.engines=[
> ForceResetter(),
> InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.2,
> label="bols"),Bo1_Facet_Aabb()], allowBiggerThanPeriod=True),
> InteractionLoop(
> [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.2,
> label="ig2s"),Ig2_Facet_Sphere_ScGeom()],
> [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1)],
>
> [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,label='lawFunctor')]
> ),
> GlobalStiffnessTimeStepper(timestepSafetyCoefficient=0.7, label = "ts"),
> PyRunner(command='sim.deletePar()',iterPeriod=30,
> label='checkPar'),
> fluidCoupling, #to be called after timestepper
> PyRunner(command='sim.printMessage()', iterPeriod= 1000,
> label='outputMessage'),
> newton,
> #PyRunner(command='sim.deletePar()',iterPeriod=50,
> label='checkPar'),
> VTKRecorder(fileName='yadep/3d-vtk-',mask =
> 2,recorders=['spheres','colors'],iterPeriod=1000)
> ]
>
>
> def printMessage(self):
>
> print("********************************YADE-ITER = " + str(O.iter) +"
> **********************************")
> if O.iter == 4000:
> maxVel = 0.05
> for b in O.bodies:
> if type(b.shape)==Sphere:
> bodyVel = abs(b.state.vel.norm())
> if bodyVel > maxVel:
> raise ValueError("Body velocity exceeds imposed shear
> velocity by ", abs(bodyVel-maxVel))
>
> def deletePar(self):
> ids = fluidCoupling.getIdList()
> for b in ids:
> temp=0
> for i in O.bodies[b].intrs():
> if i.phys.isCohesive==True:
> temp+=1
> if temp==0:
> fluidCoupling.eraseId(b)
> print('deleted '+str(b))
> print("******YADE-ITER = " + str(O.iter) +" **********")
>
> #if O.bodies[b].state.nbInitBonds==0 or
> (O.bodies[b].state.nbInitBonds==O.bodies[b].state.nbBrokenBonds):
>
>
> def steprun(self):
> O.step()
> bols.aabbEnlargeFactor = 1
> ig2s.interactionDetectionFactor = 1
> print('-------------------reset aabbEnlargeFactor------------------')
>
> def irun(self,num):
> O.run(num,1)
>
>
> if __name__=="__main__":
> sim = simulation()
> sim.steprun()
> sim.irun(5000) #run 5000 iteration and wait
> # print("body id = ", O.bodies[34].id)
> fluidCoupling.killMPI()
>
> import builtins
> builtins.sim=sim
>
> --
> 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 : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

and also in the PyRunner for deletePar, set the iterPeriod=1.

On Wed, Sep 4, 2019 at 11:19 AM Deepak Kn <email address hidden> wrote:

> Hello,
>
> There could be several reasons for this :
>
> Have you checked for the stability of the system ? Do particles 'fly' out
> of the domain within a few timesteps ? Have you made sure the courant
> number/time step in the fluid side is an acceptable value (<0.5) ?
> If all these are okay, then I suspect it is with the deletePar function as
> it may not be in sync with coupling. One way to ensure this is to add an
> if condition to check whether coupling occurs in the present time step..
>
> def deletePar(self) :
> if (O.iter%fluidCoupling.dataExchangeInterval==0) :
> ...
> try it and let me know..
>
>
>
>
>
>
>
>
> On Tue, Sep 3, 2019 at 10:33 AM Anqi H <
> <email address hidden>> wrote:
>
>> Question #683336 on Yade changed:
>> https://answers.launchpad.net/yade/+question/683336
>>
>> Status: Answered => Open
>>
>> Anqi H is still having a problem:
>> Hi Deepak,
>>
>> Thank you so much for your help so far. I've modified my script based on
>> your comment to check if a body does not have any cohesive bonds in its
>> interactions, if true then delete this sphere from fluidcoupling IDs. In
>> the first openFoam time step, some spheres were deleted as expected,
>> however, after the deletions and in the same foam timestep, it seems the
>> program hangs after 500 yade iterations and in the fluidCoupling step (I
>> have put the particle deletion pyrunner before calling fluidCoupling). I
>> have copied my python script below. I am sorry it's probably difficult
>> to run this script as it requires two other txt files that stored the
>> sphere locations. I was wondering if it was due to some obvious errors
>> in my script that I haven't noticed?
>>
>>
>> from __future__ import print_function
>> import sys
>> from yadeimport import *
>> from yade.utils import *
>> from yade import ymport
>>
>> initMPI() #Initialize the mpi environment,
>> always required.
>> fluidCoupling = yade.FoamCoupling(); #Initialize the engine
>> fluidCoupling.getRank(); #part of Initialization.
>>
>>
>> #example of spheres in shear flow : two-way point force coupling
>> class simulation():
>>
>> def __init__(self):
>>
>> O.periodic = True;
>>
>> #proppant properties
>> FrictAng_p = 0.9
>> Density_p = 2650
>> Young_p = 100e6
>> TensileStr_p=3000
>> Cohesion_p=3000
>>
>> Young = 57e8 #apparent modulus
>> FrictAng = 0.5
>> Density = 2650
>> Poisson = 0.28
>> Cohesion = 38e6 # pa
>> TensileStr = 38e6 # pa
>>
>> rock =
>> JCFpmMat(type=1,young=Young,frictionAngle=FrictAng,density=Density,poisson=Poisson,tensileStrength=TensileStr,cohesion=Cohesion,label='rock')
>> proppant =
>> JCFpmMat(type=2,young=Young_p,frictionAngle=FrictAng_p,density=Density_p,tensileStrength=TensileStr_p,cohesion=Cohesion_p,label='proppant')
>>
>> O.materials.append(JCFpmMat(young=Young_p,frictionAngle=0,density=0,label='wallmat'))
>>
>> O.materials.append(proppant)
>> O.materials.append(rock)
>>
>> rock_assembly =
>>
>> O.bodies.append(ymport.textExt('new_rock.txt','x_y_z_r',shift=Vector3(0,0,0),scale=1,material=rock))
>>
>> for b in O.bodies:
>> if b.state.pos[1]>0.017 or b.state.pos[1]<0.0135:
>> O.bodies.erase(b.id)
>> #if (b.state.pos[1]>0.0135 and b.state.pos[1]<0.015) :
>> # b.groupMask=2
>> else:
>> b.dynamic = False
>> b.state.vel=(0, 0, 0)
>>
>> proppant_assembly =
>> O.bodies.append(ymport.textExt('new_prop.txt','x_y_z_r',shift=Vector3(0,0,0),scale=1,material=proppant,color=(1.00,0.67,0.50)))
>> print ('length of bodies proppant '+str(len(O.bodies)))
>>
>> for b in proppant_assembly:
>> O.bodies[b].groupMask=2
>> if O.bodies[b].state.pos[0] < 0:
>> print("found it ")
>> O.bodies.erase(b)
>>
>>
>> sphereIDs = [b.id for b in O.bodies if (type(b.shape)==Sphere and
>> b.material.label=='proppant')]
>>
>> #coupling engine settings
>> fluidCoupling.setNumParticles(len(sphereIDs))
>> fluidCoupling.setIdList(sphereIDs)
>> fluidCoupling.isGaussianInterp=False; #use pimpleFoamYade for
>> gaussianInterp
>>
>> # Integrator
>> newton=NewtonIntegrator(damping=0.0, gravity = (0.0 ,0.0, 0.0))
>> # add small damping in case of stability issues.. ~ 0.1 max, also
>> note : If gravity is needed, set it in constant/g dir.
>>
>>
>> O.engines=[
>> ForceResetter(),
>> InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.2,
>> label="bols"),Bo1_Facet_Aabb()], allowBiggerThanPeriod=True),
>> InteractionLoop(
>> [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.2,
>> label="ig2s"),Ig2_Facet_Sphere_ScGeom()],
>> [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1)],
>>
>> [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,label='lawFunctor')]
>> ),
>> GlobalStiffnessTimeStepper(timestepSafetyCoefficient=0.7, label = "ts"),
>> PyRunner(command='sim.deletePar()',iterPeriod=30,
>> label='checkPar'),
>> fluidCoupling, #to be called after timestepper
>> PyRunner(command='sim.printMessage()', iterPeriod= 1000,
>> label='outputMessage'),
>> newton,
>> #PyRunner(command='sim.deletePar()',iterPeriod=50,
>> label='checkPar'),
>> VTKRecorder(fileName='yadep/3d-vtk-',mask =
>> 2,recorders=['spheres','colors'],iterPeriod=1000)
>> ]
>>
>>
>> def printMessage(self):
>>
>> print("********************************YADE-ITER = " + str(O.iter)
>> +" **********************************")
>> if O.iter == 4000:
>> maxVel = 0.05
>> for b in O.bodies:
>> if type(b.shape)==Sphere:
>> bodyVel = abs(b.state.vel.norm())
>> if bodyVel > maxVel:
>> raise ValueError("Body velocity exceeds imposed
>> shear velocity by ", abs(bodyVel-maxVel))
>>
>> def deletePar(self):
>> ids = fluidCoupling.getIdList()
>> for b in ids:
>> temp=0
>> for i in O.bodies[b].intrs():
>> if i.phys.isCohesive==True:
>> temp+=1
>> if temp==0:
>> fluidCoupling.eraseId(b)
>> print('deleted '+str(b))
>> print("******YADE-ITER = " + str(O.iter) +" **********")
>>
>> #if O.bodies[b].state.nbInitBonds==0 or
>> (O.bodies[b].state.nbInitBonds==O.bodies[b].state.nbBrokenBonds):
>>
>>
>> def steprun(self):
>> O.step()
>> bols.aabbEnlargeFactor = 1
>> ig2s.interactionDetectionFactor = 1
>> print('-------------------reset aabbEnlargeFactor------------------')
>>
>> def irun(self,num):
>> O.run(num,1)
>>
>>
>> if __name__=="__main__":
>> sim = simulation()
>> sim.steprun()
>> sim.irun(5000) #run 5000 iteration and wait
>> # print("body id = ", O.bodies[34].id)
>> fluidCoupling.killMPI()
>>
>> import builtins
>> builtins.sim=sim
>>
>> --
>> 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 : <email address hidden>
>> Unsubscribe : https://launchpad.net/~yade-users
>> More help : https://help.launchpad.net/ListHelp
>>
>

Hi Deepak, you are right. The system is not stable as the courant number increases to a large number very fast, it may be due to my fluid domain was set up too short. Some particles fly out that hang the simulation as well, for which I'm not too sure if it was because of icoFoam being a transient solver, and the particles were placed very close to the inlet patch.

I have another question, as I'm thinking of using pimpleFoam solver for the same problem (laminar flow through a square pipe), I was wondering if I disable paticle motion by setting b.dynamics = false and b.state.vel = (0,0,0) but still include them in fluidCoupling ID list, would this setup be suitable for studying the fluid flow in the porous media? Thank you.

Hi Deepka,

I've decided to stick with icoFoam as my fluid flow is laminar and they are both transient solvers. I've made my fluid domain a longer section so icoFoam solver can now converge with a small Courant number at the end (0.115). I'm still trying to disable the particle motion to just measure the fluid conductivity after adding the particles, however the program hangs after 1 time step. I am not too sure if disabling the particle motion (b.dynamic = False) but including the particles in the fluidCoupling ID list is not allowed by the coupling engine?

The part of the script where I disabled particle motions is

    proppant_assembly = O.bodies.append(ymport.textExt('new_prop.txt','x_y_z_r',shift=Vector3(0,0,0),scale=1,material=proppant,color=(1.00,0.67,0.50)))
    print ('length of bodies proppant '+str(len(O.bodies)))

    for b in proppant_assembly:
      O.bodies[b].groupMask=2
      O.bodies[b].dynamic = False
      O.bodies[b].state.vel=(0, 0, 0)
      if O.bodies[b].state.pos[0] < 0:
        print("found it ")
        O.bodies.erase(b)

    sphereIDs = [b.id for b in O.bodies if (type(b.shape)==Sphere and b.material.label=='proppant')]

Hi,

>I've decided to stick with icoFoam as my fluid flow is laminar and they are both transient solvers

There is more to than just the problem being laminar and transient, the choice of the solver also depends on the volume fraction of the solid particle and the Reynolds number at the particle scale ...

you need to set the ids to fluidCoupling even if the particles are static, so in your script after the sphereIDs you need to add :

fluidCoupling.setIdList(sphereIDs)

Hi Deepak,

I'm running the same simulation with icoFoam for 0.05s, the particles are defined in fluidCouplingIDList and dynamic is disabled. The fluid case without the particles converges in openFoam, but after the particles are included, the simulation hangs after 0.01275s. I've been having trouble understanding why the program hangs here as the courant number is smaller than 0.5. The terminal output is copied to this post. Another question is about inspecting p,U properties of the fluid cells, since they are saved under two directories processor0,1, does it mean paraFoam cannot be used in this case? Thank you.

Time = 0.01275

Courant Number mean: 0.142852 max: 3.7703
smoothSolver: Solving for Ux, Initial residual = 0.560774, Final residual = 7.25301e-06, No Iterations 13
smoothSolver: Solving for Uy, Initial residual = 0.210176, Final residual = 8.97415e-06, No Iterations 10
smoothSolver: Solving for Uz, Initial residual = 0.580007, Final residual = 8.33788e-06, No Iterations 13
DICPCG: Solving for p, Initial residual = 0.87801, Final residual = 0.0268803, No Iterations 11
DICPCG: Solving for p, Initial residual = 0.0229963, Final residual = 0.00113944, No Iterations 9
DICPCG: Solving for p, Initial residual = 0.00111865, Final residual = 4.86087e-05, No Iterations 11
time step continuity errors : sum local = 2.63805e-06, global = -2.6854e-07, cumulative = 2.84037e-06
DICPCG: Solving for p, Initial residual = 0.216622, Final residual = 0.00789831, No Iterations 10
DICPCG: Solving for p, Initial residual = 0.00936633, Final residual = 0.000371535, No Iterations 11
DICPCG: Solving for p, Initial residual = 0.000378234, Final residual = 9.57322e-07, No Iterations 17
time step continuity errors : sum local = 3.7861e-08, global = 3.12631e-09, cumulative = 2.8435e-06
ExecutionTime = 99.51 s ClockTime = 100 s

Hi,
I need to have a look at the case, boundary conditions and things to
identify why this happens (my best bet is bad physics/numerical settings),
paraFoam can be used to view parallel cases, did you have a look into the
openfoam documentation/coupling documentation? The methods to post process
parallel
cases are there ..

https://cfd.direct/openfoam/user-guide/v6-paraview/
https://yade-dem.org/doc/FoamCoupling.html#post-processing

On Fri, Sep 13, 2019 at 11:09 AM Anqi H <
<email address hidden>> wrote:

> Question #683336 on Yade changed:
> https://answers.launchpad.net/yade/+question/683336
>
> Status: Answered => Open
>
> Anqi H is still having a problem:
> Hi Deepak,
>
> I'm running the same simulation with icoFoam for 0.05s, the particles
> are defined in fluidCouplingIDList and dynamic is disabled. The fluid
> case without the particles converges in openFoam, but after the
> particles are included, the simulation hangs after 0.01275s. I've been
> having trouble understanding why the program hangs here as the courant
> number is smaller than 0.5. The terminal output is copied to this post.
> Another question is about inspecting p,U properties of the fluid cells,
> since they are saved under two directories processor0,1, does it mean
> paraFoam cannot be used in this case? Thank you.
>
> Time = 0.01275
>
> Courant Number mean: 0.142852 max: 3.7703
> smoothSolver: Solving for Ux, Initial residual = 0.560774, Final residual
> = 7.25301e-06, No Iterations 13
> smoothSolver: Solving for Uy, Initial residual = 0.210176, Final residual
> = 8.97415e-06, No Iterations 10
> smoothSolver: Solving for Uz, Initial residual = 0.580007, Final residual
> = 8.33788e-06, No Iterations 13
> DICPCG: Solving for p, Initial residual = 0.87801, Final residual =
> 0.0268803, No Iterations 11
> DICPCG: Solving for p, Initial residual = 0.0229963, Final residual =
> 0.00113944, No Iterations 9
> DICPCG: Solving for p, Initial residual = 0.00111865, Final residual =
> 4.86087e-05, No Iterations 11
> time step continuity errors : sum local = 2.63805e-06, global =
> -2.6854e-07, cumulative = 2.84037e-06
> DICPCG: Solving for p, Initial residual = 0.216622, Final residual =
> 0.00789831, No Iterations 10
> DICPCG: Solving for p, Initial residual = 0.00936633, Final residual =
> 0.000371535, No Iterations 11
> DICPCG: Solving for p, Initial residual = 0.000378234, Final residual =
> 9.57322e-07, No Iterations 17
> time step continuity errors : sum local = 3.7861e-08, global =
> 3.12631e-09, cumulative = 2.8435e-06
> ExecutionTime = 99.51 s ClockTime = 100 s
>
> --
> 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 : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

 I've noticed that the max Courant No. is 3.7703, that doesn't look
good...have you tried decreasing the time step in Yade?

On Fri, Sep 13, 2019 at 11:36 AM Deepak Kn <email address hidden> wrote:

> Hi,
> I need to have a look at the case, boundary conditions and things to
> identify why this happens (my best bet is bad physics/numerical settings),
> paraFoam can be used to view parallel cases, did you have a look into the
> openfoam documentation/coupling documentation? The methods to post process
> parallel
> cases are there ..
>
> https://cfd.direct/openfoam/user-guide/v6-paraview/
> https://yade-dem.org/doc/FoamCoupling.html#post-processing
>
> On Fri, Sep 13, 2019 at 11:09 AM Anqi H <
> <email address hidden>> wrote:
>
>> Question #683336 on Yade changed:
>> https://answers.launchpad.net/yade/+question/683336
>>
>> Status: Answered => Open
>>
>> Anqi H is still having a problem:
>> Hi Deepak,
>>
>> I'm running the same simulation with icoFoam for 0.05s, the particles
>> are defined in fluidCouplingIDList and dynamic is disabled. The fluid
>> case without the particles converges in openFoam, but after the
>> particles are included, the simulation hangs after 0.01275s. I've been
>> having trouble understanding why the program hangs here as the courant
>> number is smaller than 0.5. The terminal output is copied to this post.
>> Another question is about inspecting p,U properties of the fluid cells,
>> since they are saved under two directories processor0,1, does it mean
>> paraFoam cannot be used in this case? Thank you.
>>
>> Time = 0.01275
>>
>> Courant Number mean: 0.142852 max: 3.7703
>> smoothSolver: Solving for Ux, Initial residual = 0.560774, Final
>> residual = 7.25301e-06, No Iterations 13
>> smoothSolver: Solving for Uy, Initial residual = 0.210176, Final
>> residual = 8.97415e-06, No Iterations 10
>> smoothSolver: Solving for Uz, Initial residual = 0.580007, Final
>> residual = 8.33788e-06, No Iterations 13
>> DICPCG: Solving for p, Initial residual = 0.87801, Final residual =
>> 0.0268803, No Iterations 11
>> DICPCG: Solving for p, Initial residual = 0.0229963, Final residual =
>> 0.00113944, No Iterations 9
>> DICPCG: Solving for p, Initial residual = 0.00111865, Final residual =
>> 4.86087e-05, No Iterations 11
>> time step continuity errors : sum local = 2.63805e-06, global =
>> -2.6854e-07, cumulative = 2.84037e-06
>> DICPCG: Solving for p, Initial residual = 0.216622, Final residual =
>> 0.00789831, No Iterations 10
>> DICPCG: Solving for p, Initial residual = 0.00936633, Final residual =
>> 0.000371535, No Iterations 11
>> DICPCG: Solving for p, Initial residual = 0.000378234, Final residual =
>> 9.57322e-07, No Iterations 17
>> time step continuity errors : sum local = 3.7861e-08, global =
>> 3.12631e-09, cumulative = 2.8435e-06
>> ExecutionTime = 99.51 s ClockTime = 100 s
>>
>> --
>> 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 : <email address hidden>
>> Unsubscribe : https://launchpad.net/~yade-users
>> More help : https://help.launchpad.net/ListHelp
>>
>

Hi Deepak,

Thank you for the links. I've made deltaT 1/10 of the previous value, It now stalls at 0.01468s while both courant and max courant number is less than 0.5.

Time = 0.01478

Courant Number mean: 0.00188364 max: 0.00654247
smoothSolver: Solving for Ux, Initial residual = 0.000895276, Final residual = 3.99666e-07, No Iterations 1
smoothSolver: Solving for Uy, Initial residual = 0.0012464, Final residual = 5.33346e-07, No Iterations 1
smoothSolver: Solving for Uz, Initial residual = 0.00136997, Final residual = 3.99407e-07, No Iterations 1
DICPCG: Solving for p, Initial residual = 0.000356892, Final residual = 1.69725e-05, No Iterations 9
DICPCG: Solving for p, Initial residual = 1.69772e-05, Final residual = 7.49094e-07, No Iterations 10
DICPCG: Solving for p, Initial residual = 7.49094e-07, Final residual = 7.49094e-07, No Iterations 0
time step continuity errors : sum local = 6.37045e-12, global = -2.646e-14, cumulative = 1.25224e-08
DICPCG: Solving for p, Initial residual = 2.68288e-06, Final residual = 6.77363e-07, No Iterations 2
DICPCG: Solving for p, Initial residual = 6.77362e-07, Final residual = 6.77362e-07, No Iterations 0
DICPCG: Solving for p, Initial residual = 6.77362e-07, Final residual = 6.77362e-07, No Iterations 0
time step continuity errors : sum local = 5.76044e-12, global = -8.4556e-13, cumulative = 1.25216e-08
ExecutionTime = 128.46 s ClockTime = 129 s

ControlDict:
/*--------------------------------*- C++ -*----------------------------------*\
  ========= |
  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
   \\ / O peration | Website: https://openfoam.org
    \\ / A nd | Version: dev
     \\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version 2.0;
    format ascii;
    class dictionary;
    location "system";
    object controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application icoFoam;

startFrom latestTime;

startTime 0;

stopAt endTime;

endTime 0.05;

deltaT 0.000005;

writeControl timeStep;

writeInterval 20;

purgeWrite 0;

writeFormat ascii;

writePrecision 6;

writeCompression off;

timeFormat general;

timePrecision 6;

runTimeModifiable true;

// ************************************************************************* //

fv solution:
/*--------------------------------*- C++ -*----------------------------------*\
  ========= |
  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
   \\ / O peration | Website: https://openfoam.org
    \\ / A nd | Version: dev
     \\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version 2.0;
    format ascii;
    class dictionary;
    location "system";
    object fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    p
    {
        solver PCG;
        preconditioner DIC;
        tolerance 1e-06;
        relTol 0.05;
    }

    pFinal
    {
        $p;
        relTol 0;
    }

    U
    {
        solver smoothSolver;
        smoother symGaussSeidel;
        tolerance 1e-05;
        relTol 0;
    }
}

PISO
{
    nCorrectors 2;
    nNonOrthogonalCorrectors 2;
}

// ************************************************************************* //

fv scheme:
/*--------------------------------*- C++ -*----------------------------------*\
  ========= |
  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
   \\ / O peration | Website: https://openfoam.org
    \\ / A nd | Version: dev
     \\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version 2.0;
    format ascii;
    class dictionary;
    location "system";
    object fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default Euler;
}

gradSchemes
{
    default Gauss linear;
}

divSchemes
{
    default none;
    div(phi,U) Gauss limitedLinearV 1;
}

laplacianSchemes
{
    default Gauss linear corrected;
}

interpolationSchemes
{
    default linear;
}

snGradSchemes
{
    default corrected;
}

// ************************************************************************* //

0/p

/*--------------------------------*- C++ -*----------------------------------*\
  ========= |
  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
   \\ / O peration | Website: https://openfoam.org
    \\ / A nd | Version: dev
     \\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version 2.0;
    format ascii;
    class volScalarField;
    object p;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions [0 2 -2 0 0 0 0];

internalField uniform 0;

boundaryField
{
    inlet
    {
        type fixedValue;
     value uniform .2;
    }

    outlet
    {
        type fixedValue;
        value uniform 0;
    }

    walls
    {
        type zeroGradient;
    }

}

// ************************************************************************* //

0/U
/*--------------------------------*- C++ -*----------------------------------*\
  ========= |
  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
   \\ / O peration | Website: https://openfoam.org
    \\ / A nd | Version: dev
     \\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version 2.0;
    format ascii;
    class volVectorField;
    object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions [0 1 -1 0 0 0 0];

internalField uniform (1 0 0);

boundaryField
{
    inlet
    {
        type zeroGradient;
    }

    outlet
    {
        type zeroGradient;
    }

    walls
    {
        type fixedValue;
 value uniform (0 0 0);
    }
}

// ************************************************************************* //

0/uSource
/*--------------------------------*- C++ -*----------------------------------*\
  ========= |
  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
   \\ / O peration | Website: https://openfoam.org
    \\ / A nd | Version: 6
     \\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version 2.0;
    format ascii;
    class volVectorField;
    object uSource;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions [0 1 -2 0 0 0 0];

internalField uniform (0 0 0);

boundaryField
{

    inlet
    {
        type calculated;
        value uniform (0 0 0);
    }

    outlet
    {
        type calculated;
        value uniform (0 0 0);
    }

    walls
    {
        type fixedValue;
        value uniform (0 0 0);
    }

}

// ************************************************************************* //

blockMeshDict:
/*--------------------------------*- C++ -*----------------------------------*\
  ========= |
  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
   \\ / O peration | Website: https://openfoam.org
    \\ / A nd | Version: dev
     \\/ M anipulation |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version 2.0;
    format ascii;
    class dictionary;
    object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

convertToMeters 0.001;

vertices
(
    (-2 10 -2)
    (22 10 -2)
    (22 22 -2)
    (-2 22 -2)
    (-2 10 22)
    (22 10 22)
    (22 22 22)
    (-2 22 22)
);

blocks
(
    hex (0 1 2 3 4 5 6 7) (20 4 20) simpleGrading (1 1 1)
);

edges
(
);

boundary
(
    inlet
    {
        type patch;
        faces
        (
            (0 3 7 4)
        );
    }
    outlet
    {
        type patch;
        faces
        (
            (1 2 6 5)
        );
    }
    walls
    {
        type wall;
        faces
        (
            (0 1 2 3)
        (0 1 5 4)
        (4 5 6 7)
        (3 2 6 7)
        );
    }
);

// ************************************************************************* //

hi,
can you send the yade script as well?

On Fri, Sep 13, 2019 at 3:57 PM Anqi H <email address hidden>
wrote:

> Question #683336 on Yade changed:
> https://answers.launchpad.net/yade/+question/683336
>
> Anqi H posted a new comment:
> Hi Deepak,
>
> Thank you for the links. I've made deltaT 1/10 of the previous value, It
> now stalls at 0.01468s while both courant and max courant number is less
> than 0.5.
>
> Time = 0.01478
>
> Courant Number mean: 0.00188364 max: 0.00654247
> smoothSolver: Solving for Ux, Initial residual = 0.000895276, Final
> residual = 3.99666e-07, No Iterations 1
> smoothSolver: Solving for Uy, Initial residual = 0.0012464, Final
> residual = 5.33346e-07, No Iterations 1
> smoothSolver: Solving for Uz, Initial residual = 0.00136997, Final
> residual = 3.99407e-07, No Iterations 1
> DICPCG: Solving for p, Initial residual = 0.000356892, Final residual =
> 1.69725e-05, No Iterations 9
> DICPCG: Solving for p, Initial residual = 1.69772e-05, Final residual =
> 7.49094e-07, No Iterations 10
> DICPCG: Solving for p, Initial residual = 7.49094e-07, Final residual =
> 7.49094e-07, No Iterations 0
> time step continuity errors : sum local = 6.37045e-12, global =
> -2.646e-14, cumulative = 1.25224e-08
> DICPCG: Solving for p, Initial residual = 2.68288e-06, Final residual =
> 6.77363e-07, No Iterations 2
> DICPCG: Solving for p, Initial residual = 6.77362e-07, Final residual =
> 6.77362e-07, No Iterations 0
> DICPCG: Solving for p, Initial residual = 6.77362e-07, Final residual =
> 6.77362e-07, No Iterations 0
> time step continuity errors : sum local = 5.76044e-12, global =
> -8.4556e-13, cumulative = 1.25216e-08
> ExecutionTime = 128.46 s ClockTime = 129 s
>
> ControlDict:
> /*--------------------------------*- C++
> -*----------------------------------*\
> ========= |
> \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
> \\ / O peration | Website: https://openfoam.org
> \\ / A nd | Version: dev
> \\/ M anipulation |
>
> \*---------------------------------------------------------------------------*/
> FoamFile
> {
> version 2.0;
> format ascii;
> class dictionary;
> location "system";
> object controlDict;
> }
> // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
> * //
>
> application icoFoam;
>
> startFrom latestTime;
>
> startTime 0;
>
> stopAt endTime;
>
> endTime 0.05;
>
> deltaT 0.000005;
>
> writeControl timeStep;
>
> writeInterval 20;
>
> purgeWrite 0;
>
> writeFormat ascii;
>
> writePrecision 6;
>
> writeCompression off;
>
> timeFormat general;
>
> timePrecision 6;
>
> runTimeModifiable true;
>
>
> //
> ************************************************************************* //
>
> fv solution:
> /*--------------------------------*- C++
> -*----------------------------------*\
> ========= |
> \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
> \\ / O peration | Website: https://openfoam.org
> \\ / A nd | Version: dev
> \\/ M anipulation |
>
> \*---------------------------------------------------------------------------*/
> FoamFile
> {
> version 2.0;
> format ascii;
> class dictionary;
> location "system";
> object fvSolution;
> }
> // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
> * //
>
> solvers
> {
> p
> {
> solver PCG;
> preconditioner DIC;
> tolerance 1e-06;
> relTol 0.05;
> }
>
> pFinal
> {
> $p;
> relTol 0;
> }
>
> U
> {
> solver smoothSolver;
> smoother symGaussSeidel;
> tolerance 1e-05;
> relTol 0;
> }
> }
>
> PISO
> {
> nCorrectors 2;
> nNonOrthogonalCorrectors 2;
> }
>
>
> //
> ************************************************************************* //
>
> fv scheme:
> /*--------------------------------*- C++
> -*----------------------------------*\
> ========= |
> \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
> \\ / O peration | Website: https://openfoam.org
> \\ / A nd | Version: dev
> \\/ M anipulation |
>
> \*---------------------------------------------------------------------------*/
> FoamFile
> {
> version 2.0;
> format ascii;
> class dictionary;
> location "system";
> object fvSchemes;
> }
> // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
> * //
>
> ddtSchemes
> {
> default Euler;
> }
>
> gradSchemes
> {
> default Gauss linear;
> }
>
> divSchemes
> {
> default none;
> div(phi,U) Gauss limitedLinearV 1;
> }
>
> laplacianSchemes
> {
> default Gauss linear corrected;
> }
>
> interpolationSchemes
> {
> default linear;
> }
>
> snGradSchemes
> {
> default corrected;
> }
>
>
> //
> ************************************************************************* //
>
> 0/p
>
> /*--------------------------------*- C++
> -*----------------------------------*\
> ========= |
> \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
> \\ / O peration | Website: https://openfoam.org
> \\ / A nd | Version: dev
> \\/ M anipulation |
>
> \*---------------------------------------------------------------------------*/
> FoamFile
> {
> version 2.0;
> format ascii;
> class volScalarField;
> object p;
> }
> // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
> * //
>
> dimensions [0 2 -2 0 0 0 0];
>
> internalField uniform 0;
>
> boundaryField
> {
> inlet
> {
> type fixedValue;
> value uniform .2;
> }
>
> outlet
> {
> type fixedValue;
> value uniform 0;
> }
>
> walls
> {
> type zeroGradient;
> }
>
> }
>
> //
> *************************************************************************
> //
>
>
> 0/U
> /*--------------------------------*- C++
> -*----------------------------------*\
> ========= |
> \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
> \\ / O peration | Website: https://openfoam.org
> \\ / A nd | Version: dev
> \\/ M anipulation |
>
> \*---------------------------------------------------------------------------*/
> FoamFile
> {
> version 2.0;
> format ascii;
> class volVectorField;
> object U;
> }
> // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
> * //
>
> dimensions [0 1 -1 0 0 0 0];
>
> internalField uniform (1 0 0);
>
> boundaryField
> {
> inlet
> {
> type zeroGradient;
> }
>
> outlet
> {
> type zeroGradient;
> }
>
> walls
> {
> type fixedValue;
> value uniform (0 0 0);
> }
> }
>
> //
> *************************************************************************
> //
>
> 0/uSource
> /*--------------------------------*- C++
> -*----------------------------------*\
> ========= |
> \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
> \\ / O peration | Website: https://openfoam.org
> \\ / A nd | Version: 6
> \\/ M anipulation |
>
> \*---------------------------------------------------------------------------*/
> FoamFile
> {
> version 2.0;
> format ascii;
> class volVectorField;
> object uSource;
> }
> // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
> * //
>
> dimensions [0 1 -2 0 0 0 0];
>
> internalField uniform (0 0 0);
>
> boundaryField
> {
>
> inlet
> {
> type calculated;
> value uniform (0 0 0);
> }
>
> outlet
> {
> type calculated;
> value uniform (0 0 0);
> }
>
> walls
> {
> type fixedValue;
> value uniform (0 0 0);
> }
>
> }
>
> //
> *************************************************************************
> //
>
> blockMeshDict:
> /*--------------------------------*- C++
> -*----------------------------------*\
> ========= |
> \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
> \\ / O peration | Website: https://openfoam.org
> \\ / A nd | Version: dev
> \\/ M anipulation |
>
> \*---------------------------------------------------------------------------*/
> FoamFile
> {
> version 2.0;
> format ascii;
> class dictionary;
> object blockMeshDict;
> }
> // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
> * //
>
> convertToMeters 0.001;
>
> vertices
> (
> (-2 10 -2)
> (22 10 -2)
> (22 22 -2)
> (-2 22 -2)
> (-2 10 22)
> (22 10 22)
> (22 22 22)
> (-2 22 22)
> );
>
>
> blocks
> (
> hex (0 1 2 3 4 5 6 7) (20 4 20) simpleGrading (1 1 1)
> );
>
> edges
> (
> );
>
> boundary
> (
> inlet
> {
> type patch;
> faces
> (
> (0 3 7 4)
> );
> }
> outlet
> {
> type patch;
> faces
> (
> (1 2 6 5)
> );
> }
> walls
> {
> type wall;
> faces
> (
> (0 1 2 3)
> (0 1 5 4)
> (4 5 6 7)
> (3 2 6 7)
> );
> }
> );
>
> //
> *************************************************************************
> //
>
> --
> 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 : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

Hi Deepak,

I have attached the script and the txt file that has the particle positions and radius. The original txt file has about 3900 particles, I've tried running with 20 particles and the program hangs at around the same timestep. I also added the code to reset the interaction radius, now the program hangs at 0.027s still with small courant number. Can you please help me have a look? Thank you.

from __future__ import print_function
import sys
from yadeimport import *
from yade.utils import *
from yade import ymport

initMPI() #Initialize the mpi environment, always required.
fluidCoupling = yade.FoamCoupling(); #Initialize the engine
fluidCoupling.getRank(); #part of Initialization.

#example of spheres in shear flow : two-way point force coupling
class simulation():

  def __init__(self):

    O.periodic = True;

    #proppant properties
    FrictAng_p = 0.9
    Density_p = 2650
    Young_p = 100e6
    TensileStr_p=3000
    Cohesion_p=3000

    proppant = JCFpmMat(type=2,young=Young_p,frictionAngle=FrictAng_p,density=Density_p,tensileStrength=TensileStr_p,cohesion=Cohesion_p,label='proppant')

    O.materials.append(proppant)

    proppant_assembly = O.bodies.append(ymport.textExt('prop.txt','x_y_z_r',shift=Vector3(0,0,0),scale=1,material=proppant,color=(1.00,0.67,0.50)))
    print ('length of bodies proppant '+str(len(O.bodies)))

    for b in proppant_assembly:
      O.bodies[b].groupMask=2
      O.bodies[b].dynamic = False
      O.bodies[b].state.vel=(0, 0, 0)
      if O.bodies[b].state.pos[0] < 0:
        print("found it ")
        O.bodies.erase(b)

    sphereIDs = [b.id for b in O.bodies if (type(b.shape)==Sphere)]
    #coupling engine settings
    fluidCoupling.setNumParticles(len(sphereIDs))
    fluidCoupling.setIdList(sphereIDs)
    fluidCoupling.isGaussianInterp=False; #use pimpleFoamYade for gaussianInterp

    # Integrator
    newton=NewtonIntegrator(damping=0.0, gravity = (0.0 ,0.0, 0.0))
     # add small damping in case of stability issues.. ~ 0.1 max, also note : If gravity is needed, set it in constant/g dir.

    O.engines=[
  ForceResetter(),
  InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.2, label="bols"),Bo1_Facet_Aabb()], allowBiggerThanPeriod=True),
  InteractionLoop(
    [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.2, label="ig2s"),Ig2_Facet_Sphere_ScGeom()],
    [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1)],
        [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,label='lawFunctor')]
  ),
  GlobalStiffnessTimeStepper(timestepSafetyCoefficient=0.7, label = "ts"),
        fluidCoupling, #to be called after timestepper
        PyRunner(command='sim.printMessage()', iterPeriod= 1000, label='outputMessage'),
  newton,
        VTKRecorder(fileName='yadep/3d-vtk-',mask = 2,recorders=['spheres','colors'],iterPeriod=1000)
    ]

  def printMessage(self):
     print("********************************YADE-ITER = " + str(O.iter) +" **********************************")
     if O.iter == 4000:
         maxVel = 0.05
         for b in O.bodies:
             if type(b.shape)==Sphere:
                 bodyVel = abs(b.state.vel.norm())
                 if bodyVel > maxVel:
                     raise ValueError("Body velocity exceeds imposed shear velocity by ", abs(bodyVel-maxVel))

  def steprun(self):
    O.step()
    bols.aabbEnlargeFactor = 1
    ig2s.interactionDetectionFactor = 1
    print('-------------------reset aabbEnlargeFactor------------------')

  def irun(self,num):
      O.run(num,1)

if __name__=="__main__":
    sim = simulation()
    sim.steprun()
    sim.irun(5000) #run 5000 iteration and wait
    fluidCoupling.killMPI()

import builtins
builtins.sim=sim

prop.txt:
0.00852590934236 0.0153639385798 0.00176255575986 0.000217761439712
0.00764955868029 0.0152116159226 0.00436816482926 0.000217732009598
0.00560124234128 0.0151519211654 0.00842422207547 0.000217290557888
0.0063453810585 0.0151831422756 0.00480549526363 0.000216584235153
0.00960891441177 0.0149854671521 0.00879075377879 0.000216554805039
0.00483649588532 0.0153771873566 0.00930391144151 0.000216525374925
0.00842596919234 0.0153067625031 0.00235637778066 0.000216407654469
0.0040183447717 0.0152218324011 0.00854134696363 0.000216319364127
0.00730305871618 0.0150765919403 0.00722600431368 0.000216113353329
0.00601879303512 0.0150916069316 0.00798938563763 0.000215730761848
0.00963465709772 0.0149722137776 0.00959788700642 0.000215642471506
0.00678080489014 0.015448409339 0.00425480343164 0.000215407030594
0.00245574035206 0.015231337524 0.0095993867908 0.000215201019796
0.00647161027674 0.0151548077837 0.00635258123028 0.000215142159568
0.00809997911319 0.0149366515564 0.00924837631883 0.000214406406719
0.00225750713248 0.0153179180278 0.00995450815658 0.000214347546491
0.0082797736993 0.0151490590369 0.00890198988867 0.000213788374325
0.00776231509244 0.0151470998716 0.00874378437775 0.000213758944211
0.00905245969947 0.0151690471826 0.00509795120937 0.000213435212957
0.00240727047207 0.0152602638837 0.00837260770547 0.000213170341932
0.00397163696595 0.0152952188157 0.00621218690826 0.00021293490102

Hello,

The simulation hangs because Yade runs for 5000 timesteps and it exits :

if __name__=="__main__":
    sim = simulation()
    sim.steprun()
    sim.irun(5000) #run 5000 iteration and wait ---> this line.
    fluidCoupling.killMPI()

You need to set the same end time for the yade part as well (0.05), you can
do it like this in the irun function :
  def irun(self,time):
      while (O.time <= time):
            O.step()

sim.irun(0.05)

or you can set the O.dt and the corresponding number of iterations.

We do not have the interactivity of Yade in the coupling, i,e. to run for
few timesteps, pause and inspect and change stuff..

Note : even after the simulation has finished running , there is a chance
of hanging, if you see the message "Finalising parallel run" it means
everything is okay.

On Sat, Sep 14, 2019 at 7:17 AM Anqi H <email address hidden>
wrote:

> Question #683336 on Yade changed:
> https://answers.launchpad.net/yade/+question/683336
>
> Anqi H posted a new comment:
> Hi Deepak,
>
> I have attached the script and the txt file that has the particle
> positions and radius. The original txt file has about 3900 particles,
> I've tried running with 20 particles and the program hangs at around the
> same timestep. I also added the code to reset the interaction radius,
> now the program hangs at 0.027s still with small courant number. Can you
> please help me have a look? Thank you.
>
> from __future__ import print_function
> import sys
> from yadeimport import *
> from yade.utils import *
> from yade import ymport
>
> initMPI() #Initialize the mpi environment,
> always required.
> fluidCoupling = yade.FoamCoupling(); #Initialize the engine
> fluidCoupling.getRank(); #part of Initialization.
>
>
> #example of spheres in shear flow : two-way point force coupling
> class simulation():
>
> def __init__(self):
>
> O.periodic = True;
>
> #proppant properties
> FrictAng_p = 0.9
> Density_p = 2650
> Young_p = 100e6
> TensileStr_p=3000
> Cohesion_p=3000
>
> proppant =
> JCFpmMat(type=2,young=Young_p,frictionAngle=FrictAng_p,density=Density_p,tensileStrength=TensileStr_p,cohesion=Cohesion_p,label='proppant')
>
> O.materials.append(proppant)
>
> proppant_assembly =
> O.bodies.append(ymport.textExt('prop.txt','x_y_z_r',shift=Vector3(0,0,0),scale=1,material=proppant,color=(1.00,0.67,0.50)))
> print ('length of bodies proppant '+str(len(O.bodies)))
>
> for b in proppant_assembly:
> O.bodies[b].groupMask=2
> O.bodies[b].dynamic = False
> O.bodies[b].state.vel=(0, 0, 0)
> if O.bodies[b].state.pos[0] < 0:
> print("found it ")
> O.bodies.erase(b)
>
> sphereIDs = [b.id for b in O.bodies if (type(b.shape)==Sphere)]
> #coupling engine settings
> fluidCoupling.setNumParticles(len(sphereIDs))
> fluidCoupling.setIdList(sphereIDs)
> fluidCoupling.isGaussianInterp=False; #use pimpleFoamYade for
> gaussianInterp
>
> # Integrator
> newton=NewtonIntegrator(damping=0.0, gravity = (0.0 ,0.0, 0.0))
> # add small damping in case of stability issues.. ~ 0.1 max, also
> note : If gravity is needed, set it in constant/g dir.
>
>
> O.engines=[
> ForceResetter(),
> InsertionSortCollider([Bo1_Sphere_Aabb(aabbEnlargeFactor=1.2,
> label="bols"),Bo1_Facet_Aabb()], allowBiggerThanPeriod=True),
> InteractionLoop(
> [Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=1.2,
> label="ig2s"),Ig2_Facet_Sphere_ScGeom()],
> [Ip2_JCFpmMat_JCFpmMat_JCFpmPhys(cohesiveTresholdIteration=1)],
>
> [Law2_ScGeom_JCFpmPhys_JointedCohesiveFrictionalPM(recordCracks=True,label='lawFunctor')]
> ),
> GlobalStiffnessTimeStepper(timestepSafetyCoefficient=0.7, label = "ts"),
> fluidCoupling, #to be called after timestepper
> PyRunner(command='sim.printMessage()', iterPeriod= 1000,
> label='outputMessage'),
> newton,
> VTKRecorder(fileName='yadep/3d-vtk-',mask =
> 2,recorders=['spheres','colors'],iterPeriod=1000)
> ]
>
>
> def printMessage(self):
> print("********************************YADE-ITER = " + str(O.iter) +"
> **********************************")
> if O.iter == 4000:
> maxVel = 0.05
> for b in O.bodies:
> if type(b.shape)==Sphere:
> bodyVel = abs(b.state.vel.norm())
> if bodyVel > maxVel:
> raise ValueError("Body velocity exceeds imposed shear
> velocity by ", abs(bodyVel-maxVel))
>
>
> def steprun(self):
> O.step()
> bols.aabbEnlargeFactor = 1
> ig2s.interactionDetectionFactor = 1
> print('-------------------reset aabbEnlargeFactor------------------')
>
> def irun(self,num):
> O.run(num,1)
>
>
> if __name__=="__main__":
> sim = simulation()
> sim.steprun()
> sim.irun(5000) #run 5000 iteration and wait
> fluidCoupling.killMPI()
>
> import builtins
> builtins.sim=sim
>
> prop.txt:
> 0.00852590934236 0.0153639385798 0.00176255575986 0.000217761439712
> 0.00764955868029 0.0152116159226 0.00436816482926 0.000217732009598
> 0.00560124234128 0.0151519211654 0.00842422207547 0.000217290557888
> 0.0063453810585 0.0151831422756 0.00480549526363 0.000216584235153
> 0.00960891441177 0.0149854671521 0.00879075377879 0.000216554805039
> 0.00483649588532 0.0153771873566 0.00930391144151 0.000216525374925
> 0.00842596919234 0.0153067625031 0.00235637778066 0.000216407654469
> 0.0040183447717 0.0152218324011 0.00854134696363 0.000216319364127
> 0.00730305871618 0.0150765919403 0.00722600431368 0.000216113353329
> 0.00601879303512 0.0150916069316 0.00798938563763 0.000215730761848
> 0.00963465709772 0.0149722137776 0.00959788700642 0.000215642471506
> 0.00678080489014 0.015448409339 0.00425480343164 0.000215407030594
> 0.00245574035206 0.015231337524 0.0095993867908 0.000215201019796
> 0.00647161027674 0.0151548077837 0.00635258123028 0.000215142159568
> 0.00809997911319 0.0149366515564 0.00924837631883 0.000214406406719
> 0.00225750713248 0.0153179180278 0.00995450815658 0.000214347546491
> 0.0082797736993 0.0151490590369 0.00890198988867 0.000213788374325
> 0.00776231509244 0.0151470998716 0.00874378437775 0.000213758944211
> 0.00905245969947 0.0151690471826 0.00509795120937 0.000213435212957
> 0.00240727047207 0.0152602638837 0.00837260770547 0.000213170341932
> 0.00397163696595 0.0152952188157 0.00621218690826 0.00021293490102
>
> --
> 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 : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

That allows me alot. Also, I intend to enhance the triax.Py script (with facets) with the aid of practice the wall at the top and bottom as opposed to practice velocity for debris, wherein I ought to manage the pressure-->>deviatoric stress. Is that an excellent idea? You can see the data here https://plushiesshop.com/product-category/countryball-plush/