1x1i?(L1i/2)x2i?(L2i/2)S(x1,x2)dx1dx2 (7.3.11) 定义为 Si?SAi(x1i,x2i)??Ax1i?(L1i/2)?x2i?(L2i/2)如果S(x1,x2)是一个均匀随机场，则可用均值m、方差σ及归一化协方差函数ρ（r1,r2）近似描述，r1和r2分别为沿两个方向的距离。对应的局部平均随机场可用E(Si)、Var(Si)
2
及互协方差Cov(Si,Sj)近似描述。E(Yi)、Var(Yi)及Cov(Yi,Yj)可由m、σ及ρ(r1，r2)计算获得。第i和第j单元局部平均Si、Sj的互协方差可表示为
2
1?2 Cov(Si,Sj)?TxiTyiTxjTyj4??(?1)k?0l?033k?l(T1kT2l)2?2(T1k,T2l) (7.3.12)
式中：T1k、T2l(k,l?0,1,2,3)的约定如图7.3.1所示。
yT1lT1lT1lT1lTyjT1lT1lT1lTyiTxiTxjT1lx 图7.3.1 二维局部平均单元
如果有限元的网格已划分，且单元总数为n，随机场实际上被离散成n个随机变量，这n个随机变量的统计特性可由E[Si]、Var[Si]及互协方差Cov[Si，Sj]反映。由于有限元网格的疏密是由应力梯度决定的，与随机场无关，当单元数很多时，随机场可另划网格，网格的疏密可由相关距离δS决定。当然，随机场网格越密精度越高。随机场的局部平均法由于对原始数据的要求低、收敛快、精度高，是随机有限元计算中最常用的方法。 ● 随机场的插值
Liu W. K.提出随机场的插值法。该法将随机场在单元内的值用单元结点处值的插值函数来表示，于是随机场的统计特性可由各单元结点处随机变量间的统计特性近似反映。
利用形函数Ni(X)，随机场b(X)离散式表示为
b(X)??N(X)b (7.3.13)
iii?1q式中：X表示空间位置；bi为随机场在结点i处的值（i=1，2，?，q）；q为单元结点数。随机场b(X)在单元内的均值和方差可表示为 E(b(X))??N(X)E(b) (7.3.14)
iii?1q Var(b(X))?
i,j?1?N(X)Niqj(X)Cov(bi,bj) (7.3.15)
6
随机场的插值法将原连续状态的随机场仍离散成一个连续函数，未直接计算随机场引起的单元间的相关性，只需给定随机场在各结点上的值，计算相对简单，易于考虑非线性和非均匀随机场问题。但需要已知相关函数，并且要求随机场对空间参数具有较高的连续性。 ● 随机场的加权积分法
Takada、Shinozuka及Deodatis提出随机场加权积分方法。该法在单元刚度矩阵的推导过程中采用随机场在单元高斯点上的加权积分，以表征单元上的随机场。
(e)假设单元e的弹性模量为 E(e)(X)?E0[1?f(e)
(e)
(e)(X)] (7.3.16)
式中，E0(X)为弹性模量的均值；f(X)为一维零均值均匀随机场，其值域为
(e)
-1+η≤f(X)≤-1-η (0＜η＜1) 两端铰接杆单元的刚度矩阵可近似表示为
(e)(e)(e)(e)
K≈K0+XΔK0 (7.3.17)
(e)
式中，K0为弹性模量取均值时的单元刚度矩阵，而
X(e)??f(e)(x)dx （L为杆长） (7.3.18)
0L固结杆单元的刚度矩阵可表示为
(e)(e)(e)(e)(e)(e)(e)(e)
K≈K0+X0ΔK0+X1ΔK1+X2ΔK2 (7.3.19) 其中： X1(e)(e)?x?f(e)(x)dx,X2?x2?f(e)(x)dx
00LL式（7.3.17）实质上是在考虑弹性模量随机场的情况下，将单元刚度矩阵分解成确定性部分
和随机部分。式（7.3.18）等表征随机场在各单元的均值和方差，而单元间弹性模量的相关性则由下式表示 Cov(X(e)1,X)??(e)2L(e1)0?L(e2)0x1x2f(e1)(x1)f(e2)(x2)dx1dx2 (7.3.20)
不难看出，局部平均法是加权积分法的特例（即权系数全部相同）。由于采用加权积分，其计算精度相对较高，而且该法积分只需一次进行，刚度矩阵的波动性也由此得出，因此计算效率也较高。
● 随机场的正交展开
Spanos和Ghanem提出的随机场正交展开法，将材料特性参数随机场进行Karhumen-Loeve正交展开，并由此推导出刚度矩阵的级数展开式，从而获得位移、应力的统计特性。
设随机场为 S(x)?S(x)??bn?0?n?n?n(x) (7.3.21)
(x)dx
其中 bn?1?n?S(x)?Ln式中：λn、υn(x)分别为随机场S(x)相关函数的特征值和特征函数。υn(x)具有如下正交性： ?m(x)?n(x)dx??mn（Kronecker函数） (7.3.22) 式（7.3.21）对于任何分布的S(x)均收敛，可取至第r阶以满足精度要求。对于一维杆单元，单元刚度矩阵可表示为 k?k?ee??bknn?0ren
7
其中 kn??e
e
e
enTee?(x)B(x)P(x)B(x)dx n?L式中：P(x)=D(x)/S(x), D(x)为单元弹性矩阵。集成整体刚度矩阵，得 K?K?r?bKnn?0?1nrn (7.3.23)
从上式，可得位移： U?(I??bn?0KK)KF (7.3.24)
?1n?1利用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均为正态分布，其均值和标准差分别为
?S，此时z=R-S也是正态分布的随机变量，具有均值z?R?S和标准差R、S和?R、22?z??R??S。其概率密度函数为
2???11z?z??? fz(z)?exp?????2??z?2??z???????z?? (7.4.1)
其分布图示于图7.4.1，阴影部分是结构失效概率Pf，非阴影部分面积即结构的可靠度Pr。 xfz(Z)2均值点失效边界PrPf可靠区改进一次二阶矩法均值一次二阶矩法失效区O x1 图7.4.1 正态分布概率密度函数 图7.4.2 失效边界 8
在工程实践中R和S不一定是正态分布，但可以变换成标准正态分布。他们的统计特征量是均值μ、标准差σ、相关偏度或变异系数等，变异系数ν=σ/μ，表示随机变量相对于均值的变异。目前工程上一般用无量纲的可靠指标β来反映结构的可靠度，??z/?z，β越大，失效概率Pf越小，其互补的可靠度Pr就越大。
一次二阶矩法采用只需已知均值和标准差的数学模型去求解结构的可靠度。此法将功能函数z=g(X1,X2,?,Xn)在某点用Taylor级数展开，并近似地取一次项使极限状态方程线性化，然后求得可靠指标β。
*
改进的一次二阶矩法将线性化点选在失效边界上，而且选在结构最大可能失效点P上
**
（如图7.4.2）。选择设计验算点P(Xi│i=1,2,?,n)作为线性化点X0i时，线性化的极限状态方程为 z?g(X,X,?,X)?则z的均值为
*1*2*n?(Xi?Xi*)i?1nn?g?Xi?0 (7.4.2)
X*?z?g(X,X,?,X)??(?Xi?Xi*)*1*2*ni?1*
*
*
*
?g?Xi
X*由于设计验算点P选在失效边界上，有g(X1,X2?,Xn)=0，因此μz成为 ?z??(?i?1n*?XXii)?g?Xi
X*但随机变量Xi互不相关时，z的标准差σz为 ??n??g?????Xi?i?1?????Xi??X*?2z????12????i?Xi?g???Xii?1?nX*?
???其中
?Xi?i??g?XiX* ????2?n?????Xi?g?Xi?i?1???X*????12称αi为灵敏度系数，表示第i个随机变量对标准差的相对影响。于是可靠指标β为
???z??z?(?Xi?Xi*)i?1nn?g?XiX* (7.4.3)
??g?????iXi?Xii?1?X*????变换上式为
?g?i?1?Xi*n?(?Xi?Xi*???i?Xi)?0
X*即 Xi??Xi???i?Xi (i=1，2，?，n) (7.4.4) 上式中，μXi，σXi为已知的各随机变量的均值和标准差，待求的量为Xi和β，可迭代求解。
● 随机有限元的一次二阶矩法
T
设可靠性分析中的一组基本变量X=(X1,X2,?,Xn)为相互独立的正态变量（不满足时可
*
9
作变换），并进一步变换为标准正态变量Y=(Y1,Y2,?,Yn)，其中Yi=(Xi-μXi)/σXi。于是功能函数也可转换到标准正态空间，即
-1-1
g(X)=g(R(X),S(X))=g(R(T(Y)),S(T(Y)))=G(Y) 失效概率Pf=1-υ(β),由β的几何含义可知，β值为在标准正态空间中从原点到失效面的最短距离。在标准正态空间中概率密度是关于原点（均值点）旋转对称，并且随着到原点的距离的平方呈指数下降。在一次二阶矩法中，标准正态空间中的极限状态面（失效面）被一个到原点最小距离点处的切平面代替，也即按一阶Taylor展开式把功能函数在设计点处
*
展开，使之线性化。采用迭代法可以近似确定极限状态面G(Y)=0上距原点最近的点Y,然后按距离公式确定结构的可靠指标。具体的迭代格式如下：
?G(Y(i))???Y(i)nT
Yj(i?1)??G(Y)?Yj?G(Y)j?1?Yjn?Y(i)Y(i)2??G(Y)????j?1??Yj?Y(i)??(j?1,2,?,n) (7.4.5)
??Y(i)T?Y(i) (7.4.6) 以Mises屈服准则为例，功能函数 g(S)=s0-(sxx-sxxsyy+syy+3sxy)
于是极限状态方程 g(S(X(Y)))=G(Y)=0 (7.4.7)
计算可靠指标β时要用到功能函数G(Y)的梯度向量，由求导链式法则
2
2
2
2
?G(Y)?g(S)?S(X)?X(Y)??? (7.4.8) ?Y?S?X?Y由于S(X)难于用显式表达，求?S(X)/?X存在困难。有效方法是利用随机有限元一次二阶矩法计算结构的可靠指标β和设计点Y。步骤如下：
（1） 确定随机变量转换关系Y=Y(X)；
(0)
（2） 给定初值Y=0，计算式（7.4.8）右边的三个偏导数；
*
（3） 按式（7.4.5）和（7.4.6）计算设计验算点坐标Y；
（4） 用随机有限元法计算S和(?X(Y)/?Y)Y(i)；
(i)
*
（5） 重复（3）、（4）两步骤，直至收敛（G(Y)=0）。 （6） 再按式（7.4.5）和（7.4.6）计算可靠指标β值。 ● 随机有限元的最大熵法
一次二阶矩法的计算量随迭代次数成倍增加，使该法的使用受到限制。而最大熵法用于结构的可靠性分析时，可根据随机变量的二阶矩来拟合概率密度曲线。因此随机有限元结合最大熵法可用于求结构响应的概率密度曲线，从而计算结构的失效概率。
熵被定义为信息的均值，信息是对个别X值不确定性的度量。不确定性越大，熵也越大。对于一个连续随机变量，熵为 S??n(i)
?Rf(X)ln?f?X??dX，而对于离散随机变量，熵为
S???f(xi)ln?f(xi)?，这里，f(X)是随机变量X的概率密度函数；f(xi)是离散点概率
i?1函数。最大熵法通过调整概率密度函数f(X)使熵S取得最大值，并满足约束条件：
?
Rf(X)dX?1和?RXif(X)dX?m ，利用最大熵法可得近似的概率密度函数f(X)。为
10