边界节点法模拟瞬态热传导及Matlab工具箱开发

合集下载

热传导模型方程matlab

热传导模型方程引言热传导模型方程是描述热传导过程的数学模型，它在物理学、工程学等领域中有着广泛的应用。

该方程可以用来求解物体内部的温度分布以及热量的传输过程。

本文将重点介绍热传导模型方程在Matlab中的实现方法及应用。

热传导模型方程热传导模型方程是一种偏微分方程，用于描述物体内部温度分布随时间的变化。

它的一般形式如下：∂u/∂t = κ∇²u其中，u是温度分布函数，t是时间，κ是热传导率，∇²是拉普拉斯算子。

离散化处理为了在计算机中求解热传导模型方程，我们需要对其进行离散化处理。

通常，我们将物体分成若干个小区域，每个小区域的温度可以看作是一个离散节点。

然后，我们利用数值差分方法将偏微分方程转化为差分方程。

最常用的数值差分方法是有限差分法，其中最简单的方法是二阶中心差分法。

二阶中心差分法二阶中心差分法是一种常用的数值差分方法，它能够较为准确地近似热传导模型方程。

该方法利用了节点的前后两个时间步长的温度值和空间步长的温度梯度。

具体来说，对于一个二维网格，节点(i, j)的二阶中心差分近似可以表示为：(∂u/∂t)(i, j) ≈ (u(i+1, j) - 2u(i, j) + u(i-1, j))/(Δx²) + (u(i, j+1) - 2u(i, j) + u(i, j-1))/(Δy²)其中，Δx和Δy分别是空间步长。

Matlab实现在Matlab中，我们可以通过编写代码来求解热传导模型方程。

下面是实现该方程的一个示例代码：% 定义初始条件和边界条件N = 50; % 网格节点数L = 1; % 系统尺寸T = 1; % 总时间dx = L/N;dt = dx^2/4;x = linspace(0, L, N+1);y = linspace(0, L, N+1);[X, Y] = meshgrid(x, y);u = sin(pi*X).*sin(pi*Y); % 初始温度分布% 迭代求解热传导模型方程for t = 0:dt:Tu_new = u + dt*(del2(u)/dx^2);u = u_new;end% 可视化温度分布surf(X, Y, u)xlabel('x')ylabel('y')zlabel('Temperature')在上面的代码中，我们首先定义了初始条件和边界条件，然后利用循环迭代的方式求解热传导模型方程。

边界节点法计算二维瞬态热传导问题

Vol．35，No． 2，Feb．15，2014 2014 年2 月15 日出版文章编号: 1000-0887( 2014) 02-0111-10 应用数学和力学编委会，ISSN 1000-0887* 边界节点法计算二维瞬态热传导问题师晋红，傅卓佳，陈文( 1．河海大学力学与材料学院，南京210098;2．河海大学水文水资源与水利工程科学国家重点实验室，南京210098)( 我刊编委陈文来稿)摘要: 采用边界节点法( BKM) 结合双重互易法( DＲM) 求解二维瞬态热传导问题．采用差分格式处理时间变量，可将原瞬态热传导方程转化为一系列非齐次修正的Helmholtz 方程．随后，方程的解可分为特解和齐次解两部分计算，引入双重互易法在区域内部配点求解方程的特解，采用边界节点法仅需边界配点求解方程的齐次解．给出的数值算例显示该方法计算精度高，适用性好，具有很好的稳定性和收敛性，适合求解瞬态热传导问题．关键词: 瞬态热传导; 边界节点法; 双重互易法; 差分格式; 修正的Helmholtz 方程中图分类号: O241．82; O343．6文献标志码: Adoi: 10．3879 /j．issn． 1000-0887．2014．02． 001引言近年来，功能梯度材料、形状记忆合金等许多新兴材料被广泛用于高温涡轮发动机\［1 \］及伸缩翼系统\［2 \］的热障涂层等航天工程中．这些材料经常需要在高温环境中长期连续工作，所以很有必要了解这些材料及相关结构随时间演变的热传导性能，因此亟需发展计算瞬态热传导问题的数值算法，以便为航天工程材料的研发设计提供参考与指导．目前，很多数值算法\［3-9 \］已用于求解瞬态热传导问题，例如，有限差分法、有限元法、边界元方法等．有限差分法和有限元法在求解问题时要求在整个求解域上进行离散，在处理复杂几何区域等问题时存在计算量过大、网格生成困难等难题．边界元方法是一种边界型数值算法，将计算域降低了一维，然而边界元法选用含奇异性的基本解作为插值基函数，不可避免地需要处理费时费力的奇异积分计算问题．最近几十年发展起来的无网格方法\［10 \］，基于点的近似，克服网格依赖的缺陷，减少因网格畸变引起的计算困难，已广泛用于求解瞬态热传导问题\［5-7，11-12 \］．其中，基本解方法\［13 \］是* 收稿日期:基金项目:2013-10-07; 修订日期: 2013-11-19国家重点基础研究发展计划( 973 计划)( 2010CB832702 ) ; 国家杰出青年基金( 11125208) ; 国家自然科学基金( 11372097; 11302069 ) ; 高等学校学科创新引智计划( “111”计划) ( B12032)师晋红( 1989—) ，女，山西人，硕士生( E-mail: shijhhhu@ 163． com) ;作者简介:陈文( 1967—) ，男，江苏人，教授，博士，博士生导师( 通讯作者．edu． cn) ．111E-mail: chenwen@ hhu．112师晋红傅卓佳陈文与边界元法相对应的一种无网格方法．该方法无需积分，编程容易，不需要生成网格; 并只需在计算域边界上配点，且比边界元法的数值收敛速度快得多．然而为了避免基本解的源点奇异性，该法在物理边界外引入了虚假边界．虚假边界的设置有相当大的随意性，对于复杂几何区域和多连通区域，易造成计算不稳定．为了避免基本解方法布置虚假边界的问题，陈文等提出了边界节点法\［14-15 \］( boundaryknot method ，BKM) ．该方法使用控制微分方程的非奇异径向基函数通解代替奇异基本解，在克服基本解方法虚假边界问题的同时，保留了该法其他优点．边界节点法计算二维和三维几何复杂域各类物理问题的精度和稳定性都很高\［16 \］，有关论文总共已被他人引用 400 余次．目前该方法已被广泛应用于声学和扩散问题\［17-18 \］．\［19 \］结合边界型数值算法被在过去几十年中，双重互易法( dual reciprocity method ，D ＲM) 用来处理偏微分方程的非齐次项，并取得良好效果．尤其是双重互易边界元法( dual reciprocity boundary element method ，D ＲBEM) 在边界元领域就深受大家欢迎，被广泛用于求解瞬态热传导问题\［19-21 \］．双重互易法的基本想法是用一系列的近似特解来得到问题的一个特解，使用方便、高效、灵活，是求方程特解的一个有效方法．基于边界节点法和双重互易法近年来取得的众多成果，本文将采用边界节点法结合双重互易法求解二维瞬态热传导问题．文章首先采用差分格式处理热传导方程中的时间变量，将原热传导方程转化为一系列非齐次修正的 Helmholtz 方程．方程的解可分为特解和齐次解两部分．本文采用双重互易法求解方程的特解，采用边界节点法求解方程的齐次解．文章第 1 节描述热传导问题，并采用差分格式转化原热传导方程; 第 2 节介绍双重互易法和边界节点法; 第 3 节给出 3 个二维问题的数值算例并分析讨论; 第 4 节给出结论．1 问题描述设热传导问题的考察区域为 Ω，Γ Ω 为其边界．则问题的控制方程可表示如下:= ρc T ( x ，t ) ，k 2 T ( x ，t ) + Q ( x ，t ) = ( 1)x ∈ Ω．t 方程( 1) 受以下边界条件及初始条件的约束:－T ( x ，t ) q ( x ，t ) T ( x ，0)= T ( x ，t ) ， x ∈ ΓD ， ( 2) ( 3)( 4)为内热源．T ( x ，t ) = q －( x ，t ) ，，= x ∈ Γ n ( x )N = T 0 ，x ∈ Ω，其中 T ( x ，t ) 表示 t 时刻 x 点的温度，k 为导热系数，ρ 为密度，c 为比热容，Q ( x ，t ) －T ( x ，t ) ，q －( x ，t ) 分别表示边界上的温度和热流，T 表示初始温度，Γ 为 Dirichlet 边界，Γ 为0 D N Neumann 边界，Γ ΓD + ΓN ．= 采用差分格式处理原控制方程中的时间变量，对于时间间隔［t n ，t n +1］［0，t ］，物体中任一点的温度 T ( x ，t ) 及其对时间的一阶导数和内热源可以近似为\［22 \］T ( x ，t ) = T n +1( x ) ，( 5) ( 6) ( 7)n +1 n T ( x ，t ) t T ( x ) － T ( x ) ， = τ Q ( x ，t ) = Qn +1 ( x ) ，边界节点法计算二维瞬态热传导问题 113其中，τ t n+1 － t n为迭代过程中的时间步长．= 将方程( 5) 和( 6) 代入到方程( 1) 中，整理后可得(2ρc)T n+1( x ) = － ρc T n( x ) － Q n+1 ( x ) ． 1 ( 8) － τk k τk方程( 8) 中的右端项包含前一个时间步 t = t n 的近似解 T n ( x ) ，本文选取已知的初始条件T ( x ， 0) = T 0 作为迭代的第一步．方程( 8) 可简化为( 2 － λ2 ) T ( x )= f ( x ) ，( 9)其中ρc ， ( 10)=λ 槡f ( x ) = －ρcT n( x ) 1 －Q n+1( x ) ．( 11)τk方法介绍k2 方程( 9) 的解可写成如下特解与齐次解和的形式:T ( x ) = T p ( x ) + T h ( x ) ，式中，T p ( x ) 和 T h ( x ) 分别表示方程的特解和齐次解．特解 T p ( x ) 满足( 12)( 2 － λ2 ) T ( x ) = f ( x ) ，( 13) x ∈ Ω．p 齐次解 T h ( x ) 满足 ( 2 － λ2 ) T ( x ) = 0，x ∈ Ω，x ∈ ΓD ， x ∈ ΓN ．( 14) ( 15) ( 16)h －T h ( x ) q h ( x ) = T ( x ) － T p ( x )= q ( x ) － q p ( x )= g ( x ) ， = h ( x ) ，－双重互易法采用双重互易法求解方程的特解 T p ( x ) ．首先用如下线性组合来逼近右端项 f ( x )N If ( x ) =2． 1 \［23 \］:∑αi φi( x ) ， x ∈ Ω，( 17)i = 1 其中 αi 为待求系数，φi ( x ) = φ( ‖x － x i ‖) 确定 f ( x ) 后，特解 T p ( x ) 可写为N IT p ( x ) =其中 Ψi ( x ) 满足为径向基函数，x i 为区域 Ω 内N I 个不同的数据点． ∑αi Ψi ( x ) ，i = 1( 18)( 2－ λ2 ) Ψ ( x ) φ ( x ) ．( 19) = i i 如下\［24 \］:这里选取径向基函数 φ( x ) = r 2ln r ．( 20)φ 根据方程( 19) ，可得到 Ψ( x ) 的取值:2{4K0 ( λr ) = － 4 － 4 l n r r ln rλ2 r ≠ 0，Ψ －－ λ4 λ4 4 λ ( 21)4 + 4γ λ4 ( 2 )4 λ ln ， r = 0，= － +Ψ λ4 λ4114师晋红傅卓佳陈文其中 γ ≈ 0． 577 215 664 901 532 8 为 Euler 常数，K 0 为零阶第二类修正的 Bessel 函数．边界节点法采用边界节点法求解方程的齐次解 T h ( x ) ，T h ( x ) N ST h ( x ) =2． 2 \［14 \］由非奇异通解的线性组合来逼近 : ∑βi T i *( x ) ，x ∈ Ω，( 22)i = 1* * 2 其中 βi 为待求系数，T i( x ) ( ‖x － x i ‖) 为修正的 Helmholtz 算子 ( Δ － λ 的非奇异通) = T 解，x i 为源点集合，N S 为源点总数．对二维问题，非奇异通解 T *( x ) 为 T* ( λr ) ， ( 23)= I 0 其中I 0 为零阶第一类修正的 Bessel 函数．方程( 22) 中的 T h ( x ) 自动满足偏微分方程( 14) ，选取系数使得 T h ( x ) ( 15) 和( 16) ．本文采用配置点法来达到该目的，可得到如下线性方程组:N S也满足边界条件 ∑βi T i *( x j ) = g ( x j ) ， = 1，2，…，n 1 ， ( 24) j i = 1N S T *( x j ) i ∑βi = h ( x j ) ， = n 1 + 1，n 1 + 2，…，n b ，( 25)j n i = 1其中，n 1 为 Dirichlet 边界上的配置点个数，n b 为边界上的配点个数．写成矩阵形式则为 = b ， ( 26) TA β T其中 b = ( g ( x 1 ) ，g ( x 2 ) ，…，g ( x n ) ，h ( x n +1 ) ，h ( x n +2 ) ，…，h ( x n ) ) ，β ( β1 ，β2 ，…，βN ) ， = 1 1 1bSA = ( A ij ) 为插值矩阵:{T* ( x j ) ， = 1，2，…，N S ，j = 1，2，…，n 1 ，i i ( 27)A ij = T *( x j ) ， i = 1，2，…，N ，j = n + 1，n + 2，…，n i ．S 1 1 bn因此，方程( 9) 的解可表示为N IN ST ( x ) = T p ( x ) + T h ( x ) =∑αi Ψi ( x )i = 1*( x ) ， ( 28)+ ∑βi T i i = 1x ∈ Ω．数值算例本节采用边界节点法求解 3 个数值算例，从而检验该方法在求解瞬态热传导问题中的适用性及有效性．前两个算例所考察问题不含内热源，第 3 个算例包含内热源．3 为了量化数值解的好坏，分别引入绝对误差 E err ( T ) 和平均相对误差Ｒerr ( T ) 如下: E err ( T ) ，( 29)= T i － T i N∑( T i － T i )i = 12Ｒerr ( T ) ， ( 30)=N∑T2槡 i i = 1其中T i 和 T i 分别表示 x i 点的精确解和数值解，N 为测试点总个数．例 1 考虑一个热传导问题的经典算例，考察区域为正方形区域 Ω = { ( x 1 ，x 2 ) 0 ＜ x 1 ，x 2 ＜ 3 } ，各系数取值为 ρc = 1，k = 1． 25，Q = 0．初始时刻整个区域内的温度均边界节点法计算二维瞬态热传导问题 115为 30 ℃ ，边界温度均为 0 ℃ ，则该问题的边界条件及初始条件可以表述为 T ( 0，y ，t ) = T ( x ，0，t ) = T ( 3，y ，t ) = T ( x ，3，t ) = 0，T 0 ( x ，y ) 该问题的解析解为该问题的几何区域及边界条件见图 1．= 30． ∞∞–k π2 ( i 2 + j 2) t T ( x ，y ，t ) = ∑∑A i sin sin exp [i πx j πy ]， 3 3 23 i = 1 j = 1其中［( － 1) i－ 1］［( － 1) j－ A i = 4 × 30 × ． 2ij π目前已有很多学者采用有限元法( FEM) \［3 \］、边界元法 ( BEM) \［8 \］、双重互易边界元法\［19 \］、Trefftz 有限元法( Trefftz FEM) \［7 \］等方法求解该算例．这里将采用本文的方法 ( D ＲBEM) 求解该算例并与上述几种方法的数值结果进行比较．图1 例 1 的考察区域及边界条件 Fig ． 1 The domain and boundaryconditions in example 1图2 Fig ． 2 例 1 区域的边界及内部配点 The boundary and interpolationpoints in example 1在区域边界均匀配点 40 个，在区域内部均匀配点 81 个，见图 2．在计算中，时间步长取 τ = 0． 05，求 t = 1． 2 时刻的数值解．在区域中取 5 个测试点，采用边界节点法计算结果与其他方法的数值结果进行比较，见表 1．从表 1 中可以看出，本文方法所得数值解与解析解吻合得很好，并且比其他几种方法的数值精度更高．表1 = 1．2 时刻各方法数值结果比较 t Table 1 C omparison of numerical results bewteen different methods at t = 1． 2point x y BKM FEM BEM D ＲBE M Trefftz FEM exact 1 2 3 452． 4 2． 4 1． 8 1． 8 1． 51． 5 2． 4 1． 5 1． 8 1． 51． 081 1 0． 631 2 1． 778 5 1． 690 8 1． 871 01． 139 0． 670 1． 843 1． 753 1． 9381． 114 0． 657 1． 798 1． 713 1． 8871． 099 0． 645 1． 784 1． 695 1． 8771． 103 0． 660 1． 797 1． 715 1． 8941． 065 0． 626 1． 723 1． 639 1． 812为了检验该方法的收敛性，在边界取不同个数的配点分别进行计算．在整个区域均匀分布 441 个测试点，则 t 1． 2 时刻区域的平均相对误差见图 3．从图 3 中可以看出，数值解的精度= 随边界配点个数的增加而提高，最后趋于平稳，可见，该方法对边界配点个数是收敛的．例 2 为了检验边界节点法在求解不光滑区域问题上的适用性，考虑 L 形区域见图 4，各116师晋红傅卓佳陈文图3 t = 1．2 时不同边界配点个数的数值解精度 Fig ．3 Numerical accuracy variation with respect to the number of boundary points at t = 1． 2图4 例 2 中考察区域及 A ，B 测试点 Fig ． 4 The domain and 2 test pointsA ，B in example 2图5 Fig ． 5 例 2 区域的边界及内部配点 The boundary and interpolationpoints in example 2图6 例 2 中 A 点的数值解 Fig ． 6 Numerical solution at point A in example 2图 7 例2 中 B 点的数值解 Fig ． 7 Numerical solution at point B in example 2= 1，k = 1，Q = 0．该问题的精确解为 T ( x ，y ，t ) = cos( 3x ) sinh( 2y ) e －5t．系数取值为 ρc 该问题所考察区域的边界及内部配点见图 5，边界配点个数为 80，内部配点个数为 261．在计算中，时间步长取 τ = 0． 02．在区域中取 A ，B 两个测试点见图4，A ，B 两点随时间的温度变边界节点法计算二维瞬态热传导问题 117化分别见图 6、图 7．从图中可以看出，计算所得的数值解与精确解吻合得非常好，精度很高．可见，对于区域不光滑的问题，该方法同样适用．当取不同时间步长时，A ，B 两点随时间变化的数值解的绝对误差见图 8．从图中可以看出，时间步长取值越小，数值解的精度越高，同时所需的计算时间也越多．图8 例 2 中时间步长对数值解绝对误差的影响 Fig ． 8 Effect of time step length on the absolute numerical error in example 2例 3 考虑含内热源的瞬态热传导问题，考察区域见图 9．各系数取值为 ρc = 1，k = 1， Q ( x ，y ，t ) = sin x sin 2y ( 5cos t － sin t ) ．该问题的精确解为T ( x ，y ，= sin x sin 2y cos t ．在区域边界均匀配点 108 个，在区域内部均匀配点 415 个，见图 10．在计算中，时间步长取 τ = 0． 05．在整个区域均匀分布1 643 个测试点，t =2 时刻整个区域的数值解及绝对误差分布分别见图 11 和图 12．其中，最大绝对误差为 8． 196 × 10 － 5，其相应的数值解约为 0． 19，相对误差为 4． 262 × 10 － 4．因此，可以看出，该方法同样适用于求解含内热源的瞬态热传导问题，计算结果具有很高的精度．图9 例 3 的考察区域 Fig ． 9 The domain in example 3图10 例 3 区域的边界及内部配点 Fig ． 10 The boundary and interpolation points in example 3结论4 本文将边界节点法结合双重互易技术应用于求解瞬态热传导问题．其中，采用差分格式处理热传导方程中的时间变量，将原热传导方程转化为一系列修正的 Helmholtz 方程，并采用双重互易法求解方程的特解，采用边界节点法求解方程的齐次解．文中给出 3 个数值算例，从不118师晋红傅卓佳陈文图11 t = 2 时整个区域的数值解Fig ． 11 The numerical solution in the whole domain at t = 2图 12 t = 2 时整个区域数值解的绝对误差分布Fig ． 12 The absolute numerical error distribution in the whole domain at t = 2同边界形状、是否含内热源等方面检验了该方法的适用性及有效性，并与其他方法进行比较．数值结果表明，边界节点法无需积分，编程容易、不需要生成网格和划分虚假边界，计算精度高、适用性好，且具有很好的稳定性和收敛性，适合求解瞬态热传导问题．此方法还有望高效稳定求解瞬态热传导的三维问题及非线性问题，这将是我们以后的研究方向．参考文献( Ｒeferences) :［1］Burris K W ，Beardsley M B ，C huzho y L ． Component having a functionally graded material coating for improved performance \［P \］． U S Patent: 6087022，2000-7-11．Ｒ e nner E ． Thermal engine \［P \］． U S Patent: 3937019，1976-2-10．Bruch J C ，Z yvoloski G ． T ransient tw o -dimensional heat conduction problems solved by the finite ele-m ent m etho d \［J\］． I nt ernat ional Journal f or N umerical M et hods in Engineering ，1974，8( 3) : 481-494．B lobner J ，Bia ecki Ｒ A ，Kuhn G ． T ransient non-linear heat conduction-radiation pro blems —a bo und- ary elem ent fo rmulatio n \［J\］． I nt ernational Journal for N umeric al M et hods in Engineering ，1999，46( 11) : 1865-1882．［2］［3］［4］边界节点法计算二维瞬态热传导问题 119［5］Li Q H ，C h en S S ，Kou G X ． Tra ns ient hea t c onducti on a na l ys is us ing the M L PG method a nd modi-fied precise time step integration method \［J\］． Journal of Computational Physic s ，2011，230( 7) : 2736- 2750．Valtchev S S ，Ｒ oberty N C ． A time-marching M FS scheme for heat conduction problems\［J\］． Eng i- neering A nal y sis W i t h Boundary Elements ，2008，32( 6) : 480-493．Jirousek J ，Q in Q H ． Application of hybrid-T refftz element approach to transient heat conduction anal- ys is\［J\］． C omput ers ＆ St ruc t ures ，1996，58( 1) : 195-201．Brebbia C A ，T elles J C F ，W robel L C ． Boundary Element Techniques : T heory and A pplic ations in Eng i- neering \［M \］． Berlin: Springer-Verlag ，1984．欧阳华江．广义热传导方程有限元算法的计算准则 \［J\］．应用数学和力学，1992，13( 6) : 563- 571． ( O U YAN G H ua-jiang ． Criteria for finite element algorithm o f generalized heat conduction equa- tion \［J\］． A pplied M athematic s a nd M ec hanic s ，1992，13( 6) : 563-571． ( i n C hinese) )张雄，刘岩．无网格法\［M \］．北京: 清华大学出版社，2004． ( Z HAN G Xiong ，LIU Yan ． Mesh- less Methods \［M \］． Beijing: T singhua U n iversity Press ，2004． ( in Chinese) )Wang H ，Q in Q ，Kang Y ． A meshless model for transient heat conduction in functionally graded ma- teria ls\［J\］． C omputat ional M echanic s ，2006，38( 1) : 51-60．Fu Z J ，C hen W ，Q in Q H ． Three bounda ry meshless methods for heat conduction analysis in nonlin- ear FGM s w ith Kirchhoff and Laplace transformation \［J\］． A dvances in A pplied Mathematic s and Mechan- ic s ，2012，4( 5) : 519-542．Fairw eather G ，Karageorghis A ． The method of fundamental solutions fo r elliptic boundary value pro b- lems\［J\］． A dvances i n C omputational M at hematic s ，1998，9( 1 /2) : 69-95．Chen W ，T anaka M ． A meshless ，integratio n-free ，and boundary-only Ｒ BF technique \［J\］． Comput- er s ＆ M at hematic s W i t h A ppl i cations ，2002，43( 3) : 379-391．Fu Z J ，C hen W ，Q in Q H ． Boundary knot method for heat conduction in nonlinear functionally gra- ded ma teria l\［J\］． Engineering A nal y si s W i t h Boundary Elements ，2011，35( 5) : 729-734． Chen W ，H on Y C ． N umerical investigation on conve rgence of bounda ry knot method in the analysis of homogeneous Helmho ltz ，modified Helmholtz ，and convection-diffusion problems\［J \］． C omputer M et hods in A pplied M echanics and Engineering ，2003，192( 15) : 1859-1875．Hon Y C ，C hen W ． Boundary knot method fo r 2D and 3D Helmholtz and convection-diffusion prob-lems under complica ted geometry \［J\］． I nt ernat ional Journal for N umeric al M et hods in Engineering ，2003，56( 13) : 1931-1948． Jin B ，Z heng Y ． Boundary knot method for some inverse problems associated w ith the Helmholtz e-qua tion \［J\］． I nt ernat ional Journal f or N umerical M et hods in Engineering ，2005，62( 12) : 1636-1651．P artridge P W ，Brebbia C A ，Wro bel L C ． The Dual Ｒec iproc ity Boundary Element M ethod \［M \］．Southa mpton: Computa ti ona l M e c ha nics Publ ica ti ons ，1992．Bia ecki Ｒ A ，Jurga s＇ P ，Kuhn G ． D ual reciprocity BEM w ithout m atrix inversion fo r transient heat conduction \［J\］． Engineering A nal y si s W i t h Boundary Elements ，2002，26( 3) : 227-236．Bulgakov V ， arler B ，Kuhn G ． Iterative solution of systems of equations in the dual reciprocity bound- ary element method for the diffusion equation \［J \］． International J ournal for Numeric al Methods in Eng i- neering ，1998，43( 4) : 713-732． Chapko Ｒ，Kress Ｒ．Ｒ othe ＇s metho d for the heat equation and boundary integral equations\［J \］．Journal of I ntegral Equat ions a nd A ppli c ations ，1997，9( 1) : 47-69． Golberg M ．Ｒ e cent developments in the numerical evaluation of particular solutions in the boundary el-［6］［7］［8］［9］［10］［11］［12］［13］［14］［15］［16］［17］［18］［19］［20］［21］［22］［23］120师晋红傅卓佳陈文ement method \［J\］．A ppl i ed M at hematic s and C omputat ion，1996，75( 1) : 91-101．［24］ Chen C ，Ｒ ashed Y F． Evaluation o f thin plate spline based particular solutions fo r Helm holtz-type o p-era tors for the DＲ M \［J\］．M echanic s Ｒesearch C ommunicat ions，1998，25( 2) : 195-201．Boundary Knot Method for 2D Transient HeatConduction P roblem sSHI Jin-hong， FU Zhuo -jia， CHEN Wen( 1．C oll e ge of M echanic s and M at eri a l s，H o hai U ni v ersi t y，N anjing210098，P．Ｒ． C hina;2．Sta te Key Laboratory of H ydrolog y-W a ter Ｒesourc es and H ydraul i c Eng ineering，H ohaiU niversity，N a njing 210098，P．Ｒ． C hina)( Contributed by CHEN Wen，M ．AM M Editorial Boa rd)Abstract: The bounda ry knot method ( BKM ) in c onjunction w i th t he dual reciprocity method ( DＲM) w a s i ntroduce d to solve 2D tra nsi ent he at conduction proble ms．With the fini te diffe rence scheme appli ed to dea l w ith the time deri vative term ，the t ra nsie nt heat c onduction equation w as conver-ted to a se t of nonhomogeneous modifie d Helmholtz equati ons． Then the numerical sol ution to the non-homogeneous problems w a s divide d int o tw o pa rts: the partic ula r soluti on and the homoge neoussolution．The DＲM w i th few inner int erpol a tion nodes w a s empl oye d to get the pa rticul a r s ol uti on，a nd the BKM w i th bounda ry-onl y nodes us ed to obtain the homogeneous solution．N u m erical res ults show that the prese nt c ombi ned method has t he merits of high acc urac y，w ide a pplicabi lity，good sta bil it y and ra pi d conve rge nce，w h ic h w e re appe ali ng t o sol vi ng t ra nsie nt heat conduction proble ms．Key wor ds: transie nt heat c onduction; boundary knot method; dual reci procit y method; di ffere nce scheme; modi fied Helmholtz equati onFoundation item: The N ational Basic Ｒe sear ch Program of China ( 973 Program) ( 2010CB832702) ;The N ational Scie nce Fund for Distinguishe d Young Sc hola rs of China( 11125208) ;The N a tional N a tural Science Foundation of China( 11372097; 11302069)。

无单元 Galerkin 方法数值求解二维瞬态热传导问题

无单元 Galerkin 方法数值求解二维瞬态热传导问题摘要：采用无单元 Galerkin 方法数值求解具有狄利克雷边界条件的二维瞬态热传导问题。

首先离散该问题的时间变量，将该问题转化为与时间无关的边值问题；然后采用罚函数法处理 Dirichlet 边界条件，得到数值离散方程组，再利用 Matlab 软件求解给出的算例，结果表明该方法得到的数值结果与解析解吻合较好，该方法具有较高的计算精度和较好的收敛性。

关键词：二维瞬态热传导问题；无单元 Galerkin 方法；罚函数法；Matlab 软件中图分类号：G210.7 文献标识码：A 文章编号：20200168685有限差分法（FDM）[1]、有限元法（FEM）[2]、边界元法（BEM）得[3]、无网格法[4]等数值方法是解决瞬态热传导问题的常用方T k1T k1法。

Burlayenko 等人[5]用梯度有限元法计算了梯度材料中的瞬Tk1vd k vd k vdy2y态温度。

Sutradhar 和 Paulino[6]提出了一种简单的边界元法，该法只考虑边界条件，适用于梯度材料中的三维非定常热传导问T k1vd k1Q vdk1vd题。

Sladek 等人[7]用无网格局部边界积分方程方法考虑梯度材料中的不稳定温度场。

无单元 Galerkin 方法是应用最为广泛的无网格方法之一。

本文借鉴文献[8]应用无单元 Galerkin 方法数值求解具有狄利克雷边界条件的二维瞬态热传导问题。

首先离散该问题的时间变量，将该问题转化为与时间无关的边值问题；然后采用罚函数法（7）接下来，根据高斯公式化简（7）式，于是得与混合边值问k1x H1 ，使得，有题（6）等价的变分问题：求T处理 Dirichlet 边界条件，得到数值离散方程组，再利用 Matlab 软件求解给出的算例，最后给出了一个数值算例来验证理论误差a T k1,v（8）T k1vdvQ k1dR k1vd vT k2d结果。

利用Matlab进行数值模拟的方法

利用Matlab进行数值模拟的方法引言数值模拟是现代科学领域中不可或缺的一种工具，它通过数学模型和计算机算法，模拟和预测实际系统的行为。

随着科学技术的不断发展，数值模拟方法逐渐成为各个学科的重要组成部分。

Matlab作为一种强大的科学计算工具，为数值模拟提供了丰富的函数库和易于使用的编程环境。

本文将介绍一些利用Matlab进行数值模拟的方法，以及其在不同领域的应用。

一、常微分方程的数值解法常微分方程在物理、工程、生物等领域中广泛存在。

利用Matlab进行常微分方程的数值解法，可以有效地求得方程的近似解。

Matlab中的ode45函数是常用的数值解法之一，它基于龙格-库塔算法，可以处理非刚性和刚性问题。

通过设定初始条件和方程形式，利用ode45函数可以得到系统的数值解，并绘制出相应的曲线图。

例如，考虑一个一阶常微分方程dy/dx = -2xy，初始条件为y(0) = 1。

可以通过以下代码进行数值模拟：```Matlabfun = @(x, y) -2*x*y;[x, y] = ode45(fun, [0, 10], 1);plot(x, y)xlabel('x')ylabel('y')title('Solution of dy/dx = -2xy')```运行以上代码后，可以得到方程解的图像，从而对其行为有更直观的理解。

二、偏微分方程的数值解法偏微分方程在物理、流体力学、电磁学等领域中具有重要应用。

常用的偏微分方程的数值解法有有限差分法(Finite Difference Method)和有限元法(Finite Element Method)等。

在Matlab中，可以利用pdepe函数进行偏微分方程的数值模拟，其中包含了一维和二维问题的求解算法。

以热传导方程为例，假设一个长为L的均匀杆子，其温度分布满足偏微分方程∂u/∂t = α*∂²u/∂x²，其中u(x, t)表示温度分布。

MATLAB PDE-tool在热传导问题中的应用

MATLAB PDE-tool在热传导问题中的应用
熊静;张薇
【期刊名称】《工业加热》
【年(卷),期】2009(038)004
【摘要】介绍了MATLAB软件PDE工具箱求解导热问题解题步骤,并且列举了两个热传导方面的例子来说明其应用.模拟结果在工程质量控制以及在线测试材料导热系数方面具有重要的工程使用价值.
【总页数】3页(P42-44)
【作者】熊静;张薇
【作者单位】南京工业大学材料科学与工程学院,江苏南京210009;南京工业大学材料科学与工程学院,江苏南京210009
【正文语种】中文
【中图分类】TK124
【相关文献】
1.SODF-SPH方法及其在热传导问题中的应用 [J], 热合买提江·依明;阿合买提江·伊明江;买买提明·艾尼
2.lq稀疏正则化理论及其在热传导反问题中的应用 [J], 赵光伟
3.多域边界面法在稳态热传导问题中的应用 [J], 张见明;李湘贺;陆陈俊;李光耀
4.广义边界控制法在多层热传导边界识别问题中的应用 [J], 岳俊宏;李明;牛瑞萍
5.紧算子方程的不适定性分析及其在一维热传导反问题中的应用 [J], 闵涛;葛宁国;张敏;张帆
因版权原因，仅展示原文概要，查看原文内容请购买。

matlab传热计算程序

matlab传热计算程序
传热计算在工程学和科学领域中是一个重要的应用。

Matlab是一个功能强大的工程计算软件，可以用于传热计算。

在Matlab中，你可以使用各种方法来进行传热计算，比如有限元法、差分法、有限体积法等。

以下是一些常见的传热计算程序的示例：
1. 热传导方程求解，你可以编写一个Matlab程序来求解热传导方程，根据给定的边界条件和初始条件，使用差分法或有限元法来离散方程，并进行时间步进求解，得到温度场的分布。

2. 对流换热计算，对于流体内部的对流换热问题，你可以编写一个Matlab程序来求解Navier-Stokes方程和能量方程，结合有限体积法来进行流场和温度场的耦合求解。

3. 辐射换热计算，针对辐射换热问题，你可以编写一个Matlab程序来计算辐射传热，比如使用辐射传热方程和辐射传热模型，结合离散方法进行求解。

4. 传热系统优化，除了单一的传热计算，你还可以使用Matlab进行传热系统的优化设计，比如通过建立传热模型和耦合其
他工程模型，使用优化算法来寻找最优的传热系统设计参数。

总之，Matlab提供了丰富的工具和函数，可以用于传热计算的各个方面。

通过编写程序，你可以灵活地进行传热计算，并且可以根据具体的问题需求进行定制化的计算和分析。

希望这些信息对你有所帮助。

利用matlab程序解决热传导问题-推荐下载

1、题目及要求
1. 原始题目及要求 2. 各节点的离散化的代数方程 3. 源程序 4. 不同初值时的收敛快慢 5. 上下边界的热流量（λ=1W/(m℃）） 6. 计算结果的等温线图 7. 计算小结题目：已知条件如下图所示：
二、各节点的离散化的代数方程
各温度节点的代数方程
ta=(300+b+e)/4 ; tb=(200+a+c+f)/4; tc=(200+b+d+g)/4; td=(2*c+200+h)/4 te=(100+a+f+i)/4; tf=(b+e+g+j)/4; tg=(c+f+h+k)/4 ; th=(2*g+d+l)/4 ti=(100+e+m+j)/4; tj=(f+i+k+n)/4; tk=(g+j+l+o)/4; tl=(2*k+h+q)/4
0,0,-2,4,0,0,0,-1,0,0,0,0,0,0,0,0; -1,0,0,0,4,-1,0,0,-1,0,0,0,0,0,0,0; 0,-1,0,0,-1,4,-1,0,0,-1,0,0,0,0,0,0; 0,0,-1,0,0,-1,4,-1,0,0,-1,0,0,0,0,0; 0,0,0,-1,0,0,-2,4,0,0,0,-1,0,0,0,0; 0,0,0,0,-1,0,-1,0,4,0,0,0,-1,0,0,0; 0,0,0,0,0,-1,0,0,-1,4,-1,0,0,-1,0,0; 0,0,0,0,0,0,-1,0,0,-1,4,-1,0,0,-1,0; 0,0,0,0,0,0,0,-1,0,0,-2,4,0,0,0,-1; 0,0,0,0,0,0,0,0,-2,0,0,0,24,-1,0,0; 0,0,0,0,0,0,0,0,0,-2,0,0,-1,24,-1,0; 0,0,0,0,0,0,0,0,0,0,-2,0,0,-1,24,-1; 0,0,0,0,0,0,0,0,0,0,0,-1,0,0,-1,12]; b=[300,200,200,200,100,0,0,0,100,0,0,0,300,200,200,100]'; [x,n]=gauseidel(A,b,[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]',1.0e-6) xx=1:1:4; yy=xx; [X,Y]=meshgrid(xx,yy); Z=reshape(x,4,4); Z=Z' contour(X,Y,Z,30) Z= 139.6088 150.3312 153.0517 153.5639

利用matlab程序解决热传导问题

哈佛大学能源与环境学院课程作业报告作业名称：传热学大作业——利用matlab程序解决热传导问题院系：能源与环境学院专业：建筑环境与设备工程学号：姓名：盖茨比2015年6月8日一、题目及要求1.原始题目及要求2.各节点的离散化的代数方程3.源程序4.不同初值时的收敛快慢5.上下边界的热流量（λ=1W/(m℃））6.计算结果的等温线图7.计算小结题目：已知条件如下图所示：二、各节点的离散化的代数方程各温度节点的代数方程ta=(300+b+e)/4 ; tb=(200+a+c+f)/4; tc=(200+b+d+g)/4; td=(2*c+200+h)/4 te=(100+a+f+i)/4; tf=(b+e+g+j)/4; tg=(c+f+h+k)/4 ; th=(2*g+d+l)/4ti=(100+e+m+j)/4; tj=(f+i+k+n)/4; tk=(g+j+l+o)/4; tl=(2*k+h+q)/4tm=(2*i+300+n)/24; tn=(2*j+m+p+200)/24; to=(2*k+p+n+200)/24; tp=(l+o+100)/12 三、源程序【G-S迭代程序】【方法一】函数文件为：function [y,n]=gauseidel(A,b,x0,eps)D=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);G=(D-L)\U;f=(D-L)\b;y=G*x0+f;n=1;while norm(y-x0)>=epsx0=y;y=G*x0+f;n=n+1;end命令文件为：A=[4,-1,0,0,-1,0,0,0,0,0,0,0,0,0,0,0;-1,4,-1,0,0,-1,0,0,0,0,0,0,0,0,0,0;0,-1,4,-1,0,0,-1,0,0,0,0,0,0,0,0,0;0,0,-2,4,0,0,0,-1,0,0,0,0,0,0,0,0;-1,0,0,0,4,-1,0,0,-1,0,0,0,0,0,0,0;0,-1,0,0,-1,4,-1,0,0,-1,0,0,0,0,0,0;0,0,-1,0,0,-1,4,-1,0,0,-1,0,0,0,0,0;0,0,0,-1,0,0,-2,4,0,0,0,-1,0,0,0,0;0,0,0,0,-1,0,-1,0,4,0,0,0,-1,0,0,0;0,0,0,0,0,-1,0,0,-1,4,-1,0,0,-1,0,0;0,0,0,0,0,0,-1,0,0,-1,4,-1,0,0,-1,0;0,0,0,0,0,0,0,-1,0,0,-2,4,0,0,0,-1;0,0,0,0,0,0,0,0,-2,0,0,0,24,-1,0,0;0,0,0,0,0,0,0,0,0,-2,0,0,-1,24,-1,0;0,0,0,0,0,0,0,0,0,0,-2,0,0,-1,24,-1;0,0,0,0,0,0,0,0,0,0,0,-1,0,0,-1,12];b=[300,200,200,200,100,0,0,0,100,0,0,0,300,200,200,100]';[x,n]=gauseidel(A,b,[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]',1.0e-6) xx=1:1:4;yy=xx;[X,Y]=meshgrid(xx,yy);Z=reshape(x,4,4);Z=Z'contour(X,Y,Z,30)Z =139.6088 150.3312 153.0517 153.5639108.1040 108.6641 108.3119 108.1523 84.1429 67.9096 63.3793 62.4214 20.1557 15.4521 14.8744 14.7746 【方法2】>> t=zeros(5,5);t(1,1)=100;t(1,2)=100;t(1,3)=100;t(1,4)=100;t(1,5)=100;t(2,1)=200;t(3,1)=200;t(4,1)=200;t(5,1)=200;for i=1:10t(2,2)=(300+t(3,2)+t(2,3))/4 ;t(3,2)=(200+t(2,2)+t(4,2)+t(3,3))/4;t(4,2)=(200+t(3,2)+t(5,2)+t(4,3))/4;t(5,2)=(2*t(4,2)+200+t(5,3))/4;t(2,3)=(100+t(2,2)+t(3,3)+t(2,4))/4;t(3,3)=(t(3,2)+t(2,3)+t(4,3)+t(3,4))/4; t(4,3)=(t(4,2)+t(3,3)+t(5,3)+t(4,4))/4; t(5,3)=(2*t(4,3)+t(5,2)+t(5,4))/4;t(2,4)=(100+t(2,3)+t(2,5)+t(3,4))/4;t(3,4)=(t(3,3)+t(2,4)+t(4,4)+t(3,5))/4;t(4,4)=(t(4,3)+t(4,5)+t(3,4)+t(5,4))/4;t(5,4)=(2*t(4,4)+t(5,3)+t(5,5))/4;t(2,5)=(2*t(2,4)+300+t(3,5))/24;t(3,5)=(2*t(3,4)+t(2,5)+t(4,5)+200)/24;t(4,5)=(2*t(4,4)+t(3,5)+t(5,5)+200)/24;t(5,5)=(t(5,4)+t(4,5)+100)/12;t'endcontour(t',50);ans =100.0000 200.0000 200.0000 200.0000 200.0000 100.0000 136.8905 146.9674 149.8587 150.7444 100.0000 102.3012 103.2880 103.8632 104.3496 100.0000 70.6264 61.9465 59.8018 59.6008 100.0000 19.0033 14.8903 14.5393 14.5117【Jacobi迭代程序】函数文件为：function [y,n]=jacobi(A,b,x0,eps)D=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);B=D\(L+U);f=D\b;y=B*x0+f;n=1;while norm(y-x0)>=epsx0=y;y=B*x0+f;n=n+1;end命令文件为：A=[4,-1,0,0,-1,0,0,0,0,0,0,0,0,0,0,0;-1,4,-1,0,0,-1,0,0,0,0,0,0,0,0,0,0; 0,-1,4,-1,0,0,-1,0,0,0,0,0,0,0,0,0; 0,0,-2,4,0,0,0,-1,0,0,0,0,0,0,0,0;-1,0,0,0,4,-1,0,0,-1,0,0,0,0,0,0,0; 0,-1,0,0,-1,4,-1,0,0,-1,0,0,0,0,0,0; 0,0,-1,0,0,-1,4,-1,0,0,-1,0,0,0,0,0;0,0,0,-1,0,0,-2,4,0,0,0,-1,0,0,0,0;0,0,0,0,-1,0,-1,0,4,0,0,0,-1,0,0,0;0,0,0,0,0,-1,0,0,-1,4,-1,0,0,-1,0,0;0,0,0,0,0,0,-1,0,0,-1,4,-1,0,0,-1,0;0,0,0,0,0,0,0,-1,0,0,-2,4,0,0,0,-1;0,0,0,0,0,0,0,0,-2,0,0,0,24,-1,0,0;0,0,0,0,0,0,0,0,0,-2,0,0,-1,24,-1,0;0,0,0,0,0,0,0,0,0,0,-2,0,0,-1,24,-1;0,0,0,0,0,0,0,0,0,0,0,-1,0,0,-1,12];b=[300,200,200,200,100,0,0,0,100,0,0,0,300,200,200,100]'; [x,n]=jacobi(A,b,[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]',1.0e-6); xx=1:1:4;yy=xx;[X,Y]=meshgrid(xx,yy);Z=reshape(x,4,4);Z=Z'contour(X,Y,Z,30)n =97Z =139.6088 150.3312 153.0517 153.5639108.1040 108.6641 108.3119 108.152384.1429 67.9096 63.3793 62.421420.1557 15.4521 14.8744 14.7746四、不同初值时的收敛快慢1、[方法1]在Gauss 迭代和Jacobi 迭代中，本程序应用的收敛条件均为norm(y-x0)>=eps ，即使前后所求误差达到e 的-6次方时，跳出循环得出结果。

matlab偏微分方程工具箱使用手册

MATLAB偏微分方程工具箱使用手册一、Matlab偏微分方程工具箱介绍Matlab偏微分方程工具箱是Matlab中用于求解偏微分方程（PDE）问题的工具。

它提供了一系列函数和工具，可以用于建立、求解和分析PDE问题。

PDE是许多科学和工程领域中的重要数学模型，包括热传导、扩散、波动等现象的数值模拟、分析和优化。

Matlab偏微分方程工具箱为用户提供了丰富的功能和灵活的接口，使得PDE问题的求解变得更加简单和高效。

二、使用手册1. 安装和启用在开始使用Matlab偏微分方程工具箱前，首先需要确保Matlab已经安装并且包含了PDE工具箱。

确认工具箱已经安装后，可以通过以下命令启用PDE工具箱：```pdetool```这将打开PDE工具箱的图形用户界面，用户可以通过该界面进行PDE 问题的建立、求解和分析。

2. PDE建模在PDE工具箱中，用户可以通过几何建模工具进行PDE问题的建立。

用户可以定义几何形状、边界条件、初值条件等，并选择适当的PDE方程进行描述。

PDE工具箱提供了各种几何建模和PDE方程描述的选项，用户可以根据实际问题进行相应的设置和定义。

3. 求解和分析一旦PDE问题建立完成，用户可以通过PDE工具箱提供的求解器进行求解。

PDE工具箱提供了各种数值求解方法，包括有限元法、有限差分法等。

用户可以选择适当的求解方法，并进行求解。

求解完成后，PDE工具箱还提供了丰富的分析功能，用户可以对结果进行后处理、可视化和分析。

4. 结果导出和应用用户可以将求解结果导出到Matlab环境中，并进行后续的数据处理、可视化和分析。

用户也可以将结果导出到其他软件环境中进行更进一步的处理和应用。

三、个人观点和理解Matlab偏微分方程工具箱是一个非常强大的工具，它为科学和工程领域中的PDE问题提供了简单、高效的解决方案。

通过使用PDE工具箱，用户可以快速建立、求解和分析复杂的PDE问题，从而加快科学研究和工程设计的进程。

导热的反问题matlab

导热的反问题matlab
在MATLAB中，导热问题通常涉及热传导方程的数值求解。

热传导方程描述了物体内部温度分布随时间的变化，通常采用偏微分方程来描述。

解决导热问题的一种常见方法是使用有限差分法。

在MATLAB中，可以通过编写代码来离散化热传导方程，并使用迭代方法求解离散化后的方程。

另一种常见的方法是使用MATLAB的偏微分方程工具箱（Partial Differential Equation Toolbox）。

该工具箱提供了一系列函数和工具，可以帮助用户建立和求解偏微分方程，包括热传导方程。

用户可以通过定义边界条件、初始条件和热传导方程的参数来建立模型，并使用工具箱中的函数进行数值求解。

此外，MATLAB还提供了用于可视化和分析结果的丰富工具，例如绘制温度分布图、计算热通量等。

通过这些工具，用户可以全面分析导热问题的结果，并对模型进行验证和优化。

总之，在MATLAB中，可以通过编写代码、使用偏微分方程工具箱以及可视化分析工具来解决导热问题，从而全面深入地研究热传导现象。

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

边界节点法模拟瞬态热传导及Matlab工具箱开发万长风;林继;洪永兴【摘要】针对高温条件下材料热传导问题,采用边界节点法数值模拟二维和三维瞬态热传导问题并编写了相应的计算软件。

边界节点法是无网格方法的一种,采用径向基函数的非奇异通解来满足控制方程。

对于热传导方程,采用Implicit-Euler差分离散其时间项,将不含时间项的非齐次亥姆赫兹方程以特解与通解的和的形式给出。

使用双重互易法在计算区域配点求出方程的特解,使用边界节点法在边界配点来求出方程的齐次解。

使用结果表明：该方法计算精度高,数学简单,收敛快,能很好地数值模拟二维及三维瞬态热传导问题。

基于上述的算法,使用Matlab语言开发了一个热传导问题的数值模拟工具箱,该工具箱具有简单易操作、代码开源、界面友好等特点。

【期刊名称】《能源与环保》【年(卷),期】2017(039)002【总页数】6页(P12-17)【关键词】瞬态热传导边界节点法双重互易法差分格式 Matlab工具箱【作者】万长风;林继;洪永兴【作者单位】河海大学力学与材料学院,江苏南京211100【正文语种】中文【中图分类】O241.82某些工程材料在一些热处理过程中需要在高温条件下长期工作，因此研究高温环境中材料的传热有着重要的意义[1]。

由于实验过程耗时费力而且有时对设备要求较高，因此研究一种能够简单效率地计算热传导问题，能为热处理工艺和新兴材料的研发与应用提供一定的参考价值。

国内外很多学者已经使用一些常见算法例如有限元方法来模拟热传导问题[2-7]。

有些网格类算法需要进行复杂的网格划分，在解决一些复杂问题时常会遇到网格划分困难，耗时费力的网格重构等问题，这些问题直接影响着计算的效率和计算精度。

为了克服网格类算法的这些不易解决的缺点，无网格方法受到学术界的广泛重视，文献[8]中给出了近年来无网格粒子类方法的研究进展。

相比于网格类算法，此类算法使用边界或区域内配点的信息来求得插值矩阵，不需要在区域内构建网格，因此能够快速收敛且后处理简单。

无网格算法已经被广泛地应用于各种科学和工程问题的研究[5-6，9-10]。

例如基本解方法是由Kurprade和Aleksidze在1964年提出[11]，基本解法只需得到配置点的信息并求得插值矩阵，因此比网格类算法的收敛速度快。

但是在源点上基本解会出现奇异性问题，为了避免该问题，该方法在物理边界外设置了虚假边界，而虚假边界的设置一般凭借经验，导致数值结果变得不稳定[12]。

为了避免虚假边界的困扰，陈文教授提出了仅在物理边界上离散而不设置虚假边界的边界节点法(Boundary knot method，BKM)[13-14]，该方法使用径向基函数非奇异通解来满足控制方程，所以该方法不需要在计算区域边界外设置额外的虚假边界且保持了计算的稳定性与收敛性。

本文采用边界节点法并结合双重互易法(Dual reciprocity method，DRM)[15-17]模拟瞬态热传导问题。

首先，基于Implicit-Euler差分格式离散热传导方程中的时间项，那么原热传导问题的解可写成特解和通解2部分。

齐次解可由满足控制方程的非奇异一般解的线性组合来近似，特解可以用双重互易法选取适当的径向基函数并在计算区域配点求解方程来得到，源项由径向基函数近似从而可以有效的给出热传导问题的数值解。

最后基于上述算法，本文介绍了以上述算法为原理编写的Matlab工具箱，此工具箱可解决瞬态热传导问题，为数值方法的比较研究、新的算法讨论提供参考数值结果；也可用于工程问题的计算，为各种工程问题提供较为准确的参考。

假定瞬态热传导方程的计算区域为Ω，边界为Γ=∂Ω，热传导方程可表示为式(1)：边界条件约束如下：初始条件如下：T(x，0)=T0，x∈Ω其中，T(x，t)为在t时刻x点的温度；k为材料导热系数；ρ为材料的密度；c为比热容；Q(x，t)为材料内部热源；为边界上的温度；为边界上的法向热流量；T0为任意区域内一点的初始时刻的温度；Γ=ΓD+ΓN。

使用Implicit-Euler差分格式来离散方程中的时间项，在任意一段时间间隔[tn，tn+1]⊂[0，t]内，物体任一点的温度T(x，t)由T(x，tn+1)来近似，任一点的温度对时间的一阶导数由整个时间段内温度的平均变化率来近似，任一点的内热源Q(x，t)由Q(x，tn+1)来近似[18]：T(x，t)=θT(x，tn+1)+(1-θ)T(x，tn)Q(x，t)=θQ(x，tn+1)+(1-θ)Q(x，tn)其中，n为步数；dt=tn+1-tn为步长；θ为可指定的参数，这里取θ=1，即Implicit-Euler方式。

将方程代入原始方程中：方程左边的温度t=tn+1的求解需要前一个时间步t=tn的近似解Tn(x)，所以使用方程的初始条件进行迭代。

方程也可简化为：(▽2-λ2)T(x)=f(x)，x∈Ω其中方程的解可写成以下特解与齐次解和的形式：T(x)=Tp(x)+Th(x)其中，Tp(x)为方程的特解；Th(x)为齐次解。

特解Tp(x)满足式(13)：(▽2-λ2)Tp(x)=f(x)，x∈Ω需要注意的是，Tp(x)只需要满足方程而不需要满足边界条件，而齐次解Th(x)满足边界条件而不需要满足方程：(▽2-λ2)Th(x)=0，x∈Ω2.1 双重互易法为了求解得到特解Tp(x)，需引入双重互易法用如下的线性组合来近似方程的右端项f(x)：其中，αi为待求的插值矩阵；NI为配置点个数；φi(x)为选取的径向基函数。

通过在区域内配点可以求解出未知系数αi，从而方程的特解Tp(x)可以表示为式(18)：其中，ψi(x)满足：(▽2-λ2)ψi(x)=φi(x)本文选取如下的径向基函数φ(x)：φ=r2lnr根据方程，相应的ψ(x)的取值如下：其中，γ≈0.577 215 664 901 532 8，为欧拉常数；k0为零阶第二类修正的贝塞尔函数。

2.2 边界节点法本节求解齐次解Th(x)，采用边界节点法引入控制方程的非奇异一般解的线性组合来近似[14]式中，{βi}为待求系数；为修正的Helmholtz算子(Δ-λ2)的非奇异通解；{xi}为源点集合；Ns为源点总数。

对于二维和三维问题，修正Helmholtz方程的非奇异通解T*(x)分别可表示为：T*=I0(λr)，D=2其中，I0为零阶第一类修正的Bessel函数；D为方程所描述问题的维数。

方程中的Th(x)满足方程，所以只要选取适当的系数使得Th(x)使其满足边界条件和即可。

使用强格式配点法求解未知系数，从而可用以下的线性组合来表示：其中，nb1为在ΓD即狄利克雷边界上的配点个数；Ns为整个计算区域边界上的配点总个数，改写成矩阵形式为：Aβ=b，其中，b=(g(x1)，g(x2)，…，g(xnb1)，h(xnb1+1)，h(xnb1+2)，…，h(xNs))T，β=(β1，β2，…，βNs)T，A=(Aij)为插值矩阵：最终方程的数值解可表示为式(29)：基于前面介绍的算法，采用Matlab语言编写了一个易于使用的软件包，该软件使用Matlab环境编程，以其脚本语言为核心进行开发，利用基于径向基函数的边界节点法结合基于奇异值分解的正则化方法求解控制方程为Helmholtz方程的热传导问题。

该软件可在奇异计算网(http：//www.singular.software)上下载并使用。

由于计算过程中的插值矩阵可能为高度病态矩阵，因此引用了Per Christian Hansen开发的基于奇异值分解的矩阵正则化工具箱[19-20]，即软件包子目录下名称为“regu”的文件夹。

操作界面主要分为3部分：Question，Condition和Result。

Question部分主要是对于热传导问题的描述，Condition部分是对于边界条件的描述，而Result部分的主要功能是查看程序的输出结果。

各部分含有输入窗口、图形显示窗口或者各功能按钮。

表1给出了各按钮以及对应的功能。

软件的主要操作流程包括方程的前处理，计算和结果输出3个过程。

前处理主要是对相应参数进行选择和设定，包括边界形状，求解问题的控制方程形式，正则化方法等。

用户输入所需数据后，单击Compute按钮，程序将进行计算并输出结果。

(1)单击reset按钮清除所有的输入，清除之后输入问题的控制方程以及各个参数，并选择问题的维数。

当在计算区域上内热源全部为0或者在边界形状不确定的情况下，不需要输入Q，length 和radius等参数。

(2)输入问题的边界条件与计算传热时间T。

同时可以使用Import boundary nodes 按钮来导入边界点的text数据。

该按钮也可以在边界可选择的情况下添加更多的边界点。

输入问题初始条件。

如果需要更多内部点，可以使用Import inner point 按钮。

(3)由于计算得到的插值矩阵可能是病态矩阵，所以本软件采用基于奇异值分解的正则化方法进行处理。

使用Regularization method按钮选择可行的正则化方法，最后点击compute 按钮便可计算得到相应结果。

本节采用边界节点法数值模拟热传导问题算例来说明该方法的精确度与收敛性。

前2个算例为二维瞬态热传导问题，第3个算例为多连通区域的三维瞬态热传导问题。

使用绝对误差Err(T)和平均相对误差Rerr(T)来说明数值结果的好坏：例1：使用边界节点法数值模拟一个Q(x，y，t)=0较简单的热传导问题。

问题计算区域为圆形各系数ρc=1，k=1，时间步长取dt=0.05。

该问题精确解为T(x，y，t)=cos(x)sin(y)exp(-2t)+cos(x)sin(2y)exp(-5t)。

使用软件计算得到的T=(0，1]时的相对误差曲线如图1所示，从图1中可以看出该算法求解精度较高，数值结果的平均相对误差在10-2，同时文中所提工具箱也能较好模拟算例1。

例2：考虑一个含内热源的热传导问题说明该算法的收敛性和有效性，考察区域为正方体Ω={(x1，x2，x3)|0<x1，x2，x3<3)，ρc=1，k=1，该问题的内热源可表示为Q(x1，x2，x3，t)=sin(x1)sin(x2)sin(x3)(3cos(t)-sin(t))。

该问题的精确解可表示为T(x1，x2，x3，t)=sin(x1)sin(x2)sin(x3)cos(t)。

时间步长dt分别为0.01，0.02和0.05时的相对误差曲线如图2所示，dt=0.05，t=0.8时的相对误差和布点数的关系如图3所示。

可以看出，数值解的相对误差随着的降低而减小。

在一定范围内，在边界上布置的点数越多，计算精度越高。

例3：考虑三维问题。