1、第1章 緒 論1.1 課程內容(1) 研究內容本課程主要研究工程結構計算機分析（數值分析）的常用方法有限單元法、加權殘數（余量）法和邊界單元法的基本概念、基本原理及其應用。(2) 參考書籍課程的主要參考書籍如下：唐錦春，孫炳楠，郭鼎康，計算結構力學，浙江大學出版社，1989丁皓江, 謝貽權, 何福保，彈性和塑性力學中的有限單元法，機械工業出版社，1989王勖成，有限單元法，清華大學出版社，2003王勖成，邵敏，有限單元法基本原理與數值方法，第二版，清華大學出版社，1997徐次達，固體力學加權殘數法，同濟大學出版社，1987孫炳楠，項玉寅，張永元，工程中邊界單元法及其應用，浙江大學出版社，199

2、1Bath, K. J. Finite Element Procedures, Prentice-Hall, Inc., 1996.Zienkiewicz, O. C., The Finite Element Method, 5th Edition, McGraw Hill, 2001.Brebbia, C.A., The Boundary Element Method for Engineers, Pentech Press, London, 1978.Chandrupatla, T. R., Belegundu, A.D. Introduction to Finite Elements i

3、n Engineering, Prentice-Hall, Inc., 2002.1.2 結構分析方法概述一個工程技術問題總可由一組基本方程（通常是微分方程）加一組邊界條件描述，即由下式給出：基本方程：L(u)-p=0， ÎV（域內）邊界條件：B(u)-g=0， ÎS（邊界）式中L、B為算子，p、g為已知函數。工程技術問題的常用分析方法有：(1) 解析方法只適用于少數簡單問題，即形狀規則且外部作用（如外荷載）簡單的結構分析問題。(2) 數值方法數值方法可分為區域型方法和邊界型方法。常用的區域型方法包括有限差分法、加權殘數法、里茲(Ritz)法（變分法）和有限單元法等，其中有

4、限差分法是直接對基本微分方程進行離散，再對離散后的代數方程進行求解；后幾種方法則是先建立基本方程（一般是微分方程）的等效積分表達式，再進行離散求解。邊界型方法中最典型的是邊界單元法。它是先將基本微分方程變換為等效的邊界積分方程，再在邊界上對其進行離散求解。例如，圖1.1給出了一個受復雜橫向荷載（分布荷載、集中力、集中力偶等）作用的兩端固定變截面梁。為求梁的撓度和內力，可列出梁的基本方程和邊界條件如下：lq(x)yxh(x), EI(x)PmABbh(x)圖1.1 變截面單跨梁受橫向荷載作用基本方程：L(u)-p=0，ÎV（域內）EI(x)y = -M(x), 0£x

5、3;l.邊界條件：B(u)-g=0， ÎS（邊界）(y)x=0或x=l=0，(y )x=0或x=l=0以下分別就采用加權殘數法、里茲法（位移變分法）和有限單元法的基本原理進行討論。(1) 加權殘數法為求近似解，設試探函數代入基本方程和邊界條件，得殘值：RL=L(u)-p（域內），RB=B(u)-g（邊界）迫使殘值在某種平均意義（加權積分）上等于零，則有由此可得到關于待定系數ai的代數方程組，解方程可求得待定系數及解答的近似表達式，其中的試函數可以選擇多項式、三角函數、樣條函數等。(2) 里茲法（位移變分法）里茲法的理論依據是最小勢能原理。該原理可表述為：給定外力作用下，滿足幾何條件的

6、各種可能位移中，真實的位移使總勢能取極值，據此有d(U+UR)=0假設滿足位移邊界條件的位移函數為：將其代入方程得到關于待定系數Ai的代數方程，解方程可得Ai。里茲法需要在整個計算區域上假設近似函數，很難適應形狀（邊界）較復雜或解答較難預測的問題。(3) 有限單元法有限單元法的理論依據是最小勢能原理或其他形式的變分原理。該方法與里茲法的主要區別是不在整體計算區域上假設近似函數，而是先將連續的求解區域離散為一個由有限個單元組成并按一定方式相互連接的單元集合體，再以各單元連接結點處的場量（如位移量）作為基本未知量，在各單元內假設近似函數（通過結點未知量插值得到），從而將一個無限自由度問題簡化為有限

7、自由度問題。u分段假設試函數x圖1.2 一維試函數的分段假設例如圖1.2中的曲線是某個一維問題的目標函數曲線，若采用里茲法對整個區段假設一個近似的試函數，顯然比較困難。但如果現對整個區域進行分段（如圖中短線為分段線），再對各個區段假設試函數，則要簡單和準確得多，如可將各區段均假設為二次函數。喲次可見，有限單元法可視為一種分片（或分塊、分段）形式的變分法。雖然有限單元法的理論依據和里茲法是一致的，但采用了分片（或分塊、分段）假設試函數的處理方法以后，使得該方法的具體實施變得簡便易行，具有了優越的可操作性和更為明確的物理意義，也使得該方法具有了其他方法（如里茲法）所不具備的優點：1) 概念簡單、明

8、確，易為工程人員接受；也可建立嚴格的數學分析和證明；2) 適用性十分廣泛，適應于各類復雜邊界和不同外部作用的問題；3) 求解過程程序化，易于編程和計算機實現。1.3 課時安排課程的總體課時安排如下：有限單元法部分包括概論、進展；平面三角形、矩形、等參元；桿元、板元等，共約20個課時；加權殘數法部分包括基本原理、方法分類，以及伽遼金（Galerkin）方法、最小二乘法的應用，共需約46個課時；邊界單元法主要包括基本原理（以二維勢問題為例）；梁彎曲和板彎曲問題，共需約46個課時。思考題1.1 區域型分析法和邊界型分析法在對問題的基本方程和邊界條件的處理上有何不同和相同點？試分別舉例說明。1.2 里

9、茲法和有限單元法的理論依據、基本未知量的選取、試函數的假設等方面有何異同點？1.3 與里茲法相比，有限單元法在解決復雜問題上的適應性更為廣泛，你認為主要的原因在于那些方面？第2章 有限單元法2.1 概述 發展概況有限單元法的發展概況：1943年R. Courant嘗試應用三角形區域上定義的分片連續函數和最小勢能原理解決St. Venant扭轉問題，是較早的有限元思想的體現：R. Courant, Variational Methods for the Solution of Problems of Equilibrium and Vibrations, Bulletin of the Amer

10、ican Mathematical Society, 49: 1-23, 19431956年M.J. Turner，R.W. Clough等將剛架矩陣位移法推廣到彈性力學平面問題，開始了有限元的第一個成功嘗試和應用；用直接剛度法建立單元剛度特性：M.J. Turner, R.W. Clough, H.C. Martin and L.T. Topp, Stiffness and deflection analysis of complex structures, J. Aeronaut. Sci., 25: 805-823, 19561960年Clough第一次提出“有限單元法（FEM）”的名稱

11、，沿用至今。Zienkiewicz等編寫第一本有限元方面專著：O.C. Zienkiewicz and Y.K. Cheung, The Finite Element Method in Continuum and Structural Mechanics, McGraw-Hill, New York, 19651963-1964年發現該方法是基于變分原理的里茲法的另一種形式，確立其理論基礎。我國馮康在同一時期獨立提出并證明了該方法：Melosh證明了位移法就是基于最小勢能原理的Rayleigh-Ritz法馮康，基于變分原理的差分格式，應用數學和計算數學，1965，2(4): 238-2621

12、960至今：實際工程應用：平面Þ空間Þ板殼；靜力Þ動力、波動Þ穩定；彈性Þ塑性Þ粘彈性、復合材料；固體Þ流體、傳熱等連續介質力學；計算分析Þ優化設計、與CAD技術結合。E.L. Wilson：編寫了第一個公開的有限元軟件SAP；通用有限元軟件：SAP、ADINA、NASTRAN、ANSYS、ABAQUS、MIDAS等從半個多世紀以來有限單元法的萌芽、理論依據的證明和充實及其逐步的廣泛應用可以看到，它的發展和計算機軟硬件的發展基本上是同步的。如果沒有計算機的強大軟硬件支撐，有限單元法只有其微不足道的一點理論上的意義，

13、而沒有更為重要的實際應用的意義。 有限單元法概念(1) 離散化離散化的過程是將連續體劃分為有限數目、有限大小的單元的集合體。單元與單元之間只在指定點（即結點）連接，其他位置則一般保持連續即可。單元可以具有不同的形狀，即單元外形可以不同；單元與單元之間可以有不同的連接方式，即單元的結點數目、位置可以不同。連續體rg圖2.1 連續體離散為單元集合體示例(2) 單元分析對典型單元假設位移模式（由各結點位移插值），再分析單元的力學特性，建立單元的結點力與結點位移之間的關系，即單元剛度方程：Fe=kDe）并將各類荷載變換為作用在結點上的等效結點荷載。(3) 整體分析將各單元剛度方程集成整體結構的整體剛度

14、方程：F=KD根據結點的平衡條件，得最終的有限元方程：KD=R求解該方程可得到未知的結點位移。(4) 再次單元分析求出各單元的應變和應力。2.2 彈性力學平面問題的矩陣描述兩類平面問題彈性力學的平面問題可分為平面應力問題和平面應變問題兩類。實際上，所有的彈性力學問題都是空間問題。所謂平面問題，并不是說這個問題所分析的對象本身（如形狀、荷載分布）是平面的，而是指該問題的形狀、外部作用以及問題的解答（即由此產生的效應，如位移、應力等）只在平面內有變化，而沿著平面外就保持不變了。因此可以肯定地說，所謂的平面問題就是一個特殊的空間問題。那么，是不是一個問題的形狀和外部作用（即已知的位移和應力邊界條件）

15、只在平面內發生變化，而沿著平面外保持不變了，這個問題就是平面問題呢？不是的，還必須附加其他條件，這一結論才能成立。這個附加條件就是該問題沿平面外的尺寸與平面內尺寸相比要么非常小（如無限短），要么非常大（如無限長）。如果符合前者條件，則彈性體內只存在平面內的應力，而平面外的應力均為零，故這類問題稱為平面應力問題；如果符合后者條件，則彈性體內只存在平面內的位移或平面內的應變，而平面外的位移及應變均等于零，故這類問題稱為平面應變問題。zxt/2Oyt/2y體積力表面力圖2.2 平面應力問題示例xOyz圖2.3 平面應變問題示例2.2.2 基本量和基本方程的矩陣表示(1) 基本量應力：s=sx sy

16、txyT，應變：e=ex ey gxyT，位移：f=u vT；體積力：p=px pyT，表面力：q=qx qyT。(2) 基本方程幾何方程：物理方程：s=D e，式中D：彈性矩陣。對平面應力問題，平衡方程的弱形式能量原理1) 虛功原理：彈性體在外力作用下處于平衡狀態的充分和必要條件是，外力在滿足變形協調推進的虛位移上所做的虛功等于其內力在相應的虛應變上所做的虛功，即2) 最小勢能原理：滿足幾何條件的各種可能位移中，真實位移使P取駐值，即dP=d (U+UR)=0(3) 平面問題的常用單元：三結點三角形單元 六結點三角形單元 四結點矩形單元四結點任意四邊形單元 八結點曲邊四邊形單元圖2.4 平面

17、問題的常用單元單元的加密方法可分為p型加密和h型加密。前者是保持單元的大小和形狀不變，而提高單元的插值函數的階數；后者是指采用同一類單元，但加密單元的網格，即減小單元的尺寸，增加離散單元的數目。2.3 三結點三角形單元分析 離散將一平面結構離散為ne個單元，設結點綜述為n個，則結構的基本未知量取為這n個結點的位移，即D=u1 v1 u2 v2 un vnTyui (Ui)ixvi (Vi)uj (Uj) jvj (Vj)um (Um)mvm (Vm)圖2.5 三結點三角形單元 單元位移模式單元結點力向量：Fe=Ui Vi Uj Vj Um Vm T單元結點位移向量：De=ui vi uj vj

18、 um vmT(1) 位移模式的建立所謂位移函數或稱位移模式，是指利用單元的結點位移將單元內任一點的位移用坐標的函數表示出來。對于三結點三角形單元，假設其位移模式是坐標的線性函數，則有u=a1+a2 x+a3 yv=a4+a5 x+a6y系數a1 a6由6個結點位移分量（ui、vi、uj、vj、um、vm）確定。將位移模式寫成結點位移的顯式：u= Niui+ Nj uj+ Nmumv= Nivi+ Nj vj+ NmvmNi、Nj、Nm稱為形狀函數（形函數）或插值函數，其中ai、bi、ci是分母行列式第1行各元素的代數余子式，展開后可表示為ai=xjym- xmyj，bi=yj-ym，ci=

19、-xj+xm（i, j, m輪換） jyimxP j1imNi(x,y)圖2.6 形函數的物理意義(2) 形函數性質(a) (Ni)i=1，(Ni)j=0，(Ni)m=0；(b) 單元內任一點：Ni+Nj+Nm=1。(3) 位移模式的矩陣表示其中N為形函數矩形，且。 單元應變矩陣和應力矩陣將位移模式f=NDe代入幾何方程，得單元應變為：e=BDeB稱為單元應變矩陣，B= Bi Bj Bm，而各子塊為。可見B為一常量矩陣，故三結點三角形單元為一常應變單元。將e=BDe代入物理方程，可得s=SDe=DBDeS稱為單元應力矩陣，S= S1 S2 S3，對平面應力問題。可見S也是一常量矩陣，故三結點三

20、角形單元為一常應力單元。有限元方程的建立將連續體近似看作為由一系列只在結點相連的單元組裝而成的集合體，并且對各個單元均建立了其位移模式，那么該集合體的最終位移法方程，即最終的有限元方程可由變分原理（此時為最小勢能原理）建立。位移模式確定后，離散結構的可能位移就是由不同結點位移決定的分片函數。各可能位移函數中，真實的位移函數使離散結構的總勢能取駐值（處于穩定平衡狀態時為最小值），即dP=d (U+UR)=0令：De=GD這里G是一6´2n的位置矩陣，表示單元結點位移De在整體結點位移向量D中的位置；ke為單元剛度矩陣；Re單元等效結點荷載向量。于是有令：這里K為整體剛度矩陣；R為整體等

21、效結點荷載向量，則有：由勢能駐值原理dP=0可得或 單元剛度矩陣(1) 單剛列式其中，(2) 單剛性質單元剛度元素kpq的物理意義為：第q個結點位移分量為單位位移（其它結點位移等于0），所引起的第p個結點力分量。例如要獲得元素k26的直觀解釋，可先將單元的所有6個結點位移約束住，即在每個結點處沿兩個方向同時施加鏈桿約束，如圖2.7所示。ymxk26ij1圖2.7 單元剛度元素k26的物理意義單元剛度矩陣k的性質： 1) 對稱性： kpq = kqp； 2) 奇異性：因具有剛體位移； 3) 每行（列）元素之和為零；例如，第6列元素的意義：當第六個結點位移等于1（D6=1），其他結點位移等于0時，

22、在各結點位移方向施加的結點力的大小，即圖2.7中6個附加鏈桿約束上的約束力。因單元在這些結點力作用下處于平衡，故有：k16+ k36+ k56=0；k26+ k46+ k66=0。4) 主對角元素恒大于零：kqp>0； 5) k取決于單元的形狀、方位和彈性常數，與所在位置（即平移或np 轉動）無關。(3) 推導單剛的另一種方法由單元的平衡條件導出物理方程幾何方程位移模式虛功方程f =Ndee=Bdes=Sde ，S= DBFe=kd e，k= BT D BtA結點位移 位移 應變 應力 結點力 de f e s Fe圖2.8 單元分析中各物理量的聯系圖假設單元發生了虛位移，則結點力在結點

23、虛位移上做的虛功=單元應力在虛應變上的虛功（虛變形能）：令：又因dDe是任意的，故有Fe=keDe 單元等效結點荷載在將連續體用近似的單元集合體代替以后，單元集合體中的每個結點可認為一方面受到外荷載的作用，另一方面受到連接該結點的各個單元的對其產生的約束力（即結點力）的作用。各結點在這兩組力的作用下總是處于平衡的。但是，如果外荷載不是直接作用在結點上，而是作用在單元內部，那么怎么建立結點的平衡關系呢？一種有效方法是先建立單元的等效結點荷載，再對結點建立平衡關系。yYiixXiYj jXjYmmXmPxPy建立單元的等效結點荷載，或稱將非結點荷載需等效移置到結點上，該等效移置可利用虛功原理進行。

24、圖2.9 單元的等效結點荷載分量 整體剛度矩陣(1) 建立有限元方程的方法：1) 由最小勢能原理建立：K=å(Ge)Tk e Ge2) 由結點平衡條件建立(2) 整體剛度矩陣K的集成方法：對號入座。K=å(Ge)Tk e Ge，其中(Ge)T決定kpqe在K中的行，Ge決定在K中的列位置。(3) 整體剛度矩陣的性質K元素的物理意義 Kpq：結構第q個結點位移分量為單位位移（Dq=1，其它結點位移=0），在第p個結點位移方向所施加的結點力的大小。 如 K25為結構第5個結點位移分量為單位位移，即D5=1，其它結點位移均為零時，在第2個結點位移方向所需施加的結點力的大小。K的性

25、質： 1) 對稱性： Kpq = Kqp；主對角元恒大于零； 2) 稀疏矩陣，且一般呈帶狀分布；例如，平面問題最大半帶寬=2´(單元結點號之差最大值+1)。 3) 引入位移邊界條件前為奇異矩陣，引入后為正定矩陣。 位移邊界條件的引入以上集成的總體整體剛度矩陣并不包含位移邊界條件的信息，這種剛度矩陣常稱為原始剛度矩陣，它是一個奇異矩陣。要使得最終的有限元方程可解，必須引入唯一邊界條件，排除結構的剛體位移。引入位移邊界條件的常用方法有以下幾種：(1) 直接引入法（矩陣縮小法）將KD=R改寫為如下子塊形式：其中Da、Db分別為未知結點位移向量和已知結點位移向量，Kaa、Kab、Kba、Kb

26、b、Ra、Rb分別為與之對應的剛度矩陣和荷載向量子塊。將上述方程的上半部分（對應已知位移）展開，得Kaa Da=Ra-Kab Db由該方程可解出未知結點位移。該方法將改變原始剛度矩陣的階數及元素的順序，不利于程序實現，因此在計算機編程中一般很少采用。(2) 對角元素改1法（零位移邊界）設結構的總自由度數（即結點位移總數）為N，若第i個結點位移為零（即Di=0），則將K中對角元素Kii改為1，而第i行和第i列的其他元素改為0，荷載向量R中的第i個元素也改為0。其實質是將原有限方程中的第i個方程用方程Di=0代替，而其他方程中與Di對應的系數也改為0，表明Di對其他方程沒有影響，同時保證了修改后的

27、剛度矩陣仍具有對稱性。(3) 乘大數法若第i個結點位移已知，即，則將整體剛度矩陣K中的對角元素Kii改為a Kii，其中a為一大數（如1020），而荷載向量R中的第i個元素Ri改為，原方程成為以下形式其實質是將原有限方程中的第i個方程用的以下近似方程代替：將上述方程各項同除以大數a，除第i項及方程右端項外，其他各項均趨于0，故等價于。 位移模式與解答的收斂準則(1) 有限元解答的收斂準則為使有限元解答能夠收斂于精確解，單元位移模式應滿足以下條件：1) 位移模式必須包含單元的剛體位移；對彈性力學平面問題，其剛體位移表達式為：u=u0-wy，v=v0+wx因此，平面問題的位移模式必要包含上述剛體位

28、移表達式中的各項，才能保證最終解答的收斂性。2) 位移模式必須包含單元的常量應變；條件(1)、(2)合并起來可稱為完備性要求。對平面問題來說，就是要求位移模式必須包含常數項和完整的一次項。完備性條件是解答收斂的必要條件。3) 位移模式應盡可能保證位移的連續性。該條件實際上就是協調性條件，但一般情況下并不是一個必要性條件。如果位移模式同時滿足上述完備性和協調性條件，那么就組成了解答收斂的充分條件。對一般的彈性力學平面或空間問題，只需要求單元內部以及相鄰單元的公共邊界上的位移本身連續，即位移模式具有C0連續性。對有些問題，可以放松對協調性的要求，只要通過分片試驗，那么也能保證解答的收斂性。三結點三

29、角形單元的位移模式既滿足完備性，又滿足協調性的要求（在單元邊界上呈線性分布，可由兩個結點位移唯一確定），是一種協調單元。證明如下：單元內：單值連續；yx(1)ijmp(2)相鄰單元之間：uij(1)= uij(2)？ vij(1)= vij(2)？ij邊的方程：y=ax+b，則uij=a1+ a2x+a3(ax+b)= cx+duij(1)、uij(2)均為坐標的線性函數，故可由i、j兩點的結點位移唯一確定。圖2.10 兩相鄰單元(2) 多項式位移模式的一般選擇規則位移模式應與單元局部坐標的選取無關，即滿足所謂的幾何各向同性。對于平面問題，位移模式中的x和y的各階項應保持對稱，即有了xmyn項

30、，則應同時具有xnym項。為保證位移模式關于x和y坐標的對稱性，通常從以下的Pascal三角形中選取多項式位移模式的各項：1x yx2 xy y2 x3 x2y x y2 y3 x4 x3y x2y2 xy3 y4圖2.11 Pascal三角形根據最小勢能原理建立的有限元，是以結點位移作為基本未知量，這種有限單元稱為位移元。由位移元得到的近似解答總體上是精確解的一個下限。2.4 高精度的三角形單元、矩形單元 高精度三角形單元(1) 六結點三角形單元（二次單元T6）位移模式取坐標的完整二次式：u=a1+a2 x+a3 y+a4 x2+a5 xy+a6 y2v=b1+b 2 x+b 3 y+b 4

31、 x2+b5 xy+b6 y2該位移模式包含了常數項和完整的一次項，滿足完備性要求；在單元的邊界上位移呈二次拋物線分布，可由三個結點位移唯一確定，故又滿足協調性的要求，是一種協調單元。(2) 十結點三角形單元（三次單元T10）位移模式取坐標的完整三次式：u=a1+a2 x+a3 y+a4 x2+a5 xy+a6 y2+a7x3+a8x2y+a9xy2+a10y3v=b1+b 2 x+b 3 y+b 4 x2+b5 xy+b6 y2+b7x3+b8x2y+b9xy2+b10y3該位移模式包含了坐標的完整一次式（常數項和純一次項），滿足完備性要求；在單元的邊界上位移呈三次曲線分布，可由4個結點位移

32、唯一確定，故又滿足協調性的要求，是一種協調單元。ijmijm3122.12 高精度三角形單元(3) 面積坐標表示的三角形單元形函數在推導三角形單元的列式時，直接利用整體坐標系（為直角坐標）下的位移模式將使得運算十分繁瑣和復雜。如果采用三角形單元內的一種局部坐標面積坐標作為自然坐標，則可以使列式推導大為簡化。定義單元內任一點P的無量綱面積坐標（Li, Lj, Lm）為i(1,0,0)Pj(0,1,0)m(0,0,1)AiAjAmLi= Ai/ A （i, j, m）各種單元的形函數：3結點三角形單元（線性單元T3）：Ni=Li （i, j, m）6結點三角形單元（二次單元T6）：Ni=(2Li-

33、1) Li （i, j, m）N1=4Lj Lm, （1, 2, 3；i, j, m）2.13 三個結點單元示意圖面積坐標的特點：1) 三角形三個角點的面積坐標：i(1,0,0)、j(0,1,0)、m(0,0,1)三條邊用面積坐標表示的方程為：jm邊 Li=0mi邊 Lj=0ij邊 Lm=02) 三個面積坐標不獨立，其相互關系為Li+ Lj+ Lm=1面積坐標與直角坐標之間的關系：或 四結點矩形單元（R4單元）(1) 采用雙線性的位移模式：u=a1+a2 x+a3 y+a4xy1yx234hxov=b1+b 2 x+b 3 y+b 4xy可保證位移模式的完備性和協調性。若寫成形函數形式，則為2

34、.14 四結點矩形單元這里的Ni（i=1, 2, 3, 4）可以先求出8個待定系數在獲得，也可以根據形函數的性質直接構造，例如對N1，可設該表達式滿足在結點2、3、4取值均為零的性質；再令Ni(1)=1，可得待定系數這里設矩形單元的邊長各為2a、2b。(2) 局部坐標下的形函數設矩形長和寬各為2a和2b，在矩形形心為原點建立局部坐標（自然坐標）系xoh，它與直角坐標的變換式：x4h123x= -1x=1h= -1h=1o2.15 局部坐標系下的矩形單元則四個角點的局部坐標分別為1(-1, -1)，2(1, -1)，3(1,1)，4(-1,1)。位移模式為：Ni（i=1, 2, 3, 4）可根據

35、形函數的性質直接構造出來： 或(3) 單元應變矩陣和應力矩陣單元應變為：e=BDe其中B=B1 B2 B3 B4，。單元應力為：s=SDe=DBDe其中S= S1 S2 S3 S4，Si=D Bi2.5 C0連續型單元形函數的構造在有限單元法中，根據結點布置方式的不同可將矩形單元分為兩類，一類稱為Lagrange矩形單元，另一類稱為Serendipity矩形單元。前者在單元縱橫網格線的交點上均布置結點，而后者僅在單元的邊界線上布置結點，如圖2.17所示。(a) Lagrange矩形單元 (b) Serendipity矩形單元2.16 兩類矩形單元 Lagrange矩形單元(1) 一維Lagra

36、nge插值多項式過n個結點（坐標分別為x1, x2, , xn）的一維Lagrange插值多項式為xx1x2x3xnxi12.17 一維Lagrange插值多項式l1(x)例如n=2，則有：。若令x1=0，x2= l，則。用該多項式作為形函數，可滿足形函數的性質。(2) Lagrange矩形單元的形函數在矩形的各個網格交點上均布置結點，如水平方向r+1個，豎向p+1個。于是第I列J行結點i的相應形函數：Ni= NIJ= lI(r)(x) lJ(p)(h )其中，lI(r)(x)、lJ(p)(h )均為一維Lagrange插值多項式，Ni在結點i為1，在其他結點處為零；單元邊界上的結點數=形函數

37、階數+1，能夠保證邊界位移的唯一性和協調性。這類單元包含較多的內部結點，增加了單元的自由度，而實踐證明這些自由度的增加通常并不能有效提高單元的精度。對n´n結點的Lagrange單元，其Pascal三角形中包含的項數如圖2.19所示。(I,J)(r, p)(0, 0)11lI (x)lJ (h)1xnynxnynxx2x3yy2y3圖2.18 Lagrange矩形單元及插值模式 圖2.19 Pascal三角形包含項數. Serendipity矩形單元(1) 結點布置特點這類單元只在其邊界上布置結點，但不同邊界上可布置有不同數目的結點。(2) 形函數的構造方法如果只有四個角點上布置結點

38、（R4單元），則如果增加結點5使之成為5結點矩形單元（或稱R5單元），則滿足(N5) j=d5j （j=1, 2, 3, 4, 5）。原不滿足，需對其進行修正才能滿足：(N1)5=0，(N2)5=0，且保證N1，N2在除5結點外的其他4個結點取值不變。x4h12351011/21N5圖2.20 五結點矩形單元 圖2.21 形函數N1的修正對于8結點Serendipity矩形單元（R8單元），利用上述方法，可得到前4個結點對應的形函數分別為或寫成統一表達式：后4個形函數分別為對于p次單元的邊界內結點，其相應的形函數為x（或h）的p次和h（或x）的一次Lagrange多項式的乘積；而角結點的形函數

39、為四結點單元的相應函數與各內結點形函數乘以一分數的差。利用該方法可以很方便地構造一些過渡單元的形函數。2.6 平面等參單元直接用形狀規則的單元對復雜形狀的結構進行離散繪遇到邊界難于處理的問題，因此應設法采用適當的方法對規則單元進行坐標變換，使之轉化為斜邊或曲邊的單元。有限元中最常用的變換方法是等參變換，即坐標變換采用與單元位移插值一致的形函數。通過等參變換的單元稱為等參單元。 四結點四邊形等參元對于4結點任意四邊形單元，若采用雙線性的位移模式，則位移在單元斜邊界上呈二次拋物線分布，由兩個結點位移不能唯一確定。為此在單元內設法建立一個局部坐標（自然坐標），使得單元各邊界的局部坐標分別為一常量（或

40、±1），這樣，如果位移模式采用局部坐標的雙線性函數，則能夠滿足公共邊界位移的唯一性。如果單元邊界的局部坐標分別為±1，則在局部坐標下，單元的形態就是圖示的4結點正方形單元，該單元稱為基本單元或母單元，其形函數為。位移模式為：1yx234xhx4h123x= -1x=1h= -1h=1o(a) 基本單元 (b) 實際單元2.22 四結點等參數單元如果同時利用上述形函數作為局部（自然）坐標(x, h)向整體（直角）坐標(x, y)的變換函數，則可以建立兩種坐標之間的映射關系如下：根據形函數的性質：(Ni) j=dij，SNi=1，容易驗證上述變換式在四個結點處給出了節點的整體坐

41、標，而在單元的四邊上，其中一個局部坐標為±1，另一局部坐標按線性變化，因而正好給出了整體坐標下四邊的直線方程（由兩個結點整體坐標唯一確定）。例如對12邊，因此上述變換式正確地反映了局部（自然）與整體（直角）坐標之間的映射關系。經變換后的實際單元稱為子單元。子單元的位移模式仍為：利用等參變換可以構造8結點、12結點、20結點等更高次的四邊形曲邊等參單元（參見圖2.23）。hxo15247386yx15247386圖2.23 8結點曲邊四邊形等參單元 局部與整體坐標的微分和積分變換式根據復合函數的求導法則，寫成矩陣形式，式中J稱為雅可比（Jacobi）矩陣，對于具有m個結點的平面等參單元

42、，若反過來用局部坐標表示整體坐標，則可對上式作反變換，式中，Jij*是J各元素Jji的代數余子式；½J½是J的行列式，稱為雅可比（Jacobi）行列式。面積微分的變換：dA=dxdy=½J½dxdh。 單元剛度矩陣和單元等效結點荷載向量對于一具有m個結點的平面等參單元，其應變向量可寫為e=BDe其中，單元應變矩陣：B=B1 B2 Bm。單元應力向量：s=SDe=DBDe單元剛度矩陣：單元等效結點荷載列陣：式中Rpe、Rqe分別為體積力和表面力引起的等效結點荷載，且對其中一個局部坐標（x或h）為常數的邊界，線積分ds為或上述積分式通常采用Newton-Co

43、tes或高斯（Gauss）數值積分求得。可見，等參單元的剛度矩陣、等效結點荷載列陣等的計算都在規則的母單元中進行，因此各種形式的積分運算都可以采用標準的數值積分方法進行，使得不同工程問題的有限元分析能夠采用統一的通用化程序實現。 等參變換的應用條件等參單元的應用條件是母單元與實際單元之間存在一一對應關系。具體到局部與整體坐標的變換式，要求變換矩陣即雅可比J非奇異，雅可比（Jacobi）行列式½J½¹0。為確保坐標變換一一對應，在單元劃分時應避免：(1) 使任意兩個結點退化為一個結點而使½dx½=0或½dh½=0；(2) 還應

44、避免因單元過于歪斜而導致dx與dh發生共線。xh123x= -1x=1h= -1h=1o4xh123,4x= -1x=1h= -1h=1o(a) 結點退化 (b) 單元邊界接近共線2.24 畸形等參數單元 例題分析例2.1 圖示懸臂平面結構，長2m，高1m，試用圖示單元進行分析。12346(1)(2)(3)5(4)xyP jmi(1)(3) jmi(2)(4)(a) 基本單元 (b) 實際單元2.25 例2.1圖（三角形單元）1. 三結點三角形單元分析(1) 離散化：劃分為4個單元，共6個結點。(2) 依據圖示的單元局部編號規則，4個單元的剛度矩陣均相同，為為簡便起見，設泊松比n=0，于是有(

45、3) 根據局部與整體結點編號的對應關系集成整體剛度矩陣：(4) 引進位移邊界條件2. 四結點四邊角形等參單元分析12346(1)(2)5xyPx4h1232.26 例2.1圖（四邊形單元）計算步驟：(1) 離散化：劃分為2個單元，共6個結點。單元上下短邊0.75m，長邊1.25m。(2) 求，J，|J|，J-1；(3) 求B在各積分點的數值Big；(4) 利用高斯積分計算并形成k；(5) 集成K、P；(6) 引進位移邊界條件。12P例2.2 牛腿受豎向集中力作用，且結點1、2發生已知的水平位移，求應力。2.27 例2.2圖2.7 有限元程序實現有限元程序前處理程序：生成網格及數據文件主體分析程

46、序：核心計算分析后處理程序：進行結果處理，生成可直接使用或查看的結果文件、圖形文件。 有限元程序設計的一般步驟1. 算法描述和列式推導；2. 框圖設計；3. 代碼編寫；4. 上機調試、考核；5. 編寫應用說明；6. 修改、補充、完善。程序設計的一般要求：具備較齊全的功能、較強的通用性和可移植性、較好的可擴充性、良好的可讀性、足夠的可靠性、良好的自適應性。 輸入數據及分類1. 控制數據：結點總數、單元總數、約束總數、荷載總數、問題類型數等2. 幾何數據：結點坐標、單元信息（各單元的結點編號）、約束條件、單元類型數（彈性模量E、泊松比n、厚度t不同為一類）3. 物性數據：彈性模量E、泊松比n、厚度

47、t4. 荷載數據：荷載類型（集中、分布）、位置、方向、大小等 三結點三角形單元分析平面問題的主體程序1. 程序框圖（參見圖2.28）2. 結構化程序設計方法模塊化由1個主程序和若干子程序組成。子程序由通用性，采用可調數組；主程序采用動態數組存儲技術3. 動態數組存儲技術（1）按整型和實型定義兩個大型共享數組，如A(1000000)、M(1000000)；（2）設計動態數組表：將各子程序中的變界（可調）數組按各自實際需要的大小分配一維數組的空間；（3）動態數組覆蓋技術：全部或部分覆蓋。讀入控制數據開始讀入幾何、物性、荷載數據平面應力/應變?E=E0/(1-n02),n=n0/(1-n0)E=E0

48、,n=n0計算單元剛度元素疊加到整體剛度矩陣中計算等效結點荷載引入位移約束條件解線性代數方程得結點位移計算單元應力、反力等輸出結果結束圖2.28 平面問題主體分析程序框圖2.8 平面桿系結構的有限元2.8.1 等截面直梁單元（忽略剪切變形）(1) 基本方程xyp(x)q(x)圖2.29 受任意荷載的等截面直梁幾何關系：寫成矩陣系數，f= u vT，內力-位移關系：平衡關系：邊界條件：(2) 離散將一平面桿件結構離散為ne個單元，n個結點，基本未知量為結點位移D=u1 v1 q1 u2 v2 q2 un vn qnT對應的結點力向量為F=X1 Y1 M1 X2 Y2 M2 Xn Yn MnT(3

49、) 位移模式局部坐標下的單元結點位移向量和單元結點力向量：xyxyviuivjujqiqja圖2.30 等截面直梁單元由梁的平衡方程可知，在結點荷載作用下梁的軸向位移沿梁軸呈線性分布（因只有桿端作用軸向荷載），橫向位移呈三次曲線分布（因彎矩沿桿軸成線性分布），故假設u=a1+a2 xv=a3+a4 x+a5 x2+a6 x3將結點i, j的位移代入可求出6個待定系數a1 a6。將單元位移寫成結點位移的顯式，有u=N1ui+N4 ujv= N2vi+ N3qi+ N5vj+ N6qj位移模式的矩陣表示u、v獨立插值，但q不獨立插值，故要求C1連續：不僅u、v本身連續，v的一階導數也要連續。(4)

50、 單元應變矩陣和應力矩陣將位移模式f=NDe代入幾何方程，得單元應變為：將e=BDe代入物理方程，可得S=DBDe,其中(5) 單元剛度矩陣局部坐標下的單元剛度矩陣：單元剛度方程為：Fe =kDek中的每一元素krs稱為單元剛度系數，其物理意義是：第s個桿端位移為單位位移，而其他桿端位移為零時所引起的第r個桿端力。單元剛度矩陣是一個對稱矩陣，即krs=ksr，這可由反力互等定理證明；單元剛度矩陣還是一個奇異矩陣，這是由于單元中包含剛體位移。局部向整體坐標的變換：De=TTDeFe =TTFek =TTkT式中，T為一6´6的坐標變換矩陣。假設自整體坐標Ox軸沿轉角正方向（即順時針方向

51、）轉到局部坐標Ox軸的角度為a，則F=kDe坐標變換矩陣T是一個正交矩陣，即T-1=TT，故結點位移和結點力由整體向局部坐標的變換式為De=TDe，Fe=TFe2.8.2 考慮剪切變形的梁單元基本假定：垂直于梁軸線的截面在變形后仍保持為平面，但不再垂直于變形后的軸線。撓度由兩部分組成：彎曲變形部分和剪切變形部分，即v=vb+vs；dxdxqg桿軸線斜率：dv/dx=q+g（q：截面的轉角，由彎曲變形引起；g：剪切變形，截面與桿軸夾角的改變）圖2.31 發生彎曲和剪切變形的直梁微段(1) 在經典梁單元基礎上引入剪切變形的影響vb用三次式插值，參見上一節。vs用線性插值：vs = N7vsi+ N

52、8vsj。由此可導得單元剛度矩陣為式中，為剪切影響系數，k為截面剪應力不均勻分布修正系數，對矩形截面k=1.2，圓形截面k=10/9。(2) Timoshenko梁單元對u、v、q均獨立插值：u=N1ui+N2 u2， v= N1v1+ N2v2， q=N1q1+ N2q2式中q：截面的轉角，dv/dx：桿軸線的斜率。當l/h®¥時只能得到零解，即將出現剪切自鎖（Shear locking）。這是由于v、q同階插值（實際上放大了剪切應變項的量級），故dv/dx與q不同階，從而導致剪應變g=dv/dx-q=0不能處處滿足（dv/dx為常數，q為一次式），除非q=常數，意味著梁不能發生彎曲。解決方法：減縮積分數值積分時采用比精確積分要求少的積分點數，例如對兩結點梁單元采用一點積分；減縮積分后ks中的l2/3和l2/6項均改為l2/4。假設剪切應變對剪應變g另行假定插值形式；替代插值函數計算剪應變時，對q采用低一階的插值函數，如兩結點梁單元其插值函數為常數1/2，即。兩結點Timoshenko梁單元包含橫向剛體位移（v=c）和剛體轉動（q=dv/dx=c），包含常剪切應變（q=0，dv/dx=g=c）