Python
有限要素法
FEM
浸透流解析

Pythonによる2次元飽和ー不飽和定常浸透流解析プログラム(四角形要素)の理論を少し詳しく解説する

はじめに

浸透流解析プログラムについて,三角形要素に続き,四角形要素によるものの解説を作ってみた.
三角形要素用のプログラムとの違いは,入力データやマトリックス組み立てを四角形要素用にしていること,透水性マトリックス作成に数値積分を使っていることである.

解析理論は三角形要素用のものと全く同じ,四角形要素の取扱は2次元応力解析プログラムのものと同じなので,ここでは,「プログラミングのために知っているべきこと」として,四角形要素の2次元応力解析で述べた四角形要素の取扱,三角形要素の浸透流解析で述べた有限要素式の結論のみを記載し,「四角形要素による定式化」で四角形要素を用いた透水性マトリックスの作成,要素平均流速計算方法を述べている.

このシリーズの投稿のリンク

プログラミングのために知っているべきこと

四角形要素の内挿関数関係

この部分は,四角形要素による2次元応力解析の解説で述べた部分と同じである.

以下,$(a,b)$ は,$a=[−1,1]$,$b=[−1,1]$ の範囲を有する正規化パラメータ座標である.
また,$(x,y)$ は,実モデルの $x$ および $y$ 座標である.

内挿関数の数式表現

$$
\begin{align}
N_1=&\frac{1}{4}(1-a)(1-b) \qquad N_2=\frac{1}{4}(1+a)(1-b) \\
N_3=&\frac{1}{4}(1+a)(1+b) \qquad N_4=\frac{1}{4}(1-a)(1+b)
\end{align}
$$

内挿関数の微分(正規化パラメータでの微分)

$$
\begin{align}
&\cfrac{\partial N_1}{\partial a}=-\cfrac{1}{4}(1-b) &\cfrac{\partial N_2}{\partial a}=+\cfrac{1}{4}(1-b) & &\cfrac{\partial N_3}{\partial a}=+\cfrac{1}{4}(1+b) & &\cfrac{\partial N_4}{\partial a}=-\cfrac{1}{4}(1+b) \\
&\cfrac{\partial N_1}{\partial b}=-\cfrac{1}{4}(1-a) &\cfrac{\partial N_2}{\partial b}=-\cfrac{1}{4}(1+a) & &\cfrac{\partial N_3}{\partial b}=+\cfrac{1}{4}(1+a) & &\cfrac{\partial N_4}{\partial b}=+\cfrac{1}{4}(1-a)
\end{align}
$$

ヤコビ行列と行列式

ヤコビ行列 $[J]$ の要素を示す.ここに $x_{i,j,k,l}$ および $y_{i,j,k,l}$ は要素を構成する節点座標である.

$$
\begin{align}
&J_{11}=\frac{\partial N_1}{\partial a} x_i+\frac{\partial N_2}{\partial a} x_j+\frac{\partial N_3}{\partial a} x_k+\frac{\partial N_4}{\partial a} x_l \\
&J_{12}=\frac{\partial N_1}{\partial a} y_i+\frac{\partial N_2}{\partial a} y_j+\frac{\partial N_3}{\partial a} y_k+\frac{\partial N_4}{\partial a} y_l \\
&J_{21}=\frac{\partial N_1}{\partial b} x_i+\frac{\partial N_2}{\partial b} x_j+\frac{\partial N_3}{\partial b} x_k+\frac{\partial N_4}{\partial b} x_l \\
&J_{22}=\frac{\partial N_1}{\partial b} y_i+\frac{\partial N_2}{\partial b} y_j+\frac{\partial N_3}{\partial b} y_k+\frac{\partial N_4}{\partial b} y_l
\end{align}
$$

$$
\begin{equation}
\det(J)=J_{11}\cdot J_{22}-J_{12}\cdot J_{21}
\end{equation}
$$

内挿関数の微分(実座標での微分)

$$
\begin{equation}
\cfrac{\partial N_i}{\partial x}=\cfrac{1}{\det(J)}\left\{J_{22}\cfrac{\partial N_i}{\partial a}-J_{12}\cfrac{\partial N_i}{\partial b}\right\} \qquad
\cfrac{\partial N_i}{\partial y}=\cfrac{1}{\det(J)}\left\{-J_{21}\cfrac{\partial N_i}{\partial a}+J_{11}\cfrac{\partial N_i}{\partial b}\right\}
\end{equation}
$$

ガウス積分

積分点を 4 点 ($n=2$) とすれば,$a$,$b$ および重み $H$ の値は下表の通りであり,$a$,$b$ を変化させた 4 回の値の総和が積分の近似値となる.また,重みが $H=1$ であることから,計算は更に単純になる.

$$
\begin{equation}
\int_{-1}^{1}\int_{-1}^{1} f(a,b)\cdot da\cdot db
\doteqdot \sum_{i=1}^n\sum_{j=1}^n H_i\cdot H_j\cdot f(a_i,b_j)
= \sum_{kk=1}^4 f(a_{kk}, b_{kk})
\end{equation}
$$

$$
\begin{align}
& i & j & & a & & b & & H & & kk \\
& 1 & 1 & & -0.57735 02692 & & -0.57735 02692 & & 1.00000 00000 & & 1 \\
& 1 & 2 & & +0.57735 02692 & & -0.57735 02692 & & 1.00000 00000 & & 2 \\
& 2 & 1 & & +0.57735 02692 & & +0.57735 02692 & & 1.00000 00000 & & 3 \\
& 2 & 2 & & -0.57735 02692 & & +0.57735 02692 & & 1.00000 00000 & & 4
\end{align}
$$

浸透流解析のための知識

この部分は,三角形要素の浸透流解析の解説で述べた部分と同じである.

浸透流解析のための有限要素式

$$
\begin{align}
&[\boldsymbol{k}]\{\boldsymbol{\phi}\}=\{\boldsymbol{q}\} \\
&[\boldsymbol{k}]=\int_A \left(k_x\frac{\partial [\boldsymbol{N}]^T}{\partial x}\frac{\partial [\boldsymbol{N}]}{\partial x}+k_y\frac{\partial [\boldsymbol{N}]^T}{\partial y}\frac{\partial [\boldsymbol{N}]}{\partial y}\right)dA \\
&\{\boldsymbol{q}\}=-\int_s [\boldsymbol{N}]^T \left(v_x \frac{\partial x}{\partial n}+v_y \frac{\partial y}{\partial n}\right)ds
\end{align}
$$

$$
\begin{align}
&[\boldsymbol{k}] & &\text{:透水性マトリックス}\\
&\{\boldsymbol{\phi}\} & &\text{:節点全水頭ベクトル}\\
&\{\boldsymbol{q}\} & &\text{:節点流量ベクトル}
\end{align}
$$

なお $\{\boldsymbol{q}\}$ は流速ベクトルの境界長さでの積分であり,流量を示すことは容易に理解できる.

要素流速算定式

$x$ 方向および $y$ 方向の流速をそれぞれ $v_x$,$v_y$ とすれば,ダルシーの法則より

$$
\begin{equation}
v_x=-k_x\frac{\partial \phi}{\partial x} \qquad v_y=-k_x\frac{\partial \phi}{\partial y}
\end{equation}
$$

不飽和透水特性を示す式(Van Genuchten モデル)

$$
\begin{align}
&S_e=\left\{\cfrac{1}{1+(\alpha\cdot H_s)^n}\right\}^m \qquad\qquad n=\cfrac{1}{1-m} \qquad (0<m<1)\\
&K_r=(S_e)^{0.5}\cdot\left\{1-\left(1-S_e{}^{1/m}\right)^m\right\}^2 \qquad 0<S_e, K_r\leq 1\\
&K=K_r \cdot K_0
\end{align}
$$

$$
\begin{align}
&S_e & &\text{:有効飽和度}\\
&K_r & &\text{:比透水係数}\\
&H_S & &\text{:サクション(負圧が正)}\\
&\alpha & &\text{:不飽和特性を表すパラメータ(正の実数)}\\
&m & &\text{:不飽和特性を表すパラメータ(0を超え1未満の実数)}\\
&K & &\text{:不飽和透水係数}\\
&K_0 & &\text{:飽和透水係数}
\end{align}
$$

四角形要素による定式化

この部分は,四角形要素による浸透流解析を行うために必要な説明である.

透水性マトリックス

面積積分の積分点を4点,線積分の積分点を2点とすれば,重みは1となるため,それぞれの行列・ベクトルは積分点での計算値の総和として次のように示される.

$$
\begin{align}
[\boldsymbol{k}]=&\int_A \left(k_x\frac{\partial [\boldsymbol{N}]^T}{\partial x}\frac{\partial [\boldsymbol{N}]}{\partial x}+k_y\frac{\partial [\boldsymbol{N}]^T}{\partial y}\frac{\partial [\boldsymbol{N}]}{\partial y}\right)dA \\
=&\int_{-1}^{1}\int_{-1}^{1} \left(k_x\frac{\partial [\boldsymbol{N}]^T}{\partial x}\frac{\partial [\boldsymbol{N}]}{\partial x}+k_y\frac{\partial [\boldsymbol{N}]^T}{\partial y}\frac{\partial [\boldsymbol{N}]}{\partial y}\right)\cdot det(J)\cdot da\cdot db \\
=&\sum_{kk=1}^{4} \left[\left(k_x\frac{\partial [\boldsymbol{N}]^T}{\partial x}\frac{\partial [\boldsymbol{N}]}{\partial x}+k_y\frac{\partial [\boldsymbol{N}]^T}{\partial y}\frac{\partial [\boldsymbol{N}]}{\partial y}\right)\cdot det(J) \right]_ {kk}
\end{align}
$$

等方均質材料に対しては,$k_x=k_y=k$ として,

$$
\begin{equation}
[\boldsymbol{k}]=
k\cdot \sum_{kk=1}^{4} \left[\left(\frac{\partial [\boldsymbol{N}]^T}{\partial x}\frac{\partial [\boldsymbol{N}]}{\partial x}+\frac{\partial [\boldsymbol{N}]^T}{\partial y}\frac{\partial [\boldsymbol{N}]}{\partial y}\right)\cdot det(J) \right]_ {kk}
\end{equation}
$$

要素内平均流速

要素内任意点の全水頭 $\phi$ は内挿関数 $[\boldsymbol{N}]$ と節点全水頭ベクトル ${\boldsymbol{h}}$ を用いて
ダルシーの法則より要素内任意点の流速 $v_x$ および $v_y$ は

$$
\begin{align}
&v_x=-k_x\frac{\partial \phi}{\partial x} =-k_x\frac{\partial [\boldsymbol{N}]}{\partial x}\{\boldsymbol{\phi}\} \\
&v_y=-k_y\frac{\partial \phi}{\partial y} =-k_y\frac{\partial [\boldsymbol{N}]}{\partial y}\{\boldsymbol{\phi}\}
\end{align}
$$

要素内平均流速は,要素内の4つのGauss積分点の平均とすればよいため,以下の通り求められる.

$$
\begin{align}
&v_x=-\frac{k_x}{4}\sum_{kk=1}^{4}\left(\frac{\partial [\boldsymbol{N}]}{\partial x}\{\boldsymbol{\phi}\}\right)_ {kk} \\
&v_y=-\frac{k_y}{4}\sum_{kk=1}^{4}\left(\frac{\partial [\boldsymbol{N}]}{\partial y}\{\boldsymbol{\phi}\}\right)_ {kk}
\end{align}
$$

節点流量ベクトル

節点流量ベクトル ${\boldsymbol{q}}$ については,等価節点ベクトルとして直接入力すればよい.

プログラムの構想

プログラムの機能

  • プログラミング言語はPythonとする
  • 等方均質透水地盤の2次元飽和ー不飽和定常浸透流解析を行う.等方均質地盤であることに注意.
  • 要素は,1要素4節点4Gauss積分点の四角形要素とする.地盤の奥行きは入力値の単位に応じた単位奥行きとする.
  • 出力は,節点全水頭,節点圧力水頭,節点流量および要素平均流速とする
  • 境界条件として,不透水境界,全水頭指定境界,流量指定境界,浸潤面境界が扱える.
  • 浸潤面境界への流入はないものとする(雨等の影響は考えない)
  • 座標システムは,右方向を x 軸正方向,上方向あるいは標高方向を z 軸正方向とする.
  • 連立一次方程式は,numpy.linalg.solve(A,b) で解く
  • 入力データファイル,出力データファイルの文字区切りは空白とする
  • 入力データファイル名,出力データファイル名は,コマンドライン引数として入力する

扱える境界条件の説明

項目 説明
不透水境界 何も指定しなければ不透水境界となる
全水頭指定境界 既知全水頭を指定する節点数・節点番号・全水頭値を入力する.正の流量が節点への流入,負の流量が節点からの流出を示す.
流量指定境界 既知流量を指定する節点数・節点番号・流量値を入力する
浸潤面境界 浸潤面が現れる可能性のある領域の節点数・節点番号を入力する

プログラムの構成(プログラム名:py_fem_seep4.py)

三角形一次要素に対し,数値積分をおこなうため,5.および 6. が追加されている.
メインルーチンの内容は,三角形要素用を四角形要素用に書き換えただけで,同じである.

1. モジュール読み込み
2. 関数:入力データ書き出し (def PRINP_S4)
3. 関数:計算結果書き出し (def PROUT_S4)
4. 関数:不飽和透水特性計算 (def CAL_KRE_S4)
5. 関数:内挿関数の微分とヤコビ行列行列式 (def DNMAT_S4)
6. 関数:ガウス積分点指定 (def AB_S4)
7. 関数:要素透水マトリックス作成 (def SM_S4)
8. 関数:要素流速計算 (def CALV_S4)
9. メインルーチン  
    (1) 入出力ファイル名取得    
    (2) 基本数値読み込み  
    (3) 配列宣言  
    (4) 材料特性データ読み込み  
    (5) 要素構成節点番号および材料番号読み込み  
    (6) 節点座標および節点全水頭初期値読み込み  
    (7) 境界条件読み込み  
    (8) 入力データ書き出し    
    (9) 反復計算過程    
        a. 透水性マトリックス組み立て   
        b. 節点流量計算
        c. 作業用流量ベクトル更新    
        d. 境界条件処理(全水頭既知自由度の処理)  
        e. 連立一次方程式を解く(全水頭計算)  
        f. 境界条件処理(既知全水頭の組み込み) 
        g. 収束判定
    (10) 要素流速ベクトル計算  
    (11) 計算結果書き出し  

実行コマンドと入出力データ

解析実行コマンド

python3 py_fem_seep4.py inp.txt out.txt
変数 説明
py_fem_seep4.py Python による FEM プログラム名
inp.txt 入力データファイル名(空白区切りデータ)
out.txt 出力データファイル名(空白区切りデータ)

入力データファイルフォーマット

npoin  nele  nsec  koh  koq  kou  idan  # Basic values for analysis
Ak0  alpha  em                          # Material properties
    ..... (1 to nsec) .....             #
node-1  node-2  node-3  node-4  isec    # Element connectivity
    ..... (1 to nele) .....             #   (counter-clockwise order of node number)
x  z  hvec0                             # Coordinate, initial value of total head
    ..... (1 to npoin) .....            #
nokh  Hinp                              # Node with given total head, given total head
    ..... (1 to koh) .....              #     (omit data input if koh=0)
nokq  Qinp                              # Node with given duscharge, given discharge
    ..... (1 to koq) .....              #     (omit data input if koq=0)
noku                                    # Node with saturation point
    ..... (1 to kou) .....              #     (omit data input if kou=0)
変数 説明
npoin, nele, nsec 節点数,要素数,材料特性数
koh, koq, kou 全水頭指定節点数,流量指定節点数,浸潤面境界節点数
idan 解析断面(0:鉛直断面,1:水平断面)
Ak0 要素の飽和透水係数
alpha, em 要素の不飽和透水特性
node1, node2, node3, node-4, isec 節点1,節点2,節点3,節点4,材料特性番号
x, y, hvec0 節点x座標,節点y座標,全水頭初期値
nokh, Hinp 全水頭指定節点番号,節点全水頭値
nokq, Qinp 流量指定節点番号,節点流量値
noku 浸潤面境界節点番号

(コメント)

  • 変数 idan は,圧力水頭計算時に,鉛直断面では(全水頭マイナス位置水頭),水平断面では(イコール全水頭)で計算する必要があるため,この識別を行うために導入している.
  • 要素の不飽和透水特性 alpha および em は,飽和解析を行いたい場合は,負値(例えば alpha=-1, em=-1)を入力する.これにより,全要素が飽和状態であることを明確にしている.
  • 節点全水頭初期値 hvec0 は,全節点同一の値を入れて構わない.節点全水頭指定境界では,プログラムの中で指定値に入れ替えている.

出力データファイルフォーマット

npoin  nele  nsec   koh   koq   kou  idan
    (Each value of above)
sec  Ak0  alpha  em
    Ak0   : Permeability coefficient (saturated)
    alpha : Un-saturated parameter
    em    : Un-saturated parameter
    ..... (1 to nsec) .....
node  x  z  hvec  qvec  koh  koq  kou
    node : Node number
    x    : x-coordinate
    z    : z-coordinate
    hvec : Initial total head
    qvec : Initial discharge
    koh  : Index of node with given total head (0: free, 1: fixed)
    koq  : Index of node with given discharge  (0: free, 1: fixed)
    kou  : Index of node with saturation point (0: free, 1: fixed)
    ..... (1 to npoin) .....
node Hinp
    node : Node number
    Hinp : Given total head
    ..... (1 to koh) .....
node  Qinp
    node : Node number
    Qinp : Given discharge
    ..... (1 to koq) .....
elem  i  j  k  l  sec
    elem : Element number
    i    : Node number of start point
    j    : Node number of second point
    k    : Node number of third point
    l    : Node number of end point
    sec  : Material number
    ..... (1 to nele) .....
node  hvec  pvec  qvec  koh  koq  kou
    node : Node number
    hvec : Total head
    pvec : Pressure head (Total head minus elevation of the node)
    qvec : Discharge (+: inflow, -: outflow)
    koh  : Index of node with given total head (0: free, 1: fixed)
    koq  : Index of node with given discharge  (0: free, 1: fixed)
    kou  : Index of node with saturation point (0: free, 1: fixed)
    ..... (1 to npoin) .....
elem  vx  vz  vm  kr
    elem : Element number
    vx   : Flow velocity in x-direction
    vz   : Flow velocity in z-direction
    vm   : Composed flow velocity
    kr   : Parameter of saturation (1: saturated, less than 1: unsaturated)
    ..... (1 to nele) .....
Total inflow = (total inflow discharge: positive sign)
Total outflow= (total outflow discharge: negative sign)
Max.velocity in all area     =  (maximum velocity in all area)
Max.velocity in inflow area  =  (maximum velocity in inflow area)
Max.velocity in outflow area =  (maximum velocity in outflow area)
iii=(number of itaration)
icount=(converged number of degrees of freedom)
kop=(number of pressured nodes)
n=(total degrees of freedom)  time=(calculation time)
変数 説明
hvec, pvec, qvec 全水頭,圧力水頭,節点流量
vx, vzm vm 要素 x 方向流速,要素 y 方向流速,合成流速
kr 飽和度パラメータ(1:飽和,1未満:不飽和)

プログラム本体

プログラム本体を,「プログラムの構成」に沿って以下に示します.

1. モジュール読み込み

############################
# 2DFEM Seepage Flow Analysis (4-nodes)
############################
import numpy as np
import sys
import time

2. 関数:入力データ書き出し (def PRINP_S4)

############################
# Write input data
############################
def PRINP_S4(fnameW,npoin,nele,nsec,koh,koq,kou,idan,ae,x,node,nokh,nokq,noku,hvec,qvec):
    fout=open(fnameW,'w')
    # print out of input data
    print('{0:>5s} {1:>5s} {2:>5s} {3:>5s} {4:>5s} {5:>5s} {6:>5s}'
    .format('npoin','nele','nsec','koh','koq','kou','idan'),file=fout)
    print('{0:5d} {1:5d} {2:5d} {3:5d} {4:5d} {5:5d} {6:5d}'
    .format(npoin,nele,nsec,koh,koq,kou,idan),file=fout)
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s}'
    .format('sec','Ak0','alpha','em'),file=fout)
    for i in range(0,nsec):
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e}'
        .format(i+1,ae[0,i],ae[1,i],ae[2,i]),file=fout)
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s} {5:>5s} {6:>5s} {7:>5s}'
    .format('node','x','z','hvec','qvec','koh','koq','kou'),file=fout)
    for i in range(0,npoin):
        lp=i+1
        kh=0; kq=0; ku=0
        if 0<koh:
            for k in range(0,koh):
                if i==nokh[k]-1: kh=1
        if 0<koq:
            for k in range(0,koq):
                if i==nokq[k]-1: kq=1
        if 0<kou:
            for k in range(0,kou):
                if i==noku[k]-1: ku=1
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e} {5:5d} {6:5d} {7:5d}'
        .format(lp,x[0,i],x[1,i],hvec[i],qvec[i],kh,kq,ku),file=fout)
    if 0<koh:
        print('{0:>5s} {1:>5s}'.format('node','Hinp'),file=fout)
        for k in range(0,koh):
            print('{0:5d} {1:15.7e}'
            .format(nokh[k],Hinp[k]),file=fout)
    if 0<koq:
        print('{0:>5s} {1:>5s}'.format('node','Qinp'),file=fout)
        for k in range(0,koq):
            print('{0:5d} {1:15.7e}'
            .format(nokq[k],Qinp[k]),file=fout)
    print('{0:>5s} {1:>5s} {2:>5s} {3:>5s} {4:>5s} {5:>5s}'
    .format('elem','i','j','k','l','sec'),file=fout)
    for ne in range(0,nele):
        print('{0:5d} {1:5d} {2:5d} {3:5d} {4:5d} {5:5d}'
        .format(ne+1,node[0,ne],node[1,ne],node[2,ne],node[3,ne],node[4,ne]),file=fout)
    fout.close()

3. 関数:計算結果書き出し (def PROUT_S4)

############################
# Write calculation result
############################
def PROUT_S4(fnameW,npoin,nele,nsec,koh,koq,kou,idan,ae,x,node,nokh,nokq,noku,hvec,qvec,evx,evz,evm,ekr):
    fout=open(fnameW,'a')
    # total head, pressure head, discharge
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>5s} {5:>5s} {6:>5s}'.format('node','hvec','pvec','qvec','koh','koq','kou'),file=fout)
    for i in range(0,npoin):
        lp=i+1
        kh=0; kq=0; ku=0
        if 0<koh:
            for k in range(0,koh):
                if i==nokh[k]-1: kh=1
        if 0<koq:
            for k in range(0,koq):
                if i==nokq[k]-1: kq=1
        if 0<kou:
            for k in range(0,kou):
                if i==noku[k]-1: ku=1
        if idan==0: pp=hvec[i]-x[1,i]
        if idan==1: pp=hvec[i]
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:5d} {5:5d} {6:5d}'.format(lp,hvec[i],pp,qvec[i],kh,kq,ku),file=fout)
    # velocity
    print('{0:>5s} {1:>15s} {2:>15s} {3:>15s} {4:>15s}'
    .format('elem','vx','vz','vm','kr'),file=fout)
    for ne in range(0,nele):
        print('{0:5d} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e}'
        .format(ne+1,evx[ne],evz[ne],evm[ne],ekr[ne]),file=fout)
    fout.close()

4. 関数:不飽和透水特性計算 (def CAL_KRE_S4)

まず,要素飽和していることを条件に比透水係数を1にセット(Akre=1.0
次に,不飽和透水特性パラメータが,ゼロ以上であること,サクション(負圧)が生じていることから要素は不飽和であると判断し,飽和度(se)および比透水係数(Akre)を算定している.
戻り値は,比透水係数 Akre のみである.これを要素の飽和透水係数に乗じることにより,飽和ー不飽和透水係数を算定している.

############################
# Calculation of degree of saturation
############################
def CAL_KRE_S4(pp,alpha,em):
    Akre=1.0
    sucp=-pp
    if 0<alpha and 0<em:
        if 0.0<sucp:
            en=1.0/(1.0-em)
            Se=(1.0/(1.0+(alpha*sucp)**en))**em
            Akre=np.sqrt(Se)*(1.0-(1.0-Se**(1.0/em))**em)**2
    return Akre

5. 関数:内挿関数の微分とヤコビ行列行列式 (def DNMAT_S4)

############################
# Differential of shape function [N]
############################
def DNMAT_S4(ne,node,x,a,b):
    dnx=np.zeros(4,dtype=np.float64)
    dnz=np.zeros(4,dtype=np.float64)
    i=node[0,ne]-1
    j=node[1,ne]-1
    k=node[2,ne]-1
    l=node[3,ne]-1
    #[dN/dadN/db]
    dn1a=-0.25*(1.0-b)
    dn2a= 0.25*(1.0-b)
    dn3a= 0.25*(1.0+b)
    dn4a=-0.25*(1.0+b)
    dn1b=-0.25*(1.0-a)
    dn2b=-0.25*(1.0+a)
    dn3b= 0.25*(1.0+a)
    dn4b= 0.25*(1.0-a)
    #Jacobi matrix and det(J)
    J11=dn1a*x[0,i]+dn2a*x[0,j]+dn3a*x[0,k]+dn4a*x[0,l]
    J12=dn1a*x[1,i]+dn2a*x[1,j]+dn3a*x[1,k]+dn4a*x[1,l]
    J21=dn1b*x[0,i]+dn2b*x[0,j]+dn3b*x[0,k]+dn4b*x[0,l]
    J22=dn1b*x[1,i]+dn2b*x[1,j]+dn3b*x[1,k]+dn4b*x[1,l]
    detJ=J11*J22-J12*J21
    #[dN/dx][dN/dy]
    dnx[0]=( J22*dn1a-J12*dn1b)/detJ
    dnx[1]=( J22*dn2a-J12*dn2b)/detJ
    dnx[2]=( J22*dn3a-J12*dn3b)/detJ
    dnx[3]=( J22*dn4a-J12*dn4b)/detJ
    dnz[0]=(-J21*dn1a+J11*dn1b)/detJ
    dnz[1]=(-J21*dn2a+J11*dn2b)/detJ
    dnz[2]=(-J21*dn3a+J11*dn3b)/detJ
    dnz[3]=(-J21*dn4a+J11*dn4b)/detJ
    return dnx,dnz,detJ

6. 関数:ガウス積分点指定 (def AB_S4)

############################
# Gauss points
############################
def AB_S4(kk):
    if kk==0:
        a=-0.5773502692
        b=-0.5773502692
    if kk==1:
        a= 0.5773502692
        b=-0.5773502692
    if kk==2:
        a= 0.5773502692
        b= 0.5773502692
    if kk==3:
        a=-0.5773502692
        b= 0.5773502692
    return a,b

7. 関数:要素透水マトリックス作成 (def SM_S4)

要素の比透水係数は,各節点での圧力水頭から計算した比透水係数の平均値としている.

############################
# Element permeability matrix
############################
def SM_S4(ne,node,x,ae,hvec):
    sm=np.zeros([4,4],dtype=np.float64)
    k1=node[0,ne]-1
    k2=node[1,ne]-1
    k3=node[2,ne]-1
    k4=node[3,ne]-1
    m=node[4,ne]-1
    Ak0  =ae[0,m]
    alpha=ae[1,m]
    em   =ae[2,m]
    p1=hvec[k1]-x[1,k1]
    p2=hvec[k2]-x[1,k2]
    p3=hvec[k3]-x[1,k3]
    p4=hvec[k4]-x[1,k4]
    for kk in range(0,4):
        a,b=AB_S4(kk)
        dnx,dnz,detJ=DNMAT_S4(ne,node,x,a,b)
        nm1=0.25*(1.0-a)*(1.0-b)
        nm2=0.25*(1.0+a)*(1.0-b)
        nm3=0.25*(1.0+a)*(1.0+b)
        nm4=0.25*(1.0-a)*(1.0+b)
        pp=nm1*p1+nm2*p2+nm3*p3+nm4*p4
        Akre=CAL_KRE_S4(pp,alpha,em)
        wxz=np.zeros([4,4],dtype=np.float64)
        for i in range(0,4):
            for j in range(0,4):
                wxz[i,j]=wxz[i,j]+dnx[i]*dnx[j]+dnz[i]*dnz[j]
        sm=sm+Ak0*Akre*wxz*detJ
        del wxz
    return sm

8. 関数:要素流速計算 (def CALV_S4)

透水マトリックス作成時と同様の扱いで,要素の比透水係数は,各節点での圧力水頭から計算した比透水係数の平均値としている.

############################
# Velocity calculation
############################
def CALV_S4(ne,node,x,hvec,ae):
    k1=node[0,ne]-1
    k2=node[1,ne]-1
    k3=node[2,ne]-1
    k4=node[3,ne]-1
    m=node[4,ne]-1
    Ak0  =ae[0,m]
    alpha=ae[1,m]
    em   =ae[2,m]
    p1=hvec[k1]-x[1,k1]
    p2=hvec[k2]-x[1,k2]
    p3=hvec[k3]-x[1,k3]
    p4=hvec[k4]-x[1,k4]
    vx=0.0
    vz=0.0
    Akr=0.0
    for kk in range(0,4):
        a,b=AB_S4(kk)
        dnx,dnz,detJ=DNMAT_S4(ne,node,x,a,b)
        nm1=0.25*(1.0-a)*(1.0-b)
        nm2=0.25*(1.0+a)*(1.0-b)
        nm3=0.25*(1.0+a)*(1.0+b)
        nm4=0.25*(1.0-a)*(1.0+b)
        pp=nm1*p1+nm2*p2+nm3*p3+nm4*p4
        Akre=CAL_KRE_S4(pp,alpha,em)
        Akr=Akr+Akre
        wvx=dnx[0]*hvec[k1]+dnx[1]*hvec[k2]+dnx[2]*hvec[k3]+dnx[3]*hvec[k4]
        wvz=dnz[0]*hvec[k1]+dnz[1]*hvec[k2]+dnz[2]*hvec[k3]+dnz[3]*hvec[k4]
        vx=vx+wvx
        vz=vz+wvz
    Akr=0.25*Akr
    vx=-0.25*Ak0*Akr*vx
    vz=-0.25*Ak0*Akr*vz
    vm=np.sqrt(vx**2+vz**2)
    return Akr,vx,vz,vm

9. メインルーチン

(1) 入出力ファイル名取得

############################
# Main routine
############################
start=time.time()
args = sys.argv
fnameR=args[1] # input data file
fnameW=args[2] # output data file

(2) 基本数値読み込み

f=open(fnameR,'r')
text=f.readline()
text=text.strip()
text=text.split()
npoin =int(text[0]) # Number of nodes
nele  =int(text[1]) # Number of elements
nsec  =int(text[2]) # Number of sections
koh   =int(text[3]) # Number of nodes with given total head
koq   =int(text[4]) # Number of nodes with given discharge
kou   =int(text[5]) # Number of nodes with saturation point
idan  =int(text[6]) # Section (0: vertical, 1: horizontal)
nod=4               # Number of nodes per element
nfree=1             # Degree of freedom per node
n=nfree*npoin

(3) 配列宣言

x     =np.zeros([2,npoin],dtype=np.float64)     # Coordinates of nodes
ae    =np.zeros([3,nsec],dtype=np.float64)      # Section characteristics
node  =np.zeros([nod+1,nele],dtype=np.int)      # Node-element relationship
hvec  =np.zeros(nfree*npoin,dtype=np.float64)   # Total head vector
qvec  =np.zeros(nfree*npoin,dtype=np.float64)   # Discharge vector
nokh  =np.zeros(koh,dtype=np.int)               # Nodes with given total head
nokq  =np.zeros(koq,dtype=np.int)               # Nodes with given discharge
noku  =np.zeros(kou,dtype=np.int)               # Nodes with saturation point
nokp  =np.zeros(kou,dtype=np.int)               # Nodes with saturation point
Hinp  =np.zeros(koh,dtype=np.float64)           # Total head at specified nodes
Qinp  =np.zeros(koq,dtype=np.float64)           # Discharge at specified nodes (+: inflow, -: outflow)
ir    =np.zeros(nod*nfree,dtype=np.int)         # Work vector for matrix assembly
ekr   =np.zeros(nele,dtype=np.float64)          # Memory of rerative mermeability
evx   =np.zeros(nele,dtype=np.float64)          # Memory of element velocity in x-direction
evz   =np.zeros(nele,dtype=np.float64)          # Memory of element velocity in z-direction
evm   =np.zeros(nele,dtype=np.float64)          # Memory of norm of element velocity

(4) 材料特性データ読み込み

# section characteristics
for i in range(0,nsec):
    text=f.readline()
    text=text.strip()
    text=text.split()
    ae[0,i]=float(text[0]) #Ak0   : Permeability coefficient (saturated)
    ae[1,i]=float(text[1]) #alpha : Unsaturated property
    ae[2,i]=float(text[2]) #em    : Unsaturated property

(5) 要素構成節点番号および材料番号読み込み

# element-node
# element-node
for i in range(0,nele):
    text=f.readline()
    text=text.strip()
    text=text.split()
    node[0,i]=int(text[0]) #node_1
    node[1,i]=int(text[1]) #node_2
    node[2,i]=int(text[2]) #node_3
    node[3,i]=int(text[3]) #node_4
    node[4,i]=int(text[4]) #section characteristic number

(6) 節点座標および節点全水頭初期値読み込み

# node coordinates
for i in range(0,npoin):
    text=f.readline()
    text=text.strip()
    text=text.split()
    x[0,i] =float(text[0]) # x-coordinate
    x[1,i] =float(text[1]) # y-coordinate
    hvec[i]=float(text[2]) # Initial total head

(7) 境界条件読み込み

全水頭指定境界の節点番号・全水頭値,流量指定境界の節点番号・流量値,浸潤面境界の節点番号の順で読み込んでいる.

# boundary conditions
if 0<koh:
    for i in range(0,koh):
        text=f.readline()
        text=text.strip()
        text=text.split()
        nokh[i]=int(text[0])
        Hinp[i]=float(text[1])
if 0<koq:
    for i in range(0,koq):
        text=f.readline()
        text=text.strip()
        text=text.split()
        nokq[i]=int(text[0])
        Qinp[i]=float(text[1])
if 0<kou:
    for i in range(0,kou):
        text=f.readline()
        text=text.strip()
        text=text.split()
        noku[i]=int(text[0])
f.close()

(8) 入力データ書き出し

全水頭ベクトル hvec および流量ベクトル qvec は初期化により,全要素ゼロとなっているが,入力データ書き出し前に,全水頭指定節点および流量指定節点の全水頭および流量値を指定値(Hinp, Qinp)に更新している.

if 0<koh:
    for k in range(0,koh):
        hvec[nokh[k]-1]=Hinp[k] # Set given total head
if 0<koq:
    for k in range(0,koq):
        qvec[nokq[k]-1]=Qinp[k] # Set given discharge
PRINP_S4(fnameW,npoin,nele,nsec,koh,koq,kou,idan,ae,x,node,nokh,nokq,noku,hvec,qvec)

(9) 反復計算過程

反復計算に入る前に,収束判定に用いる全水頭ベクトル h1vec の定義,および反復計算の繰り返し上限 imax の指定を行っている.

h1vec=hvec # Setting initial total head
itmax=100  # Max number of iterative calculation
for iii in range(0,itmax): # Iterative calculation

a. 透水性マトリックス組み立て

仮定した全水頭ベクトル h1vec を用い,要素透水性マトリックスを作成する.

    # Assembly of permeability matrix, [k] is calculated using {h1}
    gk=np.zeros([n,n],dtype=np.float64)  # Global stiffness matrix
    for ne in range(0,nele):
        sm=SM_S4(ne,node,x,ae,h1vec)
        i=node[0,ne]-1
        j=node[1,ne]-1
        k=node[2,ne]-1
        l=node[3,ne]-1
        ir[3]=l; ir[2]=k; ir[1]=j; ir[0]=i
        for i in range(0,nod*nfree):
            it=ir[i]
            for j in range(0,nod*nfree):
                jt=ir[j]
                gk[it,jt]=gk[it,jt]+sm[i,j]

b. 節点流量計算

組み立てた全体透水性マトリックス gk と仮定された全水頭ベクトル h1vec から,節点流量ベクトル qvec を算定する.

    # calculation of nodal discharge {q}=[k]{h1}
    qvec=np.dot(gk,h1vec)

c. 作業用流量ベクトル更新

連立一次方程式作成に先立ち,流量ベクトルを確からしい値に修正するため,以下の操作を行う.

  • 作業用流量ベクトル wvec を作成・初期化(全要素ゼロセット)する.
  • 全水頭指定境界での流量は,先に計算した節点流量ベクトルの値に更新する
  • 流量指定境界での流量は,入力データで指定した指定値に更新する
  • 浸潤面境界での流量は,先に計算した節点流量ベクトルの値に更新する.ただし値が正(流入)の場合はゼロとする.
    #Update of nodal discharge {q}->{qw}
    # (1) for node with given total head
    # (2) for node with given discharge
    # (3) for node with positive discharge in the seapage point boundaris
    #     (In-flow is not allowed at seepage point. In-flow means positive discharge.)
    #     (If discharge has positive value at seepage point, the discharge is set to zero)
    wvec=np.zeros(nfree*npoin,dtype=np.float64)   # Work vector for discharge
    if 0<koh:
        for k in range(0,koh):
            wvec[nokh[k]-1]=qvec[nokh[k]-1]
    if 0<koq:
        for k in range(0,koq):
            wvec[nokq[k]-1]=Qinp[k]
    if 0<kou:
        for k in range(0,kou):
            wvec[noku[k]-1]=qvec[noku[k]-1]
            if 0.0<wvec[noku[k]-1]: wvec[noku[k]-1]=0.0

d. 境界条件処理(全水頭既知自由度の処理)

浸潤面境界節点は,大気に接しているため,圧力水頭は少なくともゼロ以下である必要がある.
これを境界条件処理の条件に取り入れ,以下の操作を行う.

  • 浸潤面境界節点のうち,圧力水頭が正,すなわち全水頭が位置水頭より大きい(x[1,i] < h1vec[i]x[1,i] は節点 i の z 座標)節点の個数(kop)と番号(nokp)を記憶する.
  • 境界条件処理において,全水頭指定節点の流量ベクトルの補正を指定全水頭を用いて行う.( wvec=wvec-Hinp[k]*gk[:,iz]
  • 先に記憶した,浸潤面境界節点のうち圧力水頭が正となっている節点において,流量ベクトルの補正を,全水頭を位置水頭に等しいとして行う.( wvec=wvec-x[1,iz]*gk[:,iz]
    # Memory of node numbers with positive pressure head in seepage points
    kop=0
    if 0<kou:
        for k in range(0,kou):
            i=noku[k]-1
            if x[1,i]<h1vec[i]: # if h1[i] is greater than z[i] (exclude equal to)
                kop=kop+1
                nokp[kop-1]=i+1
    # Deformation of FE formula considering boundary conditions
    # (1) set given total head
    # (2) set zero pressure head at the node with positive presure head in seepage points
    #     (Total head = position head)
    # treatment of boundary conditions
    if 0<koh:
        for k in range(0,koh):
            iz=nokh[k]-1
            wvec[iz]=0.0
    if 0<kop:
        for k in range(0,kop):
            iz=nokp[k]-1
            wvec[iz]=0.0
    if 0<koh:
        for k in range(0,koh):
            iz=nokh[k]-1
            wvec=wvec-Hinp[k]*gk[:,iz]
            gk[:,iz]=0.0
            gk[iz,iz]=1.0
    if 0<kop:
        for k in range(0,kop):
            iz=nokp[k]-1
            wvec=wvec-x[1,iz]*gk[:,iz]
            gk[:,iz]=0.0
            gk[iz,iz]=1.0

e. 連立一次方程式を解く(全水頭計算)

    # solution of simultaneous linear equations
    hvec = np.linalg.solve(gk, wvec)
    del gk,wvec

f. 境界条件処理(既知全水頭の組み込み)

    # Recovery of given or known total head
    if 0<koh:
        for k in range(0,koh):
            hvec[nokh[k]-1]=Hinp[k]
    if 0<kop:
        for k in range(0,kop):
            hvec[nokp[k]-1]=x[1,nokp[k]-1]
    h2vec=hvec

g. 収束判定

前回の反復過程で得られた全水頭ベクトル h1vec と今回の反復過程で得られた全水頭ベクトル h2vec の差分により収束判定を行う.
前回と今回の全水頭ベクトルの差分が,全節点で $1 \times 10^{-6}$ 未満となった時,収束とみなし,反復計算のループを抜ける.

収束しなければ,全水頭ベクトルを更新( h1vec = h2vec )して,反復計算を繰り返す.

    # Judgement of convergence |h2-h1|<1e-6
    icount=0
    for i in range(0,n):
        if np.abs(h2vec[i]-h1vec[i])<1.0e-6: icount=icount+1
    print('nnn={0}  icount={1}  kop={2}'.format(iii,icount,kop))
    if icount==n: break
    # Update of total head {h1}={h2}
    h1vec=h2vec

(10) 要素流速ベクトル計算

このプログラムの前提では,解析領域の全流入量と全流出量が等しくなる必要があるため,これを確認するため,全流入量 QTi と全流出量 QTo を求めている.
その後の for ne in range(0,nele): のループ内が,要素流速ベクトルの計算となっている.出力の一部とするため,流入領域と流筒領域での流速の最大値検索も同時に行っている.

# Calculation of total in-flow and out-flow
QTi=0.0  #Total in-flow
QTo=0.0  #Total out-flow
for i in range(0,n):
    if 0.0<qvec[i]: QTi=QTi+qvec[i]
    if qvec[i]<0.0: QTo=QTo+qvec[i]
# Calculation of element velocity & to search max.velocity
vmmax=0.0  #Max.velocity
vimax=0.0  #Max.velocity on in-flow area
vomax=0.0  #max.velocity on out-flow area
kmmax=0    #Element number with max.velocity
kimax=0    #Element number with max.velocity on in-flow area
komax=0    #Element number with max.velocity on out-flow area
for ne in range(0,nele):
    Akr,vx,vz,vm=CALV_S4(ne,node,x,hvec,ae)
    evx[ne]=vx
    evz[ne]=vz
    evm[ne]=vm
    ekr[ne]=Akr
    if vmmax<vm:
        kmmax=ne+1
        vmmax=vm
    i=node[0,ne]-1
    j=node[1,ne]-1
    k=node[2,ne]-1
    l=node[3,ne]-1
    if 0.0<(qvec[i]+qvec[j]+qvec[k]+qvec[l]):
        if vimax<vm:
            kimax=ne+1
            vimax=vm
    if (qvec[i]+qvec[j]+qvec[k]+qvec[l])<0.0:
        if vomax<vm:
            komax=ne+1
            vomax=vm

(11) 計算結果書き出し

PROUT_S4(fnameW,npoin,nele,nsec,koh,koq,kou,idan,ae,x,node,nokh,nokq,noku,hvec,qvec,evx,evz,evm,ekr)

if idan==0:
    pmax=np.max(hvec-x[1,:])
    pmin=np.min(hvec-x[1,:])
if idan==1:
    pmax=np.max(hvec)
    pmin=np.min(hvec)

dtime=time.time()-start
print('n={0}  time={1:.3f}'.format(n,dtime)+' sec')
fout=open(fnameW,'a')
print('Total inflow ={0:15.7e}'.format(QTi),file=fout)
print('Total outflow={0:15.7e}'.format(QTo),file=fout)
print('Max.velocity in all area     ={0:15.7e} (ne={1})'.format(vmmax,kmmax),file=fout)
print('Max.velocity in inflow area  ={0:15.7e} (ne={1})'.format(vimax,kimax),file=fout)
print('Max.velocity in outflow area ={0:15.7e} (ne={1})'.format(vomax,komax),file=fout)
print('Max.Pressure head ={0:15.7e}'.format(pmax),file=fout)
print('Min.Pressure head ={0:15.7e}'.format(pmin),file=fout)
print('iii={0}  icount={1}  kop={2}'.format(iii+1,icount,kop),file=fout)
print('n={0}  time={1:.3f}'.format(n,dtime)+' sec',file=fout)
fout.close()

以 上

  • プログラミングのために知っているべきこと
  • 四角形要素による定式化
  • プログラムの構想
  • 実行コマンドと入出力データ
  • 出力データファイルフォーマット
  • プログラム本体