復合材料層合板自由振動的譜切比雪夫解法

郭琛琛， 劉 濤， 王青山， 秦 斌

(1. 中南大學 輕合金研究院， 長沙 410083； 2. 中南大學 高性能復雜制造國家重點實驗室， 長沙 410083；3. 中南大學 交通運輸工程學院， 長沙 410075)

復合材料層合板由于其較高的比強度和比剛度等力學特性，被廣泛應用于土木工程、海洋工程和航空航天等諸多工業領域。這類構件通常面臨各種極端惡劣的工況，并且可能受到形式各異的邊界載荷約束。因此，不同邊界條件下復合材料層合板的振動特性預測一直是研究人員關注的焦點。

以往，層合板的振動問題求解主要是基于理論解[1-2]。然而，由于數學問題的復雜性以及對實際問題進行建模的高要求，這些常規方法難以同時滿足計算效率和計算精度要求。隨著研究的不斷深入，人們提出了一系列的數值求解方法，如Galerkin法[3]、里茲法[4-6]、改進傅里葉級數法[7]、無網格法[8-15]等。Shi等用Galerkin法對四邊固定的任意層合板進行了自由振動分析。Ye等用Raleigh-Ritz法得到了具有一般邊界約束和內線支承的復合材料中厚層板自由振動問題的修正傅立葉解。Shi等用改進的傅里葉級數方法(IFSM)研究了彈性地基上均勻支承和多點支承的中厚復合材料層合矩形板的自由和強迫振動特性。Xiang等提出了一種基于樣條徑向基函數和高階剪切變形理論的無網格法來分析固支層合板的自由振動問題。張科使用基于移動最小二乘(MLS-L)近似函數構造的三維無網格法，對正交層合板的振動問題進行了分析。陳莘莘等提出了無網格自然鄰接點Pertov-Galerkin法分析了復合材料層合板的自由振動。使用插值型重構核粒子無網格法，李情等分析了復合材料層合板的自由振動特性。通過對上述文獻分析可知，目前對于分析層合板振動特性的數值方法種類豐富，不同方法之間的區別主要是基/試函數和求解方法的不同，但是目前絕大多數的數值方法缺陷是在計算層合板振動問題建模過程復雜，并且求解過程一般采用傳統的數值積分方法，影響計算效率。

譜切比雪夫法(ST)采用具有指數收斂特性、優越的數值穩定性和高計算精度[16-17]的切比雪夫多項式作為基函數，來分析結構的動力學問題。在文獻[18]中，Yagci等首次系統闡述了譜切比雪夫法，其核心思想在于：一是采用準確的導數和內積矩陣[19]來求解方程中導數和積分運算，以提高算法的計算效率；二是采用Gauss-Lobatto節點[20]對函數進行離散化處理以保證計算結果的準確性。基于以上優點，譜切比雪夫法已成功應用于各種結構[21-26]的動力學特性研究中，例如Bediz使用二維譜切比雪夫(2D-ST)法預測了任意幾何形狀厚板在不同邊界條件下的動力特性，然而對于復合材料層合板的自由振動問題暫無涉及。

本文將譜切比雪夫法應用于復合材料層合板的自由振動特性研究。基于一階剪切變形理論，引入邊界彈簧技術實現不同邊界條件的模擬，得到層合板的能量方程，應用二維譜切比雪夫法求解能量方程獲得結構的振動分析模型。通過數值分析算例，與其他方法的計算結果進行比較，證明了此方法的有效性，還在此基礎上研究了彈性模量比和鋪設角對復合材料層合板自由振動特性的影響。

1 理論公式

1.1 模型描述

如圖1所示為復合材料層合板的幾何形狀和坐標系。直角坐標系o-xyz位于結構的幾何中面上，層合板沿x,y和z方向的長度、寬度和厚度分別是a、b和h；α是層纖維方向和x軸之間的鋪設角度；U、V、W分別表示層合板中任意一點在x、y、z方向上的位移。在本文中，采用邊界彈簧技術實現不同類型的邊界條件，引入了三組線性約束彈簧(ku、kv、kw)和兩組轉角約束彈簧(kφ、kθ)，kt(t=u,v,w,φ,θ)的值可以從零變化到無窮大，它們的不同組合可以表示各種邊界約束條件。例如，kt= 0(t=u,v,w,φ,θ)表示自由邊界條件(F)，而kt= ∞(t=u,v,w,φ,θ)表示固支邊界條件(C)。

圖1 復合材料層合板的幾何模型與坐標系Fig.1 Geometric model and coordinate system of composite laminated plates

1.2 能量方程

應用一階剪切變形理論，復合材料層合板中任意一點的位移場的表達式為

(1)

式中：廣義位移向量u= (uvwφθ)T，u,v和w表示在層合板幾何中面上沿著x,y和z方向的位移分量；φ和θ分別表示分別為x-z平面和y-z平面的轉角位移。

由式(1)得應變-位移線性關系(幾何方程)為

(2)

(3)

根據廣義胡克定律，層合板第k層的應力-應變關系(物理方程)為

(4)

力和力矩的合力是通過積分板厚度上的應力獲得的，通過引入剪切校正因子κ，最終可以得出廣義力與應變之間的關系。其關系式表示如下

(5)

式中：N={NxNyNxy}T,M={MxMyMxy}T,QS={QxQy}T分別為合力，合力矩和等效剪力向量；A為拉伸剛度矩陣，B為耦合剛度矩陣，D為彎曲剛度矩陣，AS為剪切剛度矩陣，它們的表達式如下所示

(6)

因此，存儲在板中的應變能US為

(7)

式中，R為5×5線性對稱微分算子矩陣。其表達式為

(8)

動能T的計算表達式為

(9)

利用方程(1)定義的位移場，沿z方向求積分，動能項可改寫為

(10)

G為5×5對稱質量矩陣，其表達式為

(11)

式中：(I0,I1,I2)為慣性矩；ρk為第k層材料密度。Usprings是具有彈性均勻邊界條件的板的四邊的勢能，它是用五種彈簧均勻分布在四個邊來模擬。

(12)

式中,k=diag(kukvkwkφkθ)。

因此層合板的拉格朗日能量方程可寫為

L=T-US-Usprings

(13)

1.3 譜切比雪夫近似原理

切比雪夫多項式為遞歸和正交多項式，在區間x∈[-1,1]上，其表達式為

Cl(x)=cos(larc cos(x)),l=0,1,2,…,

(14)

應用區間變換，對于任意區間[l1,l2]的函數y(x)，可以使用有限項的切比雪夫多項式展開表示為

(15)

在離散化系統(即確定采樣點)時，采用Gauss-Lobatto采樣方法[20]

l=0,1,…,M-1

(16)

擴展系數al與采樣點yl=y(xl) 之間存在一對一的映射關系，采樣函數y={y0,y1, …,yM-1}的向量與(截斷的)切比雪夫展開系數向量a={a0,a1, …,aM-1}之間的關系可以寫成

y=ΓBa

a=ΓFy

(17)

(18)

b=a

(19)

(20)

其中s= (l2-l1)/2為區間縮放系數。

應用切比雪夫多項式的積分性質，y(x)的積分表達式為

(21)

為了得到內積矩陣，令h(x) =f(x)g(x)，假設fd=diag[f1f2…fM]，跟據式(17)和(21)有

(22)

令vTΓF= [V1…VM]，故式(22)可改寫為

(23)

對于二維問題，函數f(x,y)將產生第二階張量，其組成元素表達式為

flm=f(xl,ym)

(24)

式中，l= 1, …,Mx和m= 1, …,My是張量的指數。fkl為組成Mx×My二維采樣函數矩陣F的元素值。使用以下映射算法獲得采樣二維函數的矢量形式，為

fc=flm,c=(l-1)My+m

(25)

應用二維譜切比雪夫法，在坐標(x,y)中的每個位移點都用雙切比雪夫多項式展開表示為

(26)

式中:alm是雙展開系數;Mx和My是在每個方向上使用的切比雪夫多項式數。將式(25)張量到矢量的映射應用于展開系數，式(26)可改寫為

(27)

二維問題的展開系數與采樣函數值的關系為

f=ΨBa

a=ΨFf

(28)

(29)

對于微分運算，需要獲得二維問題中關于x和y方向的導數矩陣Dx和Dy。按照文獻[25]的方法來擴展導數矩陣，就可以分別得到相對于x和y的擴展導數矩陣Dx和Dy。使用向前和向后變換矩陣以及擴展的導數矩陣，有如下定義：

fx=ΨBbx=ΨBDxa

fy=ΨBby=ΨBDya

(30)

對于面積積分，使用內積矩陣方法

?Af(x,y)g(x,y)dxdy=fTVg

(31)

V為內積矩陣，其表達式為

V=x?y

(32)

1.4 振動求解

應用二維譜切比雪夫法求解式(13)可得到

(33)

式中，M,KU和KS分別代表質量矩陣。應變能剛度矩陣和邊界彈簧剛度矩陣，其具體表達式如下

(34)

t=u,v,w,φ,θ

Vx=x?INy×Ny,Vy=INx×Nx?y

(35)

I為單位矩陣。

對于式(33)，采用弱形式的求解方法得到對未知變量求得偏導數方程

(36)

將式(36)以矩陣的形式表達，則結構的振動特性方程可以寫為以下形式

(KU+KS-Mω2)A=0

(37)

式中:A代表未知系數的矢量;ω表示結構的固有頻率。通過直接求解式(37)，可以得到板的頻率特性和相應的特征向量。

2 數值分析與算例

基于前面得到的自由振動特征求解方程，本部分將對復合材料層合板的自由振動特性進行分析。為簡化研究，如無特別說明，在接下來的計算分析中，層合板的幾何參數為：a= 1 m,b= 1 m,h= 0.1 m；材料參數為：E2= 10 GPa、E1= 25E2、G12= 0.5E2、G13=G12、G23=0.2E2、μ12= 0.25、μ12=μ21、ρ= 1 500 kg/m3，κ= 5/6。如前所述，對于不同的邊界可以設置與之對應的彈簧約束剛度實現，符號C、S及F分別用來表示固支、簡支、自由邊界。特別地，本文所涉及的邊界組合按逆時針方向設定。此外，本文方法具有簡便性和通用性，當結構的幾何參數和邊界條件發生變化時，只需要修改相關參數即可，并不用重新推導，從而可以大大簡化計算過程。此外，本文中采用式(38)對層合板的固有頻率進行無量綱化。

(38)

2.1 數值驗證

在本節中，通過與文獻中的結果進行比較，對二維譜切比雪夫法進行了數值驗證。首先，表1為改變x和y方向上所需的多項式數Mx×My，四端固支的方形層合板的前五階無量綱頻率。層合材料的纖維鋪設角為[30°/-30°]，板厚h=0.05 m。從表1可以看出，隨著多項式數的遞增，各階頻率依次從低階到高階收斂至穩定值, 當Mx×My=16×16時，所有無量綱頻率均已達到收斂值。與Jacobi-Ritz法進行對比，兩種方法均是基于一階剪切變形理論，所以兩者前五階無量綱頻率計算結果高度吻合，相對誤差在0.001%以內。此外，本文的譜切比雪夫法達到收斂值的矩陣維數為1 280×1 280 (5MxMy× 5MxMy)，而文獻[6]的矩陣維數為5 120×5 120。本文只需要少量的多項式數即可獲得足夠準確的結果，充分說明了本文所提出的譜切比雪夫法具有收斂速度快，同時計算準確性好的特點。

表1 采用不同多項式數 Mx×My的復合材料層合板的前五階無量綱頻率對比Tab.1 Comparison of the first five-order dimensionless frequency of composite laminated plates with different Mx×My

表2所示為在固支邊界條件下，[0/α]和[α/-α]兩種鋪設方式下的復合材料層合板的前三階無量綱頻率，密度為ρ=1 000 kg/m3。將計算結果分別和Galerkin法和局部移動Kriging無網格法的計算結果進行了比較,從表2可以看出，三種方法的計算結果可以很好地匹配，微小的誤差來源于主要是因為不同數值方法，不同方法之間的計算原理在本質上存在一定差異。

表2 [0/α]和[α/-α]兩種鋪設方式下復合材料層合板的前三階無量綱頻率對比Tab.2 Comparison of the first three-order dimensionless frequencies of composite laminated plates with the two layer types [0/α] and [α/-α]

最后，由于本文方法是為了求解層合板在任意邊界條件下的自由振動問題，因此有必要驗證其在各種邊界條件下的有效性和計算精度，相關對比數據結果如表3所示，參考值為修正傅里葉解計算結果。材料密度ρ=1 000 kg/m3，纖維鋪設角為[45°/-45°]。在該算例中，考慮了四種邊界條件，分別是C-C-C-C、S-S-S-S、F-C-F-C、F-S-F-S，和兩種寬厚比，即h/a=0.05和0.2。兩種方法同時基于一階剪切變形理論，只是基函數和求解形式不同，所以計算結果非常接近。結果表明，本文方法適用于任意邊界條件。

表3 不同邊界條件下[45°/-45°]復合材料層合板的前三階無量綱頻率對比Tab.3 Comparison of the first three-order dimensionless frequencies of [45°/-45°] composite laminated plates with different boundary conditions

2.2 參數化研究

根據引言調研的文獻可知，在工程實踐中，板結構的材料屬性往往是多種多樣的，而纖維的敷設角度和纖維自身的彈性模量取值對于層合板的振動特性而言至關重要，這也直接反應在復合材料的纖維類型選擇和加工工藝上。因此，詳細研究纖維的敷設角度和纖維自身的彈性模量取值對復合材料層合板的自由振動特性的影響對指導實際工程設計和應用有著重要意義。本節以四端固支，鋪設方式為[α/-α]的層合板為例，研究了彈性模量比E1/E2和鋪設角α對復合材料層合板無量綱頻率的影響。層合板的其余結構參數參考默認。E1/E2的變化區間為[1,60]，α從0°遞增至180°，提取層合板的1階、10階和100階無量綱頻率，最終計算結果如圖2所示。需要注意的是，當計算100階固有頻率時，所采用的多項式數為Mx×My=40×40。從圖2可以看出，隨著E1/E2的增大，1階、10階和100階的無量綱頻率也隨之遞增，E1/E2取最大值60時，頻率也達到了最大值；對于鋪設角α，1階、10階和100階的無量綱頻率變化規律關于α=90°對稱，并且模態階次不同，頻率變化規律也會略有不同，1階頻率在0°、90°和180°附近達到最大值，10階和100階頻率則在區間[30°,60°]和[120°,150°]內取最大值。

1階

10階

100階圖2 彈性模量比E1/ E2和鋪設角α對層合板無量綱頻率的影響Fig.2 Effects of elastic modulus ratio E1/ E2 and laying angle α on non-dimensional natural frequencies of laminated plates

3 結 論

本文在一階剪切變形理論的基礎上，應用二維譜切比雪夫法建立了任意邊界條件下復合材料層合板的自由振動分析模型。通過數值對比分析和參數化研究，得出以下結論：

(1) 所采用的譜切比雪夫法具有收斂速度快，計算準確性好的特點，可以有效分析任意邊界條件下復合材料層合板的自由振動問題。

(2) 彈性模量比E1/E2對復合材料層合板的自由振動頻率有重要影響。隨著彈性模量比E1/E2的增加，結構的無量綱固有頻率也會顯著增加。

(3) 對于鋪設角α，復合材料層合板的無量綱頻率關于α=90°對稱分布，并且與模態階次有一定的關聯性。