考慮剪切變形的變截面歐拉梁?jiǎn)卧獎(jiǎng)偠染仃?/h1>
2020-06-17 09:12:40張軍鋒尹會(huì)娜葉雨山
結(jié)構(gòu)工程師 2020年2期
張軍鋒 李 杰 尹會(huì)娜 陳 淮 葉雨山
（1.鄭州大學(xué)土木工程學(xué)院，鄭州450001；2.中國(guó)建筑第七工程局有限公司，鄭州450001）
0 引 言
歐拉梁理論又稱為歐拉-伯努利（Euler-Bernoulli）梁理論或經(jīng)典梁理論，其忽略剪切變形對(duì)撓度的影響，但對(duì)高度較大的梁，其剪切變形對(duì)撓度的貢獻(xiàn)將不可忽略［1-2］。為提高歐拉梁理論的適用性，有限元中往往基于歐拉梁理論并計(jì)入剪切變形的影響，比如ANSYS中的Beam4∕44單元，前者針對(duì)等截面梁?jiǎn)卧笳哚槍?duì)變截面梁?jiǎn)卧?/p>
等截面梁?jiǎn)卧獎(jiǎng)偠染仃嚨耐茖?dǎo)可采用結(jié)構(gòu)力學(xué)形常數(shù)方法［2］和以形函數(shù)為基礎(chǔ)的最小勢(shì)能原理或虛功原理方法［1，3-4］，但對(duì)變截面梁?jiǎn)卧獎(jiǎng)偠染仃嚨耐茖?dǎo)卻鮮有涉及。文獻(xiàn)［5］針對(duì)兩結(jié)點(diǎn)梁?jiǎn)卧o出了變截面梁?jiǎn)卧膯蝿偩仃?，但未明確剛度矩陣的具體推導(dǎo)原理和過(guò)程，也沒(méi)有計(jì)入剪切變形，且所給單剛矩陣僅針對(duì)彎曲問(wèn)題而不涉及伸縮和扭轉(zhuǎn)問(wèn)題。文獻(xiàn)［6］采用傳遞矩陣法推導(dǎo)了任意變截面梁?jiǎn)卧膭偠染仃?，但此方法不涉及形函?shù)，與通用有限元理論不一致。文獻(xiàn)［7］則以形函數(shù)為基礎(chǔ)，借助于MATLAB的符號(hào)運(yùn)算功能建立了梁高成線性或拋物線形變化的變截面梁?jiǎn)卧獎(jiǎng)偠染仃?，但其矩陣元素頗為復(fù)雜，不便使用。ANSYS的Beam44與Beam4的單剛矩陣形式相同，只是前者對(duì)截面參數(shù)包括面積A和慣性矩I采用了“平均值”，但并未給出截面參數(shù)“平均值”的推導(dǎo)原理。
基于有限元理論中的形函數(shù)和最小勢(shì)能原理，在等截面歐拉梁?jiǎn)卧獎(jiǎng)偠染仃囃茖?dǎo)過(guò)程的基礎(chǔ)上［8］，系統(tǒng)闡述了計(jì)入剪切變形的變截面歐拉梁?jiǎn)蝿偩仃囃茖?dǎo)過(guò)程，并經(jīng)過(guò)對(duì)比分析明確了ANSYS中Beam44單元?jiǎng)偠染仃嚨耐茖?dǎo)方法。為便于闡述，僅針對(duì)平面梁?jiǎn)卧M(jìn)行推導(dǎo)。另外，由于ANSYS中的Beam44這一變截面梁?jiǎn)卧蠼孛嬖趦晒?jié)點(diǎn)之間是線性梯度變化的，下文建立的變截面梁同樣符合這一要求。
1 基本方程
以XY平面梁?jiǎn)卧獮槔?，圖1給出了右手螺旋法則的坐標(biāo)系統(tǒng)，并針對(duì)伸縮、扭轉(zhuǎn)和彎曲問(wèn)題給出了對(duì)應(yīng)的桿端荷載和分布荷載，所示荷載方向均與坐標(biāo)軸正方向一致，也即荷載的正方向。圖1（c）還給出了截面內(nèi)力正方向：軸力F以受拉為正，扭矩T的正方向與F方向一致，彎矩M以梁底（即-y方向）纖維受拉為正，剪力Q以使微段發(fā)生逆時(shí)針轉(zhuǎn)動(dòng)為正［1，3-4］。
圖1 單元坐標(biāo)系以及參量正方向示意圖Fig.1 Positive direction regulations in element coordinate system
對(duì)于所示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為單元相對(duì)位置并取-1≤s≤1（圖 1，式（2）），xj和sj為單元j節(jié)點(diǎn)的絕對(duì)坐標(biāo)和相對(duì)坐標(biāo)；L為單元長(zhǎng)度；xc為單元中點(diǎn)坐標(biāo)；Ni（s）和Ni（x）表示拉格朗日形函數(shù)或厄米特形函數(shù)，表1給出兩種形函數(shù)的形式和適用變形［1］。
需要指出，從結(jié)構(gòu)力學(xué)分析可知，變截面梁的形函數(shù)極為復(fù)雜，故在有限元中對(duì)變截面梁依然直接使用等截面梁的形函數(shù)以簡(jiǎn)化分析。
表1 不同變形的形函數(shù)Table 1 Shape functions for different deformation functions
需要說(shuō)明，剪切變形的計(jì)入使梁?jiǎn)卧此拼嬖?種受力模式，但其純剪切和純彎曲受力實(shí)際上是耦合的［1］。為便于闡述，下文依然稱為4種受力模式，并對(duì)每種受力模式分別分析其單剛矩陣，再綜合為任意荷載模式下的單剛矩陣。為節(jié)約篇幅，基于前述位移函數(shù)和形函數(shù)以及最小勢(shì)能原理，以圖1所示坐標(biāo)系和參數(shù)正方向，表2直接給出各種受力模式的單剛矩陣計(jì)算方法，最小勢(shì)能原理以及所涉及的各種受力模式的幾何方程、本構(gòu)方程、截面內(nèi)力和平衡方程詳見(jiàn)文獻(xiàn)［8］；純剪切變形剛度矩陣的推導(dǎo)中，還根據(jù)能量等效原則引入了剪應(yīng)力不均勻修正系數(shù)k（以下簡(jiǎn)稱剪切系數(shù)k）以便使用統(tǒng)一的單剛矩陣計(jì)算公式，詳見(jiàn)文獻(xiàn)［1］。
表2中，E和G為彈性模量和剪切模量，μ為泊松比；截面參數(shù)A、J、I分別為面積、抗扭慣性矩和抗彎慣性矩；由于僅針對(duì)平面梁?jiǎn)卧?，此處的I為繞z軸的抗彎慣性矩也即Iz（圖1、圖2）；y為偏離截面中性軸的坐標(biāo)值；S（y）和b（y）為截面剪應(yīng)力計(jì)算位置一側(cè)截面積對(duì)中性軸的靜矩和計(jì)算位置的截面寬度。對(duì)于變截面梁?jiǎn)卧?，截面參?shù)A、J、I均隨s變化，故在剛度矩陣推導(dǎo)過(guò)程中應(yīng)始終將其放在積分號(hào)內(nèi)；對(duì)于截面剪切系數(shù)k，從其計(jì)算公式可知，其僅與截面形狀有關(guān)，如矩形、圓形和薄壁圓環(huán)截面分別為6∕5、10∕9和2，但部分變截面梁的k將隨s位置變化（圖2），故亦置于積分號(hào)內(nèi)。當(dāng)然，如果為等截面梁，截面參數(shù)A、J、I和剪切系數(shù)k均可置于積分號(hào)外，從而得對(duì)應(yīng)的剛度矩陣如表2所列。下文則以矩形、圓形和薄壁圓環(huán)和薄壁箱形截面為例，分析變截面梁的剛度矩陣，并與ANSYS所給剛度矩陣進(jìn)行對(duì)比。
2 剛度矩陣中的截面參數(shù)分析
2.1 ANSYS中Beam44單元截面參數(shù)取值
在ANSYS中，基于歐拉梁理論提供的等∕變截面梁?jiǎn)卧謩e為Beam4和Beam44，兩者的單剛矩陣形式相同（式（3）、式（4）），只是Beam44剛度矩陣中的截面參數(shù)A、I、J均采用了“平均值”作為等效截面參數(shù)（式（5）），其中的參數(shù)下標(biāo)i和j分別表示左右節(jié)點(diǎn)。或者說(shuō)ANSYS的Beam44本質(zhì)上也是一種等截面梁，只是以左右兩端截面參數(shù)的“平均值”作為截面參數(shù)取值。在ANSYS中，Beam44單元只能通過(guò)實(shí)常數(shù)命令直接定義i和j節(jié)點(diǎn)的截面參數(shù)A、I、J，并對(duì)k只能取一個(gè)值，而無(wú)法通過(guò)命令SECTYPE和SECDATA定義i和j節(jié)點(diǎn)的截面形狀和尺寸。這也說(shuō)明，Beam44單元的單剛矩陣僅由兩端截面參數(shù)和式（3）-式（5）決定，與截面形狀沿桿件的變化規(guī)律無(wú)關(guān)，但從下文的理論分析可知，這樣所得的單剛矩陣對(duì)大部分截面只是一種近似。
表2 不同受力模式的單剛矩陣推導(dǎo)Table 2 Deviation of element stiffness matrixes
2.2 剪切系數(shù)k
對(duì)于截面剪切系數(shù)k，ANSYS中的Beam44單元和Beam4單元均只能輸入一個(gè)參數(shù)，即對(duì)變截面梁認(rèn)為其不隨位置改變。實(shí)際上，對(duì)于變截面梁，矩形、圓形和薄壁圓環(huán)截面的k本身并不隨位置變化，但箱形截面、工形和T形截面的k隨位置變化：圖2依文獻(xiàn)［10］對(duì)一箱形截面梁給出了不同截面尺寸的k值。但是，在單元兩端截面參數(shù)A、J和I比值0.5～2.0的約束下（詳見(jiàn)下文），k的變化較為有限，故下文亦對(duì)k取常數(shù)即取單元中截面的k值，從而可將表2中純剪切剛度矩陣推導(dǎo)中的k置于積分號(hào)以外，同時(shí)對(duì)Beam44單元亦有式（4）。
圖2 箱型截面的截面參數(shù)Fig.2 Section parameters of a box section
2.3 矩形截面參數(shù)分析
對(duì)截面剪切系數(shù)k取定值后，表2中各受力模式剛度矩陣的積分式內(nèi)僅有截面參數(shù)和形函數(shù)導(dǎo)數(shù)向量，但因截面形狀多樣，且截面尺寸隨桿件長(zhǎng)度變化形式多樣，直接對(duì)其進(jìn)行積分仍非常復(fù)雜，且結(jié)果與ANSYS中Beam44單元所得單剛矩陣中的部分元素仍有較大差異。下文首先以最簡(jiǎn)單的矩形截面為例進(jìn)行分析，并且此時(shí)剪切系數(shù)k本身不隨s變化。
取矩形截面變截面梁i、j兩端截面尺寸分別為bj×hj=c1bi×c2hi（式（6）），其中c1和c2分別為截面寬度和高度的變化系數(shù)（以下簡(jiǎn)稱截面變化系數(shù)），假定截面尺寸隨s線性變化（這也與ANSYS中對(duì)Beam44單元的假定一致），從而可得各截面參數(shù)如式（7）所示。
對(duì)于式（7）需要說(shuō)明：對(duì)于本文的平面梁?jiǎn)卧淇箯潙T性矩I也即Iz，但同時(shí)也給出了Iy以便對(duì)比；矩形截面抗扭慣性矩J并沒(méi)有通用的計(jì)算表達(dá)式，往往采用J=χb4計(jì)算，其中系數(shù)χ表達(dá)式極為復(fù)雜，故實(shí)際計(jì)算時(shí)往往根據(jù)高寬比h∕b直接對(duì)χ取值［11］。由于難以給出χ隨s的變化表達(dá)式，故式（7）對(duì)J的表達(dá)式僅適用于截面變化系數(shù)c1=c2，也即截面高寬比h∕b不隨s變化的情況。
將式（7）所示A、J、I帶入表2進(jìn)行積分計(jì)算即可得單一受力模式下的單剛矩陣。從表2所列剛度矩陣推導(dǎo)過(guò)程可知，對(duì)于除純彎曲以外的其他3種受力模式，積分計(jì)算僅針對(duì)截面面積A和抗扭慣矩J。對(duì)于軸向伸縮和純剪切受力，積分計(jì)算僅針對(duì)截面面積A且可得式（8）。對(duì)比式（5）可知，只有當(dāng)截面變化系數(shù)c1=c2時(shí)也即截面寬度和高度按同一規(guī)律變化時(shí)，ANSYS所給等效截面面積才與理論值一致，在c1≠c2時(shí)ANSYS所給等效截面面積只是理論值的一種近似。對(duì)于扭轉(zhuǎn)受力，由于難以給出c1≠c2時(shí)J的計(jì)算表達(dá)式，但c1=c2時(shí)則有式（9）成立，同樣說(shuō)明只有當(dāng)c1=c2時(shí)ANSYS所給等效截面抗扭慣性矩（式（5））才與理論值一致。
純彎受力的積分表達(dá)式更為復(fù)雜，其抗彎慣性矩I和形函數(shù)二階導(dǎo)數(shù)N''均有積分變量s，積分所得表達(dá)式過(guò)于復(fù)雜，且通過(guò)計(jì)算對(duì)比發(fā)現(xiàn)即使c1=c2時(shí)所得剛度矩陣?yán)碚撝狄才cANSYS所得不一致。為明確ANSYS中Beam44單元單剛矩陣的計(jì)算方法，通過(guò)嘗試發(fā)現(xiàn)，如果將分進(jìn)行計(jì)算（式（10）），也即對(duì)幾何矩陣和截面參數(shù)分別積分，則在c1=c2時(shí)所得剛度矩陣將與ANSYS結(jié)果一致（式（11）-式與ANSYS所給等效截面抗彎慣性矩（式（5））一致且所得矩陣與純彎受力模式的剛度矩陣形式一致；兩者聯(lián)合并考慮系數(shù)8E∕L3后即得ANSYS變截面梁純彎剛度矩陣。這實(shí)際也明確了ANSYS中對(duì)Beam44單元純彎模式剛度矩陣的計(jì)算方法，即將″ds進(jìn)行計(jì)算，但對(duì)矩形截面也僅當(dāng)c1=c2時(shí)所得單剛矩陣與ANSYS所給單剛矩陣一致。
明確了矩形截面單一受力模式下ANSYS中單剛矩陣后，即可將其綜合為式（3）所示任意受力模式下的單剛矩陣。而從上述推導(dǎo)可知，對(duì)于矩形截面梁，各受力模式下ANSYS中Beam44單元所得單剛矩陣僅在截面寬度和高度按同一規(guī)律變化也即截面變化系數(shù)c1=c2時(shí)與理論值一致，在c1≠c2時(shí)僅為一種近似。實(shí)際上，ANSYS中Beam44單元對(duì)矩形截面并無(wú)c1=c2的限制，但ANSYS要求Beam44 單元兩端截面參數(shù)的比值 Ai∕Aj、Ji∕Jj、Ii∕Ij均應(yīng)在0.5～2.0之間，否則會(huì)提示警告信息，如果上述比值超出0.1～10.0的范圍，則會(huì)提示錯(cuò)誤信息。這可能也是為了保證當(dāng)c1≠c2時(shí)ANSYS所得等效截面參數(shù)與理論結(jié)果的差異不至過(guò)大而提出的要求：兩端截面參數(shù)比值越接近，兩者的差異也就越小。
2.4 其他常用截面參數(shù)分析
除矩形截面外，圓形、圓環(huán)和箱形截面也是常用的截面形式。對(duì)于圓形截面，因其截面尺寸僅有1個(gè)參數(shù)即半徑R，截面變化系數(shù)同樣僅有1個(gè)即c（式（14））：因圓形截面的Iz=Iy=0.5J=I故僅給出I的表達(dá)式，并且對(duì)比式（7）可知，只需將式（7）的c1和c2取為c即可的圓形截面A、J、I的表達(dá)式，且圓截面的k值亦不隨半徑變化。對(duì)其分析結(jié)果表明，在采用式（10）后，其理論推導(dǎo)所得剛度矩陣與ANSYS采用等效截面參數(shù)的表達(dá)式完全一致。
對(duì)于薄壁圓環(huán)截面，其截面尺寸有2個(gè)參數(shù)即外徑R與壁厚t，相應(yīng)的截面變化系數(shù)同樣為2個(gè)，即 c1和 c2。故其截面尺寸和 A、J、I的表達(dá)式（式（15））與式（7）一致，并且由于是圓環(huán)截面，其J值表達(dá)式不再有c1=c2的限制。對(duì)于箱形截面，其截面尺寸參數(shù)更多，理論上有6個(gè)，但一般頂?shù)装搴穸群妥笥腋拱搴穸确謩e相等故有4個(gè)參數(shù)（圖2），相應(yīng)的截面變化系數(shù)同樣為4個(gè)。對(duì)圓環(huán)和箱形其分析結(jié)果表明，在假定k值不隨桿件位置變化且使用式（10）后，其理論推導(dǎo)所得剛度矩陣同樣在c1=c2=c和c1=c2=c3=c4=c時(shí)與ANSYS采用等效截面參數(shù)的表達(dá)式一致。
需要說(shuō)明，式（15）對(duì)A、J、I的表達(dá)式僅適用于薄壁圓環(huán)，并非精確表達(dá)式（式（16））。如果采用式（16）的精確表達(dá)式，則在c1≠c2時(shí)對(duì)A和I積分不再有式（8）和式（11）的結(jié)果，而式（15）的近似表達(dá)式在c1≠c2時(shí)同樣有式（8）和式（11）成立，但精確表達(dá)式在c1=c2=c時(shí)同樣可得式（5）。
綜上可推知，ANSYS對(duì)Beam44變截面梁?jiǎn)卧獑蝿偩仃囃茖?dǎo)時(shí)：對(duì)純彎模式采用進(jìn)行計(jì)算，同時(shí)對(duì)截面剪切系數(shù)k取定值；所給等效截面參數(shù)A、J和I表達(dá)式源于截面各個(gè)尺寸的截面變化系數(shù)相同時(shí)的情況，或者梁?jiǎn)卧獌啥私孛嫘螤顬橄嗨菩吻医孛娉叽珉S桿件線性變化的情況，其他情況ANSYS所得等效截面參數(shù)與理論值會(huì)有一定偏差。
3 結(jié) 論
基于有限元理論系統(tǒng)闡述了計(jì)入剪切變形的變截面歐拉梁?jiǎn)蝿偩仃囃茖?dǎo)過(guò)程，并以矩形、圓形、圓環(huán)和箱形截面梁為例與ANSYS中Beam44單元?jiǎng)偠染仃囘M(jìn)行了對(duì)比，明確了Beam44單元?jiǎng)偠染仃嚨耐茖?dǎo)方法。研究發(fā)現(xiàn)：
（1）變截面梁因截面形狀多樣且截面尺寸隨桿件長(zhǎng)度變化形式多樣，即使對(duì)于線性梯度變化的變截面梁，對(duì)部分截面也難以給出截面參數(shù)，如抗扭剛度J和剪切系數(shù)k隨桿件位置的表達(dá)式，且變截面梁使用等截面梁的形函數(shù)本身也是一種近似，故難以對(duì)其進(jìn)行積分計(jì)算獲得普適性的剛度矩陣?yán)碚摫磉_(dá)式。
（2）ANSYS對(duì)Beam44變截面梁?jiǎn)卧獑蝿偩仃囃茖?dǎo)時(shí)，純彎模式的積分計(jì)算中對(duì)幾何矩陣和截面參數(shù)分別積分，即以代替進(jìn)行計(jì)算，這樣既可簡(jiǎn)化計(jì)算，又可使各種受力模式的計(jì)算方式一致。
（3）Beam44單元對(duì)截面剪切系數(shù)k取定值以簡(jiǎn)化計(jì)算，故在使用中宜取單元中截面的剪切系數(shù)k或左右兩端截面的均值以減少誤差。
（4）ANSYS所得Beam44單元的單剛矩陣形式與等截面梁完全一致，只是對(duì)截面面積A、抗扭慣矩J和抗彎慣矩I采用了等效值。ANSYS所給等效截面參數(shù)A、J和I表達(dá)式，源于梁?jiǎn)卧獌啥私孛嫘螤顬橄嗨菩吻医孛娉叽珉S桿件線性梯度變化的情況，其他情況的等效截面參數(shù)與理論值會(huì)有一定偏差。ANSYS對(duì)Beam44單元僅要求截面在兩節(jié)點(diǎn)之間按線性梯度變化，并無(wú)要求兩端截面形狀為相似形，這實(shí)際也是將簡(jiǎn)單模式所得表達(dá)式直接推廣應(yīng)用至各種情況。但ANSYS對(duì)Beam44單元要求兩端截面參數(shù)的比率盡可能接近1.0，這可能也是為了控制等效截面參數(shù)與理論值的偏差，故對(duì)變截面梁在建模時(shí)應(yīng)盡量減小單元長(zhǎng)度以提高計(jì)算精度。
張軍鋒 李 杰 尹會(huì)娜 陳 淮 葉雨山
（1.鄭州大學(xué)土木工程學(xué)院，鄭州450001；2.中國(guó)建筑第七工程局有限公司，鄭州450001）
0 引 言
歐拉梁理論又稱為歐拉-伯努利（Euler-Bernoulli）梁理論或經(jīng)典梁理論，其忽略剪切變形對(duì)撓度的影響，但對(duì)高度較大的梁，其剪切變形對(duì)撓度的貢獻(xiàn)將不可忽略［1-2］。為提高歐拉梁理論的適用性，有限元中往往基于歐拉梁理論并計(jì)入剪切變形的影響，比如ANSYS中的Beam4∕44單元，前者針對(duì)等截面梁?jiǎn)卧笳哚槍?duì)變截面梁?jiǎn)卧?/p>
等截面梁?jiǎn)卧獎(jiǎng)偠染仃嚨耐茖?dǎo)可采用結(jié)構(gòu)力學(xué)形常數(shù)方法［2］和以形函數(shù)為基礎(chǔ)的最小勢(shì)能原理或虛功原理方法［1，3-4］，但對(duì)變截面梁?jiǎn)卧獎(jiǎng)偠染仃嚨耐茖?dǎo)卻鮮有涉及。文獻(xiàn)［5］針對(duì)兩結(jié)點(diǎn)梁?jiǎn)卧o出了變截面梁?jiǎn)卧膯蝿偩仃?，但未明確剛度矩陣的具體推導(dǎo)原理和過(guò)程，也沒(méi)有計(jì)入剪切變形，且所給單剛矩陣僅針對(duì)彎曲問(wèn)題而不涉及伸縮和扭轉(zhuǎn)問(wèn)題。文獻(xiàn)［6］采用傳遞矩陣法推導(dǎo)了任意變截面梁?jiǎn)卧膭偠染仃?，但此方法不涉及形函?shù)，與通用有限元理論不一致。文獻(xiàn)［7］則以形函數(shù)為基礎(chǔ)，借助于MATLAB的符號(hào)運(yùn)算功能建立了梁高成線性或拋物線形變化的變截面梁?jiǎn)卧獎(jiǎng)偠染仃?，但其矩陣元素頗為復(fù)雜，不便使用。ANSYS的Beam44與Beam4的單剛矩陣形式相同，只是前者對(duì)截面參數(shù)包括面積A和慣性矩I采用了“平均值”，但并未給出截面參數(shù)“平均值”的推導(dǎo)原理。
基于有限元理論中的形函數(shù)和最小勢(shì)能原理，在等截面歐拉梁?jiǎn)卧獎(jiǎng)偠染仃囃茖?dǎo)過(guò)程的基礎(chǔ)上［8］，系統(tǒng)闡述了計(jì)入剪切變形的變截面歐拉梁?jiǎn)蝿偩仃囃茖?dǎo)過(guò)程，并經(jīng)過(guò)對(duì)比分析明確了ANSYS中Beam44單元?jiǎng)偠染仃嚨耐茖?dǎo)方法。為便于闡述，僅針對(duì)平面梁?jiǎn)卧M(jìn)行推導(dǎo)。另外，由于ANSYS中的Beam44這一變截面梁?jiǎn)卧蠼孛嬖趦晒?jié)點(diǎn)之間是線性梯度變化的，下文建立的變截面梁同樣符合這一要求。
1 基本方程
以XY平面梁?jiǎn)卧獮槔?，圖1給出了右手螺旋法則的坐標(biāo)系統(tǒng)，并針對(duì)伸縮、扭轉(zhuǎn)和彎曲問(wèn)題給出了對(duì)應(yīng)的桿端荷載和分布荷載，所示荷載方向均與坐標(biāo)軸正方向一致，也即荷載的正方向。圖1（c）還給出了截面內(nèi)力正方向：軸力F以受拉為正，扭矩T的正方向與F方向一致，彎矩M以梁底（即-y方向）纖維受拉為正，剪力Q以使微段發(fā)生逆時(shí)針轉(zhuǎn)動(dòng)為正［1，3-4］。
圖1 單元坐標(biāo)系以及參量正方向示意圖Fig.1 Positive direction regulations in element coordinate system
對(duì)于所示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為單元相對(duì)位置并取-1≤s≤1（圖 1，式（2）），xj和sj為單元j節(jié)點(diǎn)的絕對(duì)坐標(biāo)和相對(duì)坐標(biāo)；L為單元長(zhǎng)度；xc為單元中點(diǎn)坐標(biāo)；Ni（s）和Ni（x）表示拉格朗日形函數(shù)或厄米特形函數(shù)，表1給出兩種形函數(shù)的形式和適用變形［1］。
需要指出，從結(jié)構(gòu)力學(xué)分析可知，變截面梁的形函數(shù)極為復(fù)雜，故在有限元中對(duì)變截面梁依然直接使用等截面梁的形函數(shù)以簡(jiǎn)化分析。
表1 不同變形的形函數(shù)Table 1 Shape functions for different deformation functions
需要說(shuō)明，剪切變形的計(jì)入使梁?jiǎn)卧此拼嬖?種受力模式，但其純剪切和純彎曲受力實(shí)際上是耦合的［1］。為便于闡述，下文依然稱為4種受力模式，并對(duì)每種受力模式分別分析其單剛矩陣，再綜合為任意荷載模式下的單剛矩陣。為節(jié)約篇幅，基于前述位移函數(shù)和形函數(shù)以及最小勢(shì)能原理，以圖1所示坐標(biāo)系和參數(shù)正方向，表2直接給出各種受力模式的單剛矩陣計(jì)算方法，最小勢(shì)能原理以及所涉及的各種受力模式的幾何方程、本構(gòu)方程、截面內(nèi)力和平衡方程詳見(jiàn)文獻(xiàn)［8］；純剪切變形剛度矩陣的推導(dǎo)中，還根據(jù)能量等效原則引入了剪應(yīng)力不均勻修正系數(shù)k（以下簡(jiǎn)稱剪切系數(shù)k）以便使用統(tǒng)一的單剛矩陣計(jì)算公式，詳見(jiàn)文獻(xiàn)［1］。
表2中，E和G為彈性模量和剪切模量，μ為泊松比；截面參數(shù)A、J、I分別為面積、抗扭慣性矩和抗彎慣性矩；由于僅針對(duì)平面梁?jiǎn)卧?，此處的I為繞z軸的抗彎慣性矩也即Iz（圖1、圖2）；y為偏離截面中性軸的坐標(biāo)值；S（y）和b（y）為截面剪應(yīng)力計(jì)算位置一側(cè)截面積對(duì)中性軸的靜矩和計(jì)算位置的截面寬度。對(duì)于變截面梁?jiǎn)卧?，截面參?shù)A、J、I均隨s變化，故在剛度矩陣推導(dǎo)過(guò)程中應(yīng)始終將其放在積分號(hào)內(nèi)；對(duì)于截面剪切系數(shù)k，從其計(jì)算公式可知，其僅與截面形狀有關(guān)，如矩形、圓形和薄壁圓環(huán)截面分別為6∕5、10∕9和2，但部分變截面梁的k將隨s位置變化（圖2），故亦置于積分號(hào)內(nèi)。當(dāng)然，如果為等截面梁，截面參數(shù)A、J、I和剪切系數(shù)k均可置于積分號(hào)外，從而得對(duì)應(yīng)的剛度矩陣如表2所列。下文則以矩形、圓形和薄壁圓環(huán)和薄壁箱形截面為例，分析變截面梁的剛度矩陣，并與ANSYS所給剛度矩陣進(jìn)行對(duì)比。
2 剛度矩陣中的截面參數(shù)分析
2.1 ANSYS中Beam44單元截面參數(shù)取值
在ANSYS中，基于歐拉梁理論提供的等∕變截面梁?jiǎn)卧謩e為Beam4和Beam44，兩者的單剛矩陣形式相同（式（3）、式（4）），只是Beam44剛度矩陣中的截面參數(shù)A、I、J均采用了“平均值”作為等效截面參數(shù)（式（5）），其中的參數(shù)下標(biāo)i和j分別表示左右節(jié)點(diǎn)。或者說(shuō)ANSYS的Beam44本質(zhì)上也是一種等截面梁，只是以左右兩端截面參數(shù)的“平均值”作為截面參數(shù)取值。在ANSYS中，Beam44單元只能通過(guò)實(shí)常數(shù)命令直接定義i和j節(jié)點(diǎn)的截面參數(shù)A、I、J，并對(duì)k只能取一個(gè)值，而無(wú)法通過(guò)命令SECTYPE和SECDATA定義i和j節(jié)點(diǎn)的截面形狀和尺寸。這也說(shuō)明，Beam44單元的單剛矩陣僅由兩端截面參數(shù)和式（3）-式（5）決定，與截面形狀沿桿件的變化規(guī)律無(wú)關(guān)，但從下文的理論分析可知，這樣所得的單剛矩陣對(duì)大部分截面只是一種近似。
表2 不同受力模式的單剛矩陣推導(dǎo)Table 2 Deviation of element stiffness matrixes
2.2 剪切系數(shù)k
對(duì)于截面剪切系數(shù)k，ANSYS中的Beam44單元和Beam4單元均只能輸入一個(gè)參數(shù)，即對(duì)變截面梁認(rèn)為其不隨位置改變。實(shí)際上，對(duì)于變截面梁，矩形、圓形和薄壁圓環(huán)截面的k本身并不隨位置變化，但箱形截面、工形和T形截面的k隨位置變化：圖2依文獻(xiàn)［10］對(duì)一箱形截面梁給出了不同截面尺寸的k值。但是，在單元兩端截面參數(shù)A、J和I比值0.5～2.0的約束下（詳見(jiàn)下文），k的變化較為有限，故下文亦對(duì)k取常數(shù)即取單元中截面的k值，從而可將表2中純剪切剛度矩陣推導(dǎo)中的k置于積分號(hào)以外，同時(shí)對(duì)Beam44單元亦有式（4）。
圖2 箱型截面的截面參數(shù)Fig.2 Section parameters of a box section
2.3 矩形截面參數(shù)分析
對(duì)截面剪切系數(shù)k取定值后，表2中各受力模式剛度矩陣的積分式內(nèi)僅有截面參數(shù)和形函數(shù)導(dǎo)數(shù)向量，但因截面形狀多樣，且截面尺寸隨桿件長(zhǎng)度變化形式多樣，直接對(duì)其進(jìn)行積分仍非常復(fù)雜，且結(jié)果與ANSYS中Beam44單元所得單剛矩陣中的部分元素仍有較大差異。下文首先以最簡(jiǎn)單的矩形截面為例進(jìn)行分析，并且此時(shí)剪切系數(shù)k本身不隨s變化。
取矩形截面變截面梁i、j兩端截面尺寸分別為bj×hj=c1bi×c2hi（式（6）），其中c1和c2分別為截面寬度和高度的變化系數(shù)（以下簡(jiǎn)稱截面變化系數(shù)），假定截面尺寸隨s線性變化（這也與ANSYS中對(duì)Beam44單元的假定一致），從而可得各截面參數(shù)如式（7）所示。
對(duì)于式（7）需要說(shuō)明：對(duì)于本文的平面梁?jiǎn)卧淇箯潙T性矩I也即Iz，但同時(shí)也給出了Iy以便對(duì)比；矩形截面抗扭慣性矩J并沒(méi)有通用的計(jì)算表達(dá)式，往往采用J=χb4計(jì)算，其中系數(shù)χ表達(dá)式極為復(fù)雜，故實(shí)際計(jì)算時(shí)往往根據(jù)高寬比h∕b直接對(duì)χ取值［11］。由于難以給出χ隨s的變化表達(dá)式，故式（7）對(duì)J的表達(dá)式僅適用于截面變化系數(shù)c1=c2，也即截面高寬比h∕b不隨s變化的情況。
將式（7）所示A、J、I帶入表2進(jìn)行積分計(jì)算即可得單一受力模式下的單剛矩陣。從表2所列剛度矩陣推導(dǎo)過(guò)程可知，對(duì)于除純彎曲以外的其他3種受力模式，積分計(jì)算僅針對(duì)截面面積A和抗扭慣矩J。對(duì)于軸向伸縮和純剪切受力，積分計(jì)算僅針對(duì)截面面積A且可得式（8）。對(duì)比式（5）可知，只有當(dāng)截面變化系數(shù)c1=c2時(shí)也即截面寬度和高度按同一規(guī)律變化時(shí)，ANSYS所給等效截面面積才與理論值一致，在c1≠c2時(shí)ANSYS所給等效截面面積只是理論值的一種近似。對(duì)于扭轉(zhuǎn)受力，由于難以給出c1≠c2時(shí)J的計(jì)算表達(dá)式，但c1=c2時(shí)則有式（9）成立，同樣說(shuō)明只有當(dāng)c1=c2時(shí)ANSYS所給等效截面抗扭慣性矩（式（5））才與理論值一致。
純彎受力的積分表達(dá)式更為復(fù)雜，其抗彎慣性矩I和形函數(shù)二階導(dǎo)數(shù)N''均有積分變量s，積分所得表達(dá)式過(guò)于復(fù)雜，且通過(guò)計(jì)算對(duì)比發(fā)現(xiàn)即使c1=c2時(shí)所得剛度矩陣?yán)碚撝狄才cANSYS所得不一致。為明確ANSYS中Beam44單元單剛矩陣的計(jì)算方法，通過(guò)嘗試發(fā)現(xiàn)，如果將分進(jìn)行計(jì)算（式（10）），也即對(duì)幾何矩陣和截面參數(shù)分別積分，則在c1=c2時(shí)所得剛度矩陣將與ANSYS結(jié)果一致（式（11）-式與ANSYS所給等效截面抗彎慣性矩（式（5））一致且所得矩陣與純彎受力模式的剛度矩陣形式一致；兩者聯(lián)合并考慮系數(shù)8E∕L3后即得ANSYS變截面梁純彎剛度矩陣。這實(shí)際也明確了ANSYS中對(duì)Beam44單元純彎模式剛度矩陣的計(jì)算方法，即將″ds進(jìn)行計(jì)算，但對(duì)矩形截面也僅當(dāng)c1=c2時(shí)所得單剛矩陣與ANSYS所給單剛矩陣一致。
明確了矩形截面單一受力模式下ANSYS中單剛矩陣后，即可將其綜合為式（3）所示任意受力模式下的單剛矩陣。而從上述推導(dǎo)可知，對(duì)于矩形截面梁，各受力模式下ANSYS中Beam44單元所得單剛矩陣僅在截面寬度和高度按同一規(guī)律變化也即截面變化系數(shù)c1=c2時(shí)與理論值一致，在c1≠c2時(shí)僅為一種近似。實(shí)際上，ANSYS中Beam44單元對(duì)矩形截面并無(wú)c1=c2的限制，但ANSYS要求Beam44 單元兩端截面參數(shù)的比值 Ai∕Aj、Ji∕Jj、Ii∕Ij均應(yīng)在0.5～2.0之間，否則會(huì)提示警告信息，如果上述比值超出0.1～10.0的范圍，則會(huì)提示錯(cuò)誤信息。這可能也是為了保證當(dāng)c1≠c2時(shí)ANSYS所得等效截面參數(shù)與理論結(jié)果的差異不至過(guò)大而提出的要求：兩端截面參數(shù)比值越接近，兩者的差異也就越小。
2.4 其他常用截面參數(shù)分析
除矩形截面外，圓形、圓環(huán)和箱形截面也是常用的截面形式。對(duì)于圓形截面，因其截面尺寸僅有1個(gè)參數(shù)即半徑R，截面變化系數(shù)同樣僅有1個(gè)即c（式（14））：因圓形截面的Iz=Iy=0.5J=I故僅給出I的表達(dá)式，并且對(duì)比式（7）可知，只需將式（7）的c1和c2取為c即可的圓形截面A、J、I的表達(dá)式，且圓截面的k值亦不隨半徑變化。對(duì)其分析結(jié)果表明，在采用式（10）后，其理論推導(dǎo)所得剛度矩陣與ANSYS采用等效截面參數(shù)的表達(dá)式完全一致。
對(duì)于薄壁圓環(huán)截面，其截面尺寸有2個(gè)參數(shù)即外徑R與壁厚t，相應(yīng)的截面變化系數(shù)同樣為2個(gè)，即 c1和 c2。故其截面尺寸和 A、J、I的表達(dá)式（式（15））與式（7）一致，并且由于是圓環(huán)截面，其J值表達(dá)式不再有c1=c2的限制。對(duì)于箱形截面，其截面尺寸參數(shù)更多，理論上有6個(gè)，但一般頂?shù)装搴穸群妥笥腋拱搴穸确謩e相等故有4個(gè)參數(shù)（圖2），相應(yīng)的截面變化系數(shù)同樣為4個(gè)。對(duì)圓環(huán)和箱形其分析結(jié)果表明，在假定k值不隨桿件位置變化且使用式（10）后，其理論推導(dǎo)所得剛度矩陣同樣在c1=c2=c和c1=c2=c3=c4=c時(shí)與ANSYS采用等效截面參數(shù)的表達(dá)式一致。
需要說(shuō)明，式（15）對(duì)A、J、I的表達(dá)式僅適用于薄壁圓環(huán)，并非精確表達(dá)式（式（16））。如果采用式（16）的精確表達(dá)式，則在c1≠c2時(shí)對(duì)A和I積分不再有式（8）和式（11）的結(jié)果，而式（15）的近似表達(dá)式在c1≠c2時(shí)同樣有式（8）和式（11）成立，但精確表達(dá)式在c1=c2=c時(shí)同樣可得式（5）。
綜上可推知，ANSYS對(duì)Beam44變截面梁?jiǎn)卧獑蝿偩仃囃茖?dǎo)時(shí)：對(duì)純彎模式采用進(jìn)行計(jì)算，同時(shí)對(duì)截面剪切系數(shù)k取定值；所給等效截面參數(shù)A、J和I表達(dá)式源于截面各個(gè)尺寸的截面變化系數(shù)相同時(shí)的情況，或者梁?jiǎn)卧獌啥私孛嫘螤顬橄嗨菩吻医孛娉叽珉S桿件線性變化的情況，其他情況ANSYS所得等效截面參數(shù)與理論值會(huì)有一定偏差。
3 結(jié) 論
基于有限元理論系統(tǒng)闡述了計(jì)入剪切變形的變截面歐拉梁?jiǎn)蝿偩仃囃茖?dǎo)過(guò)程，并以矩形、圓形、圓環(huán)和箱形截面梁為例與ANSYS中Beam44單元?jiǎng)偠染仃囘M(jìn)行了對(duì)比，明確了Beam44單元?jiǎng)偠染仃嚨耐茖?dǎo)方法。研究發(fā)現(xiàn)：
（1）變截面梁因截面形狀多樣且截面尺寸隨桿件長(zhǎng)度變化形式多樣，即使對(duì)于線性梯度變化的變截面梁，對(duì)部分截面也難以給出截面參數(shù)，如抗扭剛度J和剪切系數(shù)k隨桿件位置的表達(dá)式，且變截面梁使用等截面梁的形函數(shù)本身也是一種近似，故難以對(duì)其進(jìn)行積分計(jì)算獲得普適性的剛度矩陣?yán)碚摫磉_(dá)式。
（2）ANSYS對(duì)Beam44變截面梁?jiǎn)卧獑蝿偩仃囃茖?dǎo)時(shí)，純彎模式的積分計(jì)算中對(duì)幾何矩陣和截面參數(shù)分別積分，即以代替進(jìn)行計(jì)算，這樣既可簡(jiǎn)化計(jì)算，又可使各種受力模式的計(jì)算方式一致。
（3）Beam44單元對(duì)截面剪切系數(shù)k取定值以簡(jiǎn)化計(jì)算，故在使用中宜取單元中截面的剪切系數(shù)k或左右兩端截面的均值以減少誤差。
（4）ANSYS所得Beam44單元的單剛矩陣形式與等截面梁完全一致，只是對(duì)截面面積A、抗扭慣矩J和抗彎慣矩I采用了等效值。ANSYS所給等效截面參數(shù)A、J和I表達(dá)式，源于梁?jiǎn)卧獌啥私孛嫘螤顬橄嗨菩吻医孛娉叽珉S桿件線性梯度變化的情況，其他情況的等效截面參數(shù)與理論值會(huì)有一定偏差。ANSYS對(duì)Beam44單元僅要求截面在兩節(jié)點(diǎn)之間按線性梯度變化，并無(wú)要求兩端截面形狀為相似形，這實(shí)際也是將簡(jiǎn)單模式所得表達(dá)式直接推廣應(yīng)用至各種情況。但ANSYS對(duì)Beam44單元要求兩端截面參數(shù)的比率盡可能接近1.0，這可能也是為了控制等效截面參數(shù)與理論值的偏差，故對(duì)變截面梁在建模時(shí)應(yīng)盡量減小單元長(zhǎng)度以提高計(jì)算精度。
