第7章随机有限元法.docx

第7章随机有限元法.docx

《第7章随机有限元法.docx》由会员分享，可在线阅读，更多相关《第7章随机有限元法.docx（26页珍藏版）》请在冰点文库上搜索。

第7章随机有限元法

第7章随机有限元法

§7.1绪论

结构工程中存在诸多的不确定性因素，从结构材料性能参数到所承受的主要荷载，如车流、阵风或地震波，无不存在随机性。

在有限单元法已成为分析复杂结构的强有力的工具和广泛使用的数值方法的今天，人们已不满足精度越来越高的确定性有限元计算，而设法用这一强有力的工具去研究工程实践中存在的大量不确定问题。

随机有限元法（StochasticFEM）,也称概率有限元法（ProbabilisticFEM）正是随机分析理论及有限元方法相结合的产物，是在传统的有限元方法的基础上发展起来的随机的数值分析方法。

最初是Monte-Carlo法及有限元法直接结合，形成独特的统计有限元方法。

Astill和Shinozuka（1972）首先将Monte-Carlo法引入结构的随机有限元法分析。

该法通过在计算机上产生的样本函数来模拟系统的随机输入量的概率特征，并对于每个给定的样本点，对系统进行确定性的有限元分析，从而得到系统的随机响应的概率特征。

由于是直接建立在大量确定性有限元计算的基础上，计算量极大，不适用于大型结构，而且最初的直接Monte-Carlo法还不是真正意义上的随机有限元法。

但及随后的摄动随机有限元法（PSFEM）相比，当样本容量足够大时，Monte-Carlo有限元法的结果更可靠也更精确。

结构系统的随机分析一般可分为两大类：

一类是统计方法，另一类是非统计方法。

因此，随机有限元法同样也有统计逼近和非统计逼近两种类型。

前者通过样本试验收集原始的数据资料，运用概率和统计理论进行分析和整理，然后作出科学推断。

这里，样本试验和数据处理的工作量很大，随着计算机的普及和发展，数值模拟法，如蒙特卡罗（MonteCarlo）模拟，已成为最常用的统计逼近法。

后者从本质上来说是利用分析工具找出结构系统的（确定的或随机的）输出随机信号及输入随机信号之间的关系，采用随机分析及求解系统控制方程相结合的方法得到输出信号的各阶随机统计量的数字特征（如各阶原点矩或中心矩）。

在20世纪70年代初，Cambou首先采用一次二阶矩方法研究线弹性问题。

由于这种方法将随机变量的影响量进行Taylor级数展开，就称之为Taylor展开法随机有限元（TSFEM）。

Shinozuka和Astill（1972）分别独立运用摄动技术研究了随机系统的特征值问题。

随后，Handa（1975）等人在考虑随机变量波动性时采用一阶和二阶摄动技术，并将这种摄动法随机有限元成功地应用于框架结构分析。

Vanmarcke等人（1983）提出随机场的局部平均理论，并将它引入随机有限元。

局部平均理论是用随机场函数在每一个离散单元上的局部平均的随机变量来代表该单元的统计量的近似理论。

LiuW.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为整体刚度矩阵

（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关键在于对有限元方程式直接进行偏微分计算，计算出有限元输出量对随机变量的梯度，故该法也称直接偏微分法或梯度分析法。

由于一阶TSFEM只需一次形成刚度矩阵，也只需一次求刚度矩阵的逆，因此效率较高。

但由于忽略了二阶以上的高次项，使TSFEM对随机变量的变异性有所限制。

一般要求一阶TSFEM随机变量的变异系数小于0.3。

如果随机变量的变异系数较大，可以采用有限元控制方程的二阶Taylor展开：

（7.2.8）

（7.2.9）

上式可见，二阶TSFEM可以放宽随机变量变异性大小的限制，但随机变量数目较多时，计算量将十分庞大，而且一阶或二阶TSFEM均无法计算响应量三阶以上的统计特性。

由于TSFEM简单明了、效率高，为我国许多学者所采用。

●摄动法随机有限元

摄动技术最初被用于非线性力学分析。

Handa等人成功地将一阶、二阶摄动技术用于随机问题，给出摄动法有限元列式。

该法假定基本随机变量在均值点处产生微小摄动，利用Taylor级数把随机变量表示为确定部分和由摄动引起的随机部分，从而将有限元控制方程（非线性的）转化为一组线性的递推方程，求解得出位移的统计特性，进而求出应力的统计特性。

假设

为随机变量

在均值点

处的微小摄动量，即

。

于是

（7.2.10）

对于U、F，也有类似上式K的表达式，式中：

K0、U0、F0分别为K、U、F在随机变量均值点的值。

根据二阶摄动法，可得

（7.2.11）

（7.2.12）

（7.2.13）

由上式可得位移的均值和协方差：

（7.2.14）

（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=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-1ΔK（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-1Δ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）的标准相关函数（7.3.1）

式中，随机场的协方差函数R（τ）=Cov[S（t+τ）,S（t）]（7.3.2）

对于一切t，随机场的方差为Var[S（t）]=R（0）=σ2（7.3.3）

相关函数也可以用谱分解表示（即Wiener-Khintchine变换对）：

（7.3.4）

常用的标准相关函数有：

非协调阶跃型、协调阶跃型、三角型、指数型、二阶AR型、高斯型等多种形式。

3.方差折减系数Γ2（h）

设Sh（t）是随机场S（t）的局部平均随机过程，即

（7.3.5）

则方差Var[Sh（t）]=σ2Γ2（h）（7.3.6）

式中，方差折减系数

（7.3.7）

可见，方差折减系数起着使Sh（t）的方差比原来S（t）的方差缩小的作用。

4.相关距离δS

δS可以看成是任意两个相隔距离为τ的随机变量不相关的最小距离（也称相关偏度）。

，则不相关，否则完全相关。

利用相关距离δS便于对随机场作近似处理，其计算式为：

（7.3.8）

以上公式均对一维问题列出，二维、三维问题也可以类似得出。

●随机场的中心离散

中心离散是随机场最简单的一种离散方法。

该法用随机场在每个单元中心点的值来表征该随机场在每个单元的属性，因而随机场在每个单元内部都是常量，且等于它在各个单元中心的值。

该法程序简单，但精度欠佳，现较少采用。

●随机场的局部平均

该法将随机场在每个单元的属性用随机场在单元上的局部平均来表征。

Baecher、Vanmarcke首先提出随机场局部平均的离散方法，我国学者曾对该法进行了深入研究[23][24]。

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’，则其协方差为

（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的互协方差可表示为

（7.3.12）

式中：

的约定如图7.3.1所示。

图7.3.1二维局部平均单元

如果有限元的网格已划分，且单元总数为n，随机场实际上被离散成n个随机变量，这n个随机变量的统计特性可由E[Si]、Var[Si]及互协方差Cov[Si，Sj]反映。

由于有限元网格的疏密是由应力梯度决定的，及随机场无关，当单元数很多时，随机场可另划网格，网格的疏密可由相关距离δS决定。

当然，随机场网格越密精度越高。

随机场的局部平均法由于对原始数据的要求低、收敛快、精度高，是随机有限元计算中最常用的方法。

●随机场的插值

LiuW.K.提出随机场的插值法。

该法将随机场在单元内的值用单元结点处值的插值函数来表示，于是随机场的统计特性可由各单元结点处随机变量间的统计特性近似反映。

利用形函数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提出随机场加权积分方法。

该法在单元刚度矩阵的推导过程中采用随机场在单元高斯点上的加权积分，以表征单元上的随机场。

假设单元e的弹性模量为

（7.3.16）

式中，E0（e）（X）为弹性模量的均值；f（e）（X）为一维零均值均匀随机场，其值域为

-1+η≤f（e）（X）≤-1-η（0＜η＜1）

两端铰接杆单元的刚度矩阵可近似表示为

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）实质上是在考虑弹性模量随机场的情况下，将单元刚度矩阵分解成确定性部分和随机部分。

式（7.3.18）等表征随机场在各单元的均值和方差，而单元间弹性模量的相关性则由下式表示

（7.3.20）

不难看出，局部平均法是加权积分法的特例（即权系数全部相同）。

由于采用加权积分，其计算精度相对较高，而且该法积分只需一次进行，刚度矩阵的波动性也由此得出，因此计算效率也较高。

●随机场的正交展开

Spanos和Ghanem提出的随机场正交展开法，将材料特性参数随机场进行Karhumen-Loeve正交展开，并由此推导出刚度矩阵的级数展开式，从而获得位移、应力的统计特性。

设随机场为

（7.3.21）

其中

式中：

λn、φn（x）分别为随机场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随机有限元及结构可靠度

●结构可靠性分析的一次二阶矩法

结构的安全性、适用性和耐久性统称为结构的可靠性。

可靠性的数学量度为可靠度，其定义为：

在规定的条件下和规定的时间内完成预定功能的概率。

“规定条件”是指正常设计、正常施工、正常使用的条件；“规定时间”是指结构的设计基准期。

结构的使用时间超过基准期后，其失效的概率将增大。

结构的一系列基本变量都具有不确定性，因此结构的可靠性分析属于随机性分析的范畴，随机有限元法将是十分有效的工具。

在结构可靠性分析中，结构的极限状态（包括承载能力极限状态、正常使用极限状态和条件极限状态）是通过功能函数来描述的。

当有n个随机变量影响结构可靠度时，结构的功能函数为:

z=g（X1,X2,…,Xn），基本变量Xi（i=1,2,…,n）是结构上的各种外因作用、材料性能和几何参数等。

z＞0时，结构处于可靠状态；z＜0,则处于失效状态；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失效边界

在工程实践中R和S不一定是正态分布，但可以变换成标准正态分布。

他们的统计特征量是均值μ、标准差σ、相关偏度或变异系数等，变异系数ν=σ/μ，表示随机变量相对于均值的变异。

目前工程上一般用无量纲的可靠指标β来反映结构的可靠度，

，β越大，失效概率Pf越小，其互补的可靠度Pr就越大。

一次二阶矩法采用只需已知均值和标准差的数学模型去求解结构的可靠度。

此法将功能函数z=g（X1,X2,…,Xn）在某点用Taylor级数展开，并近似地取一次项使极限状态方程线性化，然后求得可靠指标β。

改进的一次二阶矩法将线性化点选在失效边界上，而且选在结构最大可能失效点P*上（如图7.4.2）。

选择设计验算点P*（Xi*│i=1,2,…,n）作为线性化点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，其中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*,然后按距离公式确定结构的可靠指标。

具体的迭代格式如下：

（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）用随机有限元法计算S（i）和

；

（5）重复（3）、（4）两步骤，直至收敛（G（Y（i））=0）。

（6）再按式（7.4.5）和（7.4.6）计算可靠指标β值。

●随机有限元的最大熵法

一次二阶矩法的计算量随迭代次数成倍增加，使该法的使用受到限制。

而最大熵法用于结构的可靠性分析时，可根据随机变量的二阶矩来拟合概率密度曲线。

因此随机有限元结合最大熵法可用于求结构响应的概率密度曲线，从而计算结构的失效概率。

熵被定义为信息的均值，信息是对个别X值不确定性的度量。

不确定性越大，熵也越大。

对于一个连续随机变量，熵为

，而对于离散随机变量，熵为

，这里，f（X）是随机变量X的概率密度函数；f（xi）是离散点概率函数。

最大熵法通过调整概率密度函数f（X）使熵S取得最大值，并满足约束条件：

，利用最大熵法可得近似的概率密度函数f（X）。

为了求得最大值，可引入拉格朗日乘子法，把问题转