随机有限元法

1、 第7章 随机有限元法§7.1 绪论结构工程中存在诸多的不确定性因素，从结构材料性能参数到所承受的主要荷载，如车流、阵风或地震波，无不存在随机性。在有限单元法已成为分析复杂结构的强有力的工具和广泛使用的数值方法的今天，人们已不满足精度越来越高的确定性有限元计算，而设法用这一强有力的工具去研究工程实践中存在的大量不确定问题。随机有限元法（Stochastic FEM）,也称概率有限元法（Probabilistic FEM）正是随机分析理论与有限元方法相结合的产物，是在传统的有限元方法的基础上发展起来的随机的数值分析方法。最初是Monte-Carlo法与有限元法直接结合，形成独特的统计有

2、限元方法。Astill和Shinozuka（1972）首先将Monte-Carlo法引入结构的随机有限元法分析。该法通过在计算机上产生的样本函数来模拟系统的随机输入量的概率特征，并对于每个给定的样本点，对系统进行确定性的有限元分析，从而得到系统的随机响应的概率特征。由于是直接建立在大量确定性有限元计算的基础上，计算量极大，不适用于大型结构，而且最初的直接Monte-Carlo法还不是真正意义上的随机有限元法。但与随后的摄动随机有限元法（PSFEM）相比，当样本容量足够大时，Monte-Carlo有限元法的结果更可靠也更精确。结构系统的随机分析一般可分为两大类：一类是统计方法，另一类是非统计方法

3、。因此，随机有限元法同样也有统计逼近和非统计逼近两种类型。前者通过样本试验收集原始的数据资料，运用概率和统计理论进行分析和整理，然后作出科学推断。这里，样本试验和数据处理的工作量很大，随着计算机的普及和发展，数值模拟法，如蒙特卡罗（Monte Carlo）模拟，已成为最常用的统计逼近法。后者从本质上来说是利用分析工具找出结构系统的（确定的或随机的）输出随机信号与输入随机信号之间的关系，采用随机分析与求解系统控制方程相结合的方法得到输出信号的各阶随机统计量的数字特征（如各阶原点矩或中心矩）。在20世纪70年代初， Cambou首先采用一次二阶矩方法研究线弹性问题。由于这种方法将随机变量的影响量进

4、行Taylor级数展开，就称之为Taylor展开法随机有限元（TSFEM）。Shinozuka和Astill（1972)分别独立运用摄动技术研究了随机系统的特征值问题。随后，Handa（1975）等人在考虑随机变量波动性时采用一阶和二阶摄动技术，并将这种摄动法随机有限元成功地应用于框架结构分析。Vanmarcke等人（1983）提出随机场的局部平均理论，并将它引入随机有限元。局部平均理论是用随机场函数在每一个离散单元上的局部平均的随机变量来代表该单元的统计量的近似理论。Liu W. K.等人（1986、1988）的系列工作，提供了一种“主模态”技术，运用随机变量的特征正交化方法，将满秩的协方差

5、矩阵变换为对角矩阵，减少计算工作量，对摄动随机有限元法的发展做出贡献，此外，提出了一个随机变分原理。Yamazaki和Shinozuka(1987)创造性地将算子的Neumann级数展开式引入随机有限元的列式工作。从本质上讲，Neumann级数展开方法也是一类正则的小参数摄动方法，正定的随机刚度矩阵和微小的随机扰动量是两个基本要求，这两个基本要求保证了摄动解的正则性和收敛性，其优点在于摄动形式较简单并可以得到近似解的高阶统计量。Shinozuka等人（1987）将随机场函数的Monte-Carlo模拟与随机刚度矩阵的Neumann级数展开式结合,得到具有较好计算精度和效率的一类Neumann随

6、机有限元列式（称NSFEM）。Benaroya等（1988）指出，将出现以随机变分原理为基础的随机有限元法来逐渐取代以摄动法为基础的随机有限元法。Spanos和Ghanem等人（1989，1991）结合随机场函数的Karhuen-Loeve展式和Galerkin（迦辽金）射影方法建立了相应的随机有限元列式，并撰写了随机有限元法领域的第一本专著随机有限元谱方法。国内对随机有限元的研究起步较晚。吴世伟等人（1988）提出随机有限元的直接偏微分法及相应的可靠度计算方法。陈虬、刘先斌等人（1989、1991）提出一种新的随机场离散模型，建立了等参局部平均单元，并基于变分原理研究了一类随机有限元法的收敛

7、性和误差界。Papadrakakis(1995)采用预处理共轭梯度法给出了空间框架的非线性随机有限元列式。Schorling和Bucher(1996)基于Monte-Carlo技术，采用响应面法研究几何非线性时的可靠度随机有限元方法。刘宁（1996）则基于偏微分法，给出了三维弹塑性随机有限元列式。随机有限元法的数学理论研究和非线性随机问题的有限元分析工作还有待深入。自20世纪80年代以来，随机有限元法已在工程结构可靠性、安全性分析领域以及在各种随机激励下结构响应变异研究领域中得到应用，如应用于大型水利工程的重力坝、拱坝的可靠度计算；应用于非线性瞬态响应分析；结构振动中随机阻尼对响应的影响；结构

8、分析的随机识别；复杂结构地震响应的随机分析和两相动力系统的随机模拟等等。随着理论研究的深入，随机有限元将得到更加广泛的应用。§7.2 随机有限元的控制方程22从随机有限元控制方程的获得来看，随机有限元可分为Taylor展开法随机有限元（TSFEM）、摄动法随机有限元(PSFEM)以及Neumann展开Monte-Carlo法随机有限元（NSFEM）。 Taylor展开法随机有限元 该随机有限元法的基本思路是将有限元格式中的控制量在随机变量均值点处进行Taylor级数展开（取一阶或二阶），经过适当的数学处理得出所需的计算方程式。有限元静力分析控制方程的矩阵形式为: KU = F (7.

9、2.1)式中,U为位移矩阵，F为等效节点荷载列阵，K为整体刚度矩阵 (7.2.2)其中，B为形变矩阵,D为材料弹性矩阵。在计算出节点位移U后，即由下式求得应力列阵= DBU (7.2.3)设基本随机变量为，将位移U在均值点处一阶Taylor级数展开，并在两边同时取均值（数学期望），得 (7.2.4)式中：符号E·表示求均值，任一结点位移U的方差可由下式计算： (7.2.5)式中：符号Var·表示求方差；Cov(Xi,Xj)为Xi和Xj的协方差。其中 (7.2.6) (7.2.7)同样将在均值点处Taylor展开，也有与上面类似的表达式。可见，TSFEM关键在于对有限元方程式

10、直接进行偏微分计算，计算出有限元输出量对随机变量的梯度，故该法也称直接偏微分法或梯度分析法。由于一阶TSFEM只需一次形成刚度矩阵，也只需一次求刚度矩阵的逆，因此效率较高。但由于忽略了二阶以上的高次项，使TSFEM对随机变量的变异性有所限制。一般要求一阶TSFEM随机变量的变异系数小于0.3。如果随机变量的变异系数较大，可以采用有限元控制方程的二阶Taylor展开： (7.2.8) (7.2.9)上式可见，二阶TSFEM可以放宽随机变量变异性大小的限制，但随机变量数目较多时，计算量将十分庞大，而且一阶或二阶TSFEM均无法计算响应量三阶以上的统计特性。由于TSFEM简单明了、效率高，为我国许多

11、学者所采用。 摄动法随机有限元摄动技术最初被用于非线性力学分析。Handa等人成功地将一阶、二阶摄动技术用于随机问题，给出摄动法有限元列式。该法假定基本随机变量在均值点处产生微小摄动，利用Taylor级数把随机变量表示为确定部分和由摄动引起的随机部分，从而将有限元控制方程（非线性的）转化为一组线性的递推方程，求解得出位移的统计特性，进而求出应力的统计特性。假设为随机变量在均值点处的微小摄动量，即。于是 (7.2.10)对于U、F，也有类似上式K的表达式，式中：K0、U0、F0分别为K、U、F在随机变量均值点的值。根据二阶摄动法，可得 (7.2.11) (7.2.12) (7.2.13)由上式可

12、得位移的均值和协方差： (7.2.14) (7.2.15)由于任何量的随机性都可以引入摄动量，而且更易于考虑非线性问题，因此PSEFEM适用范围较广，对于结构几何特性的随机性（包括随机边界问题）易得出随机有限元控制方程。一阶PSFEM和一阶TSFEM一样，只需一次形成刚度矩阵、一次对刚度矩阵求逆，计算效率较高。但PSFEM需以微小的摄动量为条件，一般应小于均值的20%或30%。 Neumann展开Monte-Carlo随机有限元20世纪80年代后期，Shinozuka等人提出基于Neumann展开式的随机有限元法，使Monte-Carlo法与有限元法较完美地结合起来。Monte-Carlo法是

13、最直观、最精确、获取信息最多、对非线性问题最有效的计算统计方法。Neumann展开式的引入是为了解决矩阵求逆的效率问题。如果对每一次随机抽样，只需形成刚度矩阵，进行前代、回代以及矩阵乘和矩阵加减，而无需矩阵分解，则可大大减少工作量。在一般有限元控制方程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为单位矩阵。对于Monte-Carlo随机抽样，刚度矩阵只改变K项，而P = K0-

14、1K (7.2.17)U0 = K0-1F (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，则得如下的递推公式： U = U0+U1+U2+ (7.2.20) Ui=-K0-1KUi-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随机

15、模拟技术，NSFEM不受随机变量波动范围的限制，当变异系数小于0.2时，NSFEM与一阶TSFEM或一阶PSFEM精度相当；当变异系数大于0.2时，后两者已不能满足精度要求，但NSFEM仍能得出满意的结果。§7.3 随机场的离散模型许多物理现象和物体系统具有随机分布特性，包括系统本身的不确定或系统的激励和响应的不确定，都可以模型化为随机空间分布的随机场或随时间分布的随机过程。随机有限元法除了必须考虑材料参数等的空间变异性，需要获得随机有限元方程列式以及解决随机算子和随机矩阵的求逆问题外，还须包含对随机场的离散处理。均匀各向同性随机场的特征量1. 随机场S(t)的均值E(S(t)为常数

16、m.2. 随机场S(t)的标准相关函数 (7.3.1)式中，随机场的协方差函数 R()= CovS(t+),S(t) (7.3.2)对于一切t，随机场的方差为 VarS(t)= R(0)= 2 (7.3.3)相关函数也可以用谱分解表示（即Wiener-Khintchine变换对）： (7.3.4) 常用的标准相关函数有：非协调阶跃型、协调阶跃型、三角型、指数型、二阶AR型、高斯型等多种形式。3. 方差折减系数2(h)设Sh(t)是随机场S(t)的局部平均随机过程，即 (7.3.5)则方差 VarSh(t)= 22(h) (7.3.6)式中，方差折减系数 (7.3.7)可见，方差折减系数起着使S

17、h(t)的方差比原来S(t)的方差缩小的作用。4. 相关距离SS可以看成是任意两个相隔距离为的随机变量不相关的最小距离（也称相关偏度）。，则不相关，否则完全相关。利用相关距离S便于对随机场作近似处理，其计算式为： (7.3.8) 以上公式均对一维问题列出，二维、三维问题也可以类似得出。 随机场的中心离散中心离散是随机场最简单的一种离散方法。该法用随机场在每个单元中心点的值来表征该随机场在每个单元的属性，因而随机场在每个单元内部都是常量，且等于它在各个单元中心的值。该法程序简单，但精度欠佳，现较少采用。 随机场的局部平均 该法将随机场在每个单元的属性用随机场在单元上的局部平均来表征。Baeche

18、r、Vanmarcke首先提出随机场局部平均的离散方法，我国学者曾对该法进行了深入研究2324。1. 一维随机场的局部平均设一个一维连续平稳随机场S(x)，其均值为m，方差为2，则随机场在一个离散单元t-T/2,t+T/2上的局部平均定义为 (7.3.9)式中：T是局部平均单元的长度，ST(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的单元上的局部平均。如果局部平均随机场分别为ST和ST，则其协方差为 (

19、7.3.10)2. 二维随机场的局部平均设S(x1，x2)为二维连续参数连续状态随机场，Ai=L1iL2i为一矩形单元中心点（x1i，x2i）为中心，边平行于坐标轴x1与x2，且长为L1i与L2i的矩形面积，随机场在Ai内的局部平均定义为 (7.3.11)如果S(x1,x2)是一个均匀随机场，则可用均值m、方差2及归一化协方差函数（r1,r2）近似描述，r1和r2分别为沿两个方向的距离。对应的局部平均随机场可用E(Si)、Var(Si)及互协方差Cov(Si,Sj)近似描述。E(Yi)、Var(Yi)及Cov(Yi,Yj)可由m、2及(r1，r2)计算获得。第i和第j单元局部平均Si、Sj的互

20、协方差可表示为 (7.3.12)式中：的约定如图7.3.1所示。 二维局部平均单元 如果有限元的网格已划分，且单元总数为n，随机场实际上被离散成n个随机变量，这n个随机变量的统计特性可由ESi、VarSi及互协方差CovSi，Sj反映。由于有限元网格的疏密是由应力梯度决定的，与随机场无关，当单元数很多时，随机场可另划网格，网格的疏密可由相关距离S决定。当然，随机场网格越密精度越高。随机场的局部平均法由于对原始数据的要求低、收敛快、精度高，是随机有限元计算中最常用的方法。 随机场的插值Liu W. K.提出随机场的插值法。该法将随机场在单元内的值用单元结点处值的插值函数来表示，于是随机场的统计特

21、性可由各单元结点处随机变量间的统计特性近似反映。利用形函数Ni(X)，随机场b(X)离散式表示为 (7.3.13)式中：X表示空间位置；bi为随机场在结点i处的值（i=1，2，q）；q为单元结点数。随机场b(X)在单元内的均值和方差可表示为 (7.3.14) (7.3.15) 随机场的插值法将原连续状态的随机场仍离散成一个连续函数，未直接计算随机场引起的单元间的相关性，只需给定随机场在各结点上的值，计算相对简单，易于考虑非线性和非均匀随机场问题。但需要已知相关函数，并且要求随机场对空间参数具有较高的连续性。 随机场的加权积分法Takada、Shinozuka及Deodatis提出随机场加权积分

22、方法。该法在单元刚度矩阵的推导过程中采用随机场在单元高斯点上的加权积分，以表征单元上的随机场。假设单元e的弹性模量为 (7.3.16)式中，E0(e)(X)为弹性模量的均值；f(e)(X)为一维零均值均匀随机场，其值域为 -1+f(e)(X)-1- (01)两端铰接杆单元的刚度矩阵可近似表示为 K(e)K0(e)+X(e)K0(e) (7.3.17)式中，K0(e)为弹性模量取均值时的单元刚度矩阵，而 （L为杆长） (7.3.18)固结杆单元的刚度矩阵可表示为 K(e)K0(e)+X0(e)K0(e)+X1(e)K1(e)+X2(e)K2(e) (7.3.19)其中： 式（7.3.17）实质上

23、是在考虑弹性模量随机场的情况下，将单元刚度矩阵分解成确定性部分和随机部分。式（7.3.18）等表征随机场在各单元的均值和方差，而单元间弹性模量的相关性则由下式表示 (7.3.20) 不难看出，局部平均法是加权积分法的特例（即权系数全部相同）。由于采用加权积分，其计算精度相对较高，而且该法积分只需一次进行，刚度矩阵的波动性也由此得出，因此计算效率也较高。 随机场的正交展开Spanos和Ghanem提出的随机场正交展开法，将材料特性参数随机场进行Karhumen-Loeve正交展开，并由此推导出刚度矩阵的级数展开式，从而获得位移、应力的统计特性。设随机场为 (7.3.21)其中 式中：n、n(x)

24、分别为随机场S(x)相关函数的特征值和特征函数。n(x)具有如下正交性： （Kronecker函数） (7.3.22)式（7.3.21）对于任何分布的S(x)均收敛，可取至第r阶以满足精度要求。对于一维杆单元，单元刚度矩阵可表示为 其中 式中：Pe(x)=De(x)/S(x), De(x)为单元弹性矩阵。集成整体刚度矩阵，得 (7.3.23)从上式，可得位移： (7.3.24)利用Neumann展开式，可进一步得出位移的统计特性。该法关键在于获得特征值和特征函数。§7.4 随机有限元与结构可靠度 结构可靠性分析的一次二阶矩法结构的安全性、适用性和耐久性统称为结构的可靠性。可靠性的数学

25、量度为可靠度，其定义为：在规定的条件下和规定的时间内完成预定功能的概率。“规定条件”是指正常设计、正常施工、正常使用的条件；“规定时间”是指结构的设计基准期。结构的使用时间超过基准期后，其失效的概率将增大。结构的一系列基本变量都具有不确定性，因此结构的可靠性分析属于随机性分析的范畴，随机有限元法将是十分有效的工具。在结构可靠性分析中，结构的极限状态（包括承载能力极限状态、正常使用极限状态和条件极限状态）是通过功能函数来描述的。当有n个随机变量影响结构可靠度时，结构的功能函数为: z=g(X1,X2,Xn)，基本变量Xi(i=1,2,n)是结构上的各种外因作用、材料性能和几何参数等。z0时，结构

26、处于可靠状态；z0,则处于失效状态；z=0称结构极限状态方程（一般难以用显式表示），结构处于极限状态。如果z的概率密度函数或概率分布函数可求得，则结构可靠度的数量指标便可基于各种状态出现的概率而确定。若功能函数仅与两个随机变量有关（如结构抵抗破坏或变形的能力R和荷载引起结构内力、应力、位移等效应S），即z=g(R,S)。假设R、S均为正态分布，其均值和标准差分别为和，此时z=R-S也是正态分布的随机变量，具有均值和标准差。其概率密度函数为 (7.4.1)其分布图示于图7.4.1，阴影部分是结构失效概率Pf，非阴影部分面积即结构的可靠度Pr。 图7.4.1 正态分布概率密度函数 图7.4.2 失

27、效边界在工程实践中R和S不一定是正态分布，但可以变换成标准正态分布。他们的统计特征量是均值、标准差、相关偏度或变异系数等，变异系数=/，表示随机变量相对于均值的变异。目前工程上一般用无量纲的可靠指标来反映结构的可靠度，越大，失效概率Pf越小，其互补的可靠度Pr就越大。一次二阶矩法采用只需已知均值和标准差的数学模型去求解结构的可靠度。此法将功能函数z=g(X1,X2,Xn)在某点用Taylor级数展开，并近似地取一次项使极限状态方程线性化，然后求得可靠指标。改进的一次二阶矩法将线性化点选在失效边界上，而且选在结构最大可能失效点P*上（如图7.4.2）。选择设计验算点P*(Xi*i=1,2,n)作

28、为线性化点X0i时，线性化的极限状态方程为 (7.4.2)则z的均值为 由于设计验算点P*选在失效边界上，有g(X1*,X2*,Xn*)=0，因此z成为 但随机变量Xi互不相关时，z的标准差z为 其中 称i为灵敏度系数，表示第i个随机变量对标准差的相对影响。于是可靠指标为 (7.4.3)变换上式为 即 (i=1，2，n) (7.4.4)上式中，Xi，Xi为已知的各随机变量的均值和标准差，待求的量为Xi*和，可迭代求解。 随机有限元的一次二阶矩法设可靠性分析中的一组基本变量X=(X1,X2,Xn)T为相互独立的正态变量（不满足时可作变换），并进一步变换为标准正态变量Y=(Y1,Y2,Yn)T，其

29、中Yi=(Xi-Xi)/Xi。于是功能函数也可转换到标准正态空间，即 g(X)=g(R(X),S(X)=g(R(T-1(Y),S(T-1(Y)=G(Y)失效概率Pf=1-(),由的几何含义可知，值为在标准正态空间中从原点到失效面的最短距离。在标准正态空间中概率密度是关于原点（均值点）旋转对称，并且随着到原点的距离的平方呈指数下降。在一次二阶矩法中，标准正态空间中的极限状态面（失效面）被一个到原点最小距离点处的切平面代替，也即按一阶Taylor展开式把功能函数在设计点处展开，使之线性化。采用迭代法可以近似确定极限状态面G(Y)=0上距原点最近的点Y*,然后按距离公式确定结构的可靠指标。具体的迭代

30、格式如下： (7.4.5) (7.4.6)以Mises屈服准则为例，功能函数 g(S)=s02-(s2xx-sxxsyy+s2yy+3s2xy)于是极限状态方程 g(S(X(Y)=G(Y)=0 (7.4.7)计算可靠指标时要用到功能函数G(Y)的梯度向量，由求导链式法则 (7.4.8)由于S(X)难于用显式表达，求存在困难。有效方法是利用随机有限元一次二阶矩法计算结构的可靠指标和设计点Y*。步骤如下：（1） 确定随机变量转换关系Y=Y(X)；（2） 给定初值Y(0)=0，计算式（7.4.8）右边的三个偏导数；（3） 按式（7.4.5）和（7.4.6）计算设计验算点坐标Y*；（4） 用随机有限元

31、法计算S(i)和；（5） 重复（3）、（4）两步骤，直至收敛（G(Y(i)=0）。（6） 再按式（7.4.5）和（7.4.6）计算可靠指标值。 随机有限元的最大熵法一次二阶矩法的计算量随迭代次数成倍增加，使该法的使用受到限制。而最大熵法用于结构的可靠性分析时，可根据随机变量的二阶矩来拟合概率密度曲线。因此随机有限元结合最大熵法可用于求结构响应的概率密度曲线，从而计算结构的失效概率。熵被定义为信息的均值，信息是对个别X值不确定性的度量。不确定性越大，熵也越大。对于一个连续随机变量，熵为 ，而对于离散随机变量，熵为，这里，f(X)是随机变量X的概率密度函数；f(xi)是离散点概率函数。最大熵法通过

32、调整概率密度函数f(X)使熵S取得最大值，并满足约束条件： ，利用最大熵法可得近似的概率密度函数f(X)。为了求得最大值，可引入拉格朗日乘子法，把问题转化为确定拉格朗日乘子。选用基于Neumann级数的随机有限元法，写出递推方程。可将局部平均理论引入Neumann随机有限元法，以进一步提高计算效率。随机有限元最大熵法计算结构失效概率的步骤如下：(1) 随机场离散，随机变量化为随机向量，并计算协方差矩阵；(2) 利用特征正交化技术，用一组统计独立的随机变量描述随机场；(3) 利用Neumann随机有限元法计算响应量的前若干阶矩；(4) 拟合所有节点或部分节点的最大熵密度函数；(5) 利用数值积分

33、计算结构的失效概率。 算例例1 如图7.4.3所示桁架下端受一集中力，假设各设计变量都是均匀分布的随机变量，三根杆弹性模量E，E=2.0×106N/cm2，Cov(E)=2.5%；截面积A，A=0.8cm2，Cov(A)=2%；杆2长度L2,L2=20cm，Cov(L2)=1.5%；杆1和杆3的长度均为L1，L1=28.284cm，Cov(L1)=2.3%；外力P，P=4000N；Cov(P)=10%；材料屈服应力y，y=3200N/cm2，Cov(y)=7%。当极限状态方程取为 y-S=0（S为构件中最大应力值）时，试求该结构的失效概率。图7.4.3 简单超静定桁架解： 用随机有限

34、元一次二阶矩法计算，假设由各结点坐标值变异引起的方向余弦的变异很小，形成确定矩阵，则由随机有限元递推方程组可得各结点位移分布特征值。计算得结点4的位移均值和标准差为4=2.2928cm和4=0.0037cm；各杆件中最大应力为杆元2的应力，其均值和标准差为2=2928N/cm2和2=377N/cm2。由一次二阶矩法可计算得该结构的失效概率为0.23（由于结构超静定，杆2失效并不意味着整个结构失效），与蒙特卡罗有限元法计算结果比较接近（蒙特卡罗取样数为40时，失效概率达0.224）。随机有限元法在分析混凝土重力坝可靠度方面已得到较多应用。研究者采用抽样方法对重力坝施行随机自动剖分，研究其可靠指标

35、的统计规律和收敛于真解的规律，并进行结构优化25。在对碾压混凝土重力坝抗滑稳定的体系可靠度研究中，还采用刚体极限平衡法、线弹性随机有限元法和三维弹塑性随机有限元法进行计算比较。§7.5 随机有限元动力分析方法 不确定结构的自振特性假设某一梁-柱结构承受随机分布的轴向静荷载N(x)，沿长度有微小扰动，同时在端部有轴向推力c为随机变量，其方差为，轴向分布荷载可表示为式中p(x)是一维、均值为零的均匀随机场，方差是P2，相关偏度是P，互相关函数是RPP()。结构弹性模量和质量密度也随机变化，分别表示为，其中分别是弹性模量E和质量密度m的均值，方差分别为2E和2m；随机变量a(x)和b(x)

36、是两个一维互不相关的、均匀的、均值为零的随机场。自相关函数分别为Raa()和Rbb()或等价的功率谱密度函数为Saa()和Sbb()，相关偏度分别是E和m。任一点横向位移的有限元模式为 w(x,t)=Neqe，其中，单元位移列阵qe=w1,1,w2,2T,形函数Ne取三次多项式插值。单元的动能表达为： (7.5.1)式中：，Ae为单元横截面面积。单元的应变能为 (7.5.2)式中：，I是惯性矩。单元的外力功为 (7.5.3)式中：，p为任一随机变量据此集成结构总能量和总外力功，利用Hamilton原理，即对总位移向量q取变分，考虑到q的任意性，可得系统运动方程，进而得到如下特征问题方程： ）式

37、中，刚度系数，其中是确定性分量，而随机扰动项：。于是可以写出刚度系数、质量系数以及几何刚度系数的均值和方差的表达式26。在总刚度矩阵中两个刚度系数kij和krs的互协方差，当单元尺寸相等时，利用方差函数可作进一步简化。利用局部平均理论，在每个单元上E、m、Q的均值为零，方差为： (7.5.5)式中，和分别为随机变量a(x)、b(x)和p(x)的方差函数和相关偏度。利用协方差函数写出两刚度系数或两质量系数的互协方差。最后，特征值问题的方程可表达为 (7.5.6)或 K*x=M*x 。相应平均问题的方程为 (7.5.7) 利用单元刚度、质量和几何刚度之间的协方差的表达式，等效刚度矩阵K*的协方差矩阵能用独立的a(x)、b(x)和p(x)来表示。解未扰动的特征值问题（平均值问题）可得特征值的均值。对于满足小扰动情况，可推导出求两个