1次元熱拡散方程式の平衡解に対する有限要素法解析(コーディング編)¶

実際に１次元熱拡散方程式の定常解を解いてみる¶

前節では，１次元熱拡散方程式の定常解を数値的に解くために，有限要素法における定式化を行った．本節では，実際にコーディングを施し，解析を行ってみる．

主なプログラムの流れ，及び，メインプログラムは以下となる．

program main
  use variablesMod, only : Kmat, phi, Fvec, nNode, MtrxFile, RsltFile
  use initiatorMod, only : initProblem
  use genMatrixMod, only : gen_Kmatrix, gen_Fvector
  use matrixInvMod, only : gaussElimin
  use ioUtilityMod, only : writeResult, writeMatrix
  implicit none

  call initProblem
  call gen_KMatrix
  call gen_fVector
  call writeMatrix( Kmat, phi, Fvec, nNode, trim(MtrxFile) )
  call gaussElimin( Kmat, phi, Fvec, nNode )
  call writeResult( trim(RsltFile) )
  
end program main

定常問題の有限要素解析は，連立方程式 $Ax=b$ に帰着する．プログラムも行列・左辺生成→行列反転が基本となる．

以後，各サブルーチンについてみていく．

変数割付・初期化¶

単純にallocate文によって変数を割付けし，ゼロクリアによる初期化を行う準備部分のサブリーチンとしている．変数を節点数分だけ用意しているので，問題領域を有限要素へ分割しているプログラム領域に相当するとも考えることができる．

module initiatorMod
contains

  ! ===================================== !
  ! === Initiatialize Variables       === !
  ! ===================================== !
  subroutine initProblem
    use variablesMod
    implicit none

    ! ------------------------------------- !
    ! --- [1] Allocate & Initialize     --- !
    ! ------------------------------------- !
    allocate( phi(nNode), Fvec(nNode), Fvec_S(nNode), Fvec_V(nNode) )
    allocate( Kmat(nNode,nNode) )
    phi(:)    = 0.d0
    Fvec(:)   = 0.d0
    Fvec_S(:) = 0.d0
    Fvec_V(:) = 0.d0
    Kmat(:,:) = 0.d0
    
    return
  end subroutine initProblem

end module initiatorMod

行列・右辺ベクトル生成¶

ここでは主に次の２つの操作を行う．

基礎方程式系からのK-Matrix, F-Vectorの生成
境界条件の適用

(1)有限要素法で解析したい基礎方程式系及び入力情報は，行列(K-Matrix)及び右辺ベクトル(F-Vector)として記述される．ここでは，K-Matrix及びF-Vectorをプログラム上に表現することで問題を記述する． (2)また，境界条件を適用することでK-Matrix, F-Vectorに変更を加える．解析したい方程式系や境界条件によってK-Matrix, F-Vectorが左右されるため，これら２つの工程が有限要素法における問題定義を行うプログラム領域となる．

実際に，それぞれの要素について，matrixを計算してみる．ここでは，例として，節点数5の系について考えてみる．

1次線要素に関するk-Matrix $\{k\}^{e}$ ，及び，f-Vector の面積分成分 $\{f_S\}^{e}$ ，体積分成分 $\{f_V\}^{e}$ は前項で導出しているため，各要素が満たすべき方程式は次のように表される．

(e)=1 に対して

$\dfrac{ \lambda A }{ L } \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} \begin{bmatrix} \phi_1 \\ \phi_2 \end{bmatrix} = \dfrac{ \hat{Q} A L }{ 2 } \begin{bmatrix} 1 \\ 1 \end{bmatrix}$

(e)=2 に対して

$\dfrac{ \lambda A }{ L } \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} \begin{bmatrix} \phi_2 \\ \phi_3 \end{bmatrix} = \dfrac{ \hat{Q} A L }{ 2 } \begin{bmatrix} 1 \\ 1 \end{bmatrix}$

(e)=3 に対して

$\dfrac{ \lambda A }{ L } \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} \begin{bmatrix} \phi_3 \\ \phi_4 \end{bmatrix} = \dfrac{ \hat{Q} A L }{ 2 } \begin{bmatrix} 1 \\ 1 \end{bmatrix}$

(e)=4 に対して

$\dfrac{ \lambda A }{ L } \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} \begin{bmatrix} \phi_4 \\ \phi_5 \end{bmatrix} = \dfrac{ \hat{Q} A L }{ 2 } \begin{bmatrix} 1 \\ 1 \end{bmatrix} - qA \begin{bmatrix} 0 \\ 1 \end{bmatrix}$

これらは各要素に関係する接点の情報のみを含んだ方程式から成り立っている．上記方程式をまとめることによって，接点温度の情報を総括したGlobal Vector $\Phi$ の形で書くと，

$\dfrac{ \lambda A }{ L } \begin{bmatrix} 1 & -1 & 0 & 0 & 0 \\ -1 & 2 & -1 & 0 & 0 \\ 0 & -1 & 2 & -1 & 0 \\ 0 & 0 & -1 & 2 & -1 \\ 0 & 0 & 0 & -1 & 1 \end{bmatrix} \begin{bmatrix} \phi_1 \\ \phi_2 \\ \phi_3 \\ \phi_4 \\ \phi_5 \end{bmatrix} = \dfrac{\hat{Q}AL}{2} \begin{bmatrix} 1 \\ 2 \\ 2 \\ 2 \\ 1 \end{bmatrix} - qA \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 1 \end{bmatrix}$

このうち，右辺第二項は面積分した熱流束の流入量を表し，自由境界条件を適用している．一方で，左端において， Dirichlet 境界条件が課されており， $T(x=x_{Min})=0$ である．これは，陽に境界条件として K-Matrix, F-vectorに盛り込む必要がある．結局，

$\dfrac{ \lambda A }{ L } \begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ -1 & 2 & -1 & 0 & 0 \\ 0 & -1 & 2 & -1 & 0 \\ 0 & 0 & -1 & 2 & -1 \\ 0 & 0 & 0 & -1 & 1 \end{bmatrix} \begin{bmatrix} \phi_1 \\ \phi_2 \\ \phi_3 \\ \phi_4 \\ \phi_5 \end{bmatrix} = \begin{bmatrix} 0 \\ \hat{Q}AL \\ \hat{Q}AL \\ \hat{Q}AL \\ \dfrac{\hat{Q}AL }{2}-qA \end{bmatrix}$

module genMatrixMod
contains

  ! ===================================== !
  ! === generate [K] Matrix           === !
  ! ===================================== !
  subroutine gen_Kmatrix
    use variablesMod
    implicit none
    integer          :: i
    double precision :: coef

    ! ------------------------------------- !
    ! --- [1] Prepartion                --- !
    ! ------------------------------------- !
    coef = lambda * Area / Lelem
    
    ! ------------------------------------- !
    ! --- [2] k-Matrix Generation       --- !
    ! ------------------------------------- !
    ! -- [2-1] k-Matrix initialize      --  !
    Kmat(:,:) = 0.d0
    ! -- [2-2] row=1     (Left Edge  )  --  !
    Kmat(1,1) = + 1.d0
    Kmat(1,2) = - 1.d0
    ! -- [2-3] intermedite Region       --  !
    do i=2, nNode-1
       Kmat(i,i-1) = -1.d0
       Kmat(i,i  ) = +2.d0
       Kmat(i,i+1) = -1.d0
    enddo
    ! -- [2-4] row=nElem (Right Edge )  --  !
    Kmat(nNode,nNode-1) = - 1.d0
    Kmat(nNode,nNode  ) = + 1.d0
    ! -- [2-5] coefficient              --  !
    Kmat(:,:) = coef * Kmat(:,:)
    
    ! ------------------------------------- !
    ! --- [3] Boundary Condition        --- !
    ! ------------------------------------- !
    !  -- [3-1] Dirichlet Boundary      --  !
    Kmat(1,1)       = 1.d0
    Kmat(1,2:nNode) = 0.d0

    return
  end subroutine gen_Kmatrix


  ! ===================================== !
  ! === generate {F} vector           === !
  ! ===================================== !
  subroutine gen_Fvector
    use variablesMod
    implicit none
    integer          :: i
    double precision :: coefS, coefV

    ! ------------------------------------- !
    ! --- [1] Preparation               --- !
    ! ------------------------------------- !
    coefV     = 0.5d0 * dQdt * Area * Lelem
    coefS     = qedge * Area
    Fvec_S(:) = 0.d0
    Fvec_V(:) = 0.d0
    Fvec(:)   = 0.d0

    ! ------------------------------------- !
    ! --- [2] {F}(surface) term         --- !
    ! ------------------------------------- !
    !  -- [2-1] @ Free Boundary Cond.   --  !
    Fvec_S(nNode) = 1.d0
    !  -- [2-2] coefficient             --  !
    Fvec_S(:)     = coefS * Fvec_S(:)
    
    ! ------------------------------------- !
    ! --- [3] {F}(volume)  term         --- !
    ! ------------------------------------- !
    !  -- [2-1] @ Free Boundary Cond.   --  !
    Fvec_V(    1) = 1.d0
    do i=2, nNode-1
       Fvec_V(i) = 2.d0
    enddo
    Fvec_V(nNode) = 1.d0
    !  -- [2-2] coefficient             --  !
    Fvec_V(:)     = coefV * Fvec_V(:)

    ! ------------------------------------- !
    ! --- [4] total R.H.S.              --- !
    ! ------------------------------------- !
    Fvec(:) = Fvec_S(:) + Fvec_V(:)

    ! ------------------------------------- !
    ! --- [5] Boundary Condition        --- !
    ! ------------------------------------- !
    !  -- [5-1] Dirichlet Boundary      --  !
    Fvec(1) = Tedge
    
    return
  end subroutine gen_Fvector

  
end module genMatrixMod

行列反転¶

行列と右辺ベクトルが設定されれば，あとは解くだけである．有限要素解析として定式化された熱拡散方程式の定常解は，連立方程式 $Ax=b$ における行列反転 $x=A^{-1}b$ によって得られる．ここでは行列反転の手法として簡便に Gauss の消去法を用いる．大学教養レベルの線形代数の基本的な知識で理解でき，そのままコーディングするだけで行列反転が行える． Gauss の消去法のアルゴリズムは以下，及び，別記事に記載しておく．

module matrixInvMod
contains

  ! ========================================================== !
  ! === Gauss Elimination Solver                           === !
  ! ========================================================== !
  subroutine gaussElimin( Amat, xvec, bvec, nSize )
    implicit none
    integer         , intent(in)  :: nSize
    double precision, intent(in)  :: Amat(nSize,nSize)
    double precision, intent(in)  :: bvec(nSize)
    double precision, intent(out) :: xvec(nSize)
    integer                       :: i, j, k
    double precision              :: Dinv
    double precision, parameter   :: eps = 1.d-10
    double precision, allocatable :: Umat(:,:), vvec(:,:)

    ! ----------------------------------------- !
    ! --- [1] Preparation                   --- !
    ! ----------------------------------------- !
    !  -- [1-1] allocate Umat               --  !
    allocate( Umat(nSize,nSize) )
    !  -- [1-2] Copy AMat & bvec            --  !
    Umat(:,:) = Amat(:,:)
    xvec(:)   = bvec(:)
    
    ! ----------------------------------------- !
    ! --- [2] Upper Triangular Matrization  --- !
    ! ----------------------------------------- !
    do k=1, nSize

       !  -- [2-1] Diagonal Component       --  !
       if ( abs( Umat(k,k) ).lt.eps ) then
          write(6,*) '[gaussElimin] Amat != Regular Matrix :: No Solution End :: @ k= ', k
          stop
       endif
       Dinv      = 1.d0 / Umat(k,k)
       Umat(k,k) = 1.d0

       !  -- [2-2] Non-Diagonal Component   --  !
       if ( k.eq.nSize ) then
          ! -- [    Last Row :: k == nSize ] -- !
          xvec(k) = Dinv * xvec(k)
       else
          ! -- [Not Last Row :: k != nSize ] -- !
          !  - Division    -  !
          Umat(k,k+1:nSize) = Dinv * Umat(k,k+1:nSize)
          xvec(k)           = Dinv * xvec(k)
          !  - subtraction -  !
          do j=k+1,nSize
             Umat(j,k+1:nSize) = Umat(j,k+1:nSize) - Umat(j,k) * Umat(k,k+1:nSize)
             xvec(j)           = xvec(j)           - Umat(j,k) * xvec(k)
             Umat(j,k)         = 0.d0
          enddo
       endif
       
    end do

    ! ----------------------------------------- !
    ! --- [3] Backward Substituition        --- !
    ! ----------------------------------------- !
    do k=nSize-1, 1, -1
       do i=nSize, k+1, -1
          xvec(k) = xvec(k) - Umat(k,i)*xvec(i)
       enddo
    enddo

    return
  end subroutine gaussElimin

end module matrixInvMod

結果出力 & プロット例¶

結果の出力ルーチンの一例を示しておく．出力形式は好みであるが，今回の例では1次元問題であり，解として得られる情報も節点温度のみであるため，テキストデータとして出力している．

module ioUtilityMod
contains

  ! ===================================== !
  ! === write out Result :: phi       === !
  ! ===================================== !
  subroutine writeResult( FileName )
    use variablesMod
    implicit none
    integer                  :: i, lun
    double precision         :: dx
    character(*), intent(in) :: FileName
    
    ! ------------------------------------- !
    ! --- [1] Preparation / Open File   --- !
    ! ------------------------------------- !
    if ( Flag__writeFile ) then
       lun = 50
       open(lun,file=trim(FileName),status='replace',form='formatted')
    else
       lun = 6
    endif
    dx = xMax / dble( nNode-1 ) + xMin
    ! ------------------------------------- !
    ! --- [2] write out                 --- !
    ! ------------------------------------- !
    do i=1, nNode
       write(lun,*) dx*dble(i-1)+xMin, phi(i)
    enddo
    ! ------------------------------------- !
    ! --- [3] Open File                 --- !
    ! ------------------------------------- !
    if ( Flag__writeFile ) close(lun)

    return
  end subroutine writeResult


  subroutine writeMatrix( Amat, xvec, bvec, nSize, MtrxFile )
    implicit none
    integer         , intent(in) :: nSize
    double precision, intent(in) :: Amat(nSize,nSize), xvec(nSize), bvec(nSize)
    character(*)    , intent(in) :: MtrxFile
    integer                      :: i, j
    integer                      :: lun = 50
    character(30)                :: fmt = '(f8.4,1x)'

    open(lun,file=trim(MtrxFile),status='replace',form='formatted')
    do i=1, nSize
       do j=1, nSize
          write(lun,trim(fmt),advance='no') Amat(i,j)
       enddo
       write(lun,trim(fmt),advance='yes') bvec(i)
    enddo
    close(lun)
    return
  end subroutine writeMatrix

  
end module ioUtilityMod

作成した有限要素法サンプルの実行結果を以下に示す． [0,1.0]の区間を31分割の有限要素法で得られた節点温度をグラフ中，赤色十字(+)マークでプロットしている．対して，前述した解析解

$T = - \dfrac{ \dot{Q} }{ \lambda } [ x ( \dfrac{1}{2} x - x_{max} ) ]$

から得られる値を青実線で示している．有限要素法による解が解析解の良い近似になっていることがわかる．ここまでが，1次元有限要素解析の定式化・コーディング・実行結果である．

1次元熱拡散方程式の平衡解に対する有限要素法解析(コーディング編)¶

実際に１次元熱拡散方程式の定常解を解いてみる¶

変数割付・初期化¶

行列・右辺ベクトル生成¶

行列反転¶

結果出力 & プロット例¶

目次

前のトピックへ

次のトピックへ

このページ