随机有限元法

Yamazaki 和 Shinozuka(1987)创造性地将算子的 Neumann 级数展开式引入随机有限元 的列式工作。从本质上讲，Neumann 级数展开方法也是一类正则的小参数摄动方法，正定的 随机刚度矩阵和微小的随机扰动量是两个基本要求，这两个基本要求保证了摄动解的正则性 和收敛性，其优点在于摄动形式较简单并可以得到近似解的高阶统计量。Shinozuka 等人 （1987）将随机场函数的 Monte-Carlo 模拟与随机刚度矩阵的 Neumann 级数展开式结合,得 到具有较好计算精度和效率的一类 Neumann 随机有限元列式（称 NSFEM）。Benaroya 等（1988） 指出，将出现以随机变分原理为基础的随机有限元法来逐渐取代以摄动法为基础的随机有限 元法。Spanos 和 Ghanem 等人（1989，1991）结合随机场函数的 Karhuen-Loeve 展式和 Galerkin （迦辽金）射影方法建立了相应的随机有限元列式，并撰写了随机有限元法领域的第一本专 著《随机有限元谱方法》。
§7.2 随机有限元的控制方程[22]
从随机有限元控制方程的获得来看，随机有限元可分为 Taylor 展开法随机有限元
（TSFEM）、摄动法随机有限元(PSFEM)以及 Neumann 展开 Monte-Carlo 法随机有限元（NSFEM）。
● Taylor 展开法随机有限元
该随机有限元法的基本思路是将有限元格式中的控制量在随机变量均值点处进行
由于 TSFEM 简单明了、效率高，为我国许多学者所采用。 ● 摄动法随机有限元
摄动技术最初被用于非线性力学分析。Handa 等人成功地将一阶、二阶摄动技术用于随 机问题，给出摄动法有限元列式。该法假定基本随机变量在均值点处产生微小摄动，利用 Taylor 级数把随机变量表示为确定部分和由摄动引起的随机部分，从而将有限元控制方程 （非线性的）转化为一组线性的递推方程，求解得出位移的统计特性，进而求出应力的统计 特性。
U0
=
K
−1 0
F0
(7.2.11)
∂U ∂α i
=
K
−1 0
⎜⎜⎝⎛
∂F ∂α i
−U0
∂K ∂α i
⎟⎟⎠⎞
(7.2.12)
∂2U ∂αi∂α j
=
K0−1⎜⎜⎝⎛
∂2F ∂αi∂α
j
−U0
∂2K ∂αi∂α j
− ∂U ∂αi
wk.baidu.com
∂K ∂α j
− ∂K ∂αi
∂U ∂α j
⎟⎟⎠⎞
由上式可得位移的均值和协方差：
一阶 PSFEM 和一阶 TSFEM 一样，只需一次形成刚度矩阵、一次对刚度矩阵求逆，计算效率较
高。但 PSFEM 需以微小的摄动量为条件，一般应小于均值的 20%或 30%。
● Neumann 展开 Monte-Carlo 随机有限元
20 世纪 80 年代后期，Shinozuka 等人提出基于 Neumann 展开式的随机有限元法，使
Monte-Carlo 随机抽样，刚度矩阵只改变ΔK 项，而
P = K0-1ΔK U0 = K0-1F
(7.2.17) (7.2.18)
将式（7.2.16）代入式（7.2.1），并利用式（7.2.17）和式（7.2.18），得
U = U0-PU0+P2U0-P3U0+…
(7.2.19)
令 Ui=PiU0，则得如下的递推公式：
自 20 世纪 80 年代以来，随机有限元法已在工程结构可靠性、安全性分析领域以及在 各种随机激励下结构响应变异研究领域中得到应用，如应用于大型水利工程的重力坝、拱坝 的可靠度计算；应用于非线性瞬态响应分析；结构振动中随机阻尼对响应的影响；结构分析 的随机识别；复杂结构地震响应的随机分析和两相动力系统的随机模拟等等。随着理论研究 的深入，随机有限元将得到更加广泛的应用。
在 20 世纪 70 年代初， Cambou 首先采用一次二阶矩方法研究线弹性问题。由于这种方 法将随机变量的影响量进行 Taylor 级数展开，就称之为 Taylor 展开法随机有限元（TSFEM）。 Shinozuka 和 Astill（1972)分别独立运用摄动技术研究了随机系统的特征值问题。随后， Handa（1975）等人在考虑随机变量波动性时采用一阶和二阶摄动技术，并将这种摄动法随 机有限元成功地应用于框架结构分析。Vanmarcke 等人（1983）提出随机场的局部平均理论， 并将它引入随机有限元。局部平均理论是用随机场函数在每一个离散单元上的局部平均的随 机变量来代表该单元的统计量的近似理论。Liu W. K.等人（1986、1988）的系列工作，提 供了一种“主模态”技术，运用随机变量的特征正交化方法，将满秩的协方差矩阵变换为对 角矩阵，减少计算工作量，对摄动随机有限元法的发展做出贡献，此外，提出了一个随机变 分原理。
国内对随机有限元的研究起步较晚。吴世伟等人（1988）提出随机有限元的直接偏微
1
分法及相应的可靠度计算方法。陈虬、刘先斌等人（1989、1991）提出一种新的随机场离散 模型，建立了等参局部平均单元，并基于变分原理研究了一类随机有限元法的收敛性和误差 界。
Papadrakakis(1995)采用预处理共轭梯度法给出了空间框架的非线性随机有限元列 式。Schorling 和 Bucher(1996)基于 Monte-Carlo 技术，采用响应面法研究几何非线性时的 可靠度随机有限元方法。刘宁（1996）则基于偏微分法，给出了三维弹塑性随机有限元列式。 随机有限元法的数学理论研究和非线性随机问题的有限元分析工作还有待深入。
+
Ui′′jUk′′l E(αiα j )E(αkαl ) + E(αiαk )E(α jαl )
i=1 j=1 k=1 l=1
(7.2.15)
由于任何量的随机性都可以引入摄动量，而且更易于考虑非线性问题，因此 PSEFEM 适
3
用范围较广，对于结构几何特性的随机性（包括随机边界问题）易得出随机有限元控制方程。
结构系统的随机分析一般可分为两大类：一类是统计方法，另一类是非统计方法。因 此，随机有限元法同样也有统计逼近和非统计逼近两种类型。前者通过样本试验收集原始的 数据资料，运用概率和统计理论进行分析和整理，然后作出科学推断。这里，样本试验和数 据处理的工作量很大，随着计算机的普及和发展，数值模拟法，如蒙特卡罗（Monte Carlo） 模拟，已成为最常用的统计逼近法。后者从本质上来说是利用分析工具找出结构系统的（确 定的或随机的）输出随机信号与输入随机信号之间的关系，采用随机分析与求解系统控制方 程相结合的方法得到输出信号的各阶随机统计量的数字特征（如各阶原点矩或中心矩）。
Monte-Carlo 法与有限元法较完美地结合起来。Monte-Carlo 法是最直观、最精确、获取信
息最多、对非线性问题最有效的计算统计方法。Neumann 展开式的引入是为了解决矩阵求逆
的效率问题。如果对每一次随机抽样，只需形成刚度矩阵，进行前代、回代以及矩阵乘和矩
阵加减，而无需矩阵分解，则可大大减少工作量。
σ= DBU
(7.2.3)
设基本随机变量为 X = ( X1 , X 2 , L, X n )T ，将位移 U 在均值点 X = ( X 1, X 2 , L, X n )T
处一阶 Taylor 级数展开，并在两边同时取均值（数学期望），得
( ) E[U ] ≈ U X
=
K
−1
F
(7.2.4)
式中：符号 E[·]表示求均值，任一结点位移 U 的方差可由下式计算：
(7.2.13)
[ ] ∑∑ E U
≈U
+1 2
n i =1
n
U i′j′Cov(α i , α j )
j =1
(7.2.14)
nn
nnn
∑∑ ∑∑∑ Cov[U] ≈
Ui′U′jCov(αi , α j ) +
Ui′U′j′k E(αiα jαk )
i=1 j=1
i=1 j=1 k=1
∑∑∑∑ [ ] n n n n
− ∂K ∂Xi
∂U ∂X j
−
∂2K ∂Xi∂X
j
U
⎟⎞ ⎟⎠
(7.2.8)
∂2σ = ∂2D BU + ∂D B ∂U + ∂D + B ∂U + DB ∂2U
∂Xi∂X j ∂Xi∂X j
∂Xi ∂X j ∂X j ∂Xi
∂Xi∂X j
(7.2.9)
上式可见，二阶 TSFEM 可以放宽随机变量变异性大小的限制，但随机变量数目较多时，计算 量将十分庞大，而且一阶或二阶 TSFEM 均无法计算响应量三阶以上的统计特性。
最初是 Monte-Carlo 法与有限元法直接结合，形成独特的统计有限元方法。Astill 和 Shinozuka（1972）首先将 Monte-Carlo 法引入结构的随机有限元法分析。该法通过在计算 机上产生的样本函数来模拟系统的随机输入量的概率特征，并对于每个给定的样本点，对系 统进行确定性的有限元分析，从而得到系统的随机响应的概率特征。由于是直接建立在大量 确定性有限元计算的基础上，计算量极大，不适用于大型结构，而且最初的直接 Monte-Carlo 法还不是真正意义上的随机有限元法。但与随后的摄动随机有限元法（PSFEM）相比，当样 本容量足够大时，Monte-Carlo 有限元法的结果更可靠也更精确。
假设α i 为随机变量 X i 在均值点 X i 处的微小摄动量，即α i = X i − X i 。于是
∑ ∑ K
≈
K0
+
n
αi
i =1
∂K ∂α i
+
1 2
i,
n
α
j =1
iα
j
∂2K ∂α i∂α j
(7.2.10) 对于 U、F，也有类似上式 K 的表达式，式中：K0、U0、F0 分别为 K、U、F 在随机变量均值点 的值。根据二阶摄动法，可得
随机有限元法
§7.1 绪论
结构工程中存在诸多的不确定性因素，从结构材料性能参数到所承受的主要荷载，如 车流、阵风或地震波，无不存在随机性。在有限单元法已成为分析复杂结构的强有力的工具 和广泛使用的数值方法的今天，人们已不满足精度越来越高的确定性有限元计算，而设法用 这一强有力的工具去研究工程实践中存在的大量不确定问题。随机有限元法（Stochastic FEM）,也称概率有限元法（Probabilistic FEM）正是随机分析理论与有限元方法相结合的 产物，是在传统的有限元方法的基础上发展起来的随机的数值分析方法。
Taylor 级数展开（取一阶或二阶），经过适当的数学处理得出所需的计算方程式。有限元静
力分析控制方程的矩阵形式为: KU = F
(7.2.1)
式中,U 为位移矩阵，F 为等效节点荷载列阵，K 为整体刚度矩阵
K = ∑ ∫∫∫ΩBe T DBdv
(7.2.2)
其中，B 为形变矩阵,D 为材料弹性矩阵。在计算出节点位移 U 后，即由下式求得应力列阵σ
在一般有限元控制方程 KU = F 中，假定荷载 F 为确定值，在随机变量波动值的影响下
刚度矩阵 K 分解为 K = K0+ΔK，根据 Neumann 级数展开，有 K-1=(K0+Δ K)-1=(I-P+P2-P3+… )K0-1
(7.2.16)
式中：K0 为随机变量均值处的刚度矩阵；ΔK 为刚度矩阵的波动量；I 为单位矩阵。对于
2
但由于忽略了二阶以上的高次项，使 TSFEM 对随机变量的变异性有所限制。一般要求一阶 TSFEM 随机变量的变异系数小于 0.3。如果随机变量的变异系数较大，可以采用有限元控制 方程的二阶 Taylor 展开：
∂2U ∂Xi∂X j
=
K
−1
⎜⎛ ⎜⎝
∂2F ∂Xi∂X
j
− ∂K ∂Xi
∂U ∂X j
[ ] ∑ ∑ Var U
n
≈
n ∂U
i=1 j=1 ∂X i
⋅ ∂U X = X ∂X j
Cov( X i , X j )
X =X
式中：符号 Var[·]表示求方差；Cov(Xi,Xj)为 Xi 和 Xj 的协方差。
其中
∂U = K −1 ( ∂F − ∂K U )
∂X i
∂X i ∂X i
(7.2.5) (7.2.6)
∂σ = ∂D BU + DB ∂U
∂X i ∂X i
∂X i
(7.2.7)
同样将σ在均值点处 Taylor 展开，也有与上面类似的表达式。可见，TSFEM 关键在于 对有限元方程式直接进行偏微分计算，计算出有限元输出量对随机变量的梯度，故该法也称 直接偏微分法或梯度分析法。
由于一阶 TSFEM 只需一次形成刚度矩阵，也只需一次求刚度矩阵的逆，因此效率较高。