功能梯度碳納米管增強復合材料板彎曲和模態的廣義有限差分法1)

李煜冬 傅卓佳 湯卓超

(河海大學力學與材料學院,南京 211100)

引言

隨著復合材料的發展,納米級材料已經成為諸多領域的研究熱點,其中表現出優越力學性能的一維碳納米管在工程結構界逐漸受到人們的極大關注[1-3].Shen[4]受到碳納米管增強復合材料(carbon nanotube-reinforced composite,CNTRC)和功能梯度材料(functionally graded materials,FGMs)概念的啟發,提出并設計了功能梯度碳納米管增強復合材料(functionally graded carbon nanotube-reinforced composite,FG-CNTRC)模型.研究表明[4],摻入少量碳納米管增強體到結構基體中能夠有效提高工程結構抗拉和抗壓能力.通過對碳納米管纖維采取一定的梯度形式排布,能夠在纖維含量較低條件下合理高效地提升結構的宏觀力學性能.于是當FG-CNTRC被提出后,有關功能梯度碳納米管增強復合材料梁、板和殼結構的力學行為研究吸引了國內外諸多學者的注意[5-7].

目前已經有很多的數值分析方法對其從不同的角度做出了研究,其中網格類數值離散[8-9]是理論和應用比較成熟的典型數值計算方法.另一方面,隨著無網格與粒子類算法的快速發展[10-11],其進一步擺脫了高質量網格的束縛和避免了高階網格所需要的繁瑣積分,此類算法也被引入到功能梯度碳納米管增強復合材料板結構的研究中.Lei 等[12]采用Kp-Ritz 無網格方法計算了4 種CNTRC 矩形板(uniform distribution (UD),functionally graded (FG)-V,FG-O,FG-X 分布型)在橫向均布載荷作用下的大撓度彎曲變形,并給出了簡支方板在不同碳納米管分布形式下的載荷—撓度曲線.Do 和Lee[13]基于高階剪切變形理論(higher order shear deformation theory,HSDT)的徑向點插值法(radial point interpolation method,RPIM)數值討論了4 種CNTRC 矩形板的非線性彎曲問題.

廣義有限差分法(generalized finite difference method,GFDM)[14-15]作為一種新型的強式無網格區域配點型方法,以計算域內任意一點(中心點)為研究對象,根據“最短距離”的準則在中心點附近形成該點的局部支撐域,最后基于函數的泰勒展開式和移動最小二乘法將中心點處函數值的各階偏導數表示成其支撐域節點上函數值的線性疊加.該方法不僅無需網格劃分和數值積分而且避免了全域無網格配點法通常遇到的病態稠密矩陣問題,使得這類方法具有形式簡單、易于應用和實現等優點.其中Urena 等[16-17]做出了突出貢獻,在2001 年對影響廣義有限差分法計算精度的各種因素,如點簇形狀、權函數選取、鄰近點個數等進行了系統分析.作為無網格粒子法中重要的數值計算方法,該方法在各種力學問題中均得到了廣泛應用.2014 年Fan 等[18]將廣義有限差分法用于求解穩定二維柯西反算問題以及雙調和方程反算問題等;后來也應用于求解復雜的瞬態熱傳導問題[19-20].2018 年傅卓佳等[21-22]將廣義有限差分法用于求解腫瘤熱分析及Winkler板彎曲問題;隨后又將其應用于板結構聲振耦合以及功能梯度材料板結構的彎曲問題[23-25].

目前功能梯度碳納米管增強復合材料結構中涉及到碳納米管轉向問題的研究較少,本文旨在應用廣義有限差分法數值研究考慮了碳納米管轉向的FG-CNTRC 板的彎曲及模態問題.首先給出基于一階剪切變形理論的功能梯度碳納米管增強復合材料板的廣義有限差分法離散模型.隨后通過基準算例,檢驗廣義有限差分法的計算精度與收斂性.最后數值分析和討論碳納米管中不同分布型、體積分數、碳納米管旋轉角度、寬厚比、板傾斜角度、長寬比等對FG-CNTRC 板結構彎曲和模態的影響.

1 數學模型

1.1 控制方程

考慮一個由FG-V 分布型的碳納米管和基體材料共同組成的功能梯度碳納米管增強復合材料矩形板,如圖1(a)所示,長、寬、厚分別為 a,b,l,以該矩形板中面建立 x-y 坐標系,z 為板厚度方向的坐標.如圖1(b)所示,常見的碳納米管有4 種分布形式,為UD 型、FG-V 型、FG-O 型、FG-X 型,其中UD 型為均勻分布型,其余為功能梯度分布型.另外,還有一種分布為FG-Λ 型,其與FG-V 型成對稱性分布,在本質上是沒有差異的,所以在本文中只研究以上4 種分布類型.碳納米管的體積分數 VCNT可以表示為

圖1 功能梯度碳納米管增強復合材料板的等效模型Fig.1 Equivalent model of functionally graded carbon nanotubereinforced composite plate

為了表征功能梯度碳納米管增強復合材料的宏觀力學性能,Shen[4]提出的廣義混合律模型在FGCNTRC 研究中廣泛應用,其考慮了碳納米管的尺寸和溫度依賴性,并引入碳納米管的效能參數作為FG-CNTRC 的等效模型即

式中,ηj為碳納米管的效能參數;為不同方向下材料的楊氏模量和剪切模量;V*為材料的體積分數(“CNT”代表“碳納米管纖維”,“m”代表“基體”).該模型存在一定的假設:碳納米管在基體中排列整齊,且完全分散;不考慮碳納米管的長徑比、波紋度等微觀結構特征.此后,Mora-Didastjerdi 等[26]引入額外的效能參數對模型進行修正,考慮了碳納米管的波紋度和長徑比因素的影響,其與Shen[4]提出的廣義混合律模型基本相同.

另外,有關FG-CNTRC 等效材料模型的其他一些物理參數定義如下[4]

式中,ρ 為FG-CNTRC 的等效密度,μ12和 μ21為FG-CNTRC 的等效泊松比,和 μm為碳納米管和基體材料的泊松比.

基于一階剪切變形理論[24],板的位移場U=(u,v,w)T為

式中,(uo,vo,wo) 是板內任意一點在中面上的投影沿著 (x,y,z) 方向的位移,(φx,φy)是變形后原中面法線的相應轉動,t 為時間,坐標 z 為所考慮點到板中平面的距離.

根據幾何方程[24],有應變和位移的關系式

對于FG-CNTRC 板結構,其物理方程[8]為

根據內力與應力的關系[24],有

將式(7)代入式(8)有

其中

式中,A,B 和 C 分別表示拉伸、耦合以及彎曲剛度矩陣,A1為剪切剛度矩陣,χ 為一階剪切因子,對于矩形截面一般有 χ=5/6 .

接下來考慮如圖2 所示的碳納米管轉向問題[27],對式(9)中的 ψij有如下坐標變換關系

圖2 含轉向的碳納米管平面布置Fig.2 Plane layout of carbon nanotubes with orientation

其中,Λ 為坐標旋轉變換矩陣;β 為碳納米管旋轉角度,如圖2 所示;ψ 為碳納米管橫向布置時的FGCNTRC 物理關系矩陣,有如下定義

式中,針對第三方向相關的剪切模量滿足

最后,根據Hamilton 變分原理[8],結合式(8)導出FG-CNTRC 板模型的運動微分方程

式中,q 為板平面上的橫向載荷集度 q(x,y),(I0,I1,I2)有如下定義

1.2 邊界條件

對于板結構,常見的邊界條件有以下3 種:

(1) 固支邊界(clamped,C)

(2) 簡支邊界(simply,S)

(3) 自由邊界(free,F)

2 廣義有限差分法

考慮一個二維問題計算域 Ω(x,y),如圖3 所示對任意一個離散節點 xo形成一個包含 s 點數的點簇(支撐域),對應函數值 wj的4 階泰勒展開為

圖3 計算域中的點簇Fig.3 Cluster of points in computational domain

對于點簇展開后的余項可定義殘差函數

其中,ωj表示在該“點簇”中點 xj處的加權系數.采用式(18)所示的四次樣條權函數,其具有3 階連續性,且滿足單位性、緊支性和對稱性條件[16]

為了方便,定義

通過求極值使得殘差函數 Θ(w) 最小化,對 Dw中的每一項進行變分,形成如下矩陣

其中,H 是由點與點之間的距離分量以及權函數所形成的具有對稱性的系數矩陣,“S,Y,M”表示其具有上三角矩陣對稱性,R 是距離分量、權函數以及未知點函數值所形成的右端矩陣.為保證矩陣 H 的可逆性,研究表明[18-19,22]在二維計算中,支撐域點數一般滿足兩倍的泰勒展開項數.

對式(19)的線性方程組進行求解,有

令計算域中的節點滿足控制方程式(12),邊界節點滿足邊界條件式(13)～ 式(15),可得到僅含離散點未知物理量的線性方程組,從而獲得問題的數值解.因權函數的連續性,每個節點攜帶了相鄰近節點的物理信息,使得對節點形成移動的帶權局部支撐域能夠保證全域數值解的協調性[16].

對于控制方程式(12)所形成的FG-CNTRC 板結構彎曲問題,最終可得到總體線性方程組

其中,K 為總體剛度矩陣,F 為載荷列陣,U 為待求的中面位移及轉角向量.

對于控制方程式(12)所形成的FG-CNTRC 板結構模態問題,運用模態分析法[23]得到振型微分方程的矩陣表達

式中,K 表示由廣義有限差分法形成的整體剛度矩陣,M 為由式(12)構建的整體質量矩陣,δ 表示大小為板結構離散節點數下的單位矩陣,? 表示FGCNTRC 板自由振動的自然頻率;Ψ 表示由每個節點組成的振型列陣.

需要注意的是,基于一階剪切變形理論的FGCNTRC 板模型控制方程是一個最高階為二階的偏微分方程組,而本文采用4 階泰勒展開式的廣義有限差分法對該控制方程進行離散,目的是更進一步地提高計算精度(減小余項殘差)和擴大廣義有限差分法在更多問題中的適用性.

3 數值結果與討論

本節首先以各向同性板結構彎曲的計算實例討論廣義有限差分法這一算法特性,接著運用廣義有限差分法對含碳納米管轉向的FG-CNTRC 板的靜態線性彎曲和自振模態問題進行數值分析與討論.

3.1 各向同性板結構模型驗證

以各向同性均質板模型為例,圖4(a)展示了在相同運行環境下用廣義有限差分法(GFDM)和傳統有限單元法(traditional finite element method,TFEM)計算厚跨比為0.1 的四周固支方形板彎曲問題.圖4(b)討論了GFDM 中不同鄰近點數對數值結果的影響.基于圖5 展示的混亂布點模型,與不同厚跨比下的解析解[28]進行了比較,表1 展示了不同厚跨比下GFDM 均勻布點與散亂布點的相對誤差.可以看出廣義有限差分法的數值計算在任意混亂布點下,計算精度稍有所降低,但仍能得到滿意的結果.其無需網格和繁瑣的數值積分,僅通過散亂的節點信息便可進行計算,具有更高的計算收斂性和穩定性.

表1 相對誤差Table 1 Relative error

圖4 四周固支板彎曲撓度Fig.4 The bending of plate with all clamped boundaries

圖5 布點方式Fig.5 Pattern of nodes

對于一階剪切變形理論的中厚板問題,當板變得比較薄時,因剪切應變能過大產生錯誤的結果,這種現象稱為剪切自鎖[29].在無網格GFDM 中,通過提高泰勒展開的階次來構造節點的支撐域,能夠簡單有效地避免剪切自鎖現象的發生.如圖6 所示,采用6 階展開項能夠有效避免剪切自鎖的發生,相較與在有限元中采用縮減積分法、假設剪切應變法等[28],其簡單高效地提高了求解方法的通用性.

圖6 高階次展開項克服剪切自鎖現象Fig.6 Higher-order expansion to overcome the phenomenon of shear locking

3.2 FG-CNTRC 板結構及等效模型參數

如圖1 所示的FG-CNTRC 板結構模型,根據Shen[4]提出的廣義混合律模型,用于計算的FG-CNTRC等效模型參數如表2 和表3 所示.

表2 材料屬性Table 2 Material properties

表3 碳納米管的效能參數Table 3 Efficiency parameters of CNTs

3.3 彎曲分析

首先在沒有考慮碳納米管轉向的情況下,表4計算了體積分數為0.11,不同分布型下的四周固支FG-CNTRC 板的彎曲變形,并與已有文獻Zhu 等[8]的數值結果進行了參考比較,兩者都能很好地近似,說明基于一階剪切變形板理論的廣義有限差分法求解FG-CNTRC 板彎曲問題的計算正確和有效性.圖7為體積分數為0.11,寬厚比 a/l=10,四周固支FGCNTRC 板中心彎曲隨離散點增加的收斂情況,在總點數達到300 左右時結果近乎收斂,能夠看出廣義有限差分法具有計算快、收斂穩定的特性.

圖7 收斂性分析Fig.7 Analysis of convergence

表4 四周固支FG-CNTRC 板的中心彎曲w/lTable 4 The center bending w/l of FG-CNTRC plate with all clamped boundaries

圖8 展示了含碳納米管轉向的FG-CNTRC 板隨傾斜角 α 變化的各向異性彎曲狀態(β=-α),分別從X 軸向和Y 軸向描述了該板的彎曲變形.碳納米管具有指向性使得在X 和Y 方向上存在著各向異性的剛度變化.圖9 展示了厚跨比為0.1 的四周固支FG-CNTRC 板在載荷 q=0.1 MPa 下中心點處的 σ 應力,可以看出UD 型、FG-O 型和FG-X 分布型下的CNTRC 板應力與 z=0 成中心對稱,即 z=0 為該FGCNTRC 板的中性層,原因是UD 型、FG-O 型和FG-X 型碳納米管分布呈對稱性分布.對于沒有橫向對稱性的FG-V 碳納米管分布型,板結構中性層大概在 z=0.15 .且 σxx與 σyy應力大小不等,沿厚度方向變化情況也不相同,正因如此使得FG-CNTRC 板在宏觀上具有如圖8 所示的各向異性彎曲狀態.

圖8 FG-CNTRC 板的各向異性彎曲Fig.8 Anisotropic bending of FG-CNTRC plate

圖9 四周固支FG-CNTRC 板的中心應力Fig.9 Central stress of FG-CNTRC plate with all clamped boundaries

3.4 自振頻率分析

下面運用廣義有限差分法對含碳納米管轉向的FG-CNTRC 板的自振模態問題進行詳細分析與討論.表5展示了不同分布型下四周固支FGCNTRC 板的前6階自振頻率且與Zhu 等[8]的數值結果進行了參考比較,說明廣義有限差分法這一算法的計算正確和有效性.然后運用廣義有限差分法數值分析了不同轉向的碳納米管對FG-CNTRC 板結構自振頻率的影響.圖10 展示了不同分布型和不同傾斜角度 α 下,FG-CNTRC 板中碳納米管的旋轉角度β對結構第一階自振頻率的影響,可以看出碳納米管在橫向布置 β=0°和縱向布置 β=90°時其結構具有相同的自振效果,在 β=45°斜向配置時其結構抗彎剛度最小,自振頻率最小;隨著傾斜角 α 的增大,隨同碳納米管旋轉的FG-CNTRC 板數值結果不再呈周期性變化;在 α=30°時板最小自振發生在碳納米管β=30°左右;隨著傾斜角α 的繼續增大,使FGCNTRC 板產生最小自振頻率的碳納米管轉向角度從 β=45°逐漸減小.

表5 四周固支FG-CNTRC 板的前6 階自振頻率Table 5 Six natural frequencies of FG-CNTRC plate with all clamped boundaries

圖10 碳納米管的轉向對自振頻率的影響Fig.10 The effect of CNT orientation on natural vibration

圖11 展示了UD 分布型和不同轉向角度 β 下厚跨比l/a對FG-CN TRC板第一階自振頻率的影響,可以看出FG-CNTRC板厚跨比與自振頻率成正相關,厚度增加的剪切效應使碳納米管轉向對自振頻率的影響逐漸降低.圖12 討論了結構不同長寬比 a/b 所構成的FG-CNTRC 板對自振頻率的影響,發現當 a/b ≥1 碳納米管在β=-90°,90°時,即豎向放置,自振頻率大小近乎相等;在橫向放置 β=0°時,其隨不同的長寬比 a/b 配置產生不一樣的自振效果,且碳納米管的放置角度對自振頻率產生顯著影響.

圖11 CNTRC 板隨厚跨比變化的自振頻率Fig.11 The effect of CNT orientation on natural vibration of CNTRC plate ranging with thickness ratio

圖12 FG-X 型CNTRC 板隨長寬比變化的自振頻率Fig.12 The effect of CNT orientation on natural vibration of FG-XCNTRC plate ranging with length-width ratio

4 結論

本文采用廣義有限差分法求解了基于一階剪切變形理論的功能梯度碳納米管增強復合材料(FGCNTRC)板的靜態線性彎曲和自振模態問題.通過基準算例驗證了該算法的精確性和高效性,是一種基于移動最小二乘發展而來的簡單高效的區域型無網格方法,且能夠快速方便地克服一階剪切變形理論在薄板應用中的剪切自鎖現象.接著,在此基礎上進一步研究具有轉向的碳納米管、板寬厚比、傾斜角度等材料結構參數對FG-CNTRC 板彎曲和自振模態的影響.研究發現:

(1) 碳納米管排布形式對結構剛度的影響差異性大,FG-CNTRC 材料所形成的結構剛度大小依次有FG-O＜FG-V＜UD＜FG-X 型;

(2) 隨著板傾斜角的增大,同碳納米管旋轉的FG-CNTRC 板自振頻率不再呈周期性變化;

(3) 碳納米管 β=-90°,90°時,即豎向放置,結構長寬比的改變近乎不影響結構剛度;在橫向放置β=0°時,隨不同的長寬比變化產生不一樣的結構剛度;

(4) 在X 和Y 方向上存在著各向異性的剛度變化,且碳納米管的轉向角度對結構產生顯著影響.最后需要指出的是,本文主要側重于驗證廣義有限差分法在計算FG-CNTRC 板結構力學特性的有效性,后續將基于廣義有限差分法考慮更為復雜的數值計算模型.