Finite Element Formulation for Fluid

Fluid Dynamics
Finite Element Formulation for Fluid-10

■2次元流体解析のプログラムNSEQ8DD.FOR■
離散化および後処理に必要な計算方法が全て揃ったところで、流体解析のプログラムNSEQ8DD.FORを紹介します。プログラム名のNSEQは、Navier-Stokes EQuationsからとりましたが、前にも述べた通り、ベースになっている支配方程式は、Navier-Poisson eqations であることに注意して下さい。

プログラムの概略としては、これまでにも幾度となく述べてきましたが、ここでざーと特徴だけを述べることにしましょう。まず、

連続の式は、カップルしない。

支配方程式の時間項は、差分で離散化し、[N]^T[N]はlumpingしてあります。

λの項は、Reduced Integration technique で離散化する。

Convective terms は、連立方程式の右辺に足し込み、時間の不整合は、繰り返し計算で補う。

繰り返し計算のターミネーションは、計算結果の最大誤差(ERMAX)で判定する。ERMAXについては、変数名の説明を参照して下さい。

プログラムのアルゴリズムはラミナーNewtonian流体を扱えられ様になっているが、乱流や非Newtonian流体へ対応せるためにのプログラムの変更は用意出来るように設計されてある。

平均圧力は、速度u(x,y,t)とv(x,y,t)が得られた後に計算する。

では、プログラムの説明のために、Subroutines から紹介しましょう。下表にSubroutinesの名称と役割を述べておきますので、何が計算されているかを想像して下さい。

Subroutine名	役割	CALLさせるSubrountines
GRULE	Gauss Integration 座標を返す。	なし
DERIV	Gauss Integration 座標点で、形状関数の微分を差分で計算する。	SHAPEF
SHAPEF	Gauss Integration 座標点で、形状関数の値を計算する。	ISOPARA
ISOPARA	与えられた座標値での形状関数値を計算する。	なし
INPUT	解析に必要なデータをファイルから読む。	なし
GSM	Laplace式の[B]^T[B]を計算し、GlobalマトリックスRKXとRKYを生成する。	なし
GLAMBDA	支配方程式中の全てのλ項をReduced Integration Techniqueで計算し、GlobalマトリックスRKXとRKYとASYMXとASYMYへアッセンブルする。。	なし
GLAMBDA	Single Variableの[N]^T[N]を計算し、Globalマトリックス{GEOM}を生成する	なし
RHST0	連立方程式の右辺を生成する。	なし
SOLVEF	未知数u(x,y,t)とv(x,y,t)を計算する。	ASSEMB, FORM, SYSTEM,VELCHK
VELCHK	繰り返し計算のための誤差を計算する。誤差判定はMain で行う。	なし
ASSEMB	各式毎に作成されたglobal matrices を1つの剛性マトリックスを作成する。つまり、剛性マトリックスと右辺で連立方程式が成り立つ。	なし
FORM	Dirichlet境界条件を連立方程式へ埋め込む。	なし
SYSTEM	連立方程式を解く	なし
BANDWID	要素せ節点番号情報から、連立方程式のバンド幅を計算する。	なし
CONVMX	Convective Term を計算する。	なし
REWRITE	時間ステップが10進む毎に、Current time とu, v の値をファイルを書き込む。時間計算が再開されるときは、このファイルから初期値として計算が再開される。	なし
TIMESTP	過去の繰り返し回数をベースに、次の時間ステップ値(Δt)を計算する。	なし
BCASSY	未知数ごとの境界条件を、1つの変数にまとめる。	なし

BACK	NEXT