v0.13.2 |

Loading...

Searching...

No Matches

Unsaturated flow in 2D

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

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\)

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.

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 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};

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.

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).

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.

**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,**Show****Grid**button and press the gear as shown,

**Show****Grid**,**Show****Edges**and**Show****Ticks**,**Apply**->**OK**then close the dialogue box.

**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**two_rec.geo**isUnderstand what each line is about. Create the corresponding mesh and execute the code with it. Now create a geometry such as this// Gmsh project created on Mon Jan 27 21:51:41 2020lx = 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};

Generated by Doxygen 1.9.5 and hosted at