二维热传导方程数值解及MATLAB实现【总】

业 设 计 (论文 ) 专 业： 飞行器动力工程 学 号： 学生姓名： 所属学院 指导教师： 二 一一 年 六 月 中国民航大学本科生毕业设计 (论文 ) 二维热传导方程初边值问题数值求解的 有限元方法研究 值实现 A - 业 ： 飞行器动力工程 学生姓名 ： 学 号 ： 学 院 ： 指导教师 ： 2011 年 6 月 创见性声明 本人声明：所呈交的毕业论文是本人在指导教师的指导下进行的工作和取得的成果，论文中所引用的他人已经发表或撰写过的研究成果，均加以特别标注并在此表示致谢。与我一同工作的同志对本论文所做的任何贡献也已在论文中作了明确的说明并 表示谢意。 毕业论文作者签名： 签字日期： 年 月 日 本科毕业设计（论文）版权使用授权书 本毕业设计（论文）作者完全了解中国民航大学有关保留、使用毕业设计（论文）的规定。特授权中国民航大学可以将毕业设计（论文）的全部或部分内容编入有关数据库进行检索，并采用影印、缩印或扫描等复制手段保存、汇编以供查阅和借阅。同意学校向国家有关部门或机构送交毕业设计（论文）的复印件和磁盘。 （保密的毕业论文在解密后适用本授权说明） 毕业论文作者签名： 指导教师签名： 签字日期： 年 月 日 签字日期： 年 月 要 本文 采用有限元法结合不同的初边值条件解决二维热传导问题。 本课题要求用有限元方法对 三类）边界条件的 二维热传导方程初边值问题的数值求解 、 编程并进行 值实现 。 本文首先通过对时间、空间等定义域的离散将其转换成相对应的矩阵形式，然后再采用二维交叉隐格式求解方法将所要解决方程离散成对应矩阵，最后通过 解对应矩阵的线性方程组，解出特定时刻的温度随空间变化分布的矩阵并将其可视化。通过与 具箱 模拟数据 的对比，验证 本方法的可靠性 和稳定性，进而为二维热传导方程的数值解的实现及其可视化提供了另一种极为有效的方法。 关键 词 ： 二维 初边值 ； 热传导方程 ； 有限元法 ； 现 EM of to in At we in a in to a we a At we by to a of of By we of 录 第一章 绪论 . 1 课题背景和意义 . 1 课题研究现状 . 1 课题要求 . 2 题解决 . 2 论文结构安排 . 3 第二章 有限元法及偏微分方程解法理论分析 . 4 有限元基本知识介绍 . 4 求解二维热传导方程的基本思想 . 5 将定解区域离散 . 6 插值函数的选择 . 6 方程组的建立 . 6 方程组的求解 . 6 二维热传导方程 . 7 网格剖分 . 7 离散方程组的构建 . 7 稳定性分析 . 9 第三章 解二维热传导方程在 的实现 . 11 关知识简介 . 11 追赶法简介 . 11 相关函数命令简介 . 12 解二维热传导方程在 的实现方法 . 13 第四章 数值化和可视化举例分析 . 15 二维热传导方程数值求解的 现 . 15 究一种较为简单的情况： . 15 究一种较为复杂的情况 . 17 误差分析 . 19 结论与展望 . 21 致谢 . 23 参考文献 . 24 附录 A: 程序清单 . 25 附录 B: 外文翻译资料 . 25 国民航大学本科生毕业设计 （论文） 一章 绪论 课题背景和意义 热传导是 一种 普遍存在的物理现象 , 它广泛存在于目前的工程应用领域 中 ，对人们的生产和生活实践有着广泛而深刻的影响。 虽然 人们对于热传导现象的本质特征已有了一个比较完整的认识，并已通过 严密的数学逻辑公式 推导对其 物理 特征进行了精确的描述，但所得到的结果往往是复杂的 积分或级数表达式，其中还免不了使用 某些 特殊函数， 由于 现在数学上存在的困难，在工程中许多热传导问题还不能采用分析解法进行求解。因此 , 通过一种较为简便直观的方式掌握控制和改进热量传递的方法和技术措施 , 无论对国民经济建设还是改善人民生活都具有重要的意义 。 近年来，随着计算机技术迅速发展，数值方法已经得到广泛应用并成为有力的辅助求解工具， 人们在此基础上 已发展了诸如有限差分法、有限元法和边界元法等用于工程问题的求解 方法 。 我们知道，热传导问题的数学表达式一般为偏微分方程形式，但偏微分方程一般并没有固定有效的精确解 法，因而人们广泛采用数值方法来实现热传导方程的研究、模拟和仿真。对于热传导问题的数值计算及其可视化，原则上，可以用 语言来完成这个任务，可是实际上并没有人去这样做，原因是成本太高，需要消耗更多的人力和物力，达到的效果却并不能满足人们的需求。幸运的是，高性能数学软件（如 目前具有强大计算功能的个人计算机让这个问题变得简单起来。 本课题旨在尝试利用目前在理论上已经成熟的有限差分法原理结合 大的数据处理和模拟仿真功能编写特定的程序代码对给定初边值条件的二维热传导方程温度随 时间、空间的变化规律进行研究并将其可视化。 本课题将以二维热传导方程的数值解法及 现为主线，研究论证其可行性，从而发现一种较为简便且极为有效的热传导方程数值解法和可视化的方法，意在更好的解决目前在工程和研究领域中实际存在的热传导问题，进而推动其相关领域的发展和进步。 课题研究现状 近三十年来，解偏微分方程的理论和方法有了很大的发展，而且在各个学科技术的领域中应用也愈来愈广泛，在我国，偏微分方程数值解法作为一门课程，不但中国民航大学本科生毕业设计 （论文） 计算数学专业，而且也在其他理工科专业的研究生的大学生中开设 。同时，求解热传导方程的数值算法也取得巨大进展，特别是有限差分法方面，此算法的特点是在内边界处设计不同于整体的格式，将全局的隐式计算化为局部的分段隐式计算。而且精度上更好。 目前，在欧美各国 使用十分普及。在大学的数学、工程和科学系科，用作许多课程的辅助教学手段， 成为大学生们必不可少的计算工具，甚至 是一项必须掌握的基本技能。 在我国， 各大专院校的应用日益普遍，许多专业已把 为基本计算工具。在科研机构和工业界， 得到越来越广泛的应用 。 有强大的图形绘制功能，为科学计算和图形处理提供了很大的方便。我们只需制定的绘图方式，再提供绘图数据，有程序指令就可以得到形象、直观的图形结果。因此，近些年越来越多的人开始使用 求解数值计算和图形处理技术，我们也可以绘制出热传导方程数值解的二维、三维图形，从而可以更好的理解热传导方程的意义。 课题要求 用有限元方法对如下二维热传导方程初边值问题 进行 数值求解，编程并进行值实现： ( , , ) , ( , , ) ( 0 , ( , , 0 ) ( , ) , ( , )( , , ) 0 , ( , , ) ( 0 , uu f x y t x y t x y x y x yu x y t x y t 222 ， ( , , )f x y t 为热源函数； 为空间定义域； 为 的边界 。 要求求出给定时间情况下温度随空间的分布变化规律 。 题解决 目前，对于求解偏微分方程有很多方法，但差分法和有限元离散法式主要解决问题的两种方法。一般来说，用差分法来 解 偏微分方程，解得 的 结果就是方程的准确解函数 在 这 点上的近似值。而用变分近似的方法求解，是将近似解表示成有限维子空间中基函数的线性组合。有限元法也是基于变分原理，由于选择了特殊的基函数，使它能适用于一般的区域。这种基函数是与区域的剖分有关的，近似解 u 表示为基函数的线性组合，二线性组合中的系数，又是剖分节点上 u 或其导数的近似值。 有关一维热传导方程的有限差分法求解的 现 4，现在已经解决，本文中国民航大学本科生毕业设计 （论文） 鉴一维热传导有限差分法的 解思想，对二 维热传导方程进行转换，再对解法编程实现，从而进一步对热传导方程进行探讨。二维热传导方程求解在现实生活中的应用也更加广泛，所以有很好的现实意义。 论文 结构 安排 第一章， 绪论 。 主要介绍课题的研究背景和意义，分析了热传导现象在目前工程领域的主要解决方案和发展现状，最后说明论文的主要内容和组织结构。 第二章，有限元法理论研究。 主要介绍与二维热传导问题相关的 论，为有限元法在 的实现奠定基础。 第三章， 的实现。 主要介绍如何针对特定的一类二维热传导方程结合初边值条件编制 相应的程序代码。 第四章，二维热传导方程的数值化和可视化。 主要介绍如何利用编制好的程序代码带入相应的初边值条件，进行数值化和可视化研究。 中国民航大学本科生毕业设计 （论文） 二章 有限元法及偏微分方程解法理论 分析 本部分主要对二维热传导方程的有限元解法进行了完整的理论分析。首先从有限元的基本知识开始引出求解二维热传导方程的基本思想，然后根据 本思想使用交叉方向隐格式求解方法构建矩阵，并论证其稳定性。 有限元基本知识介绍 定义 含有未知函数12( , , , , )nu x x x 的方程称为偏微分方程。 定义 方程 111( ) ( ) ( , ) , u uk k f x tt x x x x 称热传导方程。其中， ( , )u u x t 是固体的传热过程中在 x 处、 t 时刻的温度。系数12 ( 0 ) nk k 程为 ( , ) , u u f x 其中 2 2 22 2 212 nx x x ， n 为维数 。 定义 在 特定条件下求解方程的解。这样的条件成为定解条件。给出了方程和定结条件，就构成了定解问题 。 定义 一般说，边界条件有下列形式 ( , ) ( , ) ( , ) ( , ) ( , ) ,ux y u x y x y x y x 其中 为边界的外法向导数。有如下几种特殊形式 （ 1） 第一类）条件： 0, 即 u 值给定 ； 即 只有 初始 条件而没有边界条件的定解问题 （ 2） 第二类）条件： 0 .即 u 的外法向导数给定 ； 即 只有边界条件而没有初值条件的定解问题 （ 3） 第三类）条件 ： 0, 0； 即 既有边值条件又有初值条件的定解问题 定义 定义在 , 上的函数 一个关系式，设 )( 有关系式 中国民航大学本科生毕业设计 （论文） )12( ) ( ) ,x v e d d 以上变换称为 换 。 其中 1i 是虚数单位。 定义 由第 n 个时间层推进到第 1n 个时间层时差分方程提供了逐点直接计算1的表达式，我们称 次差分方程为显式格式。 定义 有限差分格式在新的时间层上包含有多于一个的节点，这种有限差分格式称为隐式格式 。 定义 1 ( , ) ( , ) ( , ) , u x t u x t t u x ( , ) ( , ) ( , ) . u x t u x x t u x 称为 向前差分 。 定义 1 ( , ) ( , ) ( , ) , u x t u x t u x t ( , ) ( , ) ( , ) . u x t u x t u x x 称为 向后差分 。 定义 1 11( , ) ( , ) ( , ) ,22 u x t u x t t u x t 11( , ) ( , ) ( , ) u x t u x x t u x x 称为 中心差 分 。 定义 1 用微分方程的解代替差分方程的全部近似解，这样得到的方程两边的差就是截断误差 。 定义 给定一个适定的线性初值问题以及与其相容的差分格式，则差分格式的稳定性是差分格式收敛性的充要条件。 求解二维热传导方程的基本思想 基本思想是把连续的定解区域用有限个离散点构成网格来代替，这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数的近似；把原方程和定解条件中的微商用差商来近似，积分用积分来近似，于是原微分方程和定解条件近似的代之 以代数方程，即优先差分方程组，解此方程组就可以得到原问题在离散点上的近似解。下面是有限差分法数值计算的基本步骤 14： 中国民航大学本科生毕业设计 （论文） 将定解 区域离散 用有限差分方法求解偏微分方程问题必须把连续问题进行离散化。为此首先要对求解区域给出网格 划 分，由于求解的问题不同，因此求解区域也不尽相同。下面用例子来说明不同区域的剖分离散。并引入一些常用术语。 考虑 双曲型和抛物型方程的初值问题，求解区域是 ( , ) , 0 . D x t x t 我们在 的上半平面画出两族平 行于坐标轴的直线，把上 半 平面分成矩形网格。其交点称为节点（或网格点）。可设距离 0x ，称其为空间步长，平行线的距离按具体问题而定。可设距离 0t ，称其为时间步长。这样两族网格线可以写作 , 0 , 1 , 2 ,x x j x j h , 0 , 1 , 2t t n t n 网格节点有时记为 ),(nx 7 插值函数的选择 选择不同的插值函数对偏微分方程进行估计，可得到不同的差分方 程，进而稳定性和精度会有所不同。 用 数展开方法是最常用的方法， 本文采用的就是 数得到离散的差分方程结合差分格式构成一个稳定的差分格式。 方程组的建立 将离散后的差分方程转化为方程组的形式，便于求解。 方程组的求解 利用矩阵的解法求解方程组，再用 矩阵求解方法进行程序化，以便对以后类似的方程进行求解。隐式差分格式方程矩阵化后，得到的矩阵是严格的对角占优三对角矩阵，我们可以根据线性方程组的求解方法对其求解。其中这要应用的是追赶法，追赶法对于此类线 性方程组的求解非常方便，用 追赶法进行编程，就可以轻松实现矩阵的求解，进而解出差分方程的近似解。 中国民航大学本科生毕业设计 （论文） 二维热传导方程 网格剖分 在区域 : ( , , ) 0 , 0 , 0D x y t x X y Y t T 中，我们设二维热传导方程的初始值和边界条件如下： ( , , ) , ( , , ) ( 0 , ( , , 0 ) ( , ) , ( , , )( , , ) 0 , ( , , ) ( 0 , ua u f x y t x y t x y x y x y tu x y t x y t T（ 2 其中 a 为正常数 。 通过已知方程，建立一个关于时间和步长的函数，这样就把初始区域划分为一个网格图 。 先将定义域 : ( , , ) 0 , 0 , 0D x y t x X y Y t T 剖分为网格 ( , , ) , 0 , 1 , , , ;, 0 , 1 , , ; , 0 .h j h n x y t x j x j J J x Xy l y l J J y Y t n t n 其中 为时间步长，, 分别为 x 轴和 y 轴的空间步长。 离散方程组的构建 3 对于 二维热传导方程，利用向前差分格式 11122 n n n n nj j j j ju u u u 2 进行离散，引入记号 2 1 , 1 ,2, 1 , 12,2.n n n nx j l j l j l j ln n n ny j l j l j l j lu u u uu u u u （ 2 其中 ,( 计算值。 差分格式 1 2 222( ) ( , , ) , n n n nj l j l x j l y j lu u u ua f j x l y n tt x y（ 2 中国民航大学本科生毕业设计 （论文） 用 数展开易得差分格式（ 2截断误差为 22()t x y 。 根据稳定性分析我们知道 对于显式格式和隐式格式，在实际使用上都受到限制，因此构造每 层计算量不大的绝对稳定的格式就成为一个 较为关键 的问题。在 一 维中，隐式格式是绝对稳定并可用追赶法很容易求解。由此产生了交替方程隐式格式。它具有绝对稳定、容易求解和有相当精度的特点。 15 我们在构造微分方程（ 2隐式格式中，对 22和 22做了同样的处理，即同时在第 n 层或第 1n 层取值。为了构造一维形式的隐式格式，对二阶导数 22用 1n 层上用未知的二阶中心差分来代替，而 22则用 u 在第 n 层上用已知的二阶中心差分来代替，这样得到的方程组在仅 x 方向是隐 式的。比较容易求解，用追赶法就可以了。同理，为对称起见，在下一时间层上重复上述步骤，即又仅在 y 方向是隐式的，对 x 方向是显式的。这样相邻的两个时间层合并起来构成一个差分格式。故用多次追赶法就可以解出 注：关于追赶法的详细介绍在第 三 章中有详细介绍。 我们解向后差分格式的方程 。 令22,，则（ 2变为 : 2212222( ) ( , , )( , , ) , j l y j l j x j l y y j u a t t f j x l y n u r u t f j x l y n t（ 2 （ 2 x 方向是隐式时，代入（ 2变形为 1, 1 , 1 1 , 1 ,( 2 ) ( 2 ) ( , , ) , n n n n n n n nj l y j l j l j l j l x j l j l j lu r u u u u r u u u t f j x l y n t(2（ 2化 为矩阵形式： 中国民航大学本科生毕业设计 （论文） 2 0 0 011 2 0 020 1 2 00 0 0 1 2 ,1nr r r ux y r r r uy x y r r ry x y r ry x y （ ）（ ）（ ）（ ）1()1 1 ,1 1 ,1 01()2 1 ,1 1 ,1( , , ) , 1 1 ,1 1 ,1111 n n n nu r u u r uj x j j y jn n nu r u uj x j jf j x l y n tn n n nu r u u r uj m x j j y j 2 （ 2 y 方向是隐式时，代入（ 2可变形为 1, 1 , 1 1 , 1 ,( 2 ) ( 2 ) ( , , ) , n n n n n n n nj l x j l j l j l j l y j l j l j lu r u u u u r u u u t f j x l y n t （ 2 （ 2化 为矩阵形式 ： 1 2 0 0 011 2 0 020 1 2 00 0 0 1 2 1,nr r rx y x r r r ux x y r r rx x y r r x y （ ）（ ）（ ）（ ）1()1 1 , 1 1 , 1 01()2 1 , 1 1 , 1( , , ) . 1 , 1 , 1 1 , 1111 n n n nu r u u r ul y l l x ln n nu r u ul y l lt f j x l y n tn n n nu r u u r um l y l l y m l(2 稳定性分析 8 我们用 法 来简单 分析差分函数的稳定性。 在上式中我们使用了 中国民航大学本科生毕业设计 （论文） 2 222( ) ( , , ) , n n n nj l j l x j l y j lu u u ua f j x l y n tt x 差分格式 。 为了便于分析，我们令 ： x y h，2 为网格比，则上式改写为： 1 2 2( ) ( , , ) , n n n nj l j l x j l y j lu u a u u t f j x l y n t (2令 12 ,i k j h i k l v e e 代 入（ 2 ,11 2 ( c o s 1 ) 2 ( c o s 1 )12 v a k h a k h v 可得增长因子为： , . . . . . ( , )121( , )22121 4 ( s i n s i n )22k k kG t k k h k ( , ) 1G t k8此方法具有较高的 稳定 性 。 中国民航大学本科生毕业设计 （论文） 三章 解 二维热传导方程在 的实现 本部分主要介绍如何根据第二章所介绍二维热传导方程的理论知识进行值实现及可视化。首先对本课题中所涉及到的 令函数进行 相关介绍，然后详细介绍如何在 实现基于有限元法的二维热传导方程初边值问题的数值实现和可视化。 关知识简介 追赶法简介 2 求解三对角线性方程组是科学与工程计算中的重要问 题。 例如解常微分方程边值问题，解热传导方程以及船体数学放样中建立三次样条函数等，都会要求解系数矩阵呈三对角线形的线性方程组。而解三对角方程组的最简单方法是用追赶法，它公式简单，计算量小，所占用的存储单元少，所以在小机器上也能求解。 需要注意的是追赶法仅适用于三对角矩阵的线性方程组的求解方法，而不适用于其他类型的矩。 形如下图所示的矩阵被称为三对角矩阵 * * 0 0 0 0 0 0 0 0* * * 0 0 0 0 0 0 00 * * * 0 0 0 0 0 00 0 * * * 0 0 0 0 00 0 0 * * * 0 0 0 00 0 0 0 * * * 0 0 00 0 0 0 0 * * * 0 00 0 0 0 0 0 * * * 00 0 0 0 0 0 0 * * *00000000 *星号 （ *） 是数据，其他为零。 以按行为主序的原则转存为一维数组 Mk中，则 ,Ai j 的 对应关系为 22k i j 另一种计算方式为 ： 当 1 33 32国民航大学本科生毕业设计 （论文） 1 13 1 ( 1 , 2 , . . . , 1 ) . i i i ik i x q u x i n 思路为 ： 1 1 1 12 2 2 2 21 1 1 1 12 , , 1给 定 （ 三 对 角 方 程 组 ）且 按 行 严 格 对 角 占 优 ， ， （ ） ：n n n n nn n n i i n c x fa b c x fa b c x fa b x c b a + c b a i 过 消元，可将它化为同解方程组 1112221111111 消元过程）称为追。 1 1 1 1 1 1/ , / ,u c b q d ( ) ( 2 , . . . , 1 ) , i i i i iu c b u a i ) / ( ) ( 2 , . . . , ) . i i i i i i iq d q a b u a i , ( 1 , 2 , . . . , 1 ) . i i i ix q u x i n 回代过程）称为赶。 相关函数命令简介 1 数 13 数的最通常的应用是以下形式： 举个简单的例子：比如要计算 ),当 然可以直接用命令 y=);利 用 还可以这样来做： y 2); 另 外 这 里 的 函 数 名 字 还 可 以 是 一 个 函 数 句 柄 ， 即h=y=h,2);或者直接写成 y=); 2 数 13 在 令窗口、程序或函数中创建局部函数时，可用 点是不必将其储存为一个单独文件。在运用中有几点限制：不能调用另一个 数，中国民航大学本科生毕业设计 （论文） 能由一个 达式组成，并且只能返回一个变量 显然不允许 u,v这种形式。因而，任何要求逻辑运算或乘法运算以求得最终结果的场合，都不能应用 了这些限制，在许多情况下使用该函数非常方便