彈性力學(xué)平面問題的有限元法

1、Mmm3 彈性力學(xué)平面問題的有限元法本章包括以下的內(nèi)容：3.1彈性力學(xué)平面問題的基本方程3.2單元位移函數(shù)3.3單元載荷移置3.4單元剛度矩陣3.5單元剛度矩陣的性質(zhì)與物理意義3.6整體分析3.7約束條件的處理3.8整體剛度矩陣的特點與存儲方法3.9方程組解法3.1彈性力學(xué)平面問題的基本方程彈性力學(xué)是研究彈性體在約束和外載荷作用下應(yīng)力和變形分布規(guī)律的一門學(xué)科。在彈性力學(xué)中針對微小的單元體建立基本方程，把復(fù)雜形狀彈性體的受力和變形分析問題歸結(jié)為偏微分方程組的邊值問題。彈性力學(xué)的基本方程包括平衡方程、幾何方程、物理方程。彈性力學(xué)的基本假定如下：1）完全彈性，2）連續(xù)，3）均勻，4）各向同性，5）小

2、變形。3.1.1基本變量彈性力學(xué)中的基本變量為體力、面力、應(yīng)力、位移、應(yīng)變，各自的定義如下。體力體力是分布在物體體積內(nèi)的力，例如重力和慣性力。面力面力是分布在物體表面上的力，例如接觸壓力、流體壓力。應(yīng)力 物體受到約束和外力作用，其內(nèi)部將產(chǎn)生內(nèi)力。物體內(nèi)某一點的內(nèi)力就是應(yīng)力。圖3.1如圖3.1假想用通過物體內(nèi)任意一點p的一個截面mn將物理分為、兩部分。將部分撇開，根據(jù)力的平衡原則，部分將在截面mn上作用一定的內(nèi)力。在mn截面上取包含p點的微小面積，作用于面積上的內(nèi)力為。令無限減小而趨于p點時，的極限S就是物體在p點的應(yīng)力。 應(yīng)力S在其作用截面上的法向分量稱為正應(yīng)力，用表示；在作用截面上的切向分量

3、稱為剪應(yīng)力，用表示。顯然，點p在不同截面上的應(yīng)力是不同的。為分析點p的應(yīng)力狀態(tài)，即通過p點的各個截面上的應(yīng)力的大小和方向，在p點取出的一個平行六面體，六面體的各楞邊平行于坐標(biāo)軸。圖3.2將每個上的應(yīng)力分解為一個正應(yīng)力和兩個剪應(yīng)力，分別與三個坐標(biāo)軸平行。用六面體表面的應(yīng)力分量來表示p點的應(yīng)力狀態(tài)。應(yīng)力分量的下標(biāo)約定如下：第一個下標(biāo)表示應(yīng)力的作用面，第二個下標(biāo)表示應(yīng)力的作用方向。，第一個下標(biāo)x表示剪應(yīng)力作用在垂直于X軸的面上，第二個下標(biāo)y表示剪應(yīng)力指向Y軸方向。正應(yīng)力由于作用表面與作用方向垂直，用一個下標(biāo)。表示正應(yīng)力作用于垂直于X軸的面上，指向X軸方向。 應(yīng)力分量的方向定義如下： 如果某截面上的外

4、法線是沿坐標(biāo)軸的正方向，這個截面上的應(yīng)力分量以沿坐標(biāo)軸正方向為正；如果某截面上的外法線是沿坐標(biāo)軸的負(fù)方向，這個截面上的應(yīng)力分量以沿坐標(biāo)軸負(fù)方向為正。剪應(yīng)力互等：物體內(nèi)任意一點的應(yīng)力狀態(tài)可以用六個獨立的應(yīng)力分量、來表示。位移位移就是位置的移動。物體內(nèi)任意一點的位移，用位移在x，y，z坐標(biāo)軸上的投影u、v、w表示。應(yīng)變物體的形狀改變可以歸結(jié)為長度和角度的改變。各線段的單位長度的伸縮，稱為正應(yīng)變，用表示。兩個垂直線段之間的直角的改變，用弧度表示，稱為剪應(yīng)變，用表示。物體內(nèi)任意一點的變形，可以用六個應(yīng)變分量表示。3.1.2平面應(yīng)力和平面應(yīng)變問題彈性體在滿足一定條件時，其變形和應(yīng)力的分布規(guī)律可以用在某一

5、平面內(nèi)的變形和應(yīng)力的分布規(guī)律來代替，這類問題稱為平面問題。平面問題分為平面應(yīng)力問題和平面應(yīng)變問題。1） 平面應(yīng)力問題設(shè)有很薄的等厚薄板，只在板邊上受到平行于板面并且不沿厚度變化的面力，體力也平行于板面且不沿厚度變化。圖3.3設(shè)板的厚度為t，在板面上：， ， 由于平板很薄，外力不沿厚度變化，因此在整塊板上有， ， 剩下平行于XY平面的三個應(yīng)力分量未知。2）平面應(yīng)變問題設(shè)有很長的柱形體，支承情況不沿長度變化，在柱面上受到平行于橫截面而且不沿長度變化的面力，體力也如此分布。圖3.4以柱體的任一橫截面為XY平面，任一縱線為Z軸。假定該柱體為無限長，則任一截面都可以看作對稱面。由對稱性，由于沒有Z方向的

6、位移，Z方向的應(yīng)變。未知量為平行于XY平面的三個應(yīng)力分量，物體在Z方向處于自平衡狀態(tài)。平衡方程彈性力學(xué)中，在物體中取出一個微小單元體建立平衡方程。平衡方程代表了力的平衡關(guān)系，建立了應(yīng)力分量和體力分量之間的關(guān)系。對于平面問題，在物體內(nèi)的任意一點有，（3-1）幾何方程由幾何方程可以得到位移和變形之間的關(guān)系。對于平面問題，在物體內(nèi)的任意一點有，（3-2）剛體位移由位移u=0，v=0可以得到應(yīng)變分量為零，反過來，應(yīng)變分量為零則位移分量不為零。應(yīng)變分量為零時的位移稱為剛體位移。剛體位移代表了物體在平面內(nèi)的移動和轉(zhuǎn)動。由可以得到剛體位移為以下形式，由可得，將代入可得，積分后得到，得到位移分量，當(dāng)時，物體內(nèi)

7、任意一點都沿x方向移動相同的距離，可見代表物體在x方向上的剛體平移。當(dāng)時，物體內(nèi)任意一點都沿y方向移動相同的距離，可見代表物體在y方向上的剛體平移。當(dāng)時，可以假定，此時的物體內(nèi)任意一點P（x，y）的位移分量為，P點位移與y軸的夾角為，P點合成位移為，r為P點到原點的距離，可見代表物體繞z軸的剛體轉(zhuǎn)動。3.1.5物理方程彈性力學(xué)平面問題的物理方程由廣義虎克定律得到。1）平面應(yīng)力問題的物理方程（3-3）平面應(yīng)力問題有，2）平面應(yīng)變問題的物理方程（3-4）平面應(yīng)變問題有，在平面應(yīng)力問題的物理方程中，將E替換為、替換為，可以得到平面應(yīng)變問題的物理方程；在平面應(yīng)變問題的物理方程中，將E替換為、替換為，可

8、以得到平面應(yīng)力問題的物理方程。圖3.5求解彈性力學(xué)平面問題，可以歸結(jié)為在任意形狀的平面區(qū)域內(nèi)已知控制方程、在位移邊界上約束已知、在應(yīng)力邊界上受力條件已知的邊值問題。然后以應(yīng)力分量為基本未知量求解，或以位移作為基本未知量求解。如果以位移作為未知量求解，求出位移后，由幾何方程可以計算出應(yīng)變分量，得到物體的變形情況；再由物理方程計算出應(yīng)力分量，得到物體的內(nèi)力分布，就完成了對彈性力學(xué)平面問題的分析。3.2單元位移函數(shù)根據(jù)有限元法的基本思路，將彈性體離散成有限個單元體的組合，以結(jié)點的位移作為未知量。彈性體內(nèi)實際的位移分布可以用單元內(nèi)的位移分布函數(shù)來分塊近似地表示。在單元內(nèi)的位移變化可以假定一個函數(shù)來表示

9、，這個函數(shù)稱為單元位移函數(shù)、或單元位移模式。對于彈性力學(xué)平面問題，單元位移函數(shù)可以用多項式表示，（3-5）多項式中包含的項數(shù)越多，就越接近實際的位移分布，越精確。具體取多項，由單元形式來確定。即以結(jié)點位移來確定位移函數(shù)中的待定系數(shù)。圖3.6如圖3.6所示的3結(jié)點三角形單元，結(jié)點I、J、M的坐標(biāo)分別為、，結(jié)點位移分別為、。六個節(jié)點位移只能確定六個多項式的系數(shù)，所以3結(jié)點三角形單元的位移函數(shù)如下，（3-6）將3個結(jié)點上的坐標(biāo)和位移分量代入公式（3-6）就可以將六個待定系數(shù)用結(jié)點坐標(biāo)和位移分量表示出來。將水平位移分量和結(jié)點坐標(biāo)代入（3-6）中的第一式，寫成矩陣形式，（3-7）令，則有（3-8），A為

10、三角形單元的面積。T的伴隨矩陣為，（3-9）令（3-10）則（3-11）同樣，將垂直位移分量與結(jié)點坐標(biāo)代入公式（3-6）中的第二式，可得，（3-12）將（3-11）、（3-12）代回（3-6）整理后可得，令（下標(biāo)i，j，m輪換）可得（3-13）單元內(nèi)的位移記為單元的結(jié)點位移記為單元內(nèi)的位移函數(shù)可以簡寫成，（3-14）把N稱為形態(tài)矩陣，Ni稱為形態(tài)函數(shù)。選擇單元位移函數(shù)應(yīng)滿足以下條件：1） 反映單元的剛體位移與常量應(yīng)變。2）相鄰單元在公共邊界上的位移連續(xù)，即單元之間不能重疊，也不能脫離。由（3-6）可以將單元位移表示成以下的形式，反映了剛體位移和常應(yīng)變。 單元位移函數(shù)是線性插值函數(shù)，因此單元邊界

11、上各點的位移可以由兩個結(jié)點的位移完全確定。兩個單元的邊界共用兩個結(jié)點，所以邊界上的位移連續(xù)。 形態(tài)函數(shù)Ni具有以下性質(zhì)：1） 在單元結(jié)點上形態(tài)函數(shù)的值為1或為0。2）在單元中的任意一點上，三個形態(tài)函數(shù)之和等于1。用來計算三角形面積時，要注意單元結(jié)點的排列順序，當(dāng)三個結(jié)點i，j，m取逆時針順序時，；當(dāng)三個結(jié)點i，j，m取順時針順序時，。例題：如圖3.7所示等腰三角形單元，求其形態(tài)矩陣N。解：由在公式中輪換下標(biāo)可以計算得，三角形積為形態(tài)函數(shù)為形態(tài)矩陣為三角形面積的計算公式可得，如果把三個結(jié)點按順時針方向排列，即i（a，0），j（0，0），m（0，a）3.3單元載荷移置有限元法的求解對象是單元的組合

12、體，因此作用在彈性體上的外力，需要移置到相應(yīng)的結(jié)點上成為結(jié)點載荷。載荷移置要滿足靜力等效原則。靜力等效是指原載荷與結(jié)點載荷在任意虛位移上做的虛功相等。單元的虛位移可以用結(jié)點的虛位移表示為，（3-15）令結(jié)點載荷為1）集中力的移置如圖3.7所示，在單元內(nèi)任意一點作用集中力圖3.8由虛功相等可得，由于虛位移是任意的，則（3-16）例題1：在均質(zhì)、等厚的三角形單元ijm的任意一點p（xp，yp）上作用有集中載荷。2） 體力的移置令單元所受的均勻分布體力為 由虛功相等可得，（3-17）3） 分布面力的移置設(shè)在單元的邊上分布有面力，同樣可以得到結(jié)點載荷，（3-18）例題2：設(shè)有均質(zhì)、等厚的三角形單元ij

13、m，受到沿y方向的重力載荷qy的作用。求均布體力移置到各結(jié)點的載荷。同理，例題3：在均質(zhì)、等厚的三角形單元ijm的ij邊上作用有沿x方向按三角形分布的載荷，求移置后的結(jié)點載荷。取局部坐標(biāo)s，在i點s=0，在j點s=l，L為ij邊的長度。在ij邊上，以局部坐標(biāo)表示的插值函數(shù)為，載荷為3.4單元剛度矩陣 根據(jù)單元的位移函數(shù)，由幾何方程可以得到單元的應(yīng)變表達(dá)式，（3-19）記為，B矩陣稱為幾何矩陣。B矩陣可以表示為分塊矩陣的形式（3-20） 由物理方程，可以得到單元的應(yīng)力表達(dá)式，（3-21） D稱為彈性矩陣，對于平面應(yīng)力問題， 定義為應(yīng)力矩陣。將應(yīng)力矩陣分塊表示為，（3-22）應(yīng)用虛功原理可以建立單

14、元結(jié)點位移與結(jié)點力的關(guān)系矩陣，單元剛度矩陣。虛功原理：在外力作用下處于平衡狀態(tài)的彈性體，如果發(fā)生了虛位移，則所有外力在虛位移上做的虛功等于內(nèi)應(yīng)力在虛應(yīng)變上做的虛功。單元的結(jié)點力記為單元的虛應(yīng)變?yōu)閱卧耐饬μ摴?，單元的?nèi)力虛功為，由虛功原理可得，（3-23）（3-24）定義為單元剛度矩陣。在3結(jié)點等厚三角形單元中B和D的分量均為常量，則單元剛度矩陣可以表示為，（3-25）單元剛度矩陣表示為分塊矩陣：（3-26）（3-27）3.5單元剛度矩陣的性質(zhì)與物理意義（一）單元剛度矩陣的物理意義假設(shè)單元的結(jié)點位移如下：由，得到結(jié)點力如下：（3-28）表示i結(jié)點在水平方向產(chǎn)生單位位移時，在結(jié)點i的水平方向上

15、需要施加的結(jié)點力。表示i結(jié)點在水平方向產(chǎn)生單位位移時，在結(jié)點i的垂直方向上需要施加的結(jié)點力。選擇不同的單元結(jié)點位移，可以得到單元剛度矩陣中每個元素的物理含義：表示s結(jié)點在水平方向產(chǎn)生單位位移時，在結(jié)點r的水平方向上需要施加的結(jié)點力。表示s結(jié)點在水平方向產(chǎn)生單位位移時，在結(jié)點r的垂直方向上需要施加的結(jié)點力。表示s結(jié)點在垂直方向產(chǎn)生單位位移時，在結(jié)點r的水平方向上需要施加的結(jié)點力。表示s結(jié)點在垂直方向產(chǎn)生單位位移時，在結(jié)點r的垂直方向上需要施加的結(jié)點力。因此單元剛度矩陣中每個元素都可以理解為剛度系數(shù)，即在結(jié)點產(chǎn)生單位位移時需要施加的力。（二）單元剛度矩陣的性質(zhì)1） 對稱性利用分塊矩陣的性質(zhì)證明如下

16、：即2） 奇異性即單元剛度矩陣的行列式為零，。將定單元產(chǎn)生了x方向的剛體移動，此時對應(yīng)的單元結(jié)點力為零。可以得到，在單元剛度矩陣中1，3，5列中對應(yīng)行的系數(shù)相加為零，由行列式的性質(zhì)可知，。同樣如果假定單元產(chǎn)生了y方向上的剛體位移，可以得到，在單元剛度矩陣中2，4，6列中對應(yīng)行的系數(shù)相加為零。3.6整體分析 得到了單元剛度矩陣后，要將單元組成一個整體結(jié)構(gòu)，根據(jù)結(jié)點載荷平衡的原則進(jìn)行分析，即整體分析。整體分析包括以下4個步驟：1） 建立整體剛度矩陣，2） 根據(jù)支承條件修改整體剛度矩陣，3） 解方程組，求出結(jié)點的位移，4） 根據(jù)結(jié)點位移，求出單元的應(yīng)變和應(yīng)力。在這里把結(jié)點位移作為基本未知量求解。如何

17、得到整體剛度矩陣？基本方法是剛度集成法，即整體剛度矩陣是單元剛度矩陣的集成。圖3.9 如圖3.9所示，一個劃分為6個結(jié)點、4個單元的結(jié)構(gòu)。得到了每個單元的單元剛度矩陣后，要集成為整體剛度矩陣。剛度集成法的物理意義 由單元剛度矩陣的物理意義可知，單元剛度矩陣的系數(shù)是由單元結(jié)點產(chǎn)生單位位移時引起的單元結(jié)點力。在如圖3.9所示的結(jié)構(gòu)中，使結(jié)點3產(chǎn)生單位位移時，在單元（1）中的結(jié)點2上引起結(jié)點力。由于結(jié)點2、3同時屬于單元（1）、（3），在單元（2）中的結(jié)點2上同樣也引起結(jié)點力，因此，在整體結(jié)構(gòu)中當(dāng)結(jié)點3產(chǎn)生位移時，結(jié)點2上的結(jié)點力應(yīng)該是單元（1）、（2）在結(jié)點2上的結(jié)點力的疊加。剛體集成法即結(jié)構(gòu)中的

18、結(jié)點力是相關(guān)單元結(jié)點力的疊加，整體剛度矩陣的系數(shù)是相關(guān)單元的單元剛度矩陣系數(shù)的集成。結(jié)點3在整體剛度矩陣的對應(yīng)系數(shù)，應(yīng)該是單元（1）、（3）、（4）中對應(yīng)系數(shù)的集成。3.6.2剛度矩陣集成的規(guī)則1） 將單元剛度矩陣中的每個分塊放到在整體剛度矩陣中的對應(yīng)位置上，得到單元的擴(kuò)大剛度矩陣。單元剛度矩陣系數(shù)取決于單元結(jié)點的局部編號順序，必須知道單元結(jié)點的局部編號與該結(jié)點在整體結(jié)構(gòu)中的總體編號之間的關(guān)系，才能得到單元剛度矩陣中的每個分塊在整體剛度矩陣中的位置。將單元剛度矩陣中的每個分塊按總體編碼順序重新排列后，可以得到單元的擴(kuò)大矩陣。假定單元結(jié)點的局部編號與整體的對應(yīng)關(guān)系如下：單元編號單元結(jié)點局部編號單

19、元結(jié)點整體編號1i31j11m22i52j22m43i53j33m24i34j54m6單元（2）的單元擴(kuò)大矩陣的分塊矩陣形式如下，只列出非零的分塊：局部編號整體編號1234561234562）將全部單元的擴(kuò)大矩陣相加得到整體剛度矩陣。整體剛度矩陣如下所示：整體編號12345612+3+45+63.7約束條件的處理圖3.9所示的結(jié)構(gòu)的約束和載荷情況，如圖3.10所示。結(jié)點1、4上有水平方向的位移約束，結(jié)點4、6上有垂直方向的約束，結(jié)點3上作用有集中力（Px，Py）。圖3.10整體剛度矩陣K求出后，結(jié)構(gòu)上的結(jié)點力可以表示為：根據(jù)力的平衡，結(jié)點上的結(jié)點力與結(jié)點載荷或約束反力平衡。用表示結(jié)點載荷和支桿

20、反力，則可以得到結(jié)點的平衡方程：（3-29）這樣構(gòu)成的結(jié)點平衡方程組，在右邊向量P中存在未知量，因此在求解平衡方程之前，要根據(jù)結(jié)點的位移約束情況修改方程(3-29)。先考慮結(jié)點n有水平方向位移約束，與n結(jié)點水平方向?qū)?yīng)的平衡方程為：（3-30）根據(jù)支承情況，方程（3-30）應(yīng)該換成下面的方程：（3-31）對比公式（3-30）和（3-31），在式（3-29）中應(yīng)該做如下修正：在K矩陣中，第2n-1行的對角線元素改為1，該行中全部非對角線元素改為0；在P中，第2n-1個元素改為0。為了保持K矩陣的對稱性，將第2n-1列的全部非對角元素也改為0。同理，如果結(jié)點n在垂直方向有位移約束，則（3-29）中

21、的第2n個方程修改為， 在K矩陣中，第2n行的對角線元素改為1，該行中全部非對角線元素改為0；在P中，第2n個元素改為0。為了保持K矩陣的對稱性，將第2n列的全部非對角元素也改為0。（3-32）對圖3.9所示結(jié)構(gòu)的整體剛度在修改后可以得到以下的形式，（3-33）如果結(jié)點n處存在一個已知非零的水平方向位移，這時的約束條件為，（3-34）在K矩陣中，第2n-1行的對角線元素乘上一個大數(shù)A，向量P中的對應(yīng)換成，其余的系數(shù)保持不變。方程改為，（3-35）A的取值要足夠大，例如取1010。只有這樣，方程（3-35）才能與方程（3-34）等價。如果結(jié)點n處存在一個已知非零的垂直方向位移，這時的約束條件為，

22、。也可以采用同樣的方法修改整體剛度矩陣。3.8整體剛度矩陣的特點與存儲方法用有限元方法分析復(fù)雜工程問題時，結(jié)點的數(shù)目比較多，整體剛度矩陣的階數(shù)通常也是很高的。那么，是否在進(jìn)行計算時要保存整體剛度矩陣的全部元素？能否根據(jù)整體剛度矩陣的特點提高計算效率？整體剛度矩陣具有以下幾個顯著的特點：對稱性，稀疏性，非零系數(shù)帶形分布。1） 對稱性由單元剛度矩陣的對稱性和整體剛度矩陣的集成規(guī)則，可知整體剛度矩陣必為對稱矩陣。利用對稱性，只保存整體矩陣上三角部分的系數(shù)即可。2） 稀疏性單元剛度矩陣的多數(shù)元素為零，非零元素的個數(shù)只占較小的部分。如圖3.11所示的結(jié)構(gòu)，結(jié)點2只和通過單元聯(lián)接的1、3、4、5結(jié)點相關(guān)，

23、結(jié)點5只和通過單元聯(lián)接的2、3、4、6、8、9結(jié)點相關(guān)。由單元剛度矩陣的物理意義和整體剛度矩陣的形成方式可知，相關(guān)結(jié)點2、3、4、6、8、9及結(jié)點5本身產(chǎn)生位移時，才使結(jié)點5產(chǎn)生結(jié)點力，其余結(jié)點產(chǎn)生位移時不在該結(jié)點處引起結(jié)點力。在用分塊形式表示的整體矩陣中，與相關(guān)結(jié)點對應(yīng)的分塊矩陣具有非零的元素，其它位置上的分塊矩陣的元素為零，如圖3.12所示。圖3.11圖3.12整體剛度矩陣的分塊矩陣示意3）非零元素帶形分布在圖3.12中，明顯可以看出，整體剛度矩陣的非零元素分布在以對角線為中心的帶形區(qū)域內(nèi)，這種矩陣稱為帶形矩陣。在包括對角線元素的半個帶形區(qū)域內(nèi)，每行具有的元素個數(shù)叫做半帶寬，用hbd表示。

24、圖3.11所示結(jié)構(gòu)的相鄰結(jié)點編碼的最大差值為4，所以半帶寬為10。二維等帶寬存儲設(shè)整體剛度矩陣K為一個的矩陣，最大半帶寬為d。利用帶形矩陣的特點和對稱性，只需要保存以d為固定帶寬的上半帶的元素，稱為二維等帶寬存儲。進(jìn)行存儲時，把整體剛度矩陣K每行中的上半帶元素取出，保存在另一個矩陣K*的對應(yīng)行中，得到一個矩陣K*。把元素在K矩陣中的行、列編碼記為r、s，在矩陣K*中的行、列編碼記為r*、s*，對應(yīng)關(guān)系如下：r*=rs*=s-r+1 圖3.13(a)圖3.13(b)如圖3.13(a)所示的最大半帶寬為d的整體剛度矩陣K，采用二維等帶寬存儲后得到如圖3.13(b)所示的矩陣K*。用新的方法存儲后，

25、K矩陣中的對角線元素保存在新矩陣中的第1列中，K矩陣中的r行元素仍然保存在新矩陣的r行中，K矩陣中的s列元素則按照新的列編碼保存在新矩陣的不同列中。采用二維等帶寬存儲，需要保存的元素數(shù)量與K矩陣中的總元素數(shù)量之比為。所存儲的元素數(shù)量取決于最大半帶寬d的值，d的值則由單元結(jié)點的編碼方式?jīng)Q定。雖然在采用二維等帶寬存儲時，仍然會保存一些零元素，但是采用這種方法時元素尋址很方便。圖3.14(a)圖3.14(b)對于同樣的有限元單元網(wǎng)格，按照圖3.14(a)的結(jié)點編碼，最大的半帶寬為14；按照圖3.14(b)的結(jié)點編碼，最大的半帶寬為18；按照圖3.11的結(jié)點編碼，最大的半帶寬為10。3.9線性方程組解

26、法由于有限元分析需要使用較多的單元，線性方程組的階數(shù)很高，有限元求解的效率很大程度上取決于線性方程組的解法。利用矩陣的對稱、稀疏、帶狀分布等特點提高方程求解效率是關(guān)鍵。線性方程組的解法分為兩大類：直接解法，迭代解法。直接解法以高斯消去法為基礎(chǔ)，包括高斯消去法、等帶寬高斯消去法、三角分解法，以及適用于大型方程組求解的分塊算法和波前法等。迭代算法有高斯-賽德爾迭代、超松弛迭代和共軛梯度法等。在方程組的階數(shù)不是特別高時，通常采用直接解法。當(dāng)方程組的階數(shù)過高時，為避免舍入誤差和消元時有效數(shù)損失等對計算精度的影響，可以選擇迭代方法。這里給出了用Fortran語言編寫的等帶寬高斯消去法的代碼，其中NROW

27、為矩陣行的數(shù)目，NHBW為最大半帶寬。SUBROUTINE SOLVERB(NROW,NHBW,STIFF,DISPL)C Band elimination methodC .C : SOLVE EQUATIONS KS*H=(Q) :C :.:C KS is stored with the half band width methodC Input:CNROW - the quantity of rowsCNHBW -half band widthCSTIFFNROW,NHBWC - coefficient matrix stored with the half width methodCDISPLNROW - the right hand vectorC Output:CDISPL - results of variablesC Warning:Carray STIFF is changed.C*ADD:DPR* IMPLICIT DOUBLE PRECISION ( A-H,O-Z )C