非均質復合材料力學性能的確定性多尺度計算方法

裴世源,徐華,2

(1.西安交通大學現代設計與轉子軸承系統教育部重點實驗室, 710049, 西安;2.新疆大學機械工程學院, 830046, 烏魯木齊)

?

非均質復合材料力學性能的確定性多尺度計算方法

裴世源1,徐華1,2

(1.西安交通大學現代設計與轉子軸承系統教育部重點實驗室, 710049, 西安;2.新疆大學機械工程學院, 830046, 烏魯木齊)

針對具有多尺度特性的非均質復合材料進行結構強度分析,提出了一種確定性多尺度計算方法——有限細胞法(FCM)。該方法通過矩陣凝聚和插值近似得到粗網格的剛度矩陣,在計算得到粗網格解后,通過回代技術可獲得全域細網格解。為驗證FCM的計算精度和計算效率,構造了多個數值算例,并與理論和有限元(FEM)結果進行了對比,結果發現:相對于FEM,FCM可顯著擴大計算規模,提高計算速度,且易于并行;相對于其他多尺度計算方法,FCM在構造粗網格剛度矩陣過程中無需引入微觀邊界條件,計算精度高,易于降尺度計算。因此,FCM是分析復合材料結構強度的一種有效的多尺度算法,同時該方法易于推廣至多尺度的非線性結構分析、熱分析、地下水滲流分析等領域。

有限細胞法;多尺度計算;非均質復合材料

自然界和工程應用中的很多材料都具有多尺度特征,準確描述這些材料的整體細觀力學行為是當今科學界和工程界所面臨的挑戰之一。傳統數值方法在求解這些問題時需耗費大量的計算資源,雖然超級計算機和并行計算的出現在一定程度上緩解了計算難度,但對于很多大規模多尺度問題仍難以求解。在保證計算精度的前提下,發展高效的多尺度計算方法是當今計算力學領域的研究熱點[1-2]。

近些年來,國內外學者已提出了各種多尺度計算方法,其中具有代表性的有漸近均勻化法[3-7]、多尺度有限元法(MsFEM)[8-10]、非均質多尺度法(HMM)[11-12]、代表體元法(RVE)[13-15]、小波均勻化方法[16]等。這些方法通常是從宏觀域中取出有代表性的體元(微觀域),在此體元上構建一個局部邊值問題,在得到體元的等效力學參數后,將原結構等效成具有此參數的均勻結構求解宏觀域的力學響應。一般來說,選擇局部問題邊界條件具有較強的主觀性和經驗性,對計算精度具有顯著影響[9-10],因此構造局部問題合適的邊界條件是多尺度計算的重點和難點之一。另一方面,雖然這些多尺度方法均可用于分析結構的宏觀力學行為,但通常難以得到結構的微觀力學響應(如局部應力和應變等)[2,11],在很多情況下獲得全域精細網格解是至關重要的,如結構強度,疲勞和斷裂分析等.

作者曾基于有限元子結構方法提出了有限細胞法(FCM),用于求解表面結構的多尺度潤滑問題(標量場問題)效果良好[17-18]。該方法的基本思想與傳統多尺度方法類似,亦是先構建微觀域的局部問題,后求解宏觀域整體響應。主要不同點在于:①構建局部問題時,無需引入微觀邊界條件,避免了人為選擇的主觀性和經驗性;②直接通過矩陣運算得到粗網格(細胞)的剛度矩陣和載荷向量,而不是得到微觀等效材料參數后,再構建宏觀單元的剛度矩陣;③易于降尺度計算。

為對非均質多尺度復合材料進行結構分析,本文以二維連續體問題分析為例,介紹FCM求解矢量場問題的基本原理及其實施過程,并通過幾個有代表性的數值算例說明其有效性。限于篇幅,本文只介紹FCM的基本原理及其在復合材料力學分析中的應用,多尺度潤滑問題的分析可參見文獻[17-18],有關FCM在多孔介質滲流、傳熱等問題中的應用將另文論述。

1 有限細胞法

假定有結構如圖1所示,若將最小的重復單位定為細胞,則該結構由4個細胞構成,整體域可由方形細胞和圓形細胞兩種組合而成。通常FCM的實施包括以下4個步驟,如圖2所示。

圖1 宏觀域結構示意圖

(a)凝聚細胞內部自由度 (b)凝聚細胞邊界自由度

(c)求解宏觀域響應 (d)降尺度計算 圖2 有限細胞法關鍵步驟

1.1 凝聚細胞內部自由度

為便于說明,使用雙線性四邊形單元對細胞進行結構化網格剖分,每個方向均具有M個節點,如圖2a所示。使用常規有限元法組裝得到細胞整體剛度矩陣和載荷向量,細胞的平衡方程可以表示為

Kcuc=fc

(1)

式中:Kc表示細胞的剛度矩陣;uc表示細胞所有節點的自由度;fc表示載荷向量;下標c表示細胞。需要說明的是,在建立局部問題(細胞平衡方程)時,無需施加邊界條件。對于彈性力學問題fc是具有2M2個0元素的列向量,但對于傳熱、滲流或潤滑等問題,fc不一定是零向量。經過適當的節點編號,式(1)可用矩陣分塊表示為

(2)

式中:下標b表示細胞邊界;i表示細胞內部;ub表示細胞邊界上的自由度;ui表示細胞內部的自由度。基于消去法對ui凝聚,式(2)可表示為

(3)

(4)

(5)

1.2 凝聚細胞邊界自由度

圖3 細胞節點分類

為進一步減小細胞的等效剛度矩陣,在細胞的4條邊界上分別選取保留均布的k個主節點,邊界上的其余節點定義為從節點,如圖2b所示,節點分類如圖3所示。若假定主節點的自由度已知,則從節點的自由度可通過k個主節點的自由度插值計算得到。為方便實施,使用拉格朗日p次多項式作為插值函數,其中p≤k-1。例如,從節點中的第l個自由度可由同側邊界上的p+1個主節點的對應自由度插值得到

(6)

式中:下標s表示從節點;m表示主節點;φj表示拉格朗日插值基函數。將式(6)改寫為向量形式

(7)

通過矩陣疊加,全體從節點自由度可表示為

us=Tum

(8)

(9)

式中:T表示插值矩陣。

根據主節點與從節點的分類,式(3)可以用矩陣分塊表示為

(10)

將式(8)代入到式(10),可得

(11)

式中:I表示單位矩陣。式(11)第一行可表示為

Kmum=fm

(12)

Km=kmm+kmsT

(13)

即細胞的等效剛度矩陣和載荷向量分別可用Km和fm表示。在此凝聚過程中使用主節點的自由度插值計算從節點的自由度,忽略了部分信息,計算精度會有所損失,但通常在細胞邊界上的節點的自由度不會產生間斷或突變,當在細胞邊界上保留較多的主節點和選取高階插值函數后,這時精度損失變得可以忽略。一般來說,線性和二次插值函數通常誤差較大,當插值函數的階數高于三次后,計算誤差會顯著降低。本文將在數值算例部分詳細討論FCM插值方案對計算精度的影響。另一方面,經過插值近似,細胞待求的自由度取決于主節點的數量,即8(k-1),對于固體力學問題,當細胞邊界上保留11個主節點時,FCM與FEM精細網格間計算結果之間的誤差已經低于1%,此時細胞的等效剛度矩陣Km為80×80稠密矩陣。通過兩個凝聚步驟得到規模較小細胞等效剛度矩陣和載荷向量是有限細胞法的關鍵步驟。對于周期性結構,細胞的剛度矩陣僅需計算一次,對于非周期結構細胞剛度矩陣需要分別計算,但不同細胞的剛度矩陣的計算完全獨立,因此易于并行。

1.3 求解宏觀域響應

如圖2c所示,宏觀域可以認為由有限個細胞組成。在得到的等效細胞剛度矩陣Km和載荷向量fm后,系統的總體剛度矩陣和載荷向量可表示為

Kum=f

(14)

(15)

式(15)表示系統總體剛度矩陣由每個細胞的等效剛度矩陣組裝而成。施加宏觀邊界條件后,式(14)可解,從而得到宏觀域細胞主節點的自由度。

1.4 降尺度計算

對比于傳統多尺度方法(均勻化,MsFEM和HMM等),FCM能很容易地進行降尺度計算,且可以保證在細胞邊界處自由度連續。如圖2d所示,在得到主節點的自由度后,從節點的自由度可由式(8)計算得到,由此便得到了細胞的邊界自由度ub,而細胞內部的自由度為

(16)

至此整體域精細網格的全部自由度均已求出,應力和應變等信息可用傳統有限元方法計算得到。

從以上推導過程可以看出:有限細胞法完全基于矩陣變換,無需引入微觀邊界,易于降尺度計算;同時該方法與控制方程的無關,因此易于將其推廣至非線性結構分析、熱分析、地下水滲流等領域。

2 數值算例

本節將以幾個典型的算例來驗證有限細胞法的計算精度和計算效率等問題,所有計算均在個人臺式機上完成(CPU 3.2 GHz,RAM 16G),單核計算,編程語言為MATLAB 8.3,操作系統為Windows 7(64 bit)。多尺度精細網格解用FCMk×p表示,k表示細胞單側邊界上的節點個數,p表示拉格朗日插值函數的階次,把傳統有限元法精細網格解作為參考解,用FEM表示。算例中的所有參數都假設是無量綱量,且本文考察的是平面應力情況。

2.1 算例1

算例1為均質材料懸臂梁結構分析。為了驗證FCM的正確性,構造一個簡單的具有矩形截面的均質懸臂梁算例,如圖4所示,懸臂梁結構的右端固定支撐,左端承受P=100的分布力,懸臂梁長度L=100,高度h=10,寬度為1,材料的彈性模量為2×105,泊松比為0.3。懸臂梁中性軸撓度的理論表達式為[19]

(17)

圖4 均質懸臂梁結構示意圖

本算例中FCM使用了6×5的插值方案,即在細胞的每條邊界上保留了6個主節點,從節點的自由度使用5次拉格朗日函數由主節點自由度插值得到。FCM使用40×40的均布雙線性四邊形單元對細胞進行劃分,采用20×2的細胞對全域進行劃分,故全域細網格尺度共有8 000×800單元,有限元采用相同的細網格密度。通過L2范數[12]驗證數值方法的實際效果

(18)

式中:Q為懸臂梁中性軸相應的細網格節點數。本算例中vref取解析解va,即式(17)。懸臂梁中性軸y向位移對比如圖5所示。由FCM和FEM的L2范數分別為3.15×10-3和3.35×10-3可以發現,對于均質材料結構FCM在相同的細網格尺度上可以得到比傳統FEM更高的精度,且收斂于精確解。同時,FCM的計算時間僅為FEM的26.2%,計算效率也更高。

圖5 均質懸臂梁中性軸y向位移比較

2.2 算例2

算例2為具有周期結構的非均質懸臂梁結構分析。懸臂梁模型及其邊界條件如圖6所示。該結構由A×B個材料相同的正方形細胞組成,本算例中考慮了5種A×B的組合,即6×2、12×3、18×7、24×13、30×29和40×30,細胞中的邊長為1,包含邊長為0.3的正方形強化相。其中,基體材料的彈性模量為2×104,泊松比為0.3。強化相的彈性模量為2×105,泊松比為0.25,細胞采用40×40的均勻雙線性四邊形網格劃分,如圖7所示。懸臂梁左端固定,上端承受強度q=1的均布載荷。

圖6 復合材料結構示意圖

圖7 細胞網格劃分

本算例中考慮了FCM的兩種插值方案,分別為6×5和11×9,不同計算規模下,FEM和FCM計算結果的比較如表1所示。本算例中取相同網格密度下FEM的計算結果作為參考解,當A=12、B=3時,懸臂梁中性軸y向位移對比如圖8所示。

從表1可以看出,FCM和FEM的計算結果吻合良好,FCM 6×5方案的L2范數最大誤差為22.8×10-3,而FCM 11×9方案的L2范數最大誤差僅為2.43×10-3。總體看來,FCM 11×9方案的計算結果優于FCM 6×5方案。這是因為在邊界上保留更多的主節點和采用高次的插值函數,從而可以更精確地計算從節點的自由度。對于結構分析來說,細胞單側邊界上保留11個主節點,FCM已經具有足夠高的精度,保留更多的主節點并不能繼續有效提高精度,但會增加計算時間和內存。因此,綜合考慮精度和效率,在后續的算例中FCM將采用11×9的插值方案。

圖8 懸臂梁中性軸y向的位移比較(A=12,B=3)

A×B總自由度/103L2范數/10-3FCM6×5FCM11×9計算時間/sFCM6×5FCM11×9FEM6×238.46.262.4301.982.673.2312×311522.800.1732.063.499.7218×740311.000.3083.393.9235.1024×139986.910.3364.516.36118.8030×2927844.300.4096.9615.10851.0040×3038405.540.1857.7521.101643.50

在計算時間方面,隨著總自由度的增加,傳統有限元方法的計算時間顯著增加,當總自由度達到380萬時,FEM的計算時間接近0.5 h,且內存消耗也接近于物理內存極限(16 GB),因此我們沒有嘗試更大規模的計算,但在相同的網格密度下FCM的計算時間均不足22 s,內存需求也非常低(小于1 GB)。這是因為對于周期結構,細胞的剛度矩陣僅需構造一次,且規模較小,故對于周期性結構,有限細胞法可顯著降低計算量,節省內存,擴大了計算規模,具有明顯的優勢。

2.3 算例3

算例3為非完全周期結構的非均質懸臂梁結構分析。懸臂梁由30×6個細胞組成,如圖9所示。細胞的幾何尺寸、網格劃分、材料性質和結構的邊界條件與算例1相同,但在豎直方向上,復合材料增強相的邊長由0.1變化到0.6,步長為0.1。

圖9 非均質懸臂梁模型

(a)FEM

(b)FCM圖10 宏觀von Mises應力云圖比較

(a)FEM (b)FCM 圖11 左下角細胞微觀尺度von Mises應力云圖比較

圖12 懸臂梁中性軸y方向的位移比較

圖13 過應力最大點軸線方向von Mises應力比較

von Mises整體應力對比如圖10所示,細胞微觀尺度von Mises應力云圖如圖11所示,懸臂梁中性軸y向位移比較如圖12所示,過應力最大點軸線方向von Mises應力比較如圖13所示,發現FCM的計算結果與FEM參考解吻合良好。注意到懸臂梁底部增強相尺寸已經接近于細胞尺寸,在這種情況下一些傳統多尺度計算方法(均勻化,MsFEM等)的計算精度會急劇下降,但從圖10至圖13可以看到,FCM與FEM的結果依然高度吻合,說明FCM具有良好的適應性和靈活性。在計算時間方面,本算例中總自由度為57.6萬,FEM的計算時間為50.2 s,而FCM的計算時間僅為13.3 s,說明FCM在非完全周期情況下依然具有很高的計算精度和較高的計算效率。

若宏觀結構沒有任何的周期性,傳統有限元方法可能因為計算時間和內存需求過大而無法計算,但有限細胞法可以將原問題分解為很多小問題后分別求解,且不同細胞剛度矩陣的構造完全獨立,易于將細胞剛度矩陣的構造進行并行化處理,且不同線程之間的通信代價極低,因此FCM為分析非均質非周期多尺度材料的力學性能提供了一種方法。

3 結 論

本文介紹了有限細胞法(FCM)的基本原理,詳細闡述了該方法的理論基礎和關鍵實施步驟,展示了其在復合材料結構分析中的應用。FCM的基本思想是通過矩陣運算直接得到多節點的粗網格(細胞)剛度矩陣,在求解得到粗網格尺度的力學響應后,通過降尺度計算到全域精細網格尺度的響應。通過數值算例驗證有限細胞法的計算精度和計算效率等問題,結果表明,相對于傳統FEM,FCM可顯著提高計算效率,節省內存,擴大計算規模,尤其是對于周期材料,優勢非常明顯。相對于經典多尺度方法,FCM在構建局部問題時無需引入微觀邊界條件,易于進行降尺度計算,且具有與控制方程無關、方便實施等優點。FCM與控制方程無關的特性決定了其易于推廣至多尺度的非線性結構分析、熱分析、地下水滲流等領域。

[1] ROUET-LEDUC B, BARROS K, CIEREN E, et al. Spatial adaptive sampling in multiscale simulation [J]. Computer Physics Communications, 2014, 185(7): 1857-1864.

[2] KANOUTé P, BOSO D P, CHABOCHE J L, et al. Multiscale methods for composites: a review [J]. Archives of Computational Methods in Engineering, 2009, 16(1): 31-75.

[3] FISH J, SHEK K, PANDHEERADI M, et al. Computational plasticity for composite structures based on mathematical homogenization: theory and practice [J]. Computer Methods in Applied Mechanics and Engineering, 1997, 148(1): 53-73.

[4] 易壘, 文毅. 慣性載荷作用下火箭橇結構拓撲優化設計 [J]. 應用力學學報, 2013, 30(1): 80-85 YI Lei, WEN Yi. Topology optimization of rocket sled structure under inertial loads [J]. Chinese Journal of Applied Mechanics, 2013 , 30(1): 80-85.

[5] MIEHE C C B. On multiscale FE analyses of heterogeneous structures: from homogenization to multigrid solvers [J]. International Journal for Numerical Methods in Engineering, 2007, 71(10): 1135-1180.

[6] ZHANG H W, ZHANG S, BI J Y, et al. Thermo-mechanical analysis of periodic multiphase materials by a multiscale asymptotic homogenization approach [J]. International Journal for Numerical Methods in Engineering, 2007, 69(1): 87-113.

[7] 張洪武, 張盛, 畢金英. 周期性結構熱動力時間空間多尺度分析 [J]. 力學學報, 2006, 38(2): 226-235. ZHANG Hongwu, ZHANG Sheng, BI Jinying. Thermodynamic analysis of multiphase periodic structures based on a spatial and temporal multiple scale method [J]. Chinese Journal of Theoretical and Applied Mechanics, 2006, 38(2): 226-235.

[8] HOU T Y, WU X H. A multiscale finite element method for elliptic problems in composite materials and porous media V. Journal of Computational Physics, 1997, 134(1): 169-189.

[9] 張洪武, 吳敬凱, 劉輝, 等. 擴展的多尺度有限元法基本原理 [J]. 計算機輔助工程, 2010, 19(2): 3-9 ZHANG Hongwu, WU Jingkai, LIU Hui, et al. Basic theory of extended multiscale finite element method[J]. Computer Aided Engineering, 2010, 19(2): 3-9.

[10]ZHANG H, WU J, Lü J, et al. Extended multiscale finite element method for mechanical analysis of heterogeneous materials [J]. Acta Mechanica Sinica, 2010, 26(6): 899-920.

[11]E Weinan, ENGQUIST B, LI X, et al. Heterogeneous multiscale methods: a review [J]. Communications in computational physics, 2007, 2(3): 367-450.

[12]E Weinan, ENGQUIST B, HUANG Z. Heterogeneous multiscale method: a general methodology for multiscale modeling [J]. Physical Review: B, 2003, 67(9): 092101.

[13]XIA Z, ZHOU C, YONG Q, et al. On selection of repeated unit cell model and application of unified periodic boundary conditions in micro-mechanical analysis of composites [J]. International Journal of Solids and Structures, 2006, 43(2): 266-278.

[14]LI X, LIU Q, ZHANG J. A micro-macro homogenization approach for discrete particle assembly Cosserat continuum modeling of granular materials [J]. International Journal of Solids and Structures, 2010, 47(2): 291-303.

[15]PECULLAN S, GIBIANSKY L V, TORQUATO S. Scale effects on the elastic behavior of periodic and hierarchical two-dimensional composites [J]. Journal of the Mechanics and Physics of Solids, 1999, 47(7): 1509-1542.

[16]DOROBANTU M, ENGQUIST B. Wavelet-based numerical homogenization [J]. SIAM Journal on Numerical Analysis, 1998, 35(2): 540-559.

[17]裴世源, 徐華, 馬石磊, 等. 多尺度表面織構流體潤滑問題的快速求解方法 [J]. 西安交通大學學報, 2011, 45(5): 119-126. PEI Shiyan, XU Hua, MA Shilei, et al. Multiscale method for modeling surface texture effects in hydrodynamic lubrication regime [J]. Journal of Xi’an Jiaotong University, 2011, 45(5): 119-126.

[18]PEI S, MA S, XU H, et al. A multiscale method of modeling surface texture in hydrodynamic regime [J]. Tribology International, 2011, 44(12): 1810-1818.

[19]TIMOSHENKO S, GOODIER J N. Theory of elasticity [M]. New York, USA: McGraw-Hill, 1951: 506.

(編輯 杜秀杰)

Deterministic Multiscale Method for Heterogeneous Composite Material

PEI Shiyuan1,XU Hua1,2

(1. Key Laboratory of Education Ministry for Modern Design and Rotor-Bearing System, Xi’an Jiaotong University,Xi’an 710049, China; 2. School of Mechanical Engineering, Xinjiang University, Urumqi 830046, China)

For the strength analysis of heterogeneous materials in elasticity, a deterministic multiscale calculation method, finite cell method (FCM), is presented. The condensation and interpolation technique are adopted to obtain the coarse element stiffness matrix, then the original problem can be solved in coarse-scale, and the fine-scale solution can be sought out by back substitution. The comparison with analytic method and finite element method (FEM) on several numerical examples indicates that FCM enables to significantly expand computational scale and greatly improve computing rate. In FCM, the artificial boundary conditions are unnecessary for constructing macroscopic meshes, and the downscaled computing can be easily performed. FCM is also feasible for non-linear multiscale structure analysis, thermal analysis and porous flow analysis.

finite cell method; multiscale computation; heterogeneous material

2015-01-13。

裴世源(1983—),男,講師。

國家重點基礎研究發展計劃資助項目(2011 CB706602);陜西省科技攻關資助項目(2015GY022)。

時間:2015-06-29

10.7652/xjtuxb201510002

TP301.6

A

0253-987X(2015)10-0008-06

網絡出版地址:http://www.cnki.net/kcms/detail/61.1069.T.20150629.1137.003.html