一種多層組合梁自由振動分析的混合有限元法

蔣佳卿 王云 陳偉球 徐榮橋

摘要： 基于截面等效，提出了多層組合梁的二維分析模型。使用混合能變分原理，以包含頻率函數(shù)的節(jié)點(diǎn)位移及其能量對偶的應(yīng)力分量為單元節(jié)點(diǎn)未知量，引入混合元對梁長方向進(jìn)行離散，建立動力學(xué)狀態(tài)空間控制方程；接著引入微分求積法（Differential Quadrature Method， DQM）對控制方程進(jìn)行梁高度方向的離散，求得組合梁在不同軸力與邊界條件下的動力學(xué)方程。以混凝土?木材組合梁、波形鋼腹板梁和鋼?混組合梁為例進(jìn)行驗(yàn)證。該方法基于二維理論，可為梁理論提供假設(shè)依據(jù)和誤差分析的基準(zhǔn)。

關(guān)鍵詞： 自由振動； 組合梁； 混合有限元； DQM； 狀態(tài)空間法

中圖分類號： O327??? 文獻(xiàn)標(biāo)志碼： A??? 文章編號： 1004-4523（2024）05-0856-08

DOI：10.16385/j.cnki.issn.1004-4523.2024.05.014

引? 言

組合梁的動力特性是組合梁計(jì)算分析中的重要內(nèi)容，利用自振頻率可以識別結(jié)構(gòu)可能出現(xiàn)的損傷?；赥imoshenko梁理論，Xu等 ［1］推導(dǎo)了考慮層間滑移的組合梁解析表達(dá)式，分析了其靜力、動力與屈曲特性?？紤]高階剪切效應(yīng)的影響，Vo等［2］提出了使用高階多項(xiàng)式和非多項(xiàng)式函數(shù)的高階剪切變形理論。Carrera等［3］比較多種不同高階函數(shù)，包括多項(xiàng)式、三角級數(shù)、指數(shù)函數(shù)和Zig?zag函數(shù)，分析多層組合梁自振特性的異同?；赯ig?zag理論，胡霖遠(yuǎn)等［4］通過引入分段的位移線性分布假設(shè)與橫向剪應(yīng)力拋物線分布假設(shè)，建立了適用于分析波形鋼腹板梁的自由振動分析理論。Shen等［5?6］基于狀態(tài)空間架構(gòu)，給出了考慮部分組合作用的組合梁自由振動控制方程，并分析了軸力的影響。通過引入“動力折減系數(shù)”與“頻率折減系數(shù)”，侯忠明等［7］提出了鋼?混組合梁動力問題的簡便計(jì)算方法?？紤]到鋼?混組合梁在鋼梁翼緣存在剪力滯效應(yīng)，陳玉驥等［8］通過假設(shè)結(jié)構(gòu)的縱向位移函數(shù)，提出了考慮剪力滯與層間滑移效應(yīng)的組合梁一階自振頻率計(jì)算方法。

除了解析方法，有限元方法在組合梁動力分析中得到了較多應(yīng)用?；诮?jīng)典梁理論，Vo等［9］使用一維有限元單元分析了薄壁組合梁的動力特性。Chakrabarti等［10］提出了考慮層間滑移的高階梁單元。Chalak等［11］提出了基于Zig?zag理論的有限元單元，并用其分析了包含柔軟夾層的三明治結(jié)構(gòu)組合梁的動力特性。

以上研究都基于一維梁理論，若將組合梁近似為二維結(jié)構(gòu)，則無須假定截面上的變形或應(yīng)力的分布，適用于截面材料特性變化較大、層數(shù)較多，或是邊界條件沿梁高分布較為復(fù)雜的情況。Chen等［12］基于二維彈性理論，提出了使用狀態(tài)空間法與微分求積法的組合梁自由振動分析方法?；跔顟B(tài)空間法，Xu等［13］給出了簡支條件下考慮層間滑移的組合梁自由振動與屈曲荷載解析解，并使用微分求積法給出了在其他邊界條件下的半解析表達(dá)式。因?yàn)閷⑽灰婆c其能量對偶的應(yīng)力分量作為基本未知量，因此狀態(tài)空間法在分析多層組合結(jié)構(gòu)時具有獨(dú)特優(yōu)勢。

本文在組合梁統(tǒng)一的二維模型基礎(chǔ)上，推導(dǎo)了哈密頓體系下考慮固有頻率與軸力作用的混合能量泛函。然后通過梁長方向的有限元離散，求得沿梁高度方向的狀態(tài)空間方程。接著利用微分求積法（Differential Quadrature Method， DQM）對狀態(tài)方程在梁高方向上進(jìn)行分解。最后使用二分法解得組合梁的各階自振頻率。通過對混凝土?木材組合梁、波形鋼腹板梁和鋼?混組合梁的分析，驗(yàn)證了該方法的正確性。

1 組合梁的二維模型

不同形式的組合梁，包括鋼?混組合梁與波形鋼腹板梁，均可通過截面等效轉(zhuǎn)化為寬度為的矩形截面組合梁。參考胡霖遠(yuǎn)等［4］、Johnson等 ［14］ 以及Xu等［15］的工作，通過把波形鋼腹板等效為正交各向異性材料，可模擬波形鋼腹板橫向和縱向剛度的巨大差異。同時根據(jù)截面抗彎剛度與軸向剛度一致原則，對不同寬度的混凝土截面或鋼梁截面進(jìn)行材料參數(shù)轉(zhuǎn)換，使其統(tǒng)一等效為矩形截面。

為保持模型的一致性，組合梁的層間部分組合作用通過引入厚度極小的虛擬粘接層進(jìn)行模擬。如圖1中高度為（其值遠(yuǎn)小于其他層的高度）的薄層，其剪切模量為：

（1）

式中為剪力連接件的抗剪剛度，該式可方便地根據(jù)能量相等原則推導(dǎo)得到。如果需要考慮實(shí)際情況下剪力連接件在梁長方向的離散分布特性，那么虛擬粘接層的剪切模量也可在梁長方向作相應(yīng)的設(shè)置。通過一個等效的薄層材料模擬由于剪力連接件剛度有限產(chǎn)生的層間相對滑移，使得所有的等效截面成為統(tǒng)一的層間完美連接的二維模型。

2 混合有限元與DQM結(jié)合解法

2.1 混合能泛函

考慮組合梁在軸力作用下的自由振動問題，假設(shè)梁的振動頻率為，則二維系統(tǒng)的拉格朗日函數(shù)可表示為：

式中為位移振幅列陣，u和v分別對應(yīng)x方向和y方向的位移振幅；若變量下標(biāo)中包含逗號，則表示該變量對逗號后面坐標(biāo)的偏導(dǎo)數(shù)；表示材料的密度；矩陣的取值可參考胡霖遠(yuǎn)等 ［4］的工作。

選取方向作為傳遞方向，通過勒讓德變換［16］可得到如下自由振動問題的哈密頓函數(shù)：

式中? 為法線方向與x軸平行的平面上的應(yīng)力振幅矢量，以及：

求得哈密頓函數(shù)表達(dá)式后，混合能泛函可表示為：

式中表示平面區(qū)域。該混合能泛函以位移振幅和應(yīng)力振幅矢量為自變函數(shù)，使得該泛函取極值的可能位移振幅和應(yīng)力振幅是精確解。出于簡便性，后面的文字描述中統(tǒng)一省略“振幅”。

2.2 混合有限元離散

對結(jié)構(gòu)沿x方向進(jìn)行有限元離散，y方向保持不變，并對自由振動下的單元內(nèi)位移與應(yīng)力做如下假設(shè)：

（6a）

（6b）

式中表示位移與應(yīng)力隨時間的變化函數(shù)，i和t分別表示虛數(shù)與時間變量；和分別表示單元節(jié)點(diǎn)上的位移和應(yīng)力分量；和分別表示單元內(nèi)部位移與應(yīng)力的分布形函數(shù)。

有限元離散后，系統(tǒng)的混合能泛函也可以寫為：

式中? H為整個梁高；m為x方向單元數(shù)量；和為第e個單元所占的區(qū)域。

將式（3）代入式（7），并結(jié)合式（6），得：

式中

求得單元內(nèi)部的系數(shù)矩陣后，按節(jié)點(diǎn)號組裝全部矩陣，可得系統(tǒng)的混合能泛函為：

式中和表示系統(tǒng)全體位移與應(yīng)力分量；，，，，表示根據(jù)有限元節(jié)點(diǎn)號組裝后的各系數(shù)矩陣。

對式（10）取變分為零，可得如下控制方程：

式中

考慮到位移與應(yīng)力具有不同的量綱，為提高數(shù)值精度，可按下面方式進(jìn)行無量綱化。首先引入：

式中為特征長度，可取為；E為具有彈性模量量綱的量。

按式（13）無量綱化后，式（11）可改寫為：

式中

2.3 三節(jié)點(diǎn)單元

有限元離散采用三節(jié)點(diǎn)有限元單元，如圖2所示，并取中間節(jié)點(diǎn)k為首節(jié)點(diǎn)i與末節(jié)點(diǎn)j之間的中點(diǎn)。形函數(shù)表達(dá)式如下：

（16）

式中

表示單元內(nèi)的局部坐標(biāo)，與全局坐標(biāo)x的關(guān)系為：

由于形函數(shù)由參考坐標(biāo)給出，因此根據(jù)復(fù)合函數(shù)求導(dǎo)規(guī)則：

（19）

將形函數(shù)表達(dá)式代入式（9），可得各單元系數(shù)矩陣的表達(dá)式如下：

2.4 狀態(tài)方程的DQM求解

狀態(tài)方程（14）為微分方程，雖可用傳遞矩陣直接求解，但是由于引入有限單元后維數(shù)較大，容易導(dǎo)致數(shù)值穩(wěn)定問題。為此，可采用微分求積法［17?18］求解狀態(tài)方程（14）。首先， 在組合梁第i層上布置k個離散點(diǎn)，其局部坐標(biāo)為（）。根據(jù)DQM法，在離散點(diǎn)c上，狀態(tài)向量的導(dǎo)數(shù)滿足：

式中和表示第i層內(nèi)和點(diǎn)處的狀態(tài)量；為DQM權(quán)重系數(shù)，一種取法為［18］：

對于第i層上任一離散點(diǎn)上的狀態(tài)量，式（14）可以改寫為：

式中表示式（14）中矩陣A在第i層的取值。將式（21）代入式（23），可得在任一離散點(diǎn)c上滿足：

組合第i層上所有離散點(diǎn)之間的關(guān)系，可得：

（25）

式中

式中? I表示與維數(shù)相同的單位矩陣。

考慮到原組合結(jié)構(gòu)具有n層，將各層系數(shù)矩陣依次代入式（25），并組合可得：

（27）

式中

式中? 系數(shù)矩陣表示組合梁的自身動力特性，還需針對不同邊界條件進(jìn)行相應(yīng)修改。

2.5 邊界條件引入

邊界條件通過有限元分析方法中常用的“乘大數(shù)法”進(jìn)行引入，并記大數(shù)為。設(shè)邊界條件共約束f個狀態(tài)量，即：

，（29）

式中表示給定狀態(tài)量；表示給定值。

通過“乘大數(shù)法”引入邊界條件式（29）后，式（27）可改寫為：

（30）

式中? I表示與維數(shù)相同的單位矩陣；表示維數(shù)為的列向量，其中n表示材料層數(shù)，k表示每層的DQM離散點(diǎn)個數(shù)，M表示有限元節(jié)點(diǎn)數(shù)。滿足：

式中表示狀態(tài)量在總體狀態(tài)量矢量中的位置。

結(jié)合式（12），（15）和（26b），可知式（30）中的系數(shù)矩陣為頻率項(xiàng)和軸力的函數(shù)。為使得方程（30）有非零解，系數(shù)矩陣應(yīng)滿足：

（32）

計(jì)算組合梁在無軸力作用下的自振頻率，可通過設(shè)置式（12）中的軸力項(xiàng)，然后通過二分法求得滿足式（32）的組合梁各階自振頻率。通過設(shè)定式（12）中的頻率，可相應(yīng)求得組合梁的屈曲荷載。

3 算? 例

3.1 混凝土?木材組合梁

混凝土?木材組合梁由上層混凝土與下層木梁組成。上層混凝土的材料參數(shù)為，，下層木梁的材料參數(shù)為，，層間具有剪切剛度。表1展示了四種不同邊界條件下，即兩端簡支（SS）、一端簡支一端固支（SC）、兩端固支（CC）、一段固支一段自由（CF），組合梁的前五階自振頻率，并將計(jì)算結(jié)果與Xu等 ［1］解析解的結(jié)果進(jìn)行對比。

混凝土與木梁間通過剪力連接件相連接，連接件剛度決定了組合作用的強(qiáng)弱。當(dāng)連接件剛度較小時，組合作用較弱，此時組合梁的固有頻率較低；當(dāng)連接件剛度逐漸增大時，組合作用隨之增強(qiáng)，固有頻率也相應(yīng)提高。當(dāng)連接件剛度很大時，層間近似完美連接，此時固有頻率達(dá)到穩(wěn)定的最大值。

圖3反映了四種不同邊界條件下，組合梁一階頻率隨剪切剛度而變化的趨勢圖。不同支座條件下，一階頻率發(fā)生顯著變化時對應(yīng)的連接剛度區(qū)間不同。兩端簡支（SS）下，其頻率敏感的連接剛度變化區(qū)間為，而在兩端固支（CC）條件下，其頻率敏感的連接剛度變化區(qū)間為。

3.2 波形鋼腹板梁

基于二維分析模型，組合梁內(nèi)部位移與應(yīng)力的分布滿足傳遞方程式（14），因而無須預(yù)先引入其分布的假設(shè)函數(shù)，適用于分析具有各向異性特征的波形鋼腹板梁結(jié)構(gòu)。參考胡霖遠(yuǎn)等［4］、Johnson等［14］ 以及Xu等［15］的工作，波形鋼腹板可轉(zhuǎn)換為等效正交各向異性材料，并將波形鋼腹板梁（Corrugated Steel Web Beams， CSWB）等效為單位寬度下的矩形三層組合結(jié)構(gòu)。各層的材料參數(shù)如下：混凝土頂板的等效彈性模量為，泊松比為，密度為。波形鋼腹板的等效彈性模量為，，；泊松比為，；密度為?；炷恋装宓牡刃椥阅Ａ繛?img alt="" height="17" src="file：///C：/Users/examp/AppData/Local/Temp/msohtmlclip1/01/clip_image270.gif" width="98"/>，泊松比為，密度為。波形鋼腹板梁長，高，其中混凝土頂板與底板厚度均為0.25 m。

基于式（32）求得的各階自振頻率結(jié)果，與Zig?zag理論（Zig?zag Beam Theory， ZBT）［4］和有限元數(shù)值方法計(jì)算所得的結(jié)果對比如表2～5所示。表中f表示自振頻率，括號內(nèi)表示該方法與有限元結(jié)果的相對誤差（%）。有限元結(jié)果來自于ABAQUS軟件使用八節(jié)點(diǎn)平面應(yīng)力單元的計(jì)算結(jié)果，上、下混凝土板的單元尺寸均為0.025 m×0.450 m，波形鋼腹板層的單元尺寸為0.200 m×0.450 m，在建立波形鋼腹板梁的平面應(yīng)力有限元模型后計(jì)算其動力特性。

不同邊界條件下，基于此方法的計(jì)算結(jié)果與有限元結(jié)果吻合度較高，各階頻率值均與有限元結(jié)果接近。而Zig?zag理論的結(jié)果在低階頻率上精度較高；隨著頻率階數(shù)的升高，其精度有所下降。這是由于隨著振型階數(shù)升高，組合梁內(nèi)的剪切變形呈現(xiàn)出更顯著的非線性特征，因此使用Zig?zag理論所得結(jié)果的誤差將逐漸增大。基于二維模型的傳遞方程式（14），無須引入變形和應(yīng)力的預(yù)先分布假設(shè)，因此對于各階頻率均能保持較低的誤差水平。為了保證總體精度，本算例統(tǒng)一選用較大的DQM離散點(diǎn)數(shù)（k=20）并在材料界面附近加密布置離散點(diǎn)，偏多偏密的離散點(diǎn)布置使得低階頻率的誤差反而偏大。

3.3 鋼?混組合梁

鋼?混組合梁在橋梁工程中應(yīng)用廣泛。通過在界面處布置抗剪連接件，鋼?混組合梁將鋼材與混凝土進(jìn)行組合，發(fā)揮其各自優(yōu)勢?？紤]一跨徑的鋼?混凝土試驗(yàn)組合梁。鋼梁為的標(biāo)準(zhǔn)H型鋼，容重，彈性模量和泊松比分別為和?；炷梁?img alt="" height="17" src="file：///C：/Users/examp/AppData/Local/Temp/msohtmlclip1/01/clip_image292.gif" width="78"/>，寬，容重，彈性模量和泊松比分別為和?；炷僚c鋼梁之間通過栓釘（直徑）連接，并沿跨徑方向每間隔0.2 m布置一排，每排2個，產(chǎn)生的界面抗剪強(qiáng)度。參考胡霖遠(yuǎn)等［4］和Xu等［15］的工作，可將該組合梁等效為單位寬度下的矩形五層組合結(jié)構(gòu)。對于五層結(jié)構(gòu)的分析，體現(xiàn)了該方法可方便地處理任意多層結(jié)構(gòu)的特點(diǎn)。

由于栓釘沿梁長方向?yàn)殚g隔布置，采用兩種方式對其進(jìn)行模擬。方式（Ⅰ）首先計(jì)算每排栓釘處的界面抗剪強(qiáng)度：

式中表示梁長方向的栓釘排數(shù)；表示等效界面的寬度，即。求得抗剪強(qiáng)度后，可按式（1）求得栓釘位置處虛擬粘接層的剪切模量。而其他無栓釘處的剪切模量則取為一較小值（）。在進(jìn)行有限元離散時，通過設(shè)定單元間不同的剪切模量，體現(xiàn)了該方法可處理材料特性沿梁長變化的優(yōu)勢。方式（Ⅱ）則不考慮栓釘?shù)拈g隔布置，由式（1）直接求得對應(yīng)于的全梁長等效剪切模量。

表6給出了兩種計(jì)算方式下組合梁的固有頻率計(jì)算結(jié)果，并與有限元數(shù)值結(jié)果以及張?jiān)讫埖?［19］的室內(nèi)試驗(yàn)結(jié)果進(jìn)行對比。有限元數(shù)值計(jì)算基于ANSYS軟件，工字鋼與混凝土使用實(shí)體單元（SOLID65）進(jìn)行離散，栓釘使用梁單元 （BEAM188） 模擬。對比表明，采用方式（Ⅰ）的計(jì)算值與試驗(yàn)實(shí)測結(jié)果吻合度較高，誤差僅為0.64%。而不考慮栓釘間隔布置的全梁長等效方式（Ⅱ），其頻率計(jì)算結(jié)果較試驗(yàn)實(shí)測值偏低。這是由于將抗剪強(qiáng)度進(jìn)行全梁長均勻化等效，降低了端部處的連接剛度，使得組合梁的固有頻率下降。

4 結(jié)? 論

基于組合梁統(tǒng)一的二維分析模型，導(dǎo)出了哈密頓體系下組合梁動力分析的混合能量泛函。以包含頻率函數(shù)的節(jié)點(diǎn)位移及其能量對偶的應(yīng)力分量為單元節(jié)點(diǎn)未知量，使用有限元離散建立了狀態(tài)空間法下的混合有限元控制方程。狀態(tài)空間法建立了狀態(tài)量在傳遞方向上解析的傳遞關(guān)系，將復(fù)雜的多層結(jié)構(gòu)界面連接關(guān)系轉(zhuǎn)換為簡單的矩陣相乘，對分析任意多層組合梁具有顯著優(yōu)勢。而通過引入梁長方向的有限元離散，可方便地處理材料特性沿梁長變化的問題。最后，通過在傳遞方向上引入DQM方法，可有效防止由于引入有限元離散而導(dǎo)致的待解矩陣過大產(chǎn)生的數(shù)值穩(wěn)定性問題。

部分連接的組合梁以及具有正交各向異性特征的波形鋼腹板梁是工程中常見的組合梁形式，該方法準(zhǔn)確分析了其在多種邊界條件下的自振頻率。對比分析波形鋼腹板梁常用的Zig?zag理論，由于基于二維分析模型，該方法對高階頻率的計(jì)算十分準(zhǔn)確，且隨著頻率的升高，其優(yōu)勢更加明顯。而對于二層、三層及五層組合結(jié)構(gòu)的分析，表明分析過程具有通用性，并可方便地適用于任意多層組合結(jié)構(gòu)的自由振動分析。通過設(shè)置梁長方向上有限元單元之間不同的材料參數(shù)，可精細(xì)化處理材料特性沿梁長方向非均勻分布的情形，使計(jì)算模型及其結(jié)果與實(shí)際情況更加符合。

參考文獻(xiàn)：

［1］Xu R， Wu Y. Static， dynamic， and buckling analysis of partial interaction composite members using Timoshenko's beam theory［J］. International Journal of Mechanical Sciences， 2007， 49（10）： 1139-1155.

［2］Vo T P， Thai H T. Vibration and buckling of composite beams using refined shear deformation theory ［J］. International Journal of Mechanical Sciences， 2012， 62（1）： 67-76.

［3］Carrera E， Filippi M， Zappino E. Free vibration analysis of laminated beam by polynomial， trigonometric， exponential and Zig-zag theories［J］. Journal of Composite Materials， 2014， 48（19）： 2299-2316.

［4］胡霖遠(yuǎn)，陳偉球，張治成，等. 基于 Zig-zag 理論的波形鋼腹板梁自由振動分析［J］. 浙江大學(xué)學(xué)報(bào)（工學(xué)版），2019，53（3）： 503-511.

HU Lin-yuan， CHEN Wei-qiu， ZHANG Zhi-cheng， et al. Free vibration analysis of concrete beams with corrugated steel webs based on Zig-zag theory ［J］. Journal of Zhejiang University （Engineering Science）， 2019， 53（3）： 503-511.

［5］Shen X， Chen W， Wu Y， et al. Dynamic analysis of partial-interaction composite beams［J］. Composites Science and Technology， 2011， 71（10）： 1286-1294.

［6］沈旭棟，陳偉球，徐榮橋. 有軸力的部分作用組合梁的動力分析［J］. 振動工程學(xué)報(bào)，2012，25（5）： 514-521.

SHEN Xu-dong， CHEN Wei-qiu， XU Rong-qiao. Dynamic analysis of partial-interactive composite beam with axial force［J］. Journal of Vibration Engineering， 2012， 25（5）： 514-521.

［7］侯忠明，夏禾，王元清，等. 鋼-混凝土組合梁動力折減系數(shù)研究 ［J］. 振動與沖擊，2015，34（4）： 74-81.

HOU Zhong-ming， XIA He， WANG Yuan-qing， et al.? Dynamic reduction coefficients for a steel-concrete composite beam ［J］. Journal of Vibration and Shock， 2015， 34（4）： 74-81.

［8］陳玉驥，羅旗幟. 鋼-混凝土組合梁考慮剪力滯和滑移的一階自振頻率［J］. 煙臺大學(xué)學(xué)報(bào)（自然科學(xué)與工程版），2015，28（1）： 49-53.

CHEN Yu-ji， LUO Qi-zhi. First order natural vibration frequency for steel-concrete composite beams considering influence of shear lags and slips［J］. Journal of Yantai University （Natural Science and Engineering Edition）， 2015， 28（1）： 49-53.

［9］Vo T P， Lee J. Free vibration of axially loaded thin-walled composite Timoshenko beams［J］. Archive of Applied Mechanics， 2011， 81（9）： 1165-1180.

［10］Chakrabarti A， Sheikh A H， Griffith M， et al. Dynamic response of composite beams with partial shear interaction using a higher-order beam theory［J］. Journal of Structural Engineering， 2013， 139（1）： 47-56.

［11］Chalak H D， Chakrabarti A， Iqbal M A， et al. Vibration of laminated sandwich beams having soft core［J］. Journal of Vibration and Control， 2012， 18（10）： 1422-1435.

［12］Chen W Q， Lü C F， Bian Z G. Elasticity solution for free vibration of laminated beams［J］. Composite Structures， 2003， 62（1）： 75-82.

［13］Xu R， Wu Y F. Free vibration and buckling of composite beams with interlayer slip by two-dimensional theory ［J］. Journal of Sound and Vibration， 2008， 313（3-5）： 875-890.

［14］Johnson R P， Cafolla J， Bernard C. Corrugated webs in plate girders for bridges ［J］. Proceedings of the Institution of Civil Engineers-Structures and Buildings， 1997， 122（2）： 157-164.

［15］Xu R， Wu Y F. Two-dimensional analytical solutions of simply supported composite beams with interlayer slips［J］. International Journal of Solids and Structures， 2007， 44（1）： 165-175.

［16］Courant R， Hilbert D. Methods of Mathematical Physics： Partial Differential Equations［M］. Singapore： John Wiley & Sons， 2008： 32-34.

［17］Bellman R E， Kashef B G， Casti J. Differential quadrature： A technique for the rapid solution of nonlinear partial differential equations［J］. Journal of Computational Physics， 1972， 10： 40-52.

［18］Shu C. Differential Quadrature and Its Application in Engineering［M］. London： Springer Science & Business Media， 2012： 29-31.

［19］張?jiān)讫?，郭陽陽，王靜，等. 鋼-混凝土組合梁的固有頻率及其振型［J］. 吉林大學(xué)學(xué)報(bào)（工學(xué)版），2020，50（2）： 581-588.

ZHANG Yun-long， GUO Yang-yang， WANG Jing， et al. Natural frequency and mode of vibration of steel?concrete composite beam［J］. Journal of Jilin University （Engineering and Technology Edition）， 2020， 50（2）： 581-588.

A mixed finite element method for free vibration analysis of multilayer composite beams

Abstract： In this paper， a two-dimensional analytical model for composite beams is first proposed through the equivalent transformation of the cross-section. Based on the mixed variational principle， the dynamic state equations are derived through finite element meshing and interpolation along the length of the beam， with frequency contained nodal displacements and their energy-conjugated stresses as element nodal variables. The differential quadrature method （DQM） is introduced to discretize the equations along the height of the beam， and natural frequencies of composite beams under different axial forces and boundary conditions are obtained. This method was verified by numerical examples about natural frequencies of three beams， i.e. a concrete-wood composite beam， a concrete beam with a corrugated steel web and steel-concrete composite beam. Since the proposed method is based on the two-dimensional theory， it can provide benchmarks for beam theories and error analyses.

Key words： free vibration；composite beam；mixed finite element；DQM；state space method