随机有限元法

第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 为整体刚度矩阵∑⎰⎰⎰Ω=eT 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 XX 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 ii 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) ji 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 -=α。

于是 ∑∑==∂∂∂+∂∂+≈ni 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 ji 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 ijn i n j n i n j k j i jkn k i j i j i E E E E U U E U U Cov U U U Cov 111111111)()()()()(),(ααααααααααααα (7.2.15)由于任何量的随机性都可以引入摄动量，而且更易于考虑非线性问题，因此PSEFEM 适用范围较广，对于结构几何特性的随机性（包括随机边界问题）易得出随机有限元控制方程。

一阶PSFEM 和一阶TSFEM 一样，只需一次形成刚度矩阵、一次对刚度矩阵求逆，计算效率较高。

但PSFEM 需以微小的摄动量为条件，一般应小于均值的20%或30%。