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

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

0 引 言

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

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

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

雖然以上針對形函數(shù)的推導(dǎo)方法形式多樣，也解決了現(xiàn)有結(jié)構(gòu)分析中存在的各種隱性問題，但是仍然存在著很大的主觀性，而且只能適用于性質(zhì)相同的結(jié)構(gòu)，造成公式的適應(yīng)性很差。

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

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

1 基本原理

1.1 結(jié)構(gòu)平衡微分方程

工程結(jié)構(gòu)是一種靜定結(jié)構(gòu)，根據(jù)冗余度的多少可分為靜定結(jié)構(gòu)和超靜定結(jié)構(gòu)。無論結(jié)構(gòu)是處于靜止?fàn)顟B(tài)還是運(yùn)動(dòng)狀態(tài)，結(jié)構(gòu)在任意時(shí)刻都滿足平衡條件，或?yàn)殪o力平衡或?yàn)檫\(yùn)動(dòng)平衡。結(jié)構(gòu)構(gòu)件靜力平衡的微分方程為：

式中：D 為構(gòu)件的廣義剛度；V 為構(gòu)件的廣義位移；n為階數(shù)；p 為單位長度上的荷載。

1.2 微分方程的解析解

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

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

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

1.3 單元形函數(shù)

對于n 階微分方程所對應(yīng)的單元，必有n 個(gè)節(jié)點(diǎn)自由度與之對應(yīng)。將節(jié)點(diǎn)自由度及節(jié)點(diǎn)局部坐標(biāo)代入式（3），形成關(guān)于ai（i=0，1，2，…，n-1）的方程組：

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

根據(jù)式（4）就可以求得待定系數(shù)ai，將其代入式（3）并整理成關(guān)于節(jié)點(diǎn)廣義位移的矩陣，即可得到相應(yīng)單元的形函數(shù)，記為N（x）。即：

式中：δe為單元節(jié)點(diǎn)廣義位移向量組。

2 實(shí)例

2.1 桿系單元

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

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

令荷載p=0，式（6）變?yōu)椋?/p>

式（7）可進(jìn)一步表示為：

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

式中：C1、C2分別為待定系數(shù)。

根據(jù)桿件的邊界條件，u（0）=u1，u（l）=u2?？梢缘玫疥P(guān)于系數(shù)C1和C2的方程組：

將式（11）代入式（9），可以得到單元廣義位移與單元節(jié)點(diǎn)位移的關(guān)系式：

式（13）、式（14）均與文獻(xiàn)[8]相同。

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

式中：I 為慣性矩。

同理，式（15）可進(jìn)一步表示為：

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

式中：θ（x）為轉(zhuǎn)角，即撓度的1 階導(dǎo)數(shù)。

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

由式（19）可以得出:

將式（20）代入式（17），可以得到梁單元廣義位移與節(jié)點(diǎn)位移的關(guān)系式：

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

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

空間梁單元形函數(shù)為桿單元形函數(shù)和梁單元形函數(shù)的組合：

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

2.2 薄板單元

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

對于矩形板單元，一般取4 個(gè)角點(diǎn)為節(jié)點(diǎn)，每個(gè)節(jié)點(diǎn)有2 個(gè)自由度。此時(shí)，薄板僅存在2 個(gè)坐標(biāo)軸方向變形的相關(guān)項(xiàng)，則式（31）將簡化為：

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

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

為了與節(jié)點(diǎn)位移相適應(yīng)，式（34）可拓展為：

局部坐標(biāo)系中節(jié)點(diǎn)的坐標(biāo)分別記為：（ξ1，η1），（ξ2，η2），（ξ3，η3），（ξ4，η4）。將其代入式（36）可得：

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

3 結(jié) 語

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