Analyse de conduction thermique bidimensionnelle non stationnaire avec Python

Aperçu

J'ai réalisé un programme d'analyse de conduction thermique non stationnaire bidimensionnel avec Python. L'élément est un élément isoparamétrique avec 4 nœuds et 4 points d'intégration de Gauss, qui est le même que dans l'analyse de flux de perméation.

Formule d'élément limitée pour l'analyse par conduction thermique non stationnaire

Ici, la formule des éléments finis des éléments suivants est utilisée.

[\boldsymbol{k}]\\{\boldsymbol{\phi}\\}+[\boldsymbol{c}]\left\\{\cfrac{\partial \boldsymbol{\phi}}{\partial t}\right\\}=\\{\boldsymbol{f}\\}

\begin{align} &[\boldsymbol{k}]=\int_A \kappa\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)dA +\int_s \alpha_c[\boldsymbol{N}]^T [\boldsymbol{N}]ds \\\ &[\boldsymbol{c}]=\int_A \rho c [\boldsymbol{N}]^T [\boldsymbol{N}]dA \\\ &\{\boldsymbol{f}\}=\int_A \dot{Q}[\boldsymbol{N}]^T dA+\int_s \alpha_c T_c [\boldsymbol{N}]^T ds \end{align}

$ [\ boldsymbol {k}] $: matrice de conduction thermique. Le deuxième terme est dû à la limite de transfert de chaleur
$ [\ boldsymbol {c}] $: matrice de capacité thermique
$ \ {\ boldsymbol {\ phi} } $: vecteur de température du nœud
$ \ {\ boldsymbol {f} } $: vecteur de flux thermique. Le premier terme est dû à la génération de chaleur, et le second terme est dû au transfert de chaleur.
$ [\ boldsymbol {N}] $: matrice de fonction d'interpolation
$ \ Kappa $, $ \ rho $, $ c $: conductivité thermique, densité, chaleur spécifique
$ \ Alpha_c $, $ T_c $: Taux de transfert de chaleur, température externe de la limite de transfert de chaleur
$ \ dot {Q} $: taux de génération de chaleur
$ \ Int_A $, $ \ int_S $: Zone d'élément et intégration aux côtés de l'élément

Prise en compte de la génération de chaleur

Compte tenu de l'analyse de conduction thermique du béton, la caractéristique d'élévation de température d'isolation thermique du béton peut être incluse.

T(t)=T_k (1-e^{-\alpha t})

Ici, $ T_k $ est le montant final de l'augmentation de la température et $ \ alpha $ est le paramètre indiquant le taux de génération de chaleur. En utilisant l'équation ci-dessus, la valeur calorifique $ Q $ et la valeur calorifique $ \ dot {Q} $ peuvent être exprimées comme suit. $ \begin{equation} Q=\rho \cdot c \cdot T(t) \qquad \rightarrow \qquad \dot{Q}=\rho \cdot c \cdot T_k \cdot \alpha \cdot e^{-\alpha t} \end{equation} $

Méthode de dispersion temporelle

L'équation de différence Crank-Nicolson est utilisée comme méthode de dispersion temporelle pour résoudre des équations d'éléments finis non stationnaires. En fait, la formule de différence est utilisée, et la température du nœud à $ t + \ Delta t $ est séquentiellement obtenue en utilisant le vecteur connu au temps $ t $. Les symboles sont en majuscules pour indiquer le vecteur de matrice pour toute la plage d'analyse.

\begin{equation} \left(\frac{1}{2}[\boldsymbol{K}]+\frac{1}{\Delta t}[\boldsymbol{C}]\right)\{\boldsymbol{\Phi}(t+\Delta t)\} =\left(-\frac{1}{2}[\boldsymbol{K}]+\frac{1}{\Delta t}[\boldsymbol{C}]\right)\{\boldsymbol{\Phi}(t)\}+\frac{\{\boldsymbol{F}(t+\Delta t)\}+\{\boldsymbol{F}(t)\}}{2} \end{equation}

Étant donné que les propriétés matérielles des éléments sont définies de manière à ne pas changer avec le temps, dans le calcul du pas de temps, la matrice inverse n'est obtenue qu'une seule fois avec `` numpy.linalg.inv (A) '', et celle-ci est stockée et chacune est stockée. L'historique de température de l'heure est calculé.

condition limite

Les conditions aux limites peuvent prendre en compte l'isolation, la fixation de la température et le transfert de chaleur.
Si rien n'est entré comme condition aux limites, ce sera une limite adiabatique.
Pour les limites de température fixes, spécifiez la température du nœud.
Pour la limite de transfert de chaleur, spécifiez le côté de l'élément et la température extérieure.

Peut-il être utilisé?

Pour savoir s'il peut être utilisé, cela prend environ 38 secondes avec 400 éléments, 441 degrés de liberté et 1000 pas de temps de calcul. L'évaluation dépend du goût individuel, mais dans mon cas, ce n'est pas un moment que je ne peux pas attendre. D'ailleurs, dans Fortran, le calcul du même modèle est terminé en 1,6 seconde.

programme

J'ai fait un programme de dessin en même temps que le programme FEM. Le langage est Python, mais le graphique de l'historique des températures est dessiné avec matplotlib, et la carte de distribution de température en trois dimensions est dessinée avec GMT (Generic Mapping Tools). Le graphique 3D de matplotlib n'est pas très cool.

Des exemples de programmes et de données d'entrée sont affichés avec des liens vers ceux collés dans Gist.

Script d'exécution

Script d'exécution de calcul FEM

`a_py.txt`


python3 py_fem_heat.py inp_div20_model.txt inp_div20_thist.txt out_div20.txt

Script d'exécution de dessin

Script pour créer un graphique d'historique de température par Python et une carte de distribution de température 3D par GMT

`a_fig.txt`


python3 py_heat_2D.py out_11_div20_10.txt
python3 py_heat_3D.py out_11_div20_10.txt
bash _a_gmt_3D.txt
mogrify -trim -density 300 -bordercolor 'transparent' -border 10x10 -format png *.eps
rm *.eps

Création de PDF avec TeX et création d'images PNG avec ImageMagick

Un script qui crée une image png avec des marges supprimées (ajustées) par ImageMagick après avoir créé un pdf avec TeX

`a_tex.txt`


platex tex_fig.tex
dvipdfmx -p a3 tex_fig.dvi

convert -trim -density 400 tex_fig.pdf -bordercolor 'transparent' -border 20x20 -quality 100 tex_fig.png
tex_fig.tex

Document TeX

Document TeX pour organiser 8 images png côte à côte sur du papier A3

`tex_fig.tex`


\documentclass[english]{jsarticle}
\usepackage[a3paper,top=25mm,bottom=25mm,left=25mm,right=25mm]{geometry}
\usepackage[dvipdfmx]{graphicx}
\pagestyle{empty}

\begin{document}

\begin{center}
\begin{tabular}{cc}
\begin{minipage}{12.0cm}\vspace{0.2zh}\includegraphics[width=12.0cm,bb={0 0 1900 1122}]{_fig_gmt_xyz_0.png}\end{minipage}&
\begin{minipage}{12.0cm}\vspace{0.2zh}\includegraphics[width=12.0cm,bb={0 0 1900 1122}]{_fig_gmt_xyz_1.png}\end{minipage}\\
 & \\
 & \\
\begin{minipage}{12.0cm}\vspace{0.2zh}\includegraphics[width=12.0cm,bb={0 0 1900 1122}]{_fig_gmt_xyz_2.png}\end{minipage}&
\begin{minipage}{12.0cm}\vspace{0.2zh}\includegraphics[width=12.0cm,bb={0 0 1900 1122}]{_fig_gmt_xyz_3.png}\end{minipage}\\
 & \\
 & \\
\begin{minipage}{12.0cm}\vspace{0.2zh}\includegraphics[width=12.0cm,bb={0 0 1900 1122}]{_fig_gmt_xyz_4.png}\end{minipage}&
\begin{minipage}{12.0cm}\vspace{0.2zh}\includegraphics[width=12.0cm,bb={0 0 1900 1122}]{_fig_gmt_xyz_5.png}\end{minipage}\\
 & \\
 & \\
\begin{minipage}{12.0cm}\vspace{0.2zh}\includegraphics[width=12.0cm,bb={0 0 1900 1122}]{_fig_gmt_xyz_6.png}\end{minipage}&
\begin{minipage}{12.0cm}\vspace{0.2zh}\includegraphics[width=12.0cm,bb={0 0 1900 1122}]{_fig_gmt_xyz_7.png}\end{minipage}\\
\end{tabular}
\end{center}

\end{document}

Format des données d'entrée

Deux fichiers sont utilisés comme données d'entrée.

Le premier concerne les données de modèle. Entrez le nombre de nœuds, les quantités de base telles que le nombre d'éléments, les propriétés du matériau, la relation élément-nœud, les coordonnées des nœuds, les conditions aux limites, etc.

L'autre est un fichier de données d'historique temporel de la température externe de la limite de désignation de température / limite de transfert de chaleur. Étant donné que la quantité d'informations est importante, j'ai essayé de les saisir séparément des données du modèle.

Données du modèle

npoin  nele  nsec  kot  koc  delta          # Basic values for analysis
Ak  Ac  Arho  Tk  Al                        # Material properties
     ....(1 to nsec)....                    #
node-1  node-2  node-3  node-4  isec        # Element connectivity, Material set number
     ....(1 to nele)....                    # (counterclockwise order of node numbers)
x  y  T0                                    # Coordinates of node, initial temperature of node
     ....(1 to npoin)....                   # (x: right-direction, y: upward direction)
nokt                                        # Node number of temperature given boundary
     ....(1 to kot)....                     # (Omit data input if kot=0)
nekc_0  nekc_1  alphac                      # Element number, node number defined side, heat transfer rate
     ....(1 to koc)....                     # (Omit data input if koc=0)
n1out                                       # Number of nodes for output of temperature time history
n1node_1   ....(1 line input)....           # Node number for above
n2out                                       # frequency of output of temperatures of all nodes
n2step_1  ....(1 line input)....            # Time step number for above
                                            # (Omit data input if ntout=0)

npoin, nele, nsec	Nombre de nœuds, nombre d'éléments, nombre de propriétés du matériau
kot, koc	Nombre de nœuds limites spécifiés par la température, nombre de côtés d'élément limite de transfert thermique
delta	Intervalle de temps de calcul
Ak, Ac, Arho	Conductivité thermique, chaleur spécifique, densité
Tk, Al	Température de montée adiabatique maximale, paramètres de taux de génération de chaleur
x, y, T0	Coordonnée noeud x, coordonnée noeud y, température initiale du noeud
nokt	Numéro de nœud de limite de spécification de température
nekc_0, nekc_1, alphac	Numéro d'élément avec limite de transfert de chaleur, numéro de nœud côté élément de limite de transfert de chaleur, taux de transfert de chaleur de la limite de transfert de chaleur
n1out	Nombre de nœuds à afficher tout l'historique de temps
n1node	Numéro de nœud pour afficher tout l'historique de temps
n2out	Temps de sortie de toutes les températures de nœud Nombre d'étapes
n2step	Temps de sortie de toutes les températures de noeud Numéro d'étape

À la limite du transfert de chaleur, vous devez spécifier les côtés de l'élément. Par conséquent, spécifiez le numéro d'élément qui a le bord de la limite du transfert de chaleur et le numéro de nœud qui spécifie le bord. Dans FEM, le numéro de nœud qui spécifie l'élément est entré dans le sens inverse des aiguilles d'une montre, vous pouvez donc identifier le côté en entrant le premier numéro de nœud qui constitue le côté.

Limite de désignation de température / limite de transfert de chaleur Données de température externe

   1   TE[1] ... TE[kot] TC[1] ... TC[koc]   # data of 1st time step
   2   TE[1] ... TE[kot] TC[1] ... TC[koc]   # data of 2nd time step
   3   TE[1] ... TE[kot] TC[1] ... TC[koc]   #
    ....(repeating until end time).....      # data is read step by step
itime  TE[1] ... TE[kot] TC[1] ... TC[koc]   #

kot: nombre de nœuds limites spécifiés par la température
koc: nombre de côtés de l'élément limite de transfert de chaleur
TE [1 ~ kot]: Historique des températures à la limite de température spécifiée
TC [1 ~ koc]: historique de la température externe à la limite du transfert de chaleur

Format des données de sortie

npoin  nele  nsec   kot   koc           delta  niii n1out n2out
   25    16     1     0    16   1.0000000e+00   101     5     0
  sec              Ak              Ac            Arho              Tk              Al
    1   2.5000000e+00   2.8000000e-01   2.3500000e+03   4.0000000e+01   2.0000000e-01
 node               x               y          tempe0  Tfix
    1  -5.0000000e-01  -5.0000000e-01   2.0000000e+01     0
    2  -5.0000000e-01  -2.5000000e-01   2.0000000e+01     0
    3  -5.0000000e-01   0.0000000e+00   2.0000000e+01     0
..... (omit) .....
0
   23   5.0000000e-01   0.0000000e+00   2.0000000e+01     0
   24   5.0000000e-01   2.5000000e-01   2.0000000e+01     0
   25   5.0000000e-01   5.0000000e-01   2.0000000e+01     0
 nek0  nek1 alphac         
    1     1   1.0000000e+01
    5     6   1.0000000e+01
    9    11   1.0000000e+01
   13    16   1.0000000e+01
   13    21   1.0000000e+01
   14    22   1.0000000e+01
   15    23   1.0000000e+01
   16    24   1.0000000e+01
   16    25   1.0000000e+01
   12    20   1.0000000e+01
    8    15   1.0000000e+01
    4    10   1.0000000e+01
    4     5   1.0000000e+01
    3     4   1.0000000e+01
    2     3   1.0000000e+01
    1     2   1.0000000e+01
 elem     i     j     k     l   sec
    1     1     6     7     2     1
    2     2     7     8     3     1
    3     3     8     9     4     1
    4     4     9    10     5     1
    5     6    11    12     7     1
    6     7    12    13     8     1
    7     8    13    14     9     1
    8     9    14    15    10     1
    9    11    16    17    12     1
   10    12    17    18    13     1
   11    13    18    19    14     1
   12    14    19    20    15     1
   13    16    21    22    17     1
   14    17    22    23    18     1
   15    18    23    24    19     1
   16    19    24    25    20     1
  iii           ttime         Node_11         Node_12         Node_13         Node_14         Node_15
    0   0.0000000e+00   2.0000000e+01   2.0000000e+01   2.0000000e+01   2.0000000e+01   2.0000000e+01
    1   1.0000000e+00   2.5712484e+01   2.7416428e+01   2.7023876e+01   2.7416428e+01   2.5712484e+01
    2   2.0000000e+00   2.8833670e+01   3.3625990e+01   3.2820487e+01   3.3625990e+01   2.8833670e+01
..... (omit) .....
   98   9.8000000e+01   1.1182079e+01   1.2141098e+01   1.2494074e+01   1.2141098e+01   1.1182079e+01
   99   9.9000000e+01   1.1140143e+01   1.2065139e+01   1.2405593e+01   1.2065139e+01   1.1140143e+01
  100   1.0000000e+02   1.1099694e+01   1.1991874e+01   1.2320250e+01   1.1991874e+01   1.1099694e+01
n=25  time=0.147 sec

niii	Nombre d'étapes de calcul (y compris la température initiale)
Tfix	Indique si la température du nœud est spécifiée ou non (0: aucune température spécifiée, 1: température spécifiée)
iii, ttime, Node_x	Numéro de nœud pour afficher le pas de temps, l'heure, l'historique de la température Après cela, l'historique de la température à chaque fois est sorti