波前法及matlab实现.doc

有限元二维热传导波前法MATLAB程序 ? 二维热传导有限元 使用高斯消去法解线性方程组的二维热传导有限元程序 波前法的基本概念与算法 使用波前法解线性方程组的二维热传导有限元程序 消元过程 波前法与高斯消去法的效率之比较 小结：波前法的过去、现在和未来 波前法是求解线性方程组的一种方法，广泛用于有限元程序。它最初由英国人（？）B.M. Irons于1970在“国际工程计算方法杂志”上发表。30多年来，波前法有了不少变种。本文所用算法，采于法国人Pascal JOLY所著《Mise en Oeuvre de la Méthode des Eléments Finis》。这本书是我1993年在比利时一家书店买的，书中有一节波前法，六页纸，解释了基本概念和算法，但没有程序，也没有细节讨论。我曾花了两个半天的时间，在网上寻找波前法程序，或更详细的资料，没有找到（需要花钱才能看的文献不算）。倒是看到不少中国人，也在寻找。一些人说，波前法程序太难懂了。? 通过自己编写程序，我同意这些人的说法，确实难。我还真很少编如此耗费脑力的程序。完工之后，我曾对朋友老王说，有了算法，编程序还这么难，当初想出算法的人，真是了不起。? 现将我对波前法的理解和编程体会解说如下，供感兴趣的网友参考，也为填补网络上波前法空白。 二维热传导有限元 波前法和有限元密不可分。因而，在编写波前法程序之前，必须有个有限元程序。为了简化问题，最好是能解算一个节点上只有一个自由度的问题的有限元程序，而且要尽可能地简单。我手边现有的有限元程序都不符合这个要求。就决定先开发一个解算二维热传导问题的MATLAB有限元程序。? 二维热传导问题的微分方程是 其中 T 是温度，Kx 和 Ky 分别是 x 和 y 方向上的热传导系数，q 是热源。? 对于这样的比较经典的二阶微分方程，如何导出有限元表达式？? 这个问题，是有限元的首要问题！ 我相信，许许多多学过有限元的人，甚至每天都在用有限元的人，并不真的十分明了。? 我自己曾经就是这样。在我于2005年3月到巴西教书之前，我搞过20年有限元，其中有六年，还是在比利时一个专门开发有限元程序的公司里度过的。但是，如果您那时问我，任给一个二阶偏微分方程，例如上述热传导方程，如何导出有限元表达？说老实话，不看书，我还真就不会！? 直到2006年，我教了一遍有限元后，才弄明白（惟教书才是最好的学习）。? 其实简单极了：只需将那二阶偏微分方程，乘上一个任意标量函数，然后在某个有限单元上积分。? 请看下列推导： 即 其中， Ωe 是单元面积；Φ 是任意标量函数。? 注意在以上积分中，温度要遭受二阶偏导，这很不好。有限元的精髓，在于通过分步积分，将其中一阶偏导转移到那个任意标量函数 Φ 上，这样就可降低一阶温度偏导，改善对它的苛刻待遇。这得用到您在高等数学最后一章里学过的散度定理（Theorem of divergence）：? 其中，Г 是面积Ω的边界； （反）Δ 是梯度算子；F 和 G 是任意两个矢量函数。 ? 对于二维问题，上述散度定理可写为： 将此散度定理应用于（2）式中的第一项积分，令 得： 将此积分结果带入（3）式，得到热传导偏微分方程的弱化表达（weak form）： 所谓“弱化”，是指对温度函数的可导阶数要求降低了。在原热传导偏微分方程（1）中，温度函数必须是二阶可偏导函数，而在弱化表达（6）中，温度函数只要一阶偏导就行了。? 有限单元法就是以偏微分方程的弱化表达式为出发点的。? 现在，到了说明那个任意标量函数，Φ，是何方神圣的时候了。? 有限单元法有各种各样的变种，而最常见的，应用最广泛的，是所谓迦辽金法（Galerkin），即将这个任意标量函数，等同于，有限元的插值函数。迦辽金法的优点是可以最终形成对称劲度矩阵，从而具有较好的数值稳定性。? 我们知道，在一个有限单元中，任意一点的值，例如温度，是用节点上的温度表达的。对于一个四节点双线性单元来说，设四个节点的温度分别是T1,T2,T3,T4，则单元内任意一点的温度 T 可表达为? 其中 Ψ1，Ψ2，Ψ3，Ψ4，为插值函数，也称为形函数，定义如下： 其中 ξ 和 η 称为形坐标，取值区间为 [-1,1]。? 基于式（7），对温度的偏导数可写为： 将此二偏导数代入弱化表达式（6），该式就转化为以节点温度为变量的代数方程： 到此，为得到原偏微分方程的有限元表达，我们只需将任意标量函数，Φ，依次取为四个插值函数，Ψ1-4，就得到四个代数方程： ? 注意到式(9)-(12)中下标的规律，我们可将这四个方程合并写成矩阵的形式： 使用下标表达，式（13）可进一步缩写为： 式中对应于下标 i = 1…4 的每一个值，下标 j 取遍1…4。? 式（14）就是热传导偏微分方程（1）的四节点双线性有限