Zweidimensionale Analyse des gesättigten und ungesättigten Permeationsflusses mit Python

(Korrigiert am 19. November 2017)

Das Programm wurde aufgrund einiger Programmierfehler korrigiert

Überblick

Grober Eindruck

Da ich ein zweidimensionales FEM-Spannungsanalyseprogramm erstellt habe, habe ich auch ein zweidimensionales Programm zur Analyse des gesättigten und ungesättigten Permeationsflusses erstellt. Die Elemente sind isoparametrische Elemente mit 1 Element, 4 Knoten und 4 Gauß-Integrationspunkten, wie im Fall der 2D-Spannungsanalyse.

Eigentlich kann es verwendet werden, aber die Anzahl der Knoten beträgt 1317, die Anzahl der Elemente beträgt 1227, die Anzahl der Konvergenzberechnungen beträgt 10 und die Berechnungszeit beträgt 3,7 Sekunden. Wenn das Problem in dieser Größenordnung liegt, ist keine besondere Belastung zu spüren. Die Matrixoperationen und Routinen von Numpy zum Lösen simultaner linearer Gleichungen sind wahrscheinlich ausgezeichnet.

Folgende Punkte sind beim Programmieren zu beachten.

Alle Arrays sind Numpy-Arrays. Es wird keine Liste verwendet. Vermeiden Sie es, + so oft wie möglich zu verwenden, und verwenden Sie Numpy-Matrix-Operationen. Beachten Sie jedoch, dass aufgrund des Produkts eindimensionaler Vektoren unerwartete Ergebnisse auftreten können.
`x = numpy.linalg.solve (A, b)` reicht aus, um simultane lineare Gleichungen zu lösen. Da die Matrix aufgrund der Randbedingungsverarbeitung asymmetrisch wird, kann die Choleskey-Zerlegung in keiner Weise verwendet werden.
Wenn Sie dies tun, werden die Aussichten des Programms selbst sehr gut sein.
Der Eingabe- / Ausgabeabschnitt wird durcheinander gebracht. Ist dies unvermeidlich?

Analyse des gesättigten Permeationsflusses

Die Analyse des gesättigten stetigen Permeationsflusses ist ein Problem, das die folgenden simultanen linearen Gleichungen löst.

[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}

Hier ist $ \ {q \} $ der Knotenflussvektor, $ \ {h \} $ der Gesamtkopfvektor und $ [k] $ die Wasserdurchlässigkeitsmatrix, die die folgenden Eigenschaften aufweist.

Wasserpermeationsmatrix $ [k] $ ist eine Konstante,
Der Gesamtkopfvektor $ \ {h \} $, der der bekannten Komponente des Knotenflussvektors $ \ {q \} $ entspricht, ist unbekannt
Der Gesamtkopfvektor $ \ {h \} $, der der unbekannten Komponente des Knotenflussvektors $ \ {q \} $ entspricht, ist eine bekannte Zahl

Daher kann die Lösung erhalten werden, indem die simultanen Gleichungen einmal nach Verarbeitung der Beziehung zwischen der bekannten Zahl und der unbekannten Zahl gelöst werden.

(Verarbeitung der Beziehung zwischen bekannten und unbekannten Zahlen zur Lösung simultaner linearer Gleichungen)

Durchflussrate ohne $ q_i $, $ q_j $ und bekannte $ h_i $, $ h_j $, alle Köpfe ohne $ h_i $, $ h_j $ und $ q_i $, $ q_j $ in der Wasserdurchlässigkeitsmatrix Der Prozess wird so ausgeführt, dass k_ {ii} $ und $ k_ {jj} $ auf $ 1 $ gesetzt werden und die anderen Elemente der Spalte $ i $ und $ j $ auf Null gesetzt werden. Außerdem werden die Auswirkungen der Spalten $ i $ und $ j $ in der Wasserdurchlässigkeitsmatrix auf die rechte Seite übertragen.

\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}

Analyse des gesättigten und ungesättigten Permeationsflusses

Gleichzeitige Gleichungen bei der Analyse des gesättigten und ungesättigten Permeationsflusses weisen die folgenden Eigenschaften auf.

Es müssen drei Arten von Randbedingungen berücksichtigt werden: bekannte Durchflussrate, bekannte Gesamtförderhöhe und Grenze, an der eine infiltrierte Oberfläche auftritt.
Mit Ausnahme der Grenze, an der die infiltrierte Oberfläche erscheint, ist entweder die Durchflussrate oder die Gesamtförderhöhe bekannt.
An der Grenze, an der die infiltrierte Oberfläche erscheint, sind sowohl die Durchflussrate als auch die Gesamtförderhöhe unbekannt.
Da jedoch die Grenze, an der die infiltrierte Oberfläche erscheint, im Grunde die Grenze in Kontakt mit der Atmosphäre ist, ist die Durchflussrate nur ein Abfluss (0 oder weniger im Programm), sofern nicht der Zufluss von Regen berücksichtigt wird. Es gibt eine Einschränkung, dass die Druckhöhe ebenfalls 0 oder weniger beträgt.
Die Komponente der Wasserdurchlässigkeitsmatrix im gesättigten Bereich ist eine Konstante, die Komponente der Wasserpermeationsmatrix im ungesättigten Bereich ist jedoch eine Funktion des Saugdrucks (Unterdruckhöhe).

Daher ist es notwendig, das Unbekannte durch Konvergenzberechnung zu bestimmen. Das Verfahren der Konvergenzberechnung war ein sequentielles Substitutionsverfahren, bei dem die Anfangswerte aller Köpfe allen Knoten gegeben wurden und die erhaltene Lösung als Eingabewert für die nächste Konvergenzberechnung verwendet wurde. Da dieses Verfahren die Köpfe aller Knoten annimmt, kann es aus den Köpfen unter der Annahme der Wasserpermeationsmatrix aller Elemente eingestellt werden, und der Berechnungsfluss kann vereinfacht werden. Es scheint effektiv zu sein, den minimalen Gesamtkopf des Sättigungsbereichs als Anfangswert des Gesamtkopfs anzugeben, mit Ausnahme der Grenze, an der der Gesamtkopf bekannt ist, unter Berücksichtigung der Genauigkeit und Konvergenz der Lösung.

Ungesättigte Wasserdurchlässigkeit (Modell van Genuchten)

Für das charakteristische Modell der ungesättigten Wasserpermeation des Bodens gibt es auch eine Methode zur Eingabe der Sättigungs-Saugdruck-Beziehung und des Verhältnisses von Sättigungs-ungesättigtem Wasserpermeationskoeffizienten in eine Tabelle. Das Modell van Genuchten wurde übernommen.

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

Programm

Wenn Sie das Programm darauf setzen, wird es lang sein, also setzen Sie einen Link zu Ihrer eigenen HP.

Eingabedatenformat

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	Anzahl der Knoten, Anzahl der Elemente, Anzahl der Materialeigenschaften
koh, koq, kou	Anzahl der bezeichneten Wasserkopfknoten, Anzahl der angegebenen Durchflussknoten, Anzahl der Infiltrationsgrenzknoten
idan	Analyseabschnitt (0: vertikaler Abschnitt, 1: horizontaler Abschnitt)
Ak0	Durchlässigkeitskoeffizient für gesättigtes Wasser des Elements
alpa, em	Ungesättigte Wasserdurchlässigkeit von Elementen
Knoten-1, Knoten-2, Knoten-3, Knoten-4, isec	Knotennummer des Elements, Materialeigenschaftsnummer des Elements
x-, z-, hvec0	x- und z-Koordinaten des Knotens, anfänglicher Kopfwert für die Konvergenzberechnung
nokh, Hinp	Von Waterhead festgelegte Knotennummer und Waterhead-Wert
nokq, Qinp	Geben Sie die angegebene Knotennummer und den angegebenen Flusswert an
noku	Nummer des Invasionsgrenzknotens

Ausgabedatenformat

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

Knoten, hvec, pvec, qvec	Knotennummer, Gesamthöhe, Druckhöhe, Durchflussrate
elem, vx, vz, vm	Geschwindigkeit der X / z-Richtung und synthetische Geschwindigkeit der Elemente
Kr	Durchlässigkeit für ungesättigtes Wasser (1: gesättigt, 1>: ungesättigt)
Gesamtzufluss, Gesamtabfluss	Gesamtzufluss und -abfluss zum Modellbereich
iii, icount	Anzahl der Konvergenzberechnungen, Anzahl der Freiheitsgrade, die die Bedingung für den Abschluss der Berechnung erfüllen
kop	Die Anzahl der Knoten, bei denen der Druckkopf 0 ist, unter den infiltrierten Oberflächengrenzknoten
n, Zeit	Gesamtfreiheit, Berechnungszeit

Ausgabebeispiel

Ein Diagramm mit Konturen mit Druckköpfen von 0, 5, 10 und 15 m, die mit einem einfachen Programm zur Konturerstellung erstellt wurden. Die Wassertiefe beträgt 18 m stromaufwärts und 4 m stromabwärts. Das Konturdiagramm selbst befindet sich auf Amateurebene, aber die Berechnung scheint gut zu funktionieren. Die rot angezeigten Druckkopf- und Dammspezifikationen sind zusätzlich in der Mac-Vorschau angegeben.

Es ist besser, die Zeichnung dieses Ortes der Konturfunktion von GMT (Generic Mapping Tools) zu überlassen, und es wird für den Bericht besser aussehen.

Übrigens ist bekannt, dass bei einem tatsächlichen Damm die Druckhöhe nahe der Grenze zwischen Kern und Filter abfällt und der Druckabfall innerhalb des Kerns nicht wie in den Analyseergebnissen gezeigt auftritt. Bitte haben Sie Verständnis dafür, dass es sich um einen Analysefall handelt.

das ist alles