Processing math: 100%
流体力学

【流束ヤコビ行列の導出】1次元圧縮性性オイラー方程式

こんにちは(@t_kun_kamakiri
本記事では、流体力学の基礎式を簡単にしたオイラー方程式に関する基礎理論のうちの流束ヤコビ行列の導出を解説します。
こちらは計算力学技術者の熱流体1級(単相流)で問われる知識です。
導出過程も含め物理的意味を押さえておくと良いです。

まず、本記事で示したい内容は以下です。

ここにタイトルを入力

システム方程式から流束ヤコビ行列を導出する

システム方程式

Qt+Fx=0

流束ヤコビ行列

A=FQ=(0103γ2u2(3γ)uγ1(γ12h)uh(γ1)u2γu)

Aを導出するのが本記事の内容になります。

システム方程式とは

圧縮性流れを扱う場合、レイノルズ数Re=ULνが大きいことが多く、粘性項目を無視したオイラー方程式での数値計算法の解法を基本として考えられています。

以下に示すのが1次元圧縮性の非粘性に対するオイラー方程式です。

質量保存則

ρt+ρux=0

運動量保存則

ρut+x(p+ρu2)=0

エネルギー保存則

et+x(e+p)u=0

  • ρ:密度
  • u:流速
  • p:圧力
  • e:単位体積当たりの気体の全エネルギー

これらをまとめると以下になります。

Qt+Fx=0

ここで、

Q=[ρρue]=[q1q2q3],F=[ρup+ρu2(e+p)u]=[f1f2f3]

スカラー輸送方程式(移流方程式など)に対して、(4)を総称してシステム方程式と呼ぶことがあります。流速がマッハ数に対して大きくなる場合、圧縮性流れとして現象をモデル化する必要が出てきますが、レイノルズ数も大きくなるため、粘性項を無視したオイラー方程式の解法を基本として計算方法が展開されているということが大事な点です。

流束ヤコビ行列の導入

(4)の特徴を理解するために以下のように式を変換します。
以下Aのような行列を導入し、行列の対角化を行うと現象を理解しやすくなるというメリットがありますが、本記事ではそこまでは踏み込みません。

Qt+FQQx=0

ここで、A=FQと置きます。
Aを流束ヤコビ行列と呼びます。

Aの成分は後で具体的に計算するとして、結果だけ先に示しておきましょう。

A=FQ=(0103γ2u2(3γ)uγ1(γ12h)uh(γ1)u2γu)

どのようにして(6)を導くのかを今から見ていきます。
Qt+Ft=0のままでは少しわかりにくいので、Q,Fは3つの成分をもつ縦ベクトルなので、

Qit+Fix=0

のように書くことにします。
そうすると、(7)はQit+jFiQjQjx=0となり、Aij=FiQjと置いていることになります。
jが2回続いているのでjは省略してQit+FiQjQjx=0と書くことが多いです。

わかりやすくi成分だけを書くとj=1,2,3だけを足すので、Qit+FiQ1Q1x+FiQ2Q2x+FiQ3Q3x=0となります。

行列で書くと以下のようになっていることからも理解できることでしょう。

t[ρup+ρu2(e+p)u]+(A11A12A13A21A22A23A31A32A33)x[ρρue]=0

つまり、Aの成分は以下のようになっています。

A=FQ=(F1Q1F1Q2F1Q3F2Q1F2Q2F2Q3F3Q1F3Q2F3Q3)=(A11A12A13A21A22A23A31A32A33)

例えば、A23=F1Q2=(e+p)u(ρu)と調子で計算していくことになります。

理想気体の状態方程式

理想気体の状態方程式pV=nRTを少し変形しておきましょう。
pV=nRTp=mVRMTと置くと、p=ρ˜RTとなります。
ここで、Mは気体のモル質量、ρ=mVは密度、˜R=RMは比気体定数です。

さらに、˜R=~Cp~Cv=~Cv(γ1)と置けば、理想気体の状態方程式はp=ρ~Cv(γ1)Tとなります。

Eは単位体積あたりの内部エネルギーE=rho~CvTとも書くことができ、単位体積あたりの全エネルギーe=E+12ρu2であるので、

p=(γ1)(e12ρu2)

となります。

  • 単位体積あたりの全エネルギーe=E+12ρu2
  • エンタルピーH=ρh
    単位体積当たりの全エンタルピーh=e+p

わざわざ気体定数Rに対して、比気体定数˜R=RM~Cp,~Cvなどチルダ(~)を頭につけたかというと、きっちり区別しておいた方がわかりやすいからです。
まず、気体定数Rや定圧比熱Cp、定積比熱Cvの単位は[J/K mol]ですが、˜R=RM~Cp,~Cvの単位は[J/K kg]です。どちらの単位が好まれて使われているかは、分野によりますが、単位が違うと値が異なります。

流束ヤコビ行列の成分の導出

さて、ようやくA=FQ=(A11A12A13A21A22A23A31A32A33)の計算に移ります。

f1q1=ρuq1=q2q1=0

f2q2=(ρu2+p)m2=((q2)2ρ+(γ1)[e(q2)22ρ])q2=q2((q2)2ρ(γ1)(q2)22ρ)=2q2ρ(γ1)q2ρ=(3γ)u

f3m3=((e+p)u)m3=((e+(γ1)[e(m2)22ρ])u)m3=m3((e+(γ1)[e(m2)22ρ])u)=γu

※計算途中は読者に任せるとして結果はこちらです。

流束ヤコビ行列

A=FQ=(0103γ2u2(3γ)uγ1(γ12h)uh(γ1)u2γu)

まとめ

本記事では、オイラー方程式に基づき、流束ヤコビ行列の導出を解説しました。以下に主要なポイントをまとめます。

1. システム方程式の基本構造

  • 圧縮性流体の流れをモデル化する際に用いられる基礎方程式であり、質量保存則、運動量保存則、エネルギー保存則の統合形として表現されます。
  • システム方程式:Qt+Fx=0
  • 状態変数ベクトルQと流束ベクトルFを具体化することで、圧縮性流れの各保存則に対応。

2. 流束ヤコビ行列の定義

  • 流束Fを状態変数 Qで偏微分した行列:A=FQ=(0103γ2u2(3γ)uγ1(γ12h)uh(γ1)u2γu)
  • 流束ヤコビ行列Aは、オイラー方程式をより簡単に扱える形へ変換し、流れの特性を表すために用いられます

参考書

本記事の計算過程やプログラムを詳しく解説している参考書を挙げておきます。

まずはPythonのプログラムとともに解説している参考書。

Pythonで学ぶ流体力学の数値計算法

Pythonで学ぶ流体力学の数値計算法

藤井孝藏, 立川智章
4,180円(04/03 18:22時点)
Amazonの情報を掲載しています

C言語(タイトルはC++とあるが、ほぼC言語です)の数値流体としてはめちゃくちゃわかりやすい参考書です。プログラムコードもすべて公開されているのでお勧めです。
ただし、難しいので何度も読んで理解することになりますが。

数値流体解析の基礎 - Visual C++とgnuplotによる圧縮性・非圧縮性流体解析 -

数値流体解析の基礎 – Visual C++とgnuplotによる圧縮性・非圧縮性流体解析 –

鋒, 肖, 孝夫, 長崎
3,960円(04/03 10:42時点)
Amazonの情報を掲載しています

理論の勉強ならやはりこちらの銀本。
数値流体の理論の本としてはとても有名で何度も読むことで理解が深まっていきます。

流体力学の数値計算法

流体力学の数値計算法

藤井孝蔵
4,743円(04/03 21:27時点)
発売日: 1994/03/23
Amazonの情報を掲載しています

計算力学技術者資格のための問題アプリ

計算力学技術者固体力学2級対策アプリをリリースしました。

  • 下記をクリックしてホームページでダウンロードできます。
  • LINE公式に登録すると無料で問題の一部を閲覧できます
    ※LINEの仕様で数式がずれていますが、アプリでは問題ありません。
関連記事もどうぞ

COMMENT

目次へ