2차원 비온도 열전도 분석

개요


파이톤으로 2차원 비평온 열전도 분석 프로그램을 만들어 보았다.
원소는 삼투분석과 마찬가지로 4 노드 4 고스 적분점을 가진 동위소 파라미터 원소다.

비온도 열전도 분석의 유한원 공식


여기에는 다음과 같은 요소의 유한 요소식을 사용한다.
$$
[\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}: 열전도 매트릭스. 두 번째 항목은 열전도 경계
  • 달러: 열용량 매트릭스
  • ${\boldsymbol{\phi}달러: 노드 온도 벡터
  • ${\boldsymbol{f}달러: 열통량 벡터. 첫 번째는 발열 때문이고, 두 번째는 전열 때문이다
  • 달러: 내장 함수 매트릭스
  • $\kappa$\rho$c: 열전도율, 밀도, 비열
  • $\alpha_c$, $T_c$: 전열 계수, 전열 경계의 외부 온도
  • $\dot{Q}달러: 발열률
  • $\int_A$, $\int_S달러: 원소 면적과 원소 변의 포인트
  • 열적 고려


    의식적으로 콘크리트에 대해 열전도 분석을 하고 콘크리트의 단열 온도 상승 특성을 첨가한다.
    $$
    T(t)=T_k (1-e^{-\alpha t})
    $$
    여기,$T최종 온도 상승량 k$$\alpha는 발열 속도를 나타내는 매개 변수입니다.
    상기 공식을 사용하면 열량 $Q$와 발열률 $\dot{Q}달러를 다음과 같이 표현할 수 있다.
    $$
    \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}
    $$

    시간 이산화 방법


    Crank-Nicolson 차분식을 비정상적인 유한원식을 해결하는 시간 이산화 방법으로 채택하다.
    실제로 차분 공식을 사용하여 시간 $t의 이미 알고 있는 벡터를 사용하여 $t+\Deltat달러의 노드 온도를 순서대로 계산합니다.
    전체 해석 범위의 행렬 벡터를 표시하기 위해 기호를 대문자로 설정합니다.
    $$
    \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}
    $$
    원소의 재료 특성은 시간에 따라 변하지 않는 설정이기 때문에 시간 단계 계산에서numpy.linalg.inv(A)는 역행렬을 한 번만 계산하고 이를 보존하여 각 시간의 온도 역사를 계산한다.

    경계 조건

  • 경계 조건은 단열, 온도 고정과 전열을 고려할 수 있다
  • 경계 조건으로 아무것도 넣지 않으면 단열 경계가 된다
  • 온도가 경계를 고정한 상태에서 노드의 온도를 지정한다.
  • 전열 경계의 부속품 가장자리와 외부 온도를 지정한다.
  • 쓸 수 있어요?


    사용에 있어서는 400가지 요소, 441의 자유도, 계산 시간 순서는 1000회, 38초 정도입니다. 평가는 개인의 느낌에 따라 하지만 제 경우는 기다릴 수 없는 시간이 아닙니다. 참고로 포란 동일 모델의 계산은 1.6초에 완성되었습니다.

    절차.


    FEM 프로그램과 함께 작도 프로그램을 만들었다. 언어는 파이톤, 온도 이력도는 matplotlib, 3차원 온도 분포도는 GMT(Generic Mapping Tools)가 그렸다. matplotlib의 3차원 그림이 별로 예쁘지 않기 때문이다.
    프로그램 및 입력 데이터 사례는 Gist에 붙여진 데이터로 연결됩니다.
  • 비온도 열전도 분석 프로그램 (py fem heat.py)
  • 노드 온도 시간 그래프 제작 (py heat 2D.py)
  • GMT를 이용한 3차원 온도 분포도 작성(py heat 3D.py)
  • 모델 데이터 사례(inp div20 모델.txt)
  • 온도 지정/전달 경계의 온도 기록 데이터 인스턴스 (inp div20 thist.txt)
  • 실행할 스크립트


    FEM 컴퓨팅 실행 스크립트


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

    드로잉 실행 스크립트


    파이썬의 온도 히스토리, GMT의 3차원 온도 분포도 제작 스크립트
    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
    

    TeX를 통한 pdf 및 ImageMagick을 통한 pg 이미지 만들기


    TeX로 pdf를 제작한 후 ImageMagick을 통해 경계를 삭제하고 pg 이미지를 생성하는 스크립트
    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 파일


    A3 용지에 8개의 pg 이미지를 나란히 구성하는 TeX 파일
    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}
    

    입력 데이터 형식


    데이터를 입력하려면 두 개의 파일을 사용하십시오.
    첫째, 모델 데이터, 노드 수, 요소 수 등 기본 제량, 재료 특성, 요소 노드 관계, 노드 좌표, 경계 조건 등을 입력한다.
    다른 하나는 온도 지정 경계와 전열 경계 외부 온도의 시간력 데이터 파일입니다. 정보량이 증가하기 때문에 모델 데이터와 분리해서 입력합니다.

    모델 데이터

    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
    노드 수, 요소 수, 재료 특성 수
    kot, koc
    온도 지정 경계 노드 수, 전열 경계 컴포넌트 변수
    delta
    시간 간격 계산
    Ak, Ac, Arho
    열전도율
    Tk, Al
    최고 단열 상승 온도, 발열 속도 파라미터
    x, y, T0
    노드 x 좌표, 노드 y 좌표, 노드 초기 온도
    nokt
    온도 지정 경계 노드 번호
    nekc_0, nekc_1, alphac
    전열 경계가 있는 소자 번호, 전열 경계 소자 경계의 노드 번호, 전열 경계의 전열 계수
    n1out
    모든 시간이 경과한 노드 수를 출력합니다
    n1node
    모든 시간이 경과한 노드 번호를 출력합니다
    n2out
    모든 노드의 온도를 출력하는 시간 보폭
    n2step
    모든 노드의 온도를 출력하는 시간 단계 번호
    전열 경계에서는 요소의 모서리를 지정해야 합니다. 따라서 전열 경계가 있는 요소 번호와 지정된 모서리가 있는 노드 번호를 지정합니다. FEM에서는 지정된 요소의 노드 번호를 시계 반대 방향으로 입력하므로 모서리를 구성하는 첫 번째 노드 번호를 입력하면 모서리를 지정할 수 있습니다.

    온도 지정 경계/전열 경계 외부 온도 데이터

       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: 온도 지정 경계 노드 수
  • koc: 전열 경계 소자 변수
  • TE[1~kot]: 온도가 지정한 경계에 있는 온도 기록
  • TC[1~koc]: 전열 경계부의 외부 온도 기록
  • 출력 데이터 형식

    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
    단계 수 계산(초기 온도 포함)
    Tfix
    노드 온도 지정 없음(0: 온도 지정 없음, 1: 온도 지정)
    iii, ttime, Node_x
    출력 시간 절차, 시간, 온도 역사의 노드 번호 이후 각 시각의 온도 역사 기록

    수출 사례


    matplotlib으로 만든 온도 이력도



    GMT가 제작한 시간별 온도 분포도


    시간 단위는 hour이다.

    이상

    좋은 웹페이지 즐겨찾기