基于結(jié)合擴展精度技術(shù)的基本解方法的非線性功能梯度材料熱傳導(dǎo)問題求解

習強+傅卓佳+蔡加正

摘要： 采用基本解方法結(jié)合擴展精度技術(shù)和Kirchhoff變換求解功能梯度材料的二維熱傳導(dǎo)問題.在求解瞬態(tài)熱傳導(dǎo)問題時運用Laplace變換處理時間變量，將時域問題轉(zhuǎn)化為頻域問題求解；采用基本解方法計算得到高精度的頻域數(shù)值解，再分別采用Stehfest和Talbot這2種數(shù)值Laplace逆變換恢復(fù)原瞬態(tài)熱傳導(dǎo)問題的計算結(jié)果.通過3個非線性功能梯度材料的穩(wěn)態(tài)和瞬態(tài)熱傳導(dǎo)基準算例，分析結(jié)合擴展精度技術(shù)的基本解方法的計算精度與擴展精度位數(shù)、邊界布點數(shù)和虛擬邊界參數(shù)三者之間的關(guān)系.比較Stehfest和Talbot這2種數(shù)值Laplace逆變換算法的優(yōu)劣.采用結(jié)合擴展精度技術(shù)的基本解方法數(shù)值研究熱傳導(dǎo)系數(shù)隨位置劇烈變化的功能梯度材料熱傳導(dǎo)行為.數(shù)值結(jié)果表明該方法具有求解精度高、適用性好等特點，能高效模擬非線性功能梯度材料的二維穩(wěn)態(tài)與瞬態(tài)熱傳導(dǎo)行為.

關(guān)鍵詞： 基本解方法； 擴展精度技術(shù)； 數(shù)值Laplace逆變換； 功能梯度材料； 熱傳導(dǎo)

中圖分類號： O241.82；O343.6文獻標志碼： A

傅卓佳（1985—），男，浙江紹興人，副教授，博士，研究方向為徑向基函數(shù)、軟物質(zhì)力學建模、流固耦合、工程反問題、偏微分方程數(shù)值解、波傳播及結(jié)構(gòu)振動，（Email） paul212063@hhu.edu.cn

Abstract： The method of fundamental solution method in conjunction with the extended precision technique and Kirchhoff transformation is used to solve 2D heat conduction problem of nonlinear functionally graded materials. For transient heat conduction analysis， Laplace transform technique is applied to handle the time variable and the corresponding time domain problem can be transformed frequency domain problem. The numerical solution with highprecision in frequency domain is calculated by fundamental solution method. Two numerical Laplace inversion techniques Stehfest and Talbot are employed to retrieve the calculation results of original transient heat conduction. The effects of the extended precision digit， the boundary node number and fictitious boundary parameter on the calculation accuracy of the fundamental solution method in conjunction with the extended precision arithmetic are analyzed under three benchmark examples about steadystate and transient heat conduction problems in nonlinear functionally graded materials. The two numerical Laplace inversion techniques Stehfest and Talbot are compared. The thermal behavior of functionally graded materials with drasticchanged thermal conductivity is numerically studied by the fundamental solution method in conjunction with the extended precision arithmetic. The numerical results show that the method is with high accuracy and good applicability， and can effectively simulate 2D steadystate and transient heat conduction behavior of nonlinear functionally graded materials.

Key words： fundamental solution method； extended precision arithmetic； numerical Laplace inverse transform； functionally graded material； heat conduction

0引言

功能梯度材料是新一代復(fù)合材料，其導(dǎo)熱系數(shù)、比熱容和密度等物理參數(shù)沿著材料位置逐漸改變.[1]由于良好的熱傳導(dǎo)特性，功能梯度材料的使用環(huán)境多為高溫或超高溫[2-4]，比如：熱障涂層、返回艙的熱保護等，因此很有必要了解功能梯度材料的熱傳導(dǎo)性能，以便為實際工程中功能梯度材料的研發(fā)和設(shè)計提供參考.

在過去的十年里，國內(nèi)外學者開展大量關(guān)于功能梯度材料熱傳導(dǎo)行為數(shù)值模擬的研究工作.目前，求解此類問題的主要方法有有限元法和邊界元法.CHAROENSUK等[5]應(yīng)用有限元法研究功能梯度材料的瞬態(tài)熱傳導(dǎo)；SUTRADHAR等[6]運用邊界元法分析功能梯度材料的瞬態(tài)熱傳導(dǎo)問題.

近年來，無網(wǎng)格方法[7-11]發(fā)展迅速，由于其具有形式簡單、計算簡單、編程容易、收斂迅速，且可以很好地消除網(wǎng)格依賴缺陷，以及使求解精度變高等優(yōu)點，被廣泛應(yīng)用于求解熱傳導(dǎo)問題.基本解方法[12-14]是一種運用最為廣泛的無網(wǎng)格方法，該方法不需要對區(qū)域和邊界劃分網(wǎng)格，只需要選定一條物理求解域外的虛假邊界，并將源點布置在虛假邊界上以克服基本解的源點奇異性問題，從而可避免邊界元法及積分型邊界無網(wǎng)格方法中數(shù)學復(fù)雜且計算量大的奇異積分.

基本解方法在求解熱傳導(dǎo)問題時，理論上其精度應(yīng)是指數(shù)收斂的，但由于基本解方法通常離散得到稠密病態(tài)矩陣，當計算規(guī)模增大后，計算機本身的舍入誤差會造成數(shù)值結(jié)果的精度變差甚至出現(xiàn)錯誤結(jié)果.為避免稠密矩陣的病態(tài)性，奇異值分解技術(shù)被應(yīng)用于病態(tài)稠密矩陣的計算中，但是該技術(shù)只能在犧牲計算精度的前提下在一定程度上緩解這一問題.隨著航空航天、核能等工程領(lǐng)域的快速發(fā)展，對計算模型數(shù)值模擬功能梯度材料熱傳導(dǎo)行為的計算精度提出更高的要求，上述結(jié)合奇異值分解技術(shù)的計算模型已無法滿足這一需求.

關(guān)于瞬態(tài)熱傳導(dǎo)分析，時間方向的數(shù)值離散方法主要有時間步進方法和Laplace變換技術(shù).[15]時間步進方法雖然使用簡單，但是由于誤差累積和時間步長對計算效率和穩(wěn)定性的影響，使其在模擬長時間歷程熱傳導(dǎo)行為時計算效率較低，而Laplace變換技術(shù)則可以很好地避免時間步進方法中誤差累積和時間步長選取的問題.常用的數(shù)值Laplace變換技術(shù)有Stehfest算法[16]和Talbot算法[17]等，但研究發(fā)現(xiàn)數(shù)值Laplace變換技術(shù)通常是數(shù)學不適定性的，因此需要選取合適的頻域項數(shù)以確保Stehfest算法和Talbot算法等數(shù)值Laplace變換技術(shù)得到正確的計算結(jié)果.

隨著計算機技術(shù)的快速發(fā)展，擴展精度技術(shù)[17]得到越來越多的關(guān)注.TSAI等[18-19]將擴展精度技術(shù)運用于基本解方法，成功求解穩(wěn)態(tài)熱傳導(dǎo)問題和Helmholtz方程本征值問題，LING[20]將擴展精度技術(shù)運用于Kansa方法成功求解穩(wěn)態(tài)熱傳導(dǎo)問題.另一方面，數(shù)值Laplace變換技術(shù)結(jié)合擴展精度技術(shù)可以提高計算精度.[17]上述研究發(fā)現(xiàn)，擴展精度技術(shù)能在保證計算精度的同時避免稠密矩陣的病態(tài)性，并且可以很好地解決數(shù)值Laplace變換技術(shù)數(shù)學不適定性的問題.

此外，實際工程中的功能梯度材料熱傳導(dǎo)問題常常為一類非線性偏微分方程問題，常用的數(shù)值方法包括迭代解法[21]和Kirchhoff變換解法[22].不同于迭代解法，Kirchhoff變換解法可以將一些特定形式的非線性偏微分方程問題應(yīng)用Kirchhoff變換轉(zhuǎn)化為線性問題進行求解，可大大提高計算精度.

本文將擴展精度技術(shù)應(yīng)用于基于Laplace變換和Kirchhoff變換技術(shù)的基本解方法，用于求解二維非線性功能梯度材料熱傳導(dǎo)問題.

在求解瞬態(tài)熱傳導(dǎo)問題的過程中，首先對方程進行Laplace變換，采用結(jié)合擴展精度技術(shù)的基本解方法計算得到高精度的頻域數(shù)值解，隨后分別運用Stehfest算法和Talbot算法進行數(shù)值Laplace逆變換得到原瞬態(tài)熱傳導(dǎo)問題的計算結(jié)果.由于本文方法的計算誤差主要出現(xiàn)在基本解方法計算Laplace變換后的頻域問題和數(shù)值Laplace逆變換2個數(shù)值計算過程中，因此以基準算例2的精確解推導(dǎo)得到Laplace變換后的頻域精確解-T（x，p）為參考，確定上述2個數(shù)值計算過程分別對最終計算結(jié)果精度的影響.同時分別采用Stehfest算法和Talbot算法這2種數(shù)值Laplace逆變換算法進行計算，其中Stehfest算法采用30項頻域項數(shù)，即M=30，Talbot算法采用15項頻域項數(shù)，即M=15.

邊界布點數(shù)統(tǒng)一取N=60，虛擬邊界參數(shù)統(tǒng)一取d=8.在x1=0.5上均勻布置11個點作為測試點，分別取5個不同時刻，即t=0.002，0.010，0.020，0.050和0.100，計算結(jié)果的絕對誤差見圖5，圖中MFS表示基本解方法.

從圖5中可以看出：基本解方法結(jié)合擴展精度技術(shù)能提高計算精度，并且數(shù)值解與精確解經(jīng)過數(shù)值Laplace逆變換后的絕對誤差相比較，數(shù)值解精度基本沒有損失，這就說明誤差主要來源于數(shù)值Laplace逆變換的過程.

數(shù)值結(jié)果顯示，Talbot算法僅需采用15項頻域項數(shù)即可得到Stehfest算法采用30項頻域項數(shù)的計算精度.數(shù)值Laplace逆變換技術(shù)計算時間比較見表1.從表1可以看出：Talbot算法的計算時間少于Stehfest算法的計算時間.由此可認為，Talbot算法在求解這類瞬態(tài)熱傳導(dǎo)問題時比Stehfest算法更有效.

由圖7可知：隨著β1的增大，傳統(tǒng)基本解方法的誤差越來越大，當β1≥400時，傳統(tǒng)基本解方法在部分點處無法得到正確結(jié)果，而結(jié)合擴展精度技術(shù)的基本解方法依然適用.并且在β1相同的情況下，結(jié)合擴展精度技術(shù)的基本解方法的求解精度高于傳統(tǒng)基本解方法.

經(jīng)過大量的數(shù)值實驗發(fā)現(xiàn)，傳統(tǒng)基本解方法只能適用于β1≤800的情況（當β1=850時，邊界布點數(shù)N=320）；當β1>800時，傳統(tǒng)基本解方法在部分點處無法得到正確結(jié)果.結(jié)合擴展精度技術(shù)的基本解方法可以準確求解β1=2 000>800的情況，取邊界布點數(shù)N=160，β1=800，1 500和2 000，計算結(jié)果見圖8.

3結(jié)論

采用基本解方法結(jié)合擴展精度技術(shù)，Kirchhoff變換和坐標轉(zhuǎn)換求解非線性功能梯度材料的二維穩(wěn)態(tài)和瞬態(tài)熱傳導(dǎo)問題，運用Laplace變換處理時間變量，然后分別采用Stehfest和Talbot這2種數(shù)值Laplace逆變換用于恢復(fù)時間相關(guān)解.通過3個基準算例，得到以下結(jié)論.

（1）對比擴展精度位數(shù)、邊界布點數(shù)和虛擬邊界參數(shù)這3個因素對結(jié)合擴展精度技術(shù)的基本解方法的計算精度的影響，發(fā)現(xiàn)計算精度隨著擴展精度位數(shù)、邊界布點數(shù)以及虛擬邊界參數(shù)的增加而增加.此外，數(shù)值結(jié)果顯示當擴展精度位數(shù)G=0.6N時，計算效率最佳.

（2）比較Stehfest和Talbotp這2種數(shù)值Laplace逆變換算法發(fā)現(xiàn)，雖然Stehfest算法不需要進行復(fù)數(shù)運算，在取相同頻域項數(shù)M的前提下，計算時間更少，但Talbot算法的基本解方法只需取15項頻域項數(shù)便可以達到Stehfest算法取30項頻域項數(shù)所得到的計算精度，因此認為Talbot算法具有較高的計算效率.

（3）在數(shù)值研究熱傳導(dǎo)系數(shù)隨位置劇烈變化的功能梯度材料熱傳導(dǎo)行為時，發(fā)現(xiàn)傳統(tǒng)基本解方法只能適用于β1<800的情況，當β1再增大時，傳統(tǒng)基本解方法在部分點上已經(jīng)出現(xiàn)不正確的結(jié)果，而結(jié)合擴展精度技術(shù)的基本解方法可以準確求解β1=2 000>800的情況，本文方法具有求解精度高，適用性好等特點，能更高效地模擬非線性功能梯度材料的穩(wěn)態(tài)與瞬態(tài)熱傳導(dǎo)行為.

參考文獻：

[1]SU Z， JIN G， SHI S， et al. A unified solution for vibration analysis of functionally graded cylindrical， conical shells and annular plates with general boundary conditions[J]. International Journal of Mechanical Sciences， 2014（80）： 6280. DOI： 10.1016/j.ijmecsci.2014.01.002.

[2]劉俊俏， 苗福生， 李星. 二維各向異性功能梯度材料熱傳導(dǎo)的邊界元分析[J]. 西安交通大學學報， 2013， 47（5）： 7781. DOI： 10.7652/xjtuxb201305014.

LIU J Q， MIAO F S， LI X. Boundary element analysis of the two dimensional heat conduction equation for anisotropic functional graded materials[J]. Journal of Xian Jiaotong University， 2013， 47（5）： 7781. DOI： 10.7652/xjtuxb201305014.

[3]李慶華， 陳莘莘， 薛志清， 等. 基于自然元法和精細積分法的功能梯度材料瞬態(tài)熱傳導(dǎo)問題求解[J]. 計算機輔助工程， 2011， 20（3）： 2528. DOI： 10.3969/j.issn.10060871.2011.03.006.

LI Q H， CHEN S S， XUE Z Q， et al. Transient heat conduction solution to functionally graded materials based on natural element method and precise integration method[J]. Computer Aided Engineering， 2011， 20（3）： 2528. DOI： 10.3969/j.issn.10060871.2011.03.006.

[4]陳建橋， 丁亮. 功能梯度材料瞬態(tài)熱傳導(dǎo)問題的MLPG方法[J]. 華中科技大學學報（自然科學版）， 2007， 35（4）： 119121. DOI： 10.3321/j.issn：16714512.2007.04.036.

CHEN J Q， DING L. A MLPG method of transient heat transference in FGMs[J]. Journal of Huazhong University of Science and Technology（Nature Science）， 2007， 35（4）： 119121. DOI： 10.3321/j.issn：16714512.2007.04.036.

[5]CHAROENSUK J， VESSAKOSOL P. A high order control volume finite element procedure for transient heat conduction analysis of functionally graded materials[J]. Heat and Mass Transfer， 2010， 46（11）： 12611276. DOI： 10.1007/s0023101006498.

[6]SUTRADHAR A， PAULINO G H. The simple boundary element method for transient heat conduction in functionally graded materials[J]. Computer Methods in Applied Mechanics & Engineering， 2010， 193（42）： 45114539. DOI： 10.1016/j.cma.2004.02.018.

[7]程玉民. 科學和工程計算的新方法——無網(wǎng)格方法[J]. 計算機輔助工程， 2009， 18（1）： III. DOI： 10.3969/j.issn.10060871.2009.01.001.

CHENG Y M. New computational method in science and engineering： meshless method[J]. Computer Aided Engineering， 2009， 18（1）： III. DOI： 10.3969/j.issn.10060871.2009.01.001.

[8]MORADIDASTJERDI R， FOROUTAN M， POURASGHAR A. Dynamic analysis of functionally graded nanocomposite cylinders reinforced by carbon nanotube by a meshfree method[J]. Materials & Design， 2012， 44（3）： 256266. DOI： 10.1016/j.matdes.2012.07.069.

[9]ZHAO X， LIEW K M. A meshfree method for analysis of the thermal and mechanical buckling of functionally graded cylindrical shell panels[J]. Computational Mechanics， 2010， 45（4）： 297310. DOI： 10.1007/s0046600904468.

[10]DUAN Y. A note on the meshless method using radial basis functions[J]. Computers & Mathematics with Applications， 2008， 55（1）： 6675. DOI： 10.1016/j.camwa.2007.03.011.

[11]YANG D S， LING J. A new boundarytype meshfree method for solving the multidomain heat conduction problem[J]. Numerical Heat Transfer， Part B： Fundamentals： An International Journal of Computation and Methodology， 2016， 69（2）： 112. DOI： 10.1080/10407790.2015.1092823.

[12]BARNETT A H， BETCKE T. Stability and convergence of the method of fundamental solutions for Helmholtz problems on analytic domains[J]. Journal of Computational Physics， 2008， 227（14）： 70037026. DOI： 10.1016/J.JCP.2008.04.008.

[13]SUN L， CHEN W， CHENG H D. Method of fundamental solutions without fictitious boundary for plane time harmonic linear elastic and viscoelastic wave problems[J]. Computers & Structures， 2016（162）： 8090. DOI： 10.1016/j.compstruc.2015.08.018.

[14]YANG F L， YAN L， WEI T. Reconstruction of part of a boundary for the Laplace equation by using a regularized method of fundamental solutions[J]. Inverse Problems in Science & Engineering， 2009， 17（17）： 11131128.

[15]SUNDARARAJAN D. A practical approach to signals and systems： chapter 11： the Laplace transform[M]. Singapore： John Wiley & Sons（Asia） Pte Ltd， 2010： 301302. DOI： 10.1002/9780470823552.ch11.

[16]FU Z， CHEN W， QIN Q， et al. Three boundary meshless methods for heat conduction analysis in nonlinear FGMs with Kirchhoff and Laplace transformation[J]. Advances in Applied Mathematics & Mechanics， 2012， 4（5）： 519542. DOI： 10.4208/aamm.10m1170.

[17]ABATE J， VALKó P P. Multiprecision Laplace transform inversion[J]. International Journal for Numerical Methods in Engineering， 2004， 60（5）： 979993. DOI： 10.1002/nme.995.

[18]TSAI C C， LIN P H. On the exponential convergence of the method of fundamental solutions[J]. International Journal of Computational Methods， 2013， 10（2）： 292299. DOI： 10.1142/S0219876213410077.

[19]TSAI C C， YOUNG D L. Using the method of fundamental solutions for obtaining exponentially convergent Helmholtz eigensolutions[J]. Computer Modeling in Engineering and Sciences， 2013， 94（2）： 175205.

[20]LING L. Arbitrary precision computations of variations of Kansas method[J]. Progress on Meshless Methods， 2009（11）： 7783. DOI： 10.1007/9781402088216_5.

[21]李華剛， 石智偉. 求解非線性偏微分方程的傅里葉變換方法[J]. 廣東第二師范學院學報， 2011， 31（3）： 4447. DOI： 10.3969/j.issn.10078754.2011.03.008.

LI H G， SHI Z W. The method of solving nonlinear partial differential equation by applying fourier transform[J]. Journal of Guangdong Education Institute， 2011， 31（3）： 4447. DOI： 10.3969/j.issn.10078754.2011.03.008.

[22]BAGNALL K R， MUZYCHKA Y S， WANG E N. Application of the kirchhoff transform to thermal spreading problems with convection boundary conditions[J]. IEEE Transactions on Components Packaging & Manufacturing Technology， 2014， 4（3）： 408420. DOI： 10.1109/TCPMT.2013.2292584.（編輯武曉英）