reactingFoamにおけるdQについて

[Shintaro]

こんばんはshintaroと申します。

reactingFoamのチュートリアルcounterFlowFlame2DでmaxCoの値を変えた時の計算結果への影響を調べています。
maxCoを1, 0.4, 0.001と小さくしていったところ、温度の最大値が小さくなり、同時にdQは大きくなっていきました。

そこで疑問なのですが、dQは単位が[W]であるので発熱量だと思うのですが、なぜ発熱量が大きくなっているのに温度の最大値は小さくなっていくので
しょうか？

また、maxCoを小さくするのには限界があり、みなさんはどのような値をお使いなのでしょうか？

まだまだ知識が足りず申し訳ないのですが、ご教授いただけると幸いです。

[Shintaro]

shintaroです
追加で質問させてください。

dQの単位ですが、出力ファイルでは[W]だと思うのですが、basicChemistryModel.H中では[m2/s3]と書いてあります。
どちらがあっているのでしょうか？
この疑問にもお答えいただけたら幸いです。

[ohbuchi]

こんにちは。
私も同様の問題を確認しました。maxCo を0.4,0.2,0.1と
変化させたところ、dQ が、約２倍、4倍になり、他の変数は殆ど
変わりませんでした。本来、余り変化しないのが正解だと思います。
バグでしょうか？
dQ の単位は、W=J/s=kg m2/s3です。
ソースのコメントにkgが抜けているだけだと思います。

[Shintaro]

返信ありがとうございます。
温度はあまり変化しませんでしたか？
maxCoの変化に応じて変化してしまうのがdQだけであれば最悪よいと考えていますが、
私の場合温度が約2200Kから約2000K以下にまで変化してしまったので困っておりました。
（温度の具体的な変化については質問時に書洩らしました。申し訳ありません）

dQの単位の件についても返信ありがとうございます。
私も最初はソースのミスだと考えたのですが、
ソース basicChemistryModel.H のコメントには「enthalpy/sec」とあります。
エンタルピーの単位は[J/kg]ですので、やはりdQの単位は[J/kg /s=kg m2/s2 /kg /s=m2/s3]
となり、コメントが間違いだとは一概には言えないと思い悩んでいました。
また、発熱であればエンタルピー式 hEqn.H の生成項である Sh [kg/m/s3]が存在しており
わざわざdQを定義しているのは何か意図があるのでは？とも考えておりました。


以上の疑問を解決するにはやはりソースコードを読み解くしかなさそうですね…
できればソースコードでdQについて記述した箇所はご存じないでしょうか？
もしご存じでしたらお教えいただけると幸いです。

[ohbuchi]

チュートリアルではODEchemistryModelを使っているのではないかと思います。
ソースのdQメソッドを読むと、
　dQ=mesh_.V()*Sh()
で計算しています。それで、コメントでkgが抜けているのだと思いました。

こうなると、Shの計算が怪しく思えてきます。Sh()は
　Sh=specieThermo_.Hc()*RR_
で計算しています。RRは反応速度[kg/m3/s]の様です。
Hcは化学種のエンタルピー[J/kg]です。これで、Shの単位は[J/m3/s]
=[kg/m/s3]となります。

もしdQにバグがあるとすると、RRの計算に誤りがあることになりますが
化学種の濃度値などはmaxCoに依存していない様なので本質的な部分に
誤りはなさそうです。
dQがメインの計算ループで陽に使われていなそうなので、単純にdQを
計算する関数にバグがあるのだろうと気軽に考えていましたが、なかなか
見つかりません。

[Shintaro]

本当にご丁寧な解説ありがとうございます。

Shは hEqn.H の生成項なので単位は W/m3で間違いないと思います。
そうするとdQはW/m3×m3=Wで間違いなさそうですね。
やはりコメントが間違っているのでしょう。

reactingFoamで用いられているPaSR法についてのソースコードを見てはいませんが、
その中に根本的なバグが隠れているのかもしれませんね。
そのあたりを読み解いてみたいと思います。


また、今回の質問に関係なくて申し訳ないのですが...
ODEchemistryModel.Cの以下の部分でバグではないか？と考えている箇所があります。

template<class CompType, class ThermoType>
void Foam::ODEChemistryModel<CompType, ThermoType>::derivatives
(
const scalar time,
const scalarField &c,
scalarField& dcdt
) const
{
const scalar T = c[nSpecie_];
const scalar p = c[nSpecie_ + 1];

dcdt = omega(c, T, p);

// constant pressure
// dT/dt = ...
scalar rho = 0.0;
scalar cSum = 0.0;
for (label i = 0; i < nSpecie_; i++)
{
const scalar W = specieThermo_[i].W();
cSum += c[i];
rho += W*c[i];
}
const scalar mw = rho/cSum;
scalar cp = 0.0;
for (label i=0; i<nSpecie_; i++)
{
const scalar cpi = specieThermo_[i].cp(T);
const scalar Xi = c[i]/rho;
cp += Xi*cpi;
}
cp /= mw;

scalar dT = 0.0;
for (label i = 0; i < nSpecie_; i++)
{
const scalar hi = specieThermo_[i].h(T);
dT += hi*dcdt[i];
}
dT /= rho*cp;

// limit the time-derivative, this is more stable for the ODE
// solver when calculating the allowed time step
const scalar dTLimited = min(500.0, mag(dT));
dcdt[nSpecie_] = -dT*dTLimited/(mag(dT) + 1.0e-10);

// dp/dt = ...
dcdt[nSpecie_ + 1] = 0.0;
}


以上の部分の中盤の

scalar cp = 0.0;
for (label i=0; i<nSpecie_; i++)
{
const scalar cpi = specieThermo_[i].cp(T);
const scalar Xi = c[i]/rho;
cp += Xi*cpi;
}
cp /= mw;

で cp の単位が[J/K/mol]なのであれば cp /= mw; ではなく cp *= mw;ではないでしょうか？
Xi=mw*c[i]/rho*cpi であり、mwで割る必要はないと考えました。
（そもそもこの部分がreactingFoam上で使われているのかさえわからないほど勉強不足で申し訳ないのですが…）


以上、長文でしかも関係ない質問までして申し訳ないのですが、
間違っていればご指摘をお願いいたします。

[ohbuchi]

確かにそうですね。

mw=∑(Wi*ci)/∑ci
ですから

cp=∑(ci/∑Wi*ci)*cpi/(∑Wi*ci/∑ci)
では、おかしなことになります。

cp=∑(ci/∑Wi*ci)*cpi * (∑Wi*ci/∑ci)

であれば、Cpをモル比で重み付け平均したと見なせると思います。

ただ、今回のdQの結果と関係があるかといえばなさそうにみえますが。

[Shintaro]

計算結果のmaxCoに対する依存性ですが、以下のサイトで論じられていました。
http://www.openfoam.com/mantisbt/view.php?id=249
このサイトによると、chemistry.H中で
kappa = (runTime.deltaT() + tc)/(runTime.deltaT() + tc + tk);
となっており、結果的にはShがrunTime.deltaT()に影響を受けていたようです。

[ohbuchi]

こんにちは。
試しに、kappa = tc / ( tc + tk )　と修正してmacCo=0.4,0.1の結果を比べたところ
maxCo=0.4 dQmax = 10.685
maxCo=0.1 dQmax = 32.452
となり、まだDT依存が残っています。
以前よりは若干依存性が弱まりましたが。
どうもこれだけではない様な気がします。

[ohbuchi]

kappa≡1として計算したところ、チュートリアルが変な結果になりました。
やはり、Chalmers PaSr modelが曲者の様です。

下記のプレゼンにalternativeReactingFoamの紹介があり、
http://www.opensourcecfd.com/conference2008/media/proceedings/OSCIC-08_GschaiderRehm.pdf

DTの代わりに流れの時間スケールを用いる（定常計算向きの）モデルが示されています。
こちらの方がDT依存性は少ないと思います。

[Shintaro]

なるほど、問題はchemistry.Hだけではないのですね。
わざわざ確かめていただいてありがとうございます。

altemativeReactingFoamのことは知っていまいしたが、
DTではなく流れの時間スケールを用いているとは知りませんでした。
現在、reactingFoamがFluent等のＣＦＤと比べてどの程度使えるのか、
という問題に取り組んでおり、その仕事が一段落つきましたら是非試してみたいと思います。