材料工程基础课程大作业

安徽建筑大学材料与化学工程学院

《材料工程基础》

课程大作业

专业：无机非金属材料工程

班级：无机非一班

姓名：

学号：

2017 年 5 月

目录

1 计算模拟 (3)

1.1 导热问题数值解法 (3)

1.1.1概述 (3)

1.1.2 数值解法原理 (3)

2数值模拟案例 (9)

2.1MALTAB PDE-tool工具箱 (9)

2.2计算案例 (9)

2.3模拟具体步骤 (10)

2.3.1在MATLAB命令窗口，输入命令 >> pedtool (10)

2.3.2 建立几何模型 (11)

2.3.3 设置边界条件 (12)

2.3.4 PDE方程的选择 (13)

2.3.5 网格划分 (14)

2.3.6 绘温度场分布图 (14)

1 计算模拟

1.1 导热问题数值解法

1.1.1概述

随着计算机的普及应用和性能的不断改善，以及相关的数值计算方法的发展和应用程序的开发，传热学数值计算方法作为数值求解传热问题的有效工具也得到了相应的发展。利用计算机求解传热学问题愈来愈受到人们的普遍重视，而且在计算复杂传热问题中显示出它的优越性，因而成为传热学的一个重要的分支。数值传热的相关内容也很自然地成为工程类学生学习传热学课程的不可缺少的部分。

1.1.2 数值解法原理

(1)建立控制方程及定解条件

对于一些更为复杂的导热问题，如复杂的几何形状和边界条件以及物性变化较大的情况，分析求解往往很复杂或者根本不可能。此时求解问题的唯一途径是利用数值分析的办法获得数值解

数值求解通常是对微分方程直接进行数值积分或者把微分方程转化为一组代数方程组再求解。这里要介绍的是后一种方法。如何实现从微分方程到代数方程的转化又可以采用不同的数学方法，如有限差分法、有限元法和边界元法等。

简要地介绍用有限差分析方法从微分方程确立代数方程的处理过程。有限差分法的基本思想是把原来在时间和空间坐标中连续变化的物理量，如温度、压力、速度和热流等。用有限个离散点上的数值集合来近似表示。有限差分的数学基础是用差商代替微商，导数。而几何意义是用函数在某区域内的平均变化率代替函数的真实变化率。

在右图中可以看出有限差分

表示的温度场与真实温度场的区别。图中用T0、T1、T2…表示连续的温度场T；Δx为步长，它将区域的x方向划分为有限个数的区域，Δx0、Δx1、Δx2…，它们可以相等，也可以不相等。当Δx相等时，T1处的真实变化率a可以用平均变化率b、c或d来表示，其中b、c和d分别表示三种不同差分格式下的温度随时间的变化率，即:

b为向后差分格式：

c为向前差分格式:

d为中心差分格式:

这种差分格式也可以推广到高阶微商的情形。对于二阶微商的差分格式可以在一阶差分格式的基础上得出:

(2)区域离散化

这里我们以一个矩形长柱体的非稳态导热过程为例来讨论区域离散化问题。如果不考虑矩形长柱体长度方向上的温度变化，那么它是一个二维非稳态导热问题，如图，

上图表示了长柱体矩形截面上区域离散化的情况。图中可见，对于给定的

空间区域，在x方向上的步长为Δx，在y方向上的步长为Δy，用它们作为空间尺度可以将矩形区域划分成纵横交错的网格，交点称为节点。然后以节点为中心，在两个节点的中心处划分界限，定出节点的控制面积，对于三维情况则为控制体积或控制容积，因而常在一般意义上称之为节点的控制体。

控制体的形状是随着坐标系的不同而改变的，这里的控制体是一个个的矩形面积。网格的步长在每一个方向上可以均匀划分，也可以不均匀的划分。因此，选用不同的步长和不同的划分方法，可以将同一区域划分出不同大小、不同数目的控制区域，以及不同数目的节点数。显然，随着步长的不断减小，节点数目的不断增加，由节点温度表示的离散的温度场就会更加接近连续的温度场，但计算工作量也会随之增加。在时间方向上离散化的步长常用Δη来表示，Δη的选取也是可大可小的，也可以随时间的进程而变化。显然，无限小的时间步长Δη亦会使得离散温度变化接近连续的温度改变，但随之而来的是相应的计算工作量的增加。

（3）建立节点物理量的代数方程

控制体热平衡法建立节点方程的过程是将能量守恒方程应用于控制体，建立该节点与周围节点之间的能量平衡关系式，再利用傅立叶的导热定律，最后获得控制体节点温度与周围节点温度之间的关系式。

考察右图中的节点P及其控制体，由能量平衡关系应有

式中，Q W、Q E、Q S和Q N分别为邻近节点W、E、S和N通过传导方式传给节点P的热流量；QV为单位时间控制体内热源的发热量，ΔΕ为控制体单位时间内热

能的增加量。由导热傅立叶定律，在线性温度分布的假设下，时刻K周围节点传给节点P的热流量分别为:

以及控制体的发热流量

（q V为内热源强度，即单位时间单位体积的内热源发热量。） 控制体单位时

间的内能增加量为

，前者为时间上的向前差分，而后者为时间上向后差分。以上关系式中温度T的上标为所在时刻，下标为所在空间位置。

将以上关系式一并代入方程中，且假设Δx=Δy，经整理可以得出二维非稳态导热问题的内节点的两种差分格式的差分方程，即

比较上面两种差分格式可以看出，显示差分格式最突出的优点是节点温度表达式的右边只涉及K时刻的节点温度值，那么只要知道这一时刻周围节点的温度值就可以求出该节点的下一时刻的温度值，而隐示差分格式却相反，温度表达式的两端都是同一K时刻的节点温度值，这就意味着必须同时计算同一时刻所有节

点的温度值，即必须联立求解K时刻所有节点的差分方程组，增大计算工作量是显而易见的。

（4）设立迭代初场

对传热问题的有限差分法中主要采用迭代法，需要对被求得温度场预先假定一个解，称为初场，并在求解过程中不断改进。

（5）求解代数方程组

由上面的讨论可以看出，对应于离散温度场的每一个节点均可以列出相应的差分方程，这样就可以得出与节点数目相同的一组代数方程组。当联立求解这个代数方程组时，最后就可以得出每一个节点的温度值。一般情况下，差分方程组是线性代数方程组，而线性代数方程组是可以用直接法和迭代法求解的。常用的直接法有高斯消元法、列主元素消去法和矩阵求逆法，而迭代法常用的有高斯-赛德尔迭代和超（欠）松弛迭代。下面我们仅介绍迭代法求解代数方程组的过程。

现有一组线性方程组，

迭代求解该方程组的思路为，寻找一个由T1，T2，…，T n组成的列向量，使

其收敛于某一个极限向量（T1*, T2*

,… Tn*）且该极限向量就是该方程的精确

解。当这个线性代数方程组的系数项aii≠0（i=1,2,…,n）时，可将其改写成迭代形式，有:

以上各式可以用一个通用的形式来表示，

利用上式就可以进行迭代求解了，其步骤是，合理选择（假设）各节点的初

始温度，将其作为第零次迭代的近似温度值，记为Ti (0)

(i=1,2,…,n);将Ti(0)代入上

式的右端，得到第一次迭代的近似值Ti (l)

，之后将Ti

(l)

再代入上式的右端,则得出

第二次的近似值Ti (2)

,如此反复进行下去,直至进行到K次,使相邻的两次近似解

Ti(K+!)和Ti(K) (i=1,2,…,n)之间的偏差小于预先设定的小量ε时,即满足∣Ti(K+1)-Ti(K)∣≤ε(i=1,2,…,n)或∣(Ti(K+1)-Ti(K))/Ti(K)∣≤ε(i=1,2,…,n)。此时各节

点的温度值[T1(K)

, T2(K),…, Tn(K)]已经有足够的精度用来表示代数方程组的解,从

而可以结束方程求解的迭代过程。

从上述的迭代过程不难发现，当我们用第零次迭代值去进行第一次迭代

时,Ti (1)

的值已经不断地产生出来,当计算Tr

(1 )

时,到r-1的Ti(1)已经求出。如果此

时在计算Tr(1 )时涉及到的Ti(0) (i=1,2,…,r-1)全部用已求出的Ti (1)

(i=1,2,…,r-1)代

替,这势必会加快迭代收敛的速度。这种改进后的迭代方法被称之为高斯-赛德尔迭代法。高斯-赛德尔迭代法的方程组迭代形式为:

归纳起来，高斯-赛德尔迭代法的求解步骤可表述为，将代数方程组写成迭代形式，设初始值经迭代得出节点新值，有新值则去掉旧值，不断以新换旧，且在迭代过程中应用，在迭代获得满足给定精度的节点温度值后结束方程组的迭代。

2数值模拟案例

2.1MALTAB PDE-tool工具箱

MALTAB PDE-tool工具箱为我们提供了一个功能强大，应用灵活的二维有限元偏微分方程问题求解环境，图形用户界面简单，直观，操作方便。我们使用MATLAB PDE-tool偏微分方程工具箱来求解导热问题。

2.2计算案例

规格为360*240*190（单位mm）的煤矸石自保温砌块的俯视图，如下。

假设环境温度为35℃，用此保温砌块砌墙，室内温度为10℃，热流方向为平行纸面向下，此砌块的热工计算值如下，

λ—材料导热系数, W/(m·K)

R空=0.18----------空气热阻（厚度为40mm）

δ—厚度（mm）

R—热阻（k/w）

K—传热系数（W/m2·K）

R1=R3=R5=R7 （k/w）

R2=R4=R6 （k/w）

R1=0.0375/(0.2*0.39*0.19)=2.53 （k/w）

R2a=0.03/(0.2*0.03*0.19)=26.32 （k/w）

R2b=0.18/(0.11*0.19)=8.61 （k/w）

R2d=0.18/(0.09*0.19)=10.53 （k/w）

R2f=0.18/(0.07*0.19)=13.53 （k/w）

R2a=R2c=R2e=R2g （k/w）

1/R2=1/R2a+1/R2b+1/R2c+1/R2d+1/R2e+1/R2f+1/R2g=4/R2a+1/R2b+1/R2d+1/R2f=0.437 （k/w）R2=2.89 （k/w）

R=R1+R2+R3+R4+R5+R6+R7=18.79（k/w）

R T=R*0.39*0.19=1.39（m2k/w）

K=1/R T=0.72 （W/m2·K）

2.3模拟具体步骤

2.3.1在MATLAB命令窗口，输入命令 >> pedtool

2.3.2 建立几何模型

把偏微分方程工具箱调至，点击截图最左边

的正方形小方块，建立几何模型。如下图

2.3.3 设置边界条件

点击，设置边界条件。

上图最前边红色线条（砌块垂直于热流方向的外轮廓边界线），双击边线，设置为狄克雷Dirichlet边界条件,如图。h =1;r=25.

上图最后边红色线条（砌块垂直于热流方向的外轮廓边界线），双击边线，设置为狄克雷Dirichlet边界条件,如图。h =1;r=10

上图平行于热流方向的外轮廓边界线和孔的边界条件，双击，都设置为Neumann边界条件。热流g=0和导热系数q=0.

2.3.4 PDE方程的选择

点击图片中的PDE按钮,设置PDE方程格式。

导热系数K=0.2 W/(m·K) ,

砌块内热源Q=0W,

空气对流传热系数h=10W/(m2.℃)

外界温度Text=25℃

2.3.5 网格划分

使用这两个符号划分网格，左边三角形是初步划分，右边为再次划分。效果如图。

2.3.6 绘温度场分布图

点击图标，勾选如下参数绘二维温度场分布图，

2.3.6.1二维温度场分布图

2.3.6.2三维温度场分布情况

设置如下参数

得到如下图像

正面

背面Colormap设置为jet和Prism，如下图，

得到如下图

设置一下参数

设置如下参数

《材料工程基础》课程大作业得到以下图像。