一種三維線彈性隨機數值均化方法

王 南,馬 娟,張玉林

(西安電子科技大學機電工程學院,陜西西安 710071)

一種三維線彈性隨機數值均化方法

王 南,馬 娟,張玉林

(西安電子科技大學機電工程學院,陜西西安 710071)

對小變形下具有隨機微觀結構非均質材料的數值均化問題進行了研究.當考慮材料微觀結構形態、各組分性質的隨機性及相關性時,基于多尺度有限元方法和蒙特卡羅法,建立了非均質材料的隨機均化模型.求出了不同邊界條件下的隨機有效量及其數字特征值以及有效量之間所具有的相關性.考察了隨機微觀結構參數對隨機有效量以及表征體積單元中應力分布的影響.結果表明,在分析非均質材料的宏觀有效力學性質時,考慮材料微觀結構中客觀存在的隨機性和相關性非常重要.

均化;隨機性和相關性;三維線彈性;有限元法;蒙特卡羅法

非均質材料在持續的機械應力和熱應力的作用下,其微觀特性對材料整體性能[1]有著至關重要的影響.就非均質材料而言,均化方法作為一種將材料微觀特性和宏觀響應聯系起來的有效方法已經發展了數十年.用均化方法得到的等效性質可以把宏觀的非均質材料轉化為等效的均值材料,而且極大地減少了求解時需要的自由度的個數.此外,通過引入應用了功準則的有效特征,均化方法亦極大地減少了解決宏觀問題所需的計算量.就線彈性非均質材料的均化而言,目前已發展的較為完善,其中解析均化法可參見文獻[2],數值均化方法參見文獻[3].就其實質而言,均化是通過確定和分析具有統計意義的非均質材料的表征體積單元(Representative Volume Element,RVE)[4],獲得非均質材料彈性張量的估值或一組邊界值.

非均質材料的重要特征之一是其成分材料、缺陷分布以及這些因素間的相互作用的不確定性.近年來,對材料參數和荷載不確定性[5]的認識使得隨機非均質材料得到關注并取得了一些研究成果[6-8].文獻[6]解決了多孔領域中障礙問題的均勻化,其中孔洞是周期分布且擁有隨機的形狀和尺寸.文獻[7]基于高斯隨機場的非線性轉化重構了帶有隨機形態的兩相復合材料,并按照給定的統計特征值提出了一種高效的計算方法.文獻[8]采用廣義變分原理將一個帶有隨機微觀結構的邊值問題分解為一個慢尺度的確定性問題和一個快尺度的隨機問題.

盡管取得了一些進展,但對于非均質材料的隨機均化分析仍有待進一步深入研究.由于該類問題需從非均質材料不同尺度上的不確定變量入手來求解,已有的隨機均化方法主要是通過減少微觀結構中隨機變量的數目來達到簡化的目的的,其大致可分為兩類,一類是僅考慮微觀結構形態的不確定性[5-10];另一類是直接采用組分的楊氏模量和泊松比,而不是采用多尺度有限元法從組分的體積模量和剪切模量進行推導得出[11].此外,已有的均化模型幾乎均未涉及微觀結構參數間的相關性和均化結果之間的相關性.而微觀結構所具有的隨機性和相關性是客觀存在、不容回避的,故亟需發展充分考慮微觀結構不確定性的隨機均化分析模型,并建立非均質材料微觀結構的不確定性與材料整體宏觀性質不確定性之間的量化關系.

筆者對小變形下三維線彈性隨機非均質材料的均化問題進行了研究.當同時考慮微觀結構所具有的隨機性及相關性時,基于多尺度有限元方法以及蒙特卡羅法構建了隨機均化模型,并對不同邊界條件下材料的隨機有效量如有效模量、有效彈性參數和有效彈性張量及其數字特征值進行了求解;最后考察了微觀結構的不確定性對隨機均化結果的影響,得到了一些有意義的結論.

1　線彈性非均質材料的均化問題

1.1均化問題

設由均勻各向同性的基體(M1)和顆粒(M2)材料組成非均質材料M如圖1所示,其參考形狀為R0,本構方程σ=(X,ε),其中σ和ε分別表示應力和應變,X為表征R0位置的矢量.圖1中,L表示介觀,d表示微觀,D表示宏觀,且將宏觀非均質材料用等效均質材料替代.

圖1　二相非均質材料的均化問題及其微觀-介觀-宏觀的尺度準則

關于該非均質材料M的一個力學邊值問題為:求解位移u(X,t),以便下述方程在R0上成立:

其中,在?Ru、?Rt上的邊界條件分別是:u=和t=σn=,ρ是材料的密度,ρb和ρ分別為宏觀體力和慣性力.因M固有的非均質性所引起的求解困難,采用等效均質材料M*替代M后得到式(1)中x的近似解x*,故關于M*的均化問題成為:確定u*(X,t),以使式(2)在R0上成立:

其中,?Ru、?Rt上的邊界條件分別是:u*=和t*=σ*n=,本構方程為σ*=*(X,ε*).

1.2有效本構關系

均化結果的準確性與材料的表征體積單元直接相關.考慮圖1所示非均質材料中的一個表征體積單元,其體積為V0,即V0?R0.在微觀結構分析中,忽略宏觀體力和加速度,且邊界條件服從Hill功準則.若將表征體積單元進行有限元網格劃分,則均化過程中需要求解的邊值問題即簡化為求解u(X),以便

材料的有效彈性張量C*可由和計算得到[3]

1.3幾類特殊的邊界條件

在表征體積單元上施加的邊界條件可分為由平均應力Σ控制的邊界條件和由平均應變ε控制的邊界條件,其中平均應力Σ為給定的常量,而平均應變ε=(〈H〉+〈H〉T)/2中平均位移梯度〈H〉是給定的常量.均化理論中4種廣泛使用的由平均應變ε和平均應力Σ控制的邊界條件[12]分別為:平均應變控制的線性位移邊界條件(ε-Linear Displacement Boundary Conditions,ε-LD-BCs)、平均應變控制的周期性邊界條件(ε-PeRiodic Boundary Conditions,ε-PR-BCs)、平均應變控制的均勻應力邊界條件(ε-Uniform Traction Boundary Conditions,ε-UT-BCs)和平均應力控制的均勻應力邊界條件(Σ-UT-BCs).

欲求解C*,可通過在表征體積單元上施加邊界條件并通過體積平均獲得其有效應力后按式(6)獲得.

1.4基于有限元法的宏觀有效量

此外,將邊界條件施加于表征體積單元得到有效應力σ*和有效應變ε*后,有效體積模量和有效剪切模量分別為

其中,有效應力σ*和有效應變ε*是表征體積單元上的體積平均量和;張量A的(′)部分定義為A′=A-(tr(A)/3)I.

三維線彈性材料的有效楊氏模量E*和泊松比ν*可計算如下:

2　基于蒙特卡羅法的線彈性非均質材料的隨機均化

考慮由兩種均勻各向同性材料組成的二相非均質材料,即在基體中加入不同材質、隨機分布的圓形顆粒.基體和顆粒的體積模量和剪切模量分別為k1、u1和k2、u2,顆粒體積比V2,半徑為r.隨機均化流程如圖2所示.假設k1、 u1、k2、u2、V2和r為正態分布的隨機變量:、u1～,其中μ和dmsd分別為變量的均值和均方差.k1、 u1和k2、u2的相關系數分別為和.利用蒙特卡羅法可生成各隨機參數變量的m個樣本,其中第j組樣本為k1j、k2j、u1j、u2j、V2j和rj,且k1j、u1j和k2j、u2j是根據隨機向量協方差矩陣的喬列斯基分解[13]并按照和產生的,故滿足及

圖2　隨機均化流程圖

用隨機序列添加算法[14]生成表征體積單元并將其劃分為n個單元后,對其施加邊界條件.對第i個單元,將第j組參數樣本代入一個內部有限元均化程序中,首先可獲得第i個單元上第l個高斯點的隨機微觀應力和微觀應變.之后按照每個高斯點的權重系數ωl可得到第i個單元上的隨機應力,即

進一步由體積平均可得到表征體積單元上的隨機有效應力為

同理可由式(6)～(8)得到表征體積單元上的有效彈性張量、有效體積模量、有效剪切模量、有效楊氏模量和有效泊松比,將這些有效量統一用Q*表示.以此類推,m組參數樣本代入程序后即獲得m組有效量Q*.由數理統計即可求得有效量Q*的均值μ、均方差dmsd和變異系數CV為

3　數值算例

針對一種顆粒隨機分布的二相非均質材料,表1給出了基體及顆粒的各向同性材料參數的均值和變異系數,k1、u1和k2、u2之間的相關系數

表1　一種二相隨機非均質材料的微觀結構參數

下面均化計算中所用的微觀結構參數均為表1中的值,蒙特卡羅法模擬的次數為5000次.

表2是3種邊界條件下求解獲得的隨機有效彈性參數E*和ν*.從表2可見,ν*的均方差隨著V2的增大而增加,而E*的數字特征值和ν*的均值則相反;不同邊界條件下ν*的均方差滿足不等式(·)ε-LD-BCs＞(·)ε-PR-BCs＞(·)ε-UT-BCs;而E*的數字特征值和ν*的均值則滿足不等式(·)ε-LD-BCs＜(·)ε-PR-BCs＜(·)ε-UT-BCs.圖4則給出了±3dmsd原則下有效模量的取值范圍,可見ρ=0.6時的取值范圍大于ρ=0時的取值范圍.從表2和圖4可知,是否考慮參數間的相關性對均化結果是有影響的.

圖3　非均質材料線彈性均化的解析解和數值解的均值比較

表2　隨機有效彈性參數的均值和均方差

圖4　ε-LD-BCs下兩種相關情況時k*和u*的取值范圍

有效彈性張量C*中部分元素的數字特征值列于表3中,明顯可見其數字特征值隨著V2的增加而增加.

表3　有效彈性張量C*部分元素的數字特征值

當各隨機參數的變異系數均為0.05時,考察了有效量間的相關性與隨機參數間相關性的關系.由圖5可知:隨著隨機參數間ρ的增加,k*和u*間的正相關性顯著增強,而E*和ν*間的負相關性也隨之增強.

圖5　當V2=0.25時和隨變化的情況

圖6　ρ=0.9時3種邊界下表征體積單元中隨機應力的分布圖

從圖6可見,在ε-PR-BCs下,表征體積單元的邊界相對于ε-LD-BCs下的邊界呈周期性變化,而ε-UT-BCs下表征體積單元邊界的波動則無規律,這與各邊界條件的本質密切相關;此外,不同的邊界條件下應力分布的上界、均值和下界不相同;各圖中應力分布均不均勻,離顆粒越近的基體中單元應力梯度越大且應力也越大,反之則越小.

圖7　ε-LD-BCs下,當ρ=0時表征體積單元中的隨機應力分布圖

圖7為ε-LD-BCs下,各成分模量相互獨立時表征體積單元中的隨機應力分布圖.與圖6中第1行圖形比較可見,兩種相關情況下應力分布的上下邊界明顯不同,這是因為本質上兩種相關條件下求得的有效模量、有效彈性參數和有效彈性張量是不同的.故在非均質材料有效力學性質分析或在進行含有非均質材料且承受靜動態載荷的結構響應分析時,需要慎重考慮材料微觀結構中客觀存在的隨機性和相關性.

4　結束語

筆者提出了一種多尺度隨機均化方法來量化隨機微觀結構對材料隨機有效量的影響,通過多尺度有限元法和蒙特卡羅法解決了小變形下三維線彈性材料的隨機均化問題.數值算例采用二相非均質材料,當施加不同邊界條件時求解了各隨機有效量的數字特征如均方差、變異系數和相關系數等,并得出了一些有意義的結論.顯然,隨機數值均化模型是基于蒙特卡羅法發展而來的,因此其計算效率和結果的精確性與模擬次數直接相關.此外,也亟需采用其它隨機方法來解決隨機數值均化問題,從而使得兩種方法的結果可以相互驗證.三維數值算例表明:已建立的隨機多尺度方法在有效評估及優化新材料、新結構方面具有潛在的應用價值;為更準確預測失效概率,概率高階矩量信息和統計數據的獲取也非常重要,相關研究也在進行.

[1]BRUNO D,GRECO F,LUCIANO R,et al.Nonlinear Homogenized Properties of Defected Composite Materials[J]. Computers&Structures,2014,134:102-111.

[2]NEMAT-NASSER S,HORI M.Micromechanics:Overall Properties of Heterogeneous Materials[M].2nd.Amsterdam: North-Holland,1999:27-71.

[3]ZOHDI T I,WRIGGERS P.Introduction to Computational Micromechanics[M].Berlin:Springer Verlag,2005: 63-104.

[4]NGUYEN V P,LLOBERAS-VALLS O,STROEVEN M,et al.Homogenization-based Multiscale Crack Modelling: from Micro-diffusive Damage to Macro-cracks[J].Computer Methods in Applied Mechanics and Engineering,2011,200 (9-12):1220-1236.

[5]MA J,TEMIZER I,WRIGGERS P.Uncertain Analysis of the Homogenization in the Heterogeneous Material of Linear Elasticity[J].International Journal of Solids and Structures,2011,48:280-291.

[6]CAFFARELLI L A,MELLET A.Random Homogenization of an Obstacle Problem[J].Annales de l’Institut Henri Poincare(C)Non Linear Analysis,2009,26(2):375-395.

[7]FENG J W,LI C F,CEN S,et al.Statistical Reconstruction of Two-phase Random Media[J].Computers& Structures,2013,137:78-92.

[8]XU X F,CHEN X.Stochastic Homogenization of Random Elastic Multi-phase Composites and Size Quantification of Representative Volume Element[J].Mechanics of Materials,2009,41(2):174-186.

[9]LE GUENNEC Y,COTTEREAU R,CLOUTEAU D,et al.A Coupling Method for Stochastic Continuum Models at Different Scales[J].Probabilistic Engineering Mechanics,2014,37:138-147.

[10]SAKATA S,ASHIDA F,OHSUMIMOTO K.Stochastic Homogenization Analysis of a Porous Material with the Perturbation Method Considering a Microscopic Geometrical Random Variation[J].International Journal of Mechanical Sciences,2013,77:145-154.

[11]SAKATA S,ASHIDA F,FUJIWARA K.A Stochastic Homogenization Analysis for a Thermoelastic Problem of a Unidirectional Fiber-reinforced Composite Material with the Homogenization Theory[J].Journal of Thermal Stresses, 2013,36(5):405-425.

[12]HILL R,RICE J R.Elastic Potentials and the Structure of Inelastic Constitutive Laws[J].SIAM Journal on Applied Mathematics,1973,25(3):448-461.

[13]TOURAN A,WISER E P.Monte Carlo Technique with Correlated Random Variables[J].Journal of Construction Engineering and Management,1992,118(2):258-272.

[14]JIA X,WILLIAMS R A.A Packing Algorithm for Particles of Arbitrary Shapes[J].Powder Technology,2001,120 (3):175-186.

[15]TEMIZER I,ZOHDI T I.A Numerical Method for Homogenization in Non-linear Elasticity[J].Computational Mechanics,2007,40(2):281-298.

(編輯:王 瑞)

Random computational homogenization for three-dimensional linear elasticity

WANG Nan,MA Juan,ZHANG Yulin
(School of Mechano-electronic Engineering,Xidian Univ.,Xi’an 710071,China)

The computational homogenization of heterogeneous materials under infinitesimal deformation is addressed in the context of elasticity when the uncertainty in microstructure is fully considered.Based on the multi-scale finite element method and Monte-carlo method,the random homogenization model of heterogeneous materials is established when the randomness of microstructural morphology and of material properties of constituents as well as the correlation of material properties are accounted for simultaneously. The random effective quantities and their numerical characteristics as well as their correlations under different boundary conditions are then found.And the impacts of microstructural parameters on random effective quantities are also investigated and illustrated.Finally,the random stress distributions within a representative volume element under different boundary conditions are revealed as well.Obviously,it is necessary that the randomness and correlation existing in the microstructure should be fully considered during the solution of the effective mechanical properties of heterogeneous materials.

homogenization;randomness and correlation;three-dimensional linear elasticity;finite element method;Monte-Carlo method

TB33;O343.2

A

1001-2400(2016)04-0092-08

10.3969/j.issn.1001-2400.2016.04.017

2015-04-20 網絡出版時間：2015-10-21

國家自然科學基金青年基金資助項目(11102143);西安電子科技大學留學回國人員擇優資助項目(6450041101)

王 南(1990-),男,西安電子科技大學碩士研究生,E-mail:wangnan@stu.xidian.edu.cn.

網絡出版地址：http://www.cnki.net/kcms/detail/61.1076.TN.20151021.1046.034.html