流固耦合模擬分析之探討VIP

MSC.Dytran 流固耦合模擬分析之探討摘要：流場會驅(qū)動固體；而固體的運動也可能驅(qū)動流體，甚至引起流場振動。這稱為流固互制或流固耦合。MSC.Dytran軟件的流固耦合功能包含一般耦合、任意偶合，且采用拉格朗日法與歐拉法分別描述固體與流體的運動。Dytran中流固耦合計算種類拉格朗日的元素節(jié)點依附在材料上，節(jié)點隨著材料質(zhì)點作運動，故各物理量也作用在節(jié)點上隨材料流動而變化。相反，除任意耦合外，歐拉元素網(wǎng)格與節(jié)點不隨時間而變，其物理量雖也作用在歐拉元素節(jié)點上，但對于通過歐拉元素面的各時間的質(zhì)量、動量與能量的進與出，加之模擬，即模擬元素面的材料流，而不模擬各材料質(zhì)點的時間歷程。因為對一般固體材料，要模擬各材料質(zhì)點的時間歷程，因此大多用拉格朗日法。而對于流體不需要模擬材料質(zhì)點的時間歷程，故采用歐拉法，MSC.Dytran的歐拉法需用三維的計算域、三維的體元素與DMAT通用材料。此外，歐拉法容許一個元素內(nèi)含有兩種以上的材料，這就是模擬計算材料流的擴散與混合行為歐拉法模擬計算材料流的擴散與混合行為。欲推廣應(yīng)用型的計算軟件，需有充分的應(yīng)用范例。關(guān)于流固耦合的模擬，除MSC.Dytran中文范例手冊（1999b）與Example Problem Manual(2001a)論述到一般耦合與任意耦合的應(yīng)用范例之外，本文進一步探討其他的應(yīng)用。而搭配的MSC.Patran軟件除有Results的后處理工具外；該軟件的Insight工具，能透明地看到體元素所構(gòu)成的歐拉域內(nèi)部，因此，更需用Insight，去展示歐拉域內(nèi)部的流固耦合計算結(jié)果。一、 前言MSC.Dytran之流固耦合計算功能大致上包括一般耦合(general coupling)和任意耦合(ALE, Arbitrary Lagrangian-Eulerian)兩大類。欲推動一項泛用型計算軟件被廣泛地應(yīng)用，須有可供參考的文件及豐富的應(yīng)用案例。而MSC.Dytran之中范例手冊包括一個一般耦合和兩個任意耦合的應(yīng)用范例。MSC.Dytran 之Example Problem Manual (2001)也包括兩個一般耦合、兩個任意耦合的范例。本文之要旨是經(jīng)由案例，探討MSC.Dytran軟件對流固耦合（歐拉域）與固體（拉格朗日域）的一般耦合與任意耦合的模擬計算，以及MSC.Dytran軟件對流固耦合計算之前后處理與模型結(jié)果展示。本文所敘述的流固耦合計算應(yīng)用場合如下： 水下之固體物的告訴移動 造波板與之波浪水槽 海面上之高速物體撞擊混凝土墻 上游側(cè)與下游側(cè)水深相等情況之潛堰 上游側(cè)與下游側(cè)水深不等情況之固定堰 上游側(cè)與下游側(cè)水深不等情況之階梯式渠道 固定開度之水利閘門 隧洞內(nèi)之氣爆壓作用在鋼板上 隧洞內(nèi)之水流推動固體物二、 對運動的數(shù)學描述法拉格朗日法與歐拉法是對運動現(xiàn)象的兩類不同的數(shù)學描述，可說是分別對材料質(zhì)點流與空間流之描述。拉格朗日法與歐拉法之元素網(wǎng)格可在同一計算模型內(nèi)，但拉格朗日法的元素與歐拉法的元素分別擁有節(jié)點，只采用介面(interface)，稱為耦合面，才能將兩者連結(jié)在一起；否則，縱使兩者在空間內(nèi)相互重疊，也彼此不相干，即忽視對方之存在。2.1 拉格朗日法對固定的坐標系而言，拉格朗日元素的節(jié)點可相對地運動。因節(jié)點系附在材料上，故材料連續(xù)體之節(jié)點系一起隨著材料質(zhì)點流而運動。各拉格朗日元素的質(zhì)量是不變量(invariant)，但其元素體積可隨時間而改變。此外，速度、壓力強度或質(zhì)量密度等物理量系作用在拉格朗日元素的節(jié)點上，因此，各物理量系隨著材料流(material flow)而改變。因?qū)腆w材料之行為， 較須追蹤各材料質(zhì)點之時間歷程，故適宜采用拉格朗日法。拉格朗日法也適宜用以分析材料破壞(failure)或應(yīng)變硬化(strain hardening)問題。2.2歐拉法除任意耦合(ALE)外，所有歐拉元素的網(wǎng)格與節(jié)點均保持固定，不隨著時間或其他運動或變形。換言之，歐拉元素不隨時間而變。各時間之速度、 壓力強度或密度等物理量也是作用在歐拉節(jié)點上，系計算經(jīng)歐拉元素面的各時間的質(zhì)量、動量與能量等的進與出之量，而不追蹤各材料質(zhì)點之時間歷程。 因?qū)α黧w較不須追蹤各材料質(zhì)點之時間歷程，故一般系采用歐拉法。MSC.Dytran之歐拉法必定使用三維的計算域及三維的體元素(solid elements)，且限于通用材料 DMAT。歐拉法之特點是：須采用較大之計算域，計算結(jié)果才不會受到計算域的邊界之影響.歐拉法的特點和要求由于歐拉法系仿真經(jīng)過歐拉元素面的材料流；且在一個元素之內(nèi)，容許兩種以上之材料，因此，應(yīng)用歐拉法可模擬計算空氣或水等材料流之擴散與混合現(xiàn)象。歐拉法計算混合現(xiàn)象三、 流固耦合計算法之種類當固體影響流體之后，被改變后之流體反過來影響固體；另一方向，亦然，就是流固耦合。數(shù)值模擬計算是探討分析流固耦合問題的主要方法之一。由于流場動壓變化所產(chǎn)生之流場特性、固體之幾何形狀、振幅與振頻等互有關(guān)系，故流固耦合本質(zhì)相當復雜?？v使是單方向之流固耦合分析，對影響固體振動之水動力的附加質(zhì)量也大多是估計。同時，也不易準確地預測及量化固體系統(tǒng)內(nèi)部之阻尼與流體吸收動能之效應(yīng)。因此，相關(guān)之數(shù)值模擬計算的難度相當高。關(guān)于流固耦合之計算分析法可分為兩大類如下：3.1 單方向之流固耦合分析此類分析為簡化之流固耦合計算。即考慮固體單方向影響流體，但不考慮流體反過來影響固體；或反之。3.2 雙方向之流固互制分析此類分析包括流固耦合的順序分析法與完全的流固耦合分析法兩種如下：3.2.1 流固耦合之順序分析法此法系先分別計算流體與固體領(lǐng)域，每完成其中的一個領(lǐng)域之計算后，將計算結(jié)果作為另一個領(lǐng)域之荷載(loading)或邊界條件，進行另一個領(lǐng)域的計算；反過來，也是相同的作法。當計算軟件系由固力計算程序與流力計算軟程序結(jié)合而成時，就采用此種分析法，MSC.Dytran 就是一例。 有些文獻作較嚴格的定義，稱此法并非耦合(uncoupled)分析，僅能稱為流固互制分析。此法并非耦合分析，僅能稱為流固互制分析3.2.2 完全的流固耦合分析法將所有之流固耦合相關(guān)的參數(shù)、邊界與荷載等均融入流場與固體所共享之控制方程式組內(nèi)，再采用數(shù)值計算法，求解耦合的(coupled)聯(lián)立方程式，故作較嚴格定義的文獻認為此法才是真正的流固耦合分析。此種方式雖最完整，但難度也最高，故使用者最少。真正的流固耦合分析四、流固耦合計算功能MSC.Dytran 之流固耦合計算就是拉格朗日域(固體)與歐拉域(流場)的耦合計算，輸入檔內(nèi)容直接相關(guān)者系 COUPLE卡或 COUPLE1卡。在計算過程內(nèi)，拉格朗日域與歐拉域分別進行計算，每完成一個領(lǐng)域之計算后，再計算另一個領(lǐng)域；即前一領(lǐng)域之計算結(jié)果作為另一個領(lǐng)域計算所需之荷載或邊界條件。至于 MSC.Dytran 的流固耦合計算，則分為兩種如下：4.1 一般耦合通常是采用拉格朗日法模擬固體，以及采用歐拉法模擬流體。至于一般耦合(general coupling)，大多是拉格朗日的固體在歐拉的流場范圍內(nèi)運動，即拉格朗日域驅(qū)動歐拉域一般耦合應(yīng)用在Lagrange的固體在Euler的流場范圍內(nèi)運動。；流場雖有速度，但代表流場的歐拉格網(wǎng)系固定及不受拉格朗日的固體之影響。換言之，在流固耦合過程內(nèi)，歐拉格網(wǎng)不移動，也不變形。一般耦合計算之前處理大多用封閉之假殼(dummy shell) 耦合面隔開拉格朗日域與歐拉域；于起始計算之時，拉格朗日的固體至少有一微小量(譬如，0.001 m)重疊在歐拉域的范圍內(nèi)，且固體、流體或兩者有運動，才能啟動流固耦合之計算器制。若拉格朗日域與歐拉域毫無重疊，則無法起動流固耦合計算。當然，拉格朗日的固體可完全位于歐拉域內(nèi)，不因拉格朗日的固體運動而使歐拉格網(wǎng)移動或變形。此種流固耦合最適用在固體驅(qū)動流體現(xiàn)象之模擬，譬如，水上飛機降落在水面上之行為、水面下之物體的運動所引致之流體動力行為、隧道內(nèi)之高速車輛引致的氣動力行為、 管道內(nèi)的固體活塞推動氣體等。此外，一般耦合也可應(yīng)用在具有不規(guī)則的固定固體邊界之流場模擬計算，即將不動的固體邊界與流場之關(guān)聯(lián)視為流固耦合現(xiàn)象。一般耦合須在拉格朗日的幾何體之外露面(幾何面段或元素表面，經(jīng)由SURFACE 卡建立)，定義封閉的假殼元素群，作為流固耦合面，形成閉合之體積，完全隔開拉格朗日域與歐拉域；且對假殼元素，須進行 Equivalence 與 Verify 法線方向之前處理。耦合面雖為虛擬，但也須輸入物理性質(zhì)(Properties)為 2D/Summy Shell。然后，對耦合面的全部假殼元素， 使用 MSC.Patran 軟件的 Coupling 等，就可完成相關(guān)的前處理。在計算過程內(nèi)，耦合面是歐拉流場格網(wǎng)的邊界，即假殼元素之耦合面與仿真流體之體元素相接觸(Contact)， MSC.Dytran 先計算歐拉流體施加在耦合面上之荷載，然后，耦合面使拉格朗日格網(wǎng)產(chǎn)生應(yīng)力與變形。至于輸入檔相關(guān)的 Couple 卡上之參數(shù) COVER，則定義被耦合面所包裹者（即不含歐拉領(lǐng)域者）是屬內(nèi)部(INSIDE)，還是外部(OUTSIDE)。譬如，石頭掉入水中之模擬，其 COVER 是 INSIDE。汽車的安全氣囊的模擬之COVER 是 OUTSIDE。4.2 任意耦合任意耦合主要是用來仿真歐拉域驅(qū)動拉格朗日域，在該耦合過程內(nèi)，于計算模型承受來自歐拉域的荷載之后，拉格朗日的元素與歐拉的體元素之網(wǎng)格均可能變形或移位。任意耦合最適宜用來仿真流體驅(qū)動固體，譬如，容器內(nèi)之氣體爆炸、空中之鳥體撞擊飛機殼體、氣爆壓作用在殼形固體物上、 管路內(nèi)的流體推動固體、因壓力波之鋼管的膨脹或收縮等之動力行為。任意耦合應(yīng)用范圍任意耦合之流固界面或耦合面不須是封閉的面，可為兩片以上的不相連接之面段，也可為一個拉格朗日域?qū)?yīng)于兩個不相連接之歐拉域。因目前的 MSC.Patran 軟件尚不能用來直接進行 ALE 相關(guān)之前處理， 故于實際建置任意耦合的計算模型時，須先經(jīng)由 Master-slave surface 的接觸過程，建立流固耦合所作用的面段，再于執(zhí)行 MSC.Patran 之后，使用文本編輯器如記事本，修改 *.DAT 文件內(nèi)之接觸相關(guān)指令，將 Contact 卡 改寫成 ALE 卡及建立 ALEGRID 卡。此外，也須包含所有相關(guān)的歐拉節(jié)點，建立 ALEGRID 組 (Group/Set)，且使用文本編輯器(如記事本)去定義ALEGRID 卡。在建置主(Master)從(Slave)接觸模型的過程內(nèi)，一般是被 動者為主，主動者為從。主從關(guān)系的建立ALE 之耦合面系經(jīng)由輸入檔內(nèi)的 ALE 卡產(chǎn)生。一般而言，ALE 耦合面之歐拉與拉格朗日兩類元素節(jié)點須是一對一，位置重合。由于 ALE 耦合面之互制移動與變形，代表固體之拉格朗日格網(wǎng)隨著時間而移動或變形，且歐拉格網(wǎng)也在空間移動或變形。即固體變形時，耦合面之位置與形狀隨之改變，并 帶動鄰接耦合面的歐拉格網(wǎng)與其他部分歐拉格網(wǎng)作相應(yīng)之移動與變形。因 此，一方面，流體材料在歐拉格網(wǎng)內(nèi)移動；另一方面，歐拉格網(wǎng)本身也在移動，以致歐拉網(wǎng)格的位置與形狀系隨時間作不斷地改變。惟因歐拉網(wǎng)格一定是體元素，故 ALE 之仿真只具有 3 個自由度(DOF)。至于歐拉元素之節(jié)點可定義為 ALE 節(jié)點之必要條件，則如下： 須為單種歐拉材料； 適用于固體變形光滑者； 歐拉材料不能具有剪切剛度； 對流場，不使用 Roe 求解器 (根據(jù) Philip Roe 教授的理論之流力計算法)。五、 一般耦合計算之應(yīng)用例5.1 固體驅(qū)動流體之仿真計算5.1.1 隧道內(nèi)之高速車輛對流場的效應(yīng)高速車輛從較空曠的環(huán)境進入空間受限的隧道內(nèi)之后，行進的車體會較顯著地壓縮前方與其周圍的空氣，以致形成明顯的活塞之氣動力（aerodynamic force）效應(yīng)，如空氣壓與行車阻力均增大，以及在隧道出口產(chǎn)生微氣壓波(micropressure waves)等，而微氣壓波可使乘客感覺耳朵不適 。此種車體進入隧道時之周圍空氣壓的波動稱為Kelvin-Helmholtz 不穩(wěn)現(xiàn)象。根據(jù)日本新干線之實測結(jié)果，微氣壓波的最大波壓約與隧道內(nèi)的車速的三次方成正比，且與離開隧道出口之距離成反比。針對此現(xiàn)象，通常的高速鐵路設(shè)計基準之一系滿足隧道內(nèi)之最大空氣壓 的變動率不超過 3.0 kPa/sec。此外，隧道內(nèi)之有軌列車若有顯著的振動， 大多系因車輛系統(tǒng)之阻尼不足與隧道內(nèi)之氣動力等引致列車之鋼輪、鋼軌與空氣流場的耦合振動。于高速列車駛?cè)胨淼纼?nèi)時，因隧道內(nèi)的空氣壓與其他氣動力特性會有明顯的改變，故列車之側(cè)向振動可能加劇，使乘客感到不適，甚至造成乘客之恐懼感及影響行車之安全。若隧道內(nèi)系鋪設(shè)雙向軌道，雙向行車，列車大致占隧道之一半空間，則車體左右兩側(cè)所受之空氣壓會顯 著不同，空氣壓差更使得車體有側(cè)向振動加劇之傾向?；谏鲜鲋紤]，本例系針對高速列車進入隧道后之情況，應(yīng) 用 MSC.Dytran，進行三維的流固耦合計算，探討空氣動力行為，亦即根據(jù) 模擬計算結(jié)果，探求隧道內(nèi)的空氣壓與空氣流速等之時間與空間的變化情 形，俾利于將來進一步研析隧道內(nèi)可能引致的高速車體振動。茲考慮由 x-、y-、z-坐標軸所構(gòu)成之三維空間，z- 與 x- 軸形成水平面，y-軸之反方向是重力方向，車行方向為 z- 軸。幾何形狀包括隧道內(nèi)的空氣流場與車體兩部分。計算例的隧道長 40 m，其橫斷面系由直徑9.2 m 之半圓與高度 2 m、寬度 9.2 m 之矩形所構(gòu)成。隧道內(nèi)的空氣流場與 進入隧道之車體均是三維的幾何實體(solids)，且均由 8 個節(jié)點且 6 個面 的 CHEXA 體元素(solid elements)構(gòu)成?？諝饬鲌龉膊捎?1,120 個歐拉的體元素，材料為理想氣體之通用材料 Ideal Gas (DMAT)，元素之物理性質(zhì)(Properties)系選用Hydro (Peuler1)。車體之模擬系使用 8 個三維之拉格朗日的體元素，其 物理性質(zhì)選用彈塑性之通用材料 ElasPlas (DMAT)。而根據(jù)高速鐵路列車 之數(shù)據(jù)，車之質(zhì)量密度采用 329 kg/m3。空氣流場與車體之有限元素模型如圖1。所考慮之車體寬度 3 m、高度 3.5 m 及長度 10m，行經(jīng)隧道之左半部 內(nèi)，如圖2 所示。圖 1 交通隧道之幾何與有限元素模型 (圖標之車體在右上角，車體外部包 裹著流固耦合面，車行方向系由圖示之右上至左下；對隧道之上游與下游斷 面均給予歐拉之 flow 邊界條件)流固耦合面系密閉地包裹車體，并由 24 片二維之假殼元素(dummy shell)構(gòu)成，緊附在拉格朗日的體元素上。至于該元素材料的物理 性質(zhì)，也選用 “Dummy Shell”。 因流固耦合面是假殼元素，故不給予材料數(shù)據(jù)。空氣之材料數(shù)據(jù)系采用伽馬爾(Gamma law)之狀態(tài)方程式(EOS, equation of state)。因空氣壓是質(zhì)量密度、比熱比與單位質(zhì)量之內(nèi)能的函數(shù)，故空氣之狀態(tài)方程式如下：p = (1) e在上式中，p = 空氣壓；圖 2 隧道與車體之橫斷面 (自隧道進口往出口看之立視圖)圖 3 第 0.65 sec 時之隧道內(nèi)面的流速分布云紋圖 = 比熱比 = 氣體之等壓比熱對等容積比熱之比值 (ratio of specific heats) =C p = 等壓之比熱；Cv = 等容積之比熱； = 空氣之質(zhì)量密度 = 1.2887 gm/cm3 ；e = 空氣單位質(zhì)量之內(nèi)能。彈塑性之車體系采用多項式狀態(tài)方程式(EOSPOL)，主要性質(zhì)如下： 質(zhì)量密度 = 329 kg/m3 剪切模數(shù) = 8.181010 N/m2 降伏應(yīng)力 = 3.50108 N/m2關(guān)于隧道內(nèi)之空氣流場，元素性質(zhì)采用 Hydro (Peuler1)， 各歐拉的體元素內(nèi)之初流速為零，各體元素內(nèi)之空氣質(zhì)量密度為 1.2887gm/cm3，并經(jīng)由 MSC.Patran 之 “Load & BC/Init. Cond. Euler.”， 分別依序建立”Shape”、”Initial Values”與 ”Region Definition” 三個 option，構(gòu)成空氣流場的起始條件。在隧道入口斷面與出口斷面的空氣流場均采用 flow 邊界條件。該邊界條件是：在歐拉的體元素與外界之接觸面上，容許空氣流入或流出。至于不給 flow 條件之歐拉域的邊界，則預設(shè)(default)為歐拉材料不 穿透之剛性壁(rigid wall)，也就是固定邊界。車體強制速度的邊界條件是：于所有時候，z- 軸向之速度 60m/sec (即 216 km/hr)作用在所有之拉格朗日的體元素節(jié)點上，該速度代 表車速，并經(jīng)由 MSC.Patran 之 ”Loads/BCs/Velocity”菜單建立。流固耦合面是由 24 片假殼元素所構(gòu)成，密閉地包裹車體，也是歐拉的空氣流場之邊界?？諝鈮鹤饔迷隈詈厦嫔?，再傳至拉格朗日的體元 素節(jié)點上。耦合邊界條件系經(jīng)由 MSC.Patran 之 Load & BC/Coupling 菜單建立。本例系采用一般耦合(general coupling)計算，即歐拉的 體元素之網(wǎng)格不隨時間而改變。至于 COVER 為 “INSIDE” ，則表示假殼元 素所包裹之空間內(nèi)，無歐拉元素存在。于 開 始 執(zhí) 行 MSC.Patran 之 前 ，經(jīng)由 MSC.Patran 之 “Analysis/Analyze/Input Deck/Translate” 的 “ExecutionControl”菜單下，設(shè)定仿真計算時間步長的總個數(shù)、起始計算之時間步長(INISTEP)與容許最小之計算時間步長(MINSTEP)等參數(shù)。本例系設(shè)定時間 步長的總個數(shù)為 4 萬個、起始之時間步長為 1.010-5 秒、最小之時間步長 為 1.010-6 秒。MSC.Dytran 之歐拉域的計算結(jié)果輸出系以每個元素為單位。為存取隧道內(nèi)的空氣流場之計算結(jié)果，于計算之前，在 MSC.Patran 之 “Analysis/Analyze/Input Deck/Translate”菜單的 “OutputRequest” 下，建立 Archive 及 “Element Output”型式之輸出文件。本例系每 1,000個時間步長，輸出計算結(jié)果一次，輸出項目包括 Pressure、x-vel、y-vel、z-vel (即空氣壓與三個坐標軸向之流速分量)等。計算完 成后，經(jīng)由 MSC.Patran 之“Analysis/Read Archive Files/Object= Model and Results/Translate”菜單，讀取相關(guān)之 Archive 檔(即*.ARC 文件)的數(shù)據(jù)。圖3 顯示第 Cycle 40,000 (第0.65 sec)之隧道內(nèi)面的 空氣流速分布的云紋圖，該圖系由隧道出口左側(cè)往進口之方向看。該云紋圖所顯示者包括：隧道左側(cè)壁、隧道頂與出口斷面之各元素的空氣平均流速。 大致上，車頭附近之空氣流速較高，車尾附近之空氣流速較低。在隧道出口斷面上，車體左側(cè)的空氣元素(#574)與右側(cè)的空氣元素(#540)之平均空氣壓的歷時曲線比較如圖4，其顯示左右兩側(cè)之空氣壓顯著不同。因車體系靠近隧道左側(cè)壁，故于各時間，車體左側(cè)的平均空氣壓明顯地較車體右側(cè)為大。圖示之空氣壓最大值系發(fā)生于第 0.35 sec，即車頭行至隧道中點之時。因本例所考慮之車行時間與隧道長度均尚短，故計 算結(jié)果尚未得到隧道內(nèi)的負空氣壓與壓力波。圖 4 隧道出口斷面的車體左側(cè)與右側(cè)元素之平均空氣壓的歷時曲線未來可能應(yīng)用 MSC.Dytran 軟件，進一步研究之隧道與高速 列車的課題如下： 隧道內(nèi)的微氣壓波之模擬計算 車體因不平衡的氣動力之振動模擬計算5.1.2 水下之固體物的高速移動本計算例所用之歐拉域系長 20 m (x 向)，寬 21 m (z 向)及高10 m (y 向)的矩形立方體，在 x-y-z 向的有限元素之網(wǎng)格數(shù)分別為 20、20、21。負 y 軸向是重力的方向。歐拉域的下部 5 m 深為水域，上部 5 m 深為模擬空氣之空 (void) 域 。因水域與空域之材料均Hydro(peuler1)，故須分別在 MSC.Patran 的 Loads/BCs 菜單下，分Shape、Initial Value 與 Region Definition 的步驟，給予歐拉的起始條件。運動速度為 10 m/sec( velocity 邊界條件)的水下高速物體系由拉格朗日元素所構(gòu)成，它長 4 m，寬 1 m 及 高 1 m ，如圖5 所示，對水域底面與上游邊界之凈距均為 1 m。歐拉域之上游面(x= 0 m)與下游面采用 flow 邊界條件，針對水之材料，x、y 與 z 向的流速分量均給予零。歐拉域內(nèi)的所有節(jié)點的起始流速也給予零。圖 5 采用 FMAT 等于 0.002488 的等值面模擬水面本例系應(yīng)用 MSC.Patran 軟件，采用流固的一般耦合計算，在水內(nèi)，用假殼(dummy shell)元素包裹高速固體物的外部，形成封閉的流固耦合面，經(jīng) MSC.Patran 的有限元素之 Equivalence 及 Verify 假殼的法線方向，使所有之法線指向外，以致假殼所構(gòu)成之耦合面無破洞；然后，藉由該耦合面及 Loads/BCs/Create/Coupling 菜單之設(shè)定，執(zhí)行流固的一般耦合計算。圖5 顯示波動的自由水面與水下之流速分量；但未顯示水下固體物之移動。換言之，Insight 工具雖顯示歐拉域內(nèi)部之動態(tài)狀況，但不顯示歐拉域范圍內(nèi)之拉格朗日固體物的運動狀況。此為目前 MSC.Patran的限制。5.1.3 造波板與之波浪水槽本計算模型是 20m(x 向)10m(y 向)3m(z 向)的波浪水槽，內(nèi)含 7m 高(y 向)、3m 寬(z 向)與 0.5m 厚(x 向)之造波板。波浪水槽之內(nèi)部是歐拉域，分成 20203 個體元素；造波板是拉格朗日域，分成 1143個體元素，如圖6。如圖所示的圓點代表造波板之外部包裹假殼的耦合面，使拉格朗日域與歐拉域產(chǎn)生一般耦合計算。造波板之所有元素節(jié)點均有正弦時變性之 x 向的速(velocity BC)，最大速度為 1 m/sec，周期為 1 sec。圖 6 模擬造波板與波浪水槽之有限元素格網(wǎng)歐拉域之下部 5 m 是水，材料是 LinFluid；上部是空氣，當作空域處理。整個歐拉域之上游面(x= 0 m)與下游面(x= 20 m)均給與予flow 邊界條件。水與空域之性質(zhì)均是 Hydro(Peler1)，搭配 MSC.Patran的Loads/BCs 菜單下之歐拉起始條件(Init. Cond. Euler)，但僅給予水域 x 向之初速度 0.1 m/sec。在代表水槽之歐拉域的下游面之 flow 邊界條件，系給予流速0.1 m/sec (x 向)；若歐拉域之上游面之 flow 邊界條件，也是 x 向流速0.1 m/sec，則第 1.0 sec 時的上游邊界附近會有水面特別壅高之現(xiàn)象，如圖7 所示。若上游面之邊界流速改為0.1 m/sec，則該壅高水面就可消除。圖 7 水槽內(nèi)的造波板(在中央者)引致約第 0.8 sec 時之水面壅高及流速的空間分布云紋圖5.1.4 海面上之高速物體撞擊混凝土墻本例之目的在顯示一個計算模型可包括兩組以上的一般耦合，且被包裹耦合面的固體物可有一部分在水域內(nèi)與模擬空氣之空域內(nèi)。本例也探討一般耦合與固體對固體之撞擊同在一個計算模型內(nèi)的狀況。計算模型之幾何體是三個較小的長方體(solid)在一個較大長方體內(nèi)，如圖8。較小的長方體由左至右，分別稱為 maker、body 與 struc，均是用拉格朗日法描述之固體。尺寸 0.5m7m3m 的 maker 是造波板，用來模擬海洋環(huán)境；1m1m1m的body，用來模擬高速運動之水上物體；尺寸 1m9m3m的struc用來模擬混凝土墻。尺寸為 20m10m3m 之較大的長方體稱為 euler，系用來模擬包括水與空氣之歐拉域，其范圍的 x 向系由 0 m 至 20 m。計算起始時的 maker 與 body 之 x 向凈距是 5.5 m，body 與 struc 之 x 向凈距是2.0 m。水上物體 body 的初速度是 60 m/sec，朝 x 坐標軸的方向。圖 8 含自由水面的水域、造波板、高速固體物與混凝土墻構(gòu)成之模型注：最大之長方體是包含水與空氣之歐拉域，較小長方體是拉革朗日域。圖8 由左至右，依序為：由 1143 個三維體元素所構(gòu)成造波板 maker、由 222 個三維體元素所構(gòu)成之高速物體 body、由 4186 個三維體元素構(gòu)成之混凝土墻 struc(詳如圖9)。雖 maker、body 與 struc均是彈塑性的通用材料 ElasPlas (DMAT)，但 struc 的材料參數(shù)與前兩者不同且含材料破壞模型。對包含水與空域之歐拉域，其在 x 值等于 0 m 與 20 m 之邊界面分別全面給予水流進與流出的 flow 邊界條件，且均有 x 向之流速 0.1m/sec，材料對象僅是水。對其它未給予邊界條件之四個面，則預設(shè)為固定邊界。同時，對所有的歐拉元素，給予歐拉之起始條件，含 “shape”、 “initial values”與 “region definition”三部分。水域內(nèi)之初流速為 0.1 m/sec。至于空域，則不給予初流速。圖 9 高速固體物與混凝土墻之有限元素格網(wǎng)本例對仿真造波板之 maker 各元素節(jié)點，給予正弦運動之velocity 邊界條件，最高速度為 1 m/sec，周期為1秒。正弦運動速度數(shù)據(jù)輸入過程系經(jīng)由 MSC.Patran 的 “Loads/BCs/Create/Velocity”菜單，配合 Field 菜單所建置之各時間的正弦函數(shù)值。對高速固體物之各元素節(jié)點，系給予初速度為 x 向 60 m/sec之荷載條件(Loads/BCs)。對混凝土墻的最底面之所有元素節(jié)點，給予 x、y 與 z 三個方向之位移量為零的 displacement 邊界條件(相當于節(jié)點固定之條件)；其前處理過程，可針對最底面之整個幾何面(即 solid 3.3)，程序會自動地選定相應(yīng)的元素節(jié)點。對高速固體物與混凝土墻之接觸及撞擊的模擬，系給予主面與從屬面接觸(master-slave surface contact)之條件。主面系位于混凝土墻上，從屬面則位于高速運動的固體物上。此外，也在 “Execution Controls/Inertial Loads”之菜單下，給予重力加速度之大小(9.81m/sec2)與其作用方向。在 MSC.Patran 的 “Execution Controls/Execution Control Parameters”下之菜單內(nèi)，系給予計算終止時間(或終止步數(shù)End Step) 0.07 sec、起始時間步長t 為 3.910-5 sec、容許之最小時間步長(Min. t) 為 1.010-14 sec。于計算之時間步長將短于 Min. t時，計算即告終止。圖10 與圖11 系應(yīng)用 MSC.Patran 的 Results 工具繪出計算模型的外圍各歐拉元素之FMAT值的云紋圖，它呈現(xiàn)高速物體引致之自由水面位置的變動與水面波。兩圖之時間分別為開始計算后之第 0.07 sec與 0.20 sec 時，圖上約中央處是造波板，靠右者是混凝土墻；因高速物體系位于歐拉域的內(nèi)部，故不能經(jīng)由 Results 工具看到它。于第 0.06 sec時，初速度 60 m/sec之固體物已貫穿混凝土墻。圖12 顯現(xiàn)混凝土墻受到固體物接觸及貫入后，混凝土墻之材料破壞與塑性應(yīng)變的分布情形。因水與混凝土之阻尼效應(yīng)，在接觸及撞擊過程內(nèi)，固體物的速度會隨時間而降低。圖 10 于第 0.07 sec 時之 FMAT 云紋圖顯示自由水面波動情形圖 11 于第 0.20 sec 時之 FMAT 云紋圖顯示自由水面的波動圖 12 高速固體物貫穿混凝土墻與混凝土墻之塑性應(yīng)變分布圖 13 Insight 工具顯示計算模型內(nèi)部與 FMAT 等值面所代表之自由水面圖13 是采用 MSC.Patran 的 Insight 工具所產(chǎn)生之影像，它顯示第 0.04 sec 時的自由水面之波動與水下的流速向量。產(chǎn)生該圖像文件的過程是：于使用 MSC.Patran 的 “Insight Application”菜單之 “Action/Tool” 制 作 Insight 工 具 之 后 ， 選 取 menu bar 上之 “insight control/animation control” ( 動 畫 控 制 ) 下的 “animation setup” (動畫建置)，然后，點選 “save frames to files”，再點選 “animate”，如此，若有 8 幅動畫之畫面(frames)，就會建立且儲存 8 個圖像文件。如圖13 所示，Insight 工具不顯示歐拉域內(nèi)之拉格朗日的固體物之移動，此為目前之 Insight 工具的限制。于第 0.03 sec 時，高速固體物尚未接觸混凝土墻。于第 0.04sec時，高速固體物約貫入混凝土墻厚度之一半。于第 0.06 sec 時，高速固體物已完全貫穿混凝土墻，且離開混凝土墻。固體物之速度引致流速變化。在 0.07 sec 的過程內(nèi)，流場內(nèi)元素得最高流速變化如表1。表 1 各時間之流場內(nèi)的元素最高流速時間(sec)0.000.010.020.030.040.050.060.07流速(m/sec)0.136.660.538.128.523.722.018.15.2 應(yīng)用流固耦合功能形成流場較復雜之固體邊界由于含不規(guī)則之固體邊界之流場空間不能劃分成矩形格網(wǎng)(iso-mesh)如圖14，不利于給予含自由液面之歐拉的起始條件，為避免此項限制，可將流體空間先采用 simple solid (即四面體、五面體或六面體)定義成歐拉域；然后，在該歐拉域內(nèi)，迭上不規(guī)則之拉格朗日 solid，形成不規(guī)則之流場的固體邊界。如此，對流場空間仍可使用矩形格網(wǎng)。圖 14 有限元素的矩形格網(wǎng)之平面5.2.1 上游側(cè)與下游側(cè)水深相等情況之潛堰包含水與空氣的歐拉域之矩形體系長 20 m(x 向)、高 3 m(y向)及寬 3 m(z 向)，分成 20123 個元素。歐拉域的下部 1.5 m 是水域，上部之 1.5 m是代表空氣之空域。堰是主軸與流向垂直之結(jié)構(gòu)物，用以抬高水位，俾藉重力取水及引水。拉格朗日的混凝土堰體尺寸是 1m1m3m，其上游面位于 x= 9 m 處，堰體被劃分成 223 個元素，如圖15，堰體上的餅圖樣代表它被一般耦合的接口所包裹著經(jīng)由 MSC.Patran 的 Loads/BCs 菜單，在歐拉域的上游與下游邊界面，分別給予 flow 邊界條件；針對水之流進與流出，給予流速10 m/sec。歐拉域下部 1.5 m 深之水域也均給與初流速 10 m/sec。另由Analysis 之 “Inertial Loads”表，給予重力加速度 9.81 m/sec2，作用在負 y 軸向。相關(guān)的是須由 Loads/BCs，對堰體底面節(jié)點給予位移量為零之邊界條件；否則，堰體會因重力作用而往下掉落至歐拉域的范圍外。圖 15 拉格朗日域(含耦合面之堰體)及歐拉域(水與空氣)的元素格網(wǎng)圖 16 第 0.026762 sec 時的越過潛堰的水域之流速空間分布圖 17 采用 FMAT 為 0.002488 之等值面顯示越過潛堰之水面形狀本例之計算起始時間步長為 110-4 sec，容許最小時間步長為 110-14 sec。根據(jù)計算結(jié)果，圖16 是 Cycle 240 (第0.026762 sec)時之流速空間分布，在堰體附近之流速提高，最高達 13.1 m/sec。由于水下潛堰的效應(yīng)，越過潛堰處之水面會起伏，如圖17 所示；該圖顯示每個元素內(nèi)之 FMAT 值等于 0.002488 的等值面，其接近自由水面。應(yīng)用 MSC.Dytran 之優(yōu)點是不但可模擬分析流況，且可同時模擬計算混凝土堰體之應(yīng)力與應(yīng)變量。惟 MSC.Dytran 之缺點是針對由上游往下游之沖擊水流，適宜用在瞬時(transient)行為之模擬計算；想由瞬時計算至穩(wěn)態(tài)(steady)，往往會遇到流速變得很高(即計算發(fā)散)與計算時間步長t 變得很小，以致計算中止的情況。本例于 Cycle 570(第 0.063588 sec)以后，就有計算顯著發(fā)散的現(xiàn)象，故該時間以后之計算結(jié)果不可信。圖 18 水域與空域之歐拉元素及堰體之拉格朗日元素格網(wǎng)5.2.2 上游側(cè)與下游側(cè)水深不等情況之固定堰關(guān)于本例的整個計算模型之材料，除水與混凝土兩種外，還有空氣是性質(zhì) Hydro (Peuler1)及歐拉體(solid)的空域。本例的歐拉域系長 20 m(x 向)、寬 3 m(y 向)及高 10 m(z 向)之長方體。而采用拉格朗日法描述之堰體系高 1.5 m(y 向)、厚為 1 m(x 向)及長 3 m(z 向)。堰體完全在歐拉域之內(nèi)，且堰底位于歐拉域之底面。歐拉域之元素數(shù)目為 20203，共有1,200個矩形元素，拉格朗日之元素數(shù)目為 233，如圖18 所示。經(jīng)由 MSC.Patran 之 Loads/BCs 菜單下的 flow 邊界條件，給予整個歐拉域上游邊界面(x = 0 m) 之流速 4 m/sec 與下游邊界面(x = 20 m)之流速為零，flow 邊界條件之材料均僅針對水。至于圖18 下部中央之堰體外部符號，則顯示堰體被一般耦合面包裹著，使堰體與歐拉域有流固的一般耦合效應(yīng)。圖 19 起始 Euler1 群組(圖左)與 Euler2 群組(圖右)的水域與空域定義圖 20 歐拉域于第 0.58 sec 時之 FMAT 數(shù)值圖顯示瞬時的水面形狀在 MSC.Patran 的 Analysis 菜單下，對整個計算模型給予重力荷載 9.81 m/sec2，重力方向是負 y 軸向。以堰厚中央位置(x= 10 m)為界線，將將歐拉域分成 Euler1與 Euler2 兩個群組(group)，Euler1 由 x= 0 m 至 x= 10 m，Euler2由 x= 10 m 至 x= 20 m，俾給予兩部分不同之歐拉起始條件(Init. Cond.Euler.)，如圖19。Euler1 的下部 6 m 是水域，起始流速為 4 m/sec；上部 4 m 是空域。Euler2 的下部 1 m 是水域，起始流速為零；上部 9 m 是空域。圖20 是 Cycle 6300(第 0.58 sec)時之水流剖面，那時候，混凝土堰體之最大有效應(yīng)力為 1.58104 N/m2。本例所用之計算起始時間步長t 為110-5 sec，容許之Min.t 為 310-16 sec。計算結(jié)果，于時間超過 Cycle 6300(第 0.576261sec)，計算時間步長小于 Min. t，因而計算中止。本例之重點在于： 使用 Euler1 與 Euler2 兩個群組，去定義歐拉域的上游部分與下游部分的不同水深之邊界條件與起始條件。對此例，值得思考是： 下游部分的 Euler2 之水域的起始流速若不是零，而是 4 m/sec或其他數(shù)值，則所代表之意義及計算結(jié)果之差異程度如何。5.2.3 上游側(cè)與下游側(cè)水深不等情況之階梯式渠道如同前節(jié)之固定潛堰，所用歐拉域系長 20 m(x 向)、寬 3 m(y向)及高 10 m(z 向)之長方體。惟有四個長方體連在一起及構(gòu)成拉格朗日域。每個長方體均厚 2 m(x 向)及寬 3 m(z 向)，而高度(y 向)分別為 4 m、3m、2 m與1 m，形成階梯狀，如圖21 之左下方。歐拉域類似前節(jié)之固定潛堰，將 x等于0 m至10 m之范圍定義為 Euler1 群組，且定義下部分 5 m 深為水，6 m/sec 之上游面 flow邊界條件與該水域初流速；上部分之5 m深為代表空氣之空域。圖 21 四個幾何長方體構(gòu)成一個階梯狀之有限元素體與歐拉域之階梯狀邊界類似地，歐拉域 x 等于 10 m 至 20 m 之范圍定義為 Euler2群組，且定義下部分 1 m 深為水，3 m/sec 之下游面 flow 邊界條件與該水域初流速；上部分之 9 m 深為代表空氣之空域。前述之構(gòu)成拉格朗日域的四個長方體，材料使用混凝土，它們一起建置有限元素格網(wǎng)，直接用所形成之階梯狀體元素群之外露的元素面作為耦合面，進行 Loads/BCs 菜單下之一般耦合建置。此方式省去在拉格朗日之體元素外部建構(gòu) dummy shell 的耦合面、equivalence 與 verifynormal 等程序，仍可完成流固之一般耦合，如圖21；然而，此方式是否到處行得通(universal)仍待進一步探討。如圖21 所示者是 Cycle 780 時的水域范圍圖，大于 Cycle 780 者之時間步長小于容許之 Min.t，因而計算中止。本例之重點目的在于用階梯狀的拉格朗日之體元素，形成歐拉域之固定邊界，而歐拉域仍可采用矩形之元素格網(wǎng)，利于設(shè)定歐拉域的 flow邊界條件與內(nèi)部之起始條件。此外，如同前節(jié)，對歐拉域內(nèi)的不同范圍可分為兩個以上之群組，俾設(shè)定不同之起始條件與 flow 邊界條件。5.2.3 固定開度之水利閘門固定開度之水利閘門(hydraulic gate)也是歐拉域的固定邊界之一種。本例也是流固之一般耦合分析的應(yīng)用。流體是水與代表空氣的空域，固體是閘門。垂提式閘門開度固定為 2 m，閘門上游側(cè)之水深為 5 m，水之流速是在 x 向為 20 m/sec。本例之歐拉域的幾何是 x 向為 20 m 長，y 向為 10 m 高，z向為 3 m 寬，采用四節(jié)點六面體之有限元素共 20203 個，性質(zhì)為歐拉之體元素 Hydro (Peuler1)。在幾何體之上游與下游邊界，分別給予 flow 邊界條件，具有 x 向之流速 20 m/sec。該歐拉域在 x 向，分為各長 10 m 之兩個群組(group)，稱為 euler1 與 euler2。模擬閘門之幾何體是 x 向為 0.3 m 長，y 向為 8 m 高，z 向為 3 m 寬，采用四節(jié)點六面體之有限元素共 1163 個，性質(zhì)為拉格朗日之體元素。閘門之幾何體外部包裹假殼(dummy shell)元素所構(gòu)成之閉合面，經(jīng)有限元素之 verify，讓各法線(normal)之方向均朝外，確保無破洞，作為流固耦合面。對 euler1 群組之上部 5 m 深，定義為 void1 域之空隙；下部5 m深，定義為water1域之水，且初始有x向之流速20 m/sec。euler2 群組，則全部定義為 void2 域之空隙，不給予初始流速。圖 22 第 0.027737 sec 時之閘門附近的水流縱剖面本例目的之一是驗證：不但 water1 域之水面(即 water1 域與 void2 域之接口)會隨著時間變動，且 water1 域之水可橫向地往 void2域擴展。如圖22 是 cycle 1300 (時間為第 0.027737 sec)之水體立剖面。本例提出群組(group)之應(yīng)用；若非應(yīng)用群組，則無法對上、下游兩部分之歐拉域設(shè)定不同之起始條件，即上游之 euler1 域的水深為5 m，下游之 euler1 域的水深為零。重點之一是： 閘門體須與 water1 域重疊一部分，而與 void2域重疊數(shù)。本例系將 0.3 m 厚之閘門體完全與 water1域重疊之下游面重合。本例的難處是：因用顯式的時間積分法，計算之時間步長愈來愈小，以致算不下去，無法算至流場達穩(wěn)態(tài)之時。譬如，本例計算迄 cycle 300 (時間僅第 0.027737 sec)，接著之計算時間步長小于 1.010-14sec，而中止計算。一般而言，MSC.Dytran 對固體不動，但流體運動之模擬，其流速計算較易發(fā)散(diverge)，譬如本例，流場上游與下游的邊界流速與初始的 x 向之流速均為 20 m/sec，計算達 cycle 130 (時間僅第0.002639 sec)，最大流速即發(fā)散至 7.73103 m/sec。六、 任意耦合計算之應(yīng)用例大致上，任意耦合分析可應(yīng)用在氣體或液體的流場驅(qū)動固體物之互制狀況。執(zhí)行任意耦合計算時，拉格朗日域與歐拉域不須迭接，只須相連接。本文所述之兩個例子包括隧洞內(nèi)之氣爆壓作用在鋼板上、隧洞內(nèi)之水流推動固體物等兩項如下。6.1 隧洞內(nèi)之氣爆壓作用在鋼板上本例是進行包裹在 3m3m3m 幾何立方體(solid 1)外部的殼元素之拉格朗日域如圖23 之左半部，以及模擬空氣的 3m3m3m 幾何立方體(solid 2)之歐拉域如圖23 之右半部，兩者之任意耦合分析。兩者的相連接處是任意耦合面，即 solid 1 的左側(cè)面與 solid 2 的右側(cè)面之 x 坐標分別為 x 等于零與 9 m。初始的爆炸源之半徑為 0.05 m，中心點系位于 solid2 幾何體的中央，即 x= 4.5 m、y= 1.5 m 與 z= 1.5 m 處；該坐標數(shù)值系采用文本編輯器在 *.dat 檔內(nèi)的 TICEUL 卡上，予以鍵入。殼元素是彈塑性的通用材料 ElasPlas(DMAT)之鋼，歐拉域的材料是符合理想氣體律之空氣。圖 23 第 0.008099 sec 時的殼形固體物與爆炸之空氣外圍的位移量分布圖 24 第 0.008099 sec 時之殼體變形與位移量之分布任意耦合面之建置系經(jīng)由 MSC.Patran 的主面與從屬面之接觸(Contact/Element Uniform/Option= Master-Slave Surface)過程。主面是二維殼形體之右側(cè)面段，可選相應(yīng)之幾何面solid 1.2，予以設(shè)定。從屬面則是三維歐拉體元素群之左側(cè)外露面所構(gòu)成之面段，可選相應(yīng)之幾何面solid 2.1，予以設(shè)定。根據(jù)歷時的流固互制模擬計算結(jié)果，Cycle 200 (第 0.008099 sec)時，殼形體接觸爆炸壓的面之最大位移量是 0.372m，如圖24 之云紋圖所示。6.2 隧洞內(nèi)之水流推動固體物本例系應(yīng)用任意耦合功能，仿真隧洞內(nèi)之水流推動固體物的行為。如圖25，左邊是計算模型的上游側(cè)，右邊是下游側(cè)。計算模型系由兩個幾何長方體與一個立方體所構(gòu)成，即近乎在中央之立方體 soli