■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 | [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つの剛性マトリックスを作成する。 つまり、剛性マトリックスと右辺で連立方程式が成り立つ。 | なし |
FORM | Dirichlet境界条件を連立方程式へ埋め込む。 | なし |
SYSTEM | 連立方程式を解く | なし |
BANDWID | 要素せ節点番号情報から、連立方程式のバンド幅を計算する。 | なし |
CONVMX | Convective 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 |