微分法在有限元分析中的應用

王朝振，劉銀濤，孫建鵬，周 鵬，孫文武，張家駒

（1.中鐵建大橋工程局集團第三工程有限公司，遼寧 沈陽110000；2.西安建筑科技大學土木工程學院，陜西 西安710055）

0 引 言

有限元法自提出以來，就得到了廣泛的關注。隨著時間的發展與計算機技術的進步，有關有限元法的各種理論以及應用得到了不斷補充與完善，同時與有限元法相關的各類通用和專用的CAE 軟件也得到了廣泛開發與應用[1-3]。如今有限元法已經形成了相當完善的計算體系，可以用來計算和模擬與結構相關的各類宏觀以及微觀力學介質問題。特別是有限元法作為一種成熟的結構分析方法[4-7]，被廣泛用來求解眾多工程問題的數值解[6-8]。

在結構分析中使用有限元法處理待離散結構時，其做法是采用基本單元將待離散結構劃分為若干子單元，例如采用三角形單元，那么就會在所離散的每一個三角形上分別得到一個近似的插值近似函數，這個函數被稱為形函數。形函數是針對單元內的插值而言的，一般指的是對單元內任意位置P，可以由單元所有節點（Pi）處的值進行插值而得到P 點的值，即Vp=N1×V1+ N2×V2+···+ Ni×Vi+···，其中Ni即為對應于Pi點的形函數，在數學上就是一種插值的權函數。

作為有限元法的基本特色之一，形函數的求解往往決定著有限元數值分析的準確性。Dasgupt[9]和Malsch 等[10]采用代數方法構造多邊形單元的Wachspress 插值函數，給出了采用三角形面積表達的插值形函數公式；王兆清等[11]給出了多邊形有限單元形函數的幾何構造方法，并分析總結了Wachspress 插值、Laplace 插值和平均值插值的構造方法和性質；李術才等[12]通過采用平均值插值方法，提出了一種求解微分方程邊值問題的多邊形有限元方法；王麗等[13]采用面積坐標方法和形函數譜方法構造四邊形薄板元，從而大大簡化了新單元的推導過程，而且導出的高階形函數非常簡潔；夏曉舟等[14]基于復合函數鏈式求導法則，獲得了三維自然單元法non-Sibson 插值形函數導數的顯式格式，比現有的Lasserre法更簡潔、直觀，大大縮減了計算過程；崔夢雷等[15]采用拉格朗日插值法，針對全局坐標下結構的位移形函數進行了求解，結果表明得到的形函數求解簡單，精度與常規逆變換相當。

雖然以上針對形函數的推導方法形式多樣，也解決了現有結構分析中存在的各種隱性問題，但是仍然存在著很大的主觀性，而且只能適用于性質相同的結構，造成公式的適應性很差。

采用解析式法能夠很好地解決現有的有限元形函數適應性太差的問題。解析式法的基本原理是依靠現有的有限元原理以及結構的平衡方程構造等量關系式，并將其與現有的邊界條件相結合，從而能夠針對不同的問題采用同樣的分析方式，解決形函數的重復推導問題，例如，傅向榮等[16]將有限元的離散法與解析法成果有機融合，給出了一個在彈性力學問題中構造獨立完備解析試函數的通用方法；許晶等[17]采用解析式法重新構造了Timoshenko 梁有限元的位移形函數方程；許晶等[18]在考慮桿件扭轉以及翹曲的情況下，采用解析式法，構造了全新的桿件扭曲位移形函數方程。

綜上，本文基于有限元原理以及結構的平衡方程，提出了一種求解結構單元形函數的方法，此方法推導出的桿單元以及板單元的形函數不僅與現有公式相同，而且理論上可以用于推導任何結構的形函數，不僅大大簡化了不同結構之間公式之間相互轉換的麻煩，同時將有限元形函數的推導過程透明化，為有限元的進一步推廣和拓展提供了一條行之有效的途徑。

1 基本原理

1.1 結構平衡微分方程

工程結構是一種靜定結構，根據冗余度的多少可分為靜定結構和超靜定結構。無論結構是處于靜止狀態還是運動狀態，結構在任意時刻都滿足平衡條件，或為靜力平衡或為運動平衡。結構構件靜力平衡的微分方程為：

式中：D 為構件的廣義剛度；V 為構件的廣義位移；n為階數；p 為單位長度上的荷載。

1.2 微分方程的解析解

有限元法中的形函數反映了結構單元內部的位移模式。因此，形函數可以考慮為僅體現單元自身可能出現的變形特征。根據有限元原理，結構的荷載及變形均集中在節點上，單元內部并無與變形相耦合的荷載。此時式(1)可表示為：

微分方程（2）的解析解為：

式中：ai（i=0，1，2，…，n-1）為待定常系數；fi（i=1，2，…，n-1）為已知關于x 的函數。

1.3 單元形函數

對于n 階微分方程所對應的單元，必有n 個節點自由度與之對應。將節點自由度及節點局部坐標代入式（3），形成關于ai（i=0，1，2，…，n-1）的方程組：

式中：A 為關于ai（i=0，1，2，…，n-1）的待定向量組；F 為關于fi（x）（i=0，1，2，…，n-1；x=0，l）的已知函數值；l 為單元長度；U 為單元節點廣義位移向量組。

根據式（4）就可以求得待定系數ai，將其代入式（3）并整理成關于節點廣義位移的矩陣，即可得到相應單元的形函數，記為N（x）。即：

式中：δe為單元節點廣義位移向量組。

2 實例

2.1 桿系單元

桿系單元根據受力特點可以分為兩類，一類是只有軸向荷載和軸向變形的桿，另一類是受軸力、剪力、彎矩以及扭矩等共同作用的梁。桿單元的靜力微分方程為：

式中：E 為桿件的楊氏彈性模量；A 為桿件的截面面積；u 為桿件沿x 方向的位移。

令荷載p=0，式（6）變為：

式（7）可進一步表示為：

微分方程（8）的解析解為：

式中：C1、C2分別為待定系數。

根據桿件的邊界條件，u（0）=u1，u（l）=u2。可以得到關于系數C1和C2的方程組：

將式（11）代入式（9），可以得到單元廣義位移與單元節點位移的關系式：

式（13）、式（14）均與文獻[8]相同。

平面梁單元的靜力微分方程為：

式中：I 為慣性矩。

同理，式（15）可進一步表示為：

微分方程（16）的解析解為：

式中：θ（x）為轉角，即撓度的1 階導數。

根據梁單元的邊界條件v（0）= v1，θ（0）= θ1，v（l）=v2，θ（l）=θ2，可 以 得 到 關 于 系 數C1、C2、C3、C4的方程組：

由式（19）可以得出:

將式（20）代入式（17），可以得到梁單元廣義位移與節點位移的關系式：

空間梁單元共包含4 個微分方程，分別為：

式中：Iy、Iz和J 分別為繞y 軸、z 軸的抗彎慣性矩和繞x 軸的扭轉慣性矩；G 為剪切模量。

空間梁單元形函數為桿單元形函數和梁單元形函數的組合：

式 中：ui、vi、wi、θyi、θzi、θxi（i=1，2）分 別 為 沿x、y、z 軸方向的位移及繞y、z、x 軸方向的轉角。

2.2 薄板單元

薄板單元的平衡微分方程為：（h 為薄板厚度，ν 為泊松比）；U（x，y）為板的撓度；q 為板受到的垂直板面的荷載集度。

對于矩形板單元，一般取4 個角點為節點，每個節點有2 個自由度。此時，薄板僅存在2 個坐標軸方向變形的相關項，則式（31）將簡化為：

式(32)可進一步表示為：

微分方程(33)的解析解為：

為了與節點位移相適應，式（34）可拓展為：

局部坐標系中節點的坐標分別記為：（ξ1，η1），（ξ2，η2），（ξ3，η3），（ξ4，η4）。將其代入式（36）可得：

式中：Ni=（1+ξξi）（1+ηηi）/4 。

3 結 語

本文根據結構單元的微分方程，結合有限元原理，提出了確定單元形函數的理論解。推導了桿系單元、板單元的形函數。推導出的形函數與已有的形函數完全一致，表明本文提出的方法是有效的。該方法基于結構的微分方程，理論上可用于任何單元形函數的求解。同時為建立特殊單元的形函數，應用有限元原理求解特殊問題提供了一種新的方法。