發展一非結構性網格有限體積法對固體結構及流固耦合計算分析

國 立 交 通 大 學

機 械 工 程 學 系

碩 士 論 文

發展一非結構性網格有限體積法對固體結構及流固

耦合計算分析

An Unstructured-Grid Finite Volume Method for Calculating

Structural Dynamics and Fluid-Structure Interaction

研 究 生：黃 義 政

指 導 教 授：崔 燕 勇 博士

發展一非結構性網格有限體積法對固體結構及流固耦合計算分析

An Unstructured-Grid Finite Volume Method for Calculating

Structural Dynamics and Fluid-Structure Interaction

研 究 生：黃義政 Student：Yi-Cheng Huang 指導教授：崔燕勇 Advisor：Yeng-Yung Tsui 國 立 交 通 大 學 機械工程學系 碩士論文 A Thesis

Submitted to Department of Mechanical Engineering Collage of Engineering

National Chiao Tung University in Partial Fulfillment of the Requirements

for the degree of Master of Science

In

Mechanical Engineering July 2011

Hsinchu, Taiwan, Republic of China

發展一非結構性網格有限體積法對固體結構及流固耦合計算分析

研究生：黃義政 指導教授：崔燕勇 博士 國立交通大學機械工程學系(研究所)碩士班 摘要 本文研究主要為發展一個非結構性網格有限體積法，應用於計算結構力學及 流體力學來探討流固耦合之問題。在結構力學計算中，透過有限體積法，將方程 式離散後整理成一組聯立常微分方程式，求解方式分成直接解法與虛擬時間法兩 類，再藉由龍格庫塔法或預測修正法求解。結果發現在虛擬時間法下的求解方式 較佳，其中又以預測修正解法有較快之收斂速度和較大之時間步階。此方法能有 效的運用在典型的二維懸臂樑和固定樑問題，且分析之結果與理論之位移、自然 頻率非常相近。由於固體結構的振動，在流場計算中採用移動網格，透過空間守 衡的概念來滿足移動網格的質量守衡。方程式的離散同樣採用有限體積法，其中 對流項採用上風法和中央差分法混合的高階方法。在處理動量與連續方程式間的 耦合，則是藉由先預測再修正的 PISO 法則處理。在耦合測試問題中，第一個探 討的是流體流經一垂直彈性平板之管流，平板因受流場壓力與剪力的影響而產生 擺動；第二個探討的是一位於矩形塊體後的水平彈性平板，由於受到此塊體不穩 定渦流的影響，造成平板的振動。在此二測試問題中，闡述了結構與流場間相互 影響的關係。

An Unstructured-Grid Finite Volume Method for Calculating

Structural Dynamics and Fluid-Structure interaction

Student：Yi-Cheng Huang Advisor：Dr. Yeng-Yung Tsui

Department of Mechanical Engineering National Chiao Tung University

ABSTRACT

A new unstructured-grid finite-volume method is developed for the structural dynamics and fluid-structure interaction. In structural calculations, the governing equations are discretized using the finite-volume method to obtain a system of ordinary differential equations. These system of ODEs are solved by either a direct method or a virtual time method using Runge-Kutta or predictor-corrector methods. It is shown that the virtual time method is superior to the direct method. Incorporating the predictor-corrector scheme with the virtual time method, larger time steps and faster convergence rates are obtained. The methods are tested for cantilever and built-in beams. The results are in good agreement with theoretical analysis in terms of displacements and natural frequencies. Due to the vibration of the structure, moving grids are employed in flow calculations. By incorporating the space conservation concept, the mass is conserved in each cell with moving grids. The equations are discretized also using the finite-volume method. The convection term is approximated by a hybrid scheme which mixes the central and upwind differences. The coupling between the momentum and continuity equations is tackled using the predictor-corrector method of the PISO algorithm. In the tests of fluid-structure interaction, the first case considered is a vertical plate swinging in a channel flow. In the second case, a horizontal plate is placed behind a rectangular block. The vortices caused by the flow over the block result in vibration of the plate. The interaction between the vibration of the structure and the vortex flow of the fluid is clearly delineated.

誌謝

在這兩年來的日子中，非常感謝崔燕勇教授的指導，使我成長很多，不管 是在學業方面還是在研究方面，都給予我很大的幫助，讓我從中學習到解決問題 的態度及方法，再次感謝崔老師在論文研究上的教導，使我能完成此篇論文。 接下來是實驗室全體，感謝胡育昌、林仕文、陳信宏、郭大慶學長在研究上 給予的協助，特別是林仕文學長在我遇到問題時，都會幫我解決困難並給予鼓 勵，使我能克服重重困難。接著是感謝同窗的林子翔和賴胤男，一起渡過兩年來 的時光。感謝學弟們吳奉起、黃裕堂、鄭凱鴻，使我們實驗室更加有活力。 最後要感謝的是我的家人，爸爸、媽媽、哥哥、弟弟，有你們的支持和鼓勵， 才能使我順利完成學業。

目錄

摘要...i ABSTRACT...ii 誌謝...iii 目錄...iiiv 表目錄...vi 圖目錄...vii 符號表 ...x 第一章、緒論... 1 1.1 前言...1 1.2 流固耦合介紹...1 1.3 文獻回顧...2 1.4 研究內容...5 第二章、數學模式... 7 2.1 固體運動模式...7 2.2 流體運動模式...8 2.3 流固耦合...9 第三章、固體力學數值方法... 11 3.1 有限體積法...11 3.2 內部梯度計算方式...12 3.3 邊界梯度計算方式...13 3.4 求解方法...13 3.4.1 直接解法 ...13 3.4.2 虛擬時間法 ...15 3.5 邊界上之計算...18

第四章、流體力學數值方法... 19 4.1 空間守恆定理...19 4.2 有限體積法...20 4.3PISO演算法 ...23 4.4 邊界條件...27 4.5 網格的移動...28 4.6 流固耦合演算流程...28 第五章、結果與討論... 29 5.1 固力分析...29 5.1.1 求解方式分析 ...29 5.1.2 自然頻率和阻尼 ...31 5.1.3 網格測試 ...31 5.1.4 疊代次數和計算時間分析 ...33 5.2 與ANSYS分析比較...33 5.3 不同形式外力負載之懸臂樑...34 5.4 流固耦合分析...35 5.4.1 流經彈性平板之管流 ...35 5.4.2 渦流致動之彈性平板 ...37 第六章、結論... 40 參考文獻... 41

表目錄

表(一)結構體的材料性質和幾何尺寸【6】 ...44 表(二)不同阻尼比相對應之阻尼係數 ...44 表(三)在不同網格數下之位移和誤差百分比...44 表(四) 在不同時間步階下比較計算時間、步階次數和使用之虛擬時間步階 ...45 表(五) 流經彈性平板之管流的結構參數【14】 ...45 表(六) 流經彈性平板之管流的流體參數【14】 ...46 表(七) 在不同雷諾數和平板厚度下之位移結果 ...46 表(八) 渦流致動之彈性平板的結構參數【7】 ...46 表(九) 渦流致動之彈性平板的流體參數【7】 ...46

圖目錄

圖 3.1 以網格頂點為控制體積之示意圖... 47 圖 3.2 位移梯度計算示意圖... 48 圖 3.3 邊界位移梯度計算示意圖... 48 圖 3.4 相鄰邊界的控制體積計算示意圖... 49 圖 4.1 以網格中心為控制體積之示意圖... 50 圖 4.2 利用各面新舊點位置來計算體積改變量... 50 圖 4.3 兩個時階中的網格位置變化和各面掃過之體積... 51 圖 4.4 主格點和鄰近格點之控制面示意圖... 51 圖 4.5 流場網格分割區塊圖(一)... 52 圖 4.6 流場網格分割區塊圖(二)... 52 圖 4.7 流場網格變化示意圖... 53 圖 4.8 流固耦合演算流程圖... 54 圖 5.1 單點受力之懸臂梁示意圖... 55 圖 5.2 均佈受力之兩端固定樑示意圖... 55 圖 5.3 懸臂樑尾端頂點的位移反應... 56 圖 5.4 懸臂樑尾端頂點在不同阻尼係數下的位移反應... 56 圖 5.5 兩端固定樑中心位置的位移反應... 57 圖 5.6 懸臂樑尾端頂點在阻尼係數 12000 下的位移反應... 58 圖 5.7 兩端固定樑中心位置在阻尼係數 10000 下的位移反應... 58 圖 5.8 懸臂樑尾端頂點的位移反應... 59 圖 5.9 兩端固定樑中心位置的位移反應... 59 圖 5.10 懸臂樑尾端頂點自由振動之位移反應... 62 圖 5.11 兩端固定樑中心位置自由振動之位移反應... 62 圖 5.12 懸臂樑尾端頂點在不同阻尼比之位移反應... 63

圖 5.13 兩端固定樑中心位置在不同阻尼比之位移反應... 63 圖 5.14 五種不同網格數(上到下分別是由疏到密)... 60 圖 5.15 測試網格精準度之誤差分析... 61 圖 5.16 在真實時間步階為 0.01 sec每個時階下收斂之疊代次數... 64 圖 5.17 在真實時間步階為 0.1 sec每個時階下收斂之疊代次數... 64 圖 5.18 單點受力之懸臂樑應力圖(分別為 σxx、σ 、σyy xy)，其中左邊為ANSYS 模擬、右邊為本文模擬... 65 圖 5.19 均佈受力之兩端固定樑應力圖(分別為 σxx、σ 、σyy xy)，其中左邊為ANSYS 模擬、右邊為本文模擬 ... 66 圖 5.20 不同形式分佈受力之懸臂梁示意圖... 67 圖 5.21 不同型式受力之懸臂樑應力圖(a)均勻性的分佈力(b)線性遞減的分佈力(c) 二次曲線遞減的分佈力……….…...68 圖 5.22 流經彈性平板之管流示意圖... 69 圖 5.23 渦流致動之彈性平板示意圖... 69 圖 5.24 彈性平板(厚度 0.2m)自由端頂點的位移反應(a)不同雷諾數下之比較(b)在 雷諾數 50 的局部放大圖(c) 在雷諾數 100 的局部放大圖(c) 在雷諾數 150 的局部放大圖 ... 70 圖 5.25 彈性平板(厚度 0.1m)自由端頂點的位移反應(a)不同雷諾數下之比較(b)在 雷諾數 50 的局部放大圖(c) 在雷諾數 100 的局部放大圖(c) 在雷諾數 150 的局部放大圖 ... 71 圖 5.26 彈性平板(厚度 0.2m)自由端頂點的位移反應(a)10 個週期(b)最後一個週期 ... 72 圖 5.27 彈性平板(厚度 0.1m)自由端頂點的位移反應(a)10 個週期(b)最後一個週期 ... 72 圖 5.28 平板厚度 0.2(m)，在不同週期下之流場流線和壓力圖... 73

圖 5.29 平板厚度 0.1(m)，在不同週期下之流場流線和壓力圖... 74 圖 5.30 在不同雷洛數下彈性平板尾端頂點的位移反應... 75 圖 5.31 雷諾數 110，在時間t=20(s)時流場流線和壓力圖... 75 圖 5.32 雷諾數 204，彈性平板擺盪一個周期之流場壓力圖... 76 圖 5.33 雷諾數 204，彈性平板擺盪一個周期之流場流線圖... 77 圖 5.34 雷諾數 204，彈性平板隨時間之變化圖... 78 圖 5.35 雷諾數 300，彈性平板擺盪一個周期之流場壓力圖... 79 圖 5.36 雷諾數 300，彈性平板擺盪一個周期之流場流線圖... 80 圖 5.37 雷諾數 300，彈性平板隨時間之變化圖... 81

符號表

符號 定義 σ 應力 d 位移向量 c 阻尼係數 E 楊氏系數 υ 蒲松比 ε 應變 C f F 面上的對流通量 D f F 面上的擴散通量 m 質量流率 f m 面上的質量流率 r f m 流體相對於面上的質量流率 P 流體壓力 f S 面的法向量 t 時間 τ 虛擬時間 V 流體速度向量 g V 網格速度向量 Pnb δ P 到 nb 距離向量 ζ 阻尼比 μ 黏滯係數 ρ 密度 ω 自然角頻率 f 頻率 f Δ∀ 一個面掃過的體積 I 面積慣性矩 上標 n 現在時間的未知數 o 前一時間的已知數 (m) 疊代中的上一個疊代值 ─ 主格點與鄰近格點之內差 ′ 修正量 下標

f 面上的值

nb 鄰近格點的值

第一章、緒論

1.1 前言

從 1950 年代開始，眾多不同的方法被發展出來，藉由這些方法來處理廣泛 物理現象的微分方程式。其中有限元素法(Finite element method)最先發展出來應 用於固體力學，而這項方法也延伸應用在熱傳、流體和電磁方面。而有限體積法 (Finite volume method)是從較早的有限差分法演變發展而來，並且被廣泛的運用 在計算流體力學(CFD)方面，如燃燒、自由表面流和多相流等…。儘管一開始有 限體積法被限定在使用結構性網格，然而此方法廣泛應用於 CFD 是因為不管在 局部或整體都能守恆其物理現象。最早使用非結構性網格於有限體積法起源於 1960 年代，在往後幾年也被證明有效的運用在流體、熱傳和相變化。不過在固體 力學的問題方面，多年來一直在有限元素法中發展，也有不錯的成果出現，因此 也陸續出現很多種分析固體力學的商業軟體。因此有限體積法應用到固體力學這 方面一直沒有更進一步的突破。 直到最近幾年來，陸續有一些學者開始以有限體積法的觀點來處理計算固體 力學(CSD)的問題，像是平板彎曲的分析、受不同負載的固力分析、彈性固體的 應力分析、暫態固體力學的分析和流固耦合分析等…。透過這些學者的努力，證 明 FV 方法在 CSD 中有很好的成效。在 FV 方法中對於 CSD 可分為兩大類型， 一種是將控制體積取在網格頂點(cell-vertex)【1~7】，另一種是將控制體積取在網 格中心(cell-centered) 【8~11】。其中計算上在相同網格數時，使用 cell-vertex 方 法比 cell-centered 方法減少較多的計算量和記憶空間。在另一方面，在對於比較 不規則網格時以 cell-vertex 方法來計算應力較佳。因此一般而言，以有限體積法 在處理固體力學問題時都以 cell-vertex 方法佔大多數。 1.2 流固耦合介紹 流體和固體之間的影響在一般的工程問題中扮演非常重要的角色，而這兩種 物體之間影響稱之為流固耦合。而本文研究之流固耦合即是透過流體對固體施加

一作用力，讓固體產生擺動，同時受到固體擺動會使流場發生變化。這種現象在 生活中到處可見，像是懸吊的橋、受風力影響的高樓和血液流經血管等問題。然 而這些物體受到流體的影響是不可忽略的，因此不去考慮它將可能會產生嚴重的 後果。 在許多科學的應用和日常生活中，都會碰到流固耦合(Fluid-Structure interaction，FSI)的問題。因此在這方面的研究上，學者們提出了一些解決方法來 近似流固耦合。目前在處理流固耦合的問題可以分為兩種解法，分別是完全統一 法(monolithic approaches)和分段法(partitioned approaches)。在完全統一法方面， 是將流體力學之統御方程式和固體力學之統御方程式結合用同一套程式且同時 進行求解。在分段法中，是將流體方程式和固體方程式分開進行求解，其中固力 部分一般以有限元素法求解，不過少部分也有以有限體積法解之。另外在處理流 固耦合問題時，流體部分的網格不是靜止的，因此不能再使用一般的 Eulerian 座 標系統，而必須要以 ALE (arbitrary Lagrangian-Eulerian)座標系統之統御方程式並 且搭配移動網格求解。此方法也已經被證明有效的運用在 FSI 的問題中。

1.3 文獻回顧

(A) 固力方面(FVM)

1995 年，Bailey and Cross【2】以非結構有限體積法求解彈性力學方程式， 其中結構體所受之負載分為溫度和力兩種形式。分別以有限元素法(FE)和有限體 積法(FV)測試四種例子，並且與解析解比較。由模擬測試後發現利用 FE 和 FV 方法求解之結果非常接近。另外因為 FV 矩陣所使用的方法為 node by node 所建 立，所以計算之時間將比 FE 方法多。

1999 年，Taylor et al.【3】以有限體積法來計算固體力學(CSD)，並且使用不 同形狀的網格做測試，分別為 linear tetrahedral(LT)、bilinear pentahedral(BLP)和 trilinear hexahedral(TLH)。其中發現到在使用 TLH 網格並在較粗網格下時，由結 果觀察可知 FV 方法比 FE 方法的準確性更佳。然而這些方法中又以使用 LT 網格

時其結果最好。在使用 LT 網格時，FV 和 FE 的矩陣為對稱矩陣，所以以 conjugate gradient method(CGM)解之，而在使用 TLH 網格時 FV 為非對稱性矩陣，所以需 以 bi-conjugate gradient method(Bi-CGM)解之。在計算時間方面，FV 方法所花費 時間較 FE 多，其中在使用 TLH 網格時，FV 所花費時間幾乎是 FE 的兩倍之多。

2004 年，Fallah【12】分別以 cell vertex 和 cell centred 有限體積法來分析平 板彎曲的問題，由結果觀察出 FV 方法在網格數的收斂速度比 FE 方法快。在有 限體積法中又以 cell centred 準確度比 cell vertex 好，另外在 Mindlin-Reissner 理論 中使用有限元素法來分析薄平板問題會有 locking 現象發生，不過在使用有限體 積法時卻沒觀察到此現象。

2007 年，Xia et al.【4】以隱性非結構有限體積法來模擬分析 2D 和 3D 的固 力學問題，此篇描述一個新的 cell-vertex FVM 方法適用在三角和四面體網格。在 測試問題中，計算懸臂樑在尾端受一集中力的變形，其模擬結果與理論解比較非 常相近。在另一個問題中(四邊形平板)，分別以 simply supported 和 fully clamped 固定，其模擬結果觀察到最大的彎曲應力發生位置不一樣，並且和理論解比較也 都非常相近。因此可將此方法與計算流體力學(CFD)做結合，來模擬流固耦合(FSI) 問題。

(B)流固耦合方面

2000 年，Souli et al.【13】 利用 Lagrangian 和 Arbitrary Lagrangian Eulerian (ALE)兩種方法來探討不同問題時，對於網格幾何的影響。在測試問題中，發現 以 ALE 方式求解較不會導致網格被嚴重破壞，以至無法繼續計算，且在求解中 使用的時間步階也不會受到網格變形的影響而可繼續進行計算。另外此篇對於較 大變形量的計算時，提出了幾種不同的 smoothing 演算法來從新規劃計算區域的 網格。在耦合方面，不同的 smoothing 演算法也有效的被運用到流固耦合的問題 中。 2001年，Gluck et al.【14】此篇以兩種code分別求解流力和固力，在流力方面 CFD code為FASTEST-3D，固力方面CSD code為ASE。再以MpCCI將兩種程式結

合在一起，藉此來模擬流固耦合問題。在網格方面，流體部分使Block-structured 網格，固力部分使用非結構性網格。此種計算方法可以適用在紊流和非線性移動 的暫態流固耦合問題。

2002 年，Slone et al.【5】 在處理流固耦合的問題時用 FV-UM 法來描述，其 中使用了 PHYSICA 軟體來分析。在計算方面，將流體和固體結合用同一套程式 求解。在流力部分，使用 FVM 來解 N-S 方程式，再利用壓力修正法來修正壓力、 速度。固力部分，利用 cell-vertex FVM 和 Newmark 演算法來求解方程式。此篇 分析的例子為在三維流場中置入懸臂樑，藉由流體的作用力，讓此物體產生擺 動。透過分析之位移、壓力、應力和流場來探討流固耦合的問體。此外本篇最終 目的是發展一個三維的計算流程來處理複雜 DFSI 問題，這項工作對於未來在 DFSI 的演算方面有更進一步的發展。

2004 年，Hubner et al.【15】利用 monolithic approach 來處理流固耦合問題， 在固力和流力部分都以有限元素法(space-time finite element)來離散並且加入耦合 條件，再透過 model equation 將三個結合在一起，因此在求解過程是將流體和固 體用一套程式一起求解。藉由此方法測試流固耦合問題中，發現在耦合情況下將 會影響到流場而發生改變。另一方面耦合強度也是一個重要依據，它將會影響其 耦合行為。 2007年，Lv et al. 【6】 在處理固力方程式中，發展了一個三維matrix-free 非結構性交錯網格有限體積法，這項方法已有效運用在2維和3維懸臂樑的問題。 在網格方面中，發現使用3-level多重網格法比單一網格的收斂速度快。流固耦合 方面，結合處理流力的軟體(TETRAKE)和應用於FSI界面immersed membrane method(IMM)也被證實有非常好的成效。

2008年，Xia and Lin【7】以一個cell-vortex非結構有限體積法來解固力方程 式，透過這個方法來模擬受到流體力影響的結構體。在計算方面為了使用較大的 時間步階，將導入dual-time scheme適用於此解法中。運用此方法結合計算流體力 學來測試耦合分析，其模擬結果與解析解和學者所發表之文獻做比較皆非常符

合。其中在分析的例子中，可發現到當擺動振幅越大週期也越長，另一個值得注 意的是結構體厚度和慣性力也是影響耦合的關鍵。 2009 年，陳俊佑 【16】 研究壓電無閥式微幫浦流固耦合分析，利用 ANSYS 來計算壓電致動器之振動位移並與 CFD 程式結合，藉此來探討流固耦合之現象。 在此模擬中，是將流體的壓力視為作用在壓電致動器之外力，以探討在不同操作 頻率下微幫浦之特性。由結果得知，當負載頻率愈高，愈不可忽略流體壓力對壓 電致動器之影響。 2010 年，陳信宏【17】研究雙壓電致動器幫浦之流固耦合分析，將 ANSYS 軟體與 CFD 程式結合，藉此來分析幫浦之特性。在耦合計算中，是將幫浦內流 體之壓力和剪力視為雙壓電致動器之外力限制條件。透過模擬結果發現，幫浦之 流量同時受到電壓、擺動振幅和振動頻率所影響。而欲增加流量最有效率的方法 為增加電壓的大小。 1.4 研究內容

本文研究內容是發展一非結構性網格有限體積法(Finite Volume Method)來計 算固體力學，並與計算流體力學互相結合來探討流固耦合之問題。 在固體力學方面，我們採用非結構性網格，並利用有限體積法離散結構方程 式，而控制體積中心則是取在網格頂點(cell-vertex)。另外在計算上本文將求解方 式分成直接解法和虛擬時間法兩類，並且分別利用龍格庫塔(Runge-Kutta)法或預 測修正法求解，進而比較不同方法結果的優缺點。除了探討結構位移之外，也將 觀察結構的自然頻率和不同阻尼下的運動情形。 在流體力學方面，採用 ALE 座標系之統御方程式，利用有限體積法將方程 式離散，再透過 PISO 修正法和空間守恆定理(Space Conservation Law)來滿足連續 方程式，而控制體積中心則是取在網格中心(cell-centered)。因此即可透過計算流 體力學方法，將流場計算之流體壓力和流體剪力施加於流體與固體的作用面上， 並搭配計算固體力學來分析流固耦合之問題，而在耦合測試的例子中分別為流經

彈性平板之管流和渦流致動之彈性平板。因此將藉由上述之方法來探討流體和平 板兩者之間的互相作用並觀察其平板的振動情形和流場變化。

第二章、數學模式

本章之數學模型分為固體運動模式和流體運動模式，因此在本章節分別描述

此兩者之運動方程式，以及不同的邊界條件。

2.1 固體運動模式

固體假設

在本研究中，探討之問題為二維，假設為線性彈性應變(Linear elastic strains)

結構運動方程式 結構體運動之統御方程式為 Cauchy’s equation: 2 2 ij i i j d b t x σ ρ∂ =∂ + ∂ ∂ (2.1) 其中

σ

_ij代表應力、d_i代表位移、

ρ

代表結構密度、 代表物體力 b_i 應力和應變關係 由虎克定律(Hook’s law)知應力(stress)和應變(strain)的關係，且對於等向性均 質(isotropic homogeneous)結構而言關係式如下: 1 0 1-(1 ) 1 0 (1 )(1 2 ) 1-1-2 0 0 2(1- ) xx xx yy yy xy xy E σ υ _ε υ υ υ _ε σ _{υ υ υ} υ ε σ _υ ⎡ ⎤ ⎛ ⎞ ⎛ ⎞ ⎢ ⎥ ⎜ ⎟ ⎜ ⎟ ⎢ ⎥ ⎜ ⎟ ⎜ ⎟ − _⎢ ⎜ ⎟₌ ⎜ ⎢ ⎜ ⎟ _{+ −} ⎜ ⎢ ⎥ ⎜ ⎟ ⎜ ⎟ ⎢ ⎥ ⎜ ⎟ ⎜_⎜ ⎟_⎟ ⎜ ⎟ _{⎢ ⎥ ⎝ ⎠} ⎝ ⎠ ⎣ ⎦ 　 ⎥ ⎟ ⎥ ⎟ (2.2)

其中 E 為楊氏系數(Young’s Modulus)、

υ

為蒲松比(Poisson’s ratio)。

而對於線性(linear)的應變和位移(displacement)的關係如下: x xx y yy y x xy d x x d y d d y x ε ε ε ⎛ ⎞ ⎡∂ ⎤ ∂ ⎛ ⎞ _{⎢ ⎥ ⎜ ⎟} ⎜ ⎟ _{∂ ⎜ ∂ ⎟} ⎢ ⎥ ⎛ ⎞ ⎜ ⎟ _{⎜ ∂ ⎟} ⎢ ∂ ⎥ ⎜ ⎟ ⎜ _{⎟ = = ⎜} ⎢ ⎥ ⎜ ⎟ ⎜ ⎟ _{∂ ⎜ ⎟ ⎜ ∂ ⎟} ⎢ ⎥ ⎜ ⎟ _{⎝ ⎠ ⎜ ⎟} ⎢ ∂ ∂ ⎥ _∂ ∂ ⎜ ⎟ ⎜ ⎟ _{⎢ ⎥} ⎜_⎜ + ⎟_⎟ ⎝ ⎠ ∂ ∂ ⎟ ∂ ∂ ⎣ ⎦ _{⎝ ⎠} ｘ ｙ 　　 　 ０ ｄ ０　　 　　 ｙ ｄ 　 ｙ ｘ (2.3)

而d_x和d_y分別為 X 和 Y 方向之位移 結構體邊界條件 邊界條件可分為牽引力邊界和位移邊界 (1)牽引力邊界(traction boundary): 在受力的邊界上須滿足力平衡之限制條件，如 下所示 T⋅ = tσ p (2.4) 其中 為邊界上之牽引力(traction)裡面包含額外的力t_p f ，T 為邊界外之單位法向 量，表示如下: 0 0 x y y x n T n ⎛ ⎞ = ⎜⎜_⎝ 　 　 n 　_⎠ 　 n 　 ⎟⎟ (2.5) (2)位移邊界(prescribed displacement): 在結構體上有已知位移限制條件，如下: d =dp (2.6) 其中

d

_p為受限制已知的位移。 2.2 流體運動模式 流體假設 在本研究中，探討之問題為二維，所假設之工作流體性如下: (1) 不可壓縮流(incompressible Flow) : 工作流體之密度為定值。 (2) 忽略重力項 (3) 非穩態流場(unsteady) (4) 層流(laminar flow) ALE 座標系之統御方程式 在本文研究中因邊界會移動，所以網格不是靜止不動，因此必須使用 ALE (arbitrary Lagrangian-Eulerian)座標系來求解移動邊界問題。以下為 ALE 座標系之 連續方程式和動量方程式:

(

V V_g

)

0 t ρ _ρ ∂ _{⎡ ⎤} + ∇ ⋅_⎣ − _⎦= ∂ (2.7) (2) 動量方程式(Momentum Equation)

( )

(

g

)

(

V V V V P V t ρ ρ ∂ ⎡ ⎤ + ∇ ⋅_⎣ − _⎦= −∇ + ∇ ⋅ ∇ ∂ μ

)

(2.8) 其中

ρ

為流體密度、

μ

為黏滯係數、∇P為壓力梯度、V 為流體速度、V_g 為網格 速度。 流場邊界條件 (1)出入口邊界條件 入口: 給定入口邊界的流量。 出口: 在出口方面為流量邊界，而邊界的出口速度計算均採用對流邊界條件 (convective boundary condition )。

(2)牆邊界條件 牆上均採用無滑移邊界條件(no-slip condition)，在彈性結構體之邊界也設為 牆的邊界條件。 2.3 流固耦合 固體條件 在結構體和流體接觸面之邊界條件(即牽引力邊界條件)為受流體壓力和流體 剪力之負載限制。 流體條件 計算流體力學中，在結構體和流體接觸面的速度等於結構體移動的速度。又 結構體位置會隨時間改變，因此在結構體和流體之接觸面也必須隨時間而改變位 置。而數學式表示如下: s V =V (2.9) s d =d (2.10) 其中V 代表流體速度、V_s代表結構體速度、d代表結構體和流體接觸面之位移、

s

d 代表結構體之位移。

因此可透過上面流體和固體之間互相作用的關係，來探討固體結構運動之情 形和流場現象。

第三章、固體力學數值方法

在固力計算方面，本文將控制體積取在網格頂點(cell-vertex)上如(圖 3.1)所 示，因此所算出位移和速度的值均儲存在頂點上。 3.1 有限體積法 將上一章的統御方程式(2.1)式忽略物體力(不考慮重力項)，另外加入一項阻 尼(damping)力，對於一般黏滯阻尼的結構體而言，其假設阻尼和結構速度成正比 關係，又阻尼運動方向和結構體運動方向相反，因此可將改寫成下列式子: 2 2 ij d cV t ρ∂ = ∇⋅σ − ∂ (3.1) 其中 c 表示阻尼係數，將上式取體積分後，可表式成如下之積分型式: 2 2 dd ijd cVd t _∀ ρ _∀ σ _∀ ∂ ∀ = ∇ ⋅ ∀ − ∀ ∂

∫∫∫

∫∫∫

∫∫∫

(3.2) 【一】 【二】 【三】 再透過高斯散度定理，將(3.2)式之【二】項由體積分轉成面積分，即可得到 下式: 2 2 ij S dd dS cV d t _∀ ρ σ _∀ ∂ _{∀ = ⋅ − ∀} ∂

∫∫∫

∫∫

∫∫∫

(3.3) 離散化 【一】 項 2 2 2 d dd t _∀ ρ ρ t ∂ _{∀ = Δ∀} ∂

∫∫∫

∂ 2 ∂ (3.4) 【二】 項 將(3.2)式之【二】項離散化後，表示如下: _ij ( _ij)_{f f} (3.5) f S dS dS σ ⋅ =

∑

σ ⋅

∫∫

其中

i

=1,2、

j

=1,2。因此可知

σ

₁₁ =

σ

_xx，

σ

₂₂ =

σ

_yy稱為正應力(normal stress); 12 21 xy

σ

=

σ

=

σ

稱為剪應力(shear stress)。下標 f 代表面上的值。而以頂點取控制

體積時又分為內部點和邊界點。當為內部點時，面上的值是取面上兩頂點之平均 值如圖 3.1(a)；當為邊界點時，需特別注意到部分面中心位於邊界上，所以面上 的值則是取邊界上的值如圖 3.1(b)。 將(2.2)式展開後可分別得到

σ

_xx、

σ

_yy、

σ

_xy，如下所示 (1 ) (1 ) ( ) y x xx y x yy y x xy d d x y d d x y d d y x σ λ υ λυ σ λυ λ υ σ μ ∂ ⎧ _{= −} ∂ ₊ ⎪ _{∂ ∂} ⎪ ⎪ ∂ ∂ ⎪ _{= + −} ⎨ _∂ ⎪ ⎪ _{∂ ∂} = + ⎪ ∂ ∂ ⎪⎩ ∂ (3.6) 其中 (1 )(1 2 ) E λ υ υ = + − 、 2(1 ) E μ υ = + 【三】 項 cVd cV (3.7) ∀ −

_∫∫∫

∀ = − Δ∀ 透過上述離散化後，將(3.3)式分成 X 和 Y 方向展開，即可得: 2 2 2 2 (1 ) ( ) ( ) (1 ) y y x x x x y f _{f f} y x y x y x y f _{f f} d d d d d S S t x y y x d d d d d S S t y x x y ρ λ υ λυ μ ρ μ λυ λ υ ⎧ _∂ ⎧_{⎪⎡ ∂ ∂ ⎤ ⎡ ∂ ∂ ⎤} ⎫_⎪ ⎪ Δ∀ = _⎨⎢ − + _⎥ ⋅ Δ +_⎢ + _⎥ ⋅ Δ _⎬− Δ∀ ∂ ∂ ∂ ∂ ∂ ⎪ ⎪⎣ ⎦ ⎣ ⎦ ⎪ ⎪ ⎩ ⎭ ⎨ ⎧ ⎫ ∂ ∂ ∂ ⎪ _{Δ∀ =} ⎪⎡ ∂ ₊ ⎤ _{⋅ Δ +}⎡ ∂ _{+ −} ⎤ _{⋅ Δ} ⎪_{− Δ∀} ⎨ ⎬ ⎪ _∂ ⎢ _{∂ ∂} ⎥ ⎢ _{∂ ∂} ⎥ ⎣ ⎦ ⎣ ⎦ ⎪ ⎪ ⎪ ⎩ ⎭ ⎩

∑

∑

cu cv (3.8) 3.2 內部梯度計算方式 以下將介紹位移梯度(displacement gradient)計算方式，在此方法中是將位移 梯度的值儲存於網格中心上，與結構體之位移和速度儲存在網格頂點而有所不 同。首先將梯度取體積分後，再透過高斯散度定理轉為面積分，如下所示: 1 1 1 _{f f} f S d dS

φ

φ

φ

φ

∀ ∇ = ∇ ∀ = = Δ∀

∫∫∫

Δ∀

∫∫

Δ∀

∑

dS (3.9) 因本文探討之問題為二維，所以可將上式改寫如下: 1 1 1 N j j j S dy y x A A φ _{φ φ} = ∂ = = ∂ Δ

∫∫

Δ

∑

Δ (3.10a)

1 1 1 N j j j S dx x y A A φ _φ = ∂ _{= − = − Δ} ∂ Δ

∫∫

Δ

∑

φ (3.10b) 其中N為網格面的個數如(圖 3.2)所示，而網格面上的

φ

值是將面的兩個頂點取平 均值即 1

(

_{1 2}

)

2 f f f φ = φ +φ ， f₁ 和 f2為f 面上的兩個頂點。 3.3 邊界梯度計算方式 當網格頂點為邊界點時，在計算上將使用到邊界面上的位移梯度，因此需另 外去計算。同樣的邊界位移梯度儲存的位置是在邊界面的中心，而計算方式如下: ∇ − ∇φ_b φ_P = ∇ ∇ ⋅( φ) δ _Pb (3.11) 其中下標 b 代表邊界、δ_Pb代表 P 到 b 距離向量如(圖 3.3)。 將位移梯度的梯度項利用(3.9)式相同之方法離散，即可得: ( ) 1 _{f f} 1 _{b b f} f f S S φ φ φ φ ≠ ⎛ ⎞ ∇ ∇ = ∇ = _⎜∇ + ∇ ∀

∑

_{∀ ⎝}

∑

f _⎠ b S _⎟ (3.12) 其中 f ≠b代表不包含邊界的其餘各面。將(3.12)式代入(3.11)式可得邊界上之位移 梯度: 1 1 1 P f f Pb f b b b Pb S S φ φ δ φ δ ≠ ∇ + ∇ ⋅ ∀ ∇ = − ⋅ ∀

∑

(3.13) 3.4 求解方法 本文在求解此方程式將分為兩大類方式，第一類為直接求解此方程式，第二 類則將方程式加入一個虛擬時間概念去求解。以下將分別介紹: 3.4.1 直接解法 在(3.8)式中因位移對空間的梯度項已由前面提到之計算方式離散化，因此方 程式即可視為一個二階常微分程式(ODE) ，另外在計算上本文將多加兩條方程 式，如下所示:

x y d u t d v t ∂ ⎧ ₌ ⎪⎪ ∂ ⎨∂ ⎪ ₌ ⎪ ∂ ⎩ (3.14) 其中u、v 分別代表 X 和 Y 方向速度，並且將(3.14)關係式代入(3.8)式，改寫後 再和(3.14)聯立。整理後可得: 1 (1 ) ( ) 1 ( ) (1 ) x y y y x x x y f _{f f} y y x x x y f _{f f} d u t d v t d d d d u c S S t x y y x d d d d v c S S t y x x y λ υ λυ μ u v ρ ρ μ λυ λ υ ρ ρ ∂ ⎧ ₌ ⎪ ∂ ⎪ ∂ ⎪ ₌ ⎪ ∂ ⎪⎪ _{⎧ ⎫} ∂ ∂ ⎨∂ ₌ ⎪⎡ ₋ ∂ ₊ ⎤ _{⋅ Δ +}⎡ ∂ ₊ ⎤ _{⋅ Δ} ⎪ ⎨ ⎬ ⎪_{∂ ⋅ Δ∀} ⎢ _{∂ ∂} ⎥ ⎢ _{∂ ∂} ⎥ ⎣ ⎦ ⎣ ⎦ ⎪ ⎪ ⎪ ⎩ ⎭ ⎪ ⎧_{⎡ ∂ ∂ ⎤ ⎡ ∂ ∂ ⎤} ⎫ ∂ ⎪ ₌ ⎪ _{+ ⋅ Δ + + − ⋅ Δ} ⎪ ⎨⎢ ⎥ ⎢ ⎥ ⎬ ⎪ ∂ ⋅ Δ∀ _⎪⎣ ∂ ∂ _{⎦ ⎣} ∂ ∂ _{⎦ ⎪} ⎪ ⎩ ⎭ ⎩

∑

∑

− − (3.15) 在求解過程中，將四個變數d_x、d_y、u、 令為一個矩陣型式如下 v U= x y d d u v ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ (3.16) 因此可將(3.15)式看成一個聯立常微分方程組，表示成下式: dU F U( ) dt = (3.17) (3.17)式即是一個一階 ODE 型式，因此可使用一般常見的 ODE 數值方法求 解。在此本文將使用兩種方法來求解，第一種為龍格-庫塔(Runge-Kutta)法、第二 種為預測修正(predictor-corrector)法。但由於每一網格點的 ODE 必需與相鄰格點 的 ODE 相互聯立在一起，因此需要利用疊代的方式將所有 ODE 耦合在一起。以 下將分別介紹求解之方法: Runge-Kutta Method 一般而言以四階 Runge-Kutta【18】最為廣泛使用。因此本文即使用此方法 求解，過程步驟如下:

令 ( ) 1 ( ) 2 1 ( ) 3 2 ( ) 4 3 ( ) ( ) 2 ( ) 2 ( ) n m n m n m n m K F U t K F U K t K F U K K F U t K ⎧ ₌ ⎪ Δ ⎪ _{= +} ⎪⎪ ⎨ _Δ ⎪ _{= +} ⎪ ⎪ = + Δ ⎪⎩ (3.18) ( +1) = +Δ₆ ( 1 +2 2 +2 3 + n m o t U U K K K K₄) 其中 代表時間步階而上標 o 代表上一個時間的值、上標 n 代表現在時間的值， 上標(m)及(m+1)代表在現在時間步階 n 疊代的指標，利用舊的疊代值 t Δ ( ) n m U 求得 1 ~ 4 K K ，可得新的疊代值Un m( +1)，如此疊代到收斂，便可往下個時階繼續計算。 Predictor-Corrector Method 在眾多預測修正解法當中，本文使用的是由 Hamming 所提出的方法【18】， 此方法是一個被廣泛而有效的運用在求解 ODE 的問題。以下為求解之步驟: Predictor: − − − = + Δ − + − (1) 4 4 1 2 3 (2 2 ) 3 n n n n n U U t F F F (3.19) Corrector: + − − − − = − + Δ + − ≥ ( 1) 1 1 3 3 ( ) 1 2 (9 ) ( 2 ), 1 8 8 n m n n n m n n U U U t F F F m ) (3.20) 其中 代表時間步階，至於上標 n 代表現在時間的值，n-1、n-2、n-3、n-4 代表 上一個、上二個、上三個、上四個時間值。上標(m)及(m+1)代表在現在時間步階 n 疊代的指標，利用一次的預測值 求得舊的疊代值 ，再經由修正步驟可 得新的疊代值 ，如此在修正步驟中疊代到收斂，即往下個時階進行計算。 t Δ (1) n U Fn m( ) + ( 1 n m U 此外，在預測修正解法中，前面三次的的計算必須使用龍格-庫塔法求得舊時間 值，第四次後即可使用預測修正法進行求解。 3.4.2 虛擬時間法 將(3.15)式等號左邊的微分項利用向後差分法(backward scheme)將其離散並

移置等號右邊，另外在等號左邊加入一項虛擬時間的微分項∂_∂U_τ ，而此虛擬時間

(virtual time)方法已被 Jameson 應用於求解機翼的超音速流場【19】，因此(3.15)

式即可改寫成下式: 1 (1 ) ( ) 1 ( ) n n o n x x x n n o y y y n n n n o n y y x x x y f _{f f} n n o y x d d d u t d d d v t d d d d u u u c S S t x y y x d d v v v t y x τ τ λ υ λυ μ u τ ρ ρ μ τ ρ ⎛ ⎞ ∂ _{= −} − ₊ ⎜ ⎟ ∂ _⎝ Δ _⎠ ⎛ ⎞ ∂ − = −⎜_⎜ ⎟_⎟+ ∂ _⎝ Δ _⎠ ⎧_{⎡ ∂ ⎤ ⎡ ∂ ⎤} ⎫ ⎛ ⎞ ∂ ∂ ∂ _{= −} − ₊ ⎪ _{− + ⋅ Δ + + ⋅ Δ} ⎪ ₋ ⎨ ⎬ ⎜ ⎟ ⎢ ⎥ ⎢ ⎥ ∂ _⎝ Δ _⎠ ⋅ Δ∀ _⎪⎩⎣ ∂ ∂ _{⎦ ⎣} ∂ ∂ _{⎦ ⎪⎭} ∂ ⎡ ⎤ ⎛ ⎞ ∂ ∂ _{= −} − _{+ +} ⎜ ⎟ ⎢ ⎥ ∂ _⎝ Δ _⎠ ⋅ Δ∀ _⎣ ∂ ∂ _⎦

∑

(1 ) n n y x x y f _{f f} d d cv S S x y λυ λ υ ρ ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ _{⎧ ⎫} ∂ ⎡ ∂ ⎤ ⎪ ⎪ ⎪ _{⎨ ⋅ Δ + + − ⋅ Δ ⎬ −} ⎢ ⎥ ⎪ _{⎪ ⎣ ∂ ∂ ⎦ ⎪} ⎩ ⎭ ⎩

∑

(3.21) 同樣的也可將(3.21)式聯立常微分方程組，表示成下式: ( ) n n dU G U dτ = (3.22) 因此(3.21)式同樣使用上面提到之求解 ODE 的方法求解，而(3.21)和(3.15)式 的差異在於等號左邊項，一個是對真實時間微分、另一個是對虛擬時間微分。在 求解此虛擬時間的 ODE 時，虛擬時間的步階即代表前述直接解法中的疊代次數。 當在虛擬時間求解過程中到達穩態時，即代表達到收斂，也就是滿足G U( n)= ，0 此時就是真實時間步階 n 的解。 由於對數值法穩定性的考慮，在直接解法中的時間步階通常受限於 CFL 條 件，也就是計算的庫倫數(Courant number)需要小於 1，而在本文的計算經驗中， 為有效求得穩定的解，此庫倫數通常需在 0.1 大小左右。而虛擬時間法的目的就 在於〝欺騙〞數值求解的方法，由於過程中以虛擬時間步進，因此只需對虛擬時 間步階做嚴格限制，使其庫倫數約 0.1，而實際步進項是吸收於虛擬時間 ODE 的 源項內，因此實際時間步階可以取較大值，通常可以是虛擬時間步階的 1 至 2 個 order 以上的大小，如此實際的庫倫數必可以大為增加，此點可從下一章的結果與

討論中明確看出。 下面介紹應用 Runge-Kutta 法及預測修正法於虛擬時間的解法中: Runge-Kutta Method 同樣以四階 Runge-Kutta 求解，步驟如下: 令 ( ) 1 ( ) 2 1 ( ) 3 2 ( ) 4 3 ( ) ( ) 2 ( ) 2 ( ) n m n m n m n m K G U K G U K K G U K K G U K τ τ τ ⎧ ₌ ⎪ Δ ⎪ _{= +} ⎪⎪ ⎨ _Δ ⎪ _{= +} ⎪ ⎪ = + Δ ⎪⎩ (3.23) ( 1) ( ) 1 2 3 ( 2 2 6 n m n m U + =U + Δτ K + K + K +K₄) (3.24) 其中Δτ 代表虛擬之時間步階而上標 n 代表現在時間、(m)代表在虛擬時間下的上 一個時間值。因為在每個真實的時階下都必須要滿足G U( n)= ，所以在求解過程0 中當Gn m( )趨近於零時，表示K1 ~ K4也趨近零，因此即可求得 ( ) n m U 值逼近Un m( + )1 值，即滿足G U( n)=0，此時便可繼續往下一個真實時階計算。 Predictor-Corrector Method 同樣的將解法改成預測修正法時也須要在每個真實時階下解到滿足方程式 ( n) G U = 0 ，因此求解之步驟如下: Predictor: ( 1) ( 3) 4 ( ) ( 1) ( 2) (2 2 ) 3 n m n m n m n m n m p U + =U − + Δτ G −G − + G − (3.25) Corrector: ( 1) 1(9 ( ) ( 2)) 3 ( ( 1) 2 ( ) ( 1)) 8 8 n m n m n m n m n m n m c p U + = U −U − + Δτ G + + G −G − (3.26) 其中Δτ 代表虛擬之時間步階而上標 n 代表現在時間、p 代表預測之值、c 代表修 正後之值，至於上標(m)代表在虛擬時間下的上一個時間值、(m-1)代表在上二個 時間值、(m-2)代表在上三個時間值、(m-3)代表在上四個時間值。在求解過程中， 當Un m( )趨近於Un m( + )1 值時，代表收斂即滿足G U( n)=0_，即可往下一個真實時階 計算。

疊代收斂條件 在直接解法和虛擬時間法中都必須在每個時階下疊代至收斂，其收斂與否以 位移來判斷，數學式表示如下: ( 1) ( ) 6 ( ) 1 10 n m n m N n m i d d d + − = − ≤

∑

(3.27) 其中 N 代表網格頂點個數、n 代表現在時間、(m)代表在疊代中的上一個疊代值。 3.5 邊界上之計算 在結構體中相鄰邊界的控制體積如(圖 3.4)所示，當此邊界有外力之限制條件 時，此時處理方式跟內部點不一樣，因此必須額外處理。將邊界上控制體積的面 分為邊界面 f_b和非邊界面 f_i，在 f_b面上直接給定此面所負載之外力

F

，而在 f_i面 上給定方式和內部點一樣為面上所需的值。又外力給定方式可分為分佈力和集中 受力，因此在邊界上之計算可將(3.8)式改寫成兩種受力方式的式子。當外力負載 為分佈力 w 時，可將(3.8)式改寫成下式:

[

]

2 2 2 2 (1 ) ( ) ( ) ( ) (1 ) ( ) b i b b i b y y x x x x y x _f f y x y x y x y _f f d d d d d S S w S t x y y x d d d d d S S w t y x x y ρ λ υ λυ μ ρ μ λυ λ υ ≠ ≠ ⎧ ∂ _{Δ∀ =} ⎧⎪⎡ ₋ ∂ ₊ ∂ ⎤_{⋅ Δ +}⎡ ∂ ₊∂ ⎤_{⋅ Δ} ⎪⎫_{+ ⋅ −Δ − Δ∀} ⎪ _∂ ⎨⎢ _{∂ ∂} ⎥ ⎢ _{∂ ∂} ⎥ ⎬ ⎪⎣ ⎦ ⎣ ⎦ ⎪ ⎪ ⎩ ⎭ ⎨ ∂ _{Δ∀ =} ⎧⎪⎡ ∂ ₊∂ ⎤_{⋅ Δ +}⎡ ∂ _{+ −} ∂ ⎤_{⋅ Δ} ⎫_{⎪ ⎡+ ⋅ −Δ ⎤ − Δ∀} ⎨⎢ ⎥ ⎢ ⎥ ⎬ ⎣ ⎦ ∂ _⎪⎩⎣ ∂ ∂ _{⎦ ⎣} ∂ ∂ _{⎦ ⎪⎭}

∑

∑

⎪ ⎪ ⎩ y cu S cv (3.28) 而當外力負載為集中受力

F

時，即可將(3.8)式改寫如下: 2 2 2 2 (1 ) ( ) ( ) (1 ) i b i b y y x x x x y f y x y x y x y f d d d d d S S t x y y x d d d d d S S t y x x y ρ λ υ λυ μ ρ μ λυ λ υ ≠ ≠ ⎧ ∂ _{Δ∀ =} ⎧⎪⎡ ₋ ∂ ₊ ∂ ⎤_{⋅ Δ +}⎡ ∂ ₊∂ ⎤_{⋅ Δ} ⎪⎫_{+ − Δ∀} ⎪ _∂ ⎨⎢ _{∂ ∂} ⎥ ⎢ _{∂ ∂} ⎥ ⎬ ⎪⎣ ⎦ ⎣ ⎦ ⎪ ⎪ ⎩ ⎭ ⎨ ∂ ⎧⎪⎡ ∂ ∂ ⎤ ⎡ ∂ ∂ ⎤ ⎫⎪ ⎪ _{Δ∀ = + ⋅ Δ + + − ⋅ Δ + − Δ∀} ⎨⎢ ⎥ ⎢ ⎥ ⎬ ⎪ ∂ _⎪⎩⎣ ∂ ∂ _{⎦ ⎣} ∂ ∂ _{⎦ ⎪⎭} ⎩

∑

∑

x y F cu F cv (3.29)

第四章、流體力學數值方法

在流力計算方面，本文將控制體積取在網格中心(cell-centered)上如(圖 4.1)所 示，因此所計算出速度和壓力的值均儲存在網格中心上。 4.1 空間守恆定理 空間守恆定理介紹 在探討移動網格的問題時，有一個定理專門用來判斷網格體積隨時間變化是 否守恆，稱為空間守恆定理(Space Conservation Law)，數學表示如下:

0 g S d d V dS dt

∫∫∫

_∀ ∀ +

∫∫

⋅ = (4.1)

而此定理在 Demirdzic and Peric【20】發表的文獻中已被應用於有限體積法中。 另外在處理移動網格問題時，需用到每個時階下的網格體積改變量，因此在計算 上本文直接用新舊網格各面的節點位置來做計算。如(圖 4.2)，藉由找出網格一個 面移動前後的四個節點，利用這四個點來算出此面移動造成的體積改變量。 空間守恆定理 在每個時階中，新和舊體積的差值會等於在控制體積中各面所掃過的體積改 變量相加，如下所示: f o f t t Δ∀ ∀ − ∀ ₌ Δ Δ

∑

(4.2) 其中Δ∀_f 為一個面掃過之體積， f =e w n s, , , ，如(圖 4.3)所示。再透過(4.1)式可 知網格移動的速度等於每個時階中體積改變量，因此對於每個面上之關係式可表 示如下:

(

Vg S

)

_f f f t Δ∀ ⋅ = = ∀ Δ (4.3) 因此流過移動網格某一面的質量流率 r 可寫成下式: f

m

r

(

)

f f g _f f f f f f m =ρ V −V ⋅S =ρ V ⋅S −ρ ∀_f (4.4) 接著也必須讓連續方程式滿足空間守恆定理，將(2.7)式取體積後再將其離散 後改寫成下式: 0 o r f f m t ρ∀ − ∀ + Δ

∑

= (4.5) 將(4.4)式代入上式可得: f f f 0 (4.6) f V S ρ ⋅ =

∑

而(4.6)式即下列不可壓縮流之連續方程式的差分式子: ∇ ⋅

( )

ρV =0 (4.7) 因此本文 ALE 座標系的連續方程式(2.7)式，可簡化為和網格速度無關的(4.6) 式，再配合(2.8)式的動量方程式即可描述本文的流體運動模式。 4.2 有限體積法 由第二章中的動量方程式(2.8 式)可以得到下列通式:

( )

(

V Vg

)

(

)

q t ρφ _{ρ φ μ φ} ∂ _{⎡ ⎤} + ∇ ⋅_⎣ − _⎦= ∇ ⋅ ∇ + ∂ (4.8) 將上式積分後，可表示如下之積分型式: d

(

V Vg

)

d

(

)

d t _∀ ρφ _∀ ρ φ _∀ μ φ _∀ ∂ _{∀ + ∇ ⋅⎡ − ⎤ ∀ = ∇ ⋅ ∇ ∀ + ∀} ⎣ ⎦ ∂

∫∫∫

∫∫∫

∫∫∫

∫∫∫

qd P (4.9) 等式左邊第一項為非穩態項(Unsteady term)，第二項為對流項(Convection term)； 等式右邊第一項為擴散項(Diffusion term)，第二項為源項(Source term)。源項的來

源是壓力項，故q= −∇ 。又上面所提到之

φ

為速度分量。 再透過高斯散度定理，將(4.9)式對流項和擴散項由體積分轉換成面積分，即 可得到下式:

(

g

)

S S d V V dS dS t _∀ ρφ ρ φ μ φ _∀ ∂ _{∀ + − ⋅ = ∇ ⋅ + ∀} ∂

∫∫∫

∫∫

∫∫

∫∫∫

qd (4.10) 離散化 以下將介紹動量方程式各項的離散方法

非穩態項(Unsteady term) 對於有移動邊界的不可壓縮流，其網格的體積會隨時間而改變。

( )

(

o

)

d t t ρφ ρ _{φ φ} ∀ ∂ ⋅ Δ∀ ∀ = − ∂ Δ

∫∫∫

(4.11)

φ

是這一個時間的值、

φ

o 是上一個時間的值，因此在計算上可將

φ

0 移進源項。 對流項(Convection term) 將(4.10)式的對流項，已由體積分轉為面積分後，再將其離散化，可得:

(

g

)

( f

(

g

)

_f f) f f f C f f S V V dS V V S m F ρ − φ⋅ =

∑

ρ − ⋅ φ =

∑

φ =

∫∫

f f

∑

(4.12) 其中

S

f 代表面上之值、m_f 代表流經面上的質量流率、

F

_fC代表面上的對流通量、 下標 f 代表面上的值。而對於面上的值，可透過相鄰網格以線性內差的方式求 得。流體相對於面上的質量流率定義如下: r

(

)

f f g _f m =ρ V −V ⋅Sf (4.13) 在 ALE 座標系中，因動量方程式之對流項加進了網格速度的變化，網格移 動的速度會直接影響面上的質量流率。所以在質量流率的計算上必須配合空間守 恆定律，因此透過上小節推導可將質量流率改寫成如(4.4)式的計算方式，如下: r

(

)

f f g _f f f f f f m =ρ V −V ⋅S =ρ V ⋅S −ρ ∀_f (4.14) 本文在求解對流項所採用的方式為一階上風法(UD)和中央差分法(CD)的混 合式，因此對流通量可整理如下: FfC = FfUD + γ

(

FfCD −Ff

)

UD (4.15) 其中

γ

值介於 0~1 之間，當

γ

=0 時為一階上風法、

γ

=1 時為中央差分法。上式右 手邊第二項為前一個時間的值，故將其移至源項。 在一階上風法中

φ

_f 的值取決於上游，當m_f>0 時，則

φ

_f =

φ

_p；當m_f <0，則 f

φ

=

φ

_nb。

在中央差分法裡，

φ

_f 的值為相鄰網格的內差，表示如下: φ_f = −

(

1 w

)

φ_p+wφ_nb (4.16) 其中 為權重係數，下標w p代表主格點的值，下標nb代表鄰近格點的值。 故可將(4.15)式整理成下列型式: F_fC =max

(

m_f, 0

)

φ_p−max

(

−m_f, 0

)

φ_nb +γ

{

m_f ⎡_⎣

(

1−w

)

φ_po +wφ_nbo ⎤_⎦−max

(

m_f, 0

)

φ_po+max

(

−m_f, 0

)

φ_no_b

}

(4.17) 其中上標 o 代表上一個時間的已知數，因此在計算上將已知數移進源項。 擴散項(Diffusion term) 將(4.10)式的擴散項，已由體積分轉為面積分後，將其離散可得: D

( )

D f f f f f f S F =

∫∫

μ φ∇ ⋅dS ≈

∑

μ ∇φ ⋅S =

∑

F (4.18) 其中

μ

_f 代表面上的黏滯係數、 D f

F

代表面上的擴散通量。針對非結構性網格而 言，在擴散項的部分則利用 Over-Relax 的方法來近似其原來的面法向量。首先將 表示如下: f

S

Sf = +d

(

Sf −d

)

(4.19) 上式等號右手邊第一項為正交向量、第二項為非正交向量。而正交向量 的定義 如下: d 2 f Pnb Pnb f S d S δ δ = ⋅ (4.20) 其中

δ

_Pnb定義如(圖 4.4)所示。將(4.19)式和(4.20)式代入(4.18)整理可得:

(

)

(

2 f f D f nb P f f Pnb f S

)

o f F S d S μ φ φ μ φ δ = − + ∇ ⋅ ⋅ − (4.21) 其中上標 o 代表上一個時間的已知數。且上式右手邊第一項為正交項、第二項為

非正交項，在計算上將已知數非正交項移進源項。 源項(Source term) 源項的體積分可表示為: qd q (4.22) ∀ ∀ = ⋅ Δ∀

∫∫∫

而源項的來源是壓力項。 壓力項的面積分知化簡近似: f (4.23) f S P dS⋅ ≈

∑

P dS⋅ = ∇ ⋅ Δ∀P

∫∫

線性代數式的整理 透過前面的推導，可將非穩態項、對流項、擴散項和源項合併後，即可得到 下列的線性代數式: (4.24) P p nb nb Aφ =

∑

A φ +Q 其中AP Anb t ρΔ∀ = + Δ

∑

(4.25)

(

)

2 max , 0 f f nb f Pnb f S A S μ δ = + − ⋅ m (4.26)

{

(

CD UD

)

o o

(

)

}

o o f f f f f f Q F F S d P t γ μ φ ρ Δ∀φ ⎡ ⎤ = −_⎣ − _⎦ + ∇ − − ∇ ⋅ Δ∀ + Δ

∑

(4.27) 4.3PISO 演算法 在前面介紹了動量方程式的求解方式，但是在物理上其實未滿足連續方程 式。因此為了滿足連續方程式，必須透過壓力和速度的關係式來加以修正。到目 前為止，已發展出幾種壓力和速度的耦合方式。因本文流場為非穩態，所以採用 PISO 演算法來修正動量方程式所得到的壓力和速度，如此一來即可滿足連續方 程式。使用 PISO 演算法的好處是，此解法不需經由反覆的疊代，因此可以節省 很多的計算時間。PISO 分為一個的預測步驟和兩次修正步驟。 預測步驟 將(4.24)式中 Q 裡的壓力項提出表式成下式:

(

*

)

P P nb nb P A V∗ =

∑

A V∗ + S− ∇P Δ∀ (4.28) 利用現有之壓力 * p

P

即可得到速度

V

*p，表示如下: * ( ) * P P P P V H V P A ∗ ₌ ∗ _{− ⎜}⎛Δ∀⎞ ⎟ ⎝ ⎠ ∇ P (4.29) 其中 * nb nb nb p p A V S H A ∗ ₊ =

∑

， 為不包含壓力項的S Q。 類似的，控容面上之速度和壓力關係式可透過(4.29)式表式成: * f f P f V H P A ∗ ₌ ∗ ₋⎛Δ∀⎞ _∇ ⎜ ⎟ ⎝ ⎠ f (4.30) 其中 f f P f H V P A ∗ _{= +}⎛Δ∀⎞ _∇ ⎜ ⎟ ⎝ ⎠ f (4.31) 將(4.31)式代入(4.30)式，因此控容面上的速度和壓力關係式如下: * f f f P _f P V V P P A A ∗ _{= +}⎛Δ∀⎞ _{∇ −}⎛Δ∀⎞ _∇ ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ ⎝ ⎠f f (4.32) 其中上標”─”表示主格點p和控容面 f 相鄰之格點nb內差所得，如下: ∇ =P_f w_p∇P_nb+ −(1 w_p)∇P_p (4.33) V_f =w V_{p nb}+ −(1 w V_p) _p (4.34) 其中w_p為權重因子(weighting factor) 而 P _f A ⎛Δ∀⎞ ⎜ ⎟ ⎝ ⎠ 的值為主格點和相鄰格點取平均，如下: 1 2 P _f P _P P A A A ⎡ ⎤ ⎛Δ∀⎞ ₌ ⎛Δ∀⎞ ₊⎛Δ∀⎞ ⎢ ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎢ ⎥ ⎝ ⎠ _⎣⎝ ⎠ ⎝ ⎠nb_⎦ ⎥ (4.35) 由前面介紹的空間守恆定理知，可將 ALE 座標的連續方程式簡化為和網格 速度無關之連續方程式，因此面上的質量流率可改寫為:

mf ρfVf f ∗ ₌ ∗ _{⋅ S}

(

)

2 _f f _{f f} f _{nb P f Pnb} P f Pnb f S V S P P P A S ρ ρ δ ⎛ ⎞ ⎛_{Δ∀ ⎜}⎞ _{⎟⎡ δ ⎤} = ⋅ − _{⎜ ⎟ ⎣} − − ∇ ⋅ _⎦ ⎜ ⋅ ⎟ ⎝ _{⎠ ⎝ ⎠} (4.36) 第一次修正步驟 在前一個步驟已求得速度

V

_f∗和壓力 仍然尚未滿足連續方程式，因此接下 來必須做第一次的修正。在此假設第一次修正後的速度為 *

P

** V 、壓力為 ，故第 一次修正中，速度和壓力的關係如下: **

P

* ( ) ** P P P P V H V P A ∗∗ ₌ ∗ ₋⎛Δ∀⎞ _∇ ⎜ ⎟ ⎝ ⎠ P (4.37) 其中V∗∗ =V∗+V′；P**=P*+P'，而上標” ”即為所要修正的修正量。 ' 將(4.37)式與(4.29)式相減，可得到速度修正量和壓力修正量的關係式: P P P P V A ⎛Δ∀⎞ ′ = −_{⎜ ⎟} ∇ ⎝ ⎠ P ′ (4.38) 上述的方程式是在主格點上推出，亦其可以使用在每個主控面上。 f f P _f V A ⎛Δ∀⎞ ′ = −_{⎜ ⎟} ∇ ⎝ ⎠ P ′ (4.39) 因此第一次修正後的流量值可以表式成下列:

(

)

(

)

2 ( ) f f f f f f f f f f f f f P f P f f f f nb P f f f P _f Pnb f P _f m m m m V S m P P S d A A S m P P P A S A ρ ρ ρ ρ ρ δ ∗∗ ∗ ∗ ∗ ∗ ′ ′ = + = + ⋅ ⎛Δ∀⎞ _′ ⎛Δ∀⎞ _′ − _{⎜ ⎟} ∇ ⋅ − _{⎜ ⎟} ∇ ⋅ − ⎝ ⎠ ⎝ ⎠ ⎛Δ∀⎞ _{′ ′} ⎛Δ∀⎞ _′ − _{⎜ ⎟} − − _{⎜ ⎟} ∇ ⋅ ⎝ ⎠ ⎝ ⎠ ⋅ S −d 　 　 = d　 　 = (4.40) 為了滿足質量守恆，因此對於所有面上流量的總和應該為 0，表示如下: (4.41) + = f f f f f f m∗∗ = m∗ m′

∑

∑

∑

0

最後將(4.40)式代入(4.41)式，可整理成線性代數式: 1 1 1 2 P P nb nb P P f A P′=

∑

A P ′+S +S (4.42) 其中上標”1”代表第一次修正，而其它定義如下所示: 2 f nb f P f Pnb f S A A ρ δ ⎛Δ∀⎞ = _{⎜ ⎟} ⋅ ⎝ ⎠ S b (4.43) P n f A =

∑

A (4.44) (4.45) 1 1 P f f S = −

∑

m∗ 2

(

2 P f f f f P _f S P A ρ ⎛Δ∀⎞ ′ = _{⎜ ⎟} ∇ ⋅ ⎝ ⎠

∑

S −d

)

(4.46) 第二次修正步驟 在第一次的修正步驟中，得到了第一次修正過後的速度V 和壓力 。在第 二次修正步驟中，如同第一次修正步驟的方法，因此速度和壓力的關係式如下: ** **

P

*** ** ** *** ( ) P P P P P V P A ⎛Δ∀⎞ = −_{⎜ ⎟} ∇ ⎝ ⎠ V H (4.47) 將(4.47)與(4.37)相減，可得到速度修正量和壓力修正量的關係式: P P( ) P P P V H (4.48) V P A ⎛Δ∀⎞ ′′_{= ′ ′ −⎜ ⎟ ∇} ′′ ⎝ ⎠ P 其中V′′ =V∗∗∗−V∗∗、V′ ≡V∗∗−V∗、∇ = ∇P′′ P***− ∇P**_，另外HP′ =H**P−H* 又 f f mf ∗∗∗ ∗∗ ′′ = − m m ，因此質量流率的修正量m′′_f 可表示如下:

(

)

(

)

2 = P = nb nb f f f f f f f P f P nb nb f f f nb P f f f f f P _f Pnb f P _f P A V m S S A A A V S P P P S d S A S A A ρ ρ ρ ρ ρ δ ′ ⎛Δ∀⎞ ′′ _{− ⎜ ⎟ ∇ ⋅′′ + ⋅} ⎝ ⎠ ′ ⎛Δ∀⎞ _{′′ ′′} ⎛Δ∀⎞ _′′ − _{⎜ ⎟} − − _{⎜ ⎟} ∇ ⋅ − + ⋅ ⎝ ⎠ ⎝ ⎠

∑

∑

⋅ 　 　 　 (3.68) (4.49)