This tutorial deals with the usage of the unsaturated_2Dflow module to simulate/analyse the flow of water or chemical (or fluid in general) in a variably saturated matrix. The emphasis is not on the finite element implementation; rather it is on the preprocessor and description of how the code is run.

The presentation includes:

An overview of the mathematical equation
The creation of the geometry and physical groups in Gmsh
Meshing the geometry in Gmsh
Executing the code and parameter input arguments
Converting output into vtk
Visualizing the output in paraView

Unsaturated flow model

The mathematical equation describing the movement of water or chemicals in a variably saturated matrix is usually given by Richards' model

\[ \frac{\partial \theta(h)}{\partial t} - \mathrm{div}\,[K(h)\,\mathrm{grad}\,(h - z)] = 0 \]

where \(h\) is the pressure head, \(\theta\) the water (chemical) content is a function of \(h\). The hydraulic conductivity \(K\) is positive definite. The variable \(z\) refers to the vertical distance in the case of vertical flow against the action of gravitational force. However, when the flow is assumed to be restricted in the horizontal plane (i.e., \(z\) is constant), the contribution of \(z\) is ignored.

A full description of the above model requires constitutive laws for the fluid content \(\theta(h)\) and hydraulic conductivity \(K(h)\), and boundary and initial conditions. These constitutive laws depend on the particular material (such as the matrix and the fluid) properties under consideration.

In this example, we consider the modified van Genuchten (MVG) constitutive laws. The volumetric fluid content constitutive equation is given as

\[ \theta(h) = \begin{cases} \theta_r + \frac{\theta_m - \theta_r}{[1+(-\alpha h)^n]m}, & h < h_s\\ \theta_s, & h\leq h_s \end{cases} \]

where \(\theta_r\) and \(\theta_s\),are positive material parameters referred to as residual and saturated fluid contents. The parameter \(h_s\) is the minimal capillary pressure, and \(\theta_m\) is a fictitious modelling parameter. \(\alpha\) and \(n\) are empirical shape parameters. In addition, the constitutive equation for the hydraulic conductivity is given by

\[ \theta(h) = \begin{cases} K_s\,K_r(h), & h < h_s\\ K_s, & h\leq h_s \end{cases} \]

where \(K_s\) is the saturated hydraulic conductivity and \(K_r\) is the unsaturated hydraulic conductivity, which is related to \(h\) through the effective saturation \(S_e\) (volume fraction representing the composition of fluid in the matrix-fluid mixture) by

\[ K_r = S_e^\ell \bigg[ \frac{1-F(S_e)}{1-F(1)}\bigg]^2 \]

where \(\ell\) is the pore connectivity parameter usually assumed to be \(0.5\), and the function \(F\) is given by

\[ F = \bigg[1-\bigg(\frac{\theta_s - \theta_r}{\theta_m - \theta_r}S_e\bigg)^{1/m}\bigg]^m \]

and the effective saturation is given by

\[ S_e = \frac{\theta - \theta_r}{\theta_s - \theta_r} \]

Here \(m = 1 - 1/n\), \(n>1\).

The values of the parameters used in this problem represents \(K_s = 1.0, \ell = 0.5, \theta_r = 0.045, \theta_s = 0.43, \theta_m = 0.43, n = 5.38, \alpha = 1.8125\) and \(h_s = 0.00\).

The boundary conditions are divided into two (see Fig 1); one that prescribes the pressure head \(h\) on a given part of the boundary \(\Gamma_1\), the other one prescribes the flux on the rest of the boundary \(\Gamma_2 = \partial\Omega - \Gamma_1\). Usually the second type of boundary condition is zero flux which corresponds to an insulation of the boundary, i.e., no material is coming in and out of the boundary.

The initial condition prescribes the pressure head \(h\) at the beginning of the simulation, in this example we consider a homogeneous initial condition, i.e., \(h = h_s\)

Creating the geometry

The problem considered in this tutorial is the lateral flow of a dissolvable fluid in a simple rectangular paper slide as shown in Fig 1. Part of the boundary labeled \(\Gamma_1\) is in contact with a source of the fluid associated with pressure head \(h = C\). The rest of the boundary, i.e., \(\Gamma_2\), is insulated.

Figure 1: Domain and boundary decomposition

The geometry of the paper slide \(\Omega\) has been created in Gmsh, a open-source mesh generating software. Later, the geometry has been saved in to a .geo file that can be opened in any text editor, which in turn provides an alternative way of modifying it.

A simple two step procedure to create this rectangular geometry is as follows: first open the Gmsh software, then go to

File -> New -> Choose directory -> give name (e.g., paper_network.geo) -> save -> choose OpenCASCADE from the dialogue box

Next, on the left click on the Module; this opens three branches namely Geometry, Mesh and Solver, then follow

Geometry -> Elementary entities -> Add -> Rectangle

A dialogue page appears as shown below, press "shift" key to hold position and enter X = 0, Y = 0, Z = 0 as the coordinates of the left corner of the rectangle and enter DX = 0.02 and DY = 0.005 as the dimensions, make sure that radius rounded is set to 0. Then click Add on the dialogue box.

To finish creating the rectangle press "q" key to abort. Now the rectangle is created.

Physical groups

Physical groups refers to regions in the domain where various values are prescribed, for example, boundary conditions, initial conditions, material parameters etc. For the problem under consideration, we only need two types of physical groups, initial and boundary conditions. Now go back to our rectangle in Gmesh. Note that every Point, Line or Plane is allocated a unique ID by Gmsh as it is created, automatically. To view such ID's, under the Tools Menu choose Options, and in the dialogue box that appears press on the Geometry menu on the left then under Visibility check Curve labels and Surface labels then close the dialogue box.

To create the boundary entity where pressure head \(h\) is prescribed i.e., \(\Gamma_1\) (Note that it is Line 4), follow

Module -> Geometry -> Physical groups -> Add -> Curve

Give it a name on the dialogue box appearing, e.g., ESSENTIAL, uncheck Automatic numbering and enter 1 in the Numer field. Hover the the cursor over the rectangle and Click on Line 4 to select it (then its color is changed to red), and Press "e" to end selection and "q" to abort.

Next to choose the Physical Group where the initial condition is prescribed (in this case, only one initial condition which is prescribed over the entire rectangle, i.e. Plane 1, is considered), follow

Module -> Geometry -> Physical groups -> Add -> Surface

on the dialogue box give it a name, e.g., INITIAL. Select uncheck Automatic numbering and enter 2 on the number field, and the hover over the rectangle and select Plane 1, afterwards press "e" to end selection and "q" to abort.

This results in the file paper_network.geo

// Gmsh project created on Mon Jan 27 21:51:41 2020
SetFactory("OpenCASCADE");
Rectangle(1) = {0, 0, 0, 0.02, 0.005, 0};
Physical Curve("ESSENTIAL", 1) = {4};
Physical Surface("INITIAL", 2) = {1};

which you can open with any text editor from which you can modify the geometry, for example, for later use we can create variables lx and ly for dimensions of the rectangle as follows

// Gmsh project created on Mon Jan 27 21:51:41 2020
SetFactory("OpenCASCADE");
lx = 0.02;
ly = 0.005;
Rectangle(1) = {0, 0, 0, 0.02, 0.005, 0};
Physical Curve("ESSENTIAL", 1) = {4};
Physical Surface("INITIAL", 2) = {1};

Meshing

Creating the mesh is simple, once we have the .geo file setup correctly, you load the modified paper_network.geo file as follows

File -> Open... -> select the file you want to load (in this case is paper_network.geo)

Then to create the mesh,

Module -> Mesh -> 2D

The rectangle is the meshed into relatively coarser triangles. If you want to refine it, just press Refine by splitting under Mesh. (you can repeat pressing it as many times as required, to get to an acceptable level of refinement.) To save the resulting mesh you need to go to

File -> Export -> select a folder you want to save in -> give it a name , e.g., paper_network.msh -> save

Choose Version 2 ASCII on the Format field, and press OK.

At this point, we have everything we need to run the programme. But there is one important ingredient that needs to be considered before we use the mesh in the MoFEM context. First, note that .msh file is not directly used by MoFEM. To convert it to a mesh file readable by MoFEM, i.e., .h5m format, we need something referred to as a configuration file of .cfg format, which basically contains rules to convert the Physical groups from Gmsh into the corresponding NODESET, SIDESET and BLOCKSET that are understandable by MoFEM. For example, for the mesh generated above the configuration file bc.cfg looks like

[block_1] 
 
id=1001
add=BLOCKSET
name=ESSENTIAL
 
 
[block_2]
 
id=1002
add=BLOCKSET
name=INITIAL

Note: it is very important that you need to make sure that each Physical Group with its ID matches the block in bc.cfg. This means that, for example, the physical group "ESSENTIAL" has been assigned ID 1 which should be corresponded by [block_1], id=1001, name=ESSENTIAL. Whereas, the Physical Group "INITIAL" with ID 2 should be matched by [block_2], id=1002, name=INITIAL. If there is any mismatch it leads to an entirely wrong result.

Executing the code

All the necessary steps to run the code are encapsulated in following bash script file run.sh

#!/bin/bash
  
if  [[ $1 = "-my_file" ]]; then
    echo "Running with input file $2 ..."
    add_meshsets -my_file $2 -meshsets_config bc.cfg -output_file paper_network_pre.h5m
else
  echo "Give the input mesh file name ..."  
  exit
fi
 
rm -rf out*
 
mofem_part -my_file paper_network_pre.h5m -output_file paper_network_final.h5m -my_nparts 1 -dim 2 -adj_dim 1
if [[ -z "$3" ]]; then
    time mpirun -np 1 ./unsatu2dFlow_prob -file_name paper_network_final.h5m 2>&1 | tee log
elif [[ $3 = "-saturated_conductivity" ]]; then
    time mpirun -np 1 ./unsatu2dFlow_prob -file_name paper_network_final.h5m -saturated_conductivity $4 2>&1 | tee log
fi
 
../../basic_finite_elements/nonlinear_elasticity/do_vtk.sh out*

Among other things, the above script utilises a file named param_file.ptsc which contains input for the code such as the mesh, order of approximation and parameters for petsc's time solver, linear/non linear solvers, preconditioner, etc. It looks like as follows

 -order 1                      # order of approximation
 -saturated_conductivity 0.1   # set default value for K_s 
 -save_every_nth_step  4       # time-steps to save output at 
 
 
# petsc time solver (ts) arguments
 -ts_monitor                    # to print messages about the time step computations                  
 -ts_exact_final_time stepover  # rule when to stop at the final time
 -ts_dt 0.0000005               # time step size
 -ts_max_time 0.0005          # final time
 -ts_adapt_type none            # no time-adaptivity, i.e., fixed step size
 
 # petsc linear solver arguments, fixed
 -ksp_type gmres 
 -pc_type lu 
 -pc_factor_mat_solver_type mumps 
 
 # petsc nonlinear solver arguments, no need to change them
 -snes_atol 1e-7
 -snes_rtol 1e-12 
 -snes_monitor 
 -snes_lag_jacobian 1

Finally, the code is simply executed using the script file as

./run.sh -my_file paper_network.msh -saturated_conductivity 1.0

Here, the arguments for run.sh script are

following \(\texttt{-my_file}\) is the input mesh file, e.g., \(\texttt{paper_network.msh}\), and
following \(\texttt{-saturated_conductivity}\) is the input value for saturated conductivity \(K_s\), e.g., 1.0, (a default value 0.1, as listed in \(\texttt{param_file.petsc}\), is used if this argument is not specified).

Visualizing the output

Assuming paraView is installed on your machine, open paraView and load the file out_level*

File -> Open... -> go to the output file directory, then select the out_level* -> OK -> Apply (under Properties)

the pressure head \(h\) is output as h1, then change GLOBAL_ID to h1 in the dropdown option

press the play button to start playing through the time snapshots, for example, the snapshot at the step 106 looks

To plot the time profile of the horizontal line through the mid-section of the rectangle

Plot Over Line icon -> In Properties select X Axis button -> in Series Parameters field select only h1

Then you get

To Plot the time profile of a given point, first close the LineChartView1 window and delete the PlotOverLine1 in the Pipeline Browser window. Then press on the Select Points On (d) icon at the top the RenderView1 window to trigger the selection

hover over the domain and select a point somewhere in the middle -> press the Plot Selection Over Time button -> select Apply -> unselect all variables except h1 in the

Propertties window

paraView is a powerful visualisation and analysis tool. It has reach functionalities, you may play around to get a some understanding of it.

Exercise

Estimating speed of the front of pressure head. To do this, you need a grid over in the RenderView of VTK, in which the output is lying inside it. and then put the Time at a certain point say 56. You can get the grid in the following way, on the paraView page with the output loaded, in the Properties view under Annotations check Axis Grid and press Edit button next to it. A dialogue box appears, in it press the gear as shown,

and on the appearing dialogue box check Show Grid button and press the gear as shown,

now on the final dialogue box appearing check Show Grid, Show Edges and Show Ticks,

finally press Apply -> OK then close the dialogue box.

(since the save_every_nth_step is set at 4, it implies that output 56 means t = 56X4Xdt where dt is the time step size) then for example h1 = -0.5 is roughly at 0.008. Hence we effectively say that at t = 242Xdt the speed is 0.008/(242Xdt). Do this at various time steps and record the speeds. Now do the speed test for different values of saturated_conductivity \(K_s\). And observe that how the size of \(K_s\) influences the speed.
Add more feature to the geometry. For example, a geometry consisting of two rectangles separated by an angle \(\pi/3\), as shown

and the geo file two_rec.geo is
// Gmsh project created on Mon Jan 27 21:51:41 2020

lx = 0.02;

ly = 0.005;

SetFactory("OpenCASCADE");

Rectangle(1) = {0, 0, 0, lx, ly, 0};

Rotate {{0, 0, 1}, {lx, ly/2, 0}, Pi/3} { Duplicata { Surface{1}; }}

BooleanUnion { Surface{1}; } { Surface{2}; }

Recursive Delete { Surface{1}; Surface{2}; }

Physical Curve("ESSENTIAL", 1) = {8};

Physical Surface("INITIAL", 2) = {3};

Understand what each line is about. Create the corresponding mesh and execute the code with it. Now create a geometry such as this

Create a mesh with it and run the code using the mesh.

Table of Contents