结构可靠度分析的有限单元数值逼近法

刘旭

河南理工大学土木工程学院 焦作454000

结构可靠性分析［1，2］对工程结构的设计、稳定分析以及优化设计具有着重要意义。结构可靠性的度量称为可靠度。当今主流的可靠度理论认为结构体包含的各种参数如几何尺寸、抗力、荷载效应等，多为服从某种概率分布的随机变量，是其诸多不确定性因素的来源［3］。结构可靠度分析是将上述不确定性因素引入概率统计分析的理论［4］，描述结构体功能状态的函数表示为：

式中：X ＝（X1，X2，…，Xn）是n个影响结构体功能的基本随机变量。当Z＞0、Z＜0、Z＝0 时分别表示结构体处于可靠状态、失效状态和极限状态。设X的联合概率密度函数为fX（x），结构失效概率Pf的表达式为：

这是一个基于失效域的多维积分问题，通常由于极限状态面（Z＝0）在X 空间中难以确定，因此直接采用式（2）计算Pf是非常困难的。目前行而有效的方法大体分为两大类：近似法和模拟法。

近似法是使用平面或曲面替代原来的极限状态曲面（失效面），此类方法主要有：国际结构安全性联合委员会（JCSS）推荐的JC 法［5］，并以次为基础发展了二次二阶矩法、高次高阶矩法；
以及随着计算机技术的发展，近些年主流理论中的响应面法［6，7］、支持向量机法［8］、神经网络拟合法［9］等。也可利用函数近似理论对结构的功能函数和分布函数做近似处理的矩法，例如：孟广伟等人［10－12］提出利用降维算法，建立n个一维函数替代（1）式中的n维功能函数，并结合Edgeworth 级数拟合近似功能函数的分布函数从而直接计算失效概率。近似法的不足在于功能函数或失效面的数学特性对近似效果的影响较大，如JC法的计算精度很大程度上取决于失效面在设计验算点附近的非线性程度，且该类方法难以得到准确的误差范围。模拟法以蒙特卡洛法（MCS）［13］为主，原理为对已知概率分布的随机变量产生足够多的随机数，对功能函数进行重复模拟。虽然此类方法计算精度较高且无求解条件限制，但当处理小概率问题时，需要大量的模拟次数保证结果的精度［14］，显得有些捉襟见肘。由于是基于随机模拟的计算结果，所以该类方法同样难以对结果做出准确的误差分析（非概率误差）。

因此，本文提出一种基于有限单元法［15］处理X分布空间，计算结构失效概率的新方法。该算法分为四个部分：（1）利用功能函数中基本变量的概率集中性，在预设误差限的要求下，通过舍弃X的低概率区间保留高概率区间，建立有限单元法的计算域。（2）基于有限单元法对计算域进行均匀划分，得到若干个互不重叠、相互独立的单元，并规定以单元的几何中心点所处位置表征其是否位于失效域。（3）引入功能函数对每个单元的中心点进行筛选，保留位于失效域的单元并累加其概率得到近似的失效概率。（4）通过MATLAB编程不断提高划分单元的精细度，逐次迭代、缩小误差，实现数值逼近的计算过程，得到满足预设误差限的失效概率。

1.1 建立X空间的计算域

记Xi（Xi∈X，i＝1，2，…，n）的概率密度函数和分布函数分别为：fi（xi），Fi（xi）。由于Xi的分布具有概率集中性，为提高计算效率，可利用此性质定义Xi的高概率区间，即包含大部分概率值的区间。记Xi的高概率区间为ωi，ωi

内包含的概率为Pωi，如：

将式（5）转变为h与xi的关系，并对h取最小值得到：

因此在计算中根据要求给定Pωi后，即可利用式（6 ～8）求得对应式（3）中ωi的区间分布范围。例如：已知X1服从正态分布（μX1＝20，σX1＝4）；
X2服从对数正态分布（μX2＝22，σX2＝2）；
X3服从极值Ⅰ型分布（μX3＝14，σX3＝3.5）。在给定Pω1＝Pω2＝Pω3＝0.99 后，经上述方法分别得到X1、X2、X3的高概率区间为：ω1＝（9.6967，30.3033）、ω2＝（17.1853，27.4766）、ω3＝（7.2575，25.3613），三者的演示如图1 所示。

图1 X1、 X2、 X3 的高概率区间Fig.1 High probability interval of X1，X2，X3

综合考虑功能函数中n个基本变量的高概率区间，组成了X的高概率区间，在空间中表现为一种多维超方体，记作Ω，也即是有限单元法的计算域，形式为：

例如：将图1 中的三个高概率区间ω1、ω2、ω3共同组成Ω时，演示如图2 所示。

图2 ω1、 ω2、 ω3 组成的ΩFig.2 Ω consisting of ω1 and ω2 and ω3

Ω内包含的概率记作PΩ，PΩ是一个基于Ω内对联合概率密度函数的多重积分，由于随机变量的独立性，结合式（3）、式（4）、式（9）PΩ表示为：

在计算域Ω 内计算结构失效概率时，为达到目标精度，一般要求PΩ包含大部分概率值，即PΩ的取值近似于1。因此，为了计算方便和平衡每个ωi的压缩效率，对式（11）作如下规定：

据此结合式（11）、式（12）得到Pωi与PΩ的关系为：

因此当有n个随机变量时，可根据预设误差限的要求选取合适的PΩ，利用式（13）计算出Pωi，将其带入式（6 ～8）得到每个ωi的区间分布，最后根据式（9）建立计算域Ω，实现对X 分布空间的高效压缩。

1.2 计算域Ω内的有限单元法

压缩变量分布空间的过程实现了在区间范围较小的计算域Ω内，富集了满足计算要求的概率值。然后介绍均匀划分Ω的具体方法，采用有限单元法需要提出三点要求：（1）每个单元的相关计算需简便；
（2）每个单元内部需要有一点表征该单元是否位于失效域。（3）所有单元应毫无间隔且无重叠的排满计算域。

据此三点本文采用垂直坐标轴等距划分的方法，即将Ω中每个随机变量的高概率区间均分N份，得到Nn个尺寸相同的单元，外形上此单元可以视为等比缩小Nn倍的Ω。根据Ω空间内方向的不同分别以j、k、…、l代表单元在X1、X2、…、Xn坐标轴正方向上的位置排序（其中j、k、…、l的范围为1 到N之间的整数），将此方法得到的单元记为Aj，k，…，l。例如：对图2 中由ω1、ω2、ω3组成的Ω设置均分份数N＝10 时，A1，1，5表示X1坐标轴正方向上的第1 个、X2坐标轴正方向上的第1 个、X3坐标轴正方向上的第5 个单元，演示如图3 所示。

图3 均分方法的示意Fig.3 Schematic diagram of uniform segmentation

记Aj，k，…，l在Xi方向上的边长为hi，由上述均分方法可知其表达式为：记Aj，k，…，l的几何中心点坐标为Oj，k，…，l，结合式（16）得到其表达式为：

规定以Oj，k，…，l的坐标表征该单元是否位于失效域，因此几何尺寸hi和中心点坐标Oj，k，…，l为单元的关键参数。

1.3 由单元计算近似失效概率

记单元Aj，k，…，l内部包含的概率为Pj，k，…，l，类似地根据多重积分的性质和随机变量的独立性，使用单元的hi和Oj，k，…，l计算Pj，k，…，l的表达式为：

计算域Ω 被失效面分为可靠域与失效域，将Aj，k，…，l的中心点坐标Oj，k，…，l带入式（1）功能函数中，若G（Oj，k，…，l）＜0 则表示该单元位于失效域。因此，累加位于失效域单元的概率可以计算近似失效概率。为了实现有效的筛选，引入指示函数I（x）。如下：

结合式（16）、式（17）在Ω 内计算的近似失效概率P′f可以表示为：

综上所述，使用式（18）计算的P′f由计算域的概率PΩ和均分份数N的取值确定。

2.1 误差来源

一般来说，可靠度分析的误差源有两个［16］：一个是变量统计误差，即设计变量的模拟误差，相当于变量数学模型的可靠性问题；
另一个是数值误差，是由于建立算法模型所产生的误差。本文主要讨论算法的数值误差，记式（2）中失效概率的解析值

本文算法中使用式（18）计算P′f时，数值误差来源有两个。其一为，建立Ω 在高效压缩X分布空间的同时，忽略了Ω 以外的低概率区间，产生了数值误差。将Ω内计算的失效概率解析值记为其与的误差记作范围误差，范围误差限记为Pε，Pε是一种由PΩ 决定的误差，二者的关系为：

其二为，在Ω 内使用有限单元法产生的误差。位于失效面较近的若干单元会出现失效面穿过这些单元，使其一部分处于可靠域，另一部分处于失效域。由于使用单元的中心点来表征其是否失效是一种非此即彼的前提假设，因此产生了算法误差。图4 演示了二维变量的情况，按照中心点表征原则，该单元位于可靠域，从而忽略了位于失效域的部分，产生了相应误差。

图4 失效面穿过区域单元的情况Fig.4 The situation of failure face through the regional unit

2.2 误差分析

研究发现，第二种误差与均分份数N有关。一般N取值越大，与失效面相交的单元包含概率值越小，该误差也就越小。随着N增大，在Ω内计算的近似失效概率逐渐收敛因此，为得到一定精度下的失效概率，需要不断地增加均分份数N，代入式（18）求算新的近似失效概率，数值逼近为使每次增加均分份数N计算的近似失效概率稳健收敛，本文采用倍增均分份数的方法进行迭代。即当初次均分份数为N0时，第m（m＝1，2，…）次迭代中的均分份数为：

记第m－1 次迭代计算的近似失效概率为第m次将具体化为：

式中：α1，α2…等为0 到9 中的一个数字，且α1≠0。当迭代次数m足够大时，的数值充分接近，若其存在下式关系：

与一般近似法、模拟法等其他方法相比，本文算法的一个优点在于根据式（24）给出了准确的误差范围。

在上述计算近似失效概率的框架中，PΩ和N为主要参数：PΩ指导构建计算域Ω的空间范围，N为计算域Ω内基于有限单元法的计算过程提供划分精度，二者的取值分别对应算法误差分析中的Pε与Pε。因此该算法从预设误差限出发，得到各种参数，迭代计算满足精度要求的失效概率。

综上所述，迭代步骤如下：

（1）预设计算结果的范围误差限Pε，并给式（23）中的c赋值，设定迭代误差限Pε。

（2）根据设定的Pε利用式（19）计算PΩ。

（3）根据PΩ，利用式（13）、式（6 ～8）计算Ω的空间分布。

（4）设定初始均分份数N0。

（5）根据设定的Pε将式（22）确定为迭代终止条件。

（6）利用式（20）确定与迭代次数m对应的均分份数N。

（7）利用式（14 ～18）计算第m次迭代求算

（8）判断第m次计算的与第m－1 次计算的是否满足步骤（5）中的迭代终止条件，若不满足则令m＝m＋1 重复步骤（6 ～7），计算新的近似失效概率。

（9）若第m次计算的满足步骤（5）中的迭代终止条件，则停止迭代。

算例1 选自文献［18］，算例2 选自文献［19］，算例3 选自文献［20］。本文算法使用MT

ALB编程实现计算过程，使用高精度的MCS 作为对照，检验新算法可行性的同时，对其计算精度、效率以及稳健性进行验证。

4.1 算例1

结构中一边长为b的正方形截面轴压短柱受到的轴压为P＝1000kN，设短柱的材料强度为fc，其中b和fc服从正态分布且相互独立，它们的均值和方差分别为：μb＝300mm，σb＝6mm；
μfc＝22N／mm2，σfc＝5N／mm2。计算柱的失效概率。

此算例中结构的功能函数为：Z＝b2fc－1000。依照计算步骤，本文算法首先预设计算结果的范围误差限和迭代误差限分别为：Pε＝0.5 ×10－6，Pε＝0.5 ×10－6（c＝－6）；
设定初始均分份数N0＝5；
开始迭代。最终迭代次数m为5 次，均分份数N为80，取用区域单元的样本数为N2＝6400 个，计算结果为本文算法和MCS 计算结果与使用样本个数的分布见图5。工程中MCS被普遍认为是精确解，以其收敛时（模拟1 ×108次）的计算结果作为准确值，即因为Pε＋Pmε，佐证了该算法提出的误差范围在逻辑上是自洽的。该算法与MCS的结果对比见表1 所示。由此可见，在处理失效概率较小的算例时，相同精度下本文算法需要的样本个数远小于MCS，计算效率显著提高。

图5 算例1 计算结果的分布Fig.5 The distribution of calculation results in example 1

表1 算例1 中本文方法与MCS计算结果对比Tab.1 Comparative Method and MCS calculation results in example 1

4.2 算例2

设结构的极限状态方程为：Z＝X2－8100（X1＋其中X1服从正态分布，X2和X3服从对数正态分布，X4服从极值Ⅰ型分布，均值μX＝（60，2000，24，50）T，标准差σX＝（6.0，74.0，1.2，10.0）T。计算结构的失效概率。

该算例计算方法与算例1 相同，初始参数选取和迭代结果见表2，与MCS 的比较见表3。该算例表明本文算法在保持计算精度和效率的同时，对非正态随机变量和高次非线性功能函数有很好的兼容性。

表2 算例2 中本文算法初始参数选取和迭代结果Tab.2 The initial parameter selection and iteration of this product algorithm in example 2

表3 算例2 中本文方法与MCS计算结果对比Tab.3 Comparative Method and MCS calculation results in example 3

4.3 算例3

设结构的功能函数为：

其中，X1和X2服从标准正态分布且相互独立，P为该功能函数的参数，随着P取值的增大，极限状态方程的非线性程度也随之提高。对于算例3，选择不同的P值进行计算，同样就计算结果将本文方法与MCS 作对比，数据列入表4。为使照组的计算结果达到一定精度，MCS法对该算例采用1 ×107次模拟，平均每次模拟结果的耗时为96s。相同精度下本文方法仅需不到1s的时间。

表4 算例3 的计算结果与对比Tab.4 Calculation results and comparison in example 3

通过分析算例3 发现，本文算法除上述优点外，其对非线性程度各异的功能函数都有较好的适应性，即稳健性较强。

结构可靠度分析中，计算结果准确、工作量小是研究的主要方向。为解决MCS 面对小概率事件时求解的困难，本文基于有限单元法提出一种新算法。该算法通过压缩与均匀划分随机变量的多维积分空间实现化整为零、逐个分析；
经数值逼近计算结构失效概率。数值算例检验表明：

1.算例1 中在－6.5%的相对误差下，新算法的计算量和计算耗时可减少至MCS 的0.0001倍以下；
算例2 和算例3 中则是以－2.3%与小于6%的相对误差，分别将工作量减少至0.1 与0.01 以下。因此相同计算精度下，本文算法工作量一般可减小到MCS 方法的0.1 以下，由其当失效概率较小时，该算法的优势更明显。

2.新算法继承了MCS方法的诸多优点，即：对服从任意分布的基本变量、高次非线性、甚至不连续的功能函数都可以很好的求解，并且具有很好的稳健性。

3.由于该算法事先预设误差限，进而计算结构的失效概率，因此就计算结果可以得知其准确的误差范围。

本文方法也存在诸多不足之处，在面对“维数灾难”时，区域单元的取样将较为冗多，计算效率有待改进。但随着计算机技术的发展，该算法在高次非线性的大型复杂结构工程的可靠度分析中应有良好的应用前景。