楊 瀾，蹇開林，張 亮

(重慶大學(xué) 航空航天學(xué)院，重慶 400044)

熱力耦合現(xiàn)象存在于眾多的高新技術(shù)中，諸如航空發(fā)動(dòng)機(jī)、汽車發(fā)動(dòng)機(jī)和電子器件等一系列裝置，特別是結(jié)構(gòu)界面接觸傳熱的研究。在實(shí)際情況中，結(jié)構(gòu)部件的加工表面不可能是絕對平整光滑的，而是由許多微觀的不規(guī)則的凸峰和凹谷組成的粗糙表面，而邊界面的凹凸形狀是直接影響耦合分析的核心因素。將熱力耦合分析與計(jì)算力學(xué)相結(jié)合的數(shù)值模擬方法已經(jīng)成為重要的研究方向，尋求合適的計(jì)算方法以及如何優(yōu)化精確溫度場和位移場分布依舊是熱點(diǎn)問題。正因如此，結(jié)構(gòu)界面接觸傳熱的研究正逐步從簡單走向復(fù)雜、從宏觀走向微觀、從表象走向機(jī)理、從定性走向定量，其建模方法和分析手段面臨著巨大的挑戰(zhàn)。

由材料的本構(gòu)關(guān)系、熱力學(xué)基本定律以及自由能表達(dá)推導(dǎo)的熱傳導(dǎo)方程式可以知道溫度場的分布受位移場的影響[1-8]，而位移場的分布也受溫度的影響，這是典型的耦合問題，需將位移場和溫度場二者耦合求解；在針對耦合熱應(yīng)力問題求解方面，Nowinski[4]提出了耦合作用是將機(jī)械能轉(zhuǎn)化成熱能，繼而耗散掉；Biot等[9]提出將熱容量作為一個(gè)新的狀態(tài)變量引入虛耗散的熱力學(xué)原理中，由此導(dǎo)出的場方程將具有積分-微分形式，避免了因考慮熱源而需引入的另一變量熵；Nariboli[10]用Laplace變換法研究了球形空腔壁溫度突然有跳躍的熱彈性動(dòng)力學(xué)問題；Takeuti等[11]用熱彈性位勢和Love位移函數(shù)求解了軸對稱有限圓柱動(dòng)態(tài)耦合熱應(yīng)力問題；Oden[12]用有限單元法計(jì)算了各種物性系數(shù)下的半空間耦合問題并與精確解進(jìn)行比較；王林翔等[13]用微分代數(shù)的方法求解了耦合熱彈性動(dòng)力學(xué)問題；Gao[14]用邊界元法進(jìn)行熱彈性動(dòng)力學(xué)問題的求解。

針對于熱力耦合現(xiàn)象，求解的問題分為兩種，一種是動(dòng)態(tài)的熱力耦合問題，即耦合熱彈性動(dòng)力學(xué)問題，另一種是準(zhǔn)靜態(tài)的熱力耦合問題，即非定常擬靜態(tài)的熱力耦合問題。區(qū)分二者的關(guān)鍵在于外部機(jī)械載荷或者熱載荷隨時(shí)間的變化是否劇烈，也就是載荷是否具備沖擊的性質(zhì)，當(dāng)載荷不具備沖擊性質(zhì)時(shí)，位移產(chǎn)生的加速度項(xiàng)影響很小可忽略，這種并非靜態(tài)問題但可忽略加速度影響的計(jì)算方法，稱其為準(zhǔn)靜態(tài)的計(jì)算方法，大多數(shù)的工程實(shí)踐問題都屬于準(zhǔn)靜態(tài)的范疇，這樣一來可以極大地簡化計(jì)算過程。有限元方法是求解熱力耦合問題的主要方法，但計(jì)算精度受網(wǎng)格影響較大，在處理龐大的工程結(jié)構(gòu)時(shí)，網(wǎng)格劃分是一項(xiàng)巨大的任務(wù)，并且基于靜態(tài)形函數(shù)的變形描述很難充分逼近彈性體在熱力耦合作用下的真實(shí)位移場，因此近年來有不少學(xué)者提出了使用無網(wǎng)格方法解決耦合熱彈性問題[15]。無網(wǎng)格方法[16-18]擺脫了網(wǎng)格的束縛，在計(jì)算域內(nèi)使用節(jié)點(diǎn)分布去代替網(wǎng)格作用，計(jì)算點(diǎn)的形函數(shù)隨著計(jì)算點(diǎn)變化而變化，即使在計(jì)算節(jié)點(diǎn)分布不規(guī)則時(shí)，計(jì)算精度的損失也較小；并且由于節(jié)點(diǎn)生成較為方便，因而自適應(yīng)性能優(yōu)越，將自適應(yīng)技術(shù)引入無網(wǎng)格方法，能夠合理有效地增加節(jié)點(diǎn)密度，繼而保證計(jì)算精度，避免因?yàn)楣?jié)點(diǎn)數(shù)量和位置的不合理造成的精度失準(zhǔn)，在計(jì)算凹凸邊界時(shí)，二者的結(jié)合能夠克服凹凸邊界網(wǎng)格難以劃分問題，同時(shí)又能保證計(jì)算精度的不失準(zhǔn)。筆者以無網(wǎng)格伽遼金算法(EFG)計(jì)算了二維混合邊界條件下的非定常擬靜態(tài)熱力耦合問題，以移動(dòng)最小二乘法(MLS)構(gòu)建了形函數(shù)，應(yīng)用后驗(yàn)誤差公式結(jié)合采用H型自適應(yīng)計(jì)算，提出了適用于非定常溫度場擬靜態(tài)熱力耦合問題的EFG法計(jì)算模型，最后以數(shù)值算例驗(yàn)證了本文方法的可行性。

1 移動(dòng)最小二乘法(MLS)

考慮平面域Ω的某一場變量u(x)(x=((x，y))表示在平面域Ω內(nèi)的某一點(diǎn)的坐標(biāo))的未知標(biāo)量函數(shù)，則定義在點(diǎn)x的MLS近似表達(dá)式為

(1)

式中(p(x))是關(guān)于二維空間坐標(biāo)x=(x,y)的基函數(shù)。此基函數(shù)通常為利用Pascal三角形所決定的單項(xiàng)式以確保其最小完備性。m表示為基函數(shù)完備項(xiàng)的個(gè)數(shù)。

在二維空間中單項(xiàng)式基函數(shù)為：

線性基

m=3 (p(x))T=[1,x,y]，

(2)

二次基

m=6 (p(x))T=[1,x,y,x2,xy,y2],

(3)

(a(x))=[a1(x)a2(x) …am(x)] ，

(4)

式中a(x)為待定參數(shù)向量。

現(xiàn)將二維求解域Ω用N個(gè)節(jié)點(diǎn)進(jìn)行離散，每個(gè)節(jié)點(diǎn)的節(jié)點(diǎn)值uI已知，在每個(gè)節(jié)點(diǎn)xI(I=1,2,…,N)處定義一個(gè)權(quán)函數(shù)ωI(x)，并且該函數(shù)是緊支的，即在支撐域ΩI中權(quán)函數(shù)大于零，而在支撐域外面權(quán)函數(shù)等于零，支撐域形狀常見有矩形和圓形，如圖1所示。

為計(jì)算方便，這里采用矩形支撐域，以及3次樣條插值權(quán)函數(shù)。當(dāng)支撐域?yàn)榫匦螘r(shí)，

ωI(x)=ωI(rx)ωI(ry)，

(5)

(6)

式中：L，D分別為矩形支撐域的長和寬。為求a(x)，定義泛函J:

(7)

[A(x)]m×m(a(x))m×1=[B(x)]m×N[U]N×1，

(8)

其中

B(x)=[ω1(x)(px1)ω2(x)(px2)ω3(x)(px3)…ωN(x)(pxN)]，

(9)

U=[u1u2u3…uN]T是域Ω中所有節(jié)點(diǎn)的節(jié)點(diǎn)值所形成的向量。

繼而可得

uh(x)=(p(x))T[A(x)]-1B(x)U=N(x)U，

(10)

式中：N(x)=[N1(x)N2(x) …NN(x)]，類似于有限單元法中的形函數(shù)，不同之處在于無網(wǎng)法中的形函數(shù)是關(guān)于全求解域，而有限單元法中的形函數(shù)是關(guān)于單元區(qū)域。

2 耦合熱應(yīng)力問題的EFG方法

2.1 耦合熱應(yīng)力控制方程

設(shè)線性耦合熱彈性力學(xué)的求解區(qū)域?yàn)槎S平面Ω，Γ為求解平面的邊界。熱彈性力學(xué)中的應(yīng)力應(yīng)變關(guān)系為

σij=2μεij+λεkkδij-βθδij,

(11)

式中：σij是應(yīng)力;μ,λ代表拉梅常數(shù);εij代表應(yīng)變;β代表熱應(yīng)力系數(shù)；θ代表物體相對于初始溫度T0的溫差(θ=T-T0)。

在小變形條件下：

εij=(ui,j+uj,i)/2。

(12)

在處理耦合熱應(yīng)力問題時(shí)，采用分區(qū)變分原理，即溫度場和位移場分別采用變分離散，再并行求解，經(jīng)典熱彈性力學(xué)理論中溫度場和位移場的控制方程及邊界條件如下：

σij,j+bi=0，

(13)

(14)

(15)

(16)

(17)

2.2 耦合熱應(yīng)力控制方程的弱形式及離散化

在二維求解區(qū)域內(nèi)采用N個(gè)節(jié)點(diǎn)進(jìn)行離散，在整體區(qū)域內(nèi)采用加權(quán)余量法弱式，以溫度T的變分δT，及位移U的變分δU作為試函數(shù)進(jìn)行離散，可得式(12)的離散形式：

(18)

運(yùn)用分部積分和高斯散度定理，并施加邊界條件可得：

(19)

因?yàn)橐苿?dòng)最小二乘法(MLS)構(gòu)建的形函數(shù)并不滿足Kronecker delta特性，因此在式(18)的基礎(chǔ)上使用拉格朗日乘子法對求解域進(jìn)行本質(zhì)邊界條件(第一類邊界條件)施加。

(20)

式(18)中λ為拉格朗日乘子，其視為坐標(biāo)的位置函數(shù)，可利用本質(zhì)邊界條件上的節(jié)點(diǎn)進(jìn)行插值。

(21)

(22)

λ=φΛ，

(23)

式中：φi為拉格朗日插值函數(shù)；λi為本質(zhì)邊界條件上用于插值的節(jié)點(diǎn)值；nλ為本質(zhì)邊界上用于拉格朗日插值的節(jié)點(diǎn)數(shù)目。

同理運(yùn)用上述的變分過程外加在本質(zhì)邊界條件上施加拉格朗日乘子法，對位移場的控制方程(13)進(jìn)行離散處理可得到：

(24)

位移場近似函數(shù)以及溫度場的近似函數(shù)分別為

(25)

(26)

將式(25)和(26)代入式(24)和式(20)中可得到如下方程式：

(27)

式中：

(28)

(29)

(30)

(31)

(32)

(33)

(34)

(35)

(36)

(37)

(38)

(39)

(40)

(41)

(42)

(43)

式中：E,ν分別代表平面應(yīng)力條件下的彈性模量和泊松比，α是熱膨脹系數(shù)。

式(27)為施加邊界條件以及拉格朗日乘子法后的系統(tǒng)一階偏微分方程組，采用向后差分的形式對上述微分方程組進(jìn)行求解可得如下方程：

(44)

2.3 自適應(yīng)

自適應(yīng)分析的核心思想是根據(jù)已有的數(shù)值計(jì)算模型判斷出該計(jì)算區(qū)域內(nèi)的高誤差分布區(qū)域，然后進(jìn)行算法調(diào)整，繼而改進(jìn)數(shù)值計(jì)算結(jié)果精度。常用的自適應(yīng)方法有R型自適應(yīng)，H型自適應(yīng)，以及HP型自適應(yīng)。這里采用H型自適應(yīng)，H型自適應(yīng)在無網(wǎng)格方法中的本質(zhì)是節(jié)點(diǎn)的增加和刪除。它的實(shí)現(xiàn)主要包含以下兩個(gè)方面，一方面是后驗(yàn)誤差估計(jì)，另一方面就是節(jié)點(diǎn)的插入策略。通過引入局部誤差計(jì)算公式(45)和全局誤差計(jì)算公式(46)，利用Voronoi鄰接準(zhǔn)則去進(jìn)行節(jié)點(diǎn)的添加，分別以設(shè)定的局部誤差限εu和全局誤差限Lu作為局部節(jié)點(diǎn)的添加界限和自適應(yīng)過程終止界限。關(guān)于自適應(yīng)詳細(xì)的后驗(yàn)誤差計(jì)算公式的推導(dǎo)過程，鑒于篇幅的原因，此處不再贅述，可參考文獻(xiàn)[19-20]。有如下誤差估計(jì)式。

局部誤差估計(jì)式：

(45)

全局誤差估計(jì)式：

(46)

對公式中的參量做如下說明：設(shè)定X1=(x1,y1)，X2=(x2,y2)是平面中對應(yīng)的兩點(diǎn)，u1，u2分別是以上2個(gè)點(diǎn)在求解域內(nèi)對應(yīng)的場函數(shù)值，?代表梯度，Δx=(x2-x1,y2-y1)表示兩節(jié)點(diǎn)構(gòu)成的向量，uh代表利用EFG求得的場函數(shù)近似值。

3 數(shù)值算例

3.1 算例1

圖2 混合邊界條件下的光滑板Fig. 2 Smooth sliding plate with mixed boundary conditions

圖3 方板布點(diǎn)圖Fig. 3 Node distribution of square plate

圖4 方板網(wǎng)格圖Fig. 4 Mesh distribution of square plate

圖5 節(jié)點(diǎn)231號的時(shí)間位移歷程圖Fig. 5 Time-displacement history of node 231

圖6 節(jié)點(diǎn)306號的時(shí)間位移歷程圖Fig. 6 Time-displacement history of node 306

圖7 節(jié)點(diǎn)231號的時(shí)間溫度歷程圖Fig. 7 Time-temperature history of node 231

圖8 節(jié)點(diǎn)306號的時(shí)間溫度歷程圖Fig. 8 Time-temperature history of node 306

控制方程：

(47)

σij,i+bi=0。

(48)

邊界條件：

(49)

(50)

ux=0=0，ux=L=0，uy=L=0，

(51)

k=1，ξ2=1，ξ4=2，ρ=1，c=1，

(52)

TA=0，E=1，ν=0.3，α=0.02。

(53)

從圖5～8可以看出，無網(wǎng)格方法計(jì)算溫度值和位移值的精度較有限元法大致相同，但在時(shí)間步為0.10 s時(shí)，EFG算法的計(jì)算結(jié)果更加靠近FEM算法下0.01 s步長時(shí)的計(jì)算結(jié)果，因而無網(wǎng)格方法的收斂速度較有限元收斂速度略快。

圖9 混合邊界條件下的凹凸板Fig. 9 Concave convex plate with mixed boundary conditions

3.2 算例2

此算例是凹凸平板在混合邊界條件下的熱力耦合計(jì)算問題，如圖9所示，長度L=1 m，寬度D=1 m，厚度為1 mm，圖10和11為布點(diǎn)示意圖和網(wǎng)格示意圖。計(jì)算參數(shù)按圖示給定，左右兩端及下端施加位移固定約束，在最左端邊界是絕熱邊界條件，最上端和最右端存在對流熱交換，對流熱交換系數(shù)分別是ζ2=1，ζ4=2，最底部邊界是第一類傳熱邊界條件溫度恒為0 ℃，初始溫度場分布是T0=50 ℃，分析步長度Δt=0.001 s，共計(jì)迭代1 000步，為了更好地對比體現(xiàn)無網(wǎng)格伽遼金法的計(jì)算精度，采用在相同的時(shí)間步下進(jìn)行計(jì)算，有限元算法采用四邊形八節(jié)點(diǎn)二次單元計(jì)算，輸出頂端凹凸邊界的溫度場及位移場分布，分別截取1 s的計(jì)算結(jié)果，計(jì)算結(jié)果如圖12、13所示。

圖11 凹凸板網(wǎng)格圖Fig. 11 Concave convex plate mesh distribution

圖12 t=1 s凹凸邊界溫度場分布Fig. 12 Temperature field distribution of concave convex boundary in t=1 s

圖13 t=1 s凹凸邊界位移場分布Fig. 13 Displacement field distribution of concave convex boundary in t=1 s

控制方程：

(54)

σij,i+bi=0。

(55)

邊界條件：

(56)

(57)

ux=0=0，ux=L=0，uy=L=0。

(58)

k=1，ξ2=1，ξ4=2，ρ=1，c=1，

(59)

TA=0，E=1，ν=0.3，α=0.02。

(60)

由圖12、13可看出，在t=1 s時(shí)，兩種不同的算法下，溫度場和位移場差異較大，使用有限元網(wǎng)格節(jié)點(diǎn)替代無網(wǎng)格節(jié)點(diǎn)并不適用，由移動(dòng)最小二乘法(MLS)的推導(dǎo)過程，可以知道造成數(shù)據(jù)偏差較大的原因是無網(wǎng)格伽遼金方法的計(jì)算精度受節(jié)點(diǎn)密度、節(jié)點(diǎn)擺放位置以及權(quán)函數(shù)多重因素的影響，所以無規(guī)律地增加節(jié)點(diǎn)密度會影響無網(wǎng)格伽遼金方法計(jì)算精度，因此后續(xù)應(yīng)采用自適應(yīng)技術(shù)與之相結(jié)合，使用自適應(yīng)技術(shù)控制節(jié)點(diǎn)密度合理增加。

3.3 算例3

圖14 自適應(yīng)節(jié)點(diǎn)分布過程Fig. 14 Adaptive node distribution process

從圖15～20的計(jì)算結(jié)果分析，在局部放大圖中，可以看到隨著自適應(yīng)過程的增加，對于節(jié)點(diǎn)11和節(jié)點(diǎn)15而言，不論是溫度值還是位移值都隨著時(shí)間的變化逐漸趨于穩(wěn)定并且是收斂的，且t=1 s時(shí)整體的溫度場和位移場的分布也都是穩(wěn)定且收斂的。結(jié)果表明將后驗(yàn)誤差式引入自適應(yīng)調(diào)節(jié)之中，能夠合理有效地控制節(jié)點(diǎn)的增加，繼而保證計(jì)算精度，自適應(yīng)技術(shù)與無網(wǎng)格方法相結(jié)合能夠應(yīng)對具有凹凸邊界作用下熱力耦合問題，為后續(xù)無網(wǎng)格方法的程序化、軟件化及商業(yè)化提供了十分良好的案例。

圖15 第2、第3次自適應(yīng)溫度場分布(t=1 s)Fig. 15 Second and third adaptive temperature field distribution t=1 s

圖16 第2、第3次自適應(yīng)位移場分布(t=1 s)Fig. 16 Second and third adaptive displacement field distribution t=1 s

圖17 節(jié)點(diǎn)11號的時(shí)間溫度歷程圖Fig. 17 Time-temperature history of node 11

圖18 節(jié)點(diǎn)15號的時(shí)間溫度歷程圖Fig. 18 Time-temperature history of node 15

圖19 節(jié)點(diǎn)11號的時(shí)間位移歷程圖Fig. 19 Time-displacement history of node 11

圖20 節(jié)點(diǎn)15號的時(shí)間位移歷程圖Fig. 20 Time-displacement history of node 15

4 結(jié) 論

采用非定常溫度場擬靜態(tài)熱力耦合問題的EFG法計(jì)算模型，分別求解了光滑平板，凹凸邊界平板的熱力耦合問題的位移場和溫度場分布，并且采用H型自適應(yīng)對實(shí)際工程中的凹凸邊界熱力耦合問題進(jìn)行了算法優(yōu)化，研究表明：

1)在應(yīng)對光滑邊界的熱力耦合問題時(shí)，無網(wǎng)格方法的精度以及收斂速度略高于有限元方法，并且由于節(jié)點(diǎn)施加的便捷性，可結(jié)合有限元網(wǎng)格剖分軟件的單元節(jié)點(diǎn)去替代無網(wǎng)格節(jié)點(diǎn)，在實(shí)際工程中，可優(yōu)先考慮使用無網(wǎng)格方法對結(jié)構(gòu)溫度場和位移場隨時(shí)間的變化進(jìn)行計(jì)算，能夠獲得十分精確的數(shù)值解。

2)在應(yīng)對結(jié)構(gòu)具有凹凸邊界的熱力耦合問題時(shí)，將無網(wǎng)格方法與H型自適應(yīng)技術(shù)相結(jié)合，在確保邊界條件順利施加的情況下又能夠合理地增加節(jié)點(diǎn)密度，繼而保證計(jì)算精度，考慮到計(jì)算時(shí)間成本的問題，可用相同邊界條件下定常溫度場靜態(tài)問題的自適應(yīng)節(jié)點(diǎn)分布圖去替代計(jì)算。自適應(yīng)技術(shù)與無網(wǎng)格算法相結(jié)合，為后續(xù)無網(wǎng)格算法的程序化、軟件化及商業(yè)化提供了十分廣闊的前景。