基于無(wú)網(wǎng)格自然鄰接點(diǎn)Petrov－Galerkin法求解帶源參數(shù)瞬態(tài)熱傳導(dǎo)問(wèn)題*

李慶華 陳莘莘

(華東交通大學(xué)土木建筑學(xué)院，南昌 330013)

基于無(wú)網(wǎng)格自然鄰接點(diǎn)Petrov－Galerkin法求解帶源參數(shù)瞬態(tài)熱傳導(dǎo)問(wèn)題*

李慶華 陳莘莘?

(華東交通大學(xué)土木建筑學(xué)院，南昌 330013)

基于無(wú)網(wǎng)格自然鄰接點(diǎn)Petrov－Galerkin法，本文建立了一種求解帶源參數(shù)瞬態(tài)熱傳導(dǎo)問(wèn)題的新方法．為了克服移動(dòng)最小二乘近似難以準(zhǔn)確施加本質(zhì)邊界條件的缺點(diǎn)，采用了自然鄰接點(diǎn)插值構(gòu)造試函數(shù)．在局部多邊形子域上采用局部Petrov－Galerkin方法建立瞬態(tài)熱傳導(dǎo)問(wèn)題的積分弱形式．這些多邊形子域可由Delaunay三角形創(chuàng)建．時(shí)間域則通過(guò)傳統(tǒng)的兩點(diǎn)差分法進(jìn)行離散．最后通過(guò)算例驗(yàn)證了該數(shù)值算法的有效性和正確性．

熱傳導(dǎo)問(wèn)題， 源參數(shù)， 無(wú)網(wǎng)格法， 局部Petrov－Galerkin法， 自然鄰接點(diǎn)插值

引言

隨著時(shí)間變化，不同場(chǎng)點(diǎn)上的溫度場(chǎng)隨時(shí)間而變化，這種傳熱問(wèn)題被看作瞬態(tài)熱傳導(dǎo)問(wèn)題．瞬態(tài)熱傳導(dǎo)問(wèn)題是實(shí)際工程中十分普遍的現(xiàn)象，很多學(xué)者都對(duì)其進(jìn)行了分析．由于結(jié)構(gòu)的形狀以及變溫條件的復(fù)雜性，要精確地確定溫度場(chǎng)主要依賴于數(shù)值模擬的方法［1－2］．近年來(lái)，隨著無(wú)網(wǎng)格方法的迅速發(fā)展，無(wú)網(wǎng)格方法正以其獨(dú)特優(yōu)勢(shì)在熱傳導(dǎo)問(wèn)題的研究中得到了廣泛的應(yīng)用［3－7］．

無(wú)網(wǎng)格自然鄰接點(diǎn) Petrov－Galerkin 法［8，9］本質(zhì)上是基于自然鄰接點(diǎn)插值的無(wú)網(wǎng)格局部Petrov－Galerkin(Meshless Local Petrov－Galerkin，MLPG)法．這種方法不僅具有自然元法易于施加本質(zhì)邊界條件的優(yōu)點(diǎn)，而且還具有MLPG法的一些優(yōu)良特性．在這種方法中，任意節(jié)點(diǎn)的子域是由以該節(jié)點(diǎn)為公共頂點(diǎn)的Delaunay三角形構(gòu)成的多邊形區(qū)域．此外，采用有限元法的三節(jié)點(diǎn)三角形單元的形函數(shù)作為權(quán)函數(shù)，可以降低域積分中被積函數(shù)的階次，提高了計(jì)算效率．本文嘗試將無(wú)網(wǎng)格自然鄰接點(diǎn)Petrov－Galerkin法應(yīng)用于帶源參數(shù)瞬態(tài)熱傳導(dǎo)問(wèn)題的求解計(jì)算．首先通過(guò)變換將此類帶源參數(shù)的熱傳導(dǎo)問(wèn)題轉(zhuǎn)化為標(biāo)準(zhǔn)的熱傳導(dǎo)控制方程形式，然后基于無(wú)網(wǎng)格自然鄰接點(diǎn)Petrov－Galerkin法理論詳細(xì)推導(dǎo)了相應(yīng)的計(jì)算公式．最后，典型算例的計(jì)算和對(duì)比分析驗(yàn)證了應(yīng)用無(wú)網(wǎng)格自然鄰接點(diǎn)Petrov－Galerkin法分析帶源參數(shù)瞬態(tài)熱傳導(dǎo)問(wèn)題的有效性．

1 自然鄰接點(diǎn)插值

自然鄰接點(diǎn)插值是基于著名的Voronoi結(jié)構(gòu)和Delaunay三角形網(wǎng)格．按插值基函數(shù)的不同，自然鄰接點(diǎn)插值可分為 Sibson 插值［8－10］和 Laplace插值［11－14］．本文采用 Sibson 插值．

一階Voronoi結(jié)構(gòu)定義:到節(jié)點(diǎn)xI的距離小于到其它任何節(jié)點(diǎn)xJ的距離的區(qū)域的集合．其數(shù)學(xué)表達(dá)為

式中，d(x，xI)為點(diǎn)x與節(jié)點(diǎn)xI的距離．此外，二階Voronoi結(jié)構(gòu)TIJ的數(shù)學(xué)表述為:

在幾何意義上，TIJ是那些以xI為最近點(diǎn)，以xJ為第二近點(diǎn)的空間點(diǎn)的位置集合．

當(dāng)采用Sibson插值時(shí)，點(diǎn)x對(duì)節(jié)點(diǎn)xI的形函數(shù)φI(x)定義為

式中，AI(x)為點(diǎn)x與節(jié)點(diǎn)xI的二階Voronoi結(jié)構(gòu)TxI的面積．

定義了各節(jié)點(diǎn)的插值函數(shù)后，點(diǎn)x的溫度類似于有限元法可寫為

式中，uI(I=1，2，…，n)是點(diǎn)x周圍自然鄰節(jié)點(diǎn)I的溫度，φI(x)為對(duì)應(yīng)節(jié)點(diǎn)的形函數(shù)．

2 控制方程及其局部弱形式

考慮如下的帶源參數(shù)p(t)的二維瞬態(tài)熱傳導(dǎo)問(wèn)題

初始條件為

邊界條件為

式中，p，φ，f，g0，g1，h0，h1為已知函數(shù);u是未知函數(shù)．

為了將帶源參數(shù)p(t)的熱傳導(dǎo)問(wèn)題轉(zhuǎn)化為標(biāo)準(zhǔn)的熱傳導(dǎo)控制方程形式，作變換［5-7］

方程(5)變?yōu)?/p>

初始條件為

邊界條件為

經(jīng)過(guò)這個(gè)變換，源參數(shù)p(t)消失了，我們可以把式(13)～(18)當(dāng)成標(biāo)準(zhǔn)的瞬態(tài)熱傳導(dǎo)問(wèn)題進(jìn)行求解．

將問(wèn)題域及其邊界用N個(gè)節(jié)點(diǎn)離散，在任意節(jié)I點(diǎn)對(duì)應(yīng)的子域ΩIs上，平衡方程式(13)的局部弱形式可以寫為

式中，vI為加權(quán)函數(shù)．

對(duì)方程(19)的左邊項(xiàng)進(jìn)行分部積分并利用散度定理后，可得

式中，Q=r(t)φ(x，y，t)．在無(wú)網(wǎng)格自然鄰接點(diǎn)Petrov－Galerkin 法中，子域是由共享節(jié)點(diǎn)I的Delaunay三角形TIi構(gòu)成，如圖1所示．為了簡(jiǎn)化計(jì)算，采用有限元法的三節(jié)點(diǎn)三角形單元的形函數(shù)作為權(quán)函數(shù)，并注意到式(15)～(18)中無(wú)自然邊界條件，方程(20)可改寫為:

圖1 局部多邊形子域Fig．1 The local polygonal sub － domains

由于只對(duì)空間域進(jìn)行離散，求解域Ω內(nèi)的試函數(shù)w(x，t)可由式(4)表示為

將式(22)代入式(21)可得如下的離散方程

方程(23)是一組以時(shí)間為獨(dú)立變量的線性常微分方程組．本文采用傳統(tǒng)的兩點(diǎn)差分法求解方程(23)．對(duì)給定的時(shí)間步長(zhǎng)Δt，取插值公式

式中，θ是可以自由選擇的時(shí)間加權(quán)系數(shù)，不同的取值對(duì)應(yīng)不同的時(shí)間差分格式．本文取θ=2/3．

3 算例分析

為了驗(yàn)證本文提出的求解帶源參數(shù)瞬態(tài)熱傳導(dǎo)問(wèn)題數(shù)值方法的有效性，本文給出了具體的算例．

考慮問(wèn)題(5)～(10)，其中

圖2 節(jié)點(diǎn)布置Fig．2 Node distribution

圖 3 u(x，y，t)的數(shù)值解和解析解Fig．3 The numerical and exact solution of u(x，y，t)

圖 4 u(x，y，t)的數(shù)值解和解析解Fig．4 The numerical and exact solution of(x，y，t)

我們利用無(wú)網(wǎng)格自然鄰接點(diǎn) Petrov－Galerkin法對(duì)上述帶源參數(shù)瞬態(tài)熱傳導(dǎo)問(wèn)題進(jìn)行求解．如圖2所示，在求解域內(nèi)均勻布置了121個(gè)節(jié)點(diǎn)．進(jìn)行數(shù)值積分時(shí)，在構(gòu)成子域的每個(gè)Delaunay三角形上布置3個(gè)高斯點(diǎn)．時(shí)間步長(zhǎng)取為Δt=0．001．圖3給出了函數(shù)u(x，y，t)在t=0．1，0．5 和 0．9 時(shí)刻y=0．6處的解析解和數(shù)值解．圖4給出了函數(shù)u(x，y，t)在t=0．1，0．5 和0．9 時(shí)刻x=0．5 處的解析解和數(shù)值解．通過(guò)圖3和圖4可以看出，本文提出的無(wú)網(wǎng)格自然鄰接點(diǎn)Petrov－Galerkin法計(jì)算得到的數(shù)值解和解析解吻合很好．

4 結(jié)論

無(wú)網(wǎng)格自然鄰接點(diǎn)Petrov－Galerkin法作為一種無(wú)網(wǎng)格方法，不需要對(duì)求解域進(jìn)行網(wǎng)格劃分，前處理簡(jiǎn)單．由自然鄰接點(diǎn)插值創(chuàng)建的形函數(shù)滿足Kronecker delta函數(shù)性質(zhì)，從而能夠直接準(zhǔn)確地施加本質(zhì)邊界條件．此外，采用有限元法的三節(jié)點(diǎn)三角形單元的形函數(shù)作為權(quán)函數(shù)，可以降低域積分中被積函數(shù)的階次，提高了計(jì)算效率．本文建立了一種基于無(wú)網(wǎng)格自然鄰接點(diǎn)Petrov－Galerkin法的帶源參數(shù)瞬態(tài)熱傳導(dǎo)問(wèn)題分析方法，并推導(dǎo)了相應(yīng)的計(jì)算公式．?dāng)?shù)值算例表明，用無(wú)網(wǎng)格自然鄰接點(diǎn)Petrov－Galerkin法求解帶源參數(shù)瞬態(tài)熱傳導(dǎo)問(wèn)題是可行的，具有節(jié)點(diǎn)布置靈活、前處理簡(jiǎn)單、計(jì)算簡(jiǎn)便等優(yōu)點(diǎn)．

1 Rashid Y R．Analysis of axisymmetric composite structures by the finite element method．Nuclear Engineering and Design，1996，3(1):163～182

2 Dargush G F，Banerjee P K．Application of the boundary element method to transient heat conduction．International Journal for Numerical Methods in Engineering，1991，31:1231～1244

3 Li Q H，Chen S S，Kou G X．Transient heat conduction analysis using the MLPG method and modified precise time step integration method．Journal of Computational Physics，2011，230:2736～2750

4 Li Q H，Chen S S，Zeng J H．A meshless model for transient heat conduction analyses of 3D axisymmetric functionally graded solids．Chinese Physics B，2013，22(12):120204

5 王峰，林皋，鄭保敬，劉俊，李建設(shè)．帶源參數(shù)熱傳導(dǎo)問(wèn)題的基于滑動(dòng)Kriging插值的MLPG法．力學(xué)季刊，2013，34(2):175～180(Wang F，Lin G，Zheng B J，Liu J，Li J S．Meshless local Petrov－Galerkin method with moving Kriging interpolation for solving heat conduction problems with source parameter．Chinese Quarterly of Mechanics，2013，34(2):175～180(in Chinese))

6 程榮軍，程玉民．帶源參數(shù)的二維熱傳導(dǎo)反問(wèn)題的無(wú)網(wǎng)格方法．力學(xué)學(xué)報(bào)，2007，39(6):843～847(Cheng R J，Cheng Y M．The meshless method for a two－dimensional inverse heat conduction problem with a source parameter．Acta Mechanica Sinica，2007，39(6):843 ～ 847(in Chinese))

7 程榮軍，程玉民．帶源參數(shù)的熱傳導(dǎo)反問(wèn)題的無(wú)網(wǎng)格方法．物理學(xué)報(bào)，2007，56(10):5569～5574(Cheng R J，Cheng Y M．The meshless method for solving the inverse heat conduction problem with a source parameter．Acta Physica Sinica，2007，56(10):5569 ～ 5574(in Chinese))

8 蔡永昌，朱合華，王建華．基于Voronoi結(jié)構(gòu)的無(wú)網(wǎng)格局部Petrov－Galerkin法．力學(xué)學(xué)報(bào)，2003，35(2):187～193(Cai Y C，Zhu H H，Wang J H．The meshless local Petrov－Galerkin method based on the Voronoi cells．Acta Mechanica Sinica，2003，35(2):187 ～193(in Chinese))

9 Cai Y C，Zhu H H．A meshless local natural neighbour interpolation method for stress analysis of solids．Engineering Analysis with Boundary Elements，2004，28:607～613

10 Sukumar N，Moran B，Belytschko T．The natural element method in solid mechanics．International Journal for Numerical Methods in Engineering，1998，43:839～887

11 Sukumar N，Moran B，Semenov Y．Natural neighbour Galerkin method．International Journal for Numerical Methods in Engineering，2001，50:1～27

12 Zhu H H，Liu W J，Cai Y C，Miao Y B．A meshless local natural neighbour interpolation method for two－dimensional incompressible large deformation analysis．Engineering Analysis with Boundary Elements，2007，31，856～862

13 李曉龍，王復(fù)明，徐平．自然元與無(wú)限元耦合方法在巖土工程粘彈性分析中的應(yīng)用．振動(dòng)與沖擊，2008，27(12):117～121，184(Li X L，Wang F M，Xu P．Coupling method of natural and infinite element for visco－ elastic analysis in geotechnical engineering．Journal of Vibration and Shock，2008，27(12):117 ～ 121，184(in Chinese))

14 王衛(wèi)東，張敦福，趙國(guó)群，程鋼．基于偶應(yīng)力理論的自然單元法研究．機(jī)械強(qiáng)度，2009，31(4):634～637(Wang W D，Zhang D F，Zhao G Q，Cheng G．Study on natural element method based on couple stresses theory．Journal of Mechanical Strength，2009，31(4):634～637(in Chinese))

*The project supported by the Foundation of Hunan Educational Committee(12C0059)

? Corresponding author E－mail:chenshenshen@tsinghua．org．cn

ANALYSIS OF TRANSIENT HEAT CONDUCTION PROBLEMS WITH A SOURCE PARAMETER BASED ON MESHLESS NATURAL NEIGHBOUR PETROV-GALERKIN METHOD*

Li Qinghua Chen Shenshen?
(School of Civil Engineering and Architecture，East China Jiaotong University，Nanchang330013，China)

Based on meshless natural neighbour petrov－Galerkin method，a novel meshless method was developed to solve transient heat conduction problems with a source parameter．The essential boundary conditions cannot be enforced directly when the non－interpolative moving least squares(MLS)approximation is used．In order to overcome this difficulty，the natural neighbour interpolation was employed instead of the moving least squares approximation to construct trial functions．The local weak forms of the transient heat conduction problems were satisfied locally in a series of polygonal sub－domains，which can be constructed easily with Delaunay tessellations．The traditional two－point difference technique was selected for the time discretization scheme．A numerical example demonstrates the validity and effectiveness of the presented method．

heat conduction problems， source parameter， meshless method， local Petrov－Galerkin method，natural neighbour interpolation

12 February 2014，

25 February 2014．

10．6052/1672－6553－2014－022

2014－02－12 收到第 1 稿，2014－02－25 收到修改稿．

*湖南省教育廳科研項(xiàng)目資助(12C0059)

E－mail:chenshenshen@tsinghua．org．cn