摘要:阐述一种适用于非均质材料力学性能分析的扩展的多尺度有限元法(Extended Multiscale Finite Element Method, EMsFEM)的基本原理. 该方法的基本思想是利用数值方法构造能反映胞体单元内部材料非均质影响的多尺度基函数,在此基础上求得粗网格层次的等效单元刚度阵,从而在粗网格尺度上对原问题进行求解,很大程度地减少计算量. 以该方法进行的具有周期和随机微观结构的材料计算示例,通过与传统有限元法的结果比较,说明这一方法的有效性. EMsFEM的优势在于,能容易地进行降尺度计算,可较准确地求得单元内部的微观应力应变信息,在非均质材料强度和非线性分析中有很大的应用潜力.
关键词:扩展的多尺度有限元法; 基函数; 非均质材料; 降尺度计算
中图分类号:O241.82; TB115文献标志码:A
Basic theory of extended multiscale finite element method
ZHANG Hongwu, WU Jingkai, LIU Hui, FU Zhendong
(Dept. of Eng. Mechanics, Faculty of Vehicle Eng. & Mechanics, State Key Lab. for Structural Analysis of
Industrial Equipment, Dalian Univ. of Tech., Dalian Liaoning 116024, China)
Abstract: The basic theory of Extended Multiscale Finite Element Method(EMsFEM) for mechanical analysis of heterogeneous materials is presented. The underlying idea is to construct numerically the multiscale base functions to capture the small scale heterogeneities of coarse elements in the multiscale finite element analysis. Then the problems are solved on the coarse-grid scale, thus resulting in a reduced number of degrees of freedom in the model. Both problems with periodic and random microstructures are considered and the numerical results verify the validity and accuracy of the developed method by comparing them with the traditional finite element method. An important feature of this work is that the downscaling computation could be performed easily and the micro stress and strain in the macro elements can be obtained simultaneously in the multiscale computation. Thus, the developed method has great potential for strength and nonlinear analysis of complicated heterogeneous materials.
Key words: extended multiscale finite element method; base function; heterogeneous material; downscaling computation
收稿日期:2010-[KG*9〗04-[KG*9〗28修回日期:2010-[KG*9〗05-[KG*9〗03
基金项目:国家自然科学基金(10721062, 50679013, 90715037, 10728205);长江学者和创新团队发展计划;
国家基础性发展规划项目(2010CB832704)
作者简介: 张洪武(1964—),男,辽宁庄河人,教授,博导,博士,研究方向为计算力学与工程科学计算等,(E-mail)zhanghw@dlut.edu.cn 0引言
自然存在和人工形成的大部分材料都具有非均质性特征,例如地下岩土以及航空航天工业中广泛使用的复合材料等,因此研究非均质材料的力学性能具有非常重要的意义.当材料结构具有多尺度特征时,将整个结构体直接离散化进行分析往往要耗费巨大的计算机资源,甚至不可行.寻求既可以节省计算资源,又可以保证计算精度的多尺度数值计算方法已成为近年来的研究热点.目前国内外学者已经提出各种多尺度计算方法,其中比较常用的有均匀化方法[1-3]和代表体元法[4-5]等,但这些方法一般建立在微观结构的周期性假设基础上,对材料强度和非线性问题等的分析还有相当多的困难.因此,对于多尺度计算,研究工作还刚刚起步,许多具有挑战性的问题有待解决.
多尺度有限元法(Multiscale Finite Element Method,MsFEM)的原始思想来自于BABUSKA等[6]的工作;HOU等[7-8]和EFENDIEV等[9]对该法的发展作出重要贡献,他们通过在每个宏观单元上求解子问题,数值构造出满足局部特性微分算子的多尺度基函数(形函数),从而在粗网格尺度上对原问题进行求解就能得到较高的精度,同时,由于各个单元间基函数的构造相互独立,故该法能很容易地进行并行计算.MsFEM自提出以来已经被广泛应用于求解具有高度振荡系数的2阶椭圆型边值问题(标量场问题),但在具有矢量场特征的固体力学计算方面还少有工作报道.
①该文已被Int J Multiscale Comput Eng录用,待发表,作者是ZHANG Hongwu, WU Jingkai和FU Zhendong.ZHANG等[10]通过分别构造固相变形场和液相压力场基函数的方法首次将MsFEM用于求解非均质饱和多孔介质的耦合固结问题.文献Extended multiscale finite element method for mechanical analysis of periodic lattice truss materials①(以下简称文献①)考虑多维矢量场问题不同方向的耦合作用,通过引入基函数的耦合附加项,首次提出扩展的多尺度有限元法(Extended MsFEM,EMsFEM),能很好地在大尺度上预测周期性桁架材料的等效力学性能.文献[11]对文献①的工作作进一步发展,利用EMsFEM的优势,建立桁架类材料非线性问题多尺度分析的EMsFEM.
本文以二维连续体问题分析为例,介绍EMsFEM的基本原理及其实施过程,并通过几个有代表性的数值算例说明其有效性和精确性.可以看到,EMsFEM在显著减少计算量和计算所占内存的同时,对周期性和非周期性材料都能求得满意的结果.同时,对于均质材料,在一致粗网格下EMsFEM计算精度明显优于传统有限元法,且随着网格加密对L2范数收敛.与均匀化方法不同的是,EMsFEM能很容易地进行降尺度计算,即在宏观位移变量获得的同时直接求得单元微观尺度上真实的应力应变信息,由此为非均质材料的强度计算提供一种新的多尺度计算方法.
限于篇幅,本文只介绍EMsFEM的一般原理,关于桁架材料和多孔介质等的分析可参见文献[10-11]和文献①,关于连续介质分析更为系统的描述将由另文给出.
1EMsFEM的基本思想
传统有限元法的基函数构造一般依据单元节点量直接用多项式进行插值,所以在同一个单元内部的材料参数,如弹性模量等必须一致.正是由于这种限制,对于多尺度问题,传统有限元法只有在微观尺度上进行细致的离散求解才能获得有意义的结果,这样就要求有大量的计算机资源,对大型复杂结构体求解变得十分困难.
EMsFEM通过在单元上求解局部子问题进行基函数的数值构造,这些基函数可以准确有效地反映材料的微观非均质性,这样就可以在宏观层次上得到准确有效的解,从而避免在微观层次上求解的问题,大大节省计算资源.EMsFEM的基函数通过数值构造,故在每个单元内部的材料属性可以是非均质的,因此对复合材料和多孔介质材料的多尺度计算有很好的应用前景.
EMsFEM已被成功地应用于非均质连续体材料和桁架材料的多尺度计算分析中,本文以连续体材料为例说明该方法的实施过程,而桁架材料分析只给出算例,具体的桁架材料多尺度推导过程参考文献①和文献[11].
考虑图1所示的二维连续体材料结构,其求解域为Ω,边界为Γ,则结构满足以下平衡方程和边界条件div(D:E(u))=f-in Ω
nσ=T-or Γσ
u=u-on Γu(1)式中:D为表示材料属性的4阶刚度张量;E(u)为应变张量;u为位移向量;σ为应力向量,σ=[σxσyτxy]T;f-为体积力向量,f-=[f-xf-y]T;Γσ和Γu分别为力的边界和位移边界,并且它们满足Γσ∩Γu=;n为边界外法线的方向余弦向量.
图 1EMsFEM示意图
EMsFEM的计算流程可分为微观计算、宏观计算以及降尺度计算等3个主要部分.如图1所示,微观计算指在子网格上数值构造宏观单元(粗网格)的多尺度基函数,该基函数能反映宏观单元内部的微观非均质性,进而求得单元的等效刚度阵;宏观计算指基于微观计算求得的宏观单元等效刚度阵,在宏观尺度上对物理问题进行求解,大大减少计算自由度;降尺度计算指在宏观计算所获结果基础上,通过快速简洁的计算,获得细尺度层次的物理力学量.
2基函数构造
文献[10-11]和文献①的研究发现,由于存在体积膨胀/收缩效应(泊松效应),固体变形时各个方向间会产生耦合作用,故将MsFEM直接应用于计算固体力学中的矢量场问题一般会有较大误差,因此需要研究新的基函数构造理论.解决该问题的方法之一是在不同的坐标方向分别构造基函数,并同时引入基函数的耦合附加项.如在二维问题的EMsFEM中,需分别构造基函数Nixx和Niyy,i=1,2,3,4为宏观单元的4个节点.这样,宏观单元内部细网格(子网格)任意一节点的位移可以表示成u=4i=1Nixxu′i+4i=1Nixyv′i(2)
v=4i=1Niyyv′i+4i=1Niyxu′i(3)式中:Niyx为基函数的耦合附加项,其物理意义是指在宏观节点i发生x方向的单位位移时,单元内各微观节点产生的y方向位移值.
式(2)和(3)可统一写成u=Nu′E(4)式中:u为宏观单元内部细网格上所有节点的位移向量;u′E为宏观单元节点的位移向量.它们可写成
u=[u1v1u2v2…unvn]T(5)
N=[RT1RT2…RTn]T(6)
u′E=[u′1v′1u′2v′2u′3v′3u′4v′4]T(7)Ri=N1xx(i)N1xy(i)N2xx(i)N2xy(i)N3xx(i)N3xy(i)N4xx(i)N4xy(i)
N1yx(i)N1yy(i)N2yx(i)N2yy(i)N3yx(i)N3yy(i)N4yx(i)N4yy(i)i=1,2,…,n(8)
式中:n为宏观单元内细网格的节点总数.
如图1所示的结构内的某个粗网格单元,其内部域为K,KΩ.数值构造单元基函数,即在特定的边界条件下求解单元内部的平衡方程LNi=0in K
Ni(x)affined on K
i=1,2,…,m(9)式中:L为弹性算子,满足Lu=divD:12(u+(u)T);Ni为粗网格节点i的基函数,满足Ni|j=Ni(xj,yj)=δij,(i,j=1,2,…,m),δ为Kronecker符号;m为粗网格单元的节点数,本文取值m=4.
对于标量场问题,基函数构造时方程(9)满足如下边界条件[7-9]Ni=N0ion K(10)式中:N0i为构造基函数Ni时所施加的边界条件.
对于矢量场问题,需在不同坐标方向上分别构造基函数.以二维问题为例,Ni由Nixx和Niyy构成.下面以N1xx为例说明基函数的构造过程.
文献[7-10]和文献①中的大量数值算例表明,基函数构造时施加的边界条件N0i的不同会对结果的精度造成很大影响.这里先以简单的线性边界为例阐述基函数的构造方法,更为精确的边界条件在第3节中详细介绍.
线性边界条件求基函数N1xx,就是在边14和12上分别加上x方向的线性边界,即由N1xx(x1,y1)=1线性地变化到N1xx(x2,y2)=0和N1xx(x4,y4)=0,图 2数值基函数的构造方法边23和34上x方向位移则等于0,同时,边界上各节点的y方向都固定,约束见图2,则在上述边界约束条件下,单元内部的Nx1={N1xx,N1yx}值可由传统数值方法(如有限元法)在子网格K上对平衡方程(9)进行求解得到.依此类推,即可求得单元全部基函数N.
可以证明,上述构造的基函数满足在单元中任一点的基函数之和等于1,即4i=1Nixx=1,4i=1Niyy=1
4i=1Niyx=0,4i=1Nixy=0(11)也即保证宏观单元的刚体位移与单元间的位移协调.
3超样本技术
很显然,基函数的边界条件应能反映宏观单元边界上解的微观振荡性.当采用线性边界条件构造基函数时,宏观单元的边界强制发生线性变形,这种人工约束会在非均质单元的边界产生较强的边界层效应[7],从而使多尺度计算的刚度偏大(见文献①).为减小这种边界层效应带来的误差,这里介绍一种超样本技术以构造能更好反映单元边界变形的振荡边界条件.超样本技术首先由HOU等[7]提出,以提高MsFEM精度及收敛速度,并避免网格大小与多孔介质物理小尺度大小相近而引起共振效应.文献①和文献[11]进一步改进上述超样本技术,使之适应于计算固体力学的矢量场问题分析.然而,由文献①和文献[11]中超样本技术构造的基函数求得的单元等效刚度阵不能完全满足刚体位移,需要修正刚度阵.本文进一步改进上述超样本技术,使求得的等效刚度阵能完全满足刚体位移.数值试验表明,由于用超样本技术构造的基函数能更真实地模拟单元边界的变形,所以在很大程度上可降低子网格边界节点力不平衡造成的误差,使结果更加准确.
图 3超样本技术构造基函数如图3所示,某个较大的域K′包含原来的子网格域K,其中Δ1234为初始宏观单元,Δ1′2′3′4′为超样本单元.首先构造超样本单元的临时基函数Ψj′(j′=1′,2′,3′,4′).以Ψ1′xx的构造为例,只需将边界2′3′和3′4′的x方向位移固定,点3′的x和y方向位移全固定以防止刚体位移,而在另外两边界加上线性边界.依此方法,可分别构造出Ψj′xx和Ψj′yy 的值,则宏观单元的临时基函数φi可通过Ψj′线性组合得到ixx=4j=1cxijΨj′xx,ixy=4j=1cyijΨj′yy(12)式中:常系数项cxij和cyij可由条件ixx|j=δij和iyy|j=δij分别求得.临时基函数φi仍然满足4i=1ixx=1,4i=1iyy=1.
在具体实施过程中,只需要临时基函数φi在宏观单元边界上的值,即最终基函数的振荡边界条件.通过振荡边界条件,采用第2节中介绍的方法对宏观单元重新进行1次有限元分析,同时考虑基函数的耦合附加项,以得到新的基函数.以Nxi={Nixx,Niyx}的构造为例,可以求解下列方程
LNxi=0in K
Nixx(x)=ixx,Niyx(x)=0on K
i=1,2,…,m(13)
如此,可求得最终的宏观单元基函数Nixx,Nixy,Niyy和Niyx,它们仍然满足式(11)的要求.
4宏观计算
宏观单元多尺度基函数构造完后即可求得宏观单元的等效刚度阵.在宏观单元内任取1个细网格单元e来考虑,假设该细网格单元的节点编号为(e1,e2,e3,e4),则平面线弹性变形体的应变能为
Πe=12uTeKeue,Ke=∫ΩeBTeDeBet dΩe(14)
式中:ue为单元e 的节点位移向量;Ke为单元e的刚度矩阵;Be和De分别为单元e的应变位移矩阵和材料弹性矩阵,t为平面单元的厚度.
由式(4)~(8)可得单元e节点位移与宏观单元节点位移之间的关系
ue=Geu′E,Ge=[RTe1RTe2RTe3RTe4]T(15)
式中:Ge为细网格单元转换矩阵,即表征细网格节点与宏观单元节点之间映射关系的矩阵.
将宏观单元内部的所有细网格单元应变能相加,可得宏观单元的等效刚度矩阵KE=pe=1K′e,K′e=GTeKeGe(16)式中:p为1个宏观单元内细网格的单元数.
在获得宏观单元等效刚度阵的基础上,通过宏观结构总体刚度阵组装,建立宏观层次求解方程,进而可以进行宏观尺度的求解.
5降尺度计算
对于EMsFEM,微观与宏观尺度之间通过数值构造的基函数进行联系,该基函数可以反映宏观单元内部的非均质性.这样,当宏观位移场求出以后,可以很容易地利用基函数进行降尺度计算,求出任意单元内部小尺度上节点的位移值,进而求得单元微观应力应变信息,为材料强度计算和非线性多尺度分析打下基础.
以第4节提到的细网格单元e为例,在得到宏观位移向量U后,可通过式(15)得到微观节点的位移ue,进而利用几何方程和物理方程获得单元e的应变应力信息εe=Teu′E,Te=BeGe(17)和σe=Seu′E,Se=DeBeGe(18)在实际程序执行中,可将单元e的应变映射矩阵Te和应力映射矩阵Se储存在数据库中,这样就能方便地获取微观应变应力信息,为以后的材料非线性分析打下基础.
6数值算例
通过几个有代表性的算例说明介绍算法的有效性.对于算例2~5,分别用EMsFEM和FEM-F方法计算并比较其结果.其中,FEM-F为传统有限元法在细网格尺度上的解,可以当作参考解.值得注意的是,对于EMsFEM,构造基函数时若所取超样本域的大小与粗网格大小一致,则超样本技术不发挥作用,此时相当于基函数用图2介绍的线性边界直接构造.以下算例中各参数和结果值采用无量纲描述,并假设连续体介质算例均为平面应变问题.
算例1均质材料悬臂梁结构分析.为说明基
图 4悬臂梁结构函数构造时引入耦合附加项的合理性,构造1个简单的均质悬臂梁算例.如图4所示,悬臂梁结构所有节点x方向位移固定,中点y方向位移固定,左端在其右端受剪力P=1 000的作用,梁长L=192,梁高h=32.结构的材料全部是均质的,材料常数为弹性模量E=0.5×107,泊松比μ=0.3.剪力P按抛物线分布载荷加载.悬臂梁平面应变问题的中性轴挠度理论表达式为va=P(1-μ2)6EI(3Lx2-x3)+Ph28IGx.
可以验证超样本技术对均质单元不起作用,此时EMsFEM基函数由图2所示的线性边界直接构造.
在本算例中,分别采用M×N=12×2,24×4,48×8等几种不同的单元网格密度划分悬臂梁求解区域,其中M和N分别为x和y方向单元剖分数.用EMsFEM构造基函数时的子网格一致采用5×5网格数.分别用传统有限元法和EMsFEM求解该悬臂梁问题,通过L2范数L2=1nni=1v-vamax(|va|)212(19)验证数值方法的实际效果,其中,n为悬臂梁中性轴上相应的粗网格节点数,结果见表1.表中FEM表示传统有限元法在相应的M×N网格密度下的结果.可以看出,对于均质材料结构EMsFEM,由于基函数考虑耦合附加项,在相同的网格密度下能得到比传统有限元法更高的精度,且收敛于精确解.
表 1L2范数值M×NFEMEMsFEM12×2 6.62E-22.48E-224×4 1.93E-27.66E-348×8 5.97E-312.96E-3
算例2含有孔洞的周期性非均质材料结构分析.考虑1个具有微观周期性分布的悬臂梁结构,结
图 5有孔结构的细尺度网格构由M×N=30×6个宏观单元组成(见图5),左端固定,右端受到1 000的均布力作用.
图 6子网格结构各个宏观单元内部的子网格结构见图6.细网格单元上材料常数与算例1一致,由于子网格结构内部有孔洞,故单元整体仍为非均质.EMsFEM的超样本域采用3×3个粗网格大小.
图 7结构中性轴上节点的
y方向位移值2种方法下结构中性轴上节点的y方向位移值见图7.EMsFEM得到的结果与精细有限元FEM-F吻合得很好.同时,对EMsFEM进行降尺度分析,以求得结构微观尺度上的应力信息,并与FEM-F的结果进行比较.如图8所示,EMsFEM求得的微观von Mises应力分布与图 82种方法求得的结构细尺度
von Mises等效应力分布FEM-F总体上吻合得较好.从图8取出1个宏观单元(A位置)进行放大,结果见图9,单胞内部的微观应力分布与FEM-F大致相同.
(a)FEM-F(b)EMsFEM图 9A位置宏观单元的细尺度von Mises等效应力分布
对于具有周期性微观结构的材料,MsFEM中宏观单元的基函数只需构造1次并可应用于所有的单元,分析非常方便.
算例3周期性桁架材料结构分析.上文已经提到,可以很方便地将EMsFEM应用于桁架材料的
注:各杆弹性模量为1E+6,边界杆件截面积为0.5,其他杆件截面积为1.
图 10桁架材料的
子网格结构多尺度计算分析中.这里给出1个简单的周期性桁架材料分析算例,具体推导过程可参考文献①和文献[11].对于算例2中的具有周期性微观结构的悬臂梁结构,结构的子网格胞元(见图6)被替换成图10所示的桁架材料胞元,而边界约束和载荷情况都与算例2一致.2种方法下结构中性轴上节点
图 11结构中性轴上节点的
y方向位移值的y方向位移值见图11,与算例2类似,EMsFEM的结果与FEM-F吻合得很好.
算例4随机非均质悬臂梁结构分析.上文已经提到,构造基函数时边界条件的选择对结果精度影响很大,故通过该数值算例说明这种影响并检验超样本技术的效果.某非均质悬臂梁结构见图12,其结构左端固定,右端受到大小为1的均布载荷.梁长L=600,梁高h=120,细网格尺寸为1×1.结构细网格上的弹性模量场在1~100之间随机分布,材料泊松比μ=0.3.对于EMsFEM,宏观单元的网格划分(M×N=50×10)见图13(a),每个宏观单元采用12×12的子网格(见图13(b))构造基函数.分别计算出超样本域为12×12和36×36情况下梁的右下端点y方向位移值为-10.863 0和-11.166 5,2种情况下相对于FEM-F(值为-11.385)的误差分别为4.58%和1.92%.可以说明,由超样本技术构造的基函数能很好地反映宏观单元边界层解的微观振荡性,使单元变形更加合理.
图 12悬臂梁结构非均质弹性模量分布场
图 13EMsFEM网格划分
用2种方法求得的结构细尺度von Mises等效应力分布(对于EMsFEM,超样本域采用36×36的情况)见图14.可以看出,对于随机非均质材料EMsFEM可以得到满意的微观等效应力,这是其他均匀化方法难以做到的.同时,相对于其他多尺度方法,EMsFEM能直接在宏观单元上求得真实的微观变量信息,而不是在宏观单元积分点(Gauss点)上对代表体元问题重新计算得到微观信息,这样能进一步大量减少计算所用内存和计算量,故在大尺度非均质材料的非线性分析中有很大的应用潜力.[11]
图 14结构细尺度的von Mises等效应力分布
算例5饱和多孔介质固结分析[10].已有的工作表明,EMsFEM同样可以运用于非均质含液多孔材料的力学分析中,这里给出用该方法求解非均质饱和多孔介质固结现象的流固耦合问题算例.对于耦合场问题,可以构造2类基函数,分别用来模拟流体的压力场和固相的位移场.考虑图15的饱和多孔介质结构,对于固相材料,结构底部边界2个位移方向都约束,顶部边界自由,左右边界约束x方向位移;而水压边界在顶部保持为0,其他边界加上非渗流条件.本算例的材料参数如下:泊松比ν=0.2;水的动力黏度为μw=1.0E-3;组成多孔介质(如土壤)的固体颗粒体积模量Ks=1.0E+12;水的体积模量Kw=2.0E+9.这里,渗透系数场和弹性模量都在细网格上随机分布,其中,渗透系数的变化范围为2.996 9E-9~3.434 1E-16,弹性模量变化范围为1~100.如图15所示,在底部中心2个细网格单元上加上汇项.2种方法下点A的沉降历史响应和点B的水压历史响应分别见图16和17,2种方法吻合得比较好.当时间t=50时,2种方法得到的y方向位移场分别见图18和19,也满足精度要求.可以看出,EMsFEM可以很好地用来求解耦合多尺度问题.图 15结构的细网格和粗网格描述图 16点A垂向沉降历史响应图 17点B水压历史响应
图 18由FEM-F得到的
细尺度y方向位移场图 19由EMsFEM得到的
宏观尺度y方向位移场
7结论
介绍EMsFEM并通过几个有代表性的数值算例说明其有效性与适用性.可以发现,EMsFEM充分发挥传统有限元法的优势.以此为基础进行的多尺度计算方法的实施,具有容易理解、实现方便的特点.从数值结果可以看出,由于在EMsFEM中基函数添加耦合附加项,并考虑单元内部位移矢量场不同方向间的耦合作用,故构造的基函数能很好地捕获宏观单元内部的非均质属性.特别地,对于均质材料,在一致粗网格下EMsFEM能得到比传统有限元法更精确的结果,并对L2范数收敛.
此外,与其他多尺度方法相比,EMsFEM的另一个优势是能很容易地进行降尺度计算,快速求得宏观单元内部真实的微观应力应变信息,故在大尺度非均质材料的强度和非线性分析中有很大的应用潜力.作为一种新的多尺度计算方法,EMsFEM还有大量工作有待深入研究与应用,包括方法精度的进一步提高,向其他需要构造类似映射关系的数值方法中的推广等.
参考文献:
[1]BENSSOUSAN A, LIONS J L, PAPANICOULAU G. Asymptotic analysis for periodic structures[M]. Amsterdam: North-Holland Publishing Co., 1978.
[2]MATSUI K, TERADA K, YUGE K. Two-scale finite element analysis of heterogeneous solids with periodic microstructures[J]. Comput & Struct, 2004, 82(7-8): 593-606.
[3]YUAN Zheng, FISH J. Toward realization of computational homogenization in practice[J]. Int J Numer Methods Eng, 2008, 73(3):361-380.
[4]PECULLAN S, GIBIANSKY L V, TORQUATO S. Scale effects on the elastic behavior of periodic and hierarchical two dimensional composites[J]. J Mech & Phys Solids, 1999, 47(7): 1509-1542.
[5]HUET C. Application of variational concepts to size effects in elastic heterogeneous bodies[J]. J Mech & Phys Solids, 1990, 38(6): 813-841.
[6]BABUSKA I, CALOZ G, OSBORN E. Special finite element methods for a class of second order elliptic problems with rough coefficients[J]. SIAM J Numer Anal, 1994, 31(4):945-981.
[7]HOU T Y, WU Xiaohui. A multiscale finite element method for elliptic problems in composite materials and porous media[J]. J Comput Phys, 1997, 134(1): 169-89.
[8]HOU T Y, WU Xiaohui, CAI Zhiqiang. Convergence of a multiscale finite element method for elliptic problems with rapidly oscillating coefficients[J]. Math Comput, 1999, 68(227): 913-943.
[9]EFENDIEV Y, HOU T Y. Multiscale finite element methods for porous media flows and their applications[J]. Appl Numer Math, 2007, 57(5-7): 577-596.
[10]ZHANG Hongwu, FU Zhendong, WU Jingkai. Coupling multiscale finite element method for consolidation analysis of heterogeneous saturated porous media[J]. Adv Water Resour, 2009, 32(2): 268-279.
[11]ZHANG Hongwu, WU Jingkai, FU Zhendong. Extended multiscale finite element method for elasto-plastic analysis of 2D periodic lattice truss materials[J]. Comput Mech, 2010, 45(6): 623-635.
(编辑于杰)
