第7章 随机有限元法
第7章随机有限元法
§7.1 绪论
结构工程中存在诸多的不确定性因素，从结构材料性能参数到所承受的主要荷载，如车流、阵风或地震波，无不存在随机性。在有限单元法已成为分析复杂结构的强有力的工具和广泛使用的数值方法的今天，人们已不满足精度越来越高的确定性有限元计算，而设法用这一强有力的工具去研究工程实践中存在的大量不确定问题。随机有限元法（Stochastic FEM）,也称概率有限元法（Probabilistic FEM）正是随机分析理论与有限元方法相结合的产物，是在传统的有限元方法的基础上发展起来的随机的数值分析方法。
最初是Monte-Carlo法与有限元法直接结合，形成独特的统计有限元方法。Astill和Shinozuka（1972）首先将Monte-Carlo法引入结构的随机有限元法分析。该法通过在计算机上产生的样本函数来模拟系统的随机输入量的概率特征，并对于每个给定的样本点，对系统进行确定性的有限元分析，从而得到系统的随机响应的概率特征。由于是直接建立在大量确定性有限元计算的基础上，计算量极大，不适用于大型结构，而且最初的直接Monte-Carlo 法还不是真正意义上的随机有限元法。但与随后的摄动随机有限元法（PSFEM）相比，当样本容量足够大时，Monte-Carlo有限元法的结果更可靠也更精确。
结构系统的随机分析一般可分为两大类：一类是统计方法，另一类是非统计方法。因此，随机有限元法同样也有统计逼近和非统计逼近两种类型。前者通过样本试验收集原始的数据资料，运用概率和统计理论进行分析和整理，然后作出科学推断。这里，样本试验和数据处理的工作量很大，随着计算机的普及和发展，数值模拟法，如蒙特卡罗（Monte Carlo）模拟，已成为最常用的统计逼近法。后者从本质上来说是利用分析工具找出结构系统的（确定的或随机的）输出随机信号与输入随机信号之间的关系，采用随机分析与求解系统控制方程相结合的方法得到输出信号的各阶随机统计量的数字特征（如各阶原点矩或中心矩）。
在20世纪70年代初，Cambou首先采用一次二阶矩方法研究线弹性问题。由于这种方法将随机变量的影响量进行Taylor级数展开，就称之为Taylor展开法随机有限元（TSFEM）。Shinozuka和Astill（1972)分别独立运用摄动技术研究了随机系统的特征值问题。随后，Handa（1975）等人在考虑随机变量波动性时采用一阶和二阶摄动技术，并将这种摄动法随机有限元成功地应用于框架结构分析。Vanmarcke等人（1983）提出随机场的局部平均理论，并将它引入随机有限元。局部平均理论是用随机场函数在每一个离散单元上的局部平均的随机变量来代表该单元的统计量的近似理论。Liu W. K.等人（1986、1988）的系列工作，提供了一种“主模态”技术，运用随机变量的特征正交化方法，将满秩的协方差矩阵变换为对角矩阵，减少计算工作量，对摄动随机有限元法的发展做出贡献，此外，提出了一个随机变分原理。
Yamazaki和Shinozuka(1987)创造性地将算子的Neumann级数展开式引入随机有限元的列式工作。从本质上讲，Neumann级数展开方法也是一类正则的小参数摄动方法，正定的随机刚度矩阵和微小的随机扰动量是两个基本要求，这两个基本要求保证了摄动解的正则性和收敛性，其优点在于摄动形式较简单并可以得到近似解的高阶统计量。Shinozuka等人（1987）将随机场函数的Monte-Carlo模拟与随机刚度矩阵的Neumann级数展开式结合,得到具有较好计算精度和效率的一类Neumann随机有限元列式（称NSFEM）。Benaroya等（1988）指出，将出现以随机变分原理为基础的随机有限元法来逐渐取代以摄动法为基础的随机有限元法。Spanos和Ghanem等人（1989，1991）结合随机场函数的Karhuen-Loeve展式和Galerkin （迦辽金）射影方法建立了相应的随机有限元列式，并撰写了随机有限元法领域的第一本专著《随机有限元谱方法》。
国内对随机有限元的研究起步较晚。吴世伟等人（1988）提出随机有限元的直接偏微
分法及相应的可靠度计算方法。陈虬、刘先斌等人（1989、1991）提出一种新的随机场离散模型，建立了等参局部平均单元，并基于变分原理研究了一类随机有限元法的收敛性和误差界。
Papadrakakis(1995)采用预处理共轭梯度法给出了空间框架的非线性随机有限元列式。Schorling 和Bucher(1996)基于Monte-Carlo 技术，采用响应面法研究几何非线性时的可靠度随机有限元方法。刘宁（1996）则基于偏微分法，给出了三维弹塑性随机有限元列式。随机有限元法的数学理论研究和非线性随机问题的有限元分析工作还有待深入。
自20世纪80年代以来，随机有限元法已在工程结构可靠性、安全性分析领域以及在各种随机激励下结构响应变异研究领域中得到应用，如应用于大型水利工程的重力坝、拱坝的可靠度计算；应用于非线性瞬态响应分析；结构振动中随机阻尼对响应的影响；结构分析的随机识别；复杂结构地震响应的随机分析和两相动力系统的随机模拟等等。随着理论研究的深入，随机有限元将得到更加广泛的应用。
§7.2 随机有限元的控制方程[22]
从随机有限元控制方程的获得来看，随机有限元可分为Taylor 展开法随机有限元（TSFEM ）、摄动法随机有限元(PSFEM)以及Neumann 展开Monte-Carlo 法随机有限元（NSFEM ）。 ● Taylor 展开法随机有限元
该随机有限元法的基本思路是将有限元格式中的控制量在随机变量均值点处进行Taylor 级数展开（取一阶或二阶），经过适当的数学处理得出所需的计算方程式。有限元静力分析控制方程的矩阵形式为: KU = F (7.2.1) 式中,U 为位移矩阵，F 为等效节点荷载列阵，K 为整体刚度矩阵
∑???Ω=e
T DBdv B K (7.2.2) 其中，B 为形变矩阵,D 为材料弹性矩阵。在计算出节点位移U 后，即由下式求得应力列阵σ
σ= DBU (7.2.3)
设基本随机变量为T n X X X X ),,,(21 =，将位移U 在均值点T n X X X X ),,,(21 =处一阶Taylor 级数展开，并在两边同时取均值（数学期望），得
[]()F K X U U E 1-=≈ (7.2.4) 式中：符号E[·]表示求均值，任一结点位移U 的方差可由下式计算：
[]),(11j i X
X n i j X X n j i X X Cov X U X U U Var ====∑∑?????≈ (7.2.5) 式中：符号Var[·]表示求方差；Cov(X i ,X j )为X i 和X j 的协方差。
其中 )(1U X K X F K X U i
i i ??-??=??- (7.2.6)
i i i X U DB BU X D X ??+??=??σ (7.2.7) 同样将σ在均值点处Taylor 展开，也有与上面类似的表达式。可见，TSFEM 关键在于对有限元方程式直接进行偏微分计算，计算出有限元输出量对随机变量的梯度，故该法也称直接偏微分法或梯度分析法。
由于一阶TSFEM 只需一次形成刚度矩阵，也只需一次求刚度矩阵的逆，因此效率较高。
但由于忽略了二阶以上的高次项，使TSFEM 对随机变量的变异性有所限制。一般要求一阶TSFEM 随机变量的变异系数小于0.3。如果随机变量的变异系数较大，可以采用有限元控制方程的二阶Taylor 展开： ???
? ?????-????-????-???=???-U X X K X U X K X U X K X X F K X X U j i j i j i j i j i 2212 (7.2.8) j
i i j j i j i j i X X U DB X U B X D X U B X D BU X X D X X ???+??+??+????+???=???222σ (7.2.9) 上式可见，二阶TSFEM 可以放宽随机变量变异性大小的限制，但随机变量数目较多时，计算量将十分庞大，而且一阶或二阶TSFEM 均无法计算响应量三阶以上的统计特性。
由于TSFEM 简单明了、效率高，为我国许多学者所采用。
● 摄动法随机有限元
摄动技术最初被用于非线性力学分析。Handa 等人成功地将一阶、二阶摄动技术用于随机问题，给出摄动法有限元列式。该法假定基本随机变量在均值点处产生微小摄动，利用Taylor 级数把随机变量表示为确定部分和由摄动引起的随机部分，从而将有限元控制方程（非线性的）转化为一组线性的递推方程，求解得出位移的统计特性，进而求出应力的统计特性。
假设i α为随机变量i X 在均值点i X 处的微小摄动量，即i i i X X -=α。于是 ∑∑==???+??+≈n
i j i n j i j i i i K K K K 121,021αααααα (7.2.10)
对于U 、F ，也有类似上式K 的表达式，式中：K 0、U 0、F 0分别为K 、U 、F 在随机变量均值点的值。根据二阶摄动法，可得
0100F K U -= (7.2.11) ???
? ????-??=??-i i i K U F K U ααα010 (7.2.12) ???? ??????-????-???-???=???-j i j i j i j
i j i U K K U K U F K U αααααααααα202102 (7.2.13) 由上式可得位移的均值和协方差：
[]∑∑==''+≈n i j i n j ij Cov U U U E 11
),(21αα (7.2.14) [][]∑∑∑∑∑∑∑∑∑=========+''''+'''+''≈n i n j n k n l l j k i l k j i kl ij
n i n j n i n j k j i jk
n k i j i j i E E E E U U E U U Cov U U U Cov 111111
111)()()()()(),(ααααααααααααα (7.2.15)
由于任何量的随机性都可以引入摄动量，而且更易于考虑非线性问题，因此PSEFEM 适
用范围较广，对于结构几何特性的随机性（包括随机边界问题）易得出随机有限元控制方程。一阶PSFEM 和一阶TSFEM 一样，只需一次形成刚度矩阵、一次对刚度矩阵求逆，计算效率较高。但PSFEM 需以微小的摄动量为条件，一般应小于均值的20%或30%。
● Neumann 展开Monte-Carlo 随机有限元
20世纪80年代后期，Shinozuka 等人提出基于Neumann 展开式的随机有限元法，使Monte-Carlo 法与有限元法较完美地结合起来。Monte-Carlo 法是最直观、最精确、获取信息最多、对非线性问题最有效的计算统计方法。Neumann 展开式的引入是为了解决矩阵求逆的效率问题。如果对每一次随机抽样，只需形成刚度矩阵，进行前代、回代以及矩阵乘和矩阵加减，而无需矩阵分解，则可大大减少工作量。
在一般有限元控制方程KU = F 中，假定荷载F 为确定值，在随机变量波动值的影响下刚度矩阵K 分解为K = K 0+ΔK ，根据Neumann 级数展开，有
K -1=(K 0+ΔK )-1=(I -P +P 2-P 3+…)K 0-1 (7.2.16)
式中：K 0为随机变量均值处的刚度矩阵；ΔK 为刚度矩阵的波动量；I 为单位矩阵。对于Monte-Carlo 随机抽样，刚度矩阵只改变ΔK 项，而
P = K 0-1ΔK (7.2.17)
U 0 = K 0-1F (7.2.18)
将式（7.2.16）代入式（7.2.1），并利用式（7.2.17）和式（7.2.18），得
U = U 0-PU 0+P 2U 0-P 3U 0+… (7.2.19)
令U i =P i U 0，则得如下的递推公式：
U = U 0+U 1+U 2+… (7.2.20)
U i =-K 0-1ΔKU i-1 (i=1,2,3,…) (7.2.21)
由式（7.2.18）求出U 0后，可由式（7.2.21）求出U i (i=1,2,3,…)。
上述三种方法中，NSFEM 可以方便地调用确定性有限元的计算程序，而TSFEM 在编程上较为复杂，PSFEM 则更为复杂。由于采用Monte-Carlo 随机模拟技术，NSFEM 不受随机变量波动范围的限制，当变异系数小于0.2时，NSFEM 与一阶TSFEM 或一阶PSFEM 精度相当；当变异系数大于0.2时，后两者已不能满足精度要求，但NSFEM 仍能得出满意的结果。 §7.3 随机场的离散模型
许多物理现象和物体系统具有随机分布特性，包括系统本身的不确定或系统的激励和响应的不确定，都可以模型化为随机空间分布的随机场或随时间分布的随机过程。随机有限元法除了必须考虑材料参数等的空间变异性，需要获得随机有限元方程列式以及解决随机算子和随机矩阵的求逆问题外，还须包含对随机场的离散处理。
● 均匀各向同性随机场的特征量
1. 随机场S(t)的均值E(S(t))为常数m.
2. 随机场S(t)的标准相关函数 )
0()()(R R ττρ= (7.3.1) 式中，随机场的协方差函数 R(τ)= Cov[S(t+τ),S(t)] (7.3.2)
对于一切t ，随机场的方差为 Var[S(t)]= R(0)= σ2 (7.3.3)
相关函数也可以用谱分解表示（即Wiener-Khintchine 变换对）： ?
????==??∞+∞-∞+∞-τωτωττωττπ
ωd G R d R G cos )()(cos )(21)( (7.3.4)
常用的标准相关函数有：非协调阶跃型、协调阶跃型、三角型、指数型、二阶AR 型、
高斯型等多种形式。
3. 方差折减系数Γ2(h)
设S h (t)是随机场S(t)的局部平均随机过程，即 ?+=h t t h dt t S h
t S )(1)( (7.3.5) 则方差 Var[S h (t)]= σ2Γ2(h) (7.3.6)
式中，方差折减系数 ?-=Γh d h h h 02
)()1(2)(ττρτ (7.3.7) 可见，方差折减系数起着使S h (t)的方差比原来S(t)的方差缩小的作用。
4. 相关距离δS
δS 可以看成是任意两个相隔距离为τ的随机变量不相关的最小距离（也称相关偏度）。s δτ≥，则不相关，否则完全相关。利用相关距离δS 便于对随机场作近似处理，其计算式为：
???
????→∞→Γ=?∞0,/)()(2,)(202ωσωττρδG d h h h s (7.3.8) 以上公式均对一维问题列出，二维、三维问题也可以类似得出。
● 随机场的中心离散
中心离散是随机场最简单的一种离散方法。该法用随机场在每个单元中心点的值来表征该随机场在每个单元的属性，因而随机场在每个单元内部都是常量，且等于它在各个单元中心的值。该法程序简单，但精度欠佳，现较少采用。
● 随机场的局部平均
该法将随机场在每个单元的属性用随机场在单元上的局部平均来表征。Baecher 、
Vanmarcke 首先提出随机场局部平均的离散方法，我国学者曾对该法进行了深入研究[23][24]。
1. 一维随机场的局部平均
设一个一维连续平稳随机场S(x)，其均值为m ，方差为σ2，则随机场在一个离散单元
[t-T/2,t+T/2]上的局部平均定义为 ?+-=)2/()2/()(1)(T t T t T dx x S T t S (7.3.9)
式中：T 是局部平均单元的长度，S T (t)称为局部平均随机场，其均值也为m ，方差按（7.3.6）
式计算，h=T ，S(t)的标准相关函数为ρ(τ)=R(τ)/σ2，无量纲方差折减系数Γ2（h ）有
如下性质：①Γ2（T ）≥0；②Γ2（0）=1；③Γ2（-T ）=Γ2（T ）。
考虑随机场在两个长度分别为T 和T ’的单元上的局部平均。如果局部平均随机场分别为S T 和S T ’，则其协方差为 ∑='Γ-'=30222)()1(2),(k k k k T T T T T T S S Cov σ (7.3.10)
2. 二维随机场的局部平均
设S(x 1，x 2)为二维连续参数连续状态随机场，A i =L 1i L 2i 为一矩形单元中心点（x 1i ，x 2i ）为中心，边平行于坐标轴x 1与x 2，且长为L 1i 与L 2i 的矩形面积，随机场在A i 内的局部平均
定义为 ??+-+-==)2/()2/()2/()
2/(21212111112222),(1),(i i i i i i i i L x L x L x L x i i Ai i dx dx x x S A x x S S (7.3.11) 如果S(x 1,x 2)是一个均匀随机场，则可用均值m 、方差σ2
及归一化协方差函数ρ（r 1,r 2）近似描述，r 1和r 2分别为沿两个方向的距离。对应的局部平均随机场可用E(S i )、Var(S i )
及互协方差Cov(S i ,S j )近似描述。E(Y i )、Var(Y i )及Cov(Y i ,Y j )可由m 、σ2及ρ(r 1，r 2)计算
获得。第i 和第j 单元局部平均S i 、S j 的互协方差可表示为 ∑∑==+Γ-=30302122212
),()()
1(41),(k l l k l k l k yj xj yi xi j i T T T T T T T T S S Cov σ (7.3.12)
式中：)3,2,1,0,(21=l k T T l k 、的约定如图7.3.1所示。 x y T 1l
T 1l
T 1l
T 1l
T 1l T 1l
T 1l
T 1l T xi T xj T y j
T y i
图7.3.1 二维局部平均单元
如果有限元的网格已划分，且单元总数为n ，随机场实际上被离散成n 个随机变量，这n 个随机变量的统计特性可由E[S i ]、Var[S i ]及互协方差Cov[S i ，S j ]反映。由于有限元网格的疏密是由应力梯度决定的，与随机场无关，当单元数很多时，随机场可另划网格，网格的疏密可由相关距离δS 决定。当然，随机场网格越密精度越高。随机场的局部平均法由于对原始数据的要求低、收敛快、精度高，是随机有限元计算中最常用的方法。
● 随机场的插值
Liu W. K.提出随机场的插值法。该法将随机场在单元内的值用单元结点处值的插值函数来表示，于是随机场的统计特性可由各单元结点处随机变量间的统计特性近似反映。
利用形函数N i (X)，随机场b(X)离散式表示为
∑==q
i i
i b X N X b 1)()( (7.3.13) 式中：X 表示空间位置；b i 为随机场在结点i 处的值（i=1，2，…，q ）；q 为单元结点数。随机场b(X)在单元内的均值和方差可表示为
∑==q
i i
i b E X N X b E 1)()())(( (7.3.14) ∑==
q j i j i j i b b Cov X N X N X b Var 1,),()()())(( (7.3.15)
随机场的插值法将原连续状态的随机场仍离散成一个连续函数，未直接计算随机场引起的单元间的相关性，只需给定随机场在各结点上的值，计算相对简单，易于考虑非线性和非均匀随机场问题。但需要已知相关函数，并且要求随机场对空间参数具有较高的连续性。 ● 随机场的加权积分法
Takada 、Shinozuka 及Deodatis 提出随机场加权积分方法。该法在单元刚度矩阵的推导过程中采用随机场在单元高斯点上的加权积分，以表征单元上的随机场。
假设单元e 的弹性模量为 )](1[)()()(0)(X f
E X E e e e += (7.3.16) 式中，E 0(e)(X)为弹性模量的均值；f (e)(X)为一维零均值均匀随机场，其值域为
-1+η≤f (e)(X)≤-1-η (0＜η＜1)
两端铰接杆单元的刚度矩阵可近似表示为
K (e)≈K 0(e)+X (e)ΔK 0(e) (7.3.17)
式中，K 0(e)为弹性模量取均值时的单元刚度矩阵，而
?=L
e e dx x
f X 0)()()( （L 为杆长） (7.3.18) 固结杆单元的刚度矩阵可表示为
K (e)≈K 0(e)+X 0(e)ΔK 0(e)+X 1(e)ΔK 1(e)+X 2(e)ΔK 2(e) (7.3.19)
其中： ??==L L
e e e e dx x
f x X dx x f x X 00)(2)(2)()(1)(,)( 式（7.3.17）实质上是在考虑弹性模量随机场的情况下，将单元刚度矩阵分解成确定性部分和随机部分。式（7.3.18）等表征随机场在各单元的均值和方差，而单元间弹性模量的相关性则由下式表示
??=)(0)(0212)(1)(21)
(2)
(11221)()(),(e L e L e e e e dx dx x f x f x x X X Cov (7.3.20)
不难看出，局部平均法是加权积分法的特例（即权系数全部相同）。由于采用加权积分，其计算精度相对较高，而且该法积分只需一次进行，刚度矩阵的波动性也由此得出，因此计算效率也较高。
● 随机场的正交展开
Spanos 和Ghanem 提出的随机场正交展开法，将材料特性参数随机场进行Karhumen-Loeve 正交展开，并由此推导出刚度矩阵的级数展开式，从而获得位移、应力的统计特性。
设随机场为 )()()(0x b
x S x S n n n n ?λ∑∞=+= (7.3.21)
其中 ?=L n
n n dx x x S b )()(1
?λ 式中：λn 、υn (x)分别为随机场S(x)相关函数的特征值和特征函数。υn (x)具有如下正交性： ?=mn n m dx x x δ??)()(（Kronecker 函数） (7.3.22) 式（7.3.21）对于任何分布的S(x)均收敛，可取至第r 阶以满足精度要求。对于一维杆单元，单元刚度矩阵可表示为 ∑=+=r n e n n
e e k b k k 0
其中 ?=L e e T n n e n dx x B x P x B x k )()()()(?λ
式中：P e (x)=D e (x)/S(x), D e (x)为单元弹性矩阵。集成整体刚度矩阵，得 ∑=+
=r n n n K b K K 0 (7.3.23) 从上式，可得位移： ∑=---+=r n n n F K K K b
I U 01
11
)( (7.3.24) 利用Neumann 展开式，可进一步得出位移的统计特性。该法关键在于获得特征值和特征函数。
§7.4 随机有限元与结构可靠度
● 结构可靠性分析的一次二阶矩法
结构的安全性、适用性和耐久性统称为结构的可靠性。可靠性的数学量度为可靠度，其定义为：在规定的条件下和规定的时间内完成预定功能的概率。“规定条件”是指正常设计、正常施工、正常使用的条件；“规定时间”是指结构的设计基准期。结构的使用时间超过基准期后，其失效的概率将增大。结构的一系列基本变量都具有不确定性，因此结构的可靠性分析属于随机性分析的范畴，随机有限元法将是十分有效的工具。
在结构可靠性分析中，结构的极限状态（包括承载能力极限状态、正常使用极限状态和条件极限状态）是通过功能函数来描述的。当有n 个随机变量影响结构可靠度时，结构的功能函数为: z=g(X 1,X 2,…,X n )，基本变量X i (i=1,2,…,n)是结构上的各种外因作用、材料性能和几何参数等。z ＞0时，结构处于可靠状态；z ＜0,则处于失效状态；z=0称结构极限状态方程（一般难以用显式表示），结构处于极限状态。如果z 的概率密度函数或概率分布函数可求得，则结构可靠度的数量指标便可基于各种状态出现的概率而确定。
若功能函数仅与两个随机变量有关（如结构抵抗破坏或变形的能力R 和荷载引起结构内力、应力、位移等效应S ），即z=g(R,S)。假设R 、S 均为正态分布，其均值和标准差分别为
S R 、和S R σσ、
，此时z=R-S 也是正态分布的随机变量，具有均值S R z -=和标准差22S R z σσσ+=。其概率密度函数为 ∞∞-???????????? ??--= z z z z f z z z 221exp 21)(σσπ (7.4.1) 其分布图示于图7.4.1，阴影部分是结构失效概率P f ，非阴影部分面积即结构的可靠度P r 。
f z (Z)P f O P r
均值点可靠区改进一次二阶矩法均值一次二阶矩法
失效区
失效边界
x
2x 1
图7.4.1 正态分布概率密度函数 图7.4.2 失效边界
在工程实践中R 和S 不一定是正态分布，但可以变换成标准正态分布。他们的统计特征量是均值μ、标准差σ、相关偏度或变异系数等，变异系数ν=σ/μ，表示随机变量相对于均值的变异。目前工程上一般用无量纲的可靠指标β来反映结构的可靠度，z z σβ/=，β越大，失效概率P f 越小，其互补的可靠度P r 就越大。
一次二阶矩法采用只需已知均值和标准差的数学模型去求解结构的可靠度。此法将功能函数z=g(X 1,X 2,…,X n )在某点用Taylor 级数展开，并近似地取一次项使极限状态方程线性化，然后求得可靠指标β。
改进的一次二阶矩法将线性化点选在失效边界上，而且选在结构最大可能失效点P *上
（如图7.4.2）。选择设计验算点P *(X i *│i=1,2,…,n)作为线性化点X 0i 时，线性化的极限状
态方程为 0)(),,,(**1
**2*1=??-+
≈∑-X i i i n i n X g X X X X X g z (7.4.2) 则z 的均值为 *1
***2*
1)(),,,(X n i i i Xi n
z X g X X X X g ∑=??-+=μμ 由于设计验算点P *选在失效边界上，有g(X 1*,X 2*…,X n *)=0，因此μz 成为 ∑-??-=n i X i i Xi z X g X 1**
)(μ
μ 但随机变量X i 互不相关时，z 的标准差σz 为 ∑∑==???
? ????=???????????? ????=n i X i i X i n i i X X i z X g X g 12112**σασσ 其中 2
112**???????????? ??????=∑=n i X i i X X i i
X i X g X g σσα
称αi 为灵敏度系数，表示第i 个随机变量对标准差的相对影响。于是可靠指标β为
∑∑==???
? ??????-==n i X i i X i n i X i i i X z z X g X g X 11***)(σαμσμβ (7.4.3) 变换上式为 ∑==--???n i i X i i i X X i
X X g 1*0)(*σβαμ 即 i X i i X i X σβαμ-=* (i=1，2，…，n) (7.4.4)
上式中，μXi ，σXi 为已知的各随机变量的均值和标准差，待求的量为X i *
和β，可迭代求解。 ● 随机有限元的一次二阶矩法
设可靠性分析中的一组基本变量X=(X 1,X 2,…,X n )T 为相互独立的正态变量（不满足时可
作变换），并进一步变换为标准正态变量Y=(Y 1,Y 2,…,Y n )T ，其中Y i =(X i -μXi )/σXi 。于是功能函数也可转换到标准正态空间，即
g(X)=g(R(X),S(X))=g(R(T -1(Y)),S(T -1(Y)))=G(Y)
失效概率P f =1-υ(β),由β的几何含义可知，β值为在标准正态空间中从原点到失效面的最短距离。在标准正态空间中概率密度是关于原点（均值点）旋转对称，并且随着到原点的距离的平方呈指数下降。在一次二阶矩法中，标准正态空间中的极限状态面（失效面）被一个到原点最小距离点处的切平面代替，也即按一阶Taylor 展开式把功能函数在设计点处
展开，使之线性化。采用迭代法可以近似确定极限状态面G(Y)=0上距原点最近的点Y *,然后
按距离公式确定结构的可靠指标。具体的迭代格式如下： ),,2,1()()()()(12)(1)()1()()()(n j Y Y G Y Y Y G Y G Y Y G Y n j Y j i Y n j j i Y j i j i i i =?????????????+-?
??=∑∑
==+ (7.4.5) )()(i T i Y Y ?=β (7.4.6) 以Mises 屈服准则为例，功能函数 g(S)=s 02-(s 2xx -s xx s yy +s 2yy +3s 2xy )
于是极限状态方程 g(S(X(Y)))=G(Y)=0 (7.4.7)
计算可靠指标β时要用到功能函数G(Y)的梯度向量，由求导链式法则 Y
Y X X X S S S g Y Y G ????????=??)()()()( (7.4.8) 由于S(X)难于用显式表达，求X X S ??/)(存在困难。有效方法是利用随机有限元一次二阶矩法计算结构的可靠指标β和设计点Y *
。步骤如下：
（1） 确定随机变量转换关系Y=Y(X)；
（2） 给定初值Y (0)=0，计算式（7.4.8）右边的三个偏导数；
（3） 按式（7.4.5）和（7.4.6）计算设计验算点坐标Y *；
（4） 用随机有限元法计算S (i)和)()/)((i Y Y Y X ??； （5） 重复（3）、（4）两步骤，直至收敛（G(Y (i))=0）。
（6） 再按式（7.4.5）和（7.4.6）计算可靠指标β值。
● 随机有限元的最大熵法
一次二阶矩法的计算量随迭代次数成倍增加，使该法的使用受到限制。而最大熵法用于结构的可靠性分析时，可根据随机变量的二阶矩来拟合概率密度曲线。因此随机有限元结合最大熵法可用于求结构响应的概率密度曲线，从而计算结构的失效概率。
熵被定义为信息的均值，信息是对个别X 值不确定性的度量。不确定性越大，熵也越大。对于一个连续随机变量，熵为 ()[]?-
=R dX X f X f S ln )(，而对于离散随机变量，熵为
[]∑=-=n i i i x f x f S 1
)(ln )(，这里，f(X)是随机变量X 的概率密度函数；f(x i )是离散点概率
函数。最大熵法通过调整概率密度函数f(X)使熵S 取得最大值，并满足约束条件：??
==R i R m dX X f X dX X f )(1)(和 ，利用最大熵法可得近似的概率密度函数f(X)。为
了求得最大值，可引入拉格朗日乘子法，把问题转化为确定拉格朗日乘子。
选用基于Neumann 级数的随机有限元法，写出递推方程。可将局部平均理论引入Neumann 随机有限元法，以进一步提高计算效率。随机有限元最大熵法计算结构失效概率的步骤如下：
(1) 随机场离散，随机变量化为随机向量，并计算协方差矩阵；
(2) 利用特征正交化技术，用一组统计独立的随机变量描述随机场；
(3) 利用Neumann 随机有限元法计算响应量的前若干阶矩；
(4) 拟合所有节点或部分节点的最大熵密度函数；
(5) 利用数值积分计算结构的失效概率。
● 算例
例1 如图7.4.3所示桁架下端受一集中力，假设各设计变量都是均匀分布的随机变量，
三根杆弹性模量E ，μE =2.0×106N/cm 2，Cov(E)=2.5%；截面积A ，μA =0.8cm 2，Cov(A)=2%；
杆2长度L 2,μL2=20cm ，Cov(L 2)=1.5%；杆1和杆3的长度均为L 1，μL1=28.284cm ，
Cov(L 1)=2.3%；外力P ，μP =4000N ；Cov(P)=10%；材料屈服应力σy ，μσy =3200N/cm 2，
Cov(σy )=7%。当极限状态方程取为 σy -S=0（S 为构件中最大应力值）时，试求该结构的失效概率。
图7.4.3 简单超静定桁架
解： 用随机有限元一次二阶矩法计算，假设由各结点坐标值变异引起的方向余弦的变异很小，形成确定矩阵，则由随机有限元递推方程组可得各结点位移分布特征值。计算得结点4的位移均值和标准差为μδ4=2.2928cm 和σδ4=0.0037cm ；各杆件中最大应力为杆元2的
应力，其均值和标准差为μσ2=2928N/cm 2和σσ2=377N/cm 2。由一次二阶矩法可计算得该结构
的失效概率为0.23（由于结构超静定，杆2失效并不意味着整个结构失效），与蒙特卡罗有限元法计算结果比较接近（蒙特卡罗取样数为40时，失效概率达0.224）。
随机有限元法在分析混凝土重力坝可靠度方面已得到较多应用。研究者采用抽样方法对重力坝施行随机自动剖分，研究其可靠指标的统计规律和收敛于真解的规律，并进行结构优
化[25]。在对碾压混凝土重力坝抗滑稳定的体系可靠度研究中，还采用刚体极限平衡法、线弹性随机有限元法和三维弹塑性随机有限元法进行计算比较。
§7.5 随机有限元动力分析方法
● 不确定结构的自振特性
假设某一梁-柱结构承受随机分布的轴向静荷载N(x)，沿长度有微小扰动，同时在端部有轴向推力),1(c p p +=c 为随机变量，其方差为2
c σ，轴向分布荷载可表示为 )()(x N N x N P +=
式中p(x)是一维、均值为零的均匀随机场，方差是σP 2，相关偏度是ΘP ，互相关函数是R PP (ξ)。结构弹性模量和质量密度也随机变化，分别表示为
()[]()[]x b m m x a E x E +=+=11(和），其中m E 和分别是弹性模量E 和质量密度m 的均值，方差分别为σ2E 和σ2
m ；随机变量a(x)和b(x)是两个一维互不相关的、均匀的、均值为零的随机场。自相关函数分别为R aa (ξ)和R bb (ξ)或等价的功率谱密度函数为S aa (ω)和S bb (ω)，相关偏度分别是ΘE 和Θm 。
任一点横向位移的有限元模式为 w(x,t)=N e q e ，
其中，单元位移列阵q e =[w 1,θ1,w 2,θ2]T ,形函数N e 取三次多项式插值。单元的动能表达为： ?=??
? ????=e l e e T e e e q M q dx t w A x m T 0221)(21 (7.5.1) 式中：??+=
Ω+=e e l l e e T e e e T e e e e dx N A N x b m dx N A N m M M M 00
))(()()(，A e 为单元横截面面积。单元的应变能为 e e T e l e q K q dx x w I x E U e ?=???
? ?????=022221)(21 (7.5.2) 式中：()
()??+=
Ω+=e e l l e xx T e xx T e xx e e e dx N N x a I E dx N I E K K K 00)()(，I 是惯性矩。 单元的外力功为 ?+=??? ??????? ??+=e l e e GD T e e e Gc T e e e e
C
q K Qq q K q F dx t w l Qx F W 0121212121 (7.5.3) 式中：()()l x xp Q l x Q l Qx dx N l l L N K dx N N K e e l l e x e T
e x e GD e
x T e x e
Gc /)()/(/,)/(,002+===??，p 为任一随机变量
据此集成结构总能量和总外力功，利用Hamilton 原理??=+-2
12
10)(t t t t C dt W dt U T δδ，即对总位移向量q 取变分，考虑到q 的任意性，可得系统运动方程，进而得到如下特征问题方程： 0])([2=++-+-q QK pK K M s G D G c （7.5.4） 式中，刚度系数)(Ω+=ij ij ij k k k ，其中ij k 是确定性分量，dx x N x N I
E k j l i ij e )()(0''''=?，而随机扰动项：?''''=Ωe
l j i ij dx x N x N x a I E k 0)()()()(。于是可以写出刚度系数、质量系数以及
几何刚度系数的均值和方差的表达式[26]。在总刚度矩阵中两个刚度系数k ij 和k rs 的互协方差，
当单元尺寸相等时，利用方差函数可作进一步简化。
利用局部平均理论，在每个单元上E 、m 、Q 的均值为零，方差为：
e p p e p p e e m m e m m e e
e E e E E e l l Q Var l l m Var l l E Var /)()(/)()(/)()(22
2
222Θ==Θ==Θ==σγσσγσσγσ (7.5.5)
式中，)()()(e p e m e E l l l γγγ、、和p m E ΘΘΘ、、分别为随机变量a(x)、b(x)和p(x)的方差函数和相关偏度。利用协方差函数写出两刚度系数或两质量系数的互协方差。
最后，特征值问题的方程可表达为
x M M x K Q K Q K P K P K K G D G D G C G G )]([)]()()([Ω+=Ω++Ω++Ω+λ (7.5.6) 或 K *x=λM *x 。相应平均问题的方程为
[]
x M x K Q K P K G D G C λ=++ (7.5.7) 利用单元刚度、质量和几何刚度之间的协方差的表达式，等效刚度矩阵K *的协方差矩阵能用独立的a(x)、b(x)和p(x)来表示。解未扰动的特征值问题（平均值问题）可得特征值的均值。对于满足小扰动情况，可推导出求两个特征值之间的协方差公式。在建立了材料特性和几何尺寸有随机扰动的梁-柱的随机特征值问题的随机有限元公式后，利用局部平均理
论，随机过程由均值、方差和相关偏度来定义。上述公式可推广到二维或三维[24]。
● 线性结构系统的瞬态响应
设结构系统的质量矩阵M 是确定性的，阻尼、刚度和外力的概率分布由广义协方差矩阵Cov(b i ,b j )（i,j=1,2,…,q ）表示，线性结构系统的运动方程为
),(),()(),()(),(t b F t b x b K t b x b C t b x M =++
(7.5.8) 式中，质量阵M 、阻尼阵C(b)和刚度阵K(b)都是n 阶方阵，位移向量x 和荷载向量F 为n 阶列阵。用Taylor 级数对随机向量b 展开，并保留二阶项，则位移向量x(b,t)关于b 的二
阶摄动式为 ∑∑==??+?+=q j i j i bibj q i i bi b b t x b t x t x t b x 1,12)(21)()(),(γγ (7.5.9)
其中，γ是小参数；b 是b 的均值；)()()(t x t x t x bibi bi 和、分别表示位移向量的均值、位移向量相对于b i 的一阶差商并在b 处取值和位移向量相对于b i 、b j 的二阶差商并在b 处取值；Δb i 是b i 与i b 之差。同样可以对C(b)、K(b)和F(b,t)写出关于i b 的二阶摄动展开式。将这些展开式代入运动方程（7.5.8）式中，归并γ0，γ和γ2
项后，得到 零阶方程：)()()()(t F t x K t x C t x
M =++ （7.5.10） 一阶方程： ),,2,1(),()()()(1q i t x F t x K t x C t x
M bi bi bi bi ==++ (7.5.11) 式中，[]),,2,1()()()(),(1q i t x K t x C t F t x F bi bi bi bi =+-= (7.5.12) 二阶方程：),(~
)(~)(~)(~2222t x x F t x K t x C t x M bi ，=++ (7.5.13)
????????? ??+++-????????? ??=∑∑==),()()()(21)(21),()(21),,(~1,1,2j i bi bi bi bi bibi bibi q j i q j i j i bibi bi b b Cov t x K t x C t x K t x C b b Cov t F t x x F ：式中 (7.5.14)
∑==q j i j i bibi b b Cov t x t x 1
,2),()(21)(~ (7.5.15)
求解递归微分方程组（7.5.10）～（7.5.15）式，可同时得到)(~)()(2t x t x t x bi 、、。为节省
计算量，可利用特征正交性将协方差矩阵Cov(b i ,b j )转化成不相关的方差对角阵Var(c i )。在得到零阶方程的振型矩阵后，可使上述递归方程解耦。
假定阻尼与刚度成正比，对应于k 阶模态的阻尼比是ξk ，则各阶解耦的递归方程组为 零阶方程：n k f x x x k k k k k k k ， ,2,122==++ωωξ (7.5.16) 一阶方程：r i n k f x x x
kci kci k kci k k kci ,,2,1;,,2,122 ===++ωωξ (7.5.17) 式中，ci kci x x 是的第k 个分量；ci kci f f 是的第k 个分量，且)(1t F f ci T ci Φ=，i=1,2,…,r 二阶方程：n k f x x x k k k k k k k ,,2,1222222 ==++ωωξ (7.5.18)
式中，22x x k 是的第k 个分量；22f f k 是的第k 个分量，且)(22t F f T Φ=。
由递归方程组（7.5.16）～（7.5.18）式解得ci x x 、和2x ，再变换成)()()(2t x t x t x ci 和、，然后得到位移、应变和应力的概率分布、均值、方差和互协方差。
● 随机有限元法的计算机实施
不确定结构系统动力响应分析的随机有限元法的计算机实施过程，包含对前面导出的一系列递归方程组的时间积分，可用Newmark-β法、Wilson-θ法的直接积分法计算。为了给出均值和方差，需要积分的方程数目共有(q+2)个。所有积分过程都用了相同的有效刚度矩阵，对运用并行计算方法十分有利，也说明随机有限元法在动力分析中的效用。
对于概率非线性系统，仍可采用隐式或显式的时域直接积分法。对均值方程求解，可用Newton-Raphson 迭代进行。结构在具有时间随机性荷载作用下的动力响应分析，具有重要的现实意义，如受地震荷载、风荷载或车辆荷载作用，在某给定时刻荷载的一些物理参数还存在不确定性，因此有必要在n 维状态空间建立结构动力行为的控制方程。假定施加于系统的激励在时间历程上是零均值的正态随机向量，则响应在时间历程上也是零均值的。利用随机有限元计算的结果，均值是二阶精度，方差是一阶精度。随机有限元法可以很方便地利用已有的有限元分析程序，因此应用前景十分广阔。
多尺度方法在复合材料力学研究中的进展
多尺度方法在复合材料力学分析中的研究进展 摘要简要介绍了多尺度方法的分量及其适用范围，详细论述了多尺度分析方法在纤维增强复合材料弹性、塑性等力学性能中的研究进展，最后对多尺度分析方法的前景进行了展望。 关键词多尺度分析方法，复合材料，力学性能，细观力学，均匀化理论 1 引言 多尺度科学是一门研究不同长度尺度或时间尺度相互耦合现象的跨学科科学，是复杂系统的重要分支之一，具有丰富的科学内涵和研究价值。多尺度现象并存于生活的很多方面，它涵盖了许多领域。如介观、微观个宏观等多个物理、力学及其耦合领域[1]。空间和时间上的多尺度现象是材料科学中材料变形和失效的固有现象。 多尺度分析方法是考虑空间和时间的跨尺度与跨层次特征，并将相关尺度耦合的新方法，是求解各种复杂的计算材料科学和工程问题的重要方法和技术。对于求解与尺度相关的各种不连续问题。复合材料和异构材料的性能模拟问题，以及需要考虑材料微观或纳观物理特性，品格位错等问题，多尺度方法相当有效。 复合材料是由两种或者两种以上具有不同物理、化学性质的材料，以微观、介观或宏观等不同的结构尺度与层次，经过复杂的空间组合而形成的一个多相材料系统[2]。复合材料作为一种新型材料，由于具有较高的比强度和比刚度、低密度、强耐腐蚀性、低蠕变、高温下强度保持率高以及生物相容性好等一系列优点，越来越受到土木工程和航空航天工业等领域的重视。 复合材料是一种多相材料，其力学性能和失效机制不仅与宏观性能（如边界条件、载荷和约束等）有关，也与组分相的性能、增强相的形状、分布以及增强相与基体之间的界面特性等细观特征密切相关，为了优化复合材料和更好地开发利用复合材料，必须掌握其细观结构对材料宏观性能的影响，即应研究多尺度效应的影响。 如何建立起复合材料的有效性能和组分性能以及微观结构组织参数之间的
扩展有限元简介
扩展有限元 有限元是将一个物理实体模型离散成一组有限的相互连接的单元组合体, 该方法在考虑物体内部存在缺陷时间，单元边界与几何界面一致，会造成局部网格加密，其余区域稀疏的非均匀网格分布，在网格单元中最小的尺寸会增加计算成本，再者裂纹的扩展路径必须预先给定只能沿着单元边界发展。 1999年，美国西北大学Beleytachko 提出了扩展有限法，该方法是对传统有限元法进行了重大改进。扩展有限元法的核心思想是用扩充带有不连续性质的形函数来代表计算区域内的间断，在计算过程中，不连续场的描述完全独立于网格边界，在处理断裂问题有较好的优越性。利用扩展有限元，可以方便的模拟裂纹的任意路径，还可以模拟带有孔洞和夹杂的非均质材料。 扩展有限元是以标准有限元的理论为框架，保留传统有限元的优点，目前商业软件中如Abaqus 等都加入扩展有限元的分析模块。 扩展有限元以有限元为基本框架，主要针对不连续问题进行研究，相对于传统有限元方法，它克服了裂纹扩展问题的不足。其采用节点扩展函数，其中包括2个函数:裂纹尖端附近渐进函数表示裂纹尖端附近的应力奇异性;间断函数表示裂纹面处位移跳跃性。整体划分位移函数表示为 αααI =I I I =∑∑++=b x F a x H u x N x u N i )(])()[()('41 1 式中:)(x N I 为常用的节点位移函数；I u 为常规形状函数节点自由度，适用于模型中的所有节点；)(x H 为沿裂纹面间断跳跃函数；I a 为节点扩展自由度向量，这项只对形函数被裂纹切开的单元节点有效；)(x F α为裂纹尖端应力渐进函数；αI b 为节点扩展自由度向量，这项只对形函数被裂纹尖端切开的单元节点有效。 沿裂纹面间断跳跃函数)(x H 表达式为： otherwise n x x if x H 0)(11)(*≥-???-= 式中:x 为样本点；*x 距x 最近点；n 为单位外法线向量。 各向同性材料的裂纹尖端渐进函数)(x F α表达式为： ????? ?=2cos sin ,2sin sin ,2cos ,2sin )(θθθθθθαr r r r x F 裂纹尖端的渐进函数并不局限于各向同性弹性材料的裂纹建模。可用于弹塑性指数硬化材料，不同的裂纹尖端渐进函数的形式与裂纹位置、非线性材料变形程度有关。
多尺度有限元法在地下水模拟中的应用(1)
2004年7月 水 利 学 报SH UI LI X UE BAO 第7期 收稿日期:2003212215 基金项目:国家自然科学基金(40172082);国家自然科学重点基金(40335045);博士点基金(0284002) 作者简介:薛禹群(1931-),男,江苏无锡人,科学院院士,主要研究方向为地下水模拟。 文章编号:055929350(2004)0720007207多尺度有限元法在地下水模拟中的应用 薛禹群1,2,叶淑君1,谢春红3,张云1 (11南京大学地球科学系,江苏南京 210093;21南京大学污染控制与资源化研究国家重点实验室,江苏南京 210093; 31南京大学数学系,江苏南京 210093) 摘要:本文详细介绍了多尺度有限元法的基本原理,并将其应用于非均质多孔介质中的流动问题,对水文地质参数按函数连续变化、渐变和突变3种非均质多孔介质中的二维地下水稳定流、非稳定流分别用多尺度有限元法和传统有限元法进行了计算。计算结果的对比表明,多尺度有限元法比传统有限元法有效,既节省计算量又有较高的精度。 关键词:多尺度有限元法;非均质;多孔介质;地下水数值模拟 中图分类号:P641文献标识码:A 地下水含水系统大多是非均质的,应用传统有限元法或有限差分法解这类非均质介质中的水流问题,需要非常精细的剖分才能获得较为精确的解。若计算区面积很大,精细剖分产生的结点数过多,往往会超出普通计算机的容量,或者使计算时间过长,所以,人们一直致力于寻求某种方法[1～9],既可以减少剖分的单元数又可以保证计算精度,多尺度有限元法MsFE M (Multiscale Finite E lement Method )就可以很好实现这一目标。求解椭圆型问题的多尺度有限元方法由H ou 等[1,2]提出,它无需在小尺度上精确求解就能正确抓住解的大尺度特征,它通过基函数满足局部微分算子来实现。多尺度有限元法是数学家提出用来求解椭圆型问题的新的数值方法,在水文地质领域的应用国内外基本上未见报道。作者应用后发现,这种方法不仅对椭圆型问题(稳定流问题)有效,对抛物型问题(即非稳定流问题)也适用,在解决非均质多孔介质中的水流问题时具有既节省计算量又保证计算精度等优点。本文首先介绍多尺度有限元法的基本原理,然后对参数连续变化、参数渐变和参数突变3种非均质条件下的二维地下水稳定流、非稳定流数学模型分别应用多尺度有限元法和传统有限元法求解,将两种方法的计算结果与解析解或采用精细剖分的有限元法的数值解进行比较。H ou 等[1] 只谈及MsFE M 在参数连续变化的简单二维稳定流问题中的应用,没有涉及参数渐变、参数突变以及非稳定流问题。1 多尺度有限元法基本原理 众所周知,应用传统有限元法进行地下水模拟时通常包括以下4个步骤:(1)将研究区离散成若干个单元,假定每个单元上介质是均质;(2)采用多项式插值来表示单元内的水头分布;(3)应用迦辽金或里兹等方法建立有限元方程,再将它们集合形成整个研究区的代数方程组;(4)解该代数方程组得到各结点的水头。其中第2步常采用的是线性多项式,即单元内水头呈线性分布,这是对单元内真实流场的一种近似。如果单元内介质为非均质,则误差将增大,因此传统有限元法要求单元内介质是均质的。若含水介质非均质性明显,则必须将单元剖分得很小,以便保证每个单元的参数是常数。由此不难看出,
ABAQUS中扩展有限元(XFEM)功能简介
ABAQUS中扩展有限元（XFEM）功能简介 扩展有限元（Extended Finite Element Method）是一种解决断裂力学问题的新的有限元方法，其理论最早于1999年，由美国西北大学的教授Belyschko和Black首次提出，主要是采用独立于网格剖分的思想解决有限元中的裂纹扩展问题，在保留传统有限元所有优点的同时，并不需要对结构内部存在的裂纹等缺陷进行网格划分。 ABAQUS基于在非线性方面的突出优势，在其6.9的版本中开始加入了扩展有限元功能，到6.13做了一些修正，加入了一些可以被CAE支持的关键字。目前为止，除了手动编程，能够实现扩展有限元常用的商业软件只有ABAQUS，今天，我们就来谈谈ABAQUS 中如何实现扩展有限元。 1. XFEM理论 在XFEM理论出现之前，所有对裂纹的静态模拟（断裂）都基本上是采用预留裂缝缺角，通过细化网格仿真裂缝的轮廓。而动态的模拟（损伤）基本上都是基于统计原理的Paris 方法。然而，断裂和损伤的结合问题却一直没有得到有效的解决，究其原因，在于断裂力学认可裂纹尖端的应力奇异现象（就是在靠近裂尖的区域应力值会变无穷大），并且尽可能的绕开这个区域。而损伤力学又没有办法回避这个问题（裂纹都是从尖端开裂的）。 从理论上讲，其实单元内部的位移函数（形函数）可以是任意形状的，但大多数的计算软件都采用了多项式或者插值多项式作为手段来描述单元内部的位移场，这是因为采用这种方法更加便于在编程中进行处理。但是这种方法的缺点就是，由于形函数的连续性，导致单元内部不可能存在间断。直到Belytschko提出采用水平集函数作为手段，其基本形式为 和 上面左边的等式描述了单元内裂缝的位置，右边的等式描述了裂尖的位置。与之对应的形函数便是
裂纹扩展的扩展有限元(xfem)模拟实例详解
基于ABAQUS 扩展有限元的裂纹模拟 化工过程机械622080706010 李建 1 引言 1.1 ABAQUS 断裂力学问题模拟方法 在abaqus中求解断裂问题有两种方法（途径）：一种是基于经典断裂力学的模型；一种是基于损伤力学的模型。 断裂力学模型就是基于线弹性断裂力学及其基础上发展的弹塑性断裂力学等。如果不考虑裂纹的扩展，abaqus可采用seam型裂纹来分析(也可以不建seam，如notch型裂纹)，这就是基于断裂力学的方法。这种方法可以计算裂纹的应力强度因子，J积分及T-应力等。 损伤力学模型是指基于损伤力学发展而来的方法，单元在达到失效的条件后，刚度不断折减，并可能达到完全失效，最后形成断裂带。这两个模型是为解决不同的问题而提出来的，当然他们所处理的问题也有交叉的地方。 1.2 ABAQUS 裂纹扩展数值模拟方法 考虑模拟裂纹扩展，目前abaqus有两种技术：一种是基于debond的技术（包括VCCT）；一种是基于cohesive技术。 debond即节点松绑，或者称为节点释放，当满足一定得释放条件后（COD 等，目前abaqus提供了5种断裂准则），节点释放即裂纹扩展，采用这种方法时也可以计算出围线积分。 cohesive有人把它译为粘聚区模型，或带屈曲模型，多用于模拟film、裂纹扩展及复合材料层间开裂等。cohesive模型属于损伤力学模型，最先由Barenblatt 引入，使用拉伸-张开法则（traction-separation law）来模拟原子晶格的减聚力。这样就避免了裂纹尖端的奇异性。Cohesive 模型与有限元方法结合首先被用于混凝土计算和模拟，后来也被引入金属及复合材料。Cohesive界面单元要服从cohesive 分离法则，法则范围可包括粘塑性、粘弹性、破裂、纤维断裂、动力学失效及循环载荷失效等行为。 此外，从abaqus6.9版本开始还引入了扩展有限元法（XFEM），它既可以模拟静态裂纹，计算应力强度因子和J积分等参量，也可以模拟裂纹的开裂过程。被誉为最具有前途的裂纹数值模拟方法。本文将利用abaqus6.9版本中的扩展有限元法功能模拟常见的Ⅰ型裂纹的扩展。 2 Ⅰ型裂纹的扩展有限元分析 本文针对断裂力学中的平面Ⅰ型裂纹扩展问题用abaqus中的扩展有限元方法进行数值模拟，获得了裂纹扩展的整个过程，裂尖单元的应力变化曲线，以及裂纹尖端塑性区的形状。在此基础上绘制裂纹扩展的能量历史曲线变化趋势图。
某大跨桥梁结构一致多尺度有限元模拟
某大跨桥梁结构一致多尺度有限元模拟分析 □李晶 【摘要】大跨桥梁结构具有结构构件多、自由度数目大、连接条件复杂等特点。建立与实际桥梁结构几何构造完全一致的三维实体有限元模型操作繁复，且并无必要，这就需要建立既能准确、清晰地反映结构局部细节特性，同时又对结构整体响应没有影响的有限元模型。 【关键词】大跨桥梁结构；多尺度模拟；Ansys；有限元；耦合 【作者单位】李晶，广州番禺职业技术学院建筑工程系 结构行为一致多尺度模拟就是在同一个结构中同时使用多种尺度下的不同单元所建立的有限元模型，其中心思想是：对于结构中需要重点关注的局部构件或细节部位采用“小尺度”下的精细模型，而其余部分仍采用传统的高度简化的结构“大尺度”模型。使用多尺度模型可在结构建模和分析过程中充分考虑结构最不利部位的缺陷和其演化过程以及对结构整体响应的影响，同时实现结构整体和局部细节受力计算与应力分析。 一、大跨桥梁的一致多尺度模拟 大跨桥梁结构具有结构构件多、自由度数目大、连接条件复杂等特点。建立与实际桥梁结构几何构造完全一致的三维实体有限元模型操作繁复，且并无必要，这就需要建立既能准确、清晰地反映结构局部细节特性，同时又对结构整体响应没有影响的有限元模型。 大跨桥梁结构的有限元模拟应该是以具体的有限元分析为目标的结构行为一致多尺度模拟。针对不同的有限元分析目标，对相应的有限元分析模型也有着不同的技术要求，理应根据具体情况采取相应的有限元建模方法和策略。例如：基于桥梁设计的有限元分析只需要保证计算结果是趋于保守的就能达成设计目标，这种情况下建立高度简化的桥梁有限元模型即可满足要求。相对的，针对结构局部细节损伤或状态评估的有限元模拟，则需要建立更高精度的有限元模型，方能清晰、正确的反映出局部细节处的损伤演化等过程，否则，将会“失之毫厘、谬以千里”，从而无法实现结构损失分析和状态评估这一目标。 二、结构一致多尺度模拟的具体方法 结构一致多尺度模拟的指导思想是：针对结构中需要重点关注的关键细节部位采取小尺度建模，而对结构其他部分仍沿用传统宏观尺度结构模拟。根据大跨结构的构造特点，可以根据具体分析需要来确定所需“嵌入”的局部细节模型部位及数量。总的来说，结构行为一致多尺度模拟的关键即不同尺度的模型耦合。具体步骤总结如下： 步骤一：建立大跨结构的全尺度模型。在设计荷载或由结构健康监测系统记录的运营荷载作用下计算结构主要构件的内力，并根据计算结果确定结构关键构件和危险部位以作为结构的重点关注构件。 步骤二：建立重点关注构件的构件尺度模型。对构件进行名义应力分析，并由此确定构件中的关键焊连接细节部位，在此基础上，结合焊连接部位的细节构造以及局部几何、材料特性，建立局部细节模型。 步骤三：耦合以上三种尺度的模型。即应用某种方式将局部焊连接细节模型“嵌入”构件尺度模型或结构全尺度模型，由此可以在结构荷载的作用下，得出“热点”应力和局部损伤演化及其对构件名义应力乃至结构内力的影响。 三、某大跨桥梁结构的一致多尺度模拟 本文基于关键局部构件的受力特点分析，采用结构全尺度和局部构件尺度，对某大跨桥梁进行结构行为一致多尺度模拟。该大跨桥梁全桥孔跨布置为20+256+20+16= 312m，其中主跨为256m中承式钢管混凝土拱，边跨为2孔20m和1孔16m钢筋混凝土简支T梁，全桥桥面连续，在梁端与桥台接缝处设置伸缩缝。由于主拱拱肋采用钢管混凝土组合构件，具有钢管混凝土的优点，施工难度小，可靠性好。 （一）全桥“大尺度”模型的建立。本文采用大型通用有限元分析软件ANSYS针对某大跨桥梁结构建立空间有限元力学模型。模型建立主要包括：钢管混凝土拱肋模拟；横撑、立柱的模拟；吊杆的模拟；桥面系模拟；边界条件的处理；选取坐标系、定义各单元截面；建立全桥“大尺度”模型；局部“小尺度”模型的建立。 （二）选定局部分析部位。对于大跨度拱桥这种复杂结构，一些受力较为复杂的构件应在整体分析的基础上进行局部分析或细部分析。拱脚是整桥结构强度的关键部位之一，可能承受着自重、二期恒载、预应力和活载的作用，特别是施工阶段，拱桥结构体系和荷载状态的不断变化，结构的内力和变形也将随着不断的发生变化，尤其是拱脚在混凝土灌注时可能成为悬臂构件的固定端，受力复杂且量值很大，所以拱脚在各个状况下应力分布的有限元计算值和实验值的确定就极其重要。本文选取拱脚作为局部应力分析的关键部位并选用实体单元建立构件尺度下的“小尺度”模型。 （三）局部模型的“嵌入”。由于全桥大尺度模型中，拱脚被简化为拱肋端部的固端约束，这里首先解除拱肋端部的全部约束，并在新建的拱脚模型底部施加固定约束。将实体单 · 56 ·
扩展有限元方法和裂纹扩展
扩展有限元方法和裂纹扩展 1.1 扩展有限元方法（XFEM ）基本理论 1999年，美国Northwestern University 的Belytschko 和Black 领导的研究小 组提出了扩展有限元方法，为解决裂纹这类强不连续问题带来了曙光。他们正式 应用扩展有限元法（XFEM ）这一专业术语是在2000年，截止到目前，扩展有 限元法（XFEM ）成为我们解决强不连续力学问题的最有效的数值计算方法，也 成为计算断裂力学的重要分支。XFEM 在有限元的框架下进行求解，无需对构件 内部的物理界面进行网格划分，具有常规有限元方法的所有优点。它最明显的特 点是用已知的特征函数作为形函数来使传统有限元的位移得到逼近，进而克服了 在裂纹尖端和变形集中处进行高密度网络划分产生的困难，方便地模拟裂纹的任 意路径，而且计算精度和效率得到了显著的提高[6]。 扩展有限元方法是将已知解析解的特征函数作为插值函数增强传统有限元 的位移逼近，来使得单元内的真实位移特性得以体现，裂纹尖端和物理或几何界 面独立于有限元网格。XFEM 主要包括以下三部分内容：首先是不考虑构件的任 何内部细节，按照构件的几何外形尺寸生成有限元网格；其次，采用水平集方法 跟踪裂纹的实际位置；根据已知解，改进影响区域的单元的形函数，来反映裂纹 的扩展。最后通过引入不连续位移模式来表示不连续几何界面的演化。因为改进 的插值函数在单元内部具有单元分解的特性，其刚度矩阵的特点与常规有限元法 的刚度矩阵特性保持一致。单元分解法(Partition Of Unity Method)和水平集法 （Level Set Method ）、节点扩展函数构成了扩展有限元法的基本理论，其中，单 元分解法是通过引入加强函数计算平面裂纹扩展问题，保证了XFEM 的收敛性； 水平集法是跟踪裂纹的位置和模拟裂纹扩展的常用数值方法，任何内部几何界面 位置都可用它的零水平集函数来表示。 (1)单元分解法的基本思想是任意函数()x φ都可以用子域内一组局部函数 ()()x x N I ?表示，满足如下等式： ()()()x x N x I I ?φ∑= （1） 其中，它们满足单位分解条件：f I I ?x ()=1 ()x N I 是有限元法中的形函数，根 据上述理论，便可以根据需要对有限元的形函数进行改进。在XFEM 中，单元 分解的目的是进行数值积分，达到不引人额外的自由度的目的[7-8]。 (2)水平集法 使用水平集法来描述几何间断性。在一般情形下，多用来追踪
多尺度有限元法在地下水模拟中的应用
!""#年$月水 利学报%&’()(*’+,-.第$!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!期收稿日期：!""/01!012 基金项目：国家自然科学基金（#"1$!"3!）；国家自然科学重点基金（#"//2"#2）；博士点基金（"!3#""!） 作者简介：薛禹群（14/15），男，江苏无锡人，科学院院士，主要研究方向为地下水模拟。 文章编号："22404/2"（!""#）"$0"""$0"$ 多尺度有限元法在地下水模拟中的应用 薛禹群1，!，叶淑君1，谢春红/，张云1 （16南京大学地球科学系，江苏南京!1""4/；!6南京大学污染控制与资源化研究国家重点实验室， 江苏南京!1""4/； /6南京大学数学系，江苏南京!1""4/）摘要：本文详细介绍了多尺度有限元法的基本原理，并将其应用于非均质多孔介质中的流动问题，对水文地质参数按函数连续变化、渐变和突变/种非均质多孔介质中的二维地下水稳定流、非稳定流分别用多尺度有限元法和传统有限元法进行了计算。计算结果的对比表明，多尺度有限元法比传统有限元法有效，既节省计算量又有较高的精度。 关键词：多尺度有限元法；非均质；多孔介质；地下水数值模拟 中图分类号：78#1文献标识码：- 地下水含水系统大多是非均质的，应用传统有限元法或有限差分法解这类非均质介质中的水流问题，需要非常精细的剖分才能获得较为精确的解。若计算区面积很大，精细剖分产生的结点数过多，往 往会超出普通计算机的容量，或者使计算时间过长，所以，人们一直致力于寻求某种方法［194］，既可以减 少剖分的单元数又可以保证计算精度，多尺度有限元法:;<+:（:=>?@;AB>C <@D@?C +>CECD?:C?FGH ）就可以 很好实现这一目标。求解椭圆型问题的多尺度有限元方法由&G=等［1，!］提出，它无需在小尺度上精确 求解就能正确抓住解的大尺度特征，它通过基函数满足局部微分算子来实现。多尺度有限元法是数学家提出用来求解椭圆型问题的新的数值方法，在水文地质领域的应用国内外基本上未见报道。作者应用后发现，这种方法不仅对椭圆型问题（稳定流问题）有效，对抛物型问题（即非稳定流问题）也适用，在解决非均质多孔介质中的水流问题时具有既节省计算量又保证计算精度等优点。本文首先介绍多尺度有限元法的基本原理，然后对参数连续变化、参数渐变和参数突变/种非均质条件下的二维地下水稳定流、非稳定流数学模型分别应用多尺度有限元法和传统有限元法求解，将两种方法的计算结果与解析解 或采用精细剖分的有限元法的数值解进行比较。&G=等［1］只谈及:;<+:在参数连续变化的简单二维 稳定流问题中的应用，没有涉及参数渐变、参数突变以及非稳定流问题。 1多尺度有限元法基本原理 众所周知，应用传统有限元法进行地下水模拟时通常包括以下#个步骤：（1）将研究区离散成若干个单元，假定每个单元上介质是均质；（!）采用多项式插值来表示单元内的水头分布；（/）应用迦辽金或 里兹等方法建立有限元方程，再将它们集合形成整个研究区的代数方程组；（#） 解该代数方程组得到各结点的水头。其中第!步常采用的是线性多项式，即单元内水头呈线性分布，这是对单元内真实流场的一种近似。如果单元内介质为非均质，则误差将增大，因此传统有限元法要求单元内介质是均质的。若含水介质非均质性明显，则必须将单元剖分得很小，以便保证每个单元的参数是常数。由此不难看出，—$— 万方数据
多尺度方法应用
多尺度方法 1.多尺度方法的意义 很多自然科学和工程的问题都具有多尺度的特征。例如，高雷诺湍流的涡有大小不同的尺度，材料的微损伤有大小不同的尺度，多孔介质的孔径大小存在着不同的尺度等。然而，在实际应用中却常常忽略多尺度特征而采用经验模型。这些模型在应用中取得很大的成功，但经验模型也存在本身的局限性，主要体现在：（1）由于模型的误差大，导致很多问题求解的精度不高； （2）完全忽略细观结构的影响，不能完全反映问题本身的自然特征； （3）缺乏可靠的理论基础。 因此，对于很多问题，需要建立能反映自然属性、精度更高且具有理论基础的多尺度模型。在建立多尺度模型的同时，首先必须考虑问题自身的特征。按照问题的特征可以把多尺度问题分为以下几类： 第一类：这类多尺度问题包含了孤立的瑕点或奇异点，比如裂痕、断层、突变以及接触线。对于这类问题，只需要在孤立的瑕点火奇异点附近建立细观尺度的模型，其它区域满足某个宏观模型即可。这样细观尺度的模型只需在很小的计算区域里求解。 第二类：这类多尺度问题存在相关的宏观模型，但宏观模型不清晰，不能直接用于求解。典型的一个例子是均匀化问题，这时系数aε(x)=a(x,xε?)，其中ε表示细观尺度，虽然与宏观变量x相关的宏观模型确实存在，但宏观模型不明确。 第三类：这类问题是包含第一类和第二类特征的多尺度问题。 第四类：这类多尺度问题的习惯结构具有强烈的不规则性，难以找到相关的宏观模型。 随着多尺度模型的发展，还会出现更多类型的多尺度问题，对各类多尺度问题的求解引起了人们广泛的关注，也推动了多尺度计算方法的发展。很多科学和工程问题都存在多尺度问题，多尺度模拟是一个典型的跨学科问题，它涉及到数学、化学、物理、工程、计算机科学、环境科学等学科，越来越受到科学家的重视。目前为止，已经有一些经典的多尺度计算方法，如多重网格方法、均匀化方
扩展有限元的ABAQUS实现
扩展有限元的ABAQUS实现绪论 常规有限元方法（CEFM）和其他数值方法相比，具有一些无法比拟的优点，但仍存在一些缺陷。比如在解决类似裂纹这样的强不连续问题，由于裂纹尖端处的应力奇异性，导致计算量巨大而且精度不高。然而扩展有限元方法（extended finite element method，XFEM）的出现，和常规有限元方法相比具有显著的优势，使得我们可以在裂尖和应力、变形集中处划分高密度的网格，也可以方便的模拟裂纹的扩展，使计算量不那么巨大，保留了常规有限元法的所有优点。因此，扩展有限元得到了快速发展和应用，而且在裂纹的扩展研究中要的意义。本文开展对扩展有限元方法和裂纹问题的研究，并且基于限元ABAQUS平台，对扩展有限元方法针对裂纹扩展问题进行模拟实现。 21世纪以来，计算机硬件和数值仿真的快速发展以及工业工程实践与科学研究中存在的大量运算需求，世界上涌现出一批大型科研运算及科学模拟软件，能够极大的简化运算问题以及计算机模拟实验，使我们能够更加方便地研究虚拟工程及相关科学问题。有限元方法的出现为数值分析方法的研究带来了新的曙光，力学学科本来就是连接理工学科的桥梁，计算力学是目前力学发展的一个重要分支。有限元软件则是我们到达工程科学领域彼岸的非常重要的工具和桥梁之一。 ABAQUS软件是世界上最强大的大型有限元计算分析软件之一，具有不同种类的单元类型、材料类型和不同的分析过程，拥有很好的计算功能和模拟性能。ABAQUS软件不但可以进行一种部件和复杂物理场的分析，而且可以处理多系统的部件分析；不仅可以分析简单的线弹性问题，还可以处理复杂的非线性组合问题等，相比其它软件具有无可比拟的优势[1,2]。 固体力学中存在的两类不连续问题之一则是因为物体内部几何结构突变引起的强不连续问题，裂纹问题就是这类问题的代表。由于几何界面处的位移不连续性和裂纹尖端的应力奇异性使得这类问题的处理变得比较复杂。有限元方法、无单元方法、边界元方法等是解决不连续问题的重要的数值方法[3]。其中有限元法具有显著的优点，它可以解决各种几何形状、不同的边界条件以及材料的非线性问题和各向异性问题等，然而在处理裂纹这类强不连续问题中，传统的有限元方法有着比较大的弊端，不仅需要将裂纹几何界面设置为单元的边、裂尖设置为单元节点，而且要在裂尖和变形集中的地方进行高密度的网格划分，这些对于求解裂纹问题造成了极大的不便，使得问题变得愈来愈复杂。然而1999年，由美国西北大学以Belyschko代表的研究小组提出的扩展有限元方法（XFEM）为解决裂