基于結合擴展精度技術的基本解方法的非線性功能梯度材料熱傳導問題求解

習強+傅卓佳+蔡加正

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

關鍵詞： 基本解方法； 擴展精度技術； 數值Laplace逆變換； 功能梯度材料； 熱傳導

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

傅卓佳（1985—），男，浙江紹興人，副教授，博士，研究方向為徑向基函數、軟物質力學建模、流固耦合、工程反問題、偏微分方程數值解、波傳播及結構振動，（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引言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3結論

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

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

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

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

參考文獻：

[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]劉俊俏， 苗福生， 李星. 二維各向異性功能梯度材料熱傳導的邊界元分析[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]李慶華， 陳莘莘， 薛志清， 等. 基于自然元法和精細積分法的功能梯度材料瞬態熱傳導問題求解[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]陳建橋， 丁亮. 功能梯度材料瞬態熱傳導問題的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]程玉民. 科學和工程計算的新方法——無網格方法[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.（編輯武曉英）