2차원 포화 불포화 삼투류 분석

(2017.11.19 수정)

일부 프로그램 오류가 있어서 프로그램을 수정했습니다

개요

대강대강 한 감상

2차원 FEM 응력 분석 프로그램을 만들었기 때문에 겸사겸사 2차원 포화-불포화 침투 분석 프로그램을 만들었다.
원소는 2차원 응력 분석과 마찬가지로 원소 4개 노드 4고스 적분점을 가진 동위소 파라미터 원소다.
실제 사용 여부는 노드수 1317, 소자수 1227, 수렴 계산 횟수 10회, 계산 시간 3.7초로 이런 규모의 문제라면 특별한 압력이 없다.
Numpy의 행렬 연산과 해연립 방정식을 결합하는 프로그램은 우수하죠.
절차상 주의사항은 다음과 같다.

모든 어레이가 Numpy로 정렬됨

for의 사용은 최대한 피하고 Numpy의 행렬 연산을 사용합니다. 그러나 1차원 벡터의 적 등으로 인해 예상치 못한 결과가 발생할 수 있으므로 주의해야 합니다

연립 방정식을 한 번 풀면 충분하다x=numpy.linalg.solve(A, b). 경계 조건 처리로 행렬이 비대칭적이기 때문에 스키 분해 등은 도저히 사용할 수 없다

상기 조작을 집행함으로써 프로그램 주체의 전망이 매우 좋아졌다

입력과 출력이 엉망진창인데 피할 수 없는 건가요?

포화 삼투류 분석

포화 안정적 침류 분석은 다음과 같은 일차 방정식을 결합하는 문제를 해결하는 것이다.
$$
[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}
$$
여기서 $\{q\} 노드 유량 벡터, $\{h\} 전수두 벡터, $[k] 달러는 투명 행렬로 다음과 같은 특성을 가지고 있습니다.

투명 행렬 $[k]달러는 상수,

노드 트래픽 벡터 $\{q\}의 알려진 분량에 대응하는 온수두 벡터 $\{h\} 달러 미지수

노드 유량 벡터 $\{q\}달러의 미지의 분량과 비슷한 전수두 벡터 $\{h\}달러는 미지수

따라서 이미 알고 있는 수와 미지수 간의 관계를 처리한 후 연립 방정식을 구해하여 해를 얻는다.

(한 방정식을 연립하는 이미 알고 있는 수와 미지수 간의 관계를 해결하는 처리)

$q_i$,$q_j$데이터 및 $h 제거i$,$h_알려진 정보, $hi$,$h_j$를 제외한 모든 수두 및 $qi$,$q_만약 j$가 알 수 없다면 투명 행렬의 $k{i}달러 및 $k}$1, 이외의 $i열과 $j$j열 요소를 0으로 처리하면 투명 행렬의 $i열과 $j$j열의 효과를 오른쪽으로 옮깁니다.
$$
\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}
$$

포화-불포화 삼투류 분석

포화-불포화 삼투류 분석에서의 연립 방정식은 다음과 같은 특성을 가지고 있다.

경계 조건은 유량 기지, 전수두 기지, 침윤면의 경계가 나타나는 세 가지를 고려해야 한다.

침윤면이 나타나는 경계를 제외하고는 유량이나 전수두의 어느 하나도 이미 알고 있다.

침윤면이 나타나는 경계에서 유량과 전수두는 미지수이다

그러나 침윤면이 나타나는 경계는 기본적으로 대기와 접촉하는 경계이기 때문에 빗물의 유입 등을 고려하지 않으면 유량은 유출(절차상 0 이하)에 불과하다.
압력수두도 0 이하의 제한이 있다

포화구역 중의 투수 행렬의 성분은 상수이지만 불포화구역 중의 투수 행렬의 성분은 흡기압력(음압력수두)의 함수이다
따라서 수렴 계산을 통해 미지수를 정해야 한다.
수렴 계산 방법은 전체 노드에 전수두의 초기 값을 제공하여 얻은 해를 다음 수렴 계산의 입력 값의 순서로 대입하는 것이다. 이런 방법은 모든 노드의 수두를 가설하여 계산 절차를 간소화할 수 있다.
전수두의 초기값은 전수두가 이미 알고 있는 경계를 제외하고 해결 방안의 정밀도와 수렴성을 고려한 토대에서 포화구역의 최소 전수두를 주는 것이 효과적인 것 같다.

불포화 투수 특성(van Genuuchten 모델)

지반의 불포화 투수 특성 모델도 표에 포화도-흡압 관계, 포화도-불포화 투수 계수비를 입력하는 방법이 있지만 분석을 위해 연속 함수 처리를 사용하여 프로그램을 간단하게 할 수 있도록 van Genuuchten 모델을 사용했다.
$$
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

절차.

프로그램이 오래되어 자신의 홈페이지 링크를 붙였다.

포화-불포화 침투 분석 프로그램 (py fem seep 4.py)

간이 압력 조절기 그래픽 프로그램 (py seep4 conte.py)

영역형 잠금 필터 해석용 입력 데이터(inp zone 4.txt)

입력 데이터 형식

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
노드 수, 요소 수, 재료 특성 수
koh, koq, kou
수두 지정 노드수, 유량 지정 노드수, 침윤 경계 노드수
idan
분석 단면(0: 수직 단면, 1: 수평 단면)
Ak0
원소의 포화 투수 계수
alpa, em
원소의 불포화 투수 특성
node-1, node-2, node-3, node-4, isec
구성 요소의 노드 번호, 요소의 재료 특성 번호
x, z, hvec0
노드의 x와z 좌표, 수렴 계산용 초기 수두값
nokh, Hinp
물머리 지정 노드 번호 및 물머리 값
nokq, Qinp
트래픽 지정 노드 번호 및 트래픽 값
noku
침윤 경계 노드 번호

출력 데이터 형식

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
노드번호,온수두,압력수두,유량
elem, vx, vz, vm
원소의 x·z방향 유속 및 합성 유속
Kr
원소의 불포화 투수 특성 (1:포화, 1&:불포화)
Total inflow, Total outflow
모델 구역의 전유입량 및 전유출량
iii, icount
계산 횟수를 수렴하여 계산 완성 조건의 자유도수를 만족시키다
kop
침윤면 경계 노드 중 압력수두 0의 노드수
n, time
전자유도, 시간 계산

수출 사례

간단한 변환기 제작 프로그램을 통해 압력수두 0, 5, 10, 15미터의 변환기를 잡아당기는 그림을 시도해 본다.
상류는 높이 18m, 하류는 높이 4m의 수심을 부여한다.
컴퓨터 그래프 자체는 비록 문외한 수준이지만 계산은 매우 순조로운 것 같다.
또 압력 헤드와 댐 규격을 빨간 글자로 표시한 것은 맥의 프리뷰로 추가로 작성했다.
이곳의 작도는 GMT(Generic Mapping Tools)의 변환 기능에 맡겨 그림을 만들어 보고서로 활용하는 것이 좋다.
참고로 실제 댐에서는 심지와 여과기의 경계 부근에서 압력수두가 낮아지고 분석결과처럼 심지에서 압력이 낮아지지 않는다는 것은 잘 알려져 있다.
결론적으로 이것은 해석 사례이니 이해해 주십시오.

이상