Two-dimensional unsteady heat conduction analysis with Python

Overview

I made a two-dimensional unsteady heat conduction analysis program in Python. The element is an isoparametric element with 4 nodes and 4 Gauss integral points, which is the same as the osmotic flow analysis.

Finite element formula for unsteady heat conduction analysis

Here, the finite element formula of the following elements is used.

[\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}] $: Heat conduction matrix. The second term is due to the heat transfer boundary
$ [\ boldsymbol {c}] $: Heat capacity matrix
$ \ {\ boldsymbol {\ phi} } $: Nodal temperature vector
$ \ {\ boldsymbol {f} } $: Heat flux vector. The first term is due to heat generation, and the second term is due to heat transfer.
$ [\ boldsymbol {N}] $: Interpolated function matrix
$ \ Kappa $, $ \ rho $, $ c $: Thermal conductivity, density, specific heat
$ \ Alpha_c $, $ T_c $: Heat transfer coefficient, external temperature at the heat transfer boundary
$ \ dot {Q} $: Heat generation rate
$ \ Int_A $, $ \ int_S $: Integral on element area and element edge

Consideration of fever

Considering the heat conduction analysis of concrete, the heat insulation temperature rise characteristic of concrete can be included.

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

Here, $ T_k $ is a parameter indicating the final temperature rise, and $ \ alpha $ is a parameter indicating the heat generation rate. Using the above equation, the calorific value $ Q $ and the calorific value $ \ dot {Q} $ can be expressed as follows. $ \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} $

Time discretization method

The Crank-Nicolson finite difference method is used as a time discretization method for solving unsteady finite element equations. Actually, the difference equation is used, and the node temperature at $ t + \ Delta t $ is sequentially obtained using the known vector at time $ t $. The symbols are capitalized to indicate the matrix vector for the entire analysis range.

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

Since the material properties of the elements are set so that they do not change with time, in the time step calculation, the inverse matrix is calculated only once with `numpy.linalg.inv (A)`, and this is stored and each is stored. The temperature history of the time is calculated.

boundary condition

Boundary conditions can take into account insulation, temperature fixation, and heat transfer.
If nothing is entered as a boundary condition, it will be an adiabatic boundary.
For fixed temperature boundaries, specify the nodal temperature.
For heat transfer boundaries, specify the sides of the element and the external temperature.

Can it be used?

In terms of whether it can be used, it takes about 38 seconds with 400 elements, 441 degrees of freedom, and 1000 calculation time steps. Evaluation depends on individual taste, but in my case it is not a time I can't wait. By the way, in Fortran, the calculation of the same model is completed in 1.6 seconds.

program

At the same time as the FEM program, I made a drawing program. The language is Python, but the temperature history graph is drawn with matplotlib, and the 3D temperature distribution map is drawn with GMT (Generic Mapping Tools). The 3D graph of matplotlib is not very cool.

Program and input data examples are shown with links to those pasted in Gist.

Script for execution

FEM calculation execution script

`a_py.txt`


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

Drawing execution script

Script for creating temperature history graph by Python and 3D temperature distribution map by 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 creation with TeX and png image creation with ImageMagick

A script that creates a png image with margins removed (adjusted) by ImageMagick after creating a pdf with 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

TeX document

TeX document for arranging 8 png images side by side on A3 paper

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

Input data format

Two files are used as input data.

One is model data. Enter the number of nodes, basic quantities such as the number of elements, material properties, element-node relationship, node coordinates, boundary conditions, etc.

The other is a time history data file of the external temperature of the temperature designation boundary / heat transfer boundary. Since the amount of information is large, I tried to input it separately from the model data.

Model data

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	Number of nodes, number of elements, number of material characteristics
kot, koc	Number of temperature-specified boundary nodes, number of heat transfer boundary element sides
delta	Calculation time interval
Ak, Ac, Arho	Thermal conductivity, specific heat, density
Tk, Al	Maximum adiabatic temperature rise, heat generation rate parameters
x, y, T0	Node x coordinate, node y coordinate, node initial temperature
nokt	Temperature-specified boundary node number
nekc_0, nekc_1, alphac	Element number with heat transfer boundary, node number of heat transfer boundary element side, heat transfer coefficient of heat transfer boundary
n1out	Number of nodes to output all time history
n1node	Node number that outputs all time history
n2out	Time to output all node temperatures Number of steps
n2step	Time to output the temperature of all nodes Step number

At the heat transfer boundary, you need to specify the sides of the element. Therefore, specify the element number that has the heat transfer boundary edge and the node number that specifies the edge. In FEM, the node number that specifies the element is entered counterclockwise, so you can identify the edge by entering the first node number that makes up the edge.

Temperature designation boundary / heat transfer boundary External temperature data

   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: Number of temperature specified boundary nodes
koc: Number of heat transfer boundary element edges
TE [1 ~ kot]: Temperature history at the specified temperature boundary
TC [1 ~ koc]: External temperature history at the heat transfer boundary

Output data format

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	Number of calculation steps (including initial temperature)
Tfix	Whether or not the node temperature is specified (0: No temperature specified, 1: Temperature specified)
iii, ttime, Node_x	Node number to output time step, time, temperature history After that, temperature history at each time is output