波前法及matlab实现(3)

波前法的基本概念与算法

有限元形成的线性方程组的系数矩阵，即刚度矩阵，是稀疏矩阵，也就是说，矩阵里含有许许多多的零元素。有限元网格节点数目越巨大，非零元素与零元素的比值越小，刚度矩阵越稀疏。用普通高斯消去法求解这样的线性方程组，完全不考虑矩阵的稀疏性，对大量的零元素也进行加减乘除，结果浪费了大量时间。不仅如此，应用高斯消去法，因为需要将整个刚度矩阵存在计算机内存里，所需计算机内存量与有限元网格节点未知量总数成平方的关系。以热传导问题为例。一千个节点，光存刚度矩阵，就需内存1000 x 1000 x 8 / （1024 x 1024） = 7.8 Mb。 这还没问题。但若要计算一万个节点的问题，就需要 10 x 10 x 7.8 = 780 Mb 来存刚度矩阵。内存为 1 Gb的计算机已经不能计算这样的问题了，因为微软视窗等其它系统程序还需要内存。您当然可以开始这样的计算。微软视窗发现内存不够时，会自动启动虚拟内存，也就是把硬盘上的某一块地方，当作内存来使用。您的计算仍然能进行。但是，您一定看不见那最终的计算结果，除非几个月内不断电，计算机不出毛病。原因在于，与内存相比，虚拟内存的存取时间可认为是无限长！在这种情况下，因为普通高斯消去法需要极其频繁地使用虚拟内存，它的解算时间也就无限地延长了。

但如果您在这样的计算机上运行ANSYS，或任何需要花钱买的有限元程序，就会发现，解算同样的问题，只需几分钟。您甚至可以毫无困难地解算十万个节点

的热传导问题。 秘密就在于，这些商业有限元软件，在求解线性方程组时，都尽可能地利用刚度矩阵的稀疏性。波前法就是这样一种充分考虑了刚度矩阵的稀疏性求解方程的方

法。

刚度矩阵为什么会稀疏？

因为在有限元中，一个节点的影响，只限于它所连接的那些单元。为什么？就是因为在前面，我们推导有限元时，在式（2）中，将热传导偏微分方程乘以的那个神奇函数Φ。我们说过，Φ是任意标量函数。既然是任意的，当然可以任意选取。然而我们没有“任意”地、胡乱地选取，而是居心叵测地，让它恰恰等于有限元的插值函数！而这些插值函数，恰恰只在给定单元内部非零，在单元外边一律为零！换句话说，插值函数只是些局部函数。我们让Φ等于插值函数，它也就具有了这种局部性。正是Φ的这种局部性，使得一个节点的影响，只限于它所连接的单元。有限元方法，之所以能够在计算力学领域里令人眼花缭乱的各种各样的计算方法中，独领风骚，傲视群雄，鹤立鸡群，至今几达50年之久，其全

部奥妙，皆出于此！ 为进一步考察这些影响到底有多少，我们来看下面的例子。使用前面的MATLAB有限元程序，给 Cdiv 的值输入 8, Tdiv 输入 2 ，就会生成以下网格。它将圆

周分成8份，厚度分成2份。

图中括弧内是单元号码，其余数字为节点号码。可以看出，第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单元相连，其影响也就只限于第1和第2单元。因而，在总刚度矩阵第4行当中，第1、2、5、4、8、7列的元素非零，其余18个元素为零。第2节点影响4个单元，分别是第1、9、8、16单元，因而在总刚度矩阵第2行当中，非零元素最多，达到9个，即第1、2、3、4、5、6、22、23、24列的元素非零，其余15个元素为零。如此，可想而知，总刚度矩阵的每一行

的大部分元素都是零。 现在我们要考虑怎样利用这些零元素了。在以上网格中，共有16个单元，24个节点。使用高斯消去法，会生成24 x 24 的总刚度矩阵，即有24行，24列。而使用波前法，总刚度矩阵的行数虽然不变，也是 24， 但列数仅为11（至于为什么是11 ，下面要详细讨论）。最重要的，在本例情况下，是波前法根本就不需要将 24 x 11 的总刚度矩阵存在内存中，而是存在硬盘上的。内存里，波前

法只需要放一个 11 x 11 的所谓波前矩阵就行。

那么，什么是波前矩阵呢？

就是在某一时刻，彼此发生关系的节点的刚度系数组成的矩阵。这个矩阵是方

的，其中含有最多非零元素的那一行就叫波前。 什么叫某一时刻？就是某一单元！如前面MATLAB程序所示，计算有限元刚度矩阵有两重循环，最外面那层循环，是对单元循环，即从编号为第一的单元，到编号最大的单元。在波前法中，当循环到某一单元时，在计算该单元刚度矩阵以后，还要看看哪一个节点可以消去了，也就是消元。被消元的节点，对以后其它单元刚度矩阵就不再有影响了，该节点的刚度系数就可以存入硬盘指定文件中，而这些系数就可以从波前矩阵中清除掉，以便空出位置来，存储其它节点信息。因此，一个节点可以被消元的时间，是可以用单元循环的进度来度量的。 那么，一个节点，什么时候可以消元了？就是与该节点相连的所有单元都循环到了的时候。例如，若循环从第1单元开始，当循环完了第2单元（计算单元矩阵并装配到波前矩阵中）时，第4节点就可以消元了，因为第4节点所连接的2个单元都循环到了。同理，当循环完了第3单元时，第7节点也可以消元了。而第1节点的消元要等到循环到第8单元。而第2、3节点的消元时间最迟，要

循环到最后一个单元，第16单元。 据此，可以编制一个节点消元时间表，也就是循环到了哪个单元，该节点便可以被消元了。算法很简单，就是查找每一个节点所连接的最大单元编号。程序如下：

程序二：计算节点消元时间

以上程序中，第三行的ICO变量是个两维数组，18行，4列。它的每一行代表一个单元，该行的4列给出该单元的四个节点。这段程序执行的结果，是在一维数组变量NodeDeactiveTime中定义每个节点可以消元的时间（即单元号）。此

时，NodeDeactiveTime 的值是：

8 16 16 2 10 10 3 11 11 4 12 12 5 13 13

6 14 14 7 15 15 8 16 16。 第1个数“8”代表第1节点的消元时间是8（单元）；第2个数“16”代表第2节点的消元时间是16（单元）；余类推。请注意，消元最早的节点是第4节点，时间是“2”（单元）。其次是第7节点，时间是“3”（单元）。我们下面介绍在波前

矩阵里如何消元时要用到这两个信息。

知道了各节点的消元时间，就可以计算波前矩阵的最大阶数了。程序如下：

程序三：计算波前矩阵的最大阶数

在以上程序中，第1行开了一维辅助数组，NodeInFront，用于记录每一个节点是不是在波前矩阵中，“1”表示在，“0”表示不在。第2行将我们要计算的波前矩阵的最大阶数maxFrontWidth的值初始为零。第3行开始对单元循环。对编号为 i 的单元，第4行从单元总表ICO中取出该单元的节点（本例为四个节点）；第5行将这些（四个）节点在NodeInFront中的值定为“1”，代表它们进入了波前矩阵。注意，此时，有的节点可能已经在波前矩阵中了，即它们在NodeInFront中的值已经是“1”。但这没有关系，现在只是重新再植一次“1”，再一次表示该节点在波前矩阵中。第6行计算 maxFrontWidth。就是将NodeInFront中的“1”相加，再与当前maxFrontWidth比较，谁的值更大，maxFrontWidth就取谁

的值。也就是说，maxFrontWidth的值只增不减。 第8行对 i 单元的（四个）节点循环，查找其中每个节点是不是到了消元时间（由NodeDeactiveTime给出）。 如果是的，第10行的逻辑变量deactive的值为“1”，并在第12行将该节点在NodeInFront中的值改为零。这表示该节

点被清除出波前矩阵。 这段程序结束全部循环后，便得到maxFrontWidth的值为11，就是波前矩阵的最大阶数。前面说的波前法总刚度矩阵是24 行 x 11列。其中的11，就是这样得出的。它的含义，就是在整个计算过程中，某一时刻，同时存在于波前矩阵

的节点数，其值最大为11。 以上程序，实际上模拟了节点进入和离开波前矩阵的装配和消元过程。下表给出波前矩阵中的节点随着单元循环的整个变化过程。第1列是单元号码。第2列是将该单元的刚度矩阵装入波前矩阵后，波前矩阵中的节点。第3列是这些节点的数目。第4列是此时刻可以消元的节点。第5列是消元后，将消元节点清除

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