SmagorinskyCoeffsの物理的意味

[kazui]

皆さん　こんにちは。 kazuiです。

今、LESのSmagorinskyモデルを使った解析を行おうかと思いますが、

SmagorinskyCoeffsにある、ceとckの物理的な意味を教えて頂けないでしょうか？

下記のサイトから抜粋

http://lists.cfd-online.com/pipermail/openfoam/2011-June/003420.html

SmagorinskyCoeffs
{
ce 1.05;
ck 0.07;
}

NS式には乱流粘性としてフィードバックされるかと思いますが、

乱流粘性νtとce、ckの関係も分かる資料等があれば助かります。

よろしくお願いします。

[nakagawa]

中川です。

LES実装に関心があり，調査を始めたところです。

中途半端な情報ですが，ご参考までに投稿します。

Smagorinsky.Hに説明があります。

http://foam.sourceforge.net/docs/cpp/a08489_source.html

00024 Class
00025     Foam::incompressible::LESModels::Smagorinsky
00026 
00027 Description
00028     The Isochoric Smagorinsky Model for incompressible flows.
00029 
00030     Algebraic eddy viscosity SGS model founded on the assumption that
00031     local equilibrium prevails.
00032     Thus,
00033     \verbatim
00034         B = 2/3*k*I - 2*nuSgs*dev(D)
00035         Beff = 2/3*k*I - 2*nuEff*dev(D)
00036 
00037     where
00038 
00039         D = symm(grad(U));
00040         k = (2*ck/ce)*delta^2*||D||^2
00041         nuSgs = ck*sqrt(k)*delta
00042         nuEff = nuSgs + nu
00043     \endverbatim

下記ディスカッションも参考になるかもしれません。

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

NS式へはsgsModel->divDevBeff(U)項を通じて反映されますが，このマトリクス内でnuEffが使われています。

http://foam.sourceforge.net/docs/cpp/a08465_source.html#l00045

以上です。

より詳しい方にフォローして頂けると有難いです。

2013年1月17日木曜日 10時16分20秒 UTC+9 kazui:

[nakagawa]

中川です。先の投稿がテキストメールでは読みにくかったので再送します。

皆さん　こんにちは。 kazuiです。

[Youhei Takagi]

kazui様

横から失礼します、高木と申します。

LESに詳しくないので、CFD Onlineに載っていた論文、

1. P. P. Sullivan, J. C. McWilliams and C.-H. Moeng, "A subgrid-scale
model for large-eddy simulation of planetary boundary-layer flows",
Boundary-Layer Meteorology, 71, pp.247-276 (1994).

2. C. Fureby, G. Tabor, H. G. Weller and A. D. Gosman, "A comparative
study of subgrid scale models in homogeneous isotropic turbulence",
Phys. Fluids, 9, 1416 (1997).

を少し読んでみました。論文2の方に、kinetic SGS energy (k)の輸送方程式
が以下のように記述されています(論文中のEq.(5))。

d(k)/dt + div(k・u) = -B・D+div(nu_k*grad(k))-e (左辺第1項は時間偏微分です)

ここで、uは速度、BはSGS応力テンソル、Dはひずみ速度テンソル、nu_kは
SGS渦粘性係数、eはSGSエネルギー散逸率です。nu_kとeを局所等方性を
仮定してGS成分で表すと、

nu_k = c_k*k^(1/2)*delta, e = c_e*k^3/2*delta^(-1)

ここで、c_kとc_eが入力値として必要な値であり、deltaはフィルター幅です。
つまり、c_kとc_eはSGS運動エネルギーの輸送方程式において渦粘性と散逸率を
モデル化するのに必要な係数、と考えてよいかと思います。

c_kとc_eをどのような値にするか?ですが、論文1だと、

c_k = 0.10, c_e = 0.93

論文2ですと吉澤先生の論文から引用して、

c_k = 0.05, c_e = 1.00

としています。最適な値は流れ場によるので、何とも言えません。そもそも
局所等方性が成り立たない、という議論もありますが。

以上ご参考までに。もしかしたら私の解釈が間違っているかもしれませんので、
そのときはご容赦ください。

[kazui]

中川様　高木様

kazuiです。　お返事が遅くなり申し訳ありません。

アドバイスありがとうございます。

まずは中川様のアドバイスを元に、

DOXYGENなどで数式から理解をしてみようと試みました。

圧縮性と非圧縮性とで違いがあることも分かりました。

（圧縮性は少し複雑な感じでしたが。）

高木様からの貴重な情報ありがとうございます。

やはり実験定数みたいなものでしたか。

論文も入手してみようと思います。

背景に何があるのかで、流れ場に応じた係数の決め方などを

模索したいと思います。

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

[kazui]

中川様

kazuiです。　こんばんは。

今、圧縮性流体でのLESについて調べており、DOXYGENで追っていくと以下のソースコードに辿り着きました。

http://foam.sourceforge.net/docs/cpp/a08488_source.html

00028     The choric Smagorinsky Model for compressible flows.

00029 
00030     Algebraic eddy viscosity SGS model founded on the assumption that
00031     local equilibrium prevails.
00032     Thus,
00033     \verbatim
00034 

00035         B = 2/3*k*I - 2*nuSgs*dev(D)

00036 
00037     where
00038 
00039         D = symm(grad(U));

00040         k from rho*D:B + ce*rho*k^3/2/delta = 0
00041         muSgs = ck*rho*sqrt(k)*delta
00042     \endverbatim

ここからは数学の問題ですが、係数CkはmuSgsには何乗に比例しておりますでしょうか？

41行目を見ると、１乗に比例するようですが、この式の中にあるkを追っていきますと、

Bが関与しており、Bの中にまたnuSgsがあります。

これは（nuSgs＝muSgs/rho)かと思いますが、そうだとすると、muSgs（nuSgs）とCkの関係が分からなくなってしまいました。

知りたいのは、調整するためのパラメータCkが、乱流粘性に対し、どのようなインパクトがあるか？（何乗に比例するか？）です。

一度closeした話題ですが、フォロー頂けると幸いです。

2013年1月17日木曜日 11時32分23秒 UTC+9 nakagawa:

[nakagawa]

kazuiさん

中川です。

> 41行目を見ると、１乗に比例するようですが、この式の中にあるkを追っていきますと、

> Bが関与しており、Bの中にまたnuSgsがあります。

　この部分について教えてください。

　kにBが関与しているということですが、どこに、どのように書いてありますか？

2013年1月22日火曜日 19時08分18秒 UTC+9 kazui:

[kazui]

中川様

kazuiです。ご返信ありがとうございます。

ご質問の件、回答します。

35行目です。

00035         B = 2/3*k*I - 2*nuSgs*dev(D)
00036 
00037     where
00038 
00039         D = symm(grad(U));
00040         k from rho*D:B + ce*rho*k^3/2/delta = 0
00041         muSgs = ck*rho*sqrt(k)*delta

コードではなく、説明文ですが、

僕なりの解釈では以下の通りです。

①35行目のBについての説明があり、その下37行目以降に、式におけるD項・k項の出処について記載してある。

②B項にnuSgsがあり、これはSGS動粘性。

③41行目にmuSgs、SGS分子粘性についての記載があり、ここにもkの1/2乗が入っている。

　41行目の両辺をrho（密度）で割れば、nuSgsと等価になる。

　35行目に出てくるnuSgsとの関連性はどうなるのか、結局ckとnuSgsは1乗で比例していないのでは？

混乱している要因として、35行目のIという項目が、何の行列かよく分かりません。

40行目でBとDの内積をとっていることから、Bは行列（テンソル）と認識しております。

35行目と39行目は行列の演算をしており、

40行目はスカラーの演算になっていることから、ckとnuSgsの関係が数学的にどうなっているのか？

分からなくなってしまった次第です。

解釈に誤りがあるのでしょうか？

数学とLESをもっと勉強すればよいのですが、ご存知でしたら…と思い、

投稿させて頂きました。

よろしくお願いします。

2013年1月23日水曜日 0時15分41秒 UTC+9 nakagawa:

[nakagawa]

小野先生

　詳細，かつ，わかりやすいご解説，ありがとうございました。

kazui様

　もし，ckの効き具合が気になるようでしたら， compressibleのSmagorinsky.Cの46行目からのソースを確認して，簡易的にckとmuSgsの関係をグラフ化してはいかがでしょうか。

　流体中の任意の1点については，速度勾配を仮定すれば，手計算でkが出せ，muSgsも計算できます。表計算ソフトで計算すれば，ckだけが変化したときのmuSgsの変化を見ることもできると思います。速度場自体は変わらないので，あくまで簡易的ですが。

中川

2013年1月24日木曜日 1時15分46秒 UTC+9 ONO Hiroki:

小野です。

まずckですが、Smagorinskyモデルなどの渦粘性型SGSモデルでは、
nuSgsは「SGS速度スケール(通常はkSGS^0.5)とSGS長さスケール(通常はΔ)の積」に"比例"すると仮定します。
そのときの比例定数がckになります。

それでは、kSGS(以降、OFのコードに合わせ単にkと表記)をどのように求めるのか、が問題となるわけですが、Smagorinskyモデルでは局所平衡、
すなわち、kの生産(-rho*D:B)と散逸(-ce*rho*k^1.5/Δ)が等しいと仮定してkSGSを求めます。

(-rho*D:B)+(-ce*rho*k^1.5/Δ)=0
非圧縮(rho=一定)では、
D:B+ce*k^1.5/Δ=0


その際、非圧縮仮定(divU=0)の下では
D:Bは、
-2*nuSgs*D:dev(D)　となります。

なぜなら、Iは単位テンソル(1 0 0 0 1 0 0 0 1)ですので、
D:I = div(U) = 0　となるからです。

すなわち局所平衡の式は
-2*nuSgs*D:dev(D)+ce*k^1.5/Δ=0
となり、変形すると、
k^1.5 = 2*nuSgs*D:dev(D)*Δ /ce
nuSgs = ck*k^0.5*Δ
なので、
k^1.5 = 2*ck*k^0.5*Δ*D:dev(D)*Δ /ce
よって、
k = (2*ck /ce )  *Δ^2 * D:dev(D)
からkを求めることが出来ます。

つまり、Smagorinskyモデルでは
nuSgs = ck * [(2*ck /ce )  *Δ^2 * D:dev(D)]^0.5 * Δ
整理して、
nuSgs = (2*ck^3 /ce )^0.5 * Δ^2 * [D:dev(D)]^0.5
となります。

教科書なのでよく見かける式では、
(2*ck^3 /ce )^0.5の部分を単にCまたは(Cs)^2として、
このCまたはCsについて、この流れではこの値、という数字がのっています。
(余談ですが、CまたはCsは乱流統計理論を用いて値が定まります。ところが、この値では全ての流れ場でうまくいくわけではなく、
流れ場に応じて変化させる必要が出てきます。すごく乱暴に言うと、このCを勝手に決めてくれるのがDynamic LESモデルです。)


ところで、圧縮性の場合はD:I = 0とは限らないので、局所平衡の式は
2/3*k*D:I - 2*nuSgs*D:dev(D) + ce*k^1.5/Δ = 0
すなわち、
2/3*k*D:I - 2*ck*k^0.5*Δ*D:dev(D) + ce*k^1.5/Δ = 0
となります。
両辺をk^0.5で割ると(k^0.5)の2次方程式になるので、解の公式から(k^0.5)、
つまりそれを二乗してkが求まります。


なので、
>>結局ckとnuSgsは1乗で比例していないのでは？
については、
「nuSgsは、(k^0.5 * Δ)に対してはckの1乗を比例定数として比例している」と言えると思います。

ちなみに、テンソルに関しては、OF付属のプログラマーズガイドに目を通してみることをおすすめします。
オープンCAE学会では日本語訳も公開しております。
http://www.opencae.jp/wiki/%E3%82%BD%E3%83%95%E3%83%88%E3%82%A6%E3%82%A7%E3%82%A2%E3%83%9E%E3%83%8B%E3%83%A5%E3%82%A2%E3%83%AB%E7%BF%BB%E8%A8%B3



2013年1月23日水曜日 0時40分10秒 UTC+9 kazui:

[ONO Hiroki]

小野です。
一点誤りがありましたので修正して再投稿します。
修正個所は係数をCsとしてまとめる部分です。
-------------------------------------------

まずckですが、

nuSgs = (ck^3 /ce )^0.5 * Δ^2 * [2*D:dev(D)]^0.5
となります。

教科書なのでよく見かける式では、
(ck^3 /ce )^0.5の部分を単にCまたは(Cs)^2として、

このCまたはCsについて、この流れではこの値、という数字がのっています。
(余談ですが、CまたはCsは乱流統計理論を用いて値が定まります。ところが、この値では全ての流れ場でうまくいくわけではなく、
流れ場に応じて変化させる必要が出てきます。すごく乱暴に言うと、このCを勝手に決めてくれるのがDynamic LESモデルです。)


ところで、圧縮性の場合はD:I = 0とは限らないので、局所平衡の式は
2/3*k*D:I - 2*nuSgs*D:dev(D) + ce*k^1.5/Δ = 0
すなわち、
2/3*k*D:I - 2*ck*k^0.5*Δ*D:dev(D) + ce*k^1.5/Δ = 0
となります。
両辺をk^0.5で割ると(k^0.5)の2次方程式になるので、解の公式から(k^0.5)、
つまりそれを二乗してkが求まります。


なので、
>>結局ckとnuSgsは1乗で比例していないのでは？
については、
「nuSgsは、(k^0.5 * Δ)に対してはckの1乗を比例定数として比例している」と言えると思います。

ちなみに、テンソルに関しては、OF付属のプログラマーズガイドに目を通してみることをおすすめします。
オープンCAE学会では日本語訳も公開しております。
http://www.opencae.jp/wiki/%E3%82%BD%E3%83%95%E3%83%88%E3%82%A6%E3%82%A7%E3%82%A2%E3%83%9E%E3%83%8B%E3%83%A5%E3%82%A2%E3%83%AB%E7%BF%BB%E8%A8%B3

2013年1月24日木曜日 9時41分10秒 UTC+9 nakagawa:

[Kohaku Ginga]

Gingaです。みなさまはじめまして。

横から失礼します。

興味があり拝読させていただいておりました。

kazuiさんが問題にしていたckの聞き具合について、特に圧縮性の場合ですが。。。

（以下、長文失礼します）

> なので、
> >>結局ckとnuSgsは1乗で比例していないのでは？
> については、
> 「nuSgsは、(k^0.5 * Δ)に対してはckの1乗を比例定数として比例している」と言えると思います。

からもうちょっと考察を進めると。。。

まず、上記では

nuSgs ∝ ck (k^0.5 * Δ)

と仰っているとのことでよろしいでしょうか？その前提で考えると

> 2/3*k*D:I - 2*ck*k^0.5*Δ*D:dev(D) + ce*k^1.5/Δ = 0
> となります。
> 両辺をk^0.5で割ると(k^0.5)の2次方程式になるので、解の公式から(k^0.5)、
> つまりそれを二乗してkが求まります。

とのことなので、(k^0.5)で割った方程式を一元二次方程式とみなすと係数はそれぞれ

a= ce/Δ

b= 2/3*D:I

c=  - 2*ck*Δ*D:dev(D)

となるので、解の公式から

(k^0.5) = (-b ± √(b^2-4*a*c))/(2*a)

ですよね。ということは、ざっくりと考えて。。。

√の項を(b^2-4*a*c)^0.5 ⇒ ((1+x)^0.5)/bと考えて（ただし、x=-4*a*c/b^2）

ここで、x<<1と考えてしまってもおおよそ妥当とすると

(1+x)^0.5 ≒ 1+0.5*x

であるので

c ∝ ck　だから　⇒　x ∝ ck

と言うことができて

(k^0.5) ∝ ck

であるので、先にでてきた式

nuSgs ∝ ck (k^0.5 * Δ)

より

nuSgs ∝ ck^2

と考えることはできないでしょうか？

即ち、圧縮性の場合

『nuSgsに対して、ckは自乗で効いてくる』

との結論も結構妥当な線をいっているような気がしますが、皆様のご意見を頂戴したいです。

2013年1月24日木曜日 10時48分37秒 UTC+9 ONO Hiroki:

[ONO Hiroki]

小野です。

解の公式から、
(k^0.5) = (-b ± √(b^2-4*a*c))/(2*a)

ですが、たとえば非圧縮ではb=0であるため、
(k^0.5) = √c / √a　、つまり　(k^0.5) = (ck^0.5)*(others) となります。
よって、先の解説の通り
nuSgs = ck^1.5 *(others)
であるといえますね。

圧縮性の場合でも、あまり密度変化が大きくない場合はbが小さくなるので(x<<1が成立しない)、
上記の非圧縮性の場合にかなり近くなると思います。


圧縮性が強い場合は、そもそも「ckの何乗で」と解釈するのはすこし無理がありそうな気がします。
とはいえ、私は強い圧縮性流れを扱った事が無いので経験的な観点からはわからないのですが…。



2013年1月24日木曜日 23時56分02秒 UTC+9 Kohaku Ginga:

[Kohaku Ginga]

Gingaです。

ご意見ありがとうございます。

私は乱流モデルについては初学者なので、先日の私の考察？は各係数や項のオーダーの「感覚」が全く無い中での論理展開でした。小野さんのようなお話を聞けることはとても参考になります。

＞　圧縮性が強い場合は、そもそも「ckの何乗で」と解釈するのはすこし無理がありそうな気がします。

仰るとおり、私もこの乱流モデルの特徴をckの効きだけで解釈することは適切でないと思っています。ただ、ckがどのくらい効いているのかなという「感覚」を持っておきたいなと思っていた次第です。

ということで、もうちょっと正確な議論にするためには、一度、注目している流れ場をLESで解いてみて

＞　2/3*k*D:I - 2*ck*k^0.5*Δ*D:dev(D) + ce*k^1.5/Δ = 0

の各項（因子）のオーダー比較をすると面白いかなと思います。

例えば、xの分布を計算させておいてparaViewで眺めるのも良いかもしれません。

私は、圧縮性流れ（遷音速領域）を扱うことが多いので、気になる話題であり続けるような気がします。

ちょっと時間を作ってあれこれ試してみたいと思います。

なるほど、圧縮性の場合は、bのオーダーによって

2013年1月25日金曜日 13時47分16秒 UTC+9 ONO Hiroki:

[kazui]

中川様　小野様　高木様　Ginga様

kazuiです。　おはようございます。

皆様の活発な議論を元に、昨日プログラマーズガイドとソースコード（compressibleのSmagorinsky.C）を

私なりに読んでみました。

数学的に数式を展開して、添付ファイルにまとめてみました。

結論としては、乱流粘性μsgsはCkの1〜3/2乗に比例しているという見解になりました。

確かに、単純にμsgsはCkの何乗に比例すると、明確には言えないようです。

私の見解に誤りがありましたら、ご指摘頂ければ幸いです。

よろしくお願いします。