[執筆途中]流路網計算の代数方程式を微分方程式にして解く

ページ内Link
＊執筆途中だが公開ポリシーに則り公開。
2023.04.19

記事を通して言いたいこと（結論）

　例によって、冒頭で述べておく。
  • 流路網（FluidFlowNetwork）の代数方程式の計算を，非定常項を追加することで収束計算ではなく，微分方程式の時系列計算に置き換えてとして解いた（nodeが2つで，4未知数4方程式となるシステム）．
  • 本記事では，上記の立式過程も示す（流体は理想気体を想定）．
  • 代数方程式を収束計算手法で解くのに比べて，初期値が悪くても流量収支が釣り合うする解にたどり着く事ができ，過渡応答を評価したい場合だけでなく，計算収束性を改善したい場合にも利用できる手法である可能性が高い．

モデル化対象（とその周辺について）

    流路抵抗３つ，ボリューム（Node）２つの直列流路

    非定常項を加えて，状態方程式（微分方程式）を立式

    　式はボリューム１つに対してのみ示す．複雑な流路網ではボリュームの数だけ状態変数と状態方程式が増える．

    Mass dynamics

    dm dt = ∑ ( m flow ) − . mass flow into volume is positive sign m = ρV dm dt = dρ dt ∗ V ρ = p RT dρ dt = 1 R d dt ( p T ) = 1 R p ˙ T − p T ˙ T 2 ∑ ( m flow ) = V RT p ˙ − Vp R T 2 T ˙ {dm} over {dt} = sum (m_flow) newline -. mass flow into volume is positive sign newline m=ρV newline {dm} over {dt} = {dρ} over {dt}*V newline ρ= {p} over {RT} newline {dρ} over {dt} = {1} over {R} {d} over {dt} ({p} over {T}) newline = = {1} over {R} { dot p T – p dot T } over { T^{2} } newline sum (m_flow) = {V} over {RT} dot p – {Vp} over {R T^{2}} dot T

    Energy dynamics

    dU dt = ∑ ( mh flow ) dU dt = d dt ( mu ) = dm dt u − m du dt u = h − pv = CpT − p ρ = CpT − RT du dt = d dt ( CpT − RT ) = ( Cp − T ) dT dt = ( Cp − R ) T ˙ dU dt = V R T 2 ( p ˙ T − p T ˙ ) ( Cp − R ) T + m ( Cp − R ) T ˙ ⋮ = V ( Cp − R ) R p ˙ + [ − Vp ( Cp − R ) RT + m ( Cp − R ) ] T ˙ ∑ ( mh flow ) = V ( Cp − R ) R p ˙ + ( Cp − R ) [ − Vp RT + m ] T ˙ {dU} over {dt} = sum (mh_flow) newline {dU} over {dt} = {d} over {dt}(mu) newline == {dm} over dt u – m {du} over {dt} newline u= h – pv newline == CpT- {p} over {ρ} newline == CpT- RT newline {du} over {dt}= {d} over {dt}(CpT-RT) newline == (Cp-T) {dT} over {dt} newline == (Cp-R) dot T newline {dU} over {dt} = {V} over {R T^2}( dot p T – p dot T )(Cp-R)T + m(Cp-R) dot T newline dotsvert newline == {V(Cp-R)} over {R} dot p + [-{Vp(Cp-R)} over {RT} + m(Cp-R)] dot T newline sum (mh_flow) = {V(Cp-R)} over {R} dot p + (Cp-R)[-{Vp} over {RT} + m] dot T

    状態方程式

    [ V RT − Vp R T 2 V ( Cp − R ) R ( Cp − R ) [ − Vp RT + m ] ] [ p ˙ T ˙ ] = [ ∑ ( m flow ) ∑ ( mh flow ) ] [ p ˙ T ˙ ] = [ V RT − Vp R T 2 V ( Cp − R ) R ( Cp − R ) [ − Vp RT + m ] ] − 1 [ ∑ ( m flow ) ∑ ( mh flow ) ] left [ matrix{{V} over {RT} # -{Vp} over {R T^{2}} ## {V(Cp-R)} over {R} # (Cp-R)[-{Vp} over {RT} + m]} right ] left [ binom {dot p}{dot T} right ] = left [ binom {sum (m_flow)}{sum (mh_flow)} right ] newline left [ binom {dot p}{dot T} right ] = left [ matrix{{V} over {RT} # -{Vp} over {R T^{2}} ## {V(Cp-R)} over {R} # (Cp-R)[-{Vp} over {RT} + m]} right ]^{-1} left [ binom {sum (m_flow)}{sum (mh_flow)} right ]

    定常状態を代数方程式として解く場合の独立・従属変数

    • 独立変数
      • Volume1, 圧力
      • Volume1, 温度
      • Volume2, 圧力
      • Volume2, 温度
    • 従属変数
      • Volume1, 質量流量収支
      • Volume1, 熱量流量収支
      • Volume2, 質量流量収支
      • Volume2, 熱量流量収支

シミュレーションモデル：

    Scilab Xcosで作る：

    　Xcosなら関数を記述したブロックを繋げて複雑な計算式となるモデルを比較的容易に組立てられる．
    　一方で，XcosはModelicaとは異なり，自動で数式処理をおこなう機能は無い（Coselicaと呼ばれるModelica言語Blockは別）．代数方程式/微分方程式で記述されるモデルを解けるようにするには，Inpuot/Outputが明示された処理block，積分器block，代数方程式ソルバblockを適切に繋げて処理を作る必要が有る．
    　その為，実は裏で”数式変形して方程式が宜しく解かれている”と言うことは起きない．代数ブロックの繋ぎ方がループになっていると（代数方程式をソルバで解く必要があるモデル），エラーが出て実行がなされない．これは，今回行いたい事には適している．

    微分方程式（状態方程式）を立式して時間積分を行うモデル

    Diagram

    ●シミュレーションモデル情報

    VolumeDynamics blockの内部

    代数方程式を立式してソルバBlockで解くモデル（検証用）

    Diagram

    ●シミュレーションモデル情報

シミュレーション実行

    Input

    　無し．固定の境界条件を与える．

    Variables

      微分方程式モデル

      1. Volume圧力
      2. Volume温度
      3. Volume質量収支
      4. Volume熱量収支
      5. 最終時刻でのVolume状態

      代数方程式モデル

      　時間が経過しても値は変化しないので，シミュレーション実行時間は短く設定．
      1. Volume圧力
      2. Volume温度
      3. Volume質量収支
      4. Volume熱量収支
      5. 最終時刻でのVolume状態
以上

コメント

タイトルとURLをコピーしました