高效偉 丁金興 劉華雩

(大連理工大學(xué)航空航天學(xué)院,工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室,大連 116024)

本文創(chuàng)建一種全新的數(shù)值計算方法—有限線法,并將其用于求解流體-固體一體化耦合傳熱分析.常用的有限元法是基于體單元的離散方法,有限容積法是在控制體面上運算的方法,邊界元法是在邊界面上離散的方法,無網(wǎng)格法等是由計算點周圍的散點構(gòu)建計算格式的方法.本文提出的有限線法是一種基于有限條線段離散的方法,在每個點只需要兩條或三條直線或曲線構(gòu)成的線系,則可建立任意高階的算法格式.創(chuàng)新性思想是: 通過采用沿線段求方向?qū)?shù)的技術(shù),由一維拉格朗日插值公式,建立二維和三維問題的任意高階線系的空間導(dǎo)數(shù),并以此直接由問題的控制偏微分方程與邊界條件建立離散的系統(tǒng)方程組.有限線法理論簡單、通用性強,能以統(tǒng)一的格式求解固體與流體力學(xué)等偏微分方程邊值問題.在流體方程中,擴散項采用中心線系保證高精度計算,而對流項則采用迎風線系體現(xiàn)其方向性特征.本文將給出有限線法求解流固耦合傳熱問題的幾個算例分析,驗證其正確性與有效性.

1 引言

根據(jù)算法離散與操作維度的不同,常用的數(shù)值計算方法可歸納為4 類: 基于體單元的計算方法(如有限元法[1,2]、有限塊法[3,4]、單元微分法[5,6]等),基于面單元的方法(如有限容積法[7,8]、邊界元法[9,10]、邊界面法[11]等),基于線單元的方法(如有限元線法[12]、交叉線法[13]等),以及基于散點計算的方法(如無網(wǎng)格法[14-16]、基本解法[17]、奇異邊界點法[18]等).

在基于體單元的計算方法中最具代表性的方法是有限元法(FEM).用其分析問題時將計算域劃分成有限個體單元,物理量的表征、微分是通過與計算域維度一樣的體單元實現(xiàn)的.鑒于其強大的解題能力,FEM 已被用于解決各種工程問題,特別是固體力學(xué)問題,這主要歸因于其使用了具有良好性能的等參單元進行幾何離散和物理量插值[19],使得計算結(jié)果非常穩(wěn)定.此外,FEM 所需要的單元和節(jié)點數(shù)遠少于有限容積法(FVM)和有限差分法(FDM).FEM的弱點主要體現(xiàn)在: 1)需使用變分原理或虛功原理來建立計算格式,對于不同的問題所建立的表達式不同,從而不便于建立多場耦合問題的統(tǒng)一算法;2)用傳統(tǒng)的FEM 求解含有對流項的問題時,不易實施迎風格式與激波捕捉計算.為此,近年來間斷伽遼金有限元法(DGFEM)在流體力學(xué)中得到了快速發(fā)展[20],其在每個單元內(nèi)用獨立定義的節(jié)點進行高階函數(shù)逼近,因而可以較容易實現(xiàn)迎風計算,然而DGFEM 計算效率很低,此外如果單元內(nèi)有激波的存在,內(nèi)部間斷使得收斂性和魯棒性大幅下降.

面單元計算方法是指主要在面上進行運算的方法.將FVM 歸入面單元計算方法在維度上似乎不太協(xié)調(diào),但是科學(xué)的.這是因為FVM 是一種將問題的控制微分方程的守恒形式(全微分形式)使用高斯散度公式將控制體積分轉(zhuǎn)換成在控制體面上進行積分運算的算法[7].FVM 是一種通量守恒性最好、激波間斷面最容易捕捉的數(shù)值方法,也由此成了目前計算流體力學(xué)中的主流方法.BEM 是一種典型的面單元計算法,是通過兩次分部積分和使用高斯散度公式,并取權(quán)函數(shù)為問題基本解而形成邊界積分方程的算法.BEM 幾何離散簡單,能很好地模擬斷裂力學(xué)中的應(yīng)力集中現(xiàn)象,并能有效求解無限域與半無限域問題[21].BEM的弱點是在求解有源(熱源、體積力、化學(xué)反應(yīng)生成項等)、非線性、非均質(zhì)問題時,會有域積分出現(xiàn)在積分方程中,削弱了其只在邊界上離散的優(yōu)點.

基于線單元的計算方法籠統(tǒng)地講可以包括有限元線法[12]以及梁單元、桿單元[22]等的算法.前者是借助于能量泛函的變分,將控制偏微分方程半離散化為用結(jié)線表示的常微分方程組求解法,而后者則是通過一維離散模擬具有一維維度的構(gòu)件的方法.真正通過線單元的離散實現(xiàn)高維(二維和三維)問題的微分方程的求解方法是作者最近提出的交叉線法(CLM)[13].該方法通過在每個點建立由兩條(二維問題)或三條(三維問題)相交的線段,并通過形函數(shù)構(gòu)造法建立多維問題的形函數(shù),然后采用微分法則,導(dǎo)出各階偏導(dǎo)數(shù)計算式,最后將其代入控制微分方程與邊界條件,建立問題的離散系統(tǒng)方程組.CLM 屬于配點型的強形式計算方法,目前是通過單元的形式在自由單元法中使用[13],用到的線系可以是直線,也可以是曲線.因而,該方法具有幾何離散簡單、適應(yīng)性強、編程方便、系數(shù)矩陣稀疏等特點,是一種具有開發(fā)潛力的計算方法.然而,根據(jù)單元形函數(shù)構(gòu)造法(即滿足Delta符號性質(zhì)的解[13])很難構(gòu)造三階及其以上的形函數(shù),這給流體力學(xué)中經(jīng)常用到高階格式的算法帶來了限制.

基于散點的計算方法有無網(wǎng)格法(MLM)等[14-18],該類方法不需要有規(guī)則節(jié)點分布的單元,只需一些任意分布的散點則可建立解題算法.MLM的形函數(shù)及其偏導(dǎo)數(shù)是通過最小二乘法等技術(shù)在一定范圍(緊支域)內(nèi)的點集構(gòu)建的.基本解法與奇異邊界點法等則是通過問題微分算子的基本解在虛擬邊界或真實邊界節(jié)點上配置方程進行求解的.MLM在幾何適應(yīng)方面具有獨特的優(yōu)勢,因而在復(fù)雜幾何模擬與動邊界問題方面得到了廣泛的應(yīng)用.MLM雖然有使用靈活、幾何適應(yīng)性強的優(yōu)點,但也有計算效率低、結(jié)果穩(wěn)定性差以及不易施加邊界條件與各向異性物性的弱點.

值得介紹的一種新的配點型方法是作者近年來提出的自由單元法(FrEM)[23,24].該方法吸收了FEM 與配點型MLM的優(yōu)點,即: 在函數(shù)表征方面,采用FEM 中有規(guī)則節(jié)點分布的等參單元,因而具有幾何適應(yīng)性強、計算效率高、結(jié)果穩(wěn)定性好、容易施加各向異性物性以及邊界條件的優(yōu)點;在算法方面,采用配點法技術(shù),直接由控制方程與邊界條件逐點產(chǎn)生系統(tǒng)方程.FrEM 在每個點只需要一個和周圍自由選擇的節(jié)點形成的獨立的等參單元,因而不需要考慮物理量在單元之間的相互連接關(guān)系以及導(dǎo)數(shù)連續(xù)性問題.此外,單元的自由性也使其更容易構(gòu)造高階迎風格式計算流體力學(xué)方程中的對流項.FrEM 已在固體力學(xué)[25]、流體力學(xué)[24]與傳熱學(xué)[26]中得到了成功的應(yīng)用.然而,由于高階自由單元是由結(jié)構(gòu)化節(jié)點形成的[23],因此三維問題的高階單元需要較多的節(jié)點,這給幾何建模與數(shù)據(jù)存儲帶來了很大的不便.

本文在作者多年來研究BEM、FEM、配點型MLM 等方法的基礎(chǔ)上,提出了一種全新的數(shù)值方法—有限線法(FLM)[27].該方法在每個點只需要1—3 條交叉線段(稱為線系)則可建立一維到三維問題的算法格式.由于每條線段用拉格朗日多項式插值公式構(gòu)建函數(shù)表征關(guān)系,因而可以很方便地建立任意高階的算法.FLM 具有CLM 和FrEM的全部優(yōu)點,還有類似于FDM 容易使用的性能,并有高階線系所用節(jié)點少、容易施加不同物性與邊界條件、以及容易構(gòu)建迎風格式的特點,因此既適合于求解流體力學(xué)問題、也適合求解固體力學(xué)問題.

一些現(xiàn)代裝備結(jié)構(gòu),如發(fā)動機渦輪等,一直處于高溫服役環(huán)境,因而使用冷卻技術(shù)提高其耐高溫能力是經(jīng)常使用的重要手段,結(jié)構(gòu)內(nèi)部設(shè)置冷卻通道是最常用的熱防護形式,對其進行流體域-固體域間耦合傳熱分析是結(jié)構(gòu)設(shè)計的重要依據(jù).本文將系統(tǒng)介紹FLM 在求解冷卻通道結(jié)構(gòu)流固域間耦合傳熱方面的應(yīng)用,為現(xiàn)代裝備結(jié)構(gòu)設(shè)計提供有力的分析手段.

2 流體與固體中的傳熱偏微分方程邊值問題

流體與固體中的熱傳導(dǎo)微分方程可統(tǒng)一表示為[5,24]

式中,ρ為流體密度,cp為定壓比熱,vi為流速,T為溫度,λij為導(dǎo)熱率,s為熱源;式中的重復(fù)指標表示求和,Ω為計算域,對于固體vi為0.(1) 式中的左端為對流項,右端第一項為擴散項[7].

相關(guān)的邊界條件可表示為

其中,?Ω為Ω的邊界,ni為邊界外法線方向矢量,h為換熱系數(shù),T∞為環(huán)境溫度.

從 (1)—(4) 式可以看出,控制微分方程與邊界條件涉及到的關(guān)鍵項是溫度對空間總體坐標的一階與二階偏導(dǎo)數(shù),本文將用有限線法建立各階偏導(dǎo)數(shù)計算式.

3 計算空間偏導(dǎo)數(shù)的有限線法

作者在文獻[5]中導(dǎo)出了一般拉格朗日等參單元形函數(shù)對總體坐標的偏導(dǎo)數(shù)解析表達式,但可以證明其一階偏導(dǎo)數(shù)只與通過計算點的交叉線上的節(jié)點有關(guān)[13],因此本文采用交叉線系建立各階空間偏導(dǎo)數(shù)計算式.

3.1 交叉線系的建立

本文提出的FLM 是一種配點型的計算方法,計算時與FDM 和MLM 類似,首先通過網(wǎng)格生成技術(shù)[28],將計算域劃分成一系列配置點(CP),圖1給出了一個二維問題的節(jié)點劃分示例.

圖1 將計算域劃分成一系列配置點Fig.1.Discretizing computational domain into a series of collocation points.

然后,對每個CP,建立由兩條(二維問題)或三條(三維問題)線段組成的交叉線系.圖2 和圖3分別展示了一個內(nèi)部點的二維(d=2)和三維(d=3)的交叉線系.原理上,線系中的線段相互間越垂直計算精度越高,但只要不共線(二維線系)和不共面(三維線系)則可獲得滿意的結(jié)果.

對于位于邊界和角點上的CP,則需要將圖2和圖3 中的交叉點移到邊界和角點上,具體做法可參見文獻[13].

圖2 兩條線組成的二維(d=2)線系Fig.2.Line set of 2D (d=2) consisting of two lines.

圖3 三條線組成的三維(d=3)線系Fig.3.Line set of 3D (d=3) consisting of three lines.

3.2 基于線系的空間導(dǎo)數(shù)計算公式

在線系的每一條線上,定義m個節(jié)點(如圖2和圖3 所示),并由下列拉格朗日多項式插值公式表達坐標x與物理量u的函數(shù)變化關(guān)系[5,23]:

式中,l為從節(jié)點1 開始計算的線段長度,uα為第α個節(jié)點的u值.Lα為拉格朗日多項式插值函數(shù):

其中,lα為第α個節(jié)點的長度.以二維線系為例,標記函數(shù)u對總體坐標xi的一階偏導(dǎo)數(shù)分量為?u/?x1和?u/?x2,則沿線段l的方向?qū)?shù)?u/?l可表示為

對于每一條線可建立一個上述方程,當d條線考慮后可得到如下的矩陣方程:

對(9)式求逆,并考慮到(6)式,則可得如下的分量表達式:

為了方便使用,可將(11)式簡寫為

式中,重復(fù)指標α表示對M1個節(jié)點求和,xc表示線系中各線交叉處(即配點)的坐標稱為一階導(dǎo)數(shù)算子.需要說明的是,矩陣J是一個2×2 或3×3的矩陣,可以很容易求得其逆矩陣的解析表達式.

對于二階導(dǎo)數(shù),可以在一階導(dǎo)數(shù)公式上再次作用導(dǎo)數(shù)算子,得到:

需要說明的是,(13)式與(15)式中重復(fù)指標α求和的節(jié)點數(shù)量是不一樣的.圖4 展示了一個配置點(節(jié)點3)采用5×3 線系離散中,一階導(dǎo)數(shù)相關(guān)的7 個節(jié)點(由實心點標記的節(jié)點)與21 個二階導(dǎo)數(shù)相關(guān)的節(jié)點(由實心點+空心圓點標記的節(jié)點)的示意圖.對于更高階的導(dǎo)數(shù)算子可由類似的遞歸方法形成[27].

圖4 一階與二階導(dǎo)數(shù)相關(guān)的節(jié)點Fig.4.Related nodes of the 1st and 2nd order partial derivatives.

4 流固耦合傳熱問題中的有限線法公式

將空間導(dǎo)數(shù)計算公式(13)和(15)直接代入流固耦合傳熱方程(1)及其邊界條件(2)—(4),則可建立耦合問題的代數(shù)離散方程.不過,流體傳熱方程中存在對流項,需要用迎風格式才能精確求解,而擴散項用中心格式才能得到穩(wěn)定的計算結(jié)果.本研究通過采用兩套線系的技術(shù)解決此問題,即使用迎風線系體現(xiàn)計算點上游對其的影響,用中心線系保證擴散項的穩(wěn)定性與精度.圖5 展示了采用5×4 節(jié)點線段形成的兩套線系的節(jié)點分布,其中紅色的節(jié)點屬于中心線系(central line set),黑色的屬于迎風線系(upwind line set).

圖5 同一配置點的迎風與中心線系Fig.5.Upwind and central line sets at a same collocation point.

采用 (13) 式和(15) 式可將方程(1)—(4)離散成下列形式:

式中分別為基于迎風線系和中心線系計算的導(dǎo)數(shù)算子,其中的q為

耦合計算時,將流體域與固體域進行一體化離散.對于固體域,(17)式的第1 個方程中的第1 項(對流項)為0.

將上述方程中的配置點xc作用于流體域與固體域的所有配置點,并在流-固界面上使用下列溫度協(xié)調(diào)條件與通量平衡方程(其中上標f 與s 分別代表流體與固體):

則可形成如下的一體化分析矩陣方程:

求解上述方程則可獲得所有節(jié)點的溫度值.需要說明的是,由于每個配置點的線系所涉及到的節(jié)點數(shù)比傳統(tǒng)的數(shù)值方法FEM 和FVM 少得多,因此 (20) 式中的系數(shù)矩陣A是一個非常稀疏的非對稱矩陣.

5 算例分析

根據(jù)本文所推導(dǎo)的公式,采用FORTRAN 語言編制了有限線法程序,程序中采用了分區(qū)技術(shù)[28],可對復(fù)雜問題進行計算分析.下面對現(xiàn)代高溫裝備中經(jīng)常使用的冷卻管道結(jié)構(gòu)中的一些流固域間耦合傳熱問題進行分析,用以驗證本文所述方法的正確性與有效性.計算中,每個配置點的每條迎風線用二階精度,中心線用四階精度.

5.1 單根管道結(jié)構(gòu)流動傳熱分析

本文的第一個算例是用于驗證本文提出的FLM的正確性與有效性.算例的幾何尺寸與邊界條件如圖6 所示,管道的上下外表面處于環(huán)境溫度為900 K的對流換熱環(huán)境中,流固界面的換熱系數(shù)為1100 W·(m2·K)—1,流體域的物性為:λ=0.6 W·(m·K)—1,ρ=1000 kg/m3,cp=4200 J·(kg·K)—1;固體域為λ=200 W·(m·K)—1.

圖6 單根管道結(jié)構(gòu)尺寸與邊界條件Fig.6.Dimensions and B.C.of a single channel structure.

為了提供參考結(jié)果進行比較,本文用FLUENT軟件采用較為稠密的網(wǎng)格對該問題進行計算,其中在x方向(水平方向)和固體域,網(wǎng)格為等間距分布,流體域靠近壁面網(wǎng)格逐漸加密.在x方向劃分了151 個點,在垂直方向: 固體域20 個點,流體域40 個點.

5.1.1 FLM 使用均勻網(wǎng)格的網(wǎng)格收斂性分析

首先基于均勻網(wǎng)格對該問題進行FLM 計算分析,考察其網(wǎng)格收斂性.計算發(fā)現(xiàn),結(jié)果只對流體域中的垂直方向的網(wǎng)格敏感.因此,在FLM 計算中,在x方向流固域都使用51 個點;在垂直方向,固體域使用7 個點,流體域分別使用21、41 和71 個點進行計算,這3 種情況分別標記為FLMfine1,FLM-fine2 和 FLM-fine3.圖7 和圖8 分別給出了流速為v=1 m/s 時3 種網(wǎng)格的計算結(jié)果.可以看出,網(wǎng)格FLM-fine3 收斂到了FLUENT的計算結(jié)果.此外,變化趨勢顯示,FLM 具有很好的網(wǎng)格收斂性.

圖8 不同網(wǎng)格尺寸下流固界面溫度分布Fig.8.Temperature distribution on the fluid-solid interface under different mesh sizes.

從圖7 可以看出,溫度為300 K的流體流過冷卻通道后,可以將處于900 K 熱環(huán)境的結(jié)構(gòu)外表面溫度降低到400 K 以下,充分說明了高溫結(jié)構(gòu)中采用主動冷卻管道的冷卻效果.

圖7 不同網(wǎng)格尺寸下管道外表面溫度分布Fig.7.Temperature distribution on the channel outer surface under different mesh sizes.

5.1.2 使用比例間距網(wǎng)格計算

為了考核FLM 在非均勻網(wǎng)格中的表現(xiàn),在上述FLM-fine1的網(wǎng)格中,允許流體域與固體域的計算點從各域邊界到其中線按比例變化,相鄰間距比例因子在水平與垂直方向分別為1.2 和1.5.圖9展示了流固域左端與中間部分的FLM 線系連成的局部網(wǎng)格圖.圖10 和圖11 分別為不同流速下管道外表面和流固界面的溫度沿x方向的分布曲線.為了比較,也繪出了用同等節(jié)點數(shù)但等間隔分布的網(wǎng)格的計算結(jié)果(FLM-uniform).

圖9 管道結(jié)構(gòu)左端與中間局部網(wǎng)格放大圖Fig.9.Enhanced meshes of the left and middle parts of the channel structure.

可以看出,用比例間距的計算結(jié)果 (FLM-ratio)與Fluent的計算結(jié)果吻合的很好,而用等間隔網(wǎng)格的計算結(jié)果(FLM-uniform)在如此少的節(jié)點數(shù)下的計算結(jié)果不可接受.可見,在FLM 計算中,采用不等間距網(wǎng)格可以大量減少所需的節(jié)點數(shù).不過,從圖10 和圖11 中的v=0.1的結(jié)果可以看出,對于流速慢的情況,即使用等間距網(wǎng)格也可以獲得滿意的計算結(jié)果.

圖10 不同流速下管道外表面溫度分布Fig.10.Temperature distribution on channel outer surface under different velocities.

圖11 不同流速下流固界面溫度分布Fig.11.Temperature distribution on fluid-solid interface under different velocities.

5.2 多根管道結(jié)構(gòu)流動傳熱分析

本文的第2 個算例是含有3 根管道的冷卻結(jié)構(gòu),如圖12 所示.三根管道的直徑均為2 mm,其中,中間管道具有5.7°的傾斜角.整個冷卻結(jié)構(gòu)的左、右邊界除了流體的入口溫度Tin部分外其他部分均為絕熱邊界條件(q=0).三根管道中的流體的密度與比熱與算例5.1 中的一樣,均為ρ=1000 kg/m3與cp=4200 J·(kg·K)—1.計算中,整個冷卻結(jié)構(gòu)分為7 個計算域,其中第2,4,6為流體域,其他為固體域.幾何離散中,所有域沿水平方向(x方向)劃分為等間隔的101 個點,沿垂直方向劃分為從各域壁面到其中線間距比例因子為1.2的31 個點.圖13為FLM 分析中由所有線系連成的網(wǎng)格圖,注意到由于中間管道的傾斜使得計算域3—5的網(wǎng)格都成了非正交網(wǎng)格.

圖12 三根管道結(jié)構(gòu)尺寸與邊界條件Fig.12.Dimensions and B.C.of a three-channel structure.

圖13 三管道結(jié)構(gòu)FLM 分析線系連成的網(wǎng)格圖Fig.13.Mesh connected by line sets of all points for FLM analysis of the three channel structure.

5.2.1 流、固域采用相同條件的計算分析

首先考察整個管道結(jié)構(gòu)處于與算例5.1 所示的外部環(huán)境與邊界條件,即上下外表面所處的環(huán)境溫度為T∞=900 K,換熱系數(shù)為h=1100 W·(m2·K)—1,流體與固體的導(dǎo)熱系數(shù)λ 以及流體域的入口溫度Tin與算例5.1 中的相同.圖14為在不同流速下整個冷卻系統(tǒng)的溫度云圖,圖15為上部外表面溫度變化曲線,而圖16為斜管流體域下表面的溫度變化曲線.

從圖14 可以看出,不同流速下冷卻系統(tǒng)中間部分的溫度變化很大.當流速很低時(如v=0.002 m/s的情況),只有流體的入口部分能得到冷卻;隨著流速的增加,冷卻效果越來越明顯;當流速為v=0.5 m/s 時,整個中間結(jié)構(gòu)都能得到冷卻.另外,從圖16 和圖14(d)可以看出,當流速為v=0.5 m/s時,斜管流體域下表面的溫度基本為流體的入口溫度,說明在此流速下不需要中間的冷卻管道對結(jié)構(gòu)進行冷卻.

圖14 不同流速下的溫度云圖 (a) v=0.002 m/s;(b) v=0.005 m/s;(c) v=0.01 m/s;(d) v=0.5 m/sFig.14.Contours under different velocities: (a) v=0.002 m/s;(b) v=0.005 m/s;(c) v=0.01 m/s;(d) v=0.5 m/s.

圖15 不同流速下結(jié)構(gòu)上部外表面溫度變化曲線Fig.15.Temperature variation curve on the outer surface of upper structure under different velocities.

圖16 不同流速下斜管流體域下表面溫度變化曲線Fig.16.Temperature variation curve on the lower surface of fluid domain of the oblique channel under different velocities.

5.2.2 流、固域采用不同條件的計算分析

接下來將整個結(jié)構(gòu)的上下表面對流換熱邊界條件分別設(shè)置為T∞=1000 K,h=800 W·(m2·K)—1和T∞=700 K,h=1100 W·(m2·K)—1;各計算域的導(dǎo)熱系數(shù)λ、流體域的入口溫度Tin和流速v如表1所示.圖17為計算得到的整個冷卻系統(tǒng)的溫度云圖,圖18 給出了冷卻系統(tǒng)上表面(用Upper-surface標記)、下表面(Lower-surface 標記)、和各區(qū)域界面上的溫度沿x方向的變化曲線,其中Interface 45 表示計算域4 與5的界面,其他類推.

表1 各計算域的材料參數(shù)與入口邊界條件Table 1.Material parameters and boundary conditions of each computational domain.

從圖17 可以看出,冷卻通道附近的溫度與冷卻流體的溫度接近,說明冷卻流體帶走的熱量明顯.另外,從圖18 可以看出,入口(x=0)處上下外表面的溫度與環(huán)境溫度具有較大的差別,這說明冷卻流體可以顯著地降低結(jié)構(gòu)的表面溫度,這是高溫裝備通常使用主動冷卻熱防護的主要原因.

圖17 計算得到的管道冷卻系統(tǒng)溫度云圖Fig.17.Contours of computed temperature over the channel cooling system.

圖18 計算得到的冷卻系統(tǒng)上下表面和各區(qū)域界面的溫度變化曲線Fig.18.Variation curve of the computed temperature on the upper and lower surfaces as well as on the interfaces of the cooling system.

5.3 三維三管道結(jié)構(gòu)流固域間耦合傳熱分析

本文分析的第3 個算例是含有3 根管道的三維冷卻結(jié)構(gòu),如圖19 所示.管道的尺寸均為1.333 mm×0.667 mm.上下表面的溫度分別給定為T1=600 K 和T2=350 K.計算中,整個冷卻結(jié)構(gòu)分為6 個計算域,其中Ω4—Ω6為流體域,其他為固體域;冷卻通道內(nèi)冷卻液的入口溫度Tin和流速v見表2,流體的密度與比熱均為ρ=1000 kg/m3與cp=4200 J·(kg·K)—1.

圖19 含三根管道的三維冷卻結(jié)構(gòu)尺寸與邊界條件Fig.19.Dimensions and B.C.of a 3D cooling structure with three channels.

表2 三維流固結(jié)構(gòu)各計算域的材料參數(shù)與邊界條件Table 2.Material parameters and boundary conditions of each computational domain.

圖20為FLM 計算用的網(wǎng)格圖,總共有172789個配置點,其中沿管道走向(x方向)劃分了31 個點,靠近壁面的網(wǎng)格間距比例因子為1.2;各域沿垂直方向(z方向)均劃分為11 個點,壁面比例因子為1.5;沿橫向(y方向),除上下蓋板的外伸部分為25 個點外,其他部分均為21 個點,壁面比例因子均為1.5.圖21為橫截面網(wǎng)格與局部放大圖,可以看出下蓋板的外伸部分是曲面,而且靠近右下角的網(wǎng)格極度非正交.

圖20 所有配置點的線系連成的網(wǎng)格Fig.20.Mesh connected by line sets of all collocation points.

圖21 橫截面上的配置點線系連成的網(wǎng)格與局部放大圖Fig.21.Mesh connected by line sets on the transverse section and a locally refined part around a corner.

圖22 是計算得到的總體結(jié)構(gòu)上的溫度云圖,而圖23 是通過三根通道中線的垂直面上的溫度云圖.圖24 和圖25 分別為沿圖22 所示的線段L1與L2上的溫度分布.為了比較,也用商業(yè)有限容積法軟件FLUENT 和有限單元法軟件COMSOL 對該問題進行了計算,使用的網(wǎng)格是達到網(wǎng)格收斂下的網(wǎng)格.對比可以看出,FLM 計算結(jié)果與其他兩種方法的結(jié)果吻合的很好,證明了本文方法在使用非規(guī)則網(wǎng)格方面的正確性,這是傳統(tǒng)的有限差分法難以實現(xiàn)的.另外,表3 列出了采用3 種方法計算該問題時所用的節(jié)點數(shù)與時間的比較.其中,COMSOL只列了一種網(wǎng)格的結(jié)果,因為網(wǎng)格稍微密一點時計算不收斂,試算發(fā)現(xiàn)這主要是因為流體與固體的導(dǎo)熱系數(shù)相差較大所致.從表3 可以看出,在同等精度下FLM 需用的網(wǎng)格數(shù)最少,消耗時間也較少;COMSOL 用的時間最多,而且在流體入口處附近與其他兩種方法的結(jié)果有較大的差別,穩(wěn)定性也差,這也許是有限單元法在處理流體相關(guān)問題中存在的固有弱點[20].此外,FLM 與FLUENT 用粗、細兩種網(wǎng)格算的結(jié)果基本一樣,說明二者有很好的穩(wěn)定性.

圖22 總體結(jié)構(gòu)上的溫度云圖Fig.22.Contour plot of computed temperature over the entire structure.

圖23 管道走向垂直中面上的溫度云圖Fig.23.Contour plot of computed temperature over the vertical middle plan.

圖24 沿圖22 中所示的線段 L1 上的溫度分布Fig.24.Computed temperature along line L1 marked in Fig.22.

圖25 沿圖22 中所示的線段 L2 上的溫度分布Fig.25.Computed temperature along line L2 marked in Fig.22.

表3 三種計算方法的節(jié)點數(shù)與計算時間比較Table 3.Comparison of total number of nodes and computational time for three methods.

6 討論與結(jié)論

論文提出了一種全新的數(shù)值計算方法—有限線法,并將其應(yīng)用于求解流體域與固體域之間的耦合傳熱問題,根據(jù)算法的理論與在冷卻通道結(jié)構(gòu)中的算例分析,可以得出如下的結(jié)論:

1)本文通過線段方向?qū)?shù)建立的求導(dǎo)公式可以用統(tǒng)一的格式計算多維問題的任意高階的空間偏導(dǎo)數(shù);

2)創(chuàng)建的有限線法是一種配點型計算方法,其用法與傳統(tǒng)的有限差分法一樣方便、靈活,而且能模擬曲的和不規(guī)則的幾何形狀;

3)有限線法可以直接由問題的控制微分方程和邊界條件建立系統(tǒng)方程組,不需要變分原理、能量原理及其他需要積分的數(shù)學(xué)操作;

4)本文雖然以流固耦合傳熱問題為研究背景,但導(dǎo)出的公式也適用于求解其他工程問題的偏微分方程邊值問題,是一種通用的數(shù)值計算方法;

5)因為有限線法是一種強形式的計算方法,通常需要采用高階線系計算,經(jīng)驗表明,采用四階精度的中心線系和二階精度的迎風線系可以獲得高精度和高效率的計算結(jié)果;

本文所述方法是關(guān)于空間偏導(dǎo)數(shù)的離散的方法,考慮的控制方程是不含時間項的定常問題,但所導(dǎo)出的公式也同樣適用于含有時間項的非定常問題,此時,時間項可采用成熟的方法進行離散[24].