波前法及matlab实现(2)

式中对应于下标 i = 1…4 的每一个值，下标 j 取遍1…4。

式（14）就是热传导偏微分方程（1）的四节点双线性有限元表达，也是我们接

下来编写有限元程序的出发点。

使用高斯消去法解线性方程组的二维热传导有限元程序

这个程序是专为解一个特殊的热传导问题而设计的。所解问题是：已知一个无限长圆筒，内径100毫米，外径200毫米，筒内壁表面温度 1°C，外壁绝热，求圆筒截面上的温度分布。圆筒材料各向均质，热传导系数为 1 （单位还得查一下，但也无所谓）。问题的解答很简单：均布，截面各点温度都是 1°C。为什么？因为外部绝热，没有热量损失。温度只能是均布。而内壁温度为 1°C，所

以到处都是1°C。

因为问题的几何图形简单，有限元网格便容易自动生成。在以下程序中，第12行至第51行，生成四节点单元的有限元网格。控制变量有两个，Cdiv 和 Tdiv。 前者定义沿圆周分成多少单元；后者定义沿径向即筒壁厚度方向分成多少单元。如果 Cdiv = 8，Tdiv = 2， 则所生成有限元网格如下（由第52行子程序DrawFEM

画出）：

第64行使用MATLAB命令 syms 定义了两个符号变量 ksi(即前面公式中的 ξ),eta(η)。在MATLAB中，可对符号变量进行代数运算，例如定义公式，求导，积分等。第72行就利用代数运算定义了本文前面给出的四节点单元的形函数。第76和77行分别对形函数关于 ksi 和 eta 求导。第78至第99行，计算这些导数在高斯积分点的数值。本程序中，每个单元有四个高斯积分点，也就是说，

在 ksi 和 eta 两个方向上，各有二个高斯积分点。 第101至124行，根据式（14）计算单元劲度矩阵，Kelem，并将其装配入总劲度矩阵 K。 本问题没有热源，所以在单元循环水平上，不需计算式（14）中

的热源积分项。 第127至139行，施加边界条件。圆筒外壁是绝热条件，即法向热流等于零。在有限元中，这是自然满足的，所以式（14）中的边界热流积分项为零，不必考虑。唯一需要考虑的，是圆筒内壁温度等于 1°C 。这种温度给定的边界条件，在数学上叫第一类边界条件。在有限元技术中，有各种各样的方法施加第一类边界条件。主要是考虑提高内存效率。鉴于本程序的目的是进一步开发波前法，不需要仔细考虑如何更有效地施加这种温度给定的边界条件，因而所用的方法是最简单的一种：即将内壁边界节点的各行方程，全部换为 T = Tin (Tin = 1°C)。相应地，将这些行的主元素所占据的列，乘以Tin后，移到等号右边。这样处理后，

就得到待解的线性方程组：K x = F。 第141行，使用本人自编的高斯消去法函数，egauss，解出为未知量 x， 也就是个节点的温度，都等于 1。这一行，也可直接用MATLAB命令，x = K \\ F, 取

代。我使用egauss的目的，是为下一步与波前法解方程比较效率。

程序一：使用高斯消去法解线性方程组的二维热传导有限元程序

百度搜索“77cn”或“免费范文网”即可找到本站免费阅读全部范文。收藏本站方便下次阅读，免费范文网，提供经典小说综合文库波前法及matlab实现(2)在线全文阅读。