Two-dimensional saturated-unsaturated osmotic flow analysis with Python

(Corrected on November 19, 2017)

Corrected the program due to some program mistakes

Overview

Rough impression

Since I made a 2D FEM stress analysis program, I also made a 2D saturated-unsaturated permeation flow analysis program. The element is an isoparametric element having 1 element, 4 nodes, and 4 Gauss integral points, as in the 2D stress analysis.

Actually, it can be used, but the number of nodes is 1317, the number of elements is 1227, the number of convergence calculations is 10, and the calculation time is 3.7 seconds. If the problem is of this scale, no particular stress is felt. Numpy's matrix operations and routines for solving simultaneous linear equations are probably excellent.

The items to keep in mind when programming are as follows.

All arrays are Numpy arrays. No list is used. Avoid using + for as much as possible and use Numpy matrix operations. However, please note that unexpected results may occur due to the product of one-dimensional vectors.
`x = numpy.linalg.solve (A, b)` is enough to solve simultaneous linear equations. Boundary condition processing makes the matrix asymmetric, so Cholesky decomposition cannot be used either way.
By doing the above, the outlook of the program itself will be very good.
The input / output section gets messed up. Is this unavoidable?

Saturated osmotic flow analysis

Saturated steady-state osmotic flow analysis is a problem that solves the following simultaneous linear equations.

[k]\\{h\\}=\\{q\\}

\begin{bmatrix}k_{11} & k_{12} & \dots \\\ k_{21} & k_{22} & \dots \\\ \dots & \dots & \dots \\\ \dots & \dots & k_{nn}\end{bmatrix} \begin{Bmatrix}h_1 \\\ h_2 \\\ \dots \\\ h_n \end{Bmatrix} =\begin{Bmatrix}q_1 \\\ q_2 \\\ \dots \\\ q_n \end{Bmatrix}

\begin{bmatrix}k_{11} & k_{12} & \dots \\\ k_{21} & k_{22} & \dots \\\ \dots & \dots & \dots \\\ \dots & \dots & k_{nn}\end{bmatrix} \begin{Bmatrix}h_{\text{unknown}} \\\ h_{\text{unknown}} \\\ \dots \\\ h_{\text{given}} \end{Bmatrix} =\begin{Bmatrix}q_{\text{given}} \\\ q_{\text{given}} \\\ \dots \\\ q_{\text{unknown}} \end{Bmatrix}

Here, $ \ {q \} $ is the nodal flow vector, $ \ {h \} $ is the total head vector, and $ [k] $ is the hydraulic matrix, which has the following characteristics.

Permeability matrix $ [k] $ is a constant,
The total head vector $ \ {h \} $ corresponding to the known component of the nodal flow vector $ \ {q \} $ is unknown.
The total head vector $ \ {h \} $ corresponding to the unknown component of the nodal flow vector $ \ {q \} $ is a known number.

Therefore, the solution can be obtained by solving the simultaneous equations once after processing the relationship between the known number and the unknown number.

(Processing of the relationship between known and unknown numbers to solve simultaneous linear equations)

Flow rate excluding $ q_i $, $ q_j $ and known $ h_i $, $ h_j $, all heads excluding $ h_i $, $ h_j $ and $ q_i $, $ q_j $ in the permeable matrix The process is performed so that k_ {ii} $ and $ k_ {jj} $ are set to $ 1 $, and the other elements of the $ i $ and $ j $ columns are set to zero. It also transfers the effects of the $ i $ and $ j $ columns in the hydraulic matrix to the right side.

\begin{bmatrix} k_{11} & \ldots & 0 & \ldots & 0 & \ldots & k_{1n} \\\ \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\\ k_{i1} & \ldots & 1 & \ldots & 0 & \ldots & k_{in} \\\ \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\\ k_{j1} & \ldots & 0 & \ldots & 1 & \ldots & k_{jn} \\\ \vdots & \ddots & \vdots & \ddots & \vdots & \ddots & \vdots \\\ k_{n1} & \ldots & 0 & \ldots & 0 & \ldots & k_{nn} \end{bmatrix} \begin{Bmatrix} h_1 \\\ \vdots \\\ -q_i \\\ \vdots \\\ -q_j \\\ \vdots \\\ h_n \end{Bmatrix} =\begin{Bmatrix} q_1 \\\ \vdots \\\ 0 \\\ \vdots \\\ 0 \\\ \vdots \\\ q_n \end{Bmatrix} -h_i \begin{Bmatrix} k_{1i} \\\ \vdots \\\ k_{ii} \\\ \vdots \\\ k_{ji} \\\ \vdots \\\ k_{ni} \end{Bmatrix} -h_j \begin{Bmatrix} k_{1j} \\\ \vdots \\\ k_{ij} \\\ \vdots \\\ k_{jj} \\\ \vdots \\\ k_{nj} \end{Bmatrix}

Saturated-unsaturated osmotic flow analysis

The simultaneous equations in the saturated-unsaturated osmotic flow analysis have the following characteristics.

It is necessary to consider three types of boundary conditions: known flow rate, known total head, and boundary where infiltration surface appears.
Either the flow rate or the total head is known except for the boundary where the infiltration surface appears.
On the boundary where the infiltration surface appears, both the flow rate and the total head are unknown.
However, since the boundary where the infiltration surface appears is basically the boundary in contact with the atmosphere, the flow rate is only outflow (0 or less in the program) unless the inflow of rain is taken into consideration. The pressure head is also restricted to 0 or less.
The hydraulic conductivity component in the saturated region is constant, but the hydraulic conductivity component in the unsaturated region is a function of suction pressure (negative pressure head).

Therefore, it is necessary to determine the unknown by convergence calculation. The method of convergence calculation was a sequential substitution method in which the initial values of all heads were given to all nodes and the obtained solution was used as the input value for the next convergence calculation. Since this method assumes the heads of all nodes, it can be set from the heads assuming the hydraulic conductivity of all elements, and the calculation flow can be simplified. It seems that it is effective to give the minimum total head in the saturation region as the initial value of the total head except for the boundary where the total head is known, considering the accuracy and convergence of the solution.

Unsaturated permeability (van Genuchten model)

For the unsaturated permeability characteristic model of the ground, there is also a method of inputting the saturation-suction pressure relationship and the saturation-unsaturated permeability coefficient ratio in a table, but for analysis, it is better to handle it with a continuous function because the program can be simplified. The van Genuchten model was adopted.

S_e=\left\\{\frac{1}{1+\left(\alpha\cdot h_s\right)^n}\right\\}^m \qquad n=\frac{1}{1-m} \quad (0

K_r=(S_e)^{0.5}\cdot\left\\{1-\left(1-S_e{}^{1/m}\right)^m\right\\}^2 \qquad (0 \leqq S_r, K_r \leqq 1)

K=K_r \cdot K_0

S_e : Degree of saturation
K_r : Relative hydraulic conductivity function
h_s : Suction head
\alpha : Scaling parameter
m : Non-dimensional parameter
K : Permiability coefficient
K_0 : Saturated permiability coefficient

program

If you put the program on it, it will be long, so put a link to your own HP.

Input data format

npoin  nele  nsec  koh  koq  kou  idan
Ak0  alpha  em
..... (1~nsec) .....
node-1  node-2  node-3  node-4  isec
..... (1~nele) .....
x  z  hvec0
..... (1~npoin) .....
nokh  Hinp
..... (1~koh) .....
nokq  Qinp
..... (1~koq) .....
noku
..... (1~kou) .....

npoin, nele, nsec	Number of nodes, number of elements, number of material characteristics
koh, koq, kou	Number of head designated nodes, flow rate specified nodes, infiltration boundary nodes
idan	Analysis cross section (0: vertical cross section, 1: horizontal cross section)
Ak0	Saturated hydraulic conductivity of the element
alpa, em	Unsaturated permeability of elements
node-1, node-2, node-3, node-4, isec	Node numbers that make up the element, material property number of the element
x, z, hvec0	x and z coordinates of the node, initial head value for convergence calculation
nokh, Hinp	Head designation node number and head value
nokq, Qinp	Flow rate specified node number and flow value
noku	Infiltration boundary node number

Output data format

npoin  nele  nsec   koh   koq   kou  idan
    9     4     1     6     0     0     1
  sec             Ak0           alpha              em
    1   1.0000000e-05   0.0000000e+00   0.0000000e+00
 node               x               z            hvec            qvec   koh   koq   kou
    1   0.0000000e+00   0.0000000e+00   1.0000000e+01   0.0000000e+00     1     0     0
    2   1.0000000e+00   0.0000000e+00   1.0000000e+01   0.0000000e+00     1     0     0
    3   2.0000000e+00   0.0000000e+00   1.0000000e+01   0.0000000e+00     1     0     0
    4   0.0000000e+00   1.0000000e+00   0.0000000e+00   0.0000000e+00     0     0     0
    5   1.0000000e+00   1.0000000e+00   0.0000000e+00   0.0000000e+00     0     0     0
    6   2.0000000e+00   1.0000000e+00   0.0000000e+00   0.0000000e+00     0     0     0
    7   0.0000000e+00   2.0000000e+00   0.0000000e+00   0.0000000e+00     1     0     0
    8   1.0000000e+00   2.0000000e+00   0.0000000e+00   0.0000000e+00     1     0     0
    9   2.0000000e+00   2.0000000e+00   0.0000000e+00   0.0000000e+00     1     0     0
 node  Hinp
    1   1.0000000e+01
    2   1.0000000e+01
    3   1.0000000e+01
    7   0.0000000e+00
    8   0.0000000e+00
    9   0.0000000e+00
 elem     i     j     k     l   sec
    1     1     2     5     4     1
    2     2     3     6     5     1
    3     4     5     8     7     1
    4     5     6     9     8     1
 node            hvec            pvec            qvec   koh   koq   kou
    1   1.0000000e+01   1.0000000e+01   2.5000000e-05     1     0     0
    2   1.0000000e+01   1.0000000e+01   5.0000000e-05     1     0     0
    3   1.0000000e+01   1.0000000e+01   2.5000000e-05     1     0     0
    4   5.0000000e+00   5.0000000e+00  -6.7762636e-21     0     0     0
    5   5.0000000e+00   5.0000000e+00   2.7105054e-20     0     0     0
    6   5.0000000e+00   5.0000000e+00   0.0000000e+00     0     0     0
    7   0.0000000e+00   0.0000000e+00  -2.5000000e-05     1     0     0
    8   0.0000000e+00   0.0000000e+00  -5.0000000e-05     1     0     0
    9   0.0000000e+00   0.0000000e+00  -2.5000000e-05     1     0     0
 elem              vx              vz              vm              kr
    1  -5.5511151e-21   5.0000000e-05   5.0000000e-05   1.0000000e+00
    2   5.5511151e-21   5.0000000e-05   5.0000000e-05   1.0000000e+00
    3  -5.5511151e-21   5.0000000e-05   5.0000000e-05   1.0000000e+00
    4   5.5511151e-21   5.0000000e-05   5.0000000e-05   1.0000000e+00
Total inflow =  1.0000000e-04
Total outflow= -1.0000000e-04
Max.velocity in all area     =  5.0000000e-05 (ne=1)
Max.velocity in inflow area  =  5.0000000e-05 (ne=1)
Max.velocity in outflow area =  5.0000000e-05 (ne=3)
iii=2  icount=9  kop=0
n=9  time=0.007 sec

node, hvec, pvec, qvec	Node number, total head, pressure head, flow rate
elem, vx, vz, vm	Element x / z direction flow velocity and synthetic flow velocity
Kr	Element unsaturated permeability characteristics (1: saturated, 1>: unsaturated)
Total inflow, Total outflow	Total inflow and outflow to the model area
iii, icount	Number of convergence calculations, number of freedoms that satisfy the calculation completion condition
kop	The number of infiltration surface boundary nodes where the pressure head is 0
n, time	total degrees of freedom, calculation time

Output example

The figure which pulled the contour of the pressure head 0, 5, 10, 15m by the simple contour creation program. The water depth is 18 m on the upstream side and 4 m on the downstream side. The contour diagram itself is at the amateur level, but the calculation seems to work well. The pressure heads and dam specifications shown in red are additionally written in the Mac Preview.

It is better to leave the drawing here to the contour function of GMT (Generic Mapping Tools) and it will look better for the report.

By the way, it is a well-known fact that in an actual dam, the pressure head drops near the boundary between the core and the filter, and the pressure drop inside the core does not occur as shown in the analysis results. Please understand that this is just an analysis case.

that's all