三維球、柱坐標(biāo)系下導(dǎo)熱微分方程的離散求解

賈欣鑫，徐明海，胡國華，劉娟，路輝，梁卓

三維球、柱坐標(biāo)系下導(dǎo)熱微分方程的離散求解

賈欣鑫，徐明海，胡國華，劉娟，路輝，梁卓

(中國石油大學(xué)(華東)儲運(yùn)與建筑工程學(xué)院，山東青島266555)

根據(jù)柱坐標(biāo)系與球坐標(biāo)系的導(dǎo)熱微分方程式推出了導(dǎo)熱微分方程在球坐標(biāo)系、柱坐標(biāo)系上的三維高精度數(shù)值求解離散公式，并與解析解進(jìn)行對比，驗(yàn)證了該離散公式有較高的精確度。在球坐標(biāo)系下離散θ擴(kuò)散項(xiàng)時(shí)，運(yùn)用積分第一中值定理成功處理了復(fù)雜的θ擴(kuò)散項(xiàng)的離散系數(shù)。該離散格式為三維柱坐標(biāo)與球坐標(biāo)下導(dǎo)熱微分方程的數(shù)值求解提供了良好的借鑒作用。同時(shí)為導(dǎo)熱微分方程在工程計(jì)算中的應(yīng)用提供了精確的數(shù)值離散格式與理論依據(jù)。

導(dǎo)熱微分方程;球坐標(biāo)系;數(shù)值傳熱;積分第一中值定理

數(shù)值傳熱學(xué)在解決實(shí)際復(fù)雜傳熱問題中得到廣泛的應(yīng)用［1－3］，其對應(yīng)的導(dǎo)熱偏微分方程的離散一直都是解決數(shù)值傳熱問題的關(guān)鍵之一。在數(shù)值傳熱學(xué)中，二維柱坐標(biāo)與極坐標(biāo)導(dǎo)熱微分方程已得到了較好的推廣應(yīng)用。目前在離散與解決三維柱坐標(biāo)與球坐標(biāo)時(shí)，缺少較為明確的數(shù)值求解的離散格式［4－5］?，F(xiàn)今數(shù)值傳熱學(xué)對于柱體與球體的數(shù)值求解，只是通過簡化為徑向或者二維極坐標(biāo)的方式來解決，這給數(shù)值傳熱學(xué)在三維柱坐標(biāo)與球坐標(biāo)條件下的計(jì)算、推廣與應(yīng)用造成了諸多不便，因此合理準(zhǔn)確地得到三維圓柱體導(dǎo)熱、球體導(dǎo)熱偏微分方程顯得尤為重要。基于以上考慮本文以三維柱坐標(biāo)與球坐標(biāo)為基礎(chǔ)，從微分方程的數(shù)值解離散格式方面入手，推導(dǎo)出精度較高的柱坐標(biāo)與球坐標(biāo)導(dǎo)熱微分方程離散格式，并通過一維解析解對比驗(yàn)證其精確度。

1 導(dǎo)熱偏微分方程的提出

將導(dǎo)熱微元體置于直角坐標(biāo)系中，運(yùn)用能量守恒原理和傅里葉(Fourier)定律［6］，建立直角坐標(biāo)系下導(dǎo)熱微分方程

其中:λ為導(dǎo)熱系數(shù);c為導(dǎo)熱體熱容;S為內(nèi)熱源強(qiáng)度。

采用坐標(biāo)變換法分別得到圓柱坐標(biāo)系與球坐標(biāo)系中導(dǎo)熱微分方程

離散導(dǎo)熱微分方程離散的基本方法主要有2種:Taylor級數(shù)展開法和有限體積法［7－8］。為明確其導(dǎo)熱偏微分方程的物理概念及保證離散系數(shù)的意義，本文采用有限體積法，即控制容積法。

2 柱坐標(biāo)系下導(dǎo)熱微分方程的離散

首先給方程(2)兩邊同時(shí)乘以r，然后對偏微分方程兩邊同時(shí)在圖1所示的控制容積以及非穩(wěn)態(tài)的時(shí)間項(xiàng)中積分:

圖1 圓柱體控制容積

其中:源項(xiàng)中S表示成為未知量的線性函數(shù)［7］;SC為常數(shù)部分;SP表示S隨溫度T變化而變化的曲線在P點(diǎn)的斜率。

4)離散結(jié)果

整理以上結(jié)果可得:

3 球坐標(biāo)系下導(dǎo)熱微分方程的離散

對控制方程(3)兩邊同時(shí)乘以r2sin2θ，對導(dǎo)熱微分方程在圖2所示的控制容積在非穩(wěn)態(tài)時(shí)間項(xiàng)中積分得:

圖2 球體控制容積

4 實(shí)例驗(yàn)證

4.1 圓柱體導(dǎo)熱離散方程的驗(yàn)證及網(wǎng)格無關(guān)性分析

為驗(yàn)證其離散解數(shù)值解的準(zhǔn)確性，考察一個(gè)內(nèi)外徑分別為R1、R2的圓筒壁，圓筒的內(nèi)外兩側(cè)保持無量綱溫度T1、T2，其徑向?qū)峤馕鼋鉃?/p>

應(yīng)用本文離散方程，得到如圖3所示高精確度的數(shù)值解，誤差為8.86×10－3%。

圖3 圓柱體解析解與數(shù)值解

通過增加計(jì)算區(qū)域網(wǎng)格的的數(shù)目，計(jì)算其誤差的大小，進(jìn)行網(wǎng)格無關(guān)性的分析，誤差的計(jì)算方式為

式(9)中:error表示誤差;T表示解析解;Tnum表示本文數(shù)值解。

根據(jù)誤差分析網(wǎng)格無關(guān)解［11－12］得到如圖4所示的誤差與網(wǎng)格數(shù)的關(guān)系。通過數(shù)值計(jì)算得出:當(dāng)網(wǎng)格數(shù)增加到5 000時(shí)，計(jì)算誤差已經(jīng)控制在1%以內(nèi)，隨著網(wǎng)格的加密計(jì)算誤差逐漸減少。本文的解析解誤差為8.86×10－3%，選用計(jì)算網(wǎng)格數(shù)為37×20×16=11 840。

圖4 圓柱體網(wǎng)格數(shù)與誤差關(guān)系

4.2 球體導(dǎo)熱離散方程的驗(yàn)證

為驗(yàn)證其球殼離散解數(shù)值解的準(zhǔn)確性，選取內(nèi)外徑分別為R1、R2的球殼，球殼的內(nèi)外兩側(cè)保持恒溫T1、T2，驗(yàn)證其一維徑向?qū)釂栴}。解析解為

無熱源穩(wěn)態(tài)導(dǎo)熱圓柱邊界條件:

應(yīng)用本文離散方程，得到如圖5所示的數(shù)值解，誤差為0.53%。

圖5 球體解析解與數(shù)值解

圖6為誤差與網(wǎng)格數(shù)的關(guān)系。通過圖6可以看出:當(dāng)網(wǎng)格數(shù)增加到12 000時(shí)，計(jì)算誤差已經(jīng)控制在0.5%以內(nèi)，隨著網(wǎng)格的加密計(jì)算誤差逐漸減小。本文的解析解誤差為0.53%，選用網(wǎng)格數(shù)為47×20×12=11 280。

圖6 球體網(wǎng)格數(shù)與誤差關(guān)系

5 結(jié)束語

研究了圓柱體與球體的三維導(dǎo)熱微分方程，對其進(jìn)行有限容積法的高精度的數(shù)值計(jì)算離散格式推導(dǎo)。在離散球體的過程中，運(yùn)用積分第一中值定理從理論上處理了復(fù)雜θ擴(kuò)散項(xiàng)的離散系數(shù)。該離散格式為科研工作者進(jìn)行三維柱坐標(biāo)與球坐標(biāo)下導(dǎo)熱微分方程的數(shù)值求解提供了良好的借鑒。該離散格式的驗(yàn)證應(yīng)用FORTRAN語言編寫計(jì)算，運(yùn)行穩(wěn)定，在理論基礎(chǔ)上驗(yàn)證了解析解與數(shù)值解的誤差，將球體的誤差范圍控制在了0.5%以內(nèi)，為三維柱體與球體導(dǎo)熱偏微分方程的研究與工程應(yīng)用提供了高精度、可靠的數(shù)值計(jì)算離散格式。

［1］陶文銓.數(shù)值傳熱學(xué)［M］.2版.西安:西安交通大學(xué)出版社，2001.

［2］王輝，徐明海.邊界界面導(dǎo)熱系數(shù)插值格式的改進(jìn)［J］.燕山大學(xué)學(xué)報(bào):自然科學(xué)版，2010，34(5):411－417.

［3］黃善波，李兆敏.冪律流體同心環(huán)空內(nèi)層流脈動(dòng)流的數(shù)值分析［J］.中國石油大學(xué)學(xué)報(bào):自然科學(xué)版，2011，35(2):127－136.

［4］金有海，姬廣勤，曹晴云，等.旋風(fēng)分離器排氣管內(nèi)氣相流場的數(shù)值模擬［J］.中國石油大學(xué)學(xué)報(bào):自然科學(xué)版，2008，32(6):109－112.

［5］DING P，WU XH，HE YL，et al.A fast and efficient method for predicting fluid flow and heat transfer problems［J］.ASM E Journal of H eat Transfer，2008，130:1－17.

［6］鄧先琛.對建立球面坐標(biāo)系中的導(dǎo)熱微分方程式探討［J］.安徽工學(xué)院學(xué)報(bào)，1987，6(4):58－70.

［7］PATANKA S V，SPALDING D B.A calculation procedure for heat，mass and momentum transfer in three－dimensional parabolic flow［J］.International Journal Heat Mass Transfer，1972，15:1787－1806.

［8］帕坦卡SV.傳熱與流動(dòng)的數(shù)值計(jì)算［M］.張政，譯.北京:科學(xué)出版社，1984:40.

［9］華東師范大學(xué)數(shù)學(xué)系編數(shù)學(xué)分析(上)北京:高等教育出版社，2001:218－219.

［10］同登科，周生田，吳愛弟，等.計(jì)算方法［M］.東營:石油大學(xué)出版社，2004:59－60.

［11］徐明海，張儼彬，陶文銓.一種改進(jìn)的Delaunay三角形化剖分方法［J］.石油大學(xué)學(xué)報(bào):自然科學(xué)版，2001，25 (2):100－105.

［12］LEE C K，HOBBSR E.Automatic adaptive finite element mesh generation over arbitrary two－dimensional domain using advancing front technique［J］.Computers＆Structure，1999，71:9－34.

(責(zé)任編輯 劉舸)

Numerical M ethod for Three-dimensional Heat Conduction in
Cylindrical and Spherical Coordinates

JIA Xin－xin，XU Ming－h(huán)ai，HU Guo－h(huán)ua，LIU Juan，LU Hui，LIANG Zhuo

(College of Pipeline and Civil Engineering，China University of Petroleum，Qingdao 266555，China)

According to the differential equations of heat conduction on cylindricaland spherical coordinate system，numerical solution of the discrete formulation on cylindrical and spherical coordinate system with high accuracy has been derived.Compared with the analytical solution，this discrete formula has been verified with a high degree of accuracy.Tomake the complex dispersion coefficient of diffusion termθmore concrete in spherical coordinates，this paper derived the discretion coefficientof diffusion termθby the firstmean value theorem of integral.The accurate schemes provides a good reference for researcherswhom work in solving the equation of heat conduction of three－dimensional cylindrical coordinates and spherical coordinates，and it will provide accurate numerical schemes and the theoretical basis for solving practical engineering problems.

differential equation of heat conduction;spherical coordinate;numerical heat transfer; the firstmean value theorem for integrals

TK121

A

1674－8425(2014)01－0033－05

10.3969/j.issn.1674－8425(z).2014.01.007