Corrected the program due to some program mistakes
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.
`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.Saturated steady-state osmotic flow analysis is a problem that solves the following simultaneous linear equations.
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.
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.
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.
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.
If you put the program on it, it will be long, so put a link to your own HP.
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 td> | Number of nodes, number of elements, number of material characteristics td> tr> |
koh, koq, kou td> | Number of head designated nodes, flow rate specified nodes, infiltration boundary nodes td> tr> |
idan td> | Analysis cross section (0: vertical cross section, 1: horizontal cross section) td> tr> |
Ak0 td> | Saturated hydraulic conductivity of the element td> tr> |
alpa, em td> | Unsaturated permeability of elements td> tr> |
node-1, node-2, node-3, node-4, isec td> | Node numbers that make up the element, material property number of the element td> tr> |
x, z, hvec0 td> | x and z coordinates of the node, initial head value for convergence calculation td> tr> |
nokh, Hinp td> | Head designation node number and head value td> tr> |
nokq, Qinp td> | Flow rate specified node number and flow value td> tr> |
noku td> | Infiltration boundary node number td> tr> |
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 td> | Node number, total head, pressure head, flow rate td> tr> |
elem, vx, vz, vm td> | Element x / z direction flow velocity and synthetic flow velocity td> tr> |
Kr td> | Element unsaturated permeability characteristics (1: saturated, 1>: unsaturated) td> tr> |
Total inflow, Total outflow td> | Total inflow and outflow to the model area td> tr> |
iii, icount td> | Number of convergence calculations, number of freedoms that satisfy the calculation completion condition td> tr> |
kop td> | The number of infiltration surface boundary nodes where the pressure head is 0 td> tr> |
n, time td> | total degrees of freedom, calculation time td> tr> |
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
Recommended Posts