境界条件の設定2021年6月20日 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
はじめに境界条件の設定について。 使用バージョンOpenFOAM v2012、v8 境界タイプの設定OpenFOAM のメッシュファイル群 (constant/polyMesh) には boundary ファイルという境界情報を含むファイルがある。これは以下のような内容になっている。
5
(
inlet
{
type patch;
nFaces 30;
startFace 24170;
}
outlet
{
type patch;
nFaces 57;
startFace 24200;
}
upperWall
{
type wall;
nFaces 223;
startFace 24257;
}
lowerWall
{
type wall;
nFaces 250;
startFace 24480;
}
frontAndBack
{
type empty;
inGroups 1(empty);
nFaces 24450;
startFace 24730;
}
)
この情報はふつうメッシュ作成ソフトあるいはメッシュ変換ユーティリティが書き出すが、境界条件に合わせて境界のタイプを設定する必要がある。境界のタイプにはつぎのものがある。
パッチ (patch) というのは、入口や出口で使用されるタイプである。wall はそのまま壁条件であり、壁関数境界条件を設定する場合はこれを選択する必要がある。 symmetryPlane, symmetry は対称条件を与える面のタイプで、前者が平面、後者がそれ以外に適用されるものである (平面以外の対称面というのがよくわからんが)。基本的には前者を用いればよい。 cyclic, cyclicAMI については、対応する 2 面で設定する必要があり、それぞれ "neighbourPatch" で対応する境界を指定する。 2 次元問題の場合は、計算しない方向の両面に empty を指定する。2 次元軸対称については、実際の 3 次元モデルの 5°以下の楔形のモデルを作り、計算しない両面を wedge に設定すればよい。 v1706 以降では、オーバーセットメッシュ用の overset タイプがある。 境界条件の設定境界条件の設定は 0 ディレクトリ内の各フィールドファイルで行う。必要なフィールドファイルはソルバーごとに異なる。乱流モデルに k-εモデルを用いる場合、simpleFoam、buoyantBoussinesqSimpleFoam、buoyantSimpleFoam で必要になるフィールドファイルはそれぞれつぎのようになる。
v2.4.0 以前
※v2.4.0 以前では、nut, alphat はあってもよい (なくても自動で作られる)。ただし、buoyantSimpleFoam では nut ではなく mut。 ただし、圧力の単位に注意する。simpleFoam と buoyantBoussinesqSimpleFoam の圧力は密度で割ったものであるが、buoyantSimpleFoam の圧力は本来の圧力である。 基本設定U についてフィールドファイルの基本設定を示す。 U
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
object U;
}
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
inlet
{
type fixedValue;
value uniform (1 0 0);
}
outlet
{
type zeroGradient;
}
wall
{
type fixedValue;
value uniform (0 0 0);
}
}
ヘッダーヘッダーは「おまじない」と思えばよい。
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
object U;
}
フィールドデータの場合、class には以下のようなフィールドのクラスを書くようになっている。
dimensionsdimensions で単位を設定する。次のような形式で各単位の指数を指定する。 [kg m s K mol A Cd] "[0 1 -1 0 0 0 0]" は m/s を表している。つぎのような表現も可能である。 dimensions [m/s]; internalFieldフィールド内部の値を設定する。0 ディレクトリの場合、初期値を設定する。 boundaryField境界条件を設定する。境界条件はつぎのような形式で設定する。
inlet
{
type fixedValue;
value uniform (10 0 0);
}
メッシュの boundary に記述された境界名に対して、境界条件のタイプと設定を設定する。基本的な境界条件タイプを以下に挙げる。
特別な場合を除いて、基本的には type と value が 2 つとも必要である。value は境界条件の設定に使われる場合も使われない場合もあり、初期値として使われる場合もある。 境界名では、二重引用符で囲んで "..." のような書き方をすると、正規表現を使うことができる (v1.6 以降)。"inlet.*" と書くと "inlet-1" や "inlet-2" などに当てはまり、"(inlet-1|inlet-2)" と書くと "inlet-1" か "inlet-2" に当てはまる。".*" と書くとすべての境界に当てはまる。正規表現と直接指定が同時にある場合は、順番によらず直接指定が優先されるため、".*" で基本的な設定をしておき、それ以外は個別に設定するような書き方ができる。たとえば、次のような内容のファイルをテンプレートして用意しておけばよいかもしれない。
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
object U;
}
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
#includeEtc "caseDicts/setConstraintTypes"
".*"
{
type fixedValue;
value uniform (0 0 0);
}
}
#includeEtc 文 (v2.4.0 以降) では、symmetryPlane や empty などの拘束型 (constraint) の境界条件の設定を読み込み、これらの設定を記述しなくてもよいようにしている。これは次のように書いたものと同じである。
#include "${WM_PROJECT_DIR}/etc/caseDicts/setConstraintTypes"
代表的なフィールドの基本設定代表的なフィールドの基本設定を以下に示す。"基本設定" の中身は、type と value を表しており、追加のパラメタがあればそれを示している。 非圧縮性共通
※v4 以降は U の静止壁条件 (non-slip 条件) に noSlip が使える。 v2.4.0 以前
圧縮性共通
v2.4.0 以前
buoyantBoussinesqSimpleFoam
buoyantSimpleFoam
ここで、R (レイノルズ応力) の volSymmTensorField は対称テンソル場で、対称テンソルは 6 成分からなり、各成分を (xx xy xz yy yz zz) の順番で指定する。 0 ディレクトリのバックアップ0 ディレクトリの中身はたまに標準ユーティリティなどにより書き換えられてしまうことがあるため、基本的な設定をしたあとに、"0.org" などといった名前でオリジナルをコピーしておくとよい。 $ cp -r 0 0.org 入口・出口条件単純な条件単純な入口・出口条件は、つぎのようになる。 入口
出口
逆流を考慮した出口条件つぎのような設定により、逆流を考慮した出口条件を設定できる。ただし、逆流は計算上好ましくないので、極力逆流が起こらないように出口を設定したほうがよい。
pressureInletOutletVelocity は圧力から入出力速度を決めるもので、つぎのように設定する。
outlet
{
type pressureInletOutletVelocity;
value uniform (0 0 0);
}
totalPressure は全圧を指定するもので、つぎのように設定する。
outlet
{
type totalPressure;
value uniform 0;
gamma 1.4;
p0 uniform 0;
}
p0 で圧力を指定、gamma は比熱比である。 浮力を考慮した圧力を指定したい場合は、prghPressure を用いる。
outlet
{
type prghPressure;
p uniform 101325;
value uniform 101325;
}
浮力で流れが上部境界から抜けていく場合などは、この境界条件を用いる。 inletOutlet はつぎのように設定する。
outlet
{
type inletOutlet;
value uniform 300;
inletValue uniform 300;
}
流速の条件流入速度を指定するために fixedValue では不便な場合がある。以下に fixedValue 以外の流入境界条件を挙げる。 面の法線方向の流速の指定
inlet
{
type surfaceNormalFixedValue;
refValue uniform 10; // [m/s]
}
refValue で面の法線方向の流速の大きさを指定。 メッシュによっては面の方向が領域の内側ではなく外側に向いている場合があり、注意が必要。 体積流量の指定
inlet
{
type flowRateInletVelocity;
volumetricFlowRate constant 2.54e-4; // [m3/s]
value uniform (0 0 0);
}
質量流量の指定
inlet
{
type flowRateInletVelocity;
massFlowRate constant 2.54e-4; // [kg/s]
rhoInlet 1; // [kg/m3]
value uniform (0 0 0);
}
非圧縮性流体ソルバーの場合、rhoInlet で流入する流体の密度を指定する。圧縮性ソルバーの場合は rho の値が用いられる。 乱流流入条件
inlet
{
type turbulentInlet;
referenceField uniform (10 0 0); // 流速 [m/s]
fluctuationScale (0.02 0.01 0.01); // 変動スケール
}
変動を含んだ流入条件を設定する。 乱流諸量の条件入口の $k$ や $\varepsilon$、$\omega$ などにどのような値を指定すべきかは、一般的にははっきりしない。実験値が得られない場合、適当に見積もる必要がある。$k$、$\varepsilon$、$\omega$ はそれぞれ次式で表される。 $$k = \frac{3}{2} (UI)^2$$ $$\varepsilon = \frac{C_{\mu}^{3/4} k^{3/2}}{L_m}$$ $$\omega = \frac{\varepsilon}{C_\mu k} = \frac{k^{1/2}}{C_{\mu}^{1/4} L_m}$$ $$C_{\mu} = 0.09$$※$\omega$ の $C_{\mu}$ は、乱流粘性が $k$-$\varepsilon$ と $k$-$\omega$ で一致するとしたときに出てくる。 $I$ は乱流強度 (turbulent intensity)、$L_m$ は混合長 (mixing length) である。乱流強度は、十分に発達した流れでは数％の値をとる。混合長は、ダクト内の十分発達した流れでは次式で近似できる。 $$L_m = 0.07 L$$$L$ は代表長さで、円形のダクトであれば直径、そうでなければ水力直径 (4×断面積/断面周長) とする。 レイノルズ応力 $R$ については、対角成分に $(3/2)k$ を入れておけばよい。 $k$ や$\varepsilon$ などを直接指定するのではなく、乱流強度や混合長で指定するための特別な境界条件がある。 $k$
inlet
{
type turbulentIntensityKineticEnergyInlet;
intensity 0.05; // 5 %
value uniform 1;
}
intensity に乱流強度を設定する。乱流強度が 5 % ならば 0.05 を指定する。 $\varepsilon$
inlet
{
type turbulentMixingLengthDissipationRateInlet;
mixingLength 0.005; // [m]
value uniform 1;
}
mixingLength に混合長を設定する。 $\omega$
inlet
{
type turbulentMixingLengthFrequencyInlet;
mixingLength 0.005; // [m]
k k;
value uniform 1;
}
mixingLength に混合長を設定する。 v2.4.x 以前では、buoyantSimpleFoam などの圧縮性流体ソルバーの場合、εとωについては、タイプ名の頭に "compressible::" を付ける (たとえば "compressible::turbulentMixingLengthDissipationRateInlet")。 壁面速度の設定壁面速度の設定は fixedValue でそのまま指定すればよい。いくつか特殊な設定について以下で述べる。 スリップ壁面をスリップ条件にしたい場合は、U で slip を指定する。
wall
{
type slip;
}
それ以外の変数については zeroGradient にする。 回転速度の設定円筒の壁が回転するような場合は、つぎのようにして回転速度を設定する。
wall
{
type rotatingWallVelocity;
origin (0 0 0); // 回転軸の位置
axis (0 0 1); // 回転軸ベクトル
omega constant 10; // 角速度 [rad/s]
}
速度 v を指定したい場合、円筒半径を r とすると、omega = v/r を指定すればよい。 熱の境界条件の設定熱 (温度) の境界条件としては、断熱である zeroGradient、温度指定である fixedValue 以外に以下のようなものがある。 熱流束の設定壁面に熱流束を設定するには externalWallHeatFluxTemperature を用いる。
wall
{
type externalWallHeatFluxTemperature;
mode flux;
kappaMethod fluidThermo;
q uniform 100; // 熱流束 [W/m2]
value uniform 300; // 初期温度 [K]
}
mode で "power" を指定すれば、熱流束 [W/m2] の代わりに熱量 [W] を指定できる。
wall
{
type externalWallHeatFluxTemperature;
mode power;
kappaMethod fluidThermo;
Q 100; // 熱量 [W]
value uniform 300;
}
externalWallHeatFluxTemperature は buoyantBoussinesqSimpleFoam などでは利用できない。Boussinesq 近似と併用したい場合は、buoyantSimpleFoam などで状態方程式に Boussinesq を用いる (v4 以上)。 turbulentHeatFluxTemperaturev4 以前であれば、turbulentHeatFluxTemperature という境界条件も使える。 buoyantSimpleFoam などではつぎのように設定する。
wall
{
type compressible::turbulentHeatFluxTemperature;
heatSource flux; // flux (熱流束) または power (熱量)
q uniform 100; // 熱流束 [W/m2] または 熱量 [W]
kappaMethod fluidThermo;
kappa none;
value uniform 300;
}
v3.0 以前だと、kappaMethod, kappa の部分は以下のように書く。
kappa fluidThermo;
kappaName none;
buoyantBoussinesqSimpleFoam などではつぎのように設定する。
wall
{
type turbulentHeatFluxTemperature;
heatSource flux;
q uniform 100;
alphaEff alphaEff;
value uniform 300;
}
ただし、constant/transportProperties で Cp0 (密度×比熱) を設定する必要がある。 壁面熱伝達の設定壁面熱伝達を設定するには、externalWallHeatFluxTemperature の mode で "coefficient" を指定する。
wall
{
type externalWallHeatFluxTemperature;
mode coefficient;
kappaMethod fluidThermo;
Ta constant 500; // 外部温度 [K]
h constant 10; // 熱伝達係数 [W/m2-K]
value uniform 300;
}
外部の輻射の影響を考慮したい場合は、emissivity で放射率を指定する。
emissivity 0.8; // 放射率
また、熱伝達率の異なる板が複数枚積層している場合を表すこともできる。
thicknessLayers (0.1 0.2 0.3 0.4); // 厚み [m]
kappaLayers (1 2 3 4); // 熱伝導率 [W/m-K]
また、wallHeatTransfer という境界条件もある。
wall
{
type wallHeatTransfer;
Tinf uniform 500; // 外部温度 [K]
alphaWall uniform 10; // 熱伝達係数 [m/m2-K]
value uniform 300;
}
これは buoyantBoussinesqSimpleFoam などでは利用できない。 表による境界条件の設定fixedProfile により表形式で分布を持った境界条件を表現できる。
inlet
{
type fixedProfile;
profile csvFile;
profileCoeffs
{
file "k.csv";
nHeaderLine 0;
refColumn 0;
componentColumns (1);
separator ",";
mergeSeparators no;
}
direction (0 0 1);
origin 0;
}
ここでは profile に "csvFile" を指定することで、CSV ファイルを指定するタイプを選択している。profileCoeffs の file で CSV ファイルを指定している。CSV ファイルの中身は次のようになっている。 k.csv 0.01,0.327 0.013,0.367 0.025,0.405 ... 式による境界条件の設定OpenCFD 版では、exprFixedValue により式による境界条件を表現できる。
inlet
{
type exprFixedValue;
valueExpr "-10*face()/area()";
value uniform (0 0 0);
}
valueExpr で式を記述する。これは #eval と同じ形式を受け付けるが、ここでは境界特有の関数を用いることができる。上の例では、面の法線方向に大きさ 10 の速度を設定している。face() は境界面の面積ベクトル、area() はその面積である。マイナスが付いているのは、境界面の面積ベクトルは領域の外側を向いているためである。
inlet
{
type exprFixedValue;
h 0.0254;
valueExpr "-10*sin(pi()*pos().y()/$h)*face()/area()";
value uniform (0 0 0);
}
流速分布を与えることもできる。pos().y() は面の y 座標である。 非定常計算で速度の時間変動を与えることもできる。
inlet
{
type exprFixedValue;
T 100;
valueExpr "-10*sin(2*pi()*time()/$T)*face()/area()";
value uniform (0 0 0);
}
time() は時刻である。 | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| PENGUINITIS | ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||