功能梯度材料圆柱壳与圆筒问题的二维分析
习强+傅卓佳+蔡加正
摘要： 采用基本解方法结合扩展精度技术和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.（编辑武晓英）
