框架結構中矩形樓板動力響應求解的譜單元法

曹容寧，馬 蒙，孫曉靜，劉維寧

(北京交通大學 土木建筑工程學院，北京 100044)

環境振動影響下的建筑樓板動力響應計算對建筑室內振動和二次噪聲的預測評價起到關鍵作用。柱、梁、板是框架結構中振動傳播的主要路徑，對于此類復雜結構多采用有限元法(finite element method，FEM)建模分析。但軌道交通引起的室內振動和二次噪聲問題，關心的頻率可達到250 Hz[1]，由此帶來的過密的網格劃分不僅會造成計算效率低，還會導致計算不穩定。解析法是精確求解板動力響應的常用方法，但只能考慮一些具有特殊幾何特性和典型邊界條件的情況[2]；雙向梁法可適用于矩形板四邊為任意典型邊界條件的情況[3](如簡支、固支或自由邊界)，而無法用于具有非均勻截面特性或非典型邊界條件下的板結構。

SEM是將位移響應表示為有限個簡諧波的疊加，進而通過計算形函數和動剛度矩陣完成結構動力響應的求解。Doyle利用該方法分析波在桿、梁等結構中的傳播。在此研究基礎上，Lee等[9-11]用SEM求解周期性格柵結構、分層復合梁結構、管道結構等復雜結構。對于類一維板問題，Birgersson等[12]利用任意邊界的一組對邊，采用超譜單元法實現了矩形板結構的單方向耦合。Park等[13]采用邊界分離法和SEM結合，實現了四邊均為任意邊界條件的矩形板的動力計算，以及板結構之間的耦合。這種方法將矩形板分解為兩個具有對邊典型邊界的單向板，對于每個單向板仍按照上述降維的思路，將多項式插值法應用于一個方向，SEM應用于另一個方向，在頻域內計算板的形函數和動剛度矩陣。但是，由于該方法采用的節點自由度僅適用于板的變形，無法與梁、柱結構的節點自由度一一對應，因此只能實現板結構之間的相互耦合，而無法實現板結構與梁、柱結構的耦合。

為解決上述問題，本文在Park等的基礎上，針對框架結構提出了一種改進方法，其不僅可以求解四邊為任意邊界條件的樓板動力響應問題，還通過波的傳播對板的節點自由度進行改進，實現板與梁、柱的耦合，從而拓寬了SEM的應用范圍，實現了SEM在框架結構中的應用。該方法基本思路是將雙向板分解為兩個具有典型邊界條件的單向板問題，并應用SEM分別求解；更進一步地，采用最優化方法對板單元的節點自由度進行修正，使之與梁、柱單元一一對應，實現梁板柱耦合結構的計算。

1 矩形薄板動力響應的SEM求解方法

1.1 問題的分解

一個長、寬、厚分別為2a，2b和h的均質矩形薄板，如圖1所示，設板的中面位于xOy平面上，中面的中心點位于坐標原點上；矩形板的四邊分別為L邊、R邊、D邊和U邊，坐標方程依次為：x=-a，x=a，y=-b以及y=b。由于板的平面內剛度較大，所以本文只考慮板的平面外振動，即彎曲振動。建筑結構中板的厚度通常遠小于其長度和寬度，因此薄板彎曲振動問題可以采用Kirchhoff薄板模型進行計算。Kirchhoff薄板假定板在彎曲時，中面法線保持直線，不發生伸縮，且一直垂直于板的中面。當采用SEM求解時，板的形函數通過控制方程求出，這樣的形函數考慮了板的動力特性，較FEM多項式插值形函數更準確，計算精度受網格尺寸影響也不大，故無需劃分過密的網格單元和節點。

圖1 矩形薄板示意圖

圖2 矩形板分解

1.2 A子問題、B子問題求解

A問題的研究對象為D邊、U邊為自由邊界，L邊、R邊為任意邊界的單向板。單向板的形函數由有限條單元法[14]和SEM聯合推導而得，為了計算板的動力響應問題，還需要得到動剛度矩陣。

根據有限條單元法，沿y方向將板均勻劃分成ny個條形單元，如圖3所示，每一個條形單元內部的位移通過y方向的多項式插值函數和兩交線的位移函數表示。各個條形單元在交線上保持位移連續，依此將各個條形單元集成，單向板A的位移可以表示為

圖3 A問題有限條單元劃分

uAz(x,y)=ZA(y)·WA(x)

(1)

式中：WA為各個交線上的位移函數向量；ZA為與交線位移函數對應的y方向插值函數向量。

然后，基于SEM的思想，交線上的位移函數可以視為x方向波的疊加，即

(2)

式中：kAi為x方向的波數；φAi為y方向的截面振型向量；aAi為波的幅值。

將位移函數式(1)代入板的Hamilton方程[15]中，結合Parseval定理[16]，可以得到頻域內的特征方程。從而可以求解出兩兩互為相反數的特征值，即沿x正負方向傳播的波數kAi(ω)。每個特征值對應的特征向量即y方向的截面振型φAi。

在不同頻率點處求解上述特征值問題，繪制0～250 Hz的頻散曲線，如圖4所示(板的材料參數見3.1節)。實部為正表示波沿x軸正方向傳播，實部為負表示波沿x軸負方向傳播。在0～250 Hz內，x方向主要包括4對正負方向傳播的波，每一種波對應的y方向的截面振型如圖4所示。在截止頻率為60 Hz左右時，第3階波出現；在截止頻率為180 Hz左右時，第4階波出現。由此可知，隨著頻率的升高，參與傳播的波逐漸增多，且與之對應的y方向振型也趨于復雜；另外，每一條頻散曲線的波數大小均隨頻率增加而增加，說明隨著頻率增大，傳播波的波長逐漸變小。

圖4 頻散曲線及對應振型

為了求解波的幅值aAi，需借助各個交線兩端點處的位移值，即單向板L邊、R邊兩邊節點上的邊界條件。利用如下含有未知常數aAi的交線位移函數式表示L邊、R邊邊界節點的位移

(3)

為了分析框架結構中板的動力特性，需要將板結構與梁、柱結構耦合。本文所采用的節點自由度體系為傳統的“三平動、三轉動”自由度，即{uxuyuzθxθyθz}T。由于不考慮板的平面內變形，板的節點自由度設定為{uzθxθy}T。在前文對單向板問題進行求解時，每個節點需要提供4個自由度作為邊界條件，為了使板單元與梁柱單元的自由度一一對應，需要對板單元的節點自由度θAxy進行縮減，為此，對式(3)進行如下修正

(4)

為了解決欠定方程組問題，本文參考Birgersson等的研究，采用計權最小二乘法求解符合等式(4)的最優解。式(4)中的位移函數WA表示各個波的疊加，在進行最小二乘法求解之前，應對各個疊加波進行計權。計權的原則是：選擇對位移響應貢獻大的波進行疊加，即降低波數較大的波的權重。這樣做的原因有兩點：第一，波數較大的波波長短，沿x方向衰減很快，對位移函數的貢獻不大；第二，波數大的波y方向的截面振型復雜，y方向節點數量可能不足以描述截面振型，這會導致高階波的振型與低階波的振型向量相似，從而在進行波的疊加時代替低階波的貢獻量，使計算結果在低頻出現誤差。計權的方法是將每個波的截面振型向量進行歸一化處理，使截面振型向量的2-范數等于1/|kAi|，kAi為與各個振型向量對應的波數。這種方法實現了將各個疊加波按照波數的大小進行計權，之后采用最小二乘法求解式(4)的欠定方程組，即可得到各個疊加波的幅值aAi(注意這里易偶發病態現象)。結合式(1)和式(2)，可以得到位移函數

uA=NAdA

(5)

式中，NA為A問題的形函數。

將位移函數(5)代入板的Hamilton方程中，即可以的得到節點位移與節點力向量的關系

SAdA=fA

(6)

式中，SA為A問題的動剛度矩陣。

B問題的單向板中，L邊、R邊為固定邊界，D邊、U邊以節點位移為邊界條件。其與A問題的差異在于：①“典型對邊”由自由邊界變為固定邊界；②“任意邊界”的節點位移等于原問題與A問題的節點位移之差。沿x方向將板均勻劃分nx個條形單元，具體推導過程與A問題相似，如圖5所示。

圖5 B問題有限條單元劃分

1.3 原問題求解

由于原問題的位移場為A問題、B問題的疊加，板結構彎曲振動的位移函數則可以表示為

uz(x,y)=uAz(x,y)+uBz(x,y)=

NAdA+NBdB

(7)

如前所述，原問題與A問題、B問題的節點位移關系為

(8)

uz=Nd

(9)

式中，N為形函數。

為了能夠形象地描述矩形板四邊節點自由度對應的形函數，這里以圖6中的板構件為例，繪制幾個典型節點自由度上的形函數(材料和尺寸參數與第3章中的板結構相同)。繪制靜態(0)和動態(200 Hz)下的形函數，如圖7所示。可以看出，200 Hz下形函數較0下的形函數復雜，這是因為在形函數的計算過程中，考慮了板的動力特性。這也說明，與只能表示板的幾何特性的多項式插值函數不同，本文方法的形函數同時考慮了板結構的幾何特性和動力特性，可以用來表示較高頻率下結構的復雜振型。圖7(e)～圖7(h)中，在y=±2 m位置處，形函數在0位置上下稍有波動，這是由于板在高頻振動下，波長很短，導致即使兩相鄰節點位移均為0，但節點之間仍會出現微弱振型。這些微弱振型會使計算結果在某些頻帶產生突變，若要更準確地模擬板結構的邊界條件，則仍存在一定改進空間。

圖6 原問題單元劃分

圖7 典型節點自由度的形函數

將原問題的位移函數式(9)代入薄板彎曲振動的Hamilton方程，可以得到板的節點位移與節點力向量的關系，即動剛度的計算。

2 梁板柱的耦合

梁、柱結構與板結構具有不同的振動形式，結構間的相互連接使振動相互影響。在數值模擬中，為了模擬結構之間的耦合作用，令梁、板、柱結構在公共節點處的位移和力連續。分別計算梁、板、柱結構的動剛度矩陣，并按照FEM的集成規則集成總動剛度矩陣，即可計算梁板柱耦合結構的動力響應。

梁結構、柱結構的振動形式相同，分別為軸向振動、彎曲振動和扭轉振動，本文假設3種振動形式互不影響。應用SEM的一維求解方法，分別根據軸向振動控制方程、歐拉-伯努利梁的彎曲振動控制方程和扭轉振動控制方程，計算梁(柱)結構的動剛度矩陣，其與節點位移和節點力的關系可表示為

(10)

板結構彎曲振動的動剛度矩陣由前文計算得到，以平行于直角坐標系xOy平面的板結構為例，僅考慮彎曲振動的板結構的節點位移與節點力的關系為

(11)

本文視板的平面內變形為剛性變形，所以忽略了平面內變形的自由度ux,uy,θz。當梁、柱結構與板結構耦合時，也需將ux,uy,θz自由度對應的廣義位移設為0，即在進行動力計算時，將式(10)中與ux,uy,θz自由度對應的矩陣元素剔除。縮聚后，式(10)變為

(12)

(13)

3 算例分析

3.1 算例概況

本節基于板的A問題、B問題、原問題，以及梁板柱耦合結構，設置4種工況進行分析。板的中面均位于坐標系xOy平面，坐標原點O位于板面中點。分別采用SEM和FEM對4種工況進行模擬，比較計算結果的準確度。

在工況1中(如圖8(a)所示)，薄板的一組對邊為自由約束、另一組對邊為簡支約束。這樣設置邊界條件的目的是為了形成對邊簡支對邊自由邊界的板，進而可以通過解析方法計算各階自振頻率，作為評價數值方法準確度的依據。在A1點施加繞y軸逆時針方向的彎矩MA1=1 000 N·m/Hz，并計算A1點繞y軸轉角θA1。采用本文改進的SEM計算時，在簡支約束的兩邊分別均勻劃分5個節點；FEM模型則劃分4×4個單元共25個節點。

在工況2中(如圖8(b)所示)，薄板的一組對邊為固定約束、另一組對邊為簡支約束，形成了對邊簡支對邊固定邊界的板。在A2點施加繞x軸逆時針方向的彎矩MA2=1 000 N·m/Hz，并計算A2點繞x軸轉角θA2。采用SEM計算時，在簡支約束的兩邊分別均勻劃分3個節點；采用FEM模擬時，仍然劃分4×4個單元，共25個節點。

在工況3中(如圖8(c)所示)，固定板的4個角點，在A3點施加z方向荷載PA3=1 000 N/Hz，計算A3點的z方向位移響應。在SEM模型中，在四邊均勻劃分16個節點，形成一個16節點的單元。在FEM模型中，分別給出4種單元劃分方法：FEM 1劃分4×4個單元，共25個節點；FEM 2劃分6×6個單元，共49個節點；FEM 3劃分10×10個單元，共121個節點。

圖8 板結構概況

在工況4中(如圖9所示)，以建筑物中一個基本房間為例(不考慮墻體)，在板結構四邊耦合梁結構，在4個角點與柱結構耦合。在4個柱腳處施加固定約束，在A4點施加z方向荷載PA4=1 000 N/Hz，并計算O點的z方向位移響應。在SEM模型中，板的四邊均勻劃分16個節點，每根梁均勻劃分5個節點并與板在節點處耦合，每根柱子僅在兩端劃分節點。FEM的節點劃分見圖9(b)、圖9(c)。

圖9 梁板柱結構節點劃分

4種工況中的板均為邊長Lp=4 m，厚度h=0.2 m的正方形薄板。梁、柱結構的長度Lb=4 m，截面尺寸為0.5×0.5 m2。梁、板、柱結構的材料參數相同，彈性模量E=35 GPa；泊松比為ν=0.3；密度ρ=2 500 kg/m3，不考慮阻尼影響。

3.2 結果分析

工況1的計算結果如圖10所示，在1～250 Hz內有5階自振頻率。通過Levy方法計算自振頻率的解析解，在圖10中以灰色虛線表示。兩種數值方法與解析解對比，得到的自振頻率誤差百分比，如表 1所示。由表1可以看出，FEM計算誤差范圍為0.5%～4.1%，SEM誤差范圍為0～1.5%，總體上，前者誤差大于后者，最高約達9倍。且隨著階數的增加，兩者誤差相差越來越大。因此，SEM的計算結果更接近解析解。從圖10中可以看出，在1～70 Hz頻段內，兩種數值計算結果幾乎一致，從大約70 Hz開始，FEM計算結果與SEM結果分離，且隨著頻率增高，兩種結果偏離逐漸明顯。這是因為從70 Hz左右頻率處開始，單向板中的第3種波數的波開始傳播(見圖4)，其波長較短且結構振型復雜。而對于70 Hz以上的振動，FEM的單元尺寸已不足以擬合出板的振型，因此計算誤差越來越大。而由于SEM的形函數考慮了板的控制方程，所以即使只在單向板的一組對邊劃分節點，無需在板的內部劃分節點，也可以在高頻范圍得到精確度較高的計算結果。

圖10 工況1計算結果

表1 轉角自振頻率誤差

工況2的計算結果如圖11所示，在1～250 Hz內有三階自振頻率出現，與解析解的自振頻率(灰色虛線)相比，數值解的誤差百分比如表2所示。FEM的誤差約為SEM誤差的3～5倍，說明后者更接近解析解。與A問題的結果類似，兩種數值方法的頻譜在低頻保持一致，隨著頻率的增高，差異開始明顯。

圖11 工況2計算結果

表2 自振頻率誤差

工況3的計算結果如圖12所示，由于該工況很難得到解析解，本文認為當FEM模型的網格劃分非常密時，可以達到與解析解相近的精度。為了驗證，對網格尺寸為0.3～1.0 m(節點數為16～225)的FEM模型進行收斂性分析。由于FEM網格劃分不夠密引起的誤差主要出現在高頻范圍，因此以較高階自振頻率作為收斂性分析的評價量，本工況采用第6階自振頻率作為評價量。圖12(b)展示了不同節點數量的FEM模型的高階自振頻率，隨著節點數的增加，自振頻率值收斂，可以達到與解析解相近的精度。

圖12 工況3計算結果

從圖12(a)中可以看出，FEM 3與SEM結果最為相近，FEM 2次之，FEM 1結果較SEM結果偏差較大。這說明隨著FEM模型劃分單元數量的增多，其計算結果與SEM結果越來越相近，尤其是在較高頻率范圍內，這種規律更加明顯。這是因為SEM的位移函數由控制方程得到，考慮了板的動力特性，所以即使在較高頻段，仍能準確地描述板的振型；而FEM的位移函數只考慮了板的幾何特性，因此需要通過增加單元和節點的數量來提高高頻范圍計算結果的準確度。

工況4的計算結果如圖13所示，由于工況4中的耦合結構沒有解析解，因此仍采用高階自振頻率(第4階)作為評價量對FEM進行節點數收斂性分析。分別對網格尺寸為0.4～2.0 m(節點數為17～161)的FEM模型進行計算，如圖13(b)所示，若節點數過少，則準確度無法保證，隨著節點數的增加，高階自振頻率值收斂。

從圖13(a)中可以看出，當SEM模型與FEM模型劃分節點數量相近時，FEM計算結果從180 Hz左右開始準確度無法保證，只有當FEM節點數遠多于SEM時，二者的計算準確度相似。在50 Hz左右頻率處，SEM計算結果出現一個小突變，這是由于矩形板與梁只在耦合的節點處符合位移連續條件，但在節點之間位移并不連續，導致板在兩節點之間出現多余振型。

圖13 工況4計算結果

3.3 計算效率比較

在頻域內進行結構動力響應數值計算時，影響計算效率的主要因素有兩個：①位移求解時動剛度矩陣的求逆計算；②SEM需要在所有頻率點循環求解動剛度矩陣。由于對于材料參數不變的結構，SEM可以通過建立前期數據庫的方式求解動剛度矩陣以節省計算時間，所以對二者計算效率的比較主要聚焦在動剛度矩陣求逆的計算效率上。

以工況3中的SEM模型和FEM 3模型為例，從圖12(a)可以看出，兩者的計算準確度較為相近，在此基礎上比較兩個模型的動剛度矩陣。SEM模型由一個16節點的板單元組成，每個單元考慮3個自由度；FEM3模型由100個4節點板單元組成，共121個節點。圖14對兩個模型的動剛度矩陣進行了可視化繪制，可以看出，SEM模型的動剛度矩陣為48維數的滿布矩陣，其中非零元素共2 304個；FEM 3模型的動剛度矩陣為363維數的帶狀矩陣，共有8 637個非零元素。從計算效率上來說，滿布矩陣求逆的計算效率低于稀疏帶狀矩陣，但是，矩陣維數越大，非零元素越多求逆計算效率越低。因此，雖然從矩陣的稀疏程度來看，FEM在計算效率上優勢較大，但是，從矩陣的非零元素數目來看，SEM的計算效率則更高。由此可以推論，對于構造簡單且板的尺寸較小的小型結構，動剛度矩陣維數較低，非零元素數目較少，其對計算效率影響不大，則動剛度矩陣較為稀疏的FEM計算效率更高。而對于含有大面積矩形樓板的規則結構來說，動剛度矩陣的非零元素數目對計算效率的影響占主導，則SEM的計算效率更高。

圖14 可視化動剛度矩陣

4 結 論

本文針對框架結構提出了一種改進的譜單元法，該方法將矩形板拆分成兩個單向板的疊加并分別計算，其不僅可以求解四邊為任意邊界條件的樓板動力響應問題，還考慮了框架結構中樓板與梁、柱的耦合關系。

通過與自振頻率解析解、以及高密度網格的FEM頻譜響應對比，驗證了本文方法的準確性。

與傳統的FEM相比，本文方法的優勢在于：

(1)對于矩形薄板結構，只需在四邊邊界劃分節點，可以減少動剛度矩陣維數，對于大尺寸樓板結構，可以提高計算機的計算效率。

(2)對于框架結構，可實現振動響應的快速預測。由于只需在構件之間的耦合處劃分節點，無需在單元內部劃分節點，因此即使結構尺寸有所變化，也無需重新劃分節點，從而提高建模效率。

(3)SEM的計算結果準確度高于FEM，尤其是在較高頻段。①在計算對邊簡支對邊自由邊界的板時，采用FEM計算誤差最大約達SEM的9倍；②在計算對邊簡支對邊固定邊界的板時，采用FEM計算誤差最大約達SEM的5倍；③若二者的計算準確度相似，FEM所需劃分節點數遠多于SEM。

此外，本文討論的板結構的SEM在結構應用上存在一定局限性。①該方法僅適用于對矩形板的動力計算，無法計算其他形狀的板結構；②若板結構的尺寸很小，則SEM在計算效率上的優勢無法突顯。因此，SEM較適用于含有大尺寸矩形樓板的規則建筑結構。