二维热传导方程数值解及matlab实现.docx

资源描述

二维热传导方程数值解及matlab实现.docx

《二维热传导方程数值解及matlab实现.docx》由会员分享，可在线阅读，更多相关《二维热传导方程数值解及matlab实现.docx（30页珍藏版）》请在冰点文库上搜索。

二维热传导方程数值解及matlab实现.docx

第一章绪论 I

1.1 课题背景和意义 I

1.2 课题研究现状 I

1.3 课题要求 II

1.4课题解决 II

1.5 论文结构安排 III

第二章有限元法及偏微分方程解法理论分析 IV

2.1 有限元基本知识介绍 IV

2.2 求解二维热传导方程的基本思想 V

2.2.1 将定解区域离散 V

2.2.2 插值函数的选择 VI

2.2.3 方程组的建立 VI

2.2.4 方程组的求解 VI

2.3 二维热传导方程 VI

2.3.1 网格剖分 VI

2.3.2 离散方程组的构建 VII

2.3.3 稳定性分析 IX

第三章FEM求解二维热传导方程在MATLAB中的实现 XI

3.1 MATLAB相关知识简介 XI

3.1.1 追赶法简介 XI

3.1.2 相关函数命令简介 XII

3.2 FEM求解二维热传导方程在MATLAB中的实现方法 XIII

第四章数值化和可视化举例分析 XV

4.1 二维热传导方程数值求解的MATLAB实现 XV

4.1.1研究一种较为简单的情况：

4.1.2研究一种较为复杂的情况 XVI

4.2 误差分析 XVIII

结论与展望 XXI

致谢 XXIII

参考文献 XXIV

附录A:

程序清单 XXV

附录B:

外文翻译资料 XXV

第一章绪论

1.1 课题背景和意义

热传导是一种普遍存在的物理现象,它广泛存在于目前的工程应用领域中，对人们的生产和生活实践有着广泛而深刻的影响。

虽然人们对于热传导现象的本质特征已有了一个比较完整的认识，并已通过严密的数学逻辑公式推导对其物理特征进行了精确的描述，但所得到的结果往往是复杂的积分或级数表达式，其中还免不了使用某些特殊函数，由于现在数学上存在的困难，在工程中许多热传导问题还不能采用分析解法进行求解。

因此,通过一种较为简便直观的方式掌握控制和改进热量传递的方法和技术措施,无论对国民经济建设还是改善人民生活都具有重要的意义。

近年来，随着计算机技术迅速发展，数值方法已经得到广泛应用并成为有力的辅助求解工具，人们在此基础上已发展了诸如有限差分法、有限元法和边界元法等用于工程问题的求解方法。

我们知道，热传导问题的数学表达式一般为偏微分方程形式，但偏微分方程一般并没有固定有效的精确解法，因而人们广泛采用数值方法来实现热传导方程的研究、模拟和仿真。

对于热传导问题的数值计算及其可视化，原则上，可以用FORTRAN或C语言来完成这个任务，可是实际上并没有人去这样做，原因是成本太高，需要消耗更多的人力和物力，达到的效果却并不能满足人们的需求。

幸运的是，高性能数学软件（如MATLAB）与目前具有强大计算功能的个人计算机让这个问题变得简单起来。

本课题旨在尝试利用目前在理论上已经成熟的有限差分法原理结合MATLAB强大的数据处理和模拟仿真功能编写特定的程序代码对给定初边值条件的二维热传导方程温度随时间、空间的变化规律进行研究并将其可视化。

本课题将以二维热传导方程的数值解法及MATLAB实现为主线，研究论证其可行性，从而发现一种较为简便且极为有效的热传导方程数值解法和可视化的方法，意在更好的解决目前在工程和研究领域中实际存在的热传导问题，进而推动其相关领域的发展和进步。

1.2 课题研究现状

近三十年来，解偏微分方程的理论和方法有了很大的发展，而且在各个学科技术的领域中应用也愈来愈广泛，在我国，偏微分方程数值解法作为一门课程，不但在计算数学专业，而且也在其他理工科专业的研究生的大学生中开设。

同时，求解

热传导方程的数值算法也取得巨大进展，特别是有限差分法方面，此算法的特点是

在内边界处设计不同于整体的格式，将全局的隐式计算化为局部的分段隐式计算。

而且精度上更好。

目前，在欧美各国MATLAB的使用十分普及。

在大学的数学、工程和科学系科，MATLAB被用作许多课程的辅助教学手段，MATLAB也成为大学生们必不可少的计算工具，甚至是一项必须掌握的基本技能。

在我国，MATLAB在各大专院校的应用日益普遍，许多专业已把MATLAB作为基本计算工具。

在科研机构和工业界，MATLAB正得到越来越广泛的应用。

MATLAB具有强大的图形绘制功能，为科学计算和图形处理提供了很大的方便。

我们只需制定的绘图方式，再提供绘图数据，有程序指令就可以得到形象、直观的图形结果。

因此，近些年越来越多的人开始使用MATLAB来求解数值计算和图形处理技术，我们也可以绘制出热传导方程数值解的二维、三维图形，从而可以更好的理解热传导方程的意义。

1.3 课题要求

用有限元方法对如下二维热传导方程初边值问题进行数值求解，编程并进行MATLAB数值实现：

u f（x,y,t）, （x,y,t）（0,T]

u（x,y,0）（x,y）, （x,y）

u（x,y,t） 0, （x,y,t）（0,T]

2u 2u

其中 u

x2 y2

，f（x,y,t）为热源函数；为空间定义域；为的边界。

要求求出给定时间情况下温度随空间的分布变化规律。

1.4课题解决

目前，对于求解偏微分方程有很多方法，但差分法和有限元离散法式主要解决问题的两种方法。

一般来说，用差分法来解偏微分方程，解得的结果就是方程的准确解函数在这点上的近似值。

而用变分近似的方法求解，是将近似解表示成有限维子空间中基函数的线性组合。

有限元法也是基于变分原理，由于选择了特殊的基函数，使它能适用于一般的区域。

这种基函数是与区域的剖分有关的，近似解u表示为基函数的线性组合，二线性组合中的系数，又是剖分节点上u或其导数的近似值。

有关一维热传导方程的有限差分法求解的MATLAB实现[4]，现在已经解决，本文借鉴一维热传导有限差分法的MATLAB求解思想，对二维热传导方程进行转换，再

对解法编程实现，从而进一步对热传导方程进行探讨。

二维热传导方程求解在现实

生活中的应用也更加广泛，所以有很好的现实意义。

1.5 论文结构安排

第一章，绪论。

主要介绍课题的研究背景和意义，分析了热传导现象在目前工程领域的主要解决方案和发展现状，最后说明论文的主要内容和组织结构。

第二章，有限元法理论研究。

主要介绍与二维热传导问题相关的FEM理论，为有限元法在MATLAB中的实现奠定基础。

第三章，FEM在MATLAB中的实现。

主要介绍如何针对特定的一类二维热传导方程结合初边值条件编制相应的程序代码。

第四章，二维热传导方程的数值化和可视化。

主要介绍如何利用编制好的程序代码带入相应的初边值条件，进行数值化和可视化研究。

第二章有限元法及偏微分方程解法理论分析

本部分主要对二维热传导方程的有限元解法进行了完整的理论分析。

首先从有限元的基本知识开始引出求解二维热传导方程的基本思想，然后根据本思想使用交叉方向隐格式求解方法构建矩阵，并论证其稳定性。

2.1 有限元基本知识介绍

定义2.1[8]含有未知函数u（x,x, ,x,t）的偏导数的方程称为偏微分方程。

定义2.2[8]方程

1 2 n

u u u

（k ）（k ） f（x,t）,

t x x x x

1 1 n n

称热传导方程。

其中，u u（x,t）是固体的传热过程中在x处、t时刻的温度。

系数k

1 2

称为热传导系数，当k k k

（ 0）时，方程为

uf（x,t）,

2 2 2

其中

，n为维数。

x2 x2 x2

1 2 n

定义2.3[8]在特定条件下求解方程的解。

这样的条件成为定解条件。

给出了方程和定结条件，就构成了定解问题。

定义2.4[8] 一般说，边界条件有下列形式

（x,y）u（x,y）（x,y）（x,y）（x,y）,

其中为边界的外法向导数。

有如下几种特殊形式

（1）richlet（或第一类）条件：

0,即u值给定；即只有初始条件而没有边界条件的定解问题

（2）Neumann（或第二类）条件：

0.即u的外法向导数给定；即只有边界条件而没有初值条件的定解问题

（3）Robbins（或第三类）条件：

0, 0；

即既有边值条件又有初值条件的定解问题

定义2.5[8]定义在 , 上的函数vx

的一个关系式，设

v（x） dx

有关系式

v（x） 1

v（）ei（

x）dd,

以上变换称为Fourier变换。

其中i 是虚数单位。

定义2.6[8]由第n个时间层推进到第n

1个时间层时差分方程提供了逐点直接计算

un1的表达式，我们称次差分方程为显式格式。

定义2.7[8] 有限差分格式在新的时间层上包含有多于一个的节点，这种有限差分格式称为隐式格式。

定义2.8[11]

称为向前差分。

定义2.9[11]

tu（x,t）

xu（x,t）

u（x,t t）

u（x x,t）

u（x,t）,

u（x,t）.

称为向后差分。

定义2.10[11]

tu（x,t）

iu（x,t）

u（x,t）

u（x,t t）,

u（x x,t）.

tu（x,t）

u（x,t

1 t）

u（x,t

1 t）,

2 2

xu（x,t）

u（x

1 x,t）

u（x

1 x,t）.

2 2

称为中心差分。

定义2.11[11]用微分方程的解代替差分方程的全部近似解，这样得到的方程两边的差就是截断误差。

定义2.12[8]给定一个适定的线性初值问题以及与其相容的差分格式，则差分格式的稳定性是差分格式收敛性的充要条件。

2.2 求解二维热传导方程的基本思想

基本思想是把连续的定解区域用有限个离散点构成网格来代替，这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数的近似；把原方程和定解条件中的微商用差商来近似，积分用积分来近似，于是原微分方程和定解条件近似的代之以代数方程，即优先差分方程组，解此方程组就可以得到原问题在离散点上的近似解。

下面是有限差分法数值计算的基本步骤[14]：

2.2.1 将定解区域离散

用有限差分方法求解偏微分方程问题必须把连续问题进行离散化。

为此首先要

对求解区域给出网格划分，由于求解的问题不同，因此求解区域也不尽相同。

下面用例子来说明不同区域的剖分离散。

并引入一些常用术语。

考虑双曲型和抛物型方程的初值问题，求解区域是

D （x,t） x ,t0.

我们在xt的上半平面画出两族平行于坐标轴的直线，把上半平面分成矩形网

格。

其交点称为节点（或网格点）。

可设距离x0，称其为空间步长，平行线的距

离按具体问题而定。

可设距离 t 0，称其为时间步长。

这样两族网格线可以写作

xxj

jxjh,j

0,1,2 ,

t tn

nt n,n

0,1,2

网格节点有时记为（xx,tn）。

[7]

2.2.2 插值函数的选择

选择不同的插值函数对偏微分方程进行估计，可得到不同的差分方程，进而稳定性和精度会有所不同。

用Taylor级数展开方法是最常用的方法，本文采用的就是Taylor级数得到离散的差分方程结合差分格式构成一个稳定的差分格式。

2.2.3 方程组的建立将离散后的差分方程转化为方程组的形式，便于求解。

2.2.4 方程组的求解

利用矩阵的解法求解方程组，再用MATLAB对矩阵求解方法进行程序化，以便对以后类似的方程进行求解。

隐式差分格式方程矩阵化后，得到的矩阵是严格的对角占优三对角矩阵，我们可以根据线性方程组的求解方法对其求解。

其中这要应用的是追赶法，追赶法对于此类线性方程组的求解非常方便，用MATLAB对追赶法进行编程，就可以轻松实现矩阵的求解，进而解出差分方程的近似解。

2.3 二维热传导方程

2.3.1 网格剖分

在区域D:

（x,y,t）0

始值和边界条件如下：

xX,0

yY,0

t T中，我们设二维热传导方程的初

a u f（x,y,t）, （x,y,t）（0,T]

u（x,y,0）（x,y）, （x,y,t）

u（x,y,t） 0, （x,y,t）（0,T]

（2‐1）

其中a为正常数。

通过已知方程，建立一个关于时间和步长的函数，这样就把初始区域划分为一

个网格图。

先将定义域

（x,y,t）0

xX,0

yY,0 t T

剖分为网格

Dh {（xj,yh,tn）xj

jx,j

0,1, ,J,Jx X;

yl ly,l

0,1 ,J,Jy Y;tn

nt,n

0}.

其中t

T为时间步长，x X,y

Y分别为x轴和y轴的空间步长。

x y

M M

2.3.2 离散方程组的构建[3]

对于二维热传导方程，利用向前差分格式

n1 n n

2 n n

uj uj

uj1

uj uj1

t h2

对（2-1）式进行离散，引入记号

2un un

2un

un ,

xjl j

1,l jl j

1,l

2un un

2un

un .

（2‐2）

yjl j,l1 jl j,l1

其中un为差分方程在节点（j,l,n）的计算值。

差分格式

n1 n 2n 2n

ujl

ujl xujl yujl

a（）

f（j（

x,l（

y,n（t）,

（2‐3）

t x2 y2

利用Taylor级数展开易得差分格式（2-3）的截断误差为

（t x2

y2）。

根据稳定性分析我们知道对于显式格式和隐式格式，在实际使用上都受到限制，因此构造每层计算量不大的绝对稳定的格式就成为一个较为关键的问题。

在一维

中，隐式格式是绝对稳定并可用追赶法很容易求解。

由此产生了交替方程隐式格式。

它具有绝对稳定、容易求解和有相当精度的特点。

[1][5]

我们在构造微分方程（2-1）的隐式格式中，对

2u 2u

2和 2做了同样的处理，即

x y

同时在第n层或第n

1层取值。

为了构造一维形式的隐式格式，对二阶导数

2用u

在第n

1层上用未知的二阶中心差分来代替，而

2则用u在第n层上用已知的二阶

中心差分来代替，这样得到的方程组在仅x方向是隐式的。

比较容易求解，用追赶法就可以了。

同理，为对称起见，在下一时间层上重复上述步骤，即又仅在y方向是隐式的，对x方向是显式的。

这样相邻的两个时间层合并起来构成一个差分格式。

故用多次追赶法就可以解出ujl了。

注：

关于追赶法的详细介绍在第三章中有详细介绍。

我们解向后差分格式的方程。

ta ta

令rx

x2 y y2

，则（2-3）式变为:

2un 2un

）

un un1

a t（

x jl y jl

t f（j（

x,l（

y,n（t）

jl jl

r 2un r

2un

t f（j（

x,l（

y,n（t）,

（2‐4）

xx jl y y jl

（2-4）在x方向是隐式时，代入（2-2），变形为

n n n n n1 n n n

ujl

ry（uj,l1

2ujl

uj,l1）

ujl

rx（uj1,l

2ujl

uj1,l）

（t f（j（x,l（y,n（t）,

（2‐5）

（2-5）化为矩阵形式：

1（2rx

ry） ry

0 0 0

uj1n

ry 1

（2rx

ry） ry 0 0

uj2

0 ry

1（2rx

ry）ry 0

0 0 0

ry1

（2rx

nry）uj,m1

n1 n n n

uj1

rx（uj1,1

uj1,1）

ryuj0 1

n1 n n

uj2

rx（uj1,1

uj1,1）

（tf（j（x,l（y,n（t） .

（2‐6）

n1 n n n 1

uj,m1

rx（uj1,1

uj1,1）

ryujm

（2-4）在y方向是隐式时，代入（2-2），可变形为

n n n n n1 n n n

ujl

rx（uj,l1

2ujl

uj,l1）

ujl

ry（uj1,l

2ujl

uj1,l）

（t f（j（x,l（y,n（t）,

（2‐7）

（2-7）化为矩阵形式：

1（2rx

ry） rx

0 0 0

nu1l

r 1（2r r） r 0 0 n

x x y x

u2l

0 rx

1（2rx

ry）rx 0

0 0 0

rx1

（2rx

nry）um1,l

n1 n n n

u1l

ry（u1,l1

u1,l1）

rxu0l 1

n1 n n

u2l

ry（u1,l

r u

1u1,l1）

n n

u ru

（t f（j（x,l（y,n（t） ..

m1,l y（

1,l

1 1,l1）

yml

2.3.3 稳定性分析[8]

我们用Fourier方法来简单分析差分函数的稳定性。

在上式中我们使用了

（2‐8）

un un

1 2un

2un

a（）

jl jl x jl y jl

f（j（

x,l（

y,n（t）,

t x2 y2

的后差分格式。

为了便于分析，我们令：

为网格比，则上式改写为：

（t

（x （y h，

un un1

a （2un

2un）

（tf（j（x,l（

y,n（t）,

（2-9）

jl jl x jl y jl

令

代入（2-9）得

n vneik1jheik2lh,

vn 12a

（coskh

1） 2a

（coskh

1）vn1,

1 2

可得增长因子为：

G（t,k）

kh kh

.....k

（1,k2）

1 4a

（sin21

sin2

2）

2 2

G（t,k） 1[8]

此方法具有较高的稳定性。

第三章 FEM求解二维热传导方程在MATLAB中的实现

本部分主要介绍如何根据第二章所介绍二维热传导方程的理论知识进行MATLAB数值实现及可视化。

首先对本课题中所涉及到的MATLAB命令函数进行相关介绍，然后详细介绍如何在MATLAB中实现基于有限元法的二维热传导方程初边值问题的数值实现和可视化。

3.1 MATLAB相关知识简介

3.1.1 追赶法简介[2]

求解三对角线性方程组是科学与工程计算中的重要问题。

例如解常微分方程边值问题，解热传导方程以及船体数学放样中建立三次样条函数等，都会要求解系数矩阵呈三对角线形的线性方程组。

而解三对角方程组的最简单方法是用追赶法，它公式简单，计算量小，所占用的存储单元少，所以在小机器上也能求解。

需要注意的是追赶法仅适用于三对角矩阵的线性方程组的求解方法，而不适用于其他类型的矩。

形如下图所示的矩阵被称为三对角矩阵

展开阅读全文