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.

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.

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.

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

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.

_fig_seep4_cont.png

that's all

Recommended Posts

Two-dimensional saturated-unsaturated osmotic flow analysis with Python
Two-dimensional unsteady heat conduction analysis with Python
Voice analysis with python
Voice analysis with python
Data analysis with Python
Two-dimensional elastic skeleton geometric nonlinear analysis with Python
[Python] Morphological analysis with MeCab
[Co-occurrence analysis] Easy co-occurrence analysis with Python! [Python]
Sentiment analysis with Python (word2vec)
Planar skeleton analysis with Python
Japanese morphological analysis with Python
Muscle jerk analysis with Python
3D skeleton structure analysis with Python
Impedance analysis (EIS) with python [impedance.py]
Text mining with Python ① Morphological analysis
Data analysis starting with python (data visualization 1)
Logistic regression analysis Self-made with python
Data analysis starting with python (data visualization 2)
[In-Database Python Analysis Tutorial with SQL Server 2017]
Marketing analysis with Python ① Customer analysis (decyl analysis, RFM analysis)
Machine learning with python (2) Simple regression analysis
2D FEM stress analysis program with Python
Tweet analysis with Python, Mecab and CaboCha
Principal component analysis with Power BI + Python
Data analysis starting with python (data preprocessing-machine learning)
Python: Simplified morphological analysis with regular expressions
You can do it with Python! Structural analysis of two-dimensional colloidal crystals
Statistics with python
Scraping with Python
Python with Go
Data analysis python
[Various image analysis with plotly] Dynamic visualization with plotly [python, image]
Twilio with Python
Medical image analysis with Python 1 (Read MRI image with SimpleITK)
Integrate with Python
Play with 2016-Python
AES256 with python
Tested with Python
python starts with ()
[Python] Flow from web scraping to data analysis
with syntax (Python)
Bingo with python
Zundokokiyoshi with python
Static analysis of Python code with GitLab CI
Easy Lasso regression analysis with Python (no theory)
Flow to complete Slack authentication with Flask (Python)
Excel with Python
Microcomputer with Python
Cast with python
Flow of creating your own package with setup.py with python
Text mining with Python ① Morphological analysis (re: Linux version)
Data analysis for improving POG 1 ~ Web scraping with Python ~
Collecting information from Twitter with Python (morphological analysis with MeCab)
[OpenCV / Python] I tried image analysis of cells with OpenCV
Reading Note: An Introduction to Data Analysis with Python
Data analysis environment construction with Python (IPython notebook + Pandas)
Calculate the regression coefficient of simple regression analysis with python
3. Natural language processing with Python 4-1. Analysis for words with KWIC
Challenge principal component analysis of text data with Python
Planar skeleton analysis with Python (4) Handling of forced displacement
Summary of the basic flow of machine learning with Python