第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）提出随机有限元的直接偏微
1
分法及相应的可靠度计算方法。陈虬、刘先斌等人（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为整体刚度矩阵
K?????BTDBdv (7.2.2)
?e其中，B为形变矩阵,D为材料弹性矩阵。在计算出节点位移U后，即由下式求得应力列阵σ
σ= DBU (7.2.3) 设基本随机变量为X?(X1,X2,?,Xn)T，将位移U在均值点X?(X1,X2,?,Xn)T处一阶Taylor级数展开，并在两边同时取均值（数学期望），得 E?U??UX?K???1F (7.2.4)
式中：符号E[·]表示求均值，任一结点位移U的方差可由下式计算： Var?U???U??i?1j?1?Xinn?X?X?U?XjCov(Xi,Xj) (7.2.5)
X?X式中：符号Var[·]表示求方差；Cov(Xi,Xj)为Xi和Xj的协方差。 其中
?U?F?K?K?1(?U) (7.2.6) ?Xi?Xi?Xi???D?U?BU?DB (7.2.7) ?Xi?Xi?Xi
同样将σ在均值点处Taylor展开，也有与上面类似的表达式。可见，TSFEM关键在于
对有限元方程式直接进行偏微分计算，计算出有限元输出量对随机变量的梯度，故该法也称直接偏微分法或梯度分析法。
由于一阶TSFEM只需一次形成刚度矩阵，也只需一次求刚度矩阵的逆，因此效率较高。
2
但由于忽略了二阶以上的高次项，使TSFEM对随机变量的变异性有所限制。一般要求一阶TSFEM随机变量的变异系数小于0.3。如果随机变量的变异系数较大，可以采用有限元控制方程的二阶Taylor展开：
?2??2U?K?U?K?U?2K?1??F ?K???U? (7.2.8) ???Xi?Xj??Xi?Xj?Xi?Xj?Xi?Xj?Xi?Xj??2??2D?D?U?D?U?2U (7.2.9) ?BU?B??B?DB?Xi?Xj?Xi?Xj?Xi?Xj?Xj?Xi?Xi?Xj上式可见，二阶TSFEM可以放宽随机变量变异性大小的限制，但随机变量数目较多时，计算
量将十分庞大，而且一阶或二阶TSFEM均无法计算响应量三阶以上的统计特性。
由于TSFEM简单明了、效率高，为我国许多学者所采用。 ● 摄动法随机有限元
摄动技术最初被用于非线性力学分析。Handa等人成功地将一阶、二阶摄动技术用于随机问题，给出摄动法有限元列式。该法假定基本随机变量在均值点处产生微小摄动，利用Taylor级数把随机变量表示为确定部分和由摄动引起的随机部分，从而将有限元控制方程（非线性的）转化为一组线性的递推方程，求解得出位移的统计特性，进而求出应力的统计特性。
假设?i为随机变量Xi在均值点Xi处的微小摄动量，即?i?Xi?Xi。于是
?K1n?2K K?K0???i ???i?j??i2i,j?1??i??ji?1n(7.2.10)
对于U、F，也有类似上式K的表达式，式中：K0、U0、F0分别为K、U、F在随机变量均值点的值。根据二阶摄动法，可得
?1 U0?K0F0 (7.2.11)
??F?U?K?K0?1??U0?????i??i?i??? (7.2.12) ??? (7.2.13) ???2?2U?2K?U?K?K?U?1??F?K0?U0??
???????i??j??????????i??jijij?ij由上式可得位移的均值和协方差：
1nn??Cov(?i,?j) (7.2.14) E?U??U???Uij2i?1j?1?E(?i?j?k)Cov?U????Ui?U?jCov(?i,?j)????Ui?U?jk
i?1j?1ni?1j?1k?1nnnnn??Ukl??E(?i?j)E(?k?l)?E(?i?k)E(?j?l)?????Uiji?1j?1k?1l?1nnn?? (7.2.15)
由于任何量的随机性都可以引入摄动量，而且更易于考虑非线性问题，因此PSEFEM适
3
用范围较广，对于结构几何特性的随机性（包括随机边界问题）易得出随机有限元控制方程。一阶PSFEM和一阶TSFEM一样，只需一次形成刚度矩阵、一次对刚度矩阵求逆，计算效率较高。但PSFEM需以微小的摄动量为条件，一般应小于均值的20%或30%。 ● Neumann展开Monte-Carlo随机有限元
20世纪80年代后期，Shinozuka等人提出基于Neumann展开式的随机有限元法，使Monte-Carlo法与有限元法较完美地结合起来。Monte-Carlo法是最直观、最精确、获取信息最多、对非线性问题最有效的计算统计方法。Neumann展开式的引入是为了解决矩阵求逆的效率问题。如果对每一次随机抽样，只需形成刚度矩阵，进行前代、回代以及矩阵乘和矩阵加减，而无需矩阵分解，则可大大减少工作量。
在一般有限元控制方程KU = F中，假定荷载F为确定值，在随机变量波动值的影响下刚度矩阵K分解为K = K0+ΔK，根据Neumann级数展开，有
-1-123-1
K=(K0+ΔK)=(I-P+P-P+?)K0 (7.2.16) 式中：K0为随机变量均值处的刚度矩阵；ΔK为刚度矩阵的波动量；I为单位矩阵。对于Monte-Carlo随机抽样，刚度矩阵只改变ΔK项，而
-1
P = K0ΔK (7.2.17)
-1
U0 = K0F (7.2.18)
将式（7.2.16）代入式（7.2.1），并利用式（7.2.17）和式（7.2.18），得
23
U = U0-PU0+PU0-PU0+? (7.2.19) 令Ui=PiU0，则得如下的递推公式：
U = U0+U1+U2+? (7.2.20)
-1
Ui=-K0ΔKUi-1 (i=1,2,3,?) (7.2.21) 由式（7.2.18）求出U0后，可由式（7.2.21）求出Ui(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)的标准相关函数 ?(?)?R(?) (7.3.1) R(0)式中，随机场的协方差函数 R(τ)= Cov[S(t+τ),S(t)] (7.3.2)
2
对于一切t，随机场的方差为 Var[S(t)]= R(0)= σ (7.3.3) 相关函数也可以用谱分解表示（即Wiener-Khintchine变换对）：
?R(?)cos??d????? ? (7.3.4)
??R(?)??G(?)cos??d?????G(?)???12? 常用的标准相关函数有：非协调阶跃型、协调阶跃型、三角型、指数型、二阶AR型、
4
高斯型等多种形式。
2
3. 方差折减系数Γ(h)
设Sh(t)是随机场S(t)的局部平均随机过程，即
Sh(t)?1t?hS(t)dt (7.3.5) ?th2
2
则方差 Var[Sh(t)]= σΓ(h) (7.3.6)
2h?式中，方差折减系数 ?(h)??(1?)?(?)d? (7.3.7)
h0h2可见，方差折减系数起着使Sh(t)的方差比原来S(t)的方差缩小的作用。
4. 相关距离δS
δS可以看成是任意两个相隔距离为τ的随机变量不相关的最小距离（也称相关偏度）。
???s，则不相关，否则完全相关。利用相关距离δS便于对随机场作近似处理，其计算式
为：
?s?h?2(h),h????? (7.3.8) ??2??(?)d?0?2G(?)/?,??0?? 以上公式均对一维问题列出，二维、三维问题也可以类似得出。
● 随机场的中心离散
中心离散是随机场最简单的一种离散方法。该法用随机场在每个单元中心点的值来表征该随机场在每个单元的属性，因而随机场在每个单元内部都是常量，且等于它在各个单元中心的值。该法程序简单，但精度欠佳，现较少采用。 ● 随机场的局部平均
该法将随机场在每个单元的属性用随机场在单元上的局部平均来表征。Baecher、
[23][24]
Vanmarcke首先提出随机场局部平均的离散方法，我国学者曾对该法进行了深入研究。
1. 一维随机场的局部平均
2
设一个一维连续平稳随机场S(x)，其均值为m，方差为σ，则随机场在一个离散单元[t-T/2,t+T/2]上的局部平均定义为
ST(t)??T1t?(T/2)t?(T/2)S(x)dx (7.3.9)
式中：T是局部平均单元的长度，ST(t)称为局部平均随机场，其均值也为m，方差按（7.3.6）
22
式计算，h=T ，S(t)的标准相关函数为ρ(τ)=R(τ)/σ，无量纲方差折减系数Γ（h）有
2222
如下性质：①Γ（T）≥0；②Γ（0）=1；③Γ（-T）=Γ（T）。
考虑随机场在两个长度分别为T和T’的单元上的局部平均。如果局部平均随机场分别为ST和ST’，则其协方差为
Cov(ST,ST?)??(?1)2TT?k?0?23kTk2?2(Tk) (7.3.10)
2. 二维随机场的局部平均
设S(x1，x2)为二维连续参数连续状态随机场，Ai=L1iL2i为一矩形单元中心点（x1i，x2i）为中心，边平行于坐标轴x1与x2，且长为L1i与L2i的矩形面积，随机场在Ai内的局部平均
5