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
GRULEGauss Integration 座標を返す。なし
DERIVGauss Integration 座標点で、形状関数の微分を差分で計算する。SHAPEF
SHAPEFGauss Integration 座標点で、形状関数の値を計算する。ISOPARA
ISOPARA与えられた座標値での形状関数値を計算する。なし
INPUT解析に必要なデータをファイルから読む。なし
GSM[B]T[D][B]をμ=μ, λ=0で要素毎の剛性マトリックス計算し、RKを生成する。なし
GSMμ=0,λ=λとReduced Integration Techniqueで要素毎の剛性マトリックス計算し、RKに足しこむ。なし
GDM[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つの剛性マトリックスを作成する。 つまり、剛性マトリックスと右辺で連立方程式が成り立つ。なし
FORMDirichlet境界条件を連立方程式へ埋め込む。なし
SYSTEM連立方程式を解くなし
BANDWID要素せ節点番号情報から、連立方程式のバンド幅を計算する。なし
CONVMXConvective Term を計算する。なし
REWRITE時間ステップが10進む毎に、Current time とu, v の値をファイルを書き込む。時間計算が再開されるときは、このファイルから初期値として計算が再開される。なし
TIMESTP過去の繰り返し回数をベースに、次の時間ステップ値(Δt)を計算する。なし

上の表に示すようにSubroutine GSMは全ての項(μ項とλ項)を含んだ剛性マトリクスを計算するようになったいます。ですから、λ項のみを低減精度積分を行う 場合は、まず、μ=μとλ=0と正規積分のGauss-Legendre座標を使いSubroutine GSMで剛性マトリクスを計算し変数RKに格納します。その後、μ=0とλ=λと 低減精度積分のGauss-Legendre座標を使いSubroutine GSMで剛性マトリクスを計算し変数RKに加算しています。

BACK NEXT
Menu Govern FEM 2nd Visco DCF Cylinder /Heat Lub