2. 北京航空航天大学 航空科学与工程学院, 北京 100191;
3. 国防科学技术大学 航天科学与工程学院, 长沙 410073
2. School of Aeronautic Science and Engineering, Beijing University of Aeronautics and Astronautics, Beijing 100191, China;
3. College of Aerospace Science and Engineering, National University of Defense Technology, Changsha 410073, China
微分求积法是在近年里发展起来,并正在进一步发展的数值计算方法.它是一种高效率、高精度的微分方程的求解方法,计算工作量少,计算结果可靠,已用于运输[1]、结构力学[2, 3, 4, 5, 6, 7, 8]、流体力学[9, 10]和化工[11, 12]等多个领域.精细积分作为时域的求解方法,具有解析性好、精度高、计算稳定等特点.随着航空航天技术发展,高超声速技术已成为21世纪航空航天领域关注的热点问题,而气动加热是高超声速技术发展的重要障碍,也是飞行器关键技术问题之一,飞行器表面温度是非常复杂的,如何能简单、准确、高效地确定飞行器表面温度传播特性被人们急切关注.Tseng等[13]曾提出针对未知量为导热系数、边界温度和边界热流两两组合时的反演算法,但只有少量的数值试验.薛齐文等[14]提出了应用共轭梯度技术求解热传导边界条件反问题的一种数值模式.还有通过试验的方法根据布置在容器外壁或内部数个热电偶的温度测量值推算出内壁处温度和热流密度,如在稳态导热情况下采用外插法,但是对于瞬态导热问题该方法将会引起较大的误差[15].张驰等[16]提出了无网格边界元法求解有内热源的瞬态传热问题,研究了一种无需域内离散划分网格的纯边界元法,但该方法仅适用于二维问题,在三维域的应用有待进一步探讨研究.
本文给出一种求解瞬态传热问题的高效高精度方法,即采用微分求积法离散空间域,精细积分方法离散时间域,具有解析性好、精度高、计算稳定和能高效处理各种复杂边界条件的特性.
1 微分求积法微分求积法(Differential Quadrature Method,DQM)是一种全域逼近的数值方法,它的基本思想是以全域内各节点函数值的加权和来逼近函数偏导数在某一节点处的值.它与有限元法(Finite Element Method,FEM)有显著的区别:FEM通常采用低阶多项式逼近单元内的函数值,而DQM则是在全域内采用高阶多项式逼近域内某一连续函数;FEM是在单元内逼近函数值,由近似函数求得导数,而DQM综合全域内各节点的信息来直接逼近函数偏导数在某一点的值.因此,DQM用较少网格点就可以获得高的数值精度.
根据微分求积方法,函数f(x)在xi点的k阶导数可近似表示为
用微分求积法求解问题,必须选取合理分布的节点.均匀分布是最早被采用也是最简单的节点分布形式.但计算表明,非均匀节点具有更快的收敛速度和更高的求解精度.边界条件处理也是微分求积法中关键问题之一.对于二阶微分方程的求解每个端点只需一个条件,引入边界条件时直接将边界节点坐标代入边界条件即可.对于二阶以上高阶微分方程,每个端点的边界条件不止一个,边界节点却只有一个,因此边界条件的处理不能简单进行.解决这一问题的方法有多种,常见的有方程替代法、变量缩聚法、权系数矩阵修正法和边界自由度增添法.用这些方法可以处理复杂的边界条件.
2 精细积分法精细积分法是计算指数矩阵的一种高精度方法,其要点是利用指数函数加法定理计算指数矩阵的增量.考虑指数矩阵:
考虑隔热结构的热传导问题,如图 1所示.瞬态温度场U是时间和空间的函数,非稳态导热微分方程为
初始条件:
 |
| 图 1 隔热材料结构图Fig. 1 Thermal protection material structure |
边界条件:
在求解具体问题时,要用边界条件方程替换方程式(13)中相应的方程.本文选取使切比雪夫多项式等于零的点为节点.离散化的边界条件为
用向量和矩阵可以把域内平衡方程和边界条件方程写成
若边界条件都是时间无关向量,则式(17)为
从式(16)消去$\overline {
{U}} $ B,再用式(17)替换边界条件F,整理得
精细积分方法直接求解的方程是一阶微分方程,方程式(19)仅包含时间的一阶微分,因此不需要降阶,即不需要增加自由度,这也是本文选择精细积分方法的根据.
对方程式(19),根据李级数方法得到其解析解为
计算式(21)的关键问题就是其中指数矩阵e G h 的计算,其计算方法可参见第2节内容.进而由递推公式(21),可求得域内各节点温度随时间的变化值.
4 数值算例下面用本文方法分析某隔热材料的瞬态温度场.结构尺寸为0.2×0.2×0.025m3,热扩散系数2.47933884×10-7m2/s.
4.1 算 例 1采用三维模型,上表面1200℃恒温,下表面和四周面绝热状态,初始温度为25℃.分别采用微分求积(数值解)和NASTRAN计算它的温度响应.微分求积法中布置节点4×4×10.NASTRAN中划分20×20×10个八节点体单元.图 2给出了两种方法求解的结构下表面一点的温度响应图.图 3给出了两种方法求解的结构中间面一点的温度响应图.图 4给出了NASTRAN计算的982s状态结构的温度响应图.
由图 2~图 4可知:微分求积法计算得到的温度响应曲线和有限元计算结果吻合得很好,可见微分求积和精细积分结合的数值方法能高效高精度地求解瞬态传热问题,布置少量的节点就达到了比较高的精度.
 |
| 图 2 结构下表面一点的温度响应Fig. 2 Temperature response of one point at bottom of structure |
 |
| 图 3 结构中间面一点的温度响应Fig. 3 Temperature response of one point in the middle of structure |
 |
| 图 4 NASTRAN计算的第982s结构温度响应Fig. 4 Temperature response of structure at 982 s by NASTRAN |
考虑如下几种边界条件组合,其中环境温度为25℃.
第1类边界条件(Case1):上表面为1200℃恒温,其他面绝热.此时在式(14)中,C1jk=0,C2jk=0,C3ik=0,C4ik=0,C5ij=0,C6ij=1200.
第2类边界条件(Case2):上表面为1200℃恒温,其他面热交换,对流热交换系数为5000W/m2.C1jk=5000,C2jk=5000,C3ik=5000,C4ik=5000,C5ij=5000,C6ij=1200.
第3类边界条件(Case3):上表面温度变化规律为
第4类边界条件(Case4):上表面温度变化规律为
由此可以看出,针对不同的边界条件只需要改变相应的参数即可,求解格式是统一的.微分求积法一方面把二阶或高阶的边界条件并入加权系数中,另一方面把边界条件离散到各个空间点上,大大地简化了边界条件的处理.
图 5给出了3种不同初始温度下,图 1模型中上下表面中点处温度曲线图.由图 5可知前100s 时最下表面温度几乎没有变化,等于初始温度,主要是因为热量是从上表面传递过来的,它的传递是依赖于时间的,在最开始时热量只能传递上层位置;随着时间增长,热量传递加快,下表面温度也增加得越来越快,当达到一定时刻时,上下表面温差小到一定值时,热量传递趋于稳定,曲线又趋于平缓.
 |
| 图 5 3种不同初始温度下的温度响应Fig. 5 Temperature response with three different initial temperatures |
图 6给出绝热和热交换两种边界条件下下表面中点温度响应曲线.由图 6可知无论在何种边界下,开始 时刻温度变化缓慢,随着时间增长,温度迅速增加直至最后趋于平稳状态.但对于热交换边界条件,最开始时刻,下表面温度反而会下降,比初始时刻温度还低,之后逐渐增加.由图可知在任意时刻,绝热状态下温度总是比热交换条件高.因而在分析温度响应时正确地确定边界条件是至关重要的.
 |
| 图 6 不同边界条件下结构下表面一点温度响应Fig. 6 Temperature response of one point at bottom of structure with different boundary conditions |
本文综合利用精细积分法和微分求积法求解了瞬态温度场问题,该做法的优点是充分利用了两种方法的高效高精度特点.
1) 精细积分法对步长的依赖性很小,在有效位数范围内其结果和精确解一致，并且没有因为对方程降阶而增加自由度数的问题.
2) 微分求积法不依赖泛函和变分原理,具有数学原理简单、计算精度高、计算量和内存需求少等优点,能以较少的网格点和较小的计算量求得高精度的数值解.
3) 在处理复杂边界条件时,比如C1类边界条件,微分求积法具有普适性和灵活性.微分求积法具有精度高、效率高等特点,这是高阶方法的优点.
| [1] | Civan F, Sliepcevich C M.Application of differential quadrature to transport processes[J].Journal of Mathematical Analysis and Applications,1983,93(1):206-221. |
| Click to display the text | |
| [2] | Bert C W, Jang S K,Striz A G.Two new approximate methods for analyzing free vibration of structural components[J].AIAA Journal,1988,26(5):612-618. |
| Click to display the text | |
| [3] | Bert C W, Malik M.The differential quadrature method for irregular domains and application to plate vibration[J].International Journal of Mechanical Sciences,1996,38(6):589-606. |
| Click to display the text | |
| [4] | Jang S K, Bert C W,Striz A G.Application of differential quadrature to static analysis of structural components[J]. International Journal for Numerical Methods in Engineering,1989,28(3):561-577. |
| Click to display the text | |
| [5] | Bert C W, Wang X,Striz A G.Differential quadrature for static and free vibration analysis of anisotropic plates[J].International Journal of Solids and Structures,1993,30(13):1737-1744. |
| Click to display the text | |
| [6] | Liu F L,Liew K M. Static analysis of Reissner-Mindlin plates by differential quadrature element method[J].ASME Journal of Applied Mechanics,1998,65(3):705-710. |
| Click to display the text | |
| [7] | Wang X, Bert C W.A new approach in applying differential quadrature and free vibration analysis of beams and plates[J].Journal of Sound and Vibration,1993,162(3):566-572. |
| Click to display the text | |
| [8] | Xing Y F, Liu B.A differential quadrature analysis of dynamic and quasi-static magneto-thermo-elastic stresses in a conducting rectangular plate subjected to an arbitrary variation of magnetic field[J].International Journal of Engineering Science,2010,48(12):1944-1960. |
| Click to display the text | |
| [9] | Shu C, Richards B E.Application of generalized differential quadrature to solve two-dimensional incompressible Navier-Stokes equations[J].International Journal for Numerical Methods in Fluids,1992,15(17):791-798. |
| [10] | Malik M, Bert C W.Differential quadrature solution for steady state incompressible and compressible lubrication problems[J].ASME Journal of Teratology,1994,116(2):296-302. |
| [11] | Quan J R, Chang C T.New insights in solving distributed system equations by the quadrature methods-I.Analysis[J].Computers in Chemical Engineering,1989,13(7):779-788. |
| Click to display the text | |
| [12] | Quan J R, Chang C T.New insights in solving distributed system equations by the quadrature methods-II.Numerical experiments[J].Computers in Chemical Engineering,1989,13(9):1017-1024. |
| Click to display the text | |
| [13] | Tseng A A, Chen T C,Zhao F Z.Direct sensitivity coefficient method for solving two-dimensional inverse heat conduction problems by finite-element scheme[J].Numerical Heat Transfer,Part B:Fundamentals,1995,27(3):291-307. |
| Click to display the text | |
| [14] | 薛齐文,杨海天, 胡国俊.共轭梯度法求解瞬态传热组合边界条件多宗量反问题[J].应用基础与工程科学学报,2004,12(2):113-120. Xue Q W,Yang H T,Hu G J.Solving inverse heat conduction problems with multi-variables of boundary conditions in transient-state via conjugate gradient method[J].Journal of Basic Science and Engineering,2004,12(2):113-120(in Chinese). |
| Cited By in Cnki (9) | |
| [15] | France D M, Chiang T.Analytic solution to inverse heat conduction problem with periodicity[J].Journal of Heat Transfer,1980,102(3):579-581. |
| Click to display the text | |
| [16] | 张驰,石宏,张硕,等. 基于无网格边界元法的瞬态热传导问题研究[J].科学技术与工程,2013,13(26):7638-7643. Zhang C,Shi H,Zhang S,et al.Study on transient heat conduction by meshless boundary element method[J].Science Technology and Engineering,2013,13(26):7638-7643(in Chinese). |
| Cited By in Cnki () | |
| [17] | Xing Y F, Liu B.High-accuracy differential quadrature finite element method and its application to free vibrations of thin plate with curvilinear domain[J].International Journal for Numerical Methods in Engineering,2009,80(13):1718-1742. |
| Click to display the text | |
| [18] | 钟万勰. 结构动力方程的精细时程积分法[J].大连理工大学学报,1994,34(2):131-136. Zhong W X.On precise time-integration method for structural dynamics[J].Journal of Dalian University of Technology,1994,34(2):131-136(in Chinese). |
| Cited By in Cnki (366) | |