流れが層流になるのですが・・・

[mei]

皆様，こんにちは

以前，速度が対数則にのらないと質問した者です．確認したところ流れが層流化していました．

現在，チャネル流れにおいて速度の減衰をfvOptionsで速度指定することで防いでおります．paraviewでの速度を見ても減衰はしていませんでした．しかし，なぜか層流化しました・・・・

計算はLES(pisofoam)でモデルはSmagorinskyモデル，van Driest関数を使用しています．また，計算領域は(x,y,z)=(0.064,0.010,0.016)[m]で動粘性係数はnu=1.5*10^-5[m^2/s]です．

レイノルズ数はReb=5650(Re=Reb/2=2825)として

Reb=(2*Ub*D)/nu

からUbを算出してそれを指定しています．Dは平行平板間距離です．

境界条件は入口出口および左右を周期境界条件cyclic，上下を壁面wallに指定しています．

初期乱流場は速度場は他の研究でのチャネル流れの速度場を使用し，圧力場は初期零圧力です．圧力の境界条件は壁面零圧力勾配です．

計算式が違うのかといろいろしていましたが，根本的なところができていませんでした．

fvOptionsの使い方が間違っているのでしょうか？

見てみるとほぼ完全に層流化し，圧力場は零圧力になっていました．

流れが層流化する原因として何が考えられますでしょうか？

また，別の質問ですが領域全体に圧力勾配をかけることはできますでしょうか？ブラジウスの式より圧力損失を見積ることができるので，そう思ったのですが・・・

できればそちらで流れを維持したいと考えています・・・

よろしければ，どなたかご教授お願いいたします．

[Youhei Takagi]

高木と申します。

ソルバーは何をお使いでしょうか？周期境界条件のチャンネル流れ

でしたら、channelFoamというソルバーで解くことができます。

詳しくは過去にも同様な質問がありましたので、検索して探して

みてください。

以上ご参考までに。

[Masashi Imano]

今野です。

もしfvOptionsを使っているとすると，OpenFOAMのバージョンはv.2.2.0以降だと
思いますが，その場合にはchannelFoamはfvOptionsの実装により廃止されおりて，
channel流はpimpleFoam等の通常のソルバーで解くことができると思います。

以下はv.2.2.0以降をご使用の場合を仮定します。

> fvOptionsの使い方が間違っているのでしょうか？

* ご懸案のsystem/fvOptionsを貼ってください。

* fvOptionsを読み込んでいることを確認するために，ソルバーのログの初期部分を貼ってください。

tutorials/incompressible/pimpleFoam/channel395での例(乱流モデルはSmagorinskyに変更):
ーーー
Selecting incompressible transport model Newtonian
Selecting turbulence model type LESModel
Selecting LES turbulence model Smagorinsky
Selecting LES delta type vanDriest
Selecting LES delta type cubeRootVol
SmagorinskyCoeffs
{
ce 1.048;
ck 0.094;
}

Creating finite volume options
Creating fintite volume options from fvOptions

Selecting finite volume options model type pressureGradientExplicitSource
Source: momentumSource
- applying source for all time
- selecting all cells
- selected 60000 cell(s) with volume 16

Initial pressure gradient = 0
ーーー

> 流れが層流化する原因として何が考えられますでしょうか？

もし，Smagorinskyモデルのデフォルトの乱流モデルパラメタをお使いの場合には，
Smagorinsky係数のCsが約0.167となり，チャンネル流での推奨値である0.1より大部
大きくなっており，この場合チャンネル流の性状が正しく再現できず，場合によっては
層流に近いプロファイルになると思います。

Smagorinsky model details
http://www.cfd-online.com/Forums/openfoam-solving/65838-smagorinsky-model-details-2.html

また，そもそもOpenFOAMのようなコロケーショングリッドでは，スタッガードグリッドに比べて
エネルギー保存性が悪いために，チャンネル流中央での平均速度を過大に評価するとの指摘があります。

大岡ほか，LESにおけるコロケーショングリッドのエネルギー非保存性の検討 : channel計算によるスタッガードグリッドとコロケーショングリッドの比較
http://repository.dl.itc.u-tokyo.ac.jp/dspace/bitstream/2261/53036/1/sk049001004.pdf

> また，別の質問ですが領域全体に圧力勾配をかけることはできますでしょうか？ブラジウスの式より圧力損失を見積ることができるので，そう思ったのですが・・・

領域全体に一定の(平均)圧力勾配をかけるには，運動量の輸送方程式に圧力の空間勾配分の
ソース項を足せば良いので，fvOptionsのvectorSemiImplicitSourceでUに圧力勾配分を
一定ソース項であるSuに指定すれば良いのではないでしょうか?

system/fvOptions:
---
momentumSource
{
type vectorSemiImplicitSource;
active true; //on/off switch
selectionMode all; //cellSet // points //cellZone

vectorSemiImplicitSourceCoeffs
{
volumeMode specific; // absolute; //
injectionRateSuSp
{
U ((1 0 0) 0);
}
}
}
---

以上です。

2014年1月27日 2:42 Youhei Takagi <yotak...@gmail.com>:

> --
> このメールは Google グループのグループ「OpenFOAM」の登録者に送られています。
> このグループから退会し、メールの受信を停止するには、openfoam+u...@googlegroups.com にメールを送信します。
> このグループに投稿するには、open...@googlegroups.com にメールを送信してください。
> http://groups.google.com/group/openfoam からこのグループにアクセスしてください。
> その他のオプションについては、https://groups.google.com/groups/opt_out にアクセスしてください。



--
IMANO Masashi, Ph.D.

[mei]

高木様
今野様
ご返信ありがとうございます．

確かに，OpenFOAMはver-2.2.xでchannelFoamはありません．
また，fvOptionsとログについて添付いたしました．

また，Smagorinskyはデフォルトのまま使用しています．
とりあえずこちらを変えてみます．

よろしくおねがいします．

[Masashi Imano]

今野です。

fvOptionsはこれで良いと思います。
また，ソルバーのログを見ても正しく読み込んでます。
しかし，乱流モデルのパラメータがデフォルトのままなので，
Csが0.1なるようにパラメータを変更してみてください。

上記の変更後も非物理的な結果が得られる場合には，
離散化スキームやメッシュ分割が疑われますので，
system/fvSchemesやメッシュ図または分割数を示すと共に，
postChannelによるプロファイル図やyPlusLESのログも貼ってください。

以上です。


2014年1月27日 13:26 mei <kanade...@kke.biglobe.ne.jp>:

[mei]

今野様
ありがとうございます．

係数を変更ということでいままさにやろうとしていたのですが，Smagorinsky.Cに記述されている乱流エネルギーk，渦粘性係数nusgsおよび参考書SmagorinskySGS渦粘性係数から

Smagorinsky定数Cs=Ck*sqrt(2*Ck/Ce)

となるかと思ったのですが，この場合Ck，Ceはどのような基準で決められているのでしょうか？
よろしければご教授お願いいたします．


Smagorinsky.Cより
nusgs=Ck*sqrt(k)*delta
k=(2*Ck/Ce)*delta^2*|D|^2

よって
nusgs=Ck*sqrt(2*Ck/Ce)*delta^2*|D|^2

参考書では
nusgs=(Cs*delta)^2*|D|^2

したがって
Cs=Ck*sqrt(2*Ck/Ce)

[mei]

皆様，こんにちは

サイトの方からSmagorinsky定数が
Cs=Ck*sqrt(2*Ck/Ce)
ではなく
Cs=Ck*sqrt(Ck/Ce)
で表される理由についてはなんとなく理解出来ましたが，Ck，Ceについてはよくわかりませんでした．
教科書にはSmagorinsky定数はこのモデルに対して与える唯一の無次元定数とありますが，乱流エネルギーkに対してはCkを用いておりよくわかりません．
サイトの方のSunxingさんが最後にCs=0.1にするために定数を以下のように設定していますがこれはどのような理由からかをご教授していただけないでしょうか？

SmagorinskyCoeffs
{
ce 1.05;
ck 0.0472;
}

[mei]

皆様，こんにちは

少し進展があったので，ご報告させて頂きます．
教授に尋ねたところ，チャネルLESの対流項に対してはQUICK法を用いてはいけないという事でした．
数値粘性が効き過ぎるのが理由で層流化したんだろうとの見方です．
Smagorinsky定数の推奨値が小さいのも数値粘性を抑えるためという事で，現在対流項を中心差分で計算しています．
計算が終わった時にまた報告させて頂きますので，まだ層流になるようでしたらお力を貸していただければ幸いです．

[ONO Hiroki]

小野です。

教科書などでよくあらわれる式はnusgs=(Cs*delta)^2*|S|^2　…(1)ですが、
本来Smagorinskyモデル(非圧縮性の場合)は

"""
nusgsはsqrt(k)とdeltaに比例すると仮定して以下の式から求める(Ckは比例定数)。
nusgs=Ck*sqrt(k)*delta …(2)
なお、kは局所平衡を仮定して以下のように定める
k=(Ck/Ce)*delta^2*|S|^2 …(3)
"""

という考え方からできているモデルです。まとめられる定数はまとめてしまった方がコーディングも楽だしメモリ節約にもなることもあり、
通常は(2)に(3)を代入、整理した(1)を掲載することが多いですね。


Ceは、k(SGS)の輸送方程式を考えたとき、その散逸項に現れるモデル係数で、あまり流れ場に依存せずほぼ1前後であるといわれています。

よって、「Smagorinsky定数を調節する際には、Ceを変えずにCkで調節する」というのがある種の定番になっているようですが、
計算に影響するのは最終的に求められるnusgsの値なので、Ck*sqrt(Ck/Ce)の値が変わらなければどのようにCk,Ceをいじっても流れ場の結果には影響しないはずです(もちろんゼロ割などはアウトですが)。
ただし、kの値を出力する場合には影響があるので、やはりCeは1または1.05に固定してCkで調節するのが良いでしょう。





2014年1月28日火曜日 14時41分16秒 UTC+9 mei:

[mei]

小野様

ご返信ありがとうございます．

 

なるほど，とりあえずSunxingさんにならいCeをデフォルトでCkを小さくして計算していましたが，その理由がよく理解できました．

分かりやすいご説明ありがとうございました．

[mei]

皆様、こんにちは
いまだ計算中ですが層流にはなっていないようです．しかし，やはり対数則にのらないです．

メッシュが悪いのか，離散化法が悪いのかいろいろ変更してみます．

そこで，今野様のおっしゃられた圧力勾配を全体にかけるというのをやってみようと思ったのですが，

momentumSource
{
    type vectorSemiImplicitSource;
    active          true;            //on/off switch
    selectionMode   all;       //cellSet // points //cellZone

    vectorSemiImplicitSourceCoeffs
    {
        volumeMode      specific; // absolute; //
        injectionRateSuSp
        {
            U           ((1 0 0) 0);
        }
    }
}

ここでU((1 0 0) 0)は何を指定しているのでしょうか？(1 0 0)は見た感じ速度成分と思ったのですが，その次の0がわかりません．
CFDonlineにおいてenergySource項においてh(10 0)に対し，10はexplicit component(明白な要素)，0がinplicit component(暗に示された要素)とあるのですが，？？？状態です．(静止状態から1m/sにいたるまでの圧力？)

どなたか教えていただけませんか？

[すもも]

こんにちわ。

こちらで
http://www.openfoam.org/docs/cpp/

"SemiImplicitSource" で検索して出てくる一番上の "SemiImplicitSource" を表示して、
その Detailed Description の説明を参照してください。

数値計算の用語では、"explicit" は陽的、"implicit" は陰的です。
詳細は省きますが、ソース項の扱いには陰的な方法と陽的な方法があり、
陰的な方法では非線形ソース項を線形化して扱うために、

S = S0 + dS*X

のような表し方をします。X が方程式の変数です。

U           ((1 0 0) 0);

は、たぶん、最初が S0、次が dS に当たるのかなーと思います。
この辺のお話は、パタンカー「コンピュータによる 熱移動と流れの数値解析」に載っています。


2014年1月31日金曜日 0時45分25秒 UTC+9 mei:

[mei]

すもも様

 

ご返信ありがとうございます．勉強不足でお恥ずかしい限りです・・・

教えていただいた本が図書館にありましたので，読んでみます．

 

ありがとうございました．

 

[mei]

みなさま，こんにちは

まだ本を理解してる途中ですが，プログラム上で疑問が生じたので質問させてください．

fvOptionsというのはプログラム(計算)においてパラメータを固定するために存在するということでよろしいでしょうか？
でしたら，例えば平均圧力勾配を流れに加えることで流れを維持しますがsimpleFoamでは
    tmp<fvVectorMatrix> UEqn
    (
        fvm::div(phi, U)
      + turbulence->divDevReff(U)
      == fvOptions(U)
    );
とあり，fvOptionsは

    pressureGradientExplicitSourceCoeffs
    {
        fieldNames  (U);
        Ubar        ( 10 0 0 );
    }

このように記述し，この場合、バルク速度を一定値に設定しています．
しかし，平均圧力勾配を加える(別のパラメータを固定する)場合，fvOptions(U)はそのままの記述で大丈夫ですか？
つまり，大本のソルバーに手を加えず，fvoptionsの中身だけを変えれば流れに影響するパラメータなどは変化するのでしょうか？

よろしくおねがいします．

[すもも]

こんにちわ。

fvOptions は、MRF や多孔質体抵抗まで含めたいろんなことを、
プログラムの修正なしで (設定を変えるだけで) 行うためのもののようです。

なので、プログラムを修正する必要はないと思います。


2014年2月3日月曜日 1時51分11秒 UTC+9 mei:

[mei]

すもも様

 

ご返信ありがとうございます．

プログラムは修正する必要なかったんですね．わかりました．

ありがとうございました．

[Youhei Takagi]

mei様

高木です。

ソースコード、pressureGradientExplicitSource.C
を見ると、

    // Integrate flow variables over cell set
    scalar magUbarAve = 0.0;
    scalar rAUave = 0.0;
    const scalarField& cv = mesh_.V();
    forAll(cells_, i)
    {
        label cellI = cells_[i];

        scalar volCell = cv[cellI];
        magUbarAve += (flowDir_ & U[cellI])*volCell;
        rAUave += rAU[cellI]*volCell;
    }

    // Collect across all processors
    reduce(magUbarAve, sumOp<scalar>());
    reduce(rAUave, sumOp<scalar>());

    // Volume averages
    magUbarAve /= V_;
    rAUave /= V_;

    // Calculate the pressure gradient increment needed to adjust the average
    // flow-rate to the desired value
    scalar gradPplus = (mag(Ubar_) - magUbarAve)/rAUave;

    // Apply correction to velocity field
    forAll(cells_, i)
    {
        label cellI = cells_[i];
        U[cellI] += flowDir_*rAU[cellI]*gradPplus;
    }

    // Update pressure gradient
    gradP_.value() += gradPplus;


という計算をしており、流量一定条件を満たすように平均圧力勾配を加えていることが
わかります。

以上ご参考までに。

[mei]

高木様

 

ソースコードありがとうございます．思ったより読めてほっとしてます．

（flowDirは流れ方向のことでしょうか・・）．

 

層流にはならなくなり，あとは精度の問題のみとなりました(LES，RANSともに・・・)．

どうもutauが小さく算出されます．これによりU+は大きく，y+は小さく出ます・・・・

あとは離散化やメッシュなどいろいろ試していこうかと考えています．

 

みなさまお忙しい中，多くの質問に答えていただき，本当にありがとうございました．