Tutorial 2 - gnomeCreative/HYBIRD GitHub Wiki

Non-Newtonian Fluids with LBM

This second example introduces non-Newtonian fluids in LBM. Non-Newtonian fluids is an umbrella term for any fluid that does not exhibit a linear relationship between stress and strain rate.

The Bingham plastic.

Here, we work with a relatively simple non-Newtonian fluid called a Bingham plastic. This is a fluid that does not shear unless a minimum stress, called a yield stress, $\sigma_y$ is reached. The yield stress is therefore a material parameter. Note that this is a shear stress, but we use a $\sigma$ symbol in order to distinguish it from the relaxation time $\tau$. When the applied stress exceeds the yield stress, the Bingham plastic behaves according to the following law:

$\sigma=\sigma_y+\nu_\textrm{p}\dot{\gamma}$

where $\sigma$ is the flow resistance to shear, $\dot{\gamma}$ is the shear rate, and $\nu_\textrm{p}$ is the plastic viscosity, which is the second material parameter (the first being the yield stress). The model is depicted in the figure above.

In analogy with Newtonian fluids, we call the ratio between shear stress and shear rate an (apparent) viscosity:

$\nu=\frac{\sigma}{\dot{\gamma}}=\nu_\textrm{p}+\frac{\sigma_y}{\dot{\gamma}}$

When the shear rate approaches zero (the no-flow condition), the apparent viscosity diverges to $+\infty$.

In this tutorial, we consider a free-surface flow on an infinite incline with constant slope $\theta$, driven by gravity, see the figure below

Geometry for the Bingham incline flow.

We set a reference frame such that the $x$ direction is aligned with the incline, and the imposed pressure gradient (the external forcing $f_x$) is also aligned with $x$. The domain has periodic boundary conditions in the $x$ and $y$ directions. In the $z$ direction, a no-slip condition is applied at the bottom, and the free-surface condition is applied at the top. This type of flow is analogous to the Poiseuille flow, because the free-surface condition, being stationary, is analogous to a symmetry condition.

For this case, an analytical solution of the Navier-Stokes equations (adapted to the Bingham law) exists. Since the shear stress in the fluid must exceed the yield stress, we can only have flow if the stress induced by gravity is larger than the yield stress. This is true if

$H>\frac{\sigma_{y}}{\rho_\mathrm{f} f_x}$

or, in other terms,

$\sigma_{y}<\rho_\mathrm{f} f_x H$

If the above is true, the type of flow that is generated consists of a plug (i.e., an area with equal velocity and no internal deformation) lying on top of a shear layer, as shown in the figure above. The thickness of the plug can again be derived using simple equilibrium conditions, and is:

$H_0=\frac{\sigma_{y}}{\rho_\mathrm{f} f_x}$

In the shear layer ($z\lt H-H_0$), velocity follows a parabolic profile:

$u_x(z)=\frac{\rho_\textrm{f} f_x}{\nu_\mathrm{p}} z \left(H-H_0 -\frac{z}{2} \right)$.

From this result, the velocity of the fluid in the plug region can be calculated as

$u_x\left(H_0\right)=u_\text{B,max}=\frac{\rho_\text{f} f_x}{2\nu_\mathrm{p}}\left(H-H_0\right)^2$.

We will try to recover this result using LBM.

The main issue here is that the basic LBM formulation is suited only for Newtonian fluids, i.e., fluids with a constant finite viscosity. The LBM code cannot reproduce this result exactly, because the simulation of a Bingham model requires the treatment of the viscosity $\nu$ as a variable, which reaches its minimum value $\nu_\textrm{p}$ for $\dot{\gamma}\rightarrow\infty$ and diverges for $\dot{\gamma}\rightarrow 0$. Such a large range of viscosity variation would cause instability in the code. This is managed by bounding the viscosity range by two values, $\nu_\textrm{min}$ and $\nu_\textrm{max}$. How these values are derived is explained in the LBM background page, but it is important to highlight that they are a function of the time step $\Delta t$ and the lattice resolution $\Delta S$, and can therefore be controlled, albeit with constraints, by the user.

The trilinear model obtained by approximating the Bingham plastic.

The resulting approximated Bingham rheology is a tri-linear model, as sketched in the figure above. However, if $\nu_\textrm{p}>\nu_\text{min}$, which is usually easy to obtain by adapting the time step $\Delta t$, the rheology reduces to a bi-linear model, with a range of simulated viscosities wide enough to yield a good approximation for most cases.

In this tutorial, we impose density $\rho_\textrm{f}=1000~\textrm{kg}/\textrm{m}^3$, dynamic plastic viscosity $\nu_\textrm{p}=1.0~\textrm{kg}/\textrm{m}\textrm{s}$, external forcing $f_x=1~\textrm{m}/\textrm{s}^2$, and flow height width $H=1~\textrm{m}$. With these settings, the maximum velocity should be $u_\textrm{max}=1.25~\textrm{m}/\textrm{s}$. To speed up the computations, we discretize using a cuboid with dimensions 1×1×1201×1×120 in the three directions, therefore a 1D domain. Only the first $100$ nodes in the direction $z$, will be filled with fluid, the rest will be "gas" cells above the free surface.

Setting the Environment

All configuration files should be set inside the tutorials/Bingham_flow folder. Here as well (as in Tutorial 1) the DEM part is deactivated. The configuration file is already ready for use (bingham_flow.cfg), but the most important features are described here.

We are using only the LBM solver, and we do need the free-surface algorithm for this, so we set:

• demSolver=0 No DEM.
• lbSolver=1 LBM is on.
• freeSurfaceSolver=1 free-surface solver is on. However, the free surface will not evolve, because the initial condition corresponds already to the stationary condition.
• forceFieldSolver=1 We need an external forcing for our flow to develop.

In this case, a minimal setup is necessary for the output files:

maxTime=20.0 The code will run for $20~\textrm{s}$ simulation time. restartFluid=false Computation is from scratch, no restart file is required. fluidRecycleExpTime=false No restart files will be generated. screenExpTime=0.1 The code will show a state line on screen every $0.1~\textrm{s}$ simulation time. fluidExpTime=0 The code will not generate fluid Paraview files. fluidLagrangianExpTime=1.0 The code will generate Lagrangian fluid Paraview files every $1.0~\textrm{s}$ simulation time.

The domain boundary conditions are periodicity in xx and yy, and no-slip walls on zz. Note that in zz the upper boundary condition will not be used because the fluid will not touch this boundary. This can be set by specifying:

• boundary0=periodic
• boundary1=periodic
• boundary2=periodic
• boundary3=periodic
• boundary4=stat_wall
• boundary5=stat_wall

Domain setup:

latticeSpacing=0.002 This is $\Delta S$. See below for its influence on the flow viscosity resolution. domainSizeX=0.002 The lattice spacing in $x$ is equal to 1 lattice spacing, so that we work with a single column of fluid cells. domainSizeY=0.0002 The lattice spacing in $y$ is equal to the lattice spacing, so that we work with a single column of fluid cells. domainSizeZ=0.12 This is the flow height $H=0.1~\textrm{m}$ plus a thin layer of "gas" cells on top. fluidMinX=-100.0 and fluidMaxX=100.0 Since we set the minimum and maximum borders of the fluid to be beyond the domain borders, the fluid will fill the whole space (same for the direction $y$). `fluidMinZ=0.0

Running the Code

To run the code, the procedure is identical to the one explained in Tutorial 1. However, the configuration file is now called bingham_flow.cfg. Therefore, the command to be issued in the terminal is:

./hybird.exe -c bingham_flow.cfg -d results -n time

The code will run and print on the screen the value of the maximum velocity. You can see it increasing slowly towards convergence. You can also check the amount of material that is in a "plastic" state, i.e., unyielded. Of course, the approximated implementation of viscosity makes so that the "unyielded" material will be the one that has reached the maximum allowable viscosity $\nu_\textrm{max}$.

The code will also create a subfolder with a timestamp inside the results folder. To visualize the results, open Paraview and load the Lagrangian fluid files. Note that Lagrangian and Eulerian fluid files differ in the style of visualization, not in the information that they display. The Lagrangian files display every fluid cell as a point. This precludes the possibility to draw contours and interrogate variables over lines, but they open up other ways to visualize variables. As an exercise, try to run the code again, this time asking to generate Eulerian files as in Tutorial 1, and visualize the velocity profile.

Lagrangian visualization of the Bingham channel flow.

Here we focus on the type of representation that works well with Lagrangian files. Loading these files will result in a column of fluid nodes, displayed as points dataset. Because the nodes are aligned, Paraview will use automatically a 2D representation. However, we should switch the reference plane to $xz$, by selecting the "Set view direction to +Y" button, highlighted in red.

To visualize the flow velocity, we can use arrows, which are provided by the filter Glyph, shown in the figure below with a green square. The glyph filter needs a vector for scaling its size, and a vector for scaling its orientation. In both cases, we should choose the field "v", which contains the velocity vector. We can then press the "update" button (blue square) to rescale the arrows automatically. To activate the filter, it is necessary to press the "Apply" button. The arrow can also be coloured using the field v (Magnitude), see the yellow square. In the first time step, no arrows will be displayed, as velocity is zero everywhere. If we press "play", Paraview will automatically run through all available files. It might be necessary to stop and rescale the arrow glyphs.

Graphical visualization of the velocity field in the Bingham free-surface flow.

We can also visualize scalars using glyphs. As an example, we could use spheres. Change the glyph type from "arrow" to "sphere" (blue square). Set their radius to be equal to half the lattice spacing ($0.002/2=0.001~\textrm{m}$), then reset the orientation and scale arrays, and choose a "scale factor" of 1 (red squares), and press "Apply". At this point, a column of identical spheres should be created. We can colour these spheres based on a scalar field, for example viscosity.

Graphical visualization of the viscosity field in the Bingham free-surface flow.

Since viscosity has a large variability, we can obtain a better representation by choosing a log colour scale (yellow square). It is now wise to rescale the colours to the automatic range (purple square). Note how at the end of the simulation, the upper half of the fluid (the plug) has converged to the maximum allowable value of viscosity.