有限單元法及程序設計第一章緒論第二章平面桁架有限元分析第三章平面剛架有限元分析第四章平面問題有限元分析有限單元法及程序設計第一章緒論§1.1有限單元法的概念§1.2有限單元法的基本步驟§1.3常用有限單元分析軟件第一章緒論數值方法（模擬）有限元方法邊界元方法有限差分方法有限體積法無網格方法§1.1有限單元法的概念基本思路：借助數學和力學知識,利用計算機技術解決工程技術問題有限元分析是利用數學近似分析方法對真實物理系統（幾何、載荷工況）進行模擬，利用簡單而又相互作用的元素，即單元，用有限數量的未知量去逼近無限未知量的真實系統。有限單元法(FEM)是20世紀50年代以來隨著計算機的廣泛應用而發展起來的一種現代數值解法。該方法首先應用在連續力學領域——飛機結構靜、動態特性分析中。隨后很快就廣泛應用于求解傳導、電磁場、流體力學等連續性問題。Rxyo真實半圓周長：近似半圓周長：數學：離散化的手段，將偏微分方程、變分方程變成代數方程；計算曲線長度：將曲線離散成N段，每段用直線長度代替曲線長度，直線長度總和即可曲線長度，當N為無窮大時，即為精確解。§1.1有限單元法的概念單元數計算值理論值誤差(%)123.1415936.3422.828433.141599.9732.976073.141595.2743.035283.141593.3853.065963.141592.4163.084253.141591.8373.096193.141591.4583.10453.141591.1893.110543.141590.99103.115113.141590.841003.140763.141590.032003.14133.141590.013003.141433.141590.01§1.1有限單元法的概念圓周率的計算方法：單元數計算值單元數計算值單元數計算值121003.14076059100003.1415918220398022.8284271252003.141298567200003.1415923595933632.9760677433003.14143259300003.1415924935585943.035276184003.141488694400003.1415925496465253.0659554365003.141518269500003.1415925792140963.084252856003.141536068600003.1415925970102973.0961949547003.14154775700003.1415926086905583.1044960688003.141555901800003.1415926168403593.1105417379003.141561854900003.14159262279185103.11510595110003.1415663561000003.14159262729364小數點后面20位：3.14159265358979323846§1.1有限單元法的概念力學：離散化的手段，利用單元力學特性，將連續體離散成有限個單元的組合結構；齒輪的應力分布：將齒輪分割成很多單元，由結點相連并傳遞求解信息，用有限個連續體的組合代替原來的齒輪結構?！?.1有限單元法的概念電動鏟運機§1.1有限單元法的概念液壓挖掘機動臂§1.1有限單元法的概念在建的大連國際貿易中心大廈(78層，342米)63863個梁柱單元；34180個結點；§1.1有限單元法的概念在建的大連市體育中心§1.1有限單元法的概念§1.1有限單元法的概念溪洛渡拱壩(高297米)456個壩體單元，1040個地基單元，2163個結點；§1.1有限單元法的概念§1.1有限單元法的概念導管架碼頭結構§1.1有限單元法的概念高樁碼頭高樁碼頭下部結構高樁碼頭上部結構§1.1有限單元法的概念游艇碼頭結構§1.1有限單元法的概念重力式碼頭結構有限元模型是真實系統理想化的數學抽象。真實系統有限元模型有限元模型§1.1有限單元法的概念單元網格劃分中每一個小的塊體結點確定單元形狀、單元之間相互聯結的點結點力單元上結點處的結構內力載荷作用在單元結點上的外力（集中力、分布力）約束限制某些結點的某些自由度

單元單元載荷結點結點力約束有限元單元模型中幾個重要概念§1.1有限單元法的概念§1.2有限單元法的基本步驟一、結構離散化將求解域離散成單元表示的組合體，單元以結點相連；二、選擇插值函數（位移模式）選擇插值函數來表達單元內場變量的變化規律三、形成單元性質的矩陣方程場變量：標量向量力、位移張量應力、應變單元剛度方程：單元結點力與單元結點位移之間的關系；§1.2有限單元法的基本步驟四、形成整體系統的矩陣方程綜合所有單元性質的矩陣，形成整體系統的矩陣方程；五、約束處理，求解系統方程引入邊界條件，求解結點上的未知場變量；六、其他參數的計算利用已經求解的場變量，計算其他場變量；§1.3常用的有限元分析軟件(1)ANSYS功能強大、模塊多、比較通用；土木工程：CivilFEM是商業化比較早的一個軟件（收購了一些其他軟件公司）(2)MSC產品系列多、通用軟件、二次開發功能強；土木工程：MSC.MARC系列，Patran,MSCNastran及AdamsAdams------多體動力學；Actran------聲學仿真；Easy5-----控制仿真工具；Marc------非線性；SimXpert------多學科仿真；MSCNastran------結構化與多學科FEA；Dytran------顯式動力學與流固耦合；MSCFatigue------基于FE的耐久性仿真工具；Sinda------高級熱分析解決方案；Digimat------非線性，多尺度的材料與結構建模平臺；SimDesigner------CAD嵌入式多學科仿真；Patran------有限元分析解決方案；SimManager------仿真數據和流程管理§1.3常用的有限元分析軟件(3)ADINA復雜非線性問題與復雜動力問題；具有直接求解和迭代求解兩種流固耦合分析方法AutomaticDynamicIncrementalNonlinearAnalysis(4)ABAQUS（專注結構分析）主要應用于結構動力彈塑性分析；二次開發功能；HKS公司產品§1.3常用的有限元分析軟件(5)SAP2000SystemsApplicationsandProductsinDataProcessing通用結構分析

第二章平面桁架有限元分析及程序設計§2.1平面桁架單元的離散§2.2平面桁架單元分析§2.3結點平衡與整體剛度矩陣的集成§2.4邊界條件的處理§2.5單元內力與支座反力的計算§2.6平面桁架有限元程序設計有限單元法及程序設計解題方法ⅠⅠ方法1：節點法BaFP1ACDaaaFP方法2：截面法靜定桁架回顧第二章平面桁架有限元分析及程序設計解題方法：力法和位移法超靜定桁架P1234第二章平面桁架有限元分析及程序設計如圖所示桁架，求各桿軸力。力的平衡條件：位移的協調方程：PN1N1N21桿和3桿位移：2桿位移：超靜定桁架第二章平面桁架有限元分析及程序設計1桿軸力豎向分量：2桿軸力：式中：和為桿件的剛度系數；物理意義：4點產生單位位移，桿端產生的豎向桿端力；由桿件的物理性質和幾何性質決定；V4為第4節點豎向位移P1234第二章平面桁架有限元分析及程序設計超靜定桁架代入平衡方程：結構的整體剛度系數位移法求解超靜定結構。離散原則：每個結點離散后還是一個結點，每個桿件離散后變成一個單元1結構的離散化：盡量將結構離散成數量最少的等截面直桿單元23456⑤①②③④⑥⑦⑧⑨§2.1平面桁架單元的離散9個單元，6個結點123456⑤①②③④⑥⑦⑧⑨78⑩15121413111616個單元，8個結點§2.2.1

局部坐標系下的單元剛度矩陣局部坐標系的建立iE，A，lje■

軸：沿單元的桿軸方向；

§2.2平面桁架的單元分析■

軸：從軸逆時針旋轉90°?！?/p>

原點：以第一個結點為坐標原點；

桿端位移：ijeijeije桿端力：符號：與坐標系的方向一致為正，反之為負。單元右端桿端力：單元左端桿端力：單元應力：單元應變：右結點固定結點位移：左結點固定桿的受力分為兩種情況：§2.2平面桁架的單元分析任意情況（左右結點均有變形）即為以上兩種狀態的疊加：桿端力為：式中為單元剛度矩陣(局部坐標系)桿單元軸力為：式中為單元應力(廣義)矩陣；§2.2平面桁架的單元分析單元桿端力方程：桿端位移：ijeije桿端力：單元軸力：§2.2平面桁架的單元分析桿端位移和桿端力符號：與坐標系的方向一致為正，反之為負。桿端力：xy§2.2平面桁架的單元分析§2.2.2

整體坐標系下的單元剛度矩陣若局部坐標系與整體坐標系重合，則整體坐標系下的單元剛度矩陣與局部坐標下的單元剛度矩陣相同。若局部坐標系與整體坐標系不重合，如下圖所示：桿端位移：桿端力：桿端位移：i

結點：j

結點：§2.2平面桁架的單元分析設桿件的長度為l，則：兩邊微分：由于桿件的變形產生位移：因此，桿件應變為：桿件軸力為：符號：桿件軸力以拉為正，壓為負。桿件的結點力為：因此，桿件結點力向量為：式中是整體坐標系下的單元剛度矩陣；§2.2平面桁架的單元分析寫成分塊矩陣形式：式中：§2.2平面桁架的單元分析（1）單元剛度系數kij的意義j自由度(結點)產生的單位桿端位移引起的i自由度(結點)的桿端力（2）單元剛度矩陣是對稱矩陣反力互等定理式中：桿件單元的應力矩陣為：單元剛度矩陣的性質§2.2平面桁架的單元分析（3）單元剛度矩陣一般是不可逆的§2.2.3單元坐標轉換矩陣§2.2平面桁架的單元分析取任意桿件,建立如圖所示的局部坐標系：桿端力：桿端位移：§2.2.3單元坐標轉換矩陣xyxy§2.2平面桁架的單元分析桿端力：桿端位移：在上圖中,建立如圖所示的整體坐標系：以i結點為例：同理，對于j結點：§2.2平面桁架的單元分析寫成矩陣形式：因此：其中，[T]為轉換矩陣：轉換矩陣的性質轉換矩陣是正交矩陣；同理，位移也存在轉換關系：代入局部坐標系下的剛度方程：§2.2平面桁架的單元分析與利用微分得到的單元在總體坐標下的剛度方程相同xyF1、對總體結點位移和單元進行編碼；2、單元局部坐標系下的剛度矩陣;①②③123例：如圖所示平面桁架，桿長為l，截面積為A，求三個單元在整體坐標系下的剛度矩陣?！?.2平面桁架的單元分析3、單元①整體坐標系分析:解：單元①整體坐標系下的剛度矩陣為:4、單元②整體坐標系分析:5、單元③整體坐標系分析:xyF①②③123例：如圖所示平面桁架，桿長為l，截面積為A，求結構的剛度矩陣。§2.3結點平衡與整體剛度矩陣的集成1、單元整體坐標系下剛度矩陣分塊解：§2.3.1結點的平衡方程§2.3結點平衡與整體剛度矩陣的集成xyF1①ij③ij2、結點1的平衡方程：結點1的受力狀態為（如右圖）:結點1的平衡條件為:由單元③的剛度方程：由單元①的剛度方程：§2.3結點平衡與整體剛度矩陣的集成代入結點1的平衡條件:3、結點2的平衡方程：同理，結點2的平衡條件為:由單元①的剛度方程：由單元②的剛度方程：§2.3結點平衡與整體剛度矩陣的集成代入結點2的平衡條件:4、結點3的平衡方程：同理，結點3的平衡條件為:由單元③的剛度方程：由單元②的剛度方程：§2.3結點平衡與整體剛度矩陣的集成代入結點3的平衡條件:5、系統的平衡方程：§2.3結點平衡與整體剛度矩陣的集成6、結構整體剛度方程寫成矩陣形式，即可得到結構的整體剛度方程其中，[K]為結構的整體剛度矩陣；§2.3.2

整體剛度矩陣的集成步驟1、定位單元結點編號

2、累加整體結點編號

單元剛度系數

整體剛度系數§2.3結點平衡與整體剛度矩陣的集成單元定位向量xyF①②③123§2.3結點平衡與整體剛度矩陣的集成例：求上例平面桁架結構的整體剛度矩陣；1、定位單元①：1

2

3123(2)(1)(2)

(1)2、累加§2.3結點平衡與整體剛度矩陣的集成1、定位單元③：(2)(1)(2)

(1)2、累加1、定位單元②：2、累加1

2

3123

(1)

(2)(1)(2)1

2

3123§2.3結點平衡與整體剛度矩陣的集成結構的整體剛度矩陣為:將每個字塊展開，結構的整體剛度矩陣為:下標表示自由度編號下標表示結點編號§2.3結點平衡與整體剛度矩陣的集成1、整體剛度矩陣的性質§2.3.2整體剛度矩陣的集成方法（1）剛度系數Kij的意義分塊矩陣：（2）單元剛度矩陣是對稱矩陣反力互等定理j結點產生的單位桿端位移引起的i結點的桿端力；j自由度產生的單位桿端位移引起的i自由度的桿端力；不分塊矩陣：§2.3結點平衡與整體剛度矩陣的集成思考題：如何由剛度矩陣元素的意義確定整體剛度矩陣每個元素的組成？§2.3結點平衡與整體剛度矩陣的集成1234⑤①②③④⑥練習題：利用整體剛度矩陣的意義確定以下桁架剛度矩陣元素（分塊）的組成?！?.3.2

整體剛度矩陣的集成步驟1、定位單元自由度編號

2、累加整體自由度編號

§2.3結點平衡與整體剛度矩陣的集成自由度定位向量單元剛度系數

整體剛度系數練習題：利用整體剛度矩陣的意義確定以下桁架剛度矩陣元素（自由度）的組成?！?.4.1

結點邊界條件§2.4邊界條件的處理結點可以自由變形，整體結點力等于對應的外荷載。1、自由變形的結點2、約束結點或給定了結點位移的數值約束結點：給定結點位移：§2.4.2

邊界條件的處理方法§2.4邊界條件的處理1）劃行劃列法處理方法：若第i個自由度位移為零，則將總剛第i行和第i列劃掉，剛度矩陣相應降低一階。000優點：簡單易行，矩陣降階，減小計算工作量；缺點：矩陣行列、位移、荷載向量需重新編號，程序實現比較復雜；只適用于約束結點情況；§2.4邊界條件的處理2）0、1置換法（填0置1法）處理方法：將與約束自由度對應整體剛度矩陣對角線元素全部置換成1，相應行和列其他元素置換成0，將同一行荷載分量置換成0。000111000000000000000000000000000適用條件：只適用于約束結點，不適用給定位移邊界條件；§2.4邊界條件的處理3）乘大數法處理方法：將與約束自由度對應整體剛度矩陣對角線元素乘以一個大數N(1010-1015)，將同一行荷載分量置換成N與對角線元素的乘積與給定位移之積。NNNNK33b1NK44b2NK66b3優點：處理工作量小，適用于給定位移情況；000§2.5單元軸力及支座反力的計算單元軸力：(1)整體坐標系下的單元桿端位移(2)整體坐標系下的單元軸力§2.5.1

單元軸力的計算位移和荷載向量分解為自由結點和約束結點兩部分，剛度矩陣相應分塊，如下所示：§2.5.2

支座反力的計算其中：為自由位移;為約束位移;為外荷載;為約束自由度結點力;因此：為支座反力;為約束自由度結點荷載;其中：§2.5單元軸力及支座反力的計算1、計算分析題平面桁架結構網格如圖所示，已知EA=1500kN，采用乘大數方法引入支撐條件，試求后處理法引入支撐條件后的滿陣存貯的整體剛度矩陣[K]。結構坐標系下單元剛度矩陣計算公式如下：

例題式中：321①②③3myx4m2、計算分析題按照有限元法的計算步驟，求圖示桁架結構各桿軸力。已知：EA=10kn。結構坐標系下的單元剛度矩陣計算公式如下：：

例題式中：§2.6程序設計程序設計原則完整性擴充性兼容性邏輯性可讀性可維護性模塊化§2.6程序設計程序設計流程圖（程序框圖）什么是流程圖？

“程序流程圖”常簡稱為“流程圖”，是一種傳統的算法表示法，程序流程圖是人們對解決問題的方法、思路或算法的一種描述。它利用圖形化的符號框來代表各種不同性質的操作，并用流程線來連接這些操作。2、如何畫流程圖編碼和單元測試這個階段的任務是程序員根據目標系統的性質和實際環境，選取一種適當的高級程序設計語言(必要時用匯編語言)，把詳細設計的結果翻譯成用選定的語言書寫的程序，并且仔細測試編寫出的每一個模塊。

程序員在書寫程序模塊時，應使它的可讀性、可理解性和可維護性良好。綜合測試這個階段的任務是通過各種類型的測試，使軟件達到預定的要求。

最基本的測試是集成測試和驗收測試。集成測試是根據設計的軟件結構，把經單元測試的模塊按某種選定的策略裝配起來，在裝配過程中對程序進行必要的測試。驗收測試是按照需求規格說明書的規定，由用戶對目標系統進行驗收。

通過對軟件測試結果的分析可以預測軟件的可靠性；反之，根據對軟件可靠性的要求也可以決定測試和調試過程什么時候可以結束。

在進行測試的過程中，應該用正式的文檔把測試計劃、詳細測試方案以及實際測試結果保存下來，作為軟件配置的一部分。

§2.6程序設計§2.6.1程序框圖輸入數據單元局部剛度坐標轉換矩陣單元整體剛度集成整體剛度矩陣元素約束條件處理、解方程計算單元軸力、約束反力單元循環包括單元、結點、材料、荷載、約束數據§2.6

程序設計§2.6.2

程序說明1、總體剛度矩陣的半帶寬存儲總體剛度矩陣：對稱稀疏矩陣；00半帶寬：總體剛度矩陣集成：§2.6

程序設計行號：000列號：主對角線半帶寬存儲下三角：§2.6

程序設計2、先處理法處理邊界條件單元定位向量：xyF①②③123(1,2)(3,4)(5,6)后處理：xyF①②③123(1,2)(0,0)(3,0)劃行劃列法先處理法：單元定位向量：

voidforce(){int

i,j,ie,m;floatdx,dy,dz,l,cx,cy,cz,ea,w[7];

for(ie=1;ie<=ne;ie++){

i=jm[ie][1];

j=jm[ie][2];

m=jm[ie][0];w[1]=f[2*i-2];

w[2]=f[2*i-1];

w[3]=f[2*j-2];

w[4]=f[2*j-1];

dx=xy[j][1]-xy[i][1];

dy=xy[j][2]-xy[i][2];

例題3、程序說明下面為一個平面桁架計算程序段，試在左端有編號的程序右面寫出其注釋。（10分）

l=sqrt(dx*dx+dy*dy);

cx=dx/l;cy=dy/l;ea=EA[m]/l;

dx=w[3]-w[1];

dy=w[4]-w[2];

l=ea*(cx*dx+cy*dy