流れ解析

修士論文要旨（2006年度）

安定化有限要素法に基づく

乱流モデルの構築と適用研究

Development and Application of LES Turbulent Model Based on Stabilized Finite Element Method

土木工学専攻29号 仲嶋 晋一 NAKAJIMA Shin-ichi

1.

研究目的

工学的な流れや，自然界の流れの多くは乱流状態にあ る．乱流の予測には，実験や観測が行われてきたが，近 年計算機の発達に伴い，乱流の予測に対して計算機によ る数値解析が多く用いられるようになってきた．これま で，乱流の数値解析手法には構造格子を用いた差分法に よる方法が多く用いられていた．しかし，これらの手法 は複雑形状への適用に限界がある．このため，複雑形状 に適した精度の高い乱流の数値解析手法の必要性が高 まってきている．一方，非構造格子に基づく有限要素法 は，複雑形状への適用性に優れており，それに複雑な非 定常解析に適したLES（Large Eddy Simulation）を適 用したモデルを構築することは工学的に重要である．

そこで，本論文では安定化有限要素法に基づくLES 乱流モデルを構築し，その精度検証を行うと共に，複雑 な流れ場への適用性の検討を行ったものである．数値解 析例には，Channnel流れ解析を取り扱い，既存のDNS

（Direct Numerical Simulation）データベース¹⁾との比 較により精度検証を行うと共に，立方体周り風況解析を 取り扱い，実験値²⁾との比較により複雑な流れ場への適 用性の検討を行った．

2.

支配方程式

非圧縮性・粘性流体の流れは，Navier-Stokes方程式 及び連続式によって支配されている．これらに対して フィルター操作を施すことによりLESの支配方程式が 導出される．LESの支配方程式は式(1)，(2)で与えら れる．

∂ui

∂t +uj

∂ui

∂x_j −fi =−1 ρ

∂P

∂x_i + ∂

∂x_j

¡−τij+ 2νDij¢

(in Ω) (1)

∂ui

∂xi

= 0 (in Ω) (2)

ここで，uiはi方向の流速のGS(Grid Scale)成分，fi

はi方向の物体力，ρは流体の密度，νは流体の動粘性係 数である．P は総圧，Dij は変形速度テンソルである．

式(1)中のτij はSGS(Sub-Grid Scale)応力項であり，

式(3)で与えられる．

τ_ij=u_iu_j−u_iu_j≡L_ij+C_ij+R_ij (3) ここで，L_ij，C_ij，R_ijはそれぞれ，Leonard項，Cross 項，SGS Reynolds Stress項と呼ばれている．

3.2δ 6.4δ

2δ

x y

z

図– 1 解析領域

表– 1 メッシュ分割

メッシュ分割数 メッシュ分割幅 (x×y×z) (∆x⁺，∆y⁺，∆z⁺) Mesh1 32×40×32 30，1∼14，15 Mesh2 32×64×32 30，0.2∼9.2，15 Mesh3 64×40×64 15，1∼14，7.5 Mesh4 128×40×64 7.5，1∼14，7.5

LESでは，式(3)中のuiujを直接計算することがで きないので，何らかのモデル化を行い計算を行う．本論 文では，SGS応力項をSmagorinskyモデル³⁾によりモ デル化を行った．Smagorinskyモデルにおけるτijは式 (4)で表される．

τij=−2νeDij , νe= (Csfs∆)²¯¯D¯¯ (4) ここで，Csは，Smagorinsky定数である．∆，¯¯D¯¯はそ れぞれ式(5)で与えられる．

∆ =Ve¹³ , ¯¯D¯¯= q

2D_ijD_ij (5) ここで，Veは四面体要素の体積である．

fsは壁関数であり，本論文ではVan Driest関数を採 用した．Van Driest関数は，式(6)で表される．

fs= 1−exp−y⁺

A⁺ , A⁺= 25, y⁺ =uτy ν (6) ここで，y⁺は壁座標，yは壁面からの距離，uτは壁面 摩擦速度である．

本論文では式(1)，(2)に対して，SUPG/PSPG法⁴⁾ に基づく安定化有限要素法を適用し，複雑形状への適合 性に優れたP1/P1要素（流速・圧力1次補間四面体要 素）を用いて空間方向に離散化し，有限要素方程式を得 た．時間方向への離散化手法としては，陰的解法である Crank-Nicolson法を用いた．．

1 10 100 0

10 20

y⁺

MeanVelocity

PlandtlLaw DNS(DataBase) Mesh1 Mesh3 Mesh4

(a)平均流速図 ⁰_{0 50 100 150}

1 2 3 4

y⁺

rms

DNS(DataBase) Mesh1 Mesh3 Mesh4

(b) 乱流強度図 図– 2 主流方向のメッシュ解像度による比較

1 10 100

0 10 20

y⁺

MeanVelocity

PlandtlLaw DNS(DataBase) Mesh1 Mesh2

(a)平均流速図 ⁰_{0 50 100 150}

1 2 3 4

y⁺

rms

DNS(DataBase) Mesh1 Mesh2

(b) 乱流強度図 図– 3 鉛直方向のメッシュ解像度による比較

1 10 100

0 10 20

y⁺

MeanVelocity

PlandtlLaw DNS(DataBase) case1-1 case1-2

(a)平均流速図 0 50 100 150 0

1 2 3

y⁺

rms

DNS(DataBase) case1-1 case1-2

(b) 乱流強度図 図– 4 安定化手法による比較

1 10 100

0 10 20

y⁺

MeanVelocity

PlamdtlLaw DNS(DataBase) case2-1 case2-2

(a)平均流速図 0 50 100 150 0

1 2 3

y⁺

rms

DNS(DataBase) case2-1 case2-2

(b) 乱流強度図 図– 5 移流速度の取り扱いによる比較

3. Channel

流れ解析

数値解析例として図−１に示した2δの流路幅を有す るChannel流れ解析を取り扱った．この解析例では，す べての物理量が無次元量であるため，支配方程式をδ， u_τ，ρで無次元化して数値解析を行った．無次元化され た方程式を式(7)，(8)に示した．

∂u⁺_i

∂t^∗ +u⁺_j ∂u⁺_i

∂x^∗_j −Kpδi1

=−∂P⁺

∂x^∗_i −∂τ_ij⁺

∂x^∗_j + 1 Re_τ

∂

∂x^∗_j Ã∂u⁺_i

∂x^∗_j +∂u⁺_j

∂x^∗_i

!

(in Ω) (7)

∂u⁺_i

∂x^∗_i = 0 (in Ω) (8)

ここで，上付き文字+はν，uτ，ρによる無次元化を，上 付き文字∗^はδによる無次元化をそれぞれ表す．流れを 保持するために式(1)におけるf1 = 1.0を外力項とし て与えた．境界条件は，Channelの上下面にNon-Slip 条件，それ以外の面では周期境界条件を採用した．初期 条件は，IwamotoらのDNSデータベース¹⁾に乱れを付 加したものを使用した．この初期条件を用いて，無次元 時間で40まで解析を行い，その解析結果を初期条件と し，再度解析を行うことにより，付加した乱れを除去し た．Reynolds数は，式(9)のように定義される．

Re_τ = u_τδ

ν = 150 (9)

微小時間増分量 ∆t = 5.0×10⁻⁴，Smagorinsky定数 C_s= 0.10とした．なお本論文では，計算コスト削減の ため，領域分割法に基づく並列計算法を採用した．通信 ライブラリーにはMPIを使用し，領域分割数は主流方 向に8等分割，横方向に2等分割とした．連立一次方程 式の解法にはElement-by-Element Bi-CGSTAB2法を 用いた．

本論文では，メッシュ解像度が精度に与える影響を考 察するため，表−１に示したメッシュを用いて解析を 行った．なお，初期条件を作成する際に付加した乱れの 倍率は，Mesh1は600%，Mesh2は1000%，Mesh3は 1500%，Mesh4は1500%とした．また，安定化手法を 用いたことによる数値粘性による影響を考察するため，

SUPG，PSPG 法の両手法を用いたケース，PSPG法 のみを用いたケース，SUPG，PSPG法の両手法を用い なかったケースの3ケースの解析を行った．さらに，移 流速度の取り扱いによる影響を考察するため，2次精度 Adams-Bashforth法を用いたケース，要素を構成する 4節点の流速の平均値を用いたケースの2ケースの解析 を行った．

3.2 数値解析結果

以下に数値解析結果を示す．図２〜５の(a)，(b)は それぞれ平均流速図，乱流強度図であり，無次元時間 で40 ∼80における値を用いた．なお，安定化手法の 有無による比較及び移流速度の取り扱いによる比較では Mesh1を用いた．

3.2.1 主流方向のメッシュ解像度による比較

図−１(a)では，メッシュ解像度を主流方向へ増加す ることによる精度の向上が見られ，特にMesh4につい てはDNSデータベースと良い一致を示した．

図−１(b)では，主流方向の乱流強度について，ピー ク値に関しては，メッシュ解像度の増加による精度の向 上が見られたが，ピークを過ぎてからの値はMesh4よ りもMesh3の方がDNSデータベースに近いという結 果となった．主流方向以外の乱流強度については，主流 方向の乱流強度の増加による精度の向上が見られた．

3.2.2 鉛直方向のメッシュ解像度による比較

図−２(a)では，両解像度による異差は見られなかっ た．これは，メッシュ解像度を増加した方向が流れに鉛 直な方向であったためであると考えられる．

図−２(b)では，全方向の乱流強度においてMesh1 に比べてMesh2の方がDNSデータベースに近いとい う結果となった．これは，流れに鉛直な方向にメッシュ 解像度を増加したことにより，主流方向以外へのエネ ルギーの伝達がより適切になったためであると考えら れる．

3.2.3 安定化手法の有無による比較

各ケースについてはそれぞれ，case1-1はSUPG法，

PSPG法の両手法を用いたもの，case1-2はPSPG法の みを用いたもの，case1-3はSUPG，PSPG法の両手法 を用いなかったものとした．

図−３(a)では，case1-1，case1-2の間に大きな異差 は見られなかった．case1-3については，計算が途中で 発散してしまう結果となった．

図−３(b)では，主流方向の乱流強度について，大き な差異は見られなかった．主流方向以外の乱流強度につ いては，case1-2に比べcase1-1の方がDNSデータベー スに近い結果となった．case1-3については，計算が途 中で発散してしまう結果となった．

3.2.4 移流速度の取り扱いによる比較

各ケースについてはそれぞれ，case2-1は移流速度に 対して2次精度Adams-Bashforth法を用い，case2-2は 移流速度に要素を構成する4節点の流速を平均した平均 流速を用いた．

図−４(a)では，case2-2と比べ，ややcase2-1の方が DNSデータベースに近い結果となった．

図−４(b)では，主流方向の乱流強度について，ピー ク値はcase2-1に比べ，case2-2の方がDNSデータベー スに近い結果となったが，case2-2はピーク位置が流れ

5.2L

9.6L

15.7L 5.2L

4.8L

図– 6 解析領域

1.0L

1.0L 1.0L 1.0L

1.0L O 中心線

y x z

LineA

LineB LineC

図– 7 立方体拡大図

の中心方向へずれてしまった．主流方向以外の乱流強度 については，case2-1の方がcase2-2に比べDNSデータ ベースに近い結果となった．case2-2ではピークが現れ なかった．

4.

立方体周り風況解析

物体周りの流れ場での有効性を検討するため，数値解 析例として図−６立方体周り風況解析を取り扱い，実験 値及び持田らが行ったRANSによる解析結果²⁾との比 較を行った．立方体部分の拡大図を図−７に示した．立 方体の一辺の長さを代表長さとし，原点を立方体底面中 心に設定した．総節点数69004，総要素数381024とし た．境界条件は，流入部に式(10)に基づく流速を与え，

両側面及び上面はSlip条件とし，底面及び立方体表面は Non-Slip条件とした．流出部は自由流出とした．

u=y^1/4

v=w= 0.0 (10)

Reynolds数は8.4×10⁴とし，微少時間増分量は5.0× 10⁻⁴ とし，Smagorinsky 定数は0.1 とした．初期条 件は，全節点でu = v = w = 0.0 とした。時間方向 への離散化は，移流速度に2次精度Adams-Bashforth 法を用いた．本解析例についても，計算コスト削減の ために領域分割法に基づく並列計算法を採用した．領

域分割は greedyに基づき行い，各領域の総要素数が

等しくなるようにした．連立一次方程式の解法には，

Bi-CGSTAB2 法を用いた．測定地点は，x = 0.2L， y = 1.0L〜1.5L，z = 0.0LをLineAとし，x= 1.0L， y= 0.1L，z=−0.25L〜1.0LをLineBとし，x= 1.0L，

0 1 1

1.1 1.2 1.3 1.4 1.5

u

y

Experiment Reference(RANS) LES

図– 8 LineA

-0.5 0 0.5 1

0 0.5 1

Experiment Reference(RANS) LES

u

z

図– 9 LineB

y = 0.0L〜1.5L，z= 0.0Lとし，無次元時間10∼40 の流速uの時間平均値を測定し，比較を行った．

4.2 数値解析結果

図−８にLineAでの時間平均流速，図−９にLineB での時間平均流速，図−１０にLineCでの時間平均流 速を示した．

図−８では，yが1.1から1.15の間で不連続な挙動 が見られた．これはメッシュ解像度不足に起因して発生 したと考えられる．全体的に流速を過小評価する結果と なった．

図−９では，立方体後方(z =−0.5L∼0.5L)の流速 について，持田らの解析結果とよい一致を示した．立方 体後方の逆流も再現することができた．立方体から離れ るにつれて流速を過大評価するという結果となったが，

定性的にはよい一致を示したと言える．

図−１０では，立方体後方(y = 0.0L∼1.0L)の流速 について持田らの解析結果とよい一致を示した．底面部 付近での循環流も再現することができた．立方体天板を 超えたあたりで流速を過小評価する結果となったが，定 性的にはよい一致を示したと言える．

5.

結論と今後の課題

本論文では，安定化有限要素法に基づくLES乱流モデ ルの構築及び精度検証を行い，複雑な流れ場への適用性 の検討を行うすることを目的とし，既存のDNSデータ ベース¹⁾及び実験値²⁾との比較のもとで検討を行った．

0 1

0 0.5 1 1.5

Experiment Reference(RANS) LES

u

y

図– 10 LineC

その結果，以下の結論を得た．

• メッシュ解像度を十分に増加することにより，

DNSデータベースの値によい一致を示し，解像 度を増加する方向については主流方向がより影響 が大きいことが示された．

• 安定化手法の有無による比較において，SUPG法 の有無については大きな異差は見られなかった が，PSPG法を用いなかったケースについては計 算が発散する結果となった．

• 移流速度の取り扱いによる比較においては，大き な異差は見られなかった．

• 立方体周り風況解析において，本計算結果は実験 値及びRANSによる計算結果と定量的によい一 致を示し，複雑な流れ場での有効性を示すことが できた．

今後の課題として，立方体周り風況解析においてメッ シュ解像度による影響を考察すると共に，さらに複雑な 流れ場への適用が挙げられる．

参考文献

1) K.Iwamoto,Y.Suzuki and N.Kasagi: Reynoids number effect on wall turbulence, International Journal of Nu- merical Method for Heat and Fluid Flow,23 pp.678- 689, 2002.

2) 持田灯，村上周三，林吉彦：立方体モデル周辺の非等方乱 流場に関するk−ϵモデルとLESの比較，日本建築学会 計画系論文報告集，第423号，pp.23−31,1991.5 3) J. Smagorinsky: General circulation experiments

with the primitive equation I. the basic experiment, Monthly Weather Review,91, pp.99-164, 1963.

4) T. E. Tezduyar: Stabilized finite element formulations for incompressible flow computations, Advance in Ap- plied Mechanics,28, pp.1-44, 1992.