麻将胡了pg大奖视频

書名：材料成形過程數值模擬（第二版）
作者名：傅建肖兵
本章字數： 6982字
更新時間： 2020-07-01 17:30:30

2.1.2　有限元方程的建立與應用

針對不同工程問題構建的有限元方程（有限元模型）是有限元法應用的基礎，而變分法和加權余量法是建立有限元方程的兩類常用數學方法。本節將以求解平面彈性力學剛度問題和穩態熱傳導問題為例，分別介紹這兩類方法的實施要點。

2.1.2.1　預備知識

（1）彈性力學基本方程

①平衡方程　當彈性體中任一質點上的應力達到平衡時，有

Lσ+b=0　　在（Ω域內）　　（2-1）

式中　L——微分算子，對于平面問題；

　σ——應力，對于平面問題σ^T={σ_x　σ_y　τ_xy}；

　b——體積力（一般為重力）向量，對于平面問題b^T={b_x　b_y}。

②幾何方程　幾何方程表征質點位移與應變之間的關系

ε=Lu　　（在Ω域內）　　（2-2）

式中　ε——應變，對于平面問題ε^T={ε_x　ε_y　γ_xy}；

　u——質點位移矢量，對于平面問題u^T={u_x　u_y}。

③本構方程　即材料的應力-應變關系

σ=Dε　　（在Ω域內）　　（2-3）

式中　D——彈性矩陣，對于平面問題；

E、μ——材料的彈性模量和泊松比。

④邊界條件

Γ=Γ_F+Γ_u　　（在Ω的邊界上）　　（2-4）

式中　Γ_F——面力和集中力邊界；

　Γ_u——位移邊界。

（2）變分法

變分法的基礎是變分原理，而變分原理是求解連續介質問題的常用數學方法之一。

例如：一維穩態熱傳導的定解問題

　　（2-5）

式中　　　 ?——未知溫度場函數；

　 k——熱導率；

　Q（x）——沿x方向分布的熱載荷，。

式（2-5）中的溫度場?可以借助傅里葉積分求得。

如果采用變分原理求解這類定解問題，則應首先建立定解問題的積分形式

　　（2-6）

式中　　u——未知函數；

　F、E——特定算子；

Ω——連續求解域；

Γ——Ω的邊界；

Π——泛函數（即未知函數u的函數）。

此時，如果連續介質問題有解u，則解u必定使泛函Π對于微小變化δu取駐值（極值），即泛函的“一階變分”等于零

δΠ=0　　（2-7）

這就是所謂求解連續介質問題的變分原理。

可以證明：用微分方程加邊界條件求解連續介質問題同用約束或非約束變分法求解連續介質問題等價。即一方面滿足微分方程及其邊界條件的函數將使泛函取得駐值（極值）；另一方面，從變分角度看，使泛函取得駐值（極值）的函數正是滿足連續介質問題的微分方程及其邊界條件的解。

應用變分法建立有限元方程的目的，就是將求微分方程的定解問題轉變成求泛函的駐值問題，以方便試探函數的分片插值和分片積分（見圖2-2）。

（3）虛位移方程（位移變分方程）

虛位移方程是變分原理在求解彈性力學問題中的具體應用，它等價于幾何微分方程和應力邊界條件。所謂虛位移是指：在約束條件允許的范圍內，彈性體內質點可能發生的任意微小位移。虛位移的產生與彈性體所受外力及其時間無關。

彈性體在外力的作用下發生變形，表明外力對彈性體做了功。若不考慮變形過程中的熱損失、彈性體的動能和外界阻尼，則外力功將全部轉換為貯存于彈性體內的位能（應變能）。根據能量守恒定律，有

　　（2-8）

式中　δε——虛應變（即對應虛位移的任意可能應變），δε^T={δε_xδε_yδε_zδγ_yzδγ_zxδγ_xy}；

　δu——虛位移，δu^T={δu　δv　δw}；

　Ω——彈性體內部；

　Γ_f——Ω的面力邊界（假設邊界上無離散的集中力）；

，b——面力和體積力。

式（2-8）表明，外力（包括Ω內的體積力和邊界上的面力）使彈性體內質點產生虛位移所做的功等于彈性體內部虛應變產生的能量。

2.1.2.2　平面剛度問題

（1）離散處理

假設分析對象（構件）的厚度尺寸非常小，可以近似將其處理成厚度為常數的平面問題。用三節點三角形單元離散該對象及其邊界條件［圖2-4（a）］，并從中任取一子域進行單元分析［圖2-4（b）］。

圖2-4　離散后的分析對象（a）與單元結構（b）

圖2-4中，單元節點i、j、m按逆時針排布。

（2）單元剛度分析

①單元位移與單元形函數　一個三節點三角形單元共有12個自由度（6個節點位移分量和6個節點轉動分量），在線彈性范圍內，6個節點轉動分量可忽略不計，由此給出單元位移a^e的向量表達式如下

　　（2-9）

求解離散后的平面剛度問題存在兩種情況：a.位移分量是節點坐標的已知函數，對此可直接利用單元節點位移分量求出單元應變分量（幾何方程），再由單元應變分量求出單元應力分量（本構方程），最后綜合起來便可得到整個平面剛度問題的解；b.只有幾個節點的位移分量已知，不能直接求出單元應變分量和應力分量。對于后者，為了利用節點位移表示單元應變和應力，必須構造一個位移模式（位移函數）。理論上，定義于某一閉區域內的任意函數總可被一個多項式逼近，所以位移函數常常取為多項式。現選用一次多項式作為三節點三角形單元的位移模式（位移插值函數），于是有

　　（2-10）

式中　 u、v——單元內任意一質點的位移分量；

　a₁～a₆——待定系數。

將三角形的三個節點坐標代入式（2-10），可得

　　（2-11）

應用克萊姆法則求解式（2-11），并化簡

　　（2-12）

式中　　　　u_k、v_k——節點位移，k=i，j，m（下同）；

N_k（x，y）——三節點三角形單元的形函數（一個與單元類型和單元內部節點坐標有關的連續函數），

——三角形單元的面積；

a_k，b_k，c_k——與節點坐標有關的常數，a_k=x_jy_m-x_my_j，b_k=y_j-y_m，c_k=-x_j+x_m。

用矩陣表示式（2-12），得：

　　（2-13）

式中　I——單位矩陣；

　N——形函數矩陣；

　a^e——單元節點位移列矩陣。

式（2-12）和式（2-13）即為采用多項式插值構造的三節點三角形單元的位移模式，表明只要知道了節點位移，就能借助單元形函數求得單元內任意一質點的位移。換句話說，節點位移通過形函數控制整個單元內部的質點位移分布。

單元形函數具有以下兩條基本性質。

a.在任一節點上，形函數滿足

b.針對單元內任一質點，有

N_i（x，y）+N_j（x，y）+N_m（x，y）=1

②單元應變矩陣與應力矩陣　將式（2-13）代入幾何方程（2-2），可得單元應變表達式

　　（2-14）

式中　B——單元應變矩陣，其分塊子矩陣為：

　　（2-15）

再將單元應變表達式（2-14）代入本構方程式（2-3），可得單元應力表達式

σ=Dε=DBa^e=sa^e　　（2-16）

式中　s=DB=D［B_i　B_j　B_m］=［s_i　s_j　s_m］——單元應力矩陣，其分塊矩陣為：

　　（2-17）

式中　E₀、μ₀——材料常數，對于平面應力問題E₀=E，μ₀=μ，對于平面應變問題。

E，μ為材料的彈性模量和泊松。

由于應變矩陣B中的每一個非零元素均是由節點坐標決定的常數，而節點坐標為定值，所以應變矩陣B是常量矩陣。同理，應力矩陣s也是常量矩陣。這表明由式（2-14）和式（2-16）計算獲得的單元內各質點的應變值相同，應力值亦相同。可以證明，在彈性力學范圍內，由節點位移引發的應變和應力主要通過相鄰單元的邊界節點進行傳遞。

③單元剛度矩陣與單元剛度方程　現利用2.1.2.1小節介紹的虛位移原理建立單元剛度方程。假設：在來自相鄰單元“外力”q^e的作用下，單元節點i，j，m發生了虛位移δa^e，虛位移引起虛應變δε，由此得到單元的虛位移方程：

　　（2-18）

式中　t——單元厚度。

因虛位移是任意的，由此得單元剛度方程（單元剛度模型）

q^e=k^ea^e　　（2-19）

式中　k^e——單元剛度矩陣，表明單元e中各節點產生單位位移而引發（或需要）的節點力

　　（2-20）

其中

單元剛度矩陣k^e具有以下特征。

a.對稱性，即k_rs=k_sr；

b.奇異性，即單元剛度矩陣不存在其逆矩陣；

c.各行、各列元素之和恒為0（因在力的作用下，整個單元處于平衡狀態）；

d.主對角元素恒為正，即k_ii>0，表明節點位移與施加其上的節點力同向。

④單元等效節點力　作用在單元上的外力q^e由等效節點力構成，即

　　（2-21）

式中　，分別為體積力b、面力和集中力p_c產生的等效節點力。

所謂等效節點力是指：由作用在單元上的體積力、面力和集中力按照能量等效原則（即原載荷與等效節點載荷在虛位移上做功相等原則）移至節點上而形成的節點力。

圖2-5是作用在單元上的外力舉例，表2-1是對應于圖2-5的等效節點力表達式。

圖2-5　體積力、面力和集中力舉例

表2-1　對應于圖2-5的等效節點力

（3）整體剛度分析

①建立總剛度矩陣與總剛度方程　將式（2-19）改寫成如下形式

　　（2-22）

式（2-22）表明，當單元中任一節點發生位移時，都將在該節點處產生節點力（實際上是由節點力引起位移），且大小等于單元中各節點位移所引起的節點力之和。

因為在被離散對象的整體結構中，一個節點往往為若干比鄰單元所共有，根據線性疊加原理，作用在該節點上的力應等于所有比鄰單元的等效節點力之和，即

　　（2-23）

式中　 Q_i——作用在節點i上的合力；

——比鄰單元施加給i的等效節點力之和；

　 n_e——比鄰單元個數。

將式（2-22）代入式（2-23），并取r=i，得到當節點i處于平衡狀態時的合力表達式

　　（2-24）

假設被離散對象（區域）有n個節點，于是該對象中所有節點的平衡方程為

　　（2-25）

用矩陣表示，有

Ka=Q　　（2-26）

式中　K——總剛度矩陣，；

　Q——總節點載荷列矩陣，；

　a——總節點位移列矩陣。

式（2-26）為利用虛位移方程和節點力疊加原理建立起來的平面剛度問題的整體平衡方程（整體剛度模型），即用有限元格式表示的平面問題總剛度方程。

總剛度矩陣K具有以下性質。

a.對稱性　總剛度矩陣由單元剛度矩陣疊加形成，所以與單元剛度矩陣一樣具有對稱性，即K=K^T。

b.稀疏性　由于任一節點僅與少數節點和單元比鄰，相對于該節點，其無關節點在K中的對應元素為零。因此，存在于K中的大量零元素使總剛度矩陣具有稀疏性（稀疏矩陣）。

c.帶狀性　總剛度矩陣中的所有非零元素均集中分布在主對角線附近，從而形成所謂的帶狀矩陣。

d.奇異性　對于無任何節點約束或約束不足的結構件，其總剛度矩陣奇異（|K|=0），物理上表現為在力的作用下，結構件整體作剛性運動，即式（2-26）的解非唯一。

總剛度矩陣的上述特性，為其在計算機中的存儲與處理，以及解的唯一性判別和計算方法的選取提供了良好的數學依據。

②位移邊界條件　由于總剛度矩陣具有奇異性，因此，在求解總剛度方程式（2-26）之前必須設置一些節點位移約束，以防止結構件產生整體剛性運動。節點的位移約束一般施加在被離散對象的邊界上，故稱之為位移邊界條件，見圖2-4（a）。邊界上節點的位移約束通常分為兩大類。

a.零位移約束　設某一節點沿某一方向的位移為零（例如a_m=0），則平面構件的總剛度方程變成

　　（2-27）

觀察式（2-27）可以發現：總剛度矩陣中與零位移節點a_m對應的主對角元素k_mm為1，其余相關元素（即_足標中含有m的k元素）均為0；載荷列矩陣中與零位移節點對應的元素Q_m為0。

b.非零已知位移約束　已知某一節點沿某約束方向位移為，其平面構件的總剛度方程轉變成

　　（2-28）

此時，式（2-28）表現為：總剛度矩陣中與已知節點位移對應的主對角元素k_mm乘以一個足夠大的正數α（例如，α=10²⁰），相關列元素（即第二個下標是m的k元素）均為零，并且載荷列矩陣中的Q_m被所代替。

實際上，兩類位移邊界條件常常同時存在，所以可根據情況將式（2-27）和式（2-28）組合在一塊進行化簡。

③設置載荷邊界條件　設置載荷邊界條件實際上是為總載荷列矩陣中的各分量元素賦初值。結合位移邊界條件設置可以簡化賦值過程，即只針對未經式（2-27）和式（2-28）處理的載荷元素賦值。需要注意的是：除零位移和已知位移節點外，載荷中的等效體積力應平均加載到所有剩余節點上，等效面力應加載到面力邊界節點上，而等效集中力則應根據情況加載到相應單元的節點上。這樣處理后的總載荷列矩陣中，有的載荷分量可能為零（對應零位移節點），有的載荷分量可能是2～3種等效節點力的疊加（例如，某邊界節點既受等效面力作用，又受等效體積力作用）。

（4）求解總剛度方程

代入約束邊界條件和載荷邊界條件的總剛度方程可寫成

　　（2-29）

式中　——經約束邊界條件和載荷邊界條件處理后的總剛度矩陣與總節點載荷列矩陣。

式（2-29）是一個關于節點位移分量的線性方程組，利用迭代法、高斯法、波前法或帶寬法等數值方法可以求出節點位移列矩陣a中的各分量，然后代入幾何方程求解應變分量，再代入本構方程求解應力分量，于是便獲得已知邊界條件下平面剛度問題的全部近似解。

（5）平面剛度問題應用舉例

①離散處理　假設：將如圖2-6所示平面構件離散成8個三節點三角形單元；已知面載荷P，其等效結點力F=P/3均布在4、6、8三個節點上且垂直于x-z坐標面；各結點在z軸方向上的位移忽略不計（z方向尺寸很小，可固定設置為t），并忽略各節點的轉動；構件受力時，節點1和9的x位移分量為0，1、9、2、10的y位移分量也為0。

圖2-6　平面剛度問題舉例

②單元剛度分析　例如：針對圖2-6中的單元1建立節點1、4、2的平衡方程（單元剛度方程）

式中　k₁₁～k₆₆——2×2階子矩陣，其子矩陣中的各元素與材料的彈性模量E、泊松系數μ、節點坐標（x_i，y_i），以及單元厚度t有關，見式（2-20）。

③整體剛度分析　將所有單元分析結果依次組合，得到圖2-6平面構件的總剛度方程（有限元方程）。

{F}=［K］{δ}

式中　{F}={F₁F₂…F₁₀}^T={F₁_xF₁_yF₂_xF₂_y…F₁₀_xF₁₀_y}^T；

{δ}={δ₁δ₂…δ₁₀}^T={δ₁_xδ₁_yδ₂_xδ₂_y…δ₁₀_xδ₁₀_y}^T；

④設置邊界條件　將位移約束δ₁_x=δ₉_x=δ₁_y=δ₉_y=δ₂_y=δ₁₀_y=0和等效節點力F₄_y=F₆_y=F₈_y=P/3分別代入總剛度方程的位移列矩陣和載荷列矩陣，整理并化簡。

⑤計算求解　求解總剛度方程中剩余節點的位移分量，然后再分別利用幾何方程和本構方程求解單元應變與單元應力。

⑥分析求解結果　可將計算獲得的單元最大等效應力值與構件的許用強度進行比較，將節點最大位移值與構件的許用擾度（剛度）進行比較，以確定構件在面載荷P的作用下是否失效。取最大等效應力值和最大位移值作為失效判斷的基本依據是：單元應力和節點位移在離散區域內不會發生突變，這與實際工程問題相吻合。

2.1.2.3　平面變溫問題

變化的溫度場會在彈性體內部誘導熱應力，這種熱應力主要來自于溫度場變化所引起的彈性體不均勻熱脹冷縮。平面變溫問題的本構方程可表示為

σ=D（ε-ε₀）　在Ω域內　　（2-30）

式中　ε₀——溫度變化引起的溫度應變，ε₀=α（?-?₀）{1　1　0}^T；

　α——材料的線脹系數（假設各向同性）；

　?——溫度場；

　?₀——初始溫度場；

　ε——力載荷引起的應變；

　σ——應力；

　D——彈性矩陣。

同無溫度場變化的本構方程相比，式（2-30）用（ε-ε₀）替代了式（2-3）中的ε。

將單元應變方程ε=Ba^e代入式（2-30），得平面變溫問題的單元本構方程

σ^e=DBa^e-Dε₀　　（2-31）

若不考慮外力作用，則變溫條件下單元位能的泛函表達式為

　　　　（2-32）

式中　——單元剛度矩陣；

　——等效節點熱載荷。

總位能的泛函表達式

　　（2-33）

式中　——總剛度矩陣；

　——總變溫載荷列矩陣。

根據能量泛函取駐值條件（最小位能原理）δΠ_p=0，對式（2-33）求一階變分并化簡，得

Ka=Q　　（2-34）

式（2-34）即為求解平面變溫引起節點位移的有限元剛度方程。需要注意的是：求解變溫問題的平面應力時，應將式（2-34）的計算結果代回變溫問題的本構方程式（2-30）。利用泛函變分取駐值法建立限元方程是變分原理在彈性力學中的又一應用。同樣，方程（2-34）中的Q也可包括表面、體積、集中等載荷。

2.1.2.4　平面穩態熱傳導問題

（1）熱傳導基本方程與邊界條件

　在Ω域內　　（2-35）

　　　　　　　　　　　　　　　　　　　在Γ₁邊界上　　（2-36）

　在Γ₂邊界上　　（2-37）

　在Γ₃邊界上　　（2-38）

式中　　　　ρ——密度；

　　c——比熱；

　　t——時間；

k_x、k_y、k_z——沿x、y、z方向的傳熱系數；

　　Q——固體材料內部的熱源密度，Q=Q（x，y，z，t）；

n_x、n_y、n_z——邊界外法線的方向余弦；

　　——Γ₁邊界上給定溫度，；

　　q——Γ₂邊界上給定熱流密度（當q=0時，絕熱邊界），q=q（Γ，t）；

h（?_a-?）——Γ₃邊界上給定界面傳熱；

　　h——界面傳熱系數；

　　?——待求溫度場（即與空間和時間相關的溫度函數）；

　　?_a——環境溫度（自然對流條件下）或邊界層絕熱溫度（強制對流條件下），?_a=?_a（Γ，t）；

　　Γ——域Ω的邊界，Γ=Γ₁+Γ₂+Γ₃。

熱傳導基本方程（2-35）表示：當圖2-7所示微元體的傳熱達到平衡時，微元體升溫（或降溫）所需（或減少）熱量（方程左邊第1項）應與傳入（或傳出）微元體的熱量（方程左邊第2～4項）和微元體內部釋放（或吸收）的熱量（方程左邊第5項）相等。其中，熱源密度Q與具體工程問題有關，可以是固態相變產生的相變潛熱，也可以是壓力加工的能量轉換熱，或熱處理感應加熱產生的歐姆熱等。式（2-36）～式（2-38）分別稱為第一類（強制）邊界條件、第二類（自然）邊界條件和第三類（自然）邊界條件。