Many test cases are given to help users to understand how to use the different features offered by the MigFlow Software. Test cases are in the ~/MigFlow/testcases directory. They are given to provide users basic tools in order to build computation scripts for their own applications. Movies and figures of some applications are available on the MigFlow Project website
Test Cases Basic Structure
Let us illustrate the structure of the executable python files with a two dimensional deposit of grains in a fluid. This test case is provided in the testcases/depot-2d/depot.py file.
Some steps have to be followed in order to properly define the application you want to solve with the MigFlow Software.
- Set the problem parametres.
In an application, many variables are required to define :
g = np.array([0,-9.81,0]) # gravity
r = 1e-3 # grains radius
rhop = 1500 # grains density
rho = 1000 # fluid density
nu = 1e-6 # kinematic viscosity
outf = 5 # number of iterations between output files
dt = 1e-2 # time step
tEnd = 100 # final time
#geometrical parameters
ly = 5e-2 # grains area height
lx = 4e-1 # grains area widht
H = 1.2 # domain height
- Define the matrix of the grains positions and actually create the grains at these positions
- Setting the initial situation requires to define the initial positions of the grains. In test cases, this is often achieved by defining a function in which the positions of the centres are defined. The arguments of the function are related to the grains properties and the grains area geometry defined previously. The physical surfaces defined in your mesh.geo file are used to specify the solid boundaries for the grains. Grain objects are created locally and all the information about the initial conditions of the grains are written in an output file. The argument of the particles structure builder is the dimension of the problem.
def genInitialPosition(filename, r, H, ly, lx, rhop) :
"""Set all the particles centre positions and create the particles objects to add in the computing structure
Keyword arguments:
filename -- name of the output file
r -- max radius of the particles
H -- domain height
ly - particles area height
lx -- particles area width
rhop -- particles density
"""
# Particles structure builder
p = scontact.ParticleProblem(2)
# Load mesh.msh file specifying physical boundaries names
p.load_msh_boundaries("mesh.msh", ["Top", "Lateral","Bottom"])
#Definition of the points where the particles are located
x = np.arange(-lx/2+r, lx/2-r, 2*r)
y = np.arange(H/2-r, H/2-ly+r, -2*r)
x, y = np.meshgrid(x, y)
x = x.flat
y = y.flat
# Add a grain at each centre position
for i in range(len(x)) :
p.add_particle((x[i], y[i]), r, r**2 * np.pi * rhop)
p.write_vtk(filename,0,0)
- Use the previous function to set the initial condition and write it in the outputdir directory. Then, build the particle structure and fill the global particle structure with the information written in the initial file
genInitialPosition(outputdir, r, H, ly, lx, rhop)
p = scontact.ParticleProblem(2)
p.read_vtk(outputdir,0)
- Create the fluid structure and specify the computational domain (i.e. give the mesh to the fluid class). At this step, all the variables (velocity, pressure...) of the fluid problem and the properties of the fluid are stored in the fluid structure. The arguments for the fluid structure builder are the dimension of the problem, the gravity, the dynamic viscosity and the density of the fluid.
fluid = fluid.FluidProblem(2,g,[nu*rho],[rho])
# Set the mesh geometry for the fluid computation
fluid.load_msh("mesh.msh")
It is possible to solve simplified fluid equations using the temporal and advection parameters of the FluidProblem builder. If temporal is set to False, steady flows are computed. If advection is set to False, Stokes flows are computed. Default value for both temporal and advection is True to compute unsteady Navier-Stokes equations. 4. Once the fluid variables and properties are created, it is mandatory to set the boundary conditions of the fluid problem. Boundary conditions are of two types:
- Strong boundaries (Optional) constrain a field variable to a specified value. This is done by suppressing an equation of the linear system to replace it by this constraint :
fluid.set_strong_boundary(Physical boundary tag, field variable item, field value)
- Weak boundaries specify the flux at each boundary. The flux terms added by these conditions come from the integration by part arising from the finite element methods and can be suppressed for some weak boundary type. For example, if you consider a wall, it is quite obvious that the velocity flux across the boundary is zero but the pressure could be specified. There are two types of weak boundaries wall and open:
fluid.set_wall_boundary(tag, pressure=None)
"""Set the weak boundaries (=normal fluxes) for the fluid problem.
Keyword arguments:
tag -- physical tag (or list of tags), set in the mesh.geo file, of the wall boundaries
pressure -- the pressure value if imposed (callback or number)
"""
fluid.set_open_boundary(tag, velocity=None, pressure=None, porosity=1, concentration=1)
"""Set the weak boundaries (=normal fluxes) for the fluid problem.
Keyword arguments:
tag -- physical tag (or list of tags), set in the mesh.geo file, of the wall boundaries
velocity -- the velocity vector if imposed (callback or number))
pressure -- the pressure value if imposed (callback or number)
porosity -- the volume fraction of fluid at aperture (work in progress - should be one)
concentration -- the volume fraction of the first continuous phase
Raises:
NameError -- If the specified porosity outside the boundary is too small
NameError -- If velocity and pressure are not specified or if they are both specified at the open boundary. It should be one or the other
NameError -- if the dimension of the velocity vector is not equal to the physical dimension of the problem
"""
For now, the open boundary condition only allow to specify velocity and concentration if two continuous phases are considered (inflow) or pressure and concentration if two continuous phases are considered (outflow).
- Coupling between fluid and solid is mainly based on the detection of the mesh cells in which grains are. Computing fluid fraction in a cell or evaluating the fluid velocity at the grain position both imply to know where the grains are in the fluid. The tree structure based on mesh cells is used to locate them.
fluid.set_particles(p.mass(), p.volume(), p.position(), p.velocity(),p.contact_forces(),init=True)
- At this step, the initial state is completely defined and all the variables and parametres have been created. Then the computation can start. The computational loop contains the time integration method solving the fluid and the grains phase:
- Fluid equations solving
while t < tEnd :
time_integration.iterate(fluid, p, dt, min_nsub=5, external_particles_forces=G)
t += dt
This function can be used to solve problem with fluid only, or grains only, or fluid and grains. The external_particle_forces argument contains the external forces that are applied on the grains. For instance, if the grains are falling under falling, this force will be set as:
G = p.mass()*g
The complete documentation of the function and its parameters are given:
def iterate(fluid, particles, dt, min_nsub=1, contact_tol=1e-8, external_particles_forces=None, fixed_grains=False) :
"""Migflow solver: the solver type depends on the fluid and particles given as arguments.
Keyword arguments:
fluid -- fluid structure
particles -- particles structure
dt -- time step
min_nsub -- minimal nsub for the particles iterations
contact_tol -- tolerance for the contact solver
external_particles_forces -- vector of external forces applied on the particles
fixed_grains -- boolean variable specifying if the grains are fixed in the fluid
Raises:
ValueError -- fluid and particles cannot be both None
ValueError -- external_particles_forces must have shape (number of particles,dimension)
"""
The loop also contains instruction to write output files. It is important to note that after moving the grains in the computational loop you have to use the set_particles() function to search the new location of the grains in the mesh.
Optional Arguments
Optional options can be specified during the fluid initialization
- sigma: surface tension (only when there are two fluids see specific case in the Main Features Section)
- volume_drag: coefficient of a linear volume drag force on the entire fluid domain. It can be used to simulate the front and rear shear stress in two-dimensional simulations
- quadratic_drag: coefficient of a quadratic volume drag force on the entire fluid domain. It can be used to simulate the front and rear shear stress in two-dimensional simulations
- drag_in_stab: states if the drag force is in the stabilization term (SUPG and PSPG). If set to False, it stabilizes the computation by adding extra diffusivity.
- drag_coefficient_factor: factor multiplying the drag coefficient
- temporal: enables time derivative (False = steady equations)
- advection: enables advective terms (False = Stokes equations, True = Navier-Stokes equations)
Handy Test Cases List
Prerequisite Directory | Prerequisite Run | Test Case Directory | Test Case | State |
---|---|---|---|---|
~MigFlow/testcases/ | ~MigFlow/testcases/ | |||
avalanch/depot/ | depot.py | avalanch/avalanch1fluid/ | avalanch1fluidnofriction.py | Maintained |
avalanch/depot/ | depot.py | avalanch/avalanch2fluids/ | avalanch2fluids.py | Broken |
avalanch/avalanch2fluids/ | avalanch2fluidsnofriction.py | Broken | ||
couette-2d/ | melangeur.py | Mainained | ||
depot-2d/ | depot.py | Maintained | ||
depot-3d/ | depot.py | Broken | ||
drop-2/SimpleDrop/ | drop.py | Maintained | ||
drop-2d/InteractingDrops/Vert/ | 2VertDrops.py | Maintained | ||
drop-2d/InteractingDrops/Diag/ | 2DiagDrops.py | Maintained | ||
drop-3d/SimpleDrop/ | drop.py | Maintained | ||
drop-3d/InteractingDrops/Vert/ | 2VertDrops.py | Maintained | ||
drop-3d/InteractingDrops/Diag/ | 2DiagDrops.py | Maintained | ||
drop-3d/InteractingDrops/DeifferentDensities | 2DDDrops.py | Maintained | ||
funnel/ | depot.py | Broken | ||
sablier-2d/ | sablier.py | Maintained | ||
shaker/ | shaker.py | Maintained | ||
smallDepot-3d/ | depot.py | Broken |