用于夾層梁靜動(dòng)力及屈曲分析的新型組合結(jié)構(gòu)單元

林建平，陳昆，潘劍超，王冠楠，馮倩

(1.華僑大學(xué) 土木工程學(xué)院，福建 廈門 361021；2.福建省智慧基礎(chǔ)設(shè)施與監(jiān)測(cè)重點(diǎn)實(shí)驗(yàn)室（華僑大學(xué)），福建 廈門 361021；3.浙江大學(xué) 建筑工程學(xué)院，浙江 杭州 310058；4.公路數(shù)智養(yǎng)護(hù)浙江省工程研究中心，浙江 杭州 310058)

夾層組合梁具有可以充分利用組成材料的優(yōu)點(diǎn)，在保證強(qiáng)度與剛度的基礎(chǔ)上，節(jié)省材料、減輕重量，在正常使用狀態(tài)和承載能力極限狀態(tài)都能顯示出其優(yōu)勢(shì).組合梁的連接件往往不足以保證各層之間實(shí)現(xiàn)完美連接（full-interaction），導(dǎo)致“部分作用（partial-interaction）”[1]的現(xiàn)象，即組合梁受力時(shí)各層之間將會(huì)產(chǎn)生層間的黏結(jié)滑移作用.層間滑移的存在大大增加了組合梁計(jì)算分析的難度，國內(nèi)外學(xué)者對(duì)此展開了廣泛的研究.

以鋼-混凝土組合梁為代表的2層組合梁的研究為工程應(yīng)用提供了寶貴的解決方法.在實(shí)際工程中，3層組合梁也被廣泛應(yīng)用于工程之中，如鋼-混-鋼夾心組合梁[2]、雙鋼板-混凝土組合防護(hù)結(jié)構(gòu)[3]、夾心約束阻尼梁[4]、波紋鋼腹板混凝土梁[5]、夾層玻璃柱[6]等，因此，對(duì)3層組合梁的研究具有重要意義.

Newmark[7]通過實(shí)驗(yàn)提出層間滑移與剪力之間的線性假設(shè)，基于Euler-Bernoulli梁理論建立鋼-混凝土組合梁的控制方程，為部分作用組合梁的研究提供了值得借鑒的研究思路.目前針對(duì)3層組合梁的研究較少，Chui等[8]給出了在均布荷載和集中力荷載作用下3層簡(jiǎn)支組合梁的解析解.Sousa等[9]研究基于Euler-Bernoulli梁和Timoshenko梁理論的多層組合梁的解析解.在數(shù)值分析研究方面，Sousa[10]基于之前的解析方法提出新的有限元公式.Keo等[11]利用控制方程導(dǎo)出精確的剛度矩陣并提出精確有限元公式.Lin等[12]基于Timoshenko梁理論提出考慮獨(dú)立截面轉(zhuǎn)角和橫向剪切變形的3層部分作用組合梁有限元分析方法.

數(shù)值分析方法為組合梁變形分析帶來了便利，但目前關(guān)于3層組合梁有限元數(shù)值分析中可以顯著提高計(jì)算效率的顯式剛度矩陣的研究相對(duì)缺乏，這為其工程應(yīng)用帶來不便.另外，基于位移的組合梁有限元模擬中“滑移/曲率自鎖”[9,12-14]是常見的問題，它會(huì)引起層間滑移場(chǎng)的振蕩從而導(dǎo)致收斂性問題.

本研究在Xu等[15]工作的基礎(chǔ)上，采用Timoshenko梁理論，將高階插值函數(shù)引入內(nèi)部自由度以避免滑移/曲率鎖定問題，并利用MATLAB編譯有限元程序.為了提高計(jì)算效率，推導(dǎo)并給出了3層組合梁的8×8剛度矩陣、質(zhì)量矩陣以及幾何剛度矩陣元素的顯式表達(dá)式.進(jìn)行算例驗(yàn)證和應(yīng)用，通過與文獻(xiàn)和ABAQUS模擬結(jié)果對(duì)比分析，驗(yàn)證本研究提出的新型夾層組合梁有限元方法的準(zhǔn)確性、收斂性和穩(wěn)定性.

1 理論模型

1.1 組合梁結(jié)構(gòu)與參數(shù)

夾層組合梁剖面圖如圖1所示.圖中，Ei、Gi、Ii、Ai和κi(i=1,2,3) 分別表示每層梁的彈性模量、剪切模量、慣性矩、截面積和剪切修正系數(shù)；L0和H表示組合梁的長(zhǎng)度和高度為每層形心軸與相鄰界面的距離；h1、h2為相鄰層結(jié)構(gòu)形心軸的距離，h為頂層與底層形心軸距離，h=h1+h2.在x-z平面內(nèi)，夾層組合梁整個(gè)橫截面的形心軸與x軸重合.

圖1 3層夾層組合梁及坐標(biāo)系示意圖Fig.1 Three-layer sandwich composite beam and coordinate system

推導(dǎo)建立在以下幾個(gè)基本假設(shè)的基礎(chǔ)上[15-16]：1）各層梁均符合線性變形假設(shè)；2）夾層結(jié)構(gòu)的層間剪力與對(duì)應(yīng)的層間滑移成正比；3）忽略界面法向分離現(xiàn)象；4）根據(jù)Timoshenko梁理論，允許存在橫向剪切變形.

1.2 理論推導(dǎo)

如圖2所示為層間滑移、截面轉(zhuǎn)角和軸向位移的幾何關(guān)系，夾層梁的層間滑移與橫截面轉(zhuǎn)角ψ的關(guān)系[17]如下：

圖2 層間滑移、截面轉(zhuǎn)角和軸向位移的幾何關(guān)系Fig.2 Geometrical relationship of interlayer slips,rotation angles and longitudinal displacements

式中：u1、u2和u3為各層的軸向位移，us1、us2分別為中間層與頂層和底層之間的相對(duì)滑移.

如圖3所示為長(zhǎng)度為dx的微小單元.圖中，M、Q、F分別表示整個(gè)橫截面的彎矩、剪力和軸力，Ni、Mi和Qi(i=1,2,3) 為各層梁的軸力、彎矩和剪力，Vi為層間剪力，q、m分別表示作用在該單元上的均布載荷和彎矩.

圖3 3層組合梁的微小單元Fig.3 Infinitesimal element of three-layer composite beam

根據(jù)單元在橫截面上的平衡條件，可以得到

同時(shí)，根據(jù)x方向的平衡條件可以得到

整個(gè)單元橫截面的彎矩可以表示為

各層間縱向剪力與連接件的抗剪剛度呈線性關(guān)系：

根據(jù)Timoshenko梁理論，內(nèi)力與變形關(guān)系為

式中：C=κ1G1A1+κ2G2A2+κ3G3A3，w為梁的撓度.

各層軸向力Ni(i=1，2，3)與相應(yīng)軸向位移的關(guān)系為

1.3 最小勢(shì)能原理

為了推導(dǎo)夾層梁的有限元單元，首先須建立夾層梁的能量表達(dá)式：

式中：π1為彎曲應(yīng)變能，π2為剪切應(yīng)變能，π3為剪切連接應(yīng)變能，Wq為外力做功.表達(dá)式分別如下：

2 組合梁的有限元矩陣推導(dǎo)

2.1 單元?jiǎng)偠染仃囃茖?dǎo)

將組合梁沿軸向離散為若干個(gè)無窮小單元，如圖4所示為其中第e個(gè)單元.記單元節(jié)點(diǎn)位移為ue=[wi,ψi,usi1,usi2,wj,ψj,usj1,usj2]T.單元起點(diǎn) 和終點(diǎn)在x方向的坐標(biāo)分別為xi和xj，因此在局部單元內(nèi)任意點(diǎn)的坐標(biāo)可以表示為

圖4 第e個(gè)單元的坐標(biāo)Fig.4 Coordinate of e-th element

式中：l為該單元的長(zhǎng)度.

所假設(shè)位移的插值函數(shù)表達(dá)式如下：

式中：ξ1、ξ2、η1、η2、η3為新引入的未知參數(shù)，這些參數(shù)只影響單元的局部變形，而非外部節(jié)點(diǎn)或其他單元的變形，可稱為“內(nèi)部自由度”[18]，記ui=[ξ1,η1,ξ2,η2,η3]T.內(nèi)部自由度與節(jié)點(diǎn)自由度的關(guān)系可以通過單元的最小勢(shì)能原理來確定.

將式（18）代入式（11），可以得到新的能量表達(dá)式：

為了消除內(nèi)部的自由度，利用內(nèi)部自由度應(yīng)使單元的勢(shì)能最小化的關(guān)系：

將式（19）重新表示為

進(jìn)而可以得到內(nèi)部自由度與單元節(jié)點(diǎn)上位移的關(guān)系：

利用式(22)，位移表達(dá)式(18)可以重新表示為

將式(23)回代到能量方程（式(11)）并進(jìn)行變分可以得到8×8的對(duì)稱單元?jiǎng)偠染仃嘖e，具體表達(dá)式如附錄B所示.

廣義載荷對(duì)應(yīng)于外力做功表達(dá)式（式(15)）中的線性項(xiàng)，即

將新的插值函數(shù)（式(23)）代入式(24)得到廣義載荷向量Fe.

組合梁的整體剛度矩陣和全局載荷向量表達(dá)式為

式中：N為單位數(shù)量.

可以得到組合梁的總勢(shì)能為

采用最小勢(shì)能原理可以得到

2.2 動(dòng)力有限元分析

組合梁的動(dòng)力行為可以通過哈密頓原理獲得[19]：

式中：符號(hào)“^”表示相應(yīng)的模態(tài)形狀函數(shù).基于哈密頓原理，可以得到自振頻率的變分表達(dá)式：

式中：符號(hào)“st”表示駐值.

與自振頻率分析類似，屈曲荷載Pcr的變分方程可以寫為

將式(23)代入式(31)，沿單元長(zhǎng)度進(jìn)行積分得到不含內(nèi)部自由度的第e個(gè)單元的能量方程，根據(jù)最小勢(shì)能原理可以得到

式中：Me為質(zhì)量矩陣，其元素表達(dá)式見附錄C.

類似于自振頻率的分析，通過對(duì)式(32)進(jìn)行變分可以得到單元的幾何剛度矩陣：

3 算例驗(yàn)證和應(yīng)用

3.1 文獻(xiàn)[9]的算例

首先針對(duì)文獻(xiàn)[9]中如圖5所示的組合梁算例進(jìn)行分析.該梁跨徑L0=4 m，橫截面寬度為0.1 m，梁高為0.24 m.頂層和底層為鋼，彈性模量為200 GPa，層高為0.02 m，泊松比為0.2；中間層為混凝土，彈性模量為34.5 GPa，層高為0.2 m，泊松比為0.3.截面抗剪修正系數(shù)均為5/6.剪切連接件抗剪剛度為ks1=40 MPa，ks2=ks1/8=5 MPa.

圖5 簡(jiǎn)支和連續(xù)3層鋼-混組合梁算例Fig.5 Simply supported and continuous three-layer steel-concrete composite beams

簡(jiǎn)支和連續(xù)2種不同邊界條件下計(jì)算所得的層間最大滑移、最大撓度及連續(xù)梁跨間最大彎矩結(jié)果分別如表1、2所示.表中，δ1為相對(duì)誤差，δ1=[單元數(shù)為100時(shí)本研究結(jié)果-文獻(xiàn)[9]結(jié)果)/文獻(xiàn)[9]結(jié)果]×100%.當(dāng)單元數(shù)量為100時(shí)，計(jì)算所得相對(duì)誤差均不超過0.1%，顯示了本研究計(jì)算結(jié)果的準(zhǔn)確性.在使用100個(gè)組合結(jié)構(gòu)單元時(shí)2個(gè)算例的計(jì)算時(shí)間分別僅為0.184 7、0.309 3 s（Intel Core i5-11 400CPU @ 2.60 GHz），且對(duì)于獨(dú)立轉(zhuǎn)角假設(shè)的含三截面轉(zhuǎn)角模型[12]，計(jì)算時(shí)間分別為0.661 5、0.795 8 s，本研究方法分別節(jié)省了258%、157%的計(jì)算時(shí)間，表明本研究方法的計(jì)算效率高.對(duì)于簡(jiǎn)支梁算例，計(jì)算結(jié)果相對(duì)誤差隨著單元數(shù)的增加而變化的情況如圖6所示.圖中，Ne為單元數(shù)量；δ2為相對(duì)誤差，δ2=[(不同單元數(shù)下計(jì)算結(jié)果-單元數(shù)為500時(shí)計(jì)算結(jié)果)/單元數(shù)為500時(shí)結(jié)果]×100%.計(jì)算結(jié)果顯示，采用本研究方法建立的有限元模型收斂速度快，在4個(gè)單元時(shí)最大相對(duì)誤差為0.03%；當(dāng)單元數(shù)大于22時(shí)，相對(duì)誤差低于0.000 05%.

表1 簡(jiǎn)支邊界條件下層間滑移與最大撓度計(jì)算結(jié)果對(duì)比Tab.1 Comparison of interlayer slip and maximum deflection under simply supported condition

表2 連續(xù)邊界條件下層間滑移與最大彎矩計(jì)算結(jié)果對(duì)比Tab.2 Comparison of interlayer slip and maximum bending moment under continuous boundary condition

圖6 簡(jiǎn)支邊界條件下最大撓度、最大轉(zhuǎn)角和層間滑移相對(duì)誤差隨單元數(shù)的變化Fig.6 Variations of relative errors of maximum deflections,maximum rotation angles and interlayer slips with different element numbers under simply supported condition

如圖7(a)所示為簡(jiǎn)支邊界時(shí)跨中最大撓度隨層間抗剪剛度增大而變化的結(jié)果.有限元模型中單元數(shù)為100，ks1由0.1 MPa增至100 GPa，ks2保持為ks1/8.結(jié)果顯示，當(dāng)抗剪剛度較小時(shí)（ks1≤0.1 MPa），組合梁可以視為無抗剪剛度，跨中撓度超過14 mm；當(dāng)抗剪剛度較大時(shí)（ks1≥100 GPa），組合梁可視為完美抗剪連接，跨中撓度為2.84 mm.如圖7(b)所示為變形縮減系數(shù) ?（不同抗剪強(qiáng)度下變形計(jì)算結(jié)果與ks1=0.1 MPa時(shí)變形計(jì)算結(jié)果的比值）隨抗剪剛度增大的變化趨勢(shì).計(jì)算結(jié)果顯示存在有效抗剪剛度范圍[1 MPa,100 GPa]，可以使得夾層組合梁的變形快速降低.

圖7 變形隨抗剪剛度的變化Fig.7 Variations of deformation with shear stiffness

當(dāng)簡(jiǎn)支梁模型跨高比L0/H=5、10和20時(shí)，采用Euler-Bernoulli梁（EB）與Timoshenko梁（TB）理論所得跨中撓度結(jié)果的相對(duì)誤差隨著抗剪剛度的變化如圖8所示.圖中，δ3為相對(duì)誤差，δ3=[(wmax,EB-wmax,TB)/wmax,TB]×100%，其中，wmax,EB、wmax,TB分別為EB、TB梁理論所得跨中撓度結(jié)果的最大值.結(jié)果顯示EB梁計(jì)算結(jié)果的相對(duì)誤差隨著抗剪剛度的增大而增大.若以相對(duì)誤差3%為界限，在低抗剪剛度（ks1=0.1 MPa）情況下，當(dāng)L0/H=5時(shí)，EB梁計(jì)算結(jié)果不滿足要求；當(dāng)L0/H=10時(shí)，抗剪剛度超過6.3 GPa，EB梁計(jì)算結(jié)果誤差將超限；對(duì)于L0/H=20的細(xì)長(zhǎng)梁，EB梁計(jì)算結(jié)果相對(duì)誤差不超過1%.隨著跨高比的減小，Euler-Bernoulli組合梁將不再適用，Timoshenko梁理論組合梁適用性更強(qiáng).

圖8 最大撓度相對(duì)誤差隨抗剪剛度的變化Fig.8 Variations of relative errors of maximum deflection with shear stiffness

如圖9所示為簡(jiǎn)支梁與連續(xù)梁模型4種不同抗剪剛度時(shí)，采用Euler-Bernoulli梁（EB）與Timoshenko梁（TB）理論所得跨中撓度結(jié)果的相對(duì)誤差隨著跨高比的變化.結(jié)果顯示簡(jiǎn)支梁與連續(xù)梁在跨高比較小時(shí)，EB梁與TB梁計(jì)算結(jié)果的相對(duì)誤差較大.若以3%幅值為界限，當(dāng)ks1=1 MPa時(shí)，簡(jiǎn)支梁和連續(xù)梁的跨高比分別須大于5.4和8.8才可滿足要求；當(dāng)ks1=10 GPa時(shí)，簡(jiǎn)支梁和連續(xù)梁的跨高比則分別須大于11.7和18.8才可滿足要求.

圖9 2種梁理論下跨中最大撓度的相對(duì)誤差Fig.9 Relative errors of maximum mid span deflection obtained by two different beam theories

簡(jiǎn)支邊界時(shí)自振頻率與屈曲荷載的計(jì)算結(jié)果如表3所示，其中ABAQUS模擬時(shí)采用平面應(yīng)力有限元模型，δ4為相對(duì)誤差，δ4=[(本研究結(jié)果-ABAQUS結(jié)果)/ ABAQUS結(jié)果]×100%.結(jié)果顯示本研究方法結(jié)果可與ABAQUS數(shù)值模擬結(jié)果吻合.第1階自振頻率與屈曲荷載的相對(duì)誤差均小于1%；前4階模態(tài)時(shí)相對(duì)誤差均小于3%.如圖10所示為模態(tài)1～5的自振頻率與屈曲荷載的相對(duì)誤差δ5隨單元數(shù)增加的變化，相對(duì)誤差δ5=[計(jì)算結(jié)果-單元數(shù)為500時(shí)結(jié)果)/單元數(shù)為500時(shí)結(jié)果]×100%.對(duì)于第1階自振頻率與屈曲荷載，4個(gè)單元時(shí)相對(duì)誤差即小于0.06%.隨著自振和屈曲失穩(wěn)模態(tài)階數(shù)的提高，所需單元數(shù)量增加.對(duì)于第5階自振頻率與屈曲荷載，為了保證相對(duì)誤差小于0.5%所需的單元數(shù)量分別為16和18.

表3 自振頻率與屈曲荷載計(jì)算結(jié)果對(duì)比Tab.3 Comparison of calculation results of natural frequency and buckling load

圖10 自振頻率與屈曲荷載計(jì)算結(jié)果的相對(duì)誤差Fig.10 Relative errors of natural frequency and buckling load

3.2 文獻(xiàn)[11]的連續(xù)梁算例

本節(jié)所分析的連續(xù)邊界夾層組合梁如圖11所示，橫截面寬度為0.1 m，截面高度共為0.8 m，3層的厚度分別為0.2、0.4、0.2 m.3層材料均相同，彈性模量E=8 GPa，泊松比μ=0.3，截面抗剪修正系數(shù)κ=5/6，上下層間抗剪剛度相同，ks1=ks2.

圖11 兩跨連續(xù)夾層組合梁示意圖Fig.11 Two-span continuous sandwich composite beam

當(dāng)抗剪剛度ks=1 MPa、20 MPa、50 MPa與10 GPa時(shí)，跨中撓度計(jì)算結(jié)果及其相對(duì)誤差隨跨高比的變化如圖12所示.圖中，δ6為相對(duì)誤差，是相對(duì)于采用獨(dú)立截面轉(zhuǎn)角假設(shè)的含三截面轉(zhuǎn)角的有限元模型[12]的結(jié)果，δ6=[(本研究計(jì)算結(jié)果-文獻(xiàn)[12]計(jì)算結(jié)果)/文獻(xiàn)[12]計(jì)算結(jié)果]×100%.當(dāng)跨高比大于8時(shí)，相對(duì)誤差小于1.3%.結(jié)果顯示本研究所提單截面轉(zhuǎn)角模型與三截面轉(zhuǎn)角模型的計(jì)算結(jié)果基本吻合，對(duì)于常規(guī)跨高比的夾層組合梁，本研究方法具有足夠的計(jì)算精度.

圖12 兩跨連續(xù)夾層組合梁的最大撓度及相對(duì)誤差Fig.12 Maximum deflections and its errors of two-span continuous sandwich composite beam

4 結(jié)論

采用Timoshenko梁理論，引入含內(nèi)部變量的高階插值函數(shù)，推導(dǎo)3層組合梁的剛度矩陣、質(zhì)量矩陣以及幾何剛度矩陣元素的顯式表達(dá)式.利用MATLAB編譯有限元程序，驗(yàn)證本研究方法的準(zhǔn)確性和收斂性.

（1）通過算例分析表明，所提方法可以用于考慮界面滑移的夾層組合梁計(jì)算分析.引入含內(nèi)部變量高階插值函數(shù)，避免了滑移鎖定現(xiàn)象.所提方法收斂速度快、所需單元少，簡(jiǎn)支梁靜力算例中的單元數(shù)量取4時(shí)便有較高精度，相對(duì)誤差僅為0.03%.

（2）當(dāng)夾層梁跨高比相同時(shí)，Euler-Bernoulli梁的撓度計(jì)算結(jié)果相對(duì)于Timoshenko梁計(jì)算結(jié)果的誤差隨著抗剪剛度的增大而增大.隨著夾層梁跨高比的減小，采用基于Timoshenko梁理論的組合梁?jiǎn)卧谋匾院蛢?yōu)勢(shì)逐漸凸顯.

（3）相對(duì)于獨(dú)立轉(zhuǎn)角，假設(shè)的截面含多轉(zhuǎn)角有限元模型，所推導(dǎo)的單元自由度更少，計(jì)算速度更快.對(duì)于常規(guī)跨高比的各層材料性能沒有較大突變的夾層梁，所提方法具有足夠計(jì)算精度.

（4）所推導(dǎo)的3層組合梁的剛度矩陣、質(zhì)量矩陣以及幾何剛度矩陣元素的顯式表達(dá)式具有計(jì)算效率高的優(yōu)點(diǎn)，并可推廣用于ABAQUS或ANSYS商業(yè)有限元軟件程序中.

（5）本研究未考慮結(jié)合界面抗剪剛度的非線性力學(xué)行為，計(jì)算分析過程也未考慮材料非線性及幾何非線性，相關(guān)問題值得進(jìn)一步研究.

附錄A：含內(nèi)部自由度單元?jiǎng)偠染仃囎訅K

附錄B：?jiǎn)卧獎(jiǎng)偠染仃嘖e

附錄C：質(zhì)量矩陣 Me