陳代海， 周帥， 李銀鑫， 劉明杰， 李波

(1 鄭州大學 土木工程學院，河南，鄭州 450001; 2.中鐵第四勘察設計院集團有限公司，湖北 武漢 430063)

變截面結(jié)構(gòu)由于其能合理分配自重，充分發(fā)揮材料強度優(yōu)勢，受力合理，因此在實際工程中應用廣泛。在對其進行有限元數(shù)值分析時，需要推導變截面梁單元剛度矩陣，王勖成[1]、葉康生等[2]進行了理論基礎研究；傳光紅等[3]等采用勢能駐值的原理，推導了變截面Timoshenko梁的單元剛度矩陣；Al-Gahtani H J[4]研究了變截面Bernoulli-Euler梁的軸力、彎曲及扭轉(zhuǎn)單元剛度矩陣；張琪等[5]通過Hamilton原理，推導出變截面梁段單元的自由振動方程，得到了單元剛度矩陣；李爽等[6]利用平衡條件得到梁截面內(nèi)力的分布模式，在其變形前建立平衡方程，再由虛力原理導得剛度矩陣；馬志貴等[7]假定變截面的形式為梁高線性變化，根據(jù)有限元基本原理，推導出變截面梁單元的單元剛度矩陣公式；John S,et al[8]討論了一種新型的變截面單元的剛度矩陣，并對其進行試驗研究；陸念力[9]、武芳[10]中采用直接法或傳遞矩陣法推導單元剛度矩陣；傳光紅、陸念力、耿文賓通過將變截面梁單元進行簡化，如單元兩端截面等效為單元中間截面，將截面兩端慣性矩做平均化處理等，從而得出變截面梁單元的剛度矩陣[3,9,11]；Al-Gahtani H J、李爽、John S、武芳沒有給出變截面梁單元剛度矩陣的具體表達式或需要確定的參數(shù)較多[4,6,8,10]，不便于實際應用。因此，該文基于能量變分原理，推導變截面空間梁單元剛度矩陣，得到2結(jié)點與3結(jié)點變截面梁單元剛度矩陣的具體表達式，并通過對比相應Ansys單元類型的計算結(jié)果驗證其正確性。在此基礎上，采用Ansys建立變截面連續(xù)梁模型，分析變截面單元簡化為等截面單元、梁高變化線形及變截面單元劃分長度對變截面連續(xù)梁靜力位移的影響。

1 變截面梁單元剛度矩陣的推導

有限單元法中，常采用2結(jié)點單元或3結(jié)點單元離散變截面梁。該文基于能量變分原理，建立軸力單元、扭轉(zhuǎn)單元及彎曲單元的泛函，根據(jù)彈性力學基本方程及邊界條件，求泛函的駐值，得到變截面梁單元的剛度矩陣表達式。

1.1 2結(jié)點變截面梁單元

對變截面梁單元剛度矩陣的推導，可在局部坐標中考慮梁單元承受的橫向荷載及彎矩。以承受均布荷載q(x)和集中荷載Pj作用的懸臂梁為例，如圖1所示，其撓度函數(shù)為ω(x)。

根據(jù)彈性力學基本方程可得：

(1)

式中:k為梁中面變形后的曲率；M為截面承受的彎矩；Q為截面承受的剪力；θ為端部給定的轉(zhuǎn)角；It為截面慣性矩；E為彈性模量；x為梁端到計算截面的距離。

圖1 彎曲變形梁單元示意圖

根據(jù)能量變分原理建立泛函如下：

(2)

由有限單元法計算梁單元的彎曲剛度矩陣時，采用一維Hermite插值多項式表示撓度插值函數(shù)為：

(3)

式中：N1=1-3ξ2+2ξ3，N2=(ξ-2ξ2+ξ3)l，N3=3ξ2-2ξ3，N4=(ξ3-ξ2)l；l為懸臂梁長度；αi為結(jié)點位移；N為形函數(shù)；αe為αi的矩陣表示形式。

將式(3)代入式(2)，求解δ∏p=0，得到Kα=P，即可求得任意變截面梁單元的彎曲剛度矩陣Kω為：

(4)

推導變截面梁單元的軸力剛度矩陣及扭轉(zhuǎn)剛度矩陣時，采用拉格朗日插值多項式表示其位移函數(shù)為：

(5)

采用類似于梁單元彎曲剛度矩陣的推導方法，可得到變截面梁單元的軸力剛度矩陣Kα及扭轉(zhuǎn)剛度矩陣Kb。

以矩形變截面梁單元為例，其寬度為b，假設截面高度h=h0+mx，其中h0為起始端單元高度，m為梁高線性變化的比例系數(shù)。則沿單元局部坐標軸y軸和z軸的慣性矩分別為Iy=hb3/12和Iz=bh3/12，單元截面面積為A=bh。其軸力剛度矩陣Kα及扭轉(zhuǎn)剛度矩陣Kb為：

(6)

(7)

式中:β為截面扭轉(zhuǎn)系數(shù);G為剪切模量;則單元剛度矩陣Ke可以表示為：

(8)

在實際工程中，受到彎曲變形的影響，梁單元所受橫向剪切力將引起其附加撓度，即變形后的梁單元會發(fā)生翹曲。式(9)中并未考慮剪切變形的影響，且構(gòu)造位移函數(shù)時，截面轉(zhuǎn)動角度為：θ=dω/dx，若在此基礎上考慮剪切變形，其推導過程將較為復雜。因此，在保證精度的前提下簡化推導過程，對撓度ω與截面轉(zhuǎn)動

(9)

角度θ進行獨立插值，推導變截面梁考慮剪切變形時的單元剛度矩陣。

(10)

對式(10)求駐值，可得到梁單元彎剪剛度矩陣為：Kω=Kωy(z)+Ks。

以矩形變截面梁單元為例，參數(shù)同上，推導2結(jié)點變截面梁單元的剛度矩陣。其彎曲剛度矩陣Kωy(z)及剪切剛度矩陣Ks表示如下：

(11)

(12)

(13)

式中:k為橫截面剪應力影響系數(shù)，對于矩形截面，k=6/5。則考慮剪切變形的2結(jié)點變截面梁單元剛度矩陣如式(9)所示。

1.2 3結(jié)點變截面梁單元

3結(jié)點梁單元剛度矩陣的推導方法與2結(jié)點單元類似，其位移函數(shù)采用拉格朗日插值多項式可表示為：

(14)

3結(jié)點矩形變截面單元的軸力剛度矩陣及扭轉(zhuǎn)單元剛度矩陣可表示為：

(15)

(16)

單元彎剪剛度矩陣為：

(17)

(18)

(19)

則3結(jié)點變截面梁單元剛度矩陣可以表示為[見式(20)]：

2 變截面梁單元剛度矩陣的驗證

為驗證上述2結(jié)點單元與3結(jié)點單元剛度矩陣表達式的正確性，以如圖2所示的變截面懸臂梁為研究對象，運用Ansys有限元軟件，分別采用Beam44和Beam189單元建立其有限元模型。懸臂梁長l為1 m，i端截面高hi為0.2 m，寬bi為0.1 m；j端截面高hj為0.3 m，寬bj為0.1 m，彈性模量E為2×108Pa，泊松比μ為0.2，剪切模量G為8.333×107Pa。懸臂梁一端固定，一端自由，自由端施加荷載Fx=50 N，F(xiàn)y=50 N，F(xiàn)z=50 N，Mx=50 N·m，My=50 N·m，Mz=50 N·m。

基于上述推導過程得到單元剛度矩陣K[計算式見式(9)、(20)]，結(jié)合荷載列陣P，求解靜力平衡方程Kα=P，得到懸臂梁自由端的位移，并與Ansys計算結(jié)果進行對比，如表1所示。

由表1可以看出：根據(jù)能量變分法推導的變截面梁單元剛度矩陣，對于2結(jié)點單元與3結(jié)點單元，該文計算結(jié)果與Ansys計算結(jié)果基本一致，驗證了該文變截面梁單元剛度矩陣表達式的正確性，其可用于變截面結(jié)構(gòu)的有限元分析。

3 變截面梁單元在數(shù)值模擬中的影響因素分析

在連續(xù)梁橋的有限元數(shù)值模擬中，變截面梁單元應用較多。以往的有限元分析中，為簡化運算，常將變截面梁單元換算為等截面單元，由此帶來的計算誤差值得研究。另外，由上述變截面梁單元剛度矩陣的表達式可以看出，單元長度和截面高度是影響梁單元剛度的主要因素。因此，該文以3跨變截面連續(xù)梁橋為研究對象，采用Ansys軟件建立其有限元模型，分析變截面單元簡化為等截面單元、梁高變化線形及變截面單元劃分長度等因素對橋梁結(jié)構(gòu)靜位移的影響。

(20)

圖2 變截面懸臂梁示意圖(單位：m)

表1 懸臂梁自由端的位移對比

3.1 工程實例

某變截面連續(xù)梁橋跨度為：(40+64+40) m，橋梁寬7.1 m，變截面段長29.5 m，橋梁縱向布置如圖3所示。邊支點處截面高2.8 m，中支點處截面高5.2 m，截面尺寸如圖4所示。采用C50混凝土，其彈性模量E為3.45×1010Pa，泊松比μ為0.2，材料重度γ為25 kN/m3?？紤]橋梁結(jié)構(gòu)自重和均布荷載Q=100 kN/m的作用。

圖3 變截面連續(xù)梁縱向示意圖(單位：cm)

圖4 變截面連續(xù)梁截面示意圖(單位：cm)

3.2 變截面單元簡化為等截面單元的影響

基于上述工程實例，變截面梁高按2次拋物線變化，對于每個變截面單元，單元長度分別選取為3、4、5和6 m，采用兩種等截面分別對變截面進行簡化：一種等截面定義為等截面1，為變截面單元的中部截面；另一種等截面定義為等截面2，其截面特性在有限元模型中以實常數(shù)的形式表示，實常數(shù)中的截面面積取兩端截面面積的算術(shù)平均值，慣性矩取兩端截面的幾何平均值。用Ansys軟件中的2結(jié)點變截面單元Beam44離散變截面連續(xù)梁，建立變截面單元簡化為等截面單元的連續(xù)梁有限元模型。在均布荷載作用下，計算橋梁的豎向位移，結(jié)果如表2所示。

表2 不同計算模型的橋梁豎向位移對比

由表2可以看出：兩種等截面模型的計算結(jié)果相近，說明兩種等截面均可用于簡化變截面。變截面單元簡化為等截面單元時，當變截面單元劃分長度為3、4 m，即為變截面段長度的1/10和1/7時，等截面與變截面模型的位移值相差較小，最大相對差值小于3%；當變截面單元劃分長度為5、6 m，即為變截面段長度的1/6和1/5時，等截面與變截面模型的最大位移差值較大，最大相對差值為13.48%。因此，在建立變截面梁單元模型時，當單元劃分長度達到變截面段長度的1/6～1/5時，不宜采用等截面對其進行簡化。

3.3 梁高變化線形的影響

變截面連續(xù)梁梁高變化線形常以拋物線為主，為研究拋物線形狀引起的梁高變化對橋梁剛度的影響，采用Ansys軟件中的3結(jié)點單元Beam189離散變截面連續(xù)梁，建立梁高變化線形分別為1.5、1.8和2次拋物線時的3種有限元模型，在均布荷載作用下，計算橋梁的豎向位移，在自重作用下，計算橋梁支座反力。不同梁高變化線形下的橋梁豎向位移對比如圖5所示，支座反力對比如表3所示。

圖5 不同梁高變化線形下的橋梁豎向位移對比

表3 不同梁高變化線形下的橋梁支座反力對比

由圖5可知：在均布荷載作用下，橋梁跨中最大位移值隨著梁高變化拋物線形階次的增加而增大，最大相對差值為6.45%。由表3可知：在自重作用下，橋梁支座反力隨著梁高變化拋物線形階次的增加而減小，最大相對差值為2.5%。說明當梁高變化線形分別為1.5、1.8和2次拋物線時，梁高變化線形對變截面連續(xù)梁的質(zhì)量分布影響較小，而對其剛度分布影響較大。

3.4 變截面單元長度的影響

為分析變截面單元劃分長度對橋梁剛度的影響，采用Ansys中的3結(jié)點單元Beam189離散變截面連續(xù)梁，建立單元長度分別為3、4、5和6 m時的4種有限元模型，計算均布荷載作用下的橋梁豎向位移，得到其豎向位移對比如圖6所示。

圖6 不同變截面單元劃分長度下的橋梁豎向位移

從圖6可以看出：當變截面單元劃分長度為3、4 m，即為變截面段長度的1/10和1/7時，兩種模型的位移值基本吻合；當變截面單元劃分長度為5、6 m，即為變截面段長度的1/6和1/5時，橋梁跨中的最大豎向位移差值較大。當變截面單元劃分長度大于變截面段長度的1/7后，隨著單元長度的增大，橋梁跨中豎向位移最大值逐步減小。因此，為保證計算精度和效率，在建立變截面梁有限元模型時，宜將變截面單元劃分長度控制在變截面段長度的1/7以內(nèi)。

4 結(jié)論

(1) 基于能量變分法，推導了2結(jié)點和3結(jié)點變截面梁單元的剛度矩陣具體表達式。并通過實例計算，對比分析Ansys軟件計算結(jié)果，驗證了其正確性，說明其可用于變截面結(jié)構(gòu)的有限元分析，可為有限元計算程序的編制提供參考。

(2) 在建立變截面連續(xù)梁有限元模型時，宜將變截面單元劃分長度控制在變截面段長度的1/7以內(nèi)。當單元劃分長度達到變截面段長度的1/6～1/5時，不宜采用等截面對其進行簡化。

(3) 當梁高變化線形分別為1.5、1.8和2次拋物線時，梁高變化線形對變截面連續(xù)梁的質(zhì)量分布影響較小，而對其剛度分布影響較大。