薄壁曲線箱梁空間振動分析方法

蔡 恒，盧海林，顏燕祥，嚴若愚

(1.武漢大學土木建筑工程學院，湖北武漢 430072；2.武漢工程大學土木工程與建筑學院，湖北武漢 430074)

薄壁箱梁廣泛用于現代橋梁結構中。薄壁箱梁的空間受力比較復雜，在對稱荷載作用下其基本變形包括豎向彎曲和剪力滯變形[1]；在偏載作用下，將會產生扭轉和畸變效應[2-4]；對于薄壁箱形曲梁 ，由于存在曲率的影響，不論荷載是否對稱，將會產生相互耦聯的彎曲、扭轉、翹曲、畸變和剪力滯變形，受力復雜。

文獻[5-6]由剛度法原理推導各節點10自由度的直線箱梁空間單元剛度矩陣和剛度方程，考慮了約束扭轉、畸變和剪力滯效應。文獻[7]基于有限元理論提出同時考慮拉壓、彎曲、扭轉、翹曲、畸變、剪力滯效應及其交互作用的每節點14自由度的薄壁曲線箱梁空間分析單元剛度矩陣。文獻[8]建立薄壁曲梁撓曲扭轉分析的彈性控制微分方程組，考慮彎曲、扭翹、畸變和剪力滯效應，得到精確的解析解。文獻[9]直接從殼體結構特點出發，由8 節點曲邊四邊形膜單元和基于Reinssner 中厚板理論的彎曲單元，推導考慮剪切閉鎖效應的薄壁箱梁空間殼單元剛度矩陣及相應的剛度方程，以計算箱梁的剪力滯效應。可以看出上述研究完善了薄壁箱梁的靜力分析理論。

近年來，為了滿足大跨度橋梁結構抗風、抗震需要，學者們對箱形橋梁的動力特性也進行了研究。文獻[10]由廣義坐標法和能量原理推導曲線梁橋振動控制微分方程，給出顯式剛度矩陣和質量矩陣，但未考慮畸變和剪力滯效應。文獻[11]推導曲線橋地震分析的有限單元法，同樣未考慮畸變和剪力滯效應。而現代大跨度PC箱梁橋施工時都不再設橫隔板，使得剛性扭轉和畸變產生的縱向翹曲應力會達到恒載和活載共同作用下縱向彎曲應力的24%～26%[12]。此外還存在剪力滯效應，剪力滯不僅造成箱梁肋板交界處應力集中，引起梁體局部開裂，還會削弱彎箱梁的剛度，引起附加撓度增大[13]。基于此，本文采用動力有限元理論，通過在節點位移中引入高階位移插值函數，綜合考慮薄壁曲線箱梁的拉壓、彎曲、扭轉、翹曲、畸變和剪力滯效應，推導每個節點的11自由度單元剛度矩陣和質量矩陣，求解薄壁曲線箱梁的振動頻率和振型，并采用有限元軟件ANSYS加以驗證，為曲線箱梁橋的抗風和抗震研究提供依據。

1 單元剛度矩陣

1.1 基本假定

為簡化薄壁曲線箱梁的空間力學模型，做如下假定：

(1)材料處于線彈性工作階段；

(2)只在豎向彎曲中計入剪力滯變形；

(3)采用二次拋物線構造剪力滯翹曲位移函數；

(4)不計扭轉翹曲、畸變和剪力滯三者之間小量耦聯變形；

(5)薄壁曲梁截面尺寸相對于跨度和曲率半徑是小量級的。

1.2 剪力滯效應

圖1為薄壁曲梁示意圖，以橫向為x軸，豎向為y軸，縱向為z軸，單元節點位移可表示為

(1)

圖1 薄壁曲梁示意

式中：i、j為單元節點編號；u、v、w分別為箱梁縱向、豎向和橫向位移；v′、w′、φ分別為豎向轉角、橫向轉角和扭轉角；w″為橫向曲率；φ′、γ、γ′分別為扭轉翹曲、畸變和畸變翹曲位移；ζ為翼板剪切位移最大差值。ζ對于考慮剪力滯效應是必要的，按照假定(3)，剪力滯翹曲位移函數表示為

頂板：

(2)

底板：

(3)

懸臂板：

(4)

則箱梁縱向剪力滯翹曲位移可表示為

ξ(x,z,t)=ω(x)ζ(z,t)

(5)

式中：h1、h2分別為形心O到頂板、底板壁厚中線的距離；b1、b2分別為頂板寬度的一半和翼緣板寬度，如圖2所示。

圖2 曲線箱梁截面

圖2中，O1為扭轉中心，其到形心的距離記作e1；O2為畸變中心，是指只發生畸變而無其他位移時，各周邊切向分量對應的轉動中心，其位置可由文獻[14]確定，畸變中心到形心的距離記為e2。

1.3 剛度矩陣

根據彈性力學幾何方程，薄壁曲線箱梁廣義應變可表示為

(6)

寫成矩陣形式為

(7)

式中：等號右邊第一項為微分算子，記作P；等號右邊第二項為箱梁任一點的廣義位移列陣，記作δ；則式(7)可表示為

ε=Pδ

(8)

由文獻[15]，假設橫向位移按五次多項式插值，豎向位移和扭轉角按三次多項式插值，軸向位移按一次多項式插值，這4個插值與式(1)節點位移是相對應的，可以得到位移形函數矩陣

(9)

式中：N1～N12為位移插值函數，具體可表示為

(10)

(11)

(12)

(13)

(14)

(15)

(16)

(17)

(18)

(19)

(20)

(21)

式中：l為單元長度；z為插值點的長度。

因此，單元內任一點位移可表示為

δ=Nδe

(22)

式中：δe為單元節點位移矩陣。將式(22)代入式(8)中，得到

ε=PNδe

(23)

令B=PN，則ε=Bδe，B稱為應變矩陣，薄壁箱梁的彈性矩陣可表示為

(24)

式中：E、G分別為彈性模量和剪切模量；Ix、Iy分別為箱梁截面對x軸和y軸的慣性矩；Id、Iw分別為自由扭轉慣性矩和翹曲扭轉慣性矩；IR、Ir分別為畸變框架慣性矩和畸變翹曲慣性矩。

單元內任一點應力可表示為

σ=Dε=DBδe

(25)

令S=DB，則σ=Sδe，稱S為應力矩陣。

根據虛功原理，單元內任一點虛應變為

ε*=Bδ*e

(26)

單元內應力在虛應變上做的功為

(27)

桿端力在虛位移上做的功為

δW2=δ*eTFe

(28)

由δW1=δW2，得到

(29)

單元的剛度矩陣為

(30)

將應變矩陣B和彈性矩陣D代入式(30)，可得單元剛度矩陣K的表達式為

(31)

2 單元質量矩陣

根據圖2，考慮到質心與扭心、畸心不重合，由牛頓第二定律，箱梁任一點橫向慣性力為

(32)

式中：ρ為箱梁的密度；點表示對時間t的導數。

同樣由于偏心，橫向慣性力會引起附加扭矩和附加畸變矩，扭轉慣性力矩和畸變慣性力矩分別為

(33)

(34)

式中：Iρ為箱梁極慣性矩。箱梁任一點慣性力為

(35)

寫成矩陣形式

(36)

(37)

等效節點力為

(38)

將式(37)代入式(38)中，有

(39)

由此可得質量矩陣為

(40)

式中：l為曲線箱梁跨徑(單元長度)，將M寫成矩陣形式

(41)

3 振動頻率方程

把薄壁曲梁沿z軸方向劃分為若干個單元，每個單元有2個節點i和j，每個節點有11個自由度，如圖3所示。

圖3 薄壁曲線箱梁離散元

將單元剛度矩陣和質量矩陣按對號入座法則進行組裝形成總剛度矩陣、總質量矩陣。本文采用流動圓柱坐標系，故不需要對薄壁曲梁進行坐標轉換。

薄壁曲梁的振動方程為

(K-ω2M)δe=0

(42)

式中：ω為振動頻率，當薄壁曲梁發生自由振動時，式(42)有非零解，故有

K-ω2M=0

(43)

根據式(43)，引入相應的邊界條件，便可以求得薄壁曲梁振動頻率。

4 算例

以文獻[12]薄壁簡支箱梁為例，截面尺寸如圖4所示，跨徑為40 m，計算過程中認為它是曲線梁，即不改變截面尺寸、邊界條件和跨徑，按曲率半徑R=100 m的曲線梁考慮。

圖4 箱梁截面尺寸(單位：m)

將薄壁曲線箱梁沿軸向劃分為3個單元，采用MATLAB編制計算程序求解特征值方程。由于ANSYS中的Beam188梁單元不能反映結構的扭轉振動特性，故采用Shell63單元建立有限元模型，在腹板與底板交界處結點施加約束，兩端約束Ux、Uy、Uz。在結構的動力分析中，主要考慮低階模態，故取前5階自振頻率，結果見表1。

表1 簡支曲線箱梁自振頻率

注：誤差=|本文解-ANSYS結果|/ANSYS結果。

由表1可以看出，本文解與ANSYS結果存在一定的誤差，主要是因為在ANSYS中約束位于梁端腹板與底板交界處的結點，而理論計算中約束的是截面的中心，二者邊界條件約束存在差異，但誤差大部分不超過5%，理論計算結果與ANSYS有限元結果吻合良好，表明本文所構造的薄壁曲線箱梁單元剛度矩陣和質量矩陣是正確可靠的。此外，當曲率半徑R趨于無窮時，便可計算直線箱梁的自振特性，所提出的方法具有通用性，優于以往僅能計算直線箱梁自振特性的方法，限于篇幅，不再舉例說明。

5 結論

(1)提出一種綜合考慮曲線箱梁拉壓、彎曲、扭轉、翹曲、畸變和剪力滯變形的單元剛度矩陣和質量矩陣以計算自振頻率，經算例驗證與ANSYS有限元方法結果吻合良好，表明本文方法是正確可靠的。

(2)所提出的方法沿梁長方向僅劃分為數個單元便取得滿意的效果，表明本文方法的高效性。

(3)提出的單元剛度矩陣和質量矩陣不僅可以計算單跨和多跨曲線箱梁橋的自振頻率，還可以結合車輛模型或者地震質量矩陣，形成結構的達朗伯方程，采用Newmark-β法求解車-橋耦合振動響應或者結構的地震響應。