張軍鋒 李 杰 尹會娜 陳 淮 葉雨山

（1.鄭州大學(xué)土木工程學(xué)院，鄭州450001；2.中國建筑第七工程局有限公司，鄭州450001）

0 引 言

歐拉梁理論又稱為歐拉-伯努利（Euler-Bernoulli）梁理論或經(jīng)典梁理論，其忽略剪切變形對撓度的影響，但對高度較大的梁，其剪切變形對撓度的貢獻(xiàn)將不可忽略［1-2］。為提高歐拉梁理論的適用性，有限元中往往基于歐拉梁理論并計入剪切變形的影響，比如ANSYS中的Beam4∕44單元，前者針對等截面梁單元而后者針對變截面梁單元。

等截面梁單元剛度矩陣的推導(dǎo)可采用結(jié)構(gòu)力學(xué)形常數(shù)方法［2］和以形函數(shù)為基礎(chǔ)的最小勢能原理或虛功原理方法［1，3-4］，但對變截面梁單元剛度矩陣的推導(dǎo)卻鮮有涉及。文獻(xiàn)［5］針對兩結(jié)點(diǎn)梁單元給出了變截面梁單元的單剛矩陣，但未明確剛度矩陣的具體推導(dǎo)原理和過程，也沒有計入剪切變形，且所給單剛矩陣僅針對彎曲問題而不涉及伸縮和扭轉(zhuǎn)問題。文獻(xiàn)［6］采用傳遞矩陣法推導(dǎo)了任意變截面梁單元的剛度矩陣，但此方法不涉及形函數(shù)，與通用有限元理論不一致。文獻(xiàn)［7］則以形函數(shù)為基礎(chǔ)，借助于MATLAB的符號運(yùn)算功能建立了梁高成線性或拋物線形變化的變截面梁單元剛度矩陣，但其矩陣元素頗為復(fù)雜，不便使用。ANSYS的Beam44與Beam4的單剛矩陣形式相同，只是前者對截面參數(shù)包括面積A和慣性矩I采用了“平均值”，但并未給出截面參數(shù)“平均值”的推導(dǎo)原理。

基于有限元理論中的形函數(shù)和最小勢能原理，在等截面歐拉梁單元剛度矩陣推導(dǎo)過程的基礎(chǔ)上［8］，系統(tǒng)闡述了計入剪切變形的變截面歐拉梁單剛矩陣推導(dǎo)過程，并經(jīng)過對比分析明確了ANSYS中Beam44單元剛度矩陣的推導(dǎo)方法。為便于闡述，僅針對平面梁單元進(jìn)行推導(dǎo)。另外，由于ANSYS中的Beam44這一變截面梁單元要求截面在兩節(jié)點(diǎn)之間是線性梯度變化的，下文建立的變截面梁同樣符合這一要求。

1 基本方程

以XY平面梁單元為例，圖1給出了右手螺旋法則的坐標(biāo)系統(tǒng)，并針對伸縮、扭轉(zhuǎn)和彎曲問題給出了對應(yīng)的桿端荷載和分布荷載，所示荷載方向均與坐標(biāo)軸正方向一致，也即荷載的正方向。圖1（c）還給出了截面內(nèi)力正方向：軸力F以受拉為正，扭矩T的正方向與F方向一致，彎矩M以梁底（即-y方向）纖維受拉為正，剪力Q以使微段發(fā)生逆時針轉(zhuǎn)動為正［1，3-4］。

圖1 單元坐標(biāo)系以及參量正方向示意圖Fig.1 Positive direction regulations in element coordinate system

對于所示X-Y平面內(nèi)的兩節(jié)點(diǎn)單元，其位移函數(shù)φ也即變形包括伸縮、扭轉(zhuǎn)、彎曲和剪切位移四種變形且均可采用形函數(shù)表達(dá)為

式中：N和Φ分別為形函數(shù)向量和節(jié)點(diǎn)位移列向量；ψr表示節(jié)點(diǎn)i和j的位移，可分別為節(jié)點(diǎn)的軸向位移u或扭轉(zhuǎn)角θ、彎曲豎向位移ωb和剪切豎向位移ωs和截面轉(zhuǎn)角α；n為獨(dú)立節(jié)點(diǎn)位移數(shù)量，除彎曲豎向變形ωb的r=4外，其他三種變形均為r=2（表 1）；s為單元相對位置并取-1≤s≤1（圖 1，式（2）），xj和sj為單元j節(jié)點(diǎn)的絕對坐標(biāo)和相對坐標(biāo)；L為單元長度；xc為單元中點(diǎn)坐標(biāo)；Ni（s）和Ni（x）表示拉格朗日形函數(shù)或厄米特形函數(shù)，表1給出兩種形函數(shù)的形式和適用變形［1］。

需要指出，從結(jié)構(gòu)力學(xué)分析可知，變截面梁的形函數(shù)極為復(fù)雜，故在有限元中對變截面梁依然直接使用等截面梁的形函數(shù)以簡化分析。

表1 不同變形的形函數(shù)Table 1 Shape functions for different deformation functions

需要說明，剪切變形的計入使梁單元看似存在4種受力模式，但其純剪切和純彎曲受力實(shí)際上是耦合的［1］。為便于闡述，下文依然稱為4種受力模式，并對每種受力模式分別分析其單剛矩陣，再綜合為任意荷載模式下的單剛矩陣。為節(jié)約篇幅，基于前述位移函數(shù)和形函數(shù)以及最小勢能原理，以圖1所示坐標(biāo)系和參數(shù)正方向，表2直接給出各種受力模式的單剛矩陣計算方法，最小勢能原理以及所涉及的各種受力模式的幾何方程、本構(gòu)方程、截面內(nèi)力和平衡方程詳見文獻(xiàn)［8］；純剪切變形剛度矩陣的推導(dǎo)中，還根據(jù)能量等效原則引入了剪應(yīng)力不均勻修正系數(shù)k（以下簡稱剪切系數(shù)k）以便使用統(tǒng)一的單剛矩陣計算公式，詳見文獻(xiàn)［1］。

表2中，E和G為彈性模量和剪切模量，μ為泊松比；截面參數(shù)A、J、I分別為面積、抗扭慣性矩和抗彎慣性矩；由于僅針對平面梁單元，此處的I為繞z軸的抗彎慣性矩也即Iz（圖1、圖2）；y為偏離截面中性軸的坐標(biāo)值；S（y）和b（y）為截面剪應(yīng)力計算位置一側(cè)截面積對中性軸的靜矩和計算位置的截面寬度。對于變截面梁單元，截面參數(shù)A、J、I均隨s變化，故在剛度矩陣推導(dǎo)過程中應(yīng)始終將其放在積分號內(nèi)；對于截面剪切系數(shù)k，從其計算公式可知，其僅與截面形狀有關(guān)，如矩形、圓形和薄壁圓環(huán)截面分別為6∕5、10∕9和2，但部分變截面梁的k將隨s位置變化（圖2），故亦置于積分號內(nèi)。當(dāng)然，如果為等截面梁，截面參數(shù)A、J、I和剪切系數(shù)k均可置于積分號外，從而得對應(yīng)的剛度矩陣如表2所列。下文則以矩形、圓形和薄壁圓環(huán)和薄壁箱形截面為例，分析變截面梁的剛度矩陣，并與ANSYS所給剛度矩陣進(jìn)行對比。

2 剛度矩陣中的截面參數(shù)分析

2.1 ANSYS中Beam44單元截面參數(shù)取值

在ANSYS中，基于歐拉梁理論提供的等∕變截面梁單元分別為Beam4和Beam44，兩者的單剛矩陣形式相同（式（3）、式（4）），只是Beam44剛度矩陣中的截面參數(shù)A、I、J均采用了“平均值”作為等效截面參數(shù)（式（5）），其中的參數(shù)下標(biāo)i和j分別表示左右節(jié)點(diǎn)?；蛘哒fANSYS的Beam44本質(zhì)上也是一種等截面梁，只是以左右兩端截面參數(shù)的“平均值”作為截面參數(shù)取值。在ANSYS中，Beam44單元只能通過實(shí)常數(shù)命令直接定義i和j節(jié)點(diǎn)的截面參數(shù)A、I、J，并對k只能取一個值，而無法通過命令SECTYPE和SECDATA定義i和j節(jié)點(diǎn)的截面形狀和尺寸。這也說明，Beam44單元的單剛矩陣僅由兩端截面參數(shù)和式（3）-式（5）決定，與截面形狀沿桿件的變化規(guī)律無關(guān)，但從下文的理論分析可知，這樣所得的單剛矩陣對大部分截面只是一種近似。

表2 不同受力模式的單剛矩陣推導(dǎo)Table 2 Deviation of element stiffness matrixes

2.2 剪切系數(shù)k

對于截面剪切系數(shù)k，ANSYS中的Beam44單元和Beam4單元均只能輸入一個參數(shù)，即對變截面梁認(rèn)為其不隨位置改變。實(shí)際上，對于變截面梁，矩形、圓形和薄壁圓環(huán)截面的k本身并不隨位置變化，但箱形截面、工形和T形截面的k隨位置變化：圖2依文獻(xiàn)［10］對一箱形截面梁給出了不同截面尺寸的k值。但是，在單元兩端截面參數(shù)A、J和I比值0.5～2.0的約束下（詳見下文），k的變化較為有限，故下文亦對k取常數(shù)即取單元中截面的k值，從而可將表2中純剪切剛度矩陣推導(dǎo)中的k置于積分號以外，同時對Beam44單元亦有式（4）。

圖2 箱型截面的截面參數(shù)Fig.2 Section parameters of a box section

2.3 矩形截面參數(shù)分析

對截面剪切系數(shù)k取定值后，表2中各受力模式剛度矩陣的積分式內(nèi)僅有截面參數(shù)和形函數(shù)導(dǎo)數(shù)向量，但因截面形狀多樣，且截面尺寸隨桿件長度變化形式多樣，直接對其進(jìn)行積分仍非常復(fù)雜，且結(jié)果與ANSYS中Beam44單元所得單剛矩陣中的部分元素仍有較大差異。下文首先以最簡單的矩形截面為例進(jìn)行分析，并且此時剪切系數(shù)k本身不隨s變化。

取矩形截面變截面梁i、j兩端截面尺寸分別為bj×hj=c1bi×c2hi（式（6）），其中c1和c2分別為截面寬度和高度的變化系數(shù)（以下簡稱截面變化系數(shù)），假定截面尺寸隨s線性變化（這也與ANSYS中對Beam44單元的假定一致），從而可得各截面參數(shù)如式（7）所示。

對于式（7）需要說明：對于本文的平面梁單元，其抗彎慣性矩I也即Iz，但同時也給出了Iy以便對比；矩形截面抗扭慣性矩J并沒有通用的計算表達(dá)式，往往采用J=χb4計算，其中系數(shù)χ表達(dá)式極為復(fù)雜，故實(shí)際計算時往往根據(jù)高寬比h∕b直接對χ取值［11］。由于難以給出χ隨s的變化表達(dá)式，故式（7）對J的表達(dá)式僅適用于截面變化系數(shù)c1=c2，也即截面高寬比h∕b不隨s變化的情況。

將式（7）所示A、J、I帶入表2進(jìn)行積分計算即可得單一受力模式下的單剛矩陣。從表2所列剛度矩陣推導(dǎo)過程可知，對于除純彎曲以外的其他3種受力模式，積分計算僅針對截面面積A和抗扭慣矩J。對于軸向伸縮和純剪切受力，積分計算僅針對截面面積A且可得式（8）。對比式（5）可知，只有當(dāng)截面變化系數(shù)c1=c2時也即截面寬度和高度按同一規(guī)律變化時，ANSYS所給等效截面面積才與理論值一致，在c1≠c2時ANSYS所給等效截面面積只是理論值的一種近似。對于扭轉(zhuǎn)受力，由于難以給出c1≠c2時J的計算表達(dá)式，但c1=c2時則有式（9）成立，同樣說明只有當(dāng)c1=c2時ANSYS所給等效截面抗扭慣性矩（式（5））才與理論值一致。

純彎受力的積分表達(dá)式更為復(fù)雜，其抗彎慣性矩I和形函數(shù)二階導(dǎo)數(shù)N''均有積分變量s，積分所得表達(dá)式過于復(fù)雜，且通過計算對比發(fā)現(xiàn)即使c1=c2時所得剛度矩陣?yán)碚撝狄才cANSYS所得不一致。為明確ANSYS中Beam44單元單剛矩陣的計算方法，通過嘗試發(fā)現(xiàn)，如果將分進(jìn)行計算（式（10）），也即對幾何矩陣和截面參數(shù)分別積分，則在c1=c2時所得剛度矩陣將與ANSYS結(jié)果一致（式（11）-式與ANSYS所給等效截面抗彎慣性矩（式（5））一致且所得矩陣與純彎受力模式的剛度矩陣形式一致；兩者聯(lián)合并考慮系數(shù)8E∕L3后即得ANSYS變截面梁純彎剛度矩陣。這實(shí)際也明確了ANSYS中對Beam44單元純彎模式剛度矩陣的計算方法，即將″ds進(jìn)行計算，但對矩形截面也僅當(dāng)c1=c2時所得單剛矩陣與ANSYS所給單剛矩陣一致。

明確了矩形截面單一受力模式下ANSYS中單剛矩陣后，即可將其綜合為式（3）所示任意受力模式下的單剛矩陣。而從上述推導(dǎo)可知，對于矩形截面梁，各受力模式下ANSYS中Beam44單元所得單剛矩陣僅在截面寬度和高度按同一規(guī)律變化也即截面變化系數(shù)c1=c2時與理論值一致，在c1≠c2時僅為一種近似。實(shí)際上，ANSYS中Beam44單元對矩形截面并無c1=c2的限制，但ANSYS要求Beam44 單元兩端截面參數(shù)的比值 Ai∕Aj、Ji∕Jj、Ii∕Ij均應(yīng)在0.5～2.0之間，否則會提示警告信息，如果上述比值超出0.1～10.0的范圍，則會提示錯誤信息。這可能也是為了保證當(dāng)c1≠c2時ANSYS所得等效截面參數(shù)與理論結(jié)果的差異不至過大而提出的要求：兩端截面參數(shù)比值越接近，兩者的差異也就越小。

2.4 其他常用截面參數(shù)分析

除矩形截面外，圓形、圓環(huán)和箱形截面也是常用的截面形式。對于圓形截面，因其截面尺寸僅有1個參數(shù)即半徑R，截面變化系數(shù)同樣僅有1個即c（式（14））：因圓形截面的Iz=Iy=0.5J=I故僅給出I的表達(dá)式，并且對比式（7）可知，只需將式（7）的c1和c2取為c即可的圓形截面A、J、I的表達(dá)式，且圓截面的k值亦不隨半徑變化。對其分析結(jié)果表明，在采用式（10）后，其理論推導(dǎo)所得剛度矩陣與ANSYS采用等效截面參數(shù)的表達(dá)式完全一致。

對于薄壁圓環(huán)截面，其截面尺寸有2個參數(shù)即外徑R與壁厚t，相應(yīng)的截面變化系數(shù)同樣為2個，即 c1和 c2。故其截面尺寸和 A、J、I的表達(dá)式（式（15））與式（7）一致，并且由于是圓環(huán)截面，其J值表達(dá)式不再有c1=c2的限制。對于箱形截面，其截面尺寸參數(shù)更多，理論上有6個，但一般頂?shù)装搴穸群妥笥腋拱搴穸确謩e相等故有4個參數(shù)（圖2），相應(yīng)的截面變化系數(shù)同樣為4個。對圓環(huán)和箱形其分析結(jié)果表明，在假定k值不隨桿件位置變化且使用式（10）后，其理論推導(dǎo)所得剛度矩陣同樣在c1=c2=c和c1=c2=c3=c4=c時與ANSYS采用等效截面參數(shù)的表達(dá)式一致。

需要說明，式（15）對A、J、I的表達(dá)式僅適用于薄壁圓環(huán)，并非精確表達(dá)式（式（16））。如果采用式（16）的精確表達(dá)式，則在c1≠c2時對A和I積分不再有式（8）和式（11）的結(jié)果，而式（15）的近似表達(dá)式在c1≠c2時同樣有式（8）和式（11）成立，但精確表達(dá)式在c1=c2=c時同樣可得式（5）。

綜上可推知，ANSYS對Beam44變截面梁單元單剛矩陣推導(dǎo)時：對純彎模式采用進(jìn)行計算，同時對截面剪切系數(shù)k取定值；所給等效截面參數(shù)A、J和I表達(dá)式源于截面各個尺寸的截面變化系數(shù)相同時的情況，或者梁單元兩端截面形狀為相似形且截面尺寸隨桿件線性變化的情況，其他情況ANSYS所得等效截面參數(shù)與理論值會有一定偏差。

3 結(jié) 論

基于有限元理論系統(tǒng)闡述了計入剪切變形的變截面歐拉梁單剛矩陣推導(dǎo)過程，并以矩形、圓形、圓環(huán)和箱形截面梁為例與ANSYS中Beam44單元剛度矩陣進(jìn)行了對比，明確了Beam44單元剛度矩陣的推導(dǎo)方法。研究發(fā)現(xiàn)：

（1）變截面梁因截面形狀多樣且截面尺寸隨桿件長度變化形式多樣，即使對于線性梯度變化的變截面梁，對部分截面也難以給出截面參數(shù)，如抗扭剛度J和剪切系數(shù)k隨桿件位置的表達(dá)式，且變截面梁使用等截面梁的形函數(shù)本身也是一種近似，故難以對其進(jìn)行積分計算獲得普適性的剛度矩陣?yán)碚摫磉_(dá)式。

（2）ANSYS對Beam44變截面梁單元單剛矩陣推導(dǎo)時，純彎模式的積分計算中對幾何矩陣和截面參數(shù)分別積分，即以代替進(jìn)行計算，這樣既可簡化計算，又可使各種受力模式的計算方式一致。

（3）Beam44單元對截面剪切系數(shù)k取定值以簡化計算，故在使用中宜取單元中截面的剪切系數(shù)k或左右兩端截面的均值以減少誤差。

（4）ANSYS所得Beam44單元的單剛矩陣形式與等截面梁完全一致，只是對截面面積A、抗扭慣矩J和抗彎慣矩I采用了等效值。ANSYS所給等效截面參數(shù)A、J和I表達(dá)式，源于梁單元兩端截面形狀為相似形且截面尺寸隨桿件線性梯度變化的情況，其他情況的等效截面參數(shù)與理論值會有一定偏差。ANSYS對Beam44單元僅要求截面在兩節(jié)點(diǎn)之間按線性梯度變化，并無要求兩端截面形狀為相似形，這實(shí)際也是將簡單模式所得表達(dá)式直接推廣應(yīng)用至各種情況。但ANSYS對Beam44單元要求兩端截面參數(shù)的比率盡可能接近1.0，這可能也是為了控制等效截面參數(shù)與理論值的偏差，故對變截面梁在建模時應(yīng)盡量減小單元長度以提高計算精度。