波前法及matlab实现.

1、有限元二维热传导波前法 MATLAB程序*二维热传导有限元使用高斯消去法解线性方程组的二维热传导有限元程序*波前法的基本概念与算法使用波前法解线性方程组的二维热传导有限元程序消元过程波前法与高斯消去法的效率之比较小结：波前法的过去、现在和未来波前法是求解线性方程组的一种方法，广泛用于有限元程序。它最初由英国人(?) B.M. Iro ns于1970在国际工程计算方法杂志”上发表。30多年来，波前 法有了不少变种。本文所用算法，采于法国人Pascal JOLY所著Mise en Oeuvre de la M thode des El mtents Finis。这本书是我1993年在比利时一家书店

2、买 的，书中有一节"波前法"，六页纸，解释了基本概念和算法，但没有程序，也没 有细节讨论。我曾花了两个半天的时间，在网上寻找波前法程序，或更详细的资 料，没有找到(需要花钱才能看的文献不算)。倒是看到不少中国人，也在寻找。一些人说，波前法程序太难懂了。通过自己编写程序，我同意这些人的说法，确实难。我还真很少编如此耗费脑力 的程序。完工之后，我曾对朋友老王说，有了算法，编程序还这么难，当初想出 算法的人，真是了不起。现将我对波前法的理解和编程体会解说如下，供感兴趣的网友参考，也为填补网 络上波前法空白。二维热传导有限元波前法和有限元密不可分。因而，在编写波前法程序之前，必须有

3、个有限元程序。 为了简化问题，最好是能解算一个节点上只有一个自由度的问题的有限元程序， 而且要尽可能地简单。我手边现有的有限元程序都不符合这个要求。就决定先开 发一个解算二维热传导问题的 MATLAB有限元程序。二维热传导问题的微分方程是其中T是温度，Kx和Ky分别是x和y方向上的热传导系数，q是热源 对于这样的比较经典的二阶微分方程，如何导出有限元表达式？这个问题，是有限元的首要问题！我相信，许许多多学过有限元的人，甚至每天都在用有限元的人，并不真的十分 明了。我自己曾经就是这样。在我于 2005年3月到巴西教书之前，我搞过20年有限 元，其中有六年，还是在比利时一个专门开发有限元程序的公司

4、里度过的 。但是, 如果您那时问我，任给一个二阶偏微分方程，例如上述热传导方程，如何导出有 限元表达？说老实话，不看书，我还真就不会！直到2006年，我教了一遍有限元后，才弄明白（惟教书才是最好的学习）。其实简单极了：只需将那二阶偏微分方程，乘上一个任意标量函数，然后在某个 有限单元上积分。请看下列推导：其中，Q e是单元面积；是任意标量函数。注意在以上积分中，温度要遭受二阶偏导，这很不好。有限元的精髓，在于通过 分步积分，将其中一阶偏导转移到那个任意标量函数上，这样就可降低一阶温度偏导，改善对它的苛刻待遇。这得用到您在高等数学最后一章里学过的散度定理（Theorem of diverge n

5、ee ）:f(VF> fidJl11-f(VG)Fdft + i(n F) Ckir其中，r是面积Q的边界；（反）是梯度算子；F和G是任意两个矢量函 数。对于二维问题，上述散度定理可写为：dCi 込将此散度定理应用于（2）式中的第一项积分，令Fy "y dyGx =GV =0T.U-LL±(k 匹J dyKy dy )+砂(K匹k* rKr曲K "du、1 3x J Sy dy JIdTdx JkJ1：一 T|d<it + U dr* dnc)TaTUI Krh71' chTJ 、ddQ +fo dl> (0 qdQ =0 (6iK罟卜剽

6、uX c) ' dv .# i.z得：将此积分结果带入（3）式，得到热传导偏微分方程的弱化表达（weak form ）：所谓弱化”是指对温度函数的可导阶数要求降低了 。在原热传导偏微分方程（1） 中，温度函数必须是二阶可偏导函数，而在弱化表达（6 ）中，温度函数只要一阶偏导就行了。有限单元法就是以偏微分方程的弱化表达式为出发点的。现在，到了说明那个任意标量函数，是何方神圣的时候了。有限单元法有各种各样的变种，而最常见的，应用最广泛的，是所谓迦辽金法（Galerkin ），即将这个任意标量函数，等同于，有限元的插值函数。迦辽金法的优点是可以最终形成对称劲度矩阵，从而具有较好的数值稳定性。

7、我们知道，在一个有限单元中，任意一点的值，例如温度，是用节点上的温度表 达的。对于一个四节点双线性单元来说，设四个节点的温度分别是Tl,T2,T3,T4，则单元内任意一点的温度 T可表达为itiL1.£其中W1, W2 , W3 , W4，为插值函数，也称为形函数，定义如下:屮 I =土(1 -g)(l-】)V； = +(i + 1)0 -耳) 恢二(1十毛笊十说 斗屮 j =4(-Xi + n)斗其中E和n称为形坐标，取值区间为-1,1基于式（7），对温度的偏导数可写为：将此二偏导数代入弱化表达式（6），该式就转化为以节点温度为变量的代数方 程：到此，为得到原偏微分方程的有限元表达

8、，我们只需将任意标量函数，依次取为四个插值函数，W1-4，就得到四个代数方程：rjl|7liiL巧'Ir)Tfv0-wIIIchiJr/11dqdil+初】Itfxch屮屮加屮屮I+t):使用下标表达，式（13）可进一步缩写为:注意到式（9）-（12）中下标的规律，我们可将这四个方程合并写成矩阵的形式:Tt IfJT_ xlU, +J dr +JMH龟阳+J%q理=0rhp > rJlj出|打1湖a訂啊屮.如J州4 J'2 |t）X阪1kJdi"式中对应于下标i = 14的每一个值，下标j取遍14。式（14）就是热传导偏微分方程（1）的四节点双线性有限元表达，也

9、是我们接下来编写有限元程序的出发点。使用高斯消去法解线性方程组的二维热传导有限元程序这个程序是专为解一个特殊的热传导问题而设计的。所解问题是：已知一个无限 长圆筒，内径100毫米，外径200毫米，筒内壁表面温度1 °C，外壁绝热，求 圆筒截面上的温度分布。圆筒材料各向均质，热传导系数为1 （单位还得查一下，但也无所谓）。问题的解答很简单：均布，截面各点温度都是1°C。为什么？因为外部绝热，没有热量损失。温度只能是均布。而内壁温度为1°C，所以到处都是1°C。山址-1K> -1COG因为问题的几何图形简单，有限元网格便容易自动生成。在以下程序中，第1

10、2行至第51行，生成四节点单元的有限元网格。控制变量有两个，Cdiv和Tdiv。 前者定义沿圆周分成多少单元；后者定义沿径向即筒壁厚度方向分成多少单元。如果Cdiv = 8, Tdiv = 2，则所生成有限元网格如下（由第52行子程序DrawFEM 画出）：15"-/ / / / /X: -50>7 Z一 “第64行使用MATLAB命令syms定义了两个符号变量 ksi（即前面公式中的 E ）,eta（。在MATLAB中，可对符号变量进行代数运算，例如定义公式，求导， 积分等。第72行就利用代数运算定义了本文前面给出的四节点单元的形函数。第76和77行分别对形函数关于 ksi和

11、eta求导。第78至第99行，计算这些 导数在高斯积分点的数值。本程序中，每个单元有四个高斯积分点，也就是说， 在ksi和eta两个方向上，各有二个高斯积分点。第101至124行，根据式（14）计算单元劲度矩阵，Kelem，并将其装配入总 劲度矩阵K。本问题没有热源，所以在单元循环水平上，不需计算式（14）中 的热源积分项。第127至139行，施加边界条件。圆筒外壁是绝热条件，即法向热流等于零。在有限元中，这是自然满足的，所以式（14）中的边界热流积分项为零，不必 考虑。唯一需要考虑的，是圆筒内壁温度等于 1°C。这种温度给定的边界条件， 在数学上叫第一类边界条件。在有限元技术中，有

12、各种各样的方法施加第一类边 界条件。主要是考虑提高内存效率。鉴于本程序的目的是进一步开发波前法，不 需要仔细考虑如何更有效地施加这种温度给定的边界条件，因而所用的方法是最 简单的一种：即将内壁边界节点的各行方程，全部换为T = Tin （Tin = 1 C）°。相应地，将这些行的主元素所占据的列，乘以Tin后，移到等号右边。这样处理后，就得到待解的线性方程组：K x = F。第141行，使用本人自编的高斯消去法函数，egauss，解出为未知量x，也就 是个节点的温度，都等于1。这一行，也可直接用MATLAB命令，x = K F,取 代。我使用egauss的目的，是为下一步与波前法解方

13、程比较效率。程序一：使用高斯消去法解线性方程组的二维热传导有限元程序12% Purpoie:3ISolution of the heat conduction in a thick cyliner with FEM.4% Calling special funtions：6%dravFEILa for dravini finite tlewnt aeih6%etauM.B for the aolution of the systen of linear equations.7% Author： Yansheng Jiang8% May 20, 2007Q710%11clear all;12Cd

14、iv = inputCniaber of divisions in the circimferencial direction - ,)13Tdir = inputCnuber of divisions in the thickness direction =');14%15tie16a 100; Uimcr rtdiui. 17b - 200; MXiter radiust 18NumcroDeElcmcntos = Cdiv * Tdiv;19NumeroDeNos 二(Cdiv)*(Tdiv+l):20% Coordonnedades21X = zeros(NumeroDeNos

15、, 1):22Y = zeros(NumeroDeNos, 1):23dalpha s 2*pi/Cdiv;24dT= (b""a)/Tdiv;25delta n = Tdiv+1f26m 二 delta n (Cdiv - 1);27for i = l:Tdivil28r = a + (i-l)*dT.29nl = i;30nfin = m + nl：31X(nl:delta_n:nfin) = r*cos(0：dalpha：2*pi-dalpha);32Y(nl:delta n：nfin) = r*sin(0：dalpha:2*pidalpha):33end34% Co

16、nnect ividade de elementos35TCO = zeros(NumeroDeElementos, 4),36for i = l:Tdiv37nel = (i-l)*Cdiv * 1;38nl = i;39n2 = nl + 1;40n3 = n2 + delta n;41n4 = n3 - 1：42ICO(nel, 1:4) = nl n2 n3 n4;13ne = nel;44for j = 2:Cdiv-l4Sne = ne + 1;46IC0(ne»1:4) = IC0(ne-l, 1:4) + deltan;47end48nefin = ne + 1：49

17、IC0(nefin,l：2) = IC0(ne, 4) IC0(ne, 3)1;50ICO(nefin, 3:4) = Ln2 nlj,IFend52%Draw finite element mesh53drawFEM(NumeroDeElemento8f X Y,ICO);54%Sb Material datab6kconx1:% Coefficient ofconduction in directioo x57kcony kcxmx. % Coefficient of hmt conduction in direction y%IniciaHxar a satrix de rigid

18、71;x global Teroa independentTotalDof « NuMroDeNoa：K - 8eroo(ToUilDof9 TotalDof) #F seros(TotalDof, 1);%fundoes de interpolacaoMey*8 kai66%66% n4 oo n3%> ksi%oo70%nln27172psi = 0. 2b*( (1-kai)(lcsi) * (l-et> (lksi)<l*etA)<1-73%Dorivstiv«* of shape functions with reference to na

19、tural coordinates74nue.nodes per eloa 4.7bnuBxacsspoxnts.perelc®4:76dpaidksi » diffCpaikii)dpsidet diff(pai»«ta).G«utPointCor « 0. 577350269189626.GaussWeisht 二 1:GeuMMpoint=CaunuPointCor« GmiBsPointCor91GauasPo i ntCor, -CausaPoi ntCor92GaussPointCorf GaussPointCo

20、r93GbubhPointCor« GmjssPoinlCor,dpsidksi_nuoi zeros (nus_fBuaspoints_p«r_«】c® mai ncxiesper elm);dpsideta nuo » zeros (nunufupoi”t_per 】 nuanodetptr 】);for i = l：nun causspoints per ele«97dp«i<lksi mMB(i» ：)double (tube(dpaidkti.人(GausaPointa(i, 1). Gaus*Po

21、ints(i9 2)98ds>si<leta ma(i» ：)double (sub(dpaidetakti. eta. (GeutsPoints(£ 1) Gau«aPoxnts (i. 2);99end100%101for ttlto = 1 umeroDHEltiEtwintos102*6 口umeros global 3 do目 口口s do lesiento103Noa = ICO(eleaentoi :);104xe = X(ob);105ye - Y(Nob)；106Keletn - zeros (nunj-nodes-per-elett

22、i, nuni_nodesjpcr_eleiii).107for igauss - 1:nun saixs£point8 j)er eleii108JacohiMn109J - Ldpsidksi_niuD(iBaussi :)*xe <fps i dks i_nuiE (i gau &&, : )*ye110dps idela nunfigaussidpsideta nms(igaush ：)*yel;111invj s inv(J);112Derivatives of ahape function vith reference to global coord

23、inates113dpsidx 5 invJ(L :) * L(£tidkaLjiuio(ieau39t :)114dpsidetA nun(ifaus:);115dpidy - iiivj (2 :) * dpidki nun(igaussT :)116dDsideta nua(igauss> ：);117detj = det CD;H8Kulca 一 Kclum + (kcunx * dpidx'idx + kcony *dpsitly* +dpsidy) .1L9* det J * GaussWeight;120end121% Asaetablase to the

24、 global stiffness matrix122K(N'o»,Noe) = K(NoarXos) * Kelex;123clear J invj detJ xe ye Kelen;121end125怖 Applv essential boundarv conditions us inc ReddyJ s alcorithm, p. 173126% boudary conditions: temperature at the inner surface of cylinder = 1G127Tin = 1-128% normal thermal flux dl/dn0 a

25、t the outer aurafee.129% nodes on the inner surface *here tet£«ratur$ Tin is applied130«Mr(ii)f - 1 :del trt n：m+ I j ,131for i = 1:length(ebedof)132F F - tKebcdoffi) * Tin;J33end134F(ebcdof)二 Tin;135Ktebcdof, :) = 0;P136K(： ebedof) = 0;U371for i - l:lenith(eWof)433K (ebedof(£)te

26、bedof(£) = 1;J39end140% Solution for nodal teoiperaturewX =(K, F);442toe波前法的基本概念与算法有限元形成的线性方程组的系数矩阵，即刚度矩阵，是稀疏矩阵，也就是说，矩 阵里含有许许多多的零元素。有限元网格节点数目越巨大，非零元素与零元素的 比值越小，刚度矩阵越稀疏。用普通高斯消去法求解这样的线性方程组，完全不 考虑矩阵的稀疏性，对大量的零元素也进行加减乘除，结果浪费了大量时间。不 仅如此，应用高斯消去法，因为需要将整个刚度矩阵存在计算机内存里，所需计 算机内存量与有限元网格节点未知量总数成平方的关系。以热传导问题为

27、例。一 千个节点，光存刚度矩阵，就需内存 1000 x 1000 x 8 /( 1024 x 1024)= 7.8Mb。这还没问题。但若要计算一万个节点的问题，就需要 10 x 10 x 7.8 = 780 Mb来存刚度矩阵。内存为1 Gb的计算机已经不能计算这样的问题了，因为微 软视窗等其它系统程序还需要内存。您当然可以开始这样的计算。微软视窗发现 内存不够时，会自动启动虚拟内存，也就是把硬盘上的某一块地方，当作内存来 使用。您的计算仍然能进行。但是，您一定看不见那最终的计算结果，除非几个 月内不断电，计算机不出毛病。原因在于，与内存相比，虚拟内存的存取时间可 认为是无限长！在这种情况下，因

28、为普通高斯消去法需要极其频繁地使用虚拟内 存，它的解算时间也就无限地延长了。但如果您在这样的计算机上运行 ANSYS，或任何需要花钱买的有限元程序，就 会发现，解算同样的问题，只需几分钟。您甚至可以毫无困难地解算十万个节点 的热传导问题。秘密就在于，这些商业有限元软件，在求解线性方程组时，都尽可能地利用刚度 矩阵的稀疏性。波前法就是这样一种充分考虑了刚度矩阵的稀疏性求解方程的方 法。刚度矩阵为什么会稀疏？因为在有限元中，一个节点的影响，只限于它所连接的那些单元。为什么？就是 因为在前面，我们推导有限元时，在式(2 )中，将热传导偏微分方程乘以的那 个神奇函数。我们说过，是任意标量函数。既然是任

29、意的，当然可以任意选 取。然而我们没有任意”地、胡乱地选取，而是居心叵测地，让它恰恰等于有限 元的插值函数！而这些插值函数，恰恰只在给定单元内部非零，在单元外边一律 为零！换句话说，插值函数只是些局部函数。我们让等于插值函数，它也就具有了这种局部性。正是 的这种局部性，使得一个节点的影响，只限于它所 连接的单元。有限元方法，之所以能够在计算力学领域里令人眼花缭乱的各种各 样的计算方法中，独领风骚，傲视群雄，鹤立鸡群，至今几达50年之久，其全部奥妙，皆出于此！为进一步考察这些影响到底有多少，我们来看下面的例子。使用前面的MATLAB 有限元程序，给 Cdiv的值输入8, Tdiv输入2，就会生成

30、以下网格。它将圆 周分成8份，厚度分成2份。21图中括弧内是单元号码，其余数字为节点号码。可以看出，第1节点只与第1和第8单元相连，其影响也就只限于这两个单元。这里所说的影响，就是在刚度 矩阵中，第1和第8单元的所有节点，即第1、2、5、4、22、23节点，要发生 关系。也就是说，在总刚度矩阵（高斯消去法中的K矩阵）中，相应的行与列上的元素非零。例如在第1行当中，第1、2、5、4、22、23列的元素非零，其 余元素为零。我们知道，总刚度矩阵的列数与行数是相等的，在本列情况下，列 数等于24。在第1行当中，只有6个元素非零，其余18个元素都是零。同理， 第4节点只与第1和第2单元相连，其影响也就

31、只限于第1和第2单元。因而， 在总刚度矩阵第4行当中，第1、2、5、4、8、7列的元素非零，其余18个元 素为零。第2节点影响4个单元，分别是第1、9、& 16单元，因而在总刚度 矩阵第2行当中，非零元素最多，达到9个，即第1、2、3、4、5、6、22、23、 24列的元素非零，其余15个元素为零。如此，可想而知，总刚度矩阵的每一行 的大部分元素都是零。现在我们要考虑怎样利用这些零元素了。在以上网格中，共有16个单元，24个节点。使用高斯消去法，会生成24 x 24的总刚度矩阵，即有24 行, 24列。而 使用波前法，总刚度矩阵的行数虽然不变，也是24，但列数仅为11 （至于为什么是1

32、1，下面要详细讨论）。最重要的，在本例情况下，是波前法根本就不 需要将24 x 11的总刚度矩阵存在内存中，而是存在硬盘上的。内存里，波前 法只需要放一个11 x 11的所谓波前矩阵就行。那么，什么是波前矩阵呢？就是在某一时刻，彼此发生关系的节点的刚度系数组成的矩阵。这个矩阵是方的，其中含有最多非零元素的那一行就叫波前。什么叫某一时刻？就是某一单元！如前面 MATLAB程序所示，计算有限元刚度 矩阵有两重循环，最外面那层循环，是对单元循环，即从编号为第一的单元，到 编号最大的单元。在波前法中，当循环到某一单元时，在计算该单元刚度矩阵以 后，还要看看哪一个节点可以消去了，也就是消元。被消元的节点

33、，对以后其它 单元刚度矩阵就不再有影响了，该节点的刚度系数就可以存入硬盘指定文件中， 而这些系数就可以从波前矩阵中清除掉，以便空出位置来，存储其它节点信息。因此，一个节点可以被消元的时间，是可以用单元循环的进度来度量的。那么，一个节点，什么时候可以消元了？就是与该节点相连的所有单元都循环到 了的时候。例如，若循环从第1单元开始，当循环完了第2单元（计算单元矩 阵并装配到波前矩阵中）时，第4节点就可以消元了，因为第4节点所连接的2 个单元都循环到了。同理，当循环完了第 3单元时，第7节点也可以消元了。 而第1节点的消元要等到循环到第8单元。而第2、3节点的消元时间最迟，要 循环到最后一个单元，第

34、16单元。据此，可以编制一个节点消元时间表，也就是循环到了哪个单元，该节点便可以 被消元了。算法很简单，就是查找每一个节点所连接的最大单元编号。程序如下： 程序二：计算节点消元时间1Xudi vrUivK = /rri m u mer> dlrt:2ror i = l: N unitn J Eile nirnl o3= icotun45fnip j = 1: XhNwIics61if i > X iMk* 1ifTiinc1! t'< ns t Leri cd h chJ l' i8NodeDeactiTtTimConiHTtedNndc > - i;9

35、efidK'riidnr nd以上程序中，第三行的ICO变量是个两维数组，18行，4列。它的每一行代表 一个单元，该行的4列给出该单元的四个节点。这段程序执行的结果，是在一维 数组变量NodeDeactiveTime中定义每个节点可以消元的时间 （即单元号）。此 时，NodeDeactiveTime 的值是：8161621010311114121251313614147151581616。第1个数“8代表第1节点的消元时间是8 （单元）；第2个数“ 16代表第2节点 的消元时间是16 （单元）；余类推。请注意，消元最早的节点是第4节点，时间是“2”（单兀）。其次是第7节点，时间是“3”

36、（单兀）。我们下面介绍在波前 矩阵里如何消元时要用到这两个信息。知道了各节点的消元时间，就可以计算波前矩阵的最大阶数了。程序如下：程序三：计算波前矩阵的最大阶数I,Vmklnl iimi nriNiiiihK YrJ i;2Uhl J g汕id山i此3biri =kmtnb*丄际伙砧5i mt小1:8uud nirtUMth *Xink hi run umUtimiW nllhp；78f«j= l:hXndesTCiMttrfHlXidr = Xj|；10Lhjjiu - '“lit .kihehjin "川i|l i，h如丄 j：11il'ih iirih

37、r12、n|e liil“iiii疋Ird暫dr th13tliif吊eadlbtnd16i L.b Jr In 1i mik171MiM.iliix Diziiml： ilitUidllkiiuxl 1在以上程序中，第1行开了一维辅助数组，NodelnFront，用于记录每一个节点 是不是在波前矩阵中，“俵示在，“0”示不在。第2行将我们要计算的波前矩 阵的最大阶数maxFrontWidth 的值初始为零。第3行开始对单元循环。对编号 为i的单元，第4行从单元总表ICO中取出该单元的节点（本例为四个节点）； 第5行将这些（四个）节点在NodelnFront中的值定为“ 1”代表它们进入了波

38、前矩阵。注意，此时，有的节点可能已经在波前矩阵中了 ，即它们在NodeInFront 中的值已经是“ 1。”但这没有关系，现在只是重新再植一次“ 1”再一次表示该节点在波前矩阵中。第 6行计算 maxFrontWidth 。就是将NodeInFront中的“ 1” 相加，再与当前 maxFrontWidth 比较，谁的值更大，maxFrontWidth就取谁的值。也就是说，maxFro ntWidth 的值只增不减。第8行对i单元的（四个）节点循环，查找其中每个节点是不是到了消元时间（由NodeDeactiveTime 给出）。如果是的，第10行的逻辑变量deactive的值为“ 1,”并在第

39、12行将该节点在NodeInFront中的值改为零。这表示该节 点被清除出波前矩阵。这段程序结束全部循环后，便得到 maxFro ntWidth 的值为11，就是波前矩阵 的最大阶数。前面说的波前法总刚度矩阵是24行x 11列。其中的11，就是这样得出的。它的含义，就是在整个计算过程中，某一时刻，同时存在于波前矩阵 的节点数，其值最大为11。以上程序，实际上模拟了节点进入和离开波前矩阵的装配和消元过程。下表给出 波前矩阵中的节点随着单元循环的整个变化过程。第1列是单元号码。第2列是将该单元的刚度矩阵装入波前矩阵后，波前矩阵中的节点。第3列是这些节点的数目。第4列是此时刻可以消元的节点。第 5列

40、是消元后，将消元节点清除之后，波前矩阵里剩下的节点数，即第 3列减去第4列。可以看出，两次达到最大波前，分别当循环到第7单元和第10单元时。单朮谀机矩阵中的节点波前ft消无节点消元赶波11*544无WL2125 J X 764冒g12511 X 7 hl776A11 S It K 14 10 UBin7 5 123 It X U 171316qN6125 II H 1417 2U 16l'»10lb9?!25 11 X M 1723 艸 221110耳I 2 S IJ .s |J I?2 '22ID 1K9325 H X 1417 Hl 23 h1»无无变化

41、315 It X M 17 2U 23(>911St fi9H52 12 H H 14 J7 药 23»论*N1231 12 H 151417 Hl "9IL 1271532 IH 15 14 17 20 2314. 1514A 2 1» 21 17 戲 237门.IX§1532 24 21 HI 2J62141632 14 234儿 h 24* J3U读到这里，对照前面那个有限元网格，您要是能够理解这个表中所有数字的来龙 去脉，您就理解了波前法！剩下的，是一些技术上的细节，可以通过阅读下列全 部程序来最后彻底 抠”懂。这些技术细节里，比较难懂的，

42、是如何将单元刚度矩阵装配入波前矩阵。单元刚 度矩阵是个4阶矩阵，即4x4矩阵，代表了单元4个节点之间的相互影响。例 如，第1行里的4列元素，是单元4个节点对单元第1节点的影响；第2行里 的4列元素，是那4个节点对单元第2节点的影响；余类推。每个单元节点的局部编号都是1、2、3、4，但它们的总节点编号是不同的，而 且是唯一的。单元节点的局部编号与其总节点编号的对应关系，在本文的程序 中，由二维数组ICO给出。在用普通高斯消去法解方程时，任给一个单元刚度 矩阵，我们可以通过ICO表，查到该单元4个节点的总编号。而节点总编号与 总刚度矩阵K的行与列之间具有对应的关系。第1节点占据第1行第1列,第2节

43、点占据第2行第2列，第n节点占据第n行第n列。所以，使用 普通高斯消去法解方程时，把单元刚度矩阵装配进总刚度矩阵的算法很简单，女口 本文第一个程序的第122行所示，只一行程序便可实现。打个比方，这种情况 下的总刚度矩阵K，就像占地广阔的大观园，有许多房子，红楼梦里的所 有人等，例如金陵十二钗，各有各的住所，还有许多空地（相当于K里的零元素），可供黛玉葬花，姐妹们赏月吟诗。波前法的装配算法则要复杂得多。因为波前矩阵小，不能同时装下所有节点的所 有元素。这种情况，就如同旅店。世界上所有旅店是住不下装得下世界上所有居 民的，旅客必须有进有出才行。节点进出波前矩阵，就如同旅客进出旅馆。有的 旅客住的时

44、间长，有的短。节点在波前矩阵呆的时间也一样，有长有短。随着节 点的进进出出，波前矩阵的每一行和列都可能先后被不同节点占用，就像旅馆里 的每个房间都会被不同旅客先后住过。所以，波前矩阵的行与列，与总节点编号 之间，不再有像使用普通高斯消去法解方程时那样的一一对应的固定关系。单元 矩阵的某一行、某一列的元素，应该放到波前矩阵的哪一行、哪一列，是动态分配的。有两种可能：如果某节点已经存在于波前矩阵中，那么就把该节点在单元 矩阵中的元素加到波前矩阵的相应元素上（因而需要知道它在哪里）；反之，如果某节点还没进入过波前矩阵，那么就在波前矩阵给它分配一个自由的行与列。 在下面的程序中，我们使用一维数组 fr

45、eeLines来记录波前矩阵每一行（同时也 是每一列。我们假定行数等于列数，也就是说，一个节点占据了第n行，它也就占据了第n列）的状态：0表示自由；大于0表示已占用（即占用该行的节 点）。我们还使用一维数组GloballD2Frontld 来记录每一节点在波前矩阵占据 的行数（=列数）。也就是说，给出一个节点的总编号，就能在GlobalID2Frontld 中查到它在波前矩阵中的位置。这个位置如果是零，就表示该节点还没进入过波 前矩阵。这些算法的具体实施可见下列程序中的第148-168行。程序四：使用波前法解线性方程组的二维热传导有限元程123aSolution of the hr«

46、*t conduction in a thick cyliiier with FEM and4务a frun tai nitlhvd solver for IIk(4 Linear equnlioifih5f < Calling special funtions:6ftdra EM-ni 1'or drjwing LiihIl1 tknicnl nn sii7f<uriFrrel.ine.m tor ctlinji a fivr liz# in Ihr rrnntnl niiUiix8f< Audhiri dii>lkcn$> Jiiuii!9% June

47、 l 2UG71011dear oU;12Cdiv = inpuh1 number of dhKiutii in the ciniiEiftf*nciiiJ direclion =T;13while < Hit <4'Hdispi 'tuu Miiall number of dh ision?* in thv circunlrri'nciiil dircclwn' >15Ctlii 二 in|nit ntiniber' iTdiviMniis in Iht cirrnnlien*rinal diivrlion = " :J6e

48、nd17I'dh = inpiHrniiniber rrf divisions in rhe radial diterlion = " i：1S* hik Tdh < 219(ikpisiuall number “f dikinns in th* IhkIkJ direct20Tdii ='ninnK r nTdh inions in iht racial dirx-cUoti = h 匕21end22tic23a = 100： ft JntxT ludius mni21l> = 200： %<hiler radius nunZ5MumeroDeEl

49、einfn1<K = (."div Irtiv；26NinncroDeN<K = (Mb)*(T<lh+h；27X = zeroMMumroDeNo 1);28Y zrroyumerDeMosl);dalphj - ipi/Cdi:30<IT - i h-；ii/Tilh ；31ddtu n = 1 dh +1;m32m = rltbllji ii r < ilix - 11：for i - l:TJis+l31r = a + i-1)edT;3 bnl = i;36nftn = m + nl；37Xi nl:dcJtn n:nfin 1 - r*cck

50、<< 1:dttlpliailpUdnlplin 1;38Y(nl:dr1iM ntnl1nI = t sintltedalphiuJ11 pi-dlphau39end40R O = zrrusi iin)ciT)lkI'luntLiitip.,4);41for i = l:Tilh42nl = |bl r 17liv + 1；43ill =144n2 = n 1 + 1：45n3 = nl + ddta n;46n4 = n3 1;47ICCXnel J:4) > |nl n2 n3 n4|:48I1C = IK* 1:49forj = 2:Cdiv-I|SOne

51、= nc + 1;51ICO( ned:4) = ICO(ne-IJ:4) + delta n;52rncl53nelin = ne + 1;54ICO< ncfin,l:2) = ICO nc<4> ICO nc3>?55IC(hnelin3：4) = |n2 nl|:|56end57drawt'EMlXunwroDeElemenfosIX Y|JCOh58B(co(le = zvr<iM Xiimcrol)vN< 1);59BCaluc = zcroMunivroLkNos,l);60Tins 1：|61essential = 1:62inncr

52、>urfacc mHk*s = l:dclkcn:m-*-l;63BCcodednner suiTaw nodes) = essential;64R C al nd i n ncr.Mi rfact-.n< xics) = Tin:6bkconx =1:66kcony = kconx;67syms ki eta;68psi = 0.25*|<l-ksb*( 1-eta) (l+ksi>*( l-eta> < l4-ksi)*< I+etal (l-ksi)a<69dpsidksi =170dpsidcki =(IifT<|>shvl&

53、lt;n:71GaussPointCor = 0.577350269189626：|72(iaiixWeijihi = 1;73GaussPointss -GaussPointCor, -GaussPointCor74(laussPointCon -(aussPoinlC'or75GaussPointCon Gauss Point Cor76(ZaussPoiniCor. C>aussPoint(4>r|;77iHiin ii(><lcs|HT.ckiii = 4:78num gausspoints |)cr clem = 4;79dpsidksLnum = z

54、troM ntiin jiaussp<>iiits|M r ckm. inim n(Hlt-s|KT clcmi：80dpsideta num = zvroMnum gausspoints per elenu num nixies per ekin)：81for i = 1:num jiaiiss|H>ints prr elcin82dpsidksijHiiniL：)=doublc4subs(dp&idksiJkMjcta4GaussPoinU(i.l)tGausftPointMi«2n)h83dpsideta numli.:)二doubldsubsi dpsideta, ksiArtaJJGaussPointM i.ll.GaussPointsl84end85maxElemPerNode = 4;66NocleLk-activvl iiik* = zt ros(Xumcrul)v<A,l h87for i = 1:NumeroDe Elements88Noss ICO