李 情，陳莘莘

（華東交通大學(xué) 土木建筑學(xué)院，南昌 330013）

引言

復(fù)合材料層合結(jié)構(gòu)因其比強(qiáng)度高、耐高溫耐腐蝕、比模量高、可設(shè)計(jì)性好等優(yōu)點(diǎn)，廣泛用于航空航天、汽車、船舶海洋等工程領(lǐng)域，對(duì)其振動(dòng)特性的研究具有重要的工程意義[1-3].現(xiàn)有層合板理論的主要差別是對(duì)于橫向剪切形變的處理，經(jīng)典層合理論基于Kirchhoff-Love 假設(shè)，忽略了層間的剪切和擠壓，只適用于薄板；一階剪切變形理論(FSDT)考慮橫向剪切形變，引入了剪切修正系數(shù)，對(duì)中厚板適用結(jié)果良好；高階剪切變形板理論(HSDT)引入剪切應(yīng)力沿板厚方向的非線性變化，避免了FSDT 剪切修正系數(shù)的引入，但是需要有限元方程C1連續(xù)，在實(shí)際計(jì)算過(guò)程中較難實(shí)現(xiàn).諸多學(xué)者利用較為簡(jiǎn)潔的FSDT 對(duì)層合板進(jìn)行了自由振動(dòng)分析，得到了令人滿意的結(jié)果.

迄今為止，工程結(jié)構(gòu)分析中常見(jiàn)的數(shù)值計(jì)算方法有：有限元法、邊界元法、無(wú)網(wǎng)格法等.作為應(yīng)用最為廣泛和成熟的數(shù)值計(jì)算方法，有限元法在復(fù)合材料層合板的動(dòng)力問(wèn)題分析中應(yīng)用廣泛[4-5].然而，傳統(tǒng)的有限元法大多是基于自然坐標(biāo)的，計(jì)算時(shí)需要進(jìn)行等參變換和Jacobi 矩陣的計(jì)算，當(dāng)單元嚴(yán)重不規(guī)則或變形時(shí)，坐標(biāo)轉(zhuǎn)換的要求將顯著降低計(jì)算精度；同時(shí)，剛度矩陣存在過(guò)剛現(xiàn)象，對(duì)三角形單元等低階單元尤為突出.為了解決上述難點(diǎn)，學(xué)者們提出了不同的解決方法，Liu 等[6]將有限元法和無(wú)網(wǎng)格法中的應(yīng)變光滑技術(shù)相結(jié)合，提出了光滑有限元法(S-FEM)，展現(xiàn)了良好的計(jì)算精度、效率和穩(wěn)定性.基于光滑域構(gòu)造方式的不同，光滑有限元法包括基于單元子域的光滑有限元法(CS-FEM)[6]、基于單元節(jié)點(diǎn)的光滑有限元法(NS-FEM)[7]、基于單元邊界的光滑有限元法(ES-FEM)[8]和基于單元面的光滑有限元法(FS-FEM)[9].眾多學(xué)者將光滑有限元法拓展到彈塑性問(wèn)題[10-11]、板殼問(wèn)題[12-13]以及復(fù)合材料分析[14-15]等領(lǐng)域，并取得了良好的結(jié)果.目前Reissner-Mindlin 板殼單元依然是工程實(shí)際數(shù)值計(jì)算中常見(jiàn)的單元類型，但當(dāng)板殼厚度很薄時(shí)，常伴隨“剪切自鎖”的現(xiàn)象.為了消除自鎖，眾多學(xué)者開(kāi)展了大量的研究工作，提出了不同的解決方法，比如基于選擇性積分的板殼單元、MITC4 單元、C0 單元、EAS 單元、ANS 單元等.Bletzinger 等[16]提出了離散剪切間隙(DSG)法，構(gòu)造了DSG3 單元，利用離散單元剪切間隙得到新的剪切應(yīng)變矩陣，削弱了“自鎖”現(xiàn)象.然而，原始DSG 單元的剪切應(yīng)變會(huì)受單元節(jié)點(diǎn)編號(hào)的順序變化的影響，尤其是對(duì)于不規(guī)則網(wǎng)格和畸變網(wǎng)格，計(jì)算精度和穩(wěn)定性不高.

Nguyen-Xuan 和Nguyen-Thoi 等將光滑有限元法與DSG 法結(jié)合，提出了ES-DSG3 單元[17]、NS-DSG 單元[18]、CS-DSG 單元[12]，用于分析Reissner-Mindlin 板的問(wèn)題，得到了較好的結(jié)果.值得注意的是，在上述方法中雖然對(duì)應(yīng)變矩陣進(jìn)行了光滑處理，但是Reissner-Mindlin 板有限元方程中的荷載矢量、質(zhì)量矩陣仍需要進(jìn)行等參變換，無(wú)法完全避免坐標(biāo)映射和變換.由此，Yang 等[19]基于全局坐標(biāo)構(gòu)造了Reissner-Mindlin 板改進(jìn)的邊界光滑DSG 單元(RES-DSG3)，提出了一種非等參DSG 方法，將其與符號(hào)積分和光滑技術(shù)相結(jié)合，對(duì)矩陣進(jìn)行了光滑處理，同時(shí)避免了坐標(biāo)映射，并將此法用于分析Reissner-Mindlin 板，取得了良好的計(jì)算結(jié)果.本文借助于RES-DSG3 法，詳細(xì)推導(dǎo)了復(fù)合材料層合板自由振動(dòng)的控制方程，并編制了相應(yīng)的MATLAB 程序，通過(guò)典型算例的計(jì)算，驗(yàn)證了本文方法的可行性和有效性.

1 層合板自由振動(dòng)的有限元控制方程

基于FSDT，復(fù)合材料層合板的位移場(chǎng)可以描述為

式中，u,v,w表示層合板內(nèi)任意點(diǎn)x,y,z變形后的位移，u0,v0,w0表示層合板中面上相應(yīng)點(diǎn)的位移，βx,βy為中面法線關(guān)于y軸和x軸方向的轉(zhuǎn)角，如圖1 所示.

圖1 層合板的坐標(biāo)系Fig.1 The coordinate system for the laminate

位移應(yīng)變關(guān)系為

合力、合力矩和等效剪力向量可分別描述為

則廣義本構(gòu)關(guān)系為

式中，A為拉伸剛度矩陣，B為耦合剛度矩陣，D為彎曲剛度矩陣，As為剪切剛度矩陣，它們的元素可以分別表示為

其中

上標(biāo)(n)表示第n層，表示彈性矩陣系數(shù)，ζ 為剪切修正因子.

根據(jù)Hamiltion 原理，不考慮阻尼，可得層合板自由振動(dòng)的有限元控制方程為

式中，K為結(jié)構(gòu)的系統(tǒng)剛度矩陣，M為結(jié)構(gòu)的系統(tǒng)質(zhì)量矩陣，u為系統(tǒng)的位移矢量.同時(shí)有

且

式中，Ni為單元的形函數(shù).

單元的質(zhì)量矩陣為

式中，m為質(zhì)量密度矩陣，其表達(dá)式為

對(duì)第K層材料有

2 RES-DSG3 法[17]

2.1 基于全局坐標(biāo)系的非等參DSG3 法

有限元方程中節(jié)點(diǎn)i相關(guān)的剪切間隙可描述為全局卡氏坐標(biāo)系下的標(biāo)量形式：

同時(shí)

全局坐標(biāo)系下，三角形單元的形函數(shù)為

式中，xi，yi為三角形單元節(jié)點(diǎn)的全局坐標(biāo)，aij(i=1,2,3;j=1,2,3)是與式(21)中逆矩陣對(duì)應(yīng)的矩陣元素.

將形函數(shù)代入方程(17)、(18)，得到

式中

則剪切應(yīng)變可表示為

基于DSG 方程，復(fù)合材料層合板的單元剪切應(yīng)變可寫為矩陣形式，如下：

式中

2.2 邊界光滑技術(shù)

對(duì)于平面問(wèn)題，基于符號(hào)積分和Gauss 散度定理有

式中，Г為積分域邊界.

假設(shè)積分域包含k條邊界，且為光滑連續(xù)，則式(37)、(38)可寫為沿積分域邊界的線積分之和：

應(yīng)用近似積分技術(shù)，光滑域內(nèi)點(diǎn)xC處的任意函數(shù)可表示為

式中，?(x-xC)為光滑函數(shù)，定義為

將式(42)代入式(41)中，得

將式(43)代入式(39)、(40)中，得

3 基于REG-DSG3 的有限元方程

假設(shè)求解域Ω離散為Ne個(gè)單元，求解域這些單元網(wǎng)格共有Neg條邊，將每條邊的兩個(gè)端點(diǎn)

和這條邊相鄰的兩個(gè)三角形單元的中心相連接，這樣就在三角形單元的基礎(chǔ)上形成了Ns個(gè)基于邊的光滑域，求解域光滑域的數(shù)目和三角形單元邊的數(shù)目相等，即Neg=Ns，如圖2 所示.

圖2 三角形單元和基于邊界的光滑域Fig.2 Triangular elements and smoothing domains associated with edges

利用邊界光滑技術(shù)和DSG 單元，層合板自由振動(dòng)控制方程中的應(yīng)變矩陣可表示為

其中

式中，k表示每個(gè)光滑域的邊界段數(shù)，對(duì)于求解域內(nèi)邊界處的邊，k=3，對(duì)于求解域內(nèi)部的邊，k=4；Ni是與邊界段m相關(guān)的三角形單元節(jié)點(diǎn)i的形函數(shù).

其中

由此，光滑剛度矩陣為

同理，光滑質(zhì)量矩陣可通過(guò)以下轉(zhuǎn)換得到：

通過(guò)Gauss 積分，得

由此，層合板有限元控制方程中的矩陣都已沿光滑域的邊界段進(jìn)行光滑處理和計(jì)算，在計(jì)算過(guò)程中不需要坐標(biāo)映射和Jacobi 矩陣的計(jì)算.

4 數(shù)值算例

下面通過(guò)數(shù)值算例分析復(fù)合材料層合板的自由振動(dòng)問(wèn)題.如無(wú)特殊說(shuō)明，層合板各層具有相同的厚度，且由均勻線彈性復(fù)合材料構(gòu)成，無(wú)量綱固有頻率為.

4.1 算例1

首先以一長(zhǎng)度為a、厚度為h的四邊簡(jiǎn)支對(duì)稱層合方板(0°/90°/90°/0°)為研究對(duì)象，探討RES-DSG3 法在求解復(fù)合材料動(dòng)力學(xué)問(wèn)題時(shí)的收斂性和有效性.對(duì)于規(guī)則網(wǎng)格，α=0；不規(guī)則網(wǎng)格，α=0.4，其中α為定義在區(qū)間0 ～ 0.5 的不規(guī)則系數(shù)[20].相關(guān)的材料參數(shù)為E1/E2=10,20,30，G12=G13=0.6E2，G23=0.5E2，v12=0.25，ρ=1，剪切修正因子ζ=5/6.

由表1 可見(jiàn)，不同節(jié)點(diǎn)分布和不同彈性模量比時(shí)的四層層合板的一階無(wú)量綱固有頻率與文獻(xiàn)[21]的一階剪切理論的精確解很接近，具有較高的計(jì)算精度，表明了本文方法的有效性和收斂性.同時(shí)，網(wǎng)格“質(zhì)量”對(duì)計(jì)算精度的影響很小，即使在極不規(guī)則的網(wǎng)格下也能獲取穩(wěn)定可靠的計(jì)算結(jié)果，進(jìn)一步拓寬了該單元的使用范圍.

表1 四層復(fù)合材料簡(jiǎn)支層合板的一階無(wú)量綱固有頻率Table 1 The non-dimensional 1st fundamental frequency of the simply supported 4-layer laminated composite plate

4.2 算例2

考慮三層對(duì)稱層合方板(0°/90°/0°)，邊長(zhǎng)為a，厚度為h，材料參數(shù)為E1/E2=40.采用17 × 17 網(wǎng)格離散求解域，對(duì)不同邊界條件和不同邊厚比條件下的板進(jìn)行了數(shù)值計(jì)算.

表2 列出了通過(guò)本文方法計(jì)算的四邊簡(jiǎn)支(SS)，四邊固支(CC)，兩對(duì)邊簡(jiǎn)支、另兩對(duì)邊固支(SC)三種邊界條件和不同邊厚比的三層層合方板的一階固有頻率，并分別與文獻(xiàn)[21-22]基于FSDT 的復(fù)合二次徑向基函數(shù)法計(jì)算的結(jié)果進(jìn)行了對(duì)比，具有較好的一致性.同時(shí)，表3 列出了三種邊界條件下的前三階固有頻率，可以看出本文方法計(jì)算的無(wú)量綱固有頻率與文獻(xiàn)[23]基于高階剪切變形理論計(jì)算的結(jié)果也很接近，驗(yàn)證了本文方法的有效性.

表2 三層復(fù)合材料層合板的一階無(wú)量綱固有頻率Table 2 The non-dimensional 1st fundamental frequency of the 3-layer laminated composite plate

表3 三層復(fù)合材料層合板前三階無(wú)量綱固有頻率（a/h=10）Table 3 The non-dimensional 1st 3 fundamental frequencies of the 3-layer laminated composite plate(a/h=10)

5 結(jié)論

本文將基于全局坐標(biāo)的重構(gòu)邊界光滑DSG 單元用于復(fù)合材料層合板的自由振動(dòng)問(wèn)題分析.與原始的DSG 單元相比，該方法不需要對(duì)系統(tǒng)方程進(jìn)行坐標(biāo)映射，同時(shí)結(jié)合光滑技術(shù)將有限元矩陣的域積分簡(jiǎn)化為沿平滑單元邊界的線積分.從數(shù)值計(jì)算結(jié)果來(lái)看，RES-DSG3 法在復(fù)合材料層合板自由振動(dòng)分析時(shí)表現(xiàn)出良好的性能，具有較高的計(jì)算精度；即使是不規(guī)則網(wǎng)格，也能獲得較好的結(jié)果，是一種有效可行的方法.此外，該方法中的光滑技術(shù)也可以擴(kuò)展到其他問(wèn)題的有限元框架內(nèi).