Analyse bidimensionnelle du flux de perméation saturée-insaturée avec Python

(Corrigé le 19 novembre 2017)

Correction du programme en raison de certaines erreurs de programmation

Aperçu

Impression approximative

Depuis que j'ai réalisé un programme d'analyse de contraintes FEM bidimensionnel, j'ai également réalisé un programme d'analyse d'écoulement de perméation saturée-insaturée bidimensionnelle. L'élément est un élément isoparamétrique ayant 1 élément, 4 nœuds et 4 points d'intégration de Gauss, comme dans le cas de l'analyse de contraintes 2D.

En fait, il peut être utilisé, mais le nombre de nœuds est de 1317, le nombre d'éléments est de 1227, le nombre de calculs de convergence est de 10 et le temps de calcul est de 3,7 secondes. Si le problème est de cette échelle, aucune contrainte particulière n'est ressentie. Les opérations matricielles et les routines de Numpy pour résoudre des équations linéaires simultanées sont probablement excellentes.

Les éléments à garder à l'esprit lors de la programmation sont les suivants.

Tous les tableaux sont des tableaux Numpy. Aucune liste n'est utilisée. Évitez d'utiliser + autant que possible et utilisez les opérations matricielles Numpy. Cependant, veuillez noter que des résultats inattendus peuvent survenir en raison du produit de vecteurs unidimensionnels.
x = numpy.linalg.solve (A, b) `` est suffisant pour résoudre des équations linéaires simultanées. Étant donné que la matrice devient asymétrique en raison du traitement des conditions aux limites, la décomposition choleskey ne peut pas être utilisée dans les deux cas.
En faisant ce qui précède, les perspectives du programme lui-même seront très bonnes.
La section entrée / sortie est en panne, est-ce inévitable?

Analyse de flux de perméation saturée

L'analyse de flux de perméation stationnaire saturée est un problème qui résout les équations linéaires simultanées suivantes.

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

Ici, $ \ {q \} $ est le vecteur d'écoulement nodal, $ \ {h \} $ est le vecteur de tête totale, et $ [k] $ est la matrice de perméabilité à l'eau, qui présente les caractéristiques suivantes.

Matrice de perméation d'eau $ [k] $ est une constante,
Le vecteur de tête total $ \ {h \} $ correspondant à la composante connue du vecteur de flux nodal $ \ {q \} $ est inconnu
Le vecteur de tête total $ \ {h \} $ correspondant à la composante inconnue du vecteur de flux nodal $ \ {q \} $ est un nombre connu

Par conséquent, la solution peut être obtenue en résolvant les équations simultanées une fois après avoir traité la relation entre le nombre connu et le nombre inconnu.

(Traitement de la relation entre les nombres connus et inconnus pour résoudre des équations linéaires simultanées)

Débit hors $ q_i $, $ q_j $ et connu $ h_i $, $ h_j $, toutes les têtes hors $ h_i $, $ h_j $ et $ q_i $, $ q_j $ dans la matrice de perméabilité à l'eau Processus dans lesquels k_ {ii} $ et $ k_ {jj} $ sont mis à $ 1 $, et les autres éléments de la colonne $ i $ et de la colonne $ j $ sont mis à zéro. Il transfère également les effets des colonnes $ i $ et $ j $ de la matrice de perméabilité à l'eau sur le côté droit.

\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 flux de perméation saturés-insaturés

Les équations simultanées dans l'analyse d'écoulement de perméation saturée-insaturée ont les caractéristiques suivantes.

Il est nécessaire de considérer trois types de conditions aux limites: le débit connu, la hauteur totale connue et la limite où apparaît une surface infiltrée.
Soit le débit, soit la hauteur totale est connu sauf pour la limite où apparaît la surface infiltrée.
Sur la limite où apparaît la surface infiltrée, le débit et la hauteur totale sont inconnus.
Cependant, comme la limite où apparaît la surface infiltrée est essentiellement la limite en contact avec l'atmosphère, le débit n'est que l'écoulement sortant (0 ou moins dans le programme) à moins que l'afflux de pluie ne soit pris en compte. Il existe une restriction selon laquelle la hauteur de pression est également égale ou inférieure à 0.
Le composant de la matrice de perméation de l'eau dans la région saturée est constant, mais le composant de la matrice de perméation de l'eau dans la région non saturée est fonction de la pression d'aspiration (hauteur de pression négative).

Par conséquent, il est nécessaire de déterminer l'inconnu par calcul de convergence. La méthode de calcul de convergence était une méthode de substitution séquentielle dans laquelle les valeurs initiales de toutes les têtes étaient données à tous les nœuds et la solution obtenue était utilisée comme valeur d'entrée pour le prochain calcul de convergence. Puisque cette méthode suppose les têtes de tous les nœuds, elle peut être définie à partir des têtes en supposant la matrice de perméation d'eau de tous les éléments, et le flux de calcul peut être simplifié. Il semble qu'il soit efficace de donner la hauteur totale minimale de la région de saturation comme valeur initiale de la hauteur totale, sauf pour la limite où la hauteur totale est connue, compte tenu de la précision et de la convergence de la solution.

Perméabilité à l'eau non saturée (modèle de van Genuchten)

Pour le modèle caractéristique de perméation d'eau non saturée du sol, il existe également une méthode de saisie de la relation de pression de saturation-aspiration et du rapport de coefficient de perméation d'eau saturation-insaturée dans un tableau, mais pour l'analyse, il est préférable de gérer avec une fonction continue car le programme peut être simplifié. Le modèle van Genuchten a été adopté.

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

programme

Si vous mettez le programme dessus, ce sera long, alors mettez un lien vers votre propre HP.

Format des données d'entrée

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	Nombre de nœuds, nombre d'éléments, nombre de propriétés du matériau
koh, koq, kou	Nombre de nœuds de tête d'eau désignés, nombre de nœuds de débit spécifié, nombre de nœuds de limite d'infiltration
idan	Section d'analyse (0: coupe verticale, 1: coupe horizontale)
Ak0	Coefficient de perméabilité à l'eau saturée de l'élément
alpa, em	Perméabilité à l'eau non saturée des éléments
node-1, node-2, node-3, node-4, isec	Numéros de nœud qui composent l'élément, numéro de propriété matérielle de l'élément
x, z, hvec0	coordonnées x et z du nœud, valeur de tête initiale pour le calcul de convergence
nokh, Hinp	Numéro de nœud désigné Waterhead et valeur Waterhead
nokq, Qinp	Numéro de nœud et valeur de flux spécifiés par le flux
noku	Numéro de nœud de limite d'invasion

Format des données de sortie

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	numéro de nœud, hauteur totale, hauteur de pression, débit
elem, vx, vz, vm	Vitesse de direction X / z et vitesse synthétique des éléments
Kr	Perméabilité à l'eau insaturée de l'élément (1: saturé, 1>: insaturé)
Flux entrant total, flux sortant total	Flux entrant et sortant total vers la zone du modèle
iii, icount	Nombre de calculs de convergence, nombre de degrés de liberté satisfaisant la condition de fin de calcul
kop	Le nombre de nœuds où la hauteur de pression est 0 parmi les nœuds limites de surface infiltrés
n, temps	Liberté totale, temps de calcul

Exemple de sortie

La figure qui a tiré le contour de la tête de pression 0, 5, 10, 15m par le programme de création de contour simple. La profondeur de l'eau est de 18 m en amont et de 4 m en aval. Le diagramme de contour lui-même est au niveau amateur, mais le calcul semble bien fonctionner. Les spécifications de la tête de pression et du barrage affichées en rouge sont également écrites dans l'aperçu Mac.

Il est préférable de laisser le dessin ici à la fonction de contour de GMT (Generic Mapping Tools) et cela sera meilleur pour le rapport.

À propos, il est bien connu que dans un barrage réel, la hauteur de pression chute près de la limite entre le noyau et le filtre, et la chute de pression à l'intérieur du noyau ne se produit pas comme indiqué dans les résultats de l'analyse. Veuillez comprendre qu'il s'agit d'un cas d'analyse.

c'est tout