Zweidimensionale instationäre Wärmeleitungsanalyse mit Python

Überblick

Ich habe mit Python ein zweidimensionales Programm zur Analyse der nicht stationären Wärmeleitung erstellt. Das Element ist ein isoparametrisches Element mit 4 Knoten und 4 Gauß-Integrationspunkten, was dem in der Permeationsflussanalyse entspricht.

Begrenzte Elementformel für die instationäre Wärmeleitungsanalyse

Hier wird die Finite-Elemente-Formel der folgenden Elemente verwendet.

[\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}] $: Wärmeleitungsmatrix. Der zweite Term ist auf die Wärmeübertragungsgrenze zurückzuführen
$ [\ boldsymbol {c}] $: Wärmekapazitätsmatrix
$ \ {\ boldsymbol {\ phi} } $: Knotentemperaturvektor
$ \ {\ boldsymbol {f} } $: Wärmeflussvektor. Der erste Term ist auf die Wärmeerzeugung zurückzuführen, und der zweite Term ist auf die Wärmeübertragung zurückzuführen.
$ [\ boldsymbol {N}] $: Interpolierende Funktionsmatrix
$ \ Kappa $, $ \ rho $, $ c $: Wärmeleitfähigkeit, Dichte, spezifische Wärme
$ \ Alpha_c $, $ T_c $: Wärmeübertragungsrate, Außentemperatur der Wärmeübertragungsgrenze
$ \ dot {Q} $: Wärmeerzeugungsrate
$ \ Int_A $, $ \ int_S $: Elementbereich und Integration an den Elementseiten

Berücksichtigung der Wärmeerzeugung

Unter Berücksichtigung der Wärmeleitungsanalyse von Beton kann der für Beton charakteristische Temperaturanstieg der Wärmeisolierung einbezogen werden.

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

Hier ist $ T_k $ der endgültige Temperaturanstieg und $ \ alpha $ der Parameter, der die Wärmeerzeugungsrate angibt. Unter Verwendung der obigen Gleichung können der Heizwert $ Q $ und der Heizwert $ \ dot {Q} $ wie folgt ausgedrückt werden. $ \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} $

Zeitverteilungsmethode

Die Crank-Nicolson-Differenzgleichung wird als Zeitdispersionsmethode zum Lösen instationärer Finite-Elemente-Gleichungen verwendet. Tatsächlich wird die Differenzformel verwendet und die Knotentemperatur bei $ t + \ Delta t $ wird nacheinander unter Verwendung des bekannten Vektors zum Zeitpunkt $ t $ erhalten. Die Symbole werden großgeschrieben, um den Matrixvektor für den gesamten Analysebereich anzuzeigen.

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

Da die Materialeigenschaften der Elemente so eingestellt sind, dass sie sich nicht mit der Zeit ändern, wird bei der Zeitschrittberechnung die inverse Matrix nur einmal mit `numpy.linalg.inv (A)` berechnet, und diese wird jeweils gespeichert Der Temperaturverlauf der Zeit wird berechnet.

Randbedingung

Randbedingungen können Isolierung, Temperaturfixierung und Wärmeübertragung berücksichtigen.
Wenn nichts als Randbedingung eingegeben wird, handelt es sich um eine adiabatische Grenze.
Geben Sie für feste Temperaturgrenzen die Knotentemperatur an.
Geben Sie für die Wärmeübertragungsgrenze die Seite des Elements und die Außentemperatur an.

Kann es verwendet werden?

Mit 400 Elementen, 441 Freiheitsgraden und 1000 Berechnungszeitschritten dauert es ungefähr 38 Sekunden, ob es verwendet werden kann. Die Bewertung hängt vom individuellen Geschmack ab, aber in meinem Fall ist es keine Zeit, die ich nicht abwarten kann. In Fortran ist die Berechnung des gleichen Modells übrigens in 1,6 Sekunden abgeschlossen.

Programm

Ich habe gleichzeitig mit dem FEM-Programm ein Zeichenprogramm erstellt. Die Sprache ist Python, aber das Temperaturverlaufsdiagramm wird mit matplotlib gezeichnet, und die dreidimensionale Temperaturverteilungskarte wird mit GMT (Generic Mapping Tools) gezeichnet. Das 3D-Diagramm von Matplotlib ist nicht sehr cool.

Beispiele für Programm- und Eingabedaten werden mit Links zu den in Gist eingefügten Daten angezeigt.

Ausführungsskript

Ausführungsskript für die FEM-Berechnung

`a_py.txt`


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

Ausführungsskript zeichnen

Skript zum Erstellen eines Temperaturverlaufsdiagramms von Python und einer 3D-Temperaturverteilungskarte von 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

PDF-Erstellung mit TeX und PNG-Bilderstellung mit ImageMagick

Ein Skript, das ein PNG-Bild mit von ImageMagick entfernten (angepassten) Rändern erstellt, nachdem mit TeX ein PDF erstellt wurde

`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

TeX-Dokument

TeX-Dokument zum Anordnen von 8 PNG-Bildern nebeneinander auf A3-Papier

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

Eingabedatenformat

Zwei Dateien werden als Eingabedaten verwendet.

Eines sind Modelldaten. Geben Sie die Anzahl der Knoten, Grundgrößen wie Anzahl der Elemente, Materialeigenschaften, Element-Knoten-Beziehung, Knotenkoordinaten, Randbedingungen usw. ein.

Die andere ist eine Zeitverlaufsdatendatei der Außentemperatur der Temperaturbezeichnungsgrenze / Wärmeübertragungsgrenze. Da die Informationsmenge groß ist, habe ich versucht, sie getrennt von den Modelldaten einzugeben.

Modelldaten

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	Anzahl der Knoten, Anzahl der Elemente, Anzahl der Materialeigenschaften
kot, koc	Anzahl der temperaturspezifizierten Grenzknoten, Anzahl der Seiten der Wärmeübertragungsgrenzelemente
Delta	Berechnungszeitintervall
Ak, Ac, Arho	Wärmeleitfähigkeit, spezifische Wärme, Dichte
Tk, Al	Maximale adiabatische Anstiegstemperatur, Parameter der Wärmeerzeugungsrate
x, y, T0	Knoten x-Koordinate, Knoten y-Koordinate, Knotenanfangstemperatur
nokt	Grenzwertknotennummer der Temperaturspezifikation
nekc_0, nekc_1, alphac	Elementnummer mit Wärmeübertragungsgrenze, Knotennummer der Seite der Wärmeübertragungsgrenze, Wärmeübertragungsrate der Wärmeübertragungsgrenze
n1out	Anzahl der Knoten, die den gesamten Zeitverlauf ausgeben sollen
n1node	Knotennummer zur Ausgabe des gesamten Zeitverlaufs
n2out	Zeit für die Ausgabe aller Knotentemperaturen Anzahl der Schritte
n2step	Zeit für die Ausgabe aller Knotentemperaturen Schrittnummer

An der Wärmeübertragungsgrenze müssen Sie die Seiten des Elements angeben. Geben Sie daher die Elementnummer mit der Wärmeübertragungsgrenzkante und die Knotennummer mit der Kante an. In FEM wird die Knotennummer, die das Element angibt, gegen den Uhrzeigersinn eingegeben, sodass Sie die Seite identifizieren können, indem Sie die erste Knotennummer eingeben, aus der die Seite besteht.

Temperaturbezeichnungsgrenze / Wärmeübertragungsgrenze Externe Temperaturdaten

   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: Anzahl der temperaturspezifizierten Grenzknoten
koc: Anzahl der Seiten des Wärmeübertragungsgrenzelements
TE [1 ~ kot]: Temperaturverlauf an der angegebenen Temperaturgrenze
TC [1 ~ koc]: Externer Temperaturverlauf an der Wärmeübertragungsgrenze

Ausgabedatenformat

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	Anzahl der Berechnungsschritte (einschließlich der Anfangstemperatur)
Tfix	Gibt an, ob die Knotentemperatur angegeben ist oder nicht (0: Keine Temperatur angegeben, 1: Angabe der Temperatur)
iii, ttime, Node_x	Knotennummer zur Ausgabe von Zeitschritt, Zeit, Temperaturverlauf Danach wird jeweils der Temperaturverlauf