ファイルk出力
aki-y1985
ユーザーガイドの2.1にあるcavity flowのチュートリアルをpisoFoamを使って勉強しています．
乱流モデルはkEpsilonを使っています．
今，OpenFOAM/OpenFOAM-1.6/src/turbulenceModels/incompressible/RAS/
kEpsilon/kEpsilon.C
(下記に抜粋させて頂きます)
223 // Dissipation equation
224 tmp<fvScalarMatrix> epsEqn
225 (
226 fvm::ddt(epsilon_)
227 + fvm::div(phi_, epsilon_)
228 - fvm::Sp(fvc::div(phi_), epsilon_)
229 - fvm::laplacian(DepsilonEff(), epsilon_)
230 ==
231 C1_*G*epsilon_/k_
232 - fvm::Sp(C2_*epsilon_/k_, epsilon_)
233 );
の中のfvm::ddt(U) や，fvm::div(phi_, epsilon_)，fvm::Sp(fvc::div(phi_),
epsilon_)
といった項を各時刻で A，B，Cといったファイルで出力したいのですが，kEpsilon.Cの
ソースコードをどのように変更すれば良いかわかりません．
試しに，
volScalarField A = fvm::ddt(U);
volScalarField B = fvm::div(phi_, epsilon_);
volScalarField C = fvm::Sp(fvc::div(phi_), epsilon_);
と定義して，
tmp<fvScalarMatrix> epsEqn
(
A
+ B
- C
と変更しただけでは上手くいきませんでした．どのようにすればよろしいのでしょうか．
ご存知の方がいらっしゃいましたら，是非ともご教示のほどをお願い致します．
よろしくお願いいたします．
aki
aki-y1985
正しくは，「ファイルの出力について」です．大変申し訳ありませんでした．
ohb...@amber.plala.or.jp
なかなかResがつきませんね。
やりたい事は、Navier-Stokes式や乱流モデルの項別の履歴をファイルに記録したいこと
でしょうか？
それならば、ライブラリのソースを変更せずに、アプリケーションのソースを変更しましょう。
ライブラリソースに手を加えるのは、サポートされていない境界条件やメッシュ移動など
本当に必要な時だけにしましょう。また、手を加えるとしても派生クラスを作る形にして
元のソースには手を加えないのが基本ルールです。
オブジェクト指向の開発スタイルを勉強してください。
ソルバーの冒頭でインクルードしているcreateFields.Hで記録用のフィールドを定義して、
IOobject：：AUTO_WRITEとしておきます。
次に、ソルバの計算ループの終わりに、各フィールドデータに値を代入する文を記述
すれば、各Timeディレクトリにフィールド名のついたファイルが生成される筈です。
aki-y1985
早速ご返信頂きありがとうございます！
> やりたい事は、Navier-Stokes式や乱流モデルの項別の履歴をファイルに記録したいこと
> でしょうか？
はい．おっしゃる通り，各時刻で圧力や速度をpやUのファイルで出力しているのと同様に，
ナビエ・ストークス式やイプシロンの式の中のfvm::ddt(epsilon_)やC1_*G*epsilon_/k_
などの値を，AやBなどの適当な名前を付けて各時刻でファイルに出力したいと思っております．
> それならば、ライブラリのソースを変更せずに、アプリケーションのソースを変更しましょう。
> ライブラリソースに手を加えるのは、サポートされていない境界条件やメッシュ移動など
> 本当に必要な時だけにしましょう。また、手を加えるとしても派生クラスを作る形にして
> 元のソースには手を加えないのが基本ルールです。
> オブジェクト指向の開発スタイルを勉強してください。
ご経験を積まれた方のご助言は，大変勉強になります．
OpenFOAMを勉強するたびにC++言語の不勉強を痛感しており，早くオブジェクト指向の
考え方について理解したいと思っております．
> ソルバーの冒頭でインクルードしているcreateFields.Hで記録用のフィールドを定義して、
> IOobject：：AUTO_WRITEとしておきます。
> 次に、ソルバの計算ループの終わりに、各フィールドデータに値を代入する文を記述
> すれば、各Timeディレクトリにフィールド名のついたファイルが生成される筈です。
試しにレイノルズ応力テンソルRの項を出力することを考えまして，ご助言を参考に
createfields.Hの最後にRに関する宣言を以下のように追加し，
createfields.H
-----------------------------------------
volSymmTensorField R
(
IOobject
(
"R",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
turbulence->R()
);
-----------------------------------------
また，interFoam.Cを
interFoam.C
-----------------------------------------
turbulence->correct();
R = turbulence->R(); // 追加部分
R.write(); // 追加部分
runTime.write();
-----------------------------------------
上記のように修正したところ，Rを無事に出力できました！
しかし，どうもsystem/controlDictのwriteIntervalで指定した時間刻み
ではなく，毎時間ステップで出力してしまうようで，
0.000119904，0.00025673，0.0003475824・・・のようなフォルダがたくさん
できてしまいます．
また，任意の項，例えばC1_*G*epsilon_/k_を出力しようと思いまして，
上記の例と同様にcreatefields.Hの最後に
createfields.H
-----------------------------------------
volScalarField A
(
IOobject
(
"A",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
C1_*G*epsilon_/k_
);
-----------------------------------------
と追加し，interFoam.Cも
interFoam.C
-----------------------------------------
turbulence->correct();
A = turbulence->A(); // 追加部分
A.write(); // 追加部分
runTime.write();
-----------------------------------------
と変更したのですが，計算させる以前のinterFoamをwmakeする時点で
----------------------------------------------------------------------------------
createFields.H: In function ‘int main(int, char**)’:
createFields.H:119: error: ‘C1_’ was not declared in this scope
createFields.H:119: error: ‘G’ was not declared in this scope
createFields.H:119: error: ‘epsilon_’ was not declared in this scope
createFields.H:119: error: ‘k_’ was not declared in this scope
interFoam.C:96: error: ‘class Foam::incompressible::turbulenceModel’
has no member named ‘A’
----------------------------------------------------------------------------------
というエラーが出てしまいました．
エラーメッセージから想像するに，各定数や変数が定義されていないからなのかと思いまして，
色々といじってみたのですが，私の能力では一向に解決しませんでした．
一体どのように修正すればよろしいのでしょうか，もしよろしければお教え頂けますと大変有難いです．
どうぞよろしくお願い申し上げます．
aki
ohb...@amber.plala.or.jp
> createFields.H: In function ‘int main(int, char**)’:
> createFields.H:119: error: ‘C1_’ was not declared in this scope
> createFields.H:119: error: ‘G’ was not declared in this scope
> createFields.H:119: error: ‘epsilon_’ was not declared in this scope
> createFields.H:119: error: ‘k_’ was not declared in this scope
> interFoam.C:96: error: ‘class Foam::incompressible::turbulenceModel’
> has no member named ‘A’
interFoam.Cの中で、turbulence->Aという存在しないメソッドを呼び出していることが
エラーとなっています。
turbulence->Rは実際に存在するメソッドだったため、うまく行ったのでしょう。
IMANO Masashi
akiさん
今野です。
> しかし，どうもsystem/controlDictのwriteIntervalで指定した時間刻み
> ではなく，毎時間ステップで出力してしまうようで，
> 0.000119904，0.00025673，0.0003475824・・・のようなフォルダがたくさん
> できてしまいます．
コード追加された R.write() は、その時点でのRを出力しますので、時間ステッ
プのループ内に置いた場合には、各時点で出力されてしまいます。
akiさんはRの宣言で、元々 IOobject::AUTO_WRITE と指定しており、
runTime.write() により、writeIntervalで指定した時間刻みで出力されるの
で、R.write() は不要ということになります。
さて、本題の
> の中のfvm::ddt(U) や，fvm::div(phi_, epsilon_)，fvm::Sp(fvc::div(phi_),
> epsilon_)
> といった項を各時刻で A，B，Cといったファイルで出力したいのですが，kEpsilon.Cの
> ソースコードをどのように変更すれば良いかわかりません．
ですが、大渕さんの言われる
> それならば、ライブラリのソースを変更せずに、アプリケーションのソースを変更しましょう。
> ライブラリソースに手を加えるのは、サポートされていない境界条件やメッシュ移動など
> 本当に必要な時だけにしましょう。また、手を加えるとしても派生クラスを作る形にして
> 元のソースには手を加えないのが基本ルールです。
はごもっともで、本来そうすべきですが、どうやらpisoFoamとか単体のソルバ
だけで、各項をモニターしたいわけでもないようなので、今回はとりあえずデ
バッグ用ということで、kEpsilonのソースコードを直接いじる方法を説明しま
す。
*出力したい各項(ここではepsilonA,epsilonB,epsilonCとする)をnut同様に
kEpsilonのConstructorsでAUTO_WRITEで宣言します。初期値は各項の計算方法
により与えます(次元と初期値を陽に指定する方法でも良いです)。
例:
epsilonA_
(
IOobject
(
"epsilonA",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
fvc::ddt(epsilon_)
),
*epsilonを解いた後で、各項を毎度計算します。
例:
epsilonA_=fvc::ddt(epsilon_);
ここで、epsilonの輸送方程式では fvm::ddt(epsilon_) 等と
fvm::が使われてますが、fvm::は陰解法用にMatrixを作成するので、
volScalarFieldのような陽に決まっている場を求めるには、fvc::に変える必
要があります。
*epsilonA_等の宣言をkEpsilon.Hで行っておく必要があります。
例:
volScalarField epsilonA_;
以下、kEpsilon.[C,H]へのパッチです。
行頭が+の行を+を除いて加えてください。
今野 雅 // IMANO Masashi
--- kEpsilon.C.orig 2010-03-06 14:32:52.000000000 +0900
+++ kEpsilon.C 2010-09-15 10:38:47.000000000 +0900
@@ -126,6 +126,42 @@
IOobject::AUTO_WRITE
),
autoCreateNut("nut", mesh_)
+ ),
+ epsilonA_
+ (
+ IOobject
+ (
+ "epsilonA",
+ runTime_.timeName(),
+ mesh_,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ fvc::ddt(epsilon_)
+ ),
+ epsilonB_
+ (
+ IOobject
+ (
+ "epsilonB",
+ runTime_.timeName(),
+ mesh_,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ fvc::div(phi_, epsilon_)
+ ),
+ epsilonC_
+ (
+ IOobject
+ (
+ "epsilonC",
+ runTime_.timeName(),
+ mesh_,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE
+ ),
+ fvc::Sp(fvc::div(phi_), epsilon_)
)
{
nut_ = Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
@@ -239,6 +275,9 @@
solve(epsEqn);
bound(epsilon_, epsilon0_);
+ epsilonA_=fvc::ddt(epsilon_);
+ epsilonB_=fvc::div(phi_, epsilon_);
+ epsilonC_=fvc::Sp(fvc::div(phi_), epsilon_);
// Turbulent kinetic energy equation
tmp<fvScalarMatrix> kEqn
--- kEpsilon.H.orig 2010-03-06 14:32:52.000000000 +0900
+++ kEpsilon.H 2010-09-15 10:36:53.000000000 +0900
@@ -81,6 +81,9 @@
volScalarField k_;
volScalarField epsilon_;
volScalarField nut_;
+ volScalarField epsilonA_;
+ volScalarField epsilonB_;
+ volScalarField epsilonC_;
public:
谷藤亜希
早速ご助言を頂きましてありがとうございました！
Rができてしまいましたので，単純にそれをAに変えれば良いという安直な
発想でした・・・．今後は，もう少しコードを読めるように地力をつけようと思います．
今野さん，
コードの詳細にわたるアドバイスを頂きまして，本当にありがとうございました！
お教え頂いた通りにコードを変更してみたところ，無事に全てのファイルが
writeIntervalで指定した時間刻みで出力されました！
まだコード構造を完全には理解できてはおりませんが，お教え頂いたことを
ヒントに理解していこうと思います．お忙しい中快く答えて頂き，大変恐縮です．
回りにはこれほどOpenFOAMやC++言語に精通する方がおらず，ただただ尊敬します．
aki
谷藤亜希
先日はソースコードの変更方法をお教えくださり，ありがとうございました．
その後，先日お教え頂いた方法を参考にkEpsilonのソースコードを直接いじってみていますが，
一つだけまた疑問点が出てきましたので，質問させて頂きます．
今，下記の行をkEpsilon.Cに加えています．
epsilonA_=fvc::ddt(epsilon_);
epsilonB_=fvc::div(phi_, epsilon_);
epsilonC_=fvc::Sp(fvc::div(phi_), epsilon_);
epsilonD_=- fvc::laplacian(DepsilonEff(), epsilon_);
epsilonE_=C1_*G*epsilon_/k_;
epsilonF_=- fvc::Sp(C2_*epsilon_/k_, epsilon_);
epsilonsum_=epsilonA_+epsilonB_+epsilonC_+epsilonD_-epsilonE_-epsilonF_ ;
このとき，epsEqnの式が
// Dissipation equation
tmp<fvScalarMatrix> epsEqn
(
fvm::ddt(epsilon_)
+ fvm::div(phi_, epsilon_)
- fvm::Sp(fvc::div(phi_), epsilon_)
==
C1_*G*epsilon_/k_
- fvm::Sp(C2_*epsilon_/k_, epsilon_)
);
ですので，epsilonsum_は常に0になるべきだと思うのですが，実際に出力させてみると，
その値が0にはなりません．特に，境界付近が0とはかけ離れた値になります．
最初は空間分割数が少ないせいかと思いましたが，空間刻み幅を細かくしても同じ傾向でした．
これは，一体何が原因と考えられますでしょうか？もしよろしければご助言を頂けますと有難いです．
どうぞよろしくお願い申し上げます．
aki
IMANO Masashi
今野です。
返事が遅くなりました。
> 今野さん，
>
> 先日はソースコードの変更方法をお教えくださり，ありがとうございました．
> その後，先日お教え頂いた方法を参考にkEpsilonのソースコードを直接いじってみていますが，
> 一つだけまた疑問点が出てきましたので，質問させて頂きます．
ただ、ディスカッションボードでこのように回答者の指名のような形で書くと、
他の人からの回答が得られにくくなるので、できれば指名していないように書
いたほうが良いと思います。
> ですので，epsilonsum_は常に0になるべきだと思うのですが，実際に出力させてみると，
> その値が0にはなりません．特に，境界付近が0とはかけ離れた値になります．
壁境界付近の epsilon 輸送方程式 epsEqn の残差が0から大きく離れるのは、
壁に隣接する格子では壁関数による境界条件 (epsilonWallFunctions) を用いて、
以下のように epsilon を別途明示的に与えているからだと思います。
src/turbulenceModels/incompressible/RAS/derivedFvPatchFields/wallFunctions/epsilonWallFunctions/epsilonWallFunction/epsilonWallFunctionFvPatchScalarField.C:
void epsilonWallFunctionFvPatchScalarField::updateCoeffs()
:
// Set epsilon and G
forAll(nutw, faceI)
{
label faceCellI = patch().faceCells()[faceI];
scalar yPlus = Cmu25*y[faceI]*sqrt(k[faceCellI])/nuw[faceI];
epsilon[faceCellI] = Cmu75*pow(k[faceCellI], 1.5)/(kappa_*y[faceI]);
また、壁隣接格子以外で0にならないのは以下が原因だと思います。
*epsilon の輸送方程式を離散化した線型一次方程式は、CG法等の反復法で解か
れるが、残差が system/fvSolution の tolerance (絶対誤差の許容量)や
relTol (初期残差からの相対値の許容量) 以下になるように不完全に解かれ
るだけなので、そもそも完全には残差が0にはならない。
epsilon
{
solver PBiCG;
preconditioner DILU;
tolerance 1e-05;
relTol 0.1;
}
*epsilonの輸送方程式を反復法で解いた後に、安定性のために、epsilonを
下限値である epsilon0 以上にする bouding が行なわれるので、それが行なわれた
格子や隣接格子では、残差が増える。
solve(epsEqn);
bound(epsilon_, epsilon0_);