發展一利用非結構性網格之可壓縮/不可壓縮黏性流計算方法(III)

行政院國家科學委員會專題研究計畫 成果報告

發展一利用非結構性網格之可壓縮/不可壓縮黏性流計算方

法(III)

計畫類別： 個別型計畫

計畫編號： NSC93-2212-E-009-012-

執行期間： 93 年 08 月 01 日至 94 年 07 月 31 日

執行單位： 國立交通大學機械工程學系(所)

計畫主持人： 崔燕勇

報告類型： 精簡報告

處理方式： 本計畫可公開查詢

中 華 民 國 94 年 9 月 19 日

PDFCMD Trial Version

www.zeon.com.tw

發展一利用非結構性網格之可壓縮/不可壓縮黏性流計算方法(III)

A Method Development for Compressible/Incompressible Viscous Flow

Using Unstructured Grids (III)

計劃編號:

NSC 93-2212-E-009-012

執行期限: 93 年 8 月 1 日 94 年 7 月 31 日

主持人: 崔燕勇 國立交通大學教授

一、 中文摘要(關鍵詞: 非結構性網格, 非 線性紊流模式) 本計劃係一系列有關非結構性網格計算 方法之第三期研究,本期計劃主要是針對紊 流場計算,除了採用常用的線性渦流模式,還 推導出一通用的二次渦流模式,利用此通用 式,一般常用的非線性渦流模式皆可以此式 表示,配合前期發展的非結構性網格計算方 法,將其應用至數個基本的管道流並與實驗 值比對。

英 文 摘 要(Keywords: unstructured grids,

nonlinear turbulence models) This study is part III of a series of projects related to the method development using unstructured grids. The main concern of this year’s s project is the turbulent flow. In order to characterize turbulence a general form of second-order nonlinear eddy-viscosity model is derived. It can be seen that the linear and a number of nonlinear models can be represented by this form. A number of basic turbulent duct flows have been examined using the nonlinear models. The results have been compared with measurements.

二、計劃緣由與目的 在求解紊流流場時﹐因為對大部份的紊 流流場來說﹐直接解 Navier- Stoke 方程式 仍是不實際的﹐因此藉由紊流模式來求解紊 流流場仍是必需的工具。目前最普遍使用的 紊流模式為線性的 k-ε模式﹐主要是因為其 簡單的形式﹐計算上較穩定﹐而且在計算上 所花費時間較短﹐能有效率地得到合理的解。 但當流場的複雜性增加時﹐如流線的曲度增 加、強的渦流﹐線性 k-ε模式就不再適用 了。一般認為不適用的原因為Boussinesq 假 設雷諾應力與應變率為線性關係。使用雷諾 應力模式可以解決一些線性 k-ε模式不適 用的問題﹐但是相對的使用雷諾應力模式複 雜許多﹐且花費的計算時間也相當長。另一 個解決的方法就是使用非線性的雷諾應力與 應變率關係。最近的研究顯示採用非線性 k- ε模式﹐提供了一個兼顧效率與準確性的方 法。 在本文中使用非結構性網格﹐與結構性 網格比較﹐結構性網格因為具有自動指標性 及優先的方向﹐所以對於格點之間的相關位 置容易得知﹐使用起來也因此比較簡單﹐但 是對於較複雜的幾何形狀﹐則往往受到限 制。因此近年來發展出非結構性網格對於複 雜幾何形狀﹐因為非結構性網格可以是任意 形狀的多邊形﹐所以可以簡單的建立複雜幾 何外型的網格﹐並且也能容易的做到局部加 密。但是其在計算時需要儲存格子之間的相 對關係﹐並且在計算上也較為複雜﹐所以比 起結構性網格需要較大的記憶空間及計算速 度﹐然而因為電腦的快速發展非結構性網格 已經比結構性網格佔有更大的優勢。 三、研究方法 非線性紊流模式 在 eddy-viscosity 模式中係假設雷諾應 力與平均速度場的應變呈一多項式關系，在

PDFCMD Trial Version

www.zeon.com.tw

此假設為一二次式 : 0 1 2 ' ' ~ i i k i j j j j u u u u u a a a x x x ∂ ∂ ∂ − + + ∂ ∂ ∂ 將應變率拆成二項

(

ij ij

)

j i S x u Ω + = ∂ ∂ 2 1 其中 ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ∂ ∂ + ∂ ∂ = i j j i ij x u x u S 2 1 ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ∂ ∂ − ∂ ∂ = Ω i j j i x u x u 2 1 在此必須注意，Sij 為一對稱張量，而 Ωij 係反對稱(anti-symmetric)。利用 S_ij 與 ij Ω 可將雷諾應力關係式拆成下列數項： 0 11 12 21 22 22 23 ' ' i j ij ij ik kj ik kj ik kj ik kj u u a a S a a S S a S a S a − = + + Ω + + Ω + Ω + Ω Ω 其中每一項給予一單獨係數，雷諾應力項必 須要滿足下列二性質： (1) 對稱性(symmetry)： −ui'uj'=−uj'ui' 由此可得到： 22 22 12 0 a a a − = = (2) 當 i= j 時等式左右二邊須滿足一致性 (consistency)： 由此可推論： ij i k u a δ 3 2 '2 0 =− =− 同時其它項必須為 0，為滿足此一條件， ik

S

S

_kj及

Ω

_ik

Ω

_kj項必須由

(

S

_ik

S

_kj

)

* 及

(

Ω

ik

Ω

kj

)

* 取代，此處： ij kl lk kj ik kj ik ij kl lk kj ik kj ikS S S S S S δ δ Ω Ω − Ω Ω ≡ Ω Ω − ≡ 3 1 ) ( 3 1 ) ( * * 值得注意的是： _{( )}* kj ik kj ik kj ik kj ik S S S S Ω −Ω ≡ Ω −Ω 如此一來，雷諾應力可表示成下列形式： _{− =− + +} ij ij j iu k a S u' ' ₁₁ 3 2 δ * 23 * 22 * 21(SikSkj) a (Sik kj ikSkj) a ( ik kj) a + Ω −Ω + Ω Ω 上式即為Craft 等人之三次模式中的一次及 二次項，同時也就是Speziale 經過簡化後的 二次模式的形態，倘若只保留一次項，即為 所熟知的線性模式，另外可以經由推導將上 面的二次項改寫如下： * 3 * 2 * 1 * 23 * 22 * 21( ) ( ) ( ) ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ∂ ∂ ∂ ∂ + ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ ∂ ∂ ∂ ∂ + ∂ ∂ ∂ ∂ + ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∂ ∂ ∂ ∂ = Ω Ω + Ω − Ω + j k i k i k k j j k k i k j k i kj ik kj ik kj ik kj ik x u x u C x u x u x u x u C x u x u C a S S a S S a 此處 ij k l k l k j k i k j k i x u x u x u x u x u x u δ ∂ ∂ ∂ ∂ − ∂ ∂ ∂ ∂ = ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ ∂ ∂ ∂ ∂ 3 1 *

上面形式與Rubinstein & Barton 及 Yoshizawa

模式相同。 從以上的推導可以看出，利用此通用形 式只須改係數即可得到各種不同的二次紊流 模式，在前面所給的兩種表示法我們將採第 一種形式，因為如此在未來能較容易的將其 延伸到三次模式。 非結構性網格方法 在有限體積法中，首先將輸送方程式對 控 制 容 積 積 分 ， 利 用 高 斯 的 散 度 定 理 (divergence theorem)將對流及擴散項轉化為 面積分，再利用均值定理(mean-value theorem) 轉為差分形式。對某一面 ”e” 對流項的計 算採用下列的混合差分計算：

(

)

c UD HD UD e e e e F =F +γ F +F 此處上標UD 代表上風差分法，另一上標 HD 代表一高階差分法，

γ

乃是一介於 0 與 1 之 間的數值。 擴散項的表示可有許多表示法，但這些 方法多半在從二維問題擴展至三維時須重新 推倒，且其差分方程式變得複雜許多，下列 所用的方法卻完全不需考慮其為二維或三 維。_δr_PE為一從點P 至點 E 之向量，sr是界於 P 及 E 之間定容面”e”的表面向量，擴散通 量可由下式近似之： 2 ( ) ( ) d e e E P e e PE s F s d s φ φ φ δ Γ = − + Γ ∇ • − • r r r _r 在本方法中網格位置排成非交錯式，所 有變數儲存在控容中心，在上述對流項及在 後面所推導的壓力修正方程式中，需要計算 通過控容面上的質量流率，此質通量可由下 面方式推估。先計算面上之速度： ( ) e e e E e e Vr =Vr − D∇p −D ∇p 上式中加有上標〝—〞者代表由共用此面相 鄰兩點P 和 E 內插得之，如此面上的質通量 可由下式近似求得： P[( ) ] e e e E E P e P E m& = ρ Vr • −s A p − p − ∇p •δr 令所有流過一控容面上的質通量滿足質量守 衡，可得壓力修正方程式： 1 2 1 n P P C C P P C A p A P S S = ′ =

∑

+ +

PDFCMD Trial Version

www.zeon.com.tw

其中 * 1 1 n P f f S m = =

∑

₂ 1 ( ) n P f f f f S ρ D P s d = ′ =

∑

∇ • −r r 四、結果與討論 A.完全發展管道流 首先考慮於一管道中的完全發展流,對 於平均速度的計算結果﹐線性與非線性模式 跟實驗的值比較都相當接近。 正向應力：圖 1 為正向應力分佈圖﹐ 由此圖可看出在完全發展流中﹐線性模式所 預測的雷諾應力為等向性_{τxx =τyy =τzz}﹐但由 實驗結果發現有相當的誤差﹐而非線性模式 所預測的_τxx比起線性的模式大﹐所預測的 yy τ ﹐_τzz則比線性模式小﹐都比較符合實驗的 結果。由此可看出非線性模式中所加入的非 線性項能使預測的正向應力更準確。 B.背階管道流 下表列出利用線性模式非與線性模式所 預測的於一突張管道其背階後方迴流接觸長 度﹐由此表可以看出﹐不管線性或非線性模 式對於接觸長度的預測都比實驗值短﹐而非 線性模式預測的接觸長度比線性模式長﹐較 接近實驗值﹐其誤差少了近 5%。 紊流模式 接觸長度 誤差 實驗值 Kim(1980) 7.11H 線性 k-ε模式 6.05H -14.91% Speziale 模式 6.40H -9.98% Rubinstein 模式 6.38H -10.27% 正向應力與剪應力：圖 2 為正向應力與 剪應力分佈圖,可看出線性與非線性模式對 於

u

'u

'

預測都很接近﹐大部分與實驗值也很 吻合﹐只有在迴流區(x<7.11H)﹐線性模式與 非線性模式預測的值都略大於實驗值﹐而在 這裡非線性模式所預測的值比較接近實驗。 為

v

'v

'

的分佈圖﹐在此圖中可以發現在迴流 區中﹐線性模式明顯預測值過大﹐而非線性 模式所預測的值較小,與實驗值比較接近﹐這 與完全發展流的

τ

yy類似﹐但是非線性與線性 模式對於

v

'v

'

最大值的位置預測的都比實驗 低。而在迴流區之後逐漸恢復成完全發展流 的區域(x>8.5H)﹐線性模式則跟實驗值相當 吻合﹐而非線性模式所預測的值反而過小。 為剪應力分佈圖﹐在迴流區中線性與非線性 模式預測的剪應力明顯都過大﹐而在迴流區 之後則兩種模式都與實驗值相近。 C.軸對稱擴張管道流 下表列出對於一軸對稱擴張管在不同的 擴張角下﹐線性模式與非線性模式所預測的 迴流接觸長度﹐可看出非線性與線性模式鎖 預測的接觸長度都比實驗值來的短﹐而非線 性模式所預測的結果較佳﹐誤差比線性模式 少 3~4%左右。而當擴張角越大時﹐兩種模式 跟實驗值的誤差也越大。 (1)18o 擴張角 紊流模式 接觸長度 誤差 實驗值Stieglmeier 1.49 D2 線性 k-ε模式 1.257 D2 -15.64% Speziale 模式 1.301 D2 -12.68% Rubinstein 模式 1.307 D2 -12.28% (2)90o 擴張角 紊流模式 接觸長度 誤差 實驗值Stieglmeier 1.86 D2 線性 k-ε模式 1.347 D2 -27.58% Speziale 模式 1.425 D2 -23.38% Rubinstein 模式 1.437 D2 -22.74% 正向應力與剪應力：以 18 度擴張角為例 (圖 3),由

u

'u

'

的分佈圖發現﹐非線性與線性 模式預測的

u

'u

'

比實驗值來的小﹐而非線性 所預測的結果比線性模式好。但在擴張角為 90 度時﹐線性模式預測的仍小於實驗值﹐但 非線性模式所預測的則大於實驗值。由

v

'v

'

的分佈圖可發現﹐在 18 度中非線性模式線性 模式所預測值較小﹐與實驗值比較接預測的

'

'v

v

在 x>40mm 之後比起實驗值過小﹐反而是 線性模式較為接近實驗值。在 90 度時則不太 一樣﹐線性模式預測的都比實驗值大﹐而非 線性模式所預測的結果較佳﹐這跟完全發展 流中

τ

_yy的情形類似。由

w

'w

'

的分佈圖可發 現,在 18 度線性與非線性模式對

w

'w

'

預測都 很接近與實驗值相比﹐除了在 x=20mm 與 x=40mm 的地方稍略小於實驗值﹐其他地方於 實驗值相差不大。而在 90 度的流場中﹐線性

PDFCMD Trial Version

www.zeon.com.tw

模式預測的

w

'w

'

就稍微大於實驗值﹐而非線 性所得到的結果就比較接近實驗﹐這跟完全 發展流中

τ

zz的情形類似。

u

'v

'

的分佈圖可發 現在 18 度的流場中﹐線性模式預測值則明顯 過大﹐而非線性模式與實驗值相比則有很好 的結果。在 90 度時﹐兩種模式預測的結果相 近。 五、參考文獻

1. Craft, T. J., Launder B. E. and Suga, K. 1993 “Extending the applicability of eddy viscosity models through the use of deformation invariants and non-linear elements”, 5th IAHR conference on Refined-Flow Modeling and Turbulence Measurement Paris 7-10

2. Rubinstein, R. and Barton, J. M. 1990 “Nonlinear Reynolds stress models and the renormalization group”, Phys. Fluids A, Vol. 2, No. 8, pp. 1472-1476

3. Speziale, C. G. and Ngo, T. 1998 “Numerical Solution of Turbulent flow past a backward facing step using a nonlinear k-ε model”, Int. J. Engng Sci. Vol. 26, No. 10, pp. 1099-1112

4. Yoshizawa, A. 1984 “Statistical analysis of the deviation of the Reynolds stress from its eddy-viscosity representation”, Phys. Fluids, Vol. 27, No. 6, pp. 1377- 1387 0 0.2 0.4 0.6 0.8 1 Y/d 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 [-τxx / ρ Uma x 2] 1/2 × 10 Experiment: Laufer(19 Linear k-ε model Speziale model Rubinstein model 0 0.25 0.5 0.75 Y/d 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 [-τyy / ρ Um ax 2 ] 1/ 2 × 10 Experiment:Laufer Linear k-ε model Speziale's model Rubinstein's mode 0 0.25 0.5 0.75 1 Y/d 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 [-τzz / ρ Um ax 2] 1/2 × 10 Experiment:Laufer(1951) Linear k-ε model Speziale's model Rubinstein's model 圖一 完全發展管道流正向應力分佈圖 0 0.025 0.05 u'u'/Uo 2 0 1 2 3 y/ H Experiment:Kim(1980) Linear model Speziale model Rubinstein model X=0.975H 0 0.025 0.05 u'u'/U_o2 0 1 2 3 y/ H X=6.778H 0 0.025 0.05 v'v'/U_o2 0 1 2 3 y/ H Experiment:Kim(1980) Linear model Speziale model Rubinstein model X=0.975H 0 0.025 0.05 v'v'/Uo 2 0 1 2 3 y/ H X=6.77H 0 0.01 0.02 0.03 u'v'/Uo 2 0 1 2 3 y/H X=0.975H 0 0.01 0.02 0.03 u'v'/U_o2 0 1 2 3 y/ H X=6.77H 圖二 背階管道流正向應力分佈圖 0 0.04 0.08 u'u'/Um2 0 10 20 30 40 y (mm) X=40mm 0 0.04 0.08 v'v'/Um2 0 10 20 30 40 y( m m ) X=40mm 0 0.04 0.08 w'w'/Um2 0 10 20 30 40 y (mm) X=40mm -0.005 0 0.005 0.01 0.015 0.02 0.025 0.03 u'v'/Um2 0 10 20 30 40 y (mm) X=40mm 圖三 18o _{軸對稱擴張管道流正向應力及剪} 應力分佈圖

PDFCMD Trial Version

www.zeon.com.tw