了求得最大值，可引入拉格朗日乘子法，把问题转化为确定拉格朗日乘子。
选用基于Neumann级数的随机有限元法，写出递推方程。可将局部平均理论引入Neumann随机有限元法，以进一步提高计算效率。随机有限元最大熵法计算结构失效概率的步骤如下：
(1) 随机场离散，随机变量化为随机向量，并计算协方差矩阵； (2) 利用特征正交化技术，用一组统计独立的随机变量描述随机场； (3) 利用Neumann随机有限元法计算响应量的前若干阶矩； (4) 拟合所有节点或部分节点的最大熵密度函数； (5) 利用数值积分计算结构的失效概率。 ● 算例
例1 如图7.4.3所示桁架下端受一集中力，假设各设计变量都是均匀分布的随机变量，
622
三根杆弹性模量E，μE=2.0×10N/cm，Cov(E)=2.5%；截面积A，μA=0.8cm，Cov(A)=2%；杆2长度L2,μL2=20cm，Cov(L2)=1.5%；杆1和杆3的长度均为L1，μL1=28.284cm，
2
Cov(L1)=2.3%；外力P，μP=4000N；Cov(P)=10%；材料屈服应力σy，μσy=3200N/cm， Cov(σy)=7%。当极限状态方程取为 σy-S=0（S为构件中最大应力值）时，试求该结构的失效概率。
图7.4.3 简单超静定桁架 解： 用随机有限元一次二阶矩法计算，假设由各结点坐标值变异引起的方向余弦的变异很小，形成确定矩阵，则由随机有限元递推方程组可得各结点位移分布特征值。计算得结点4的位移均值和标准差为μδ4=2.2928cm和σδ4=0.0037cm；各杆件中最大应力为杆元2的
22
应力，其均值和标准差为μσ2=2928N/cm和σσ2=377N/cm。由一次二阶矩法可计算得该结构的失效概率为0.23（由于结构超静定，杆2失效并不意味着整个结构失效），与蒙特卡罗有限元法计算结果比较接近（蒙特卡罗取样数为40时，失效概率达0.224）。
随机有限元法在分析混凝土重力坝可靠度方面已得到较多应用。研究者采用抽样方法对重力坝施行随机自动剖分，研究其可靠指标的统计规律和收敛于真解的规律，并进行结构优[25]
化。在对碾压混凝土重力坝抗滑稳定的体系可靠度研究中，还采用刚体极限平衡法、线弹性随机有限元法和三维弹塑性随机有限元法进行计算比较。
§7.5 随机有限元动力分析方法
● 不确定结构的自振特性
假设某一梁-柱结构承受随机分布的轴向静荷载N(x)，沿长度有微小扰动，同时在端部有轴向推力p?p(1?c),c为随机变量，其方差为?c，轴向分布荷载可表示为
2N(x)?N?NP(x)
式中p(x)是一维、均值为零的均匀随机场，方差是σP，相关偏度是ΘP，互相关函数是RPP(ξ)。结构弹性模量和质量密度也随机变化，分别表示为
2
11
E(x）?E?1?a?x??和m?m?1?b?x??，其中E和m分别是弹性模量E和质量密度m的均
值，方差分别为σE和σm；随机变量a(x)和b(x)是两个一维互不相关的、均匀的、均值为零的随机场。自相关函数分别为Raa(ξ)和Rbb(ξ)或等价的功率谱密度函数为Saa(ω)和Sbb(ω)，相关偏度分别是ΘE和Θm。
ee
任一点横向位移的有限元模式为 w(x,t)=Nq，
eTe
其中，单元位移列阵q=[w1,θ1,w2,θ2],形函数N取三次多项式插值。单元的动能表达为：
2
2
1le1??w? Te??m(x)Ae?dx?qeTMeqe (7.5.1) ?202??t?式中：M?M?M(?)?面积。单元的应变能为
eee2?le0m(Ne)TAeNedx?m?b(x)(Ne)TAeNedx，Ae为单元横截面
02le??2w?1le1eTeee?dx?qKq (7.5.2) U??E(x)I????x2?202??式中：K?K?K(?)?单元的外力功为
eee?le0eeEINxxdx?EI?a(x)Nxx0??2Tle??TeNxxdx，I是惯性矩。
1le?Qx???w?1eeTee1eTee W???F1e?dx?F1qKGcq?QqKGDq (7.5.3) ???02?l???t?22eC式中：KeGc??N0le??eTxNdx,KexeGDee??Nx(L2le/l)Nxdx,Qx/l?(Qx/l)?Qxp(x)/l，p0le??T为任一随机变量
据此集成结构总能量和总外力功，利用Hamilton原理
?t2t1?(T?U)dt???WCdt?0，
t1t2即对总位移向量q取变分，考虑到q的任意性，可得系统运动方程，进而得到如下特征问题方程： [?s2M?(?K)?pKGc?QKGD]q?0 （7.5.4） 式中，刚度系数kij?kij?kij(?)，其中kij是确定性分量，kij?EI随机扰动项：kij(?)?EI?le0Ni??(x)N?j?(x)dx，而
?le0a(x)Ni??(x)N?j?(x)dx。于是可以写出刚度系数、质量系数以及
[26]
几何刚度系数的均值和方差的表达式。在总刚度矩阵中两个刚度系数kij和krs的互协方差，
当单元尺寸相等时，利用方差函数可作进一步简化。
利用局部平均理论，在每个单元上E、m、Q的均值为零，方差为：
22Var(Ee)??E?E(le)??E?e/le Var(me)??m?m(le)??m?m/le (7.5.5)
2Var(Qe)??2p?p(le)??p?p/le22式中，?E(le)、?m(le)、?p(le)和?E、?m、?p分别为随机变量a(x)、b(x)和p(x)的方差函数和相关偏度。利用协方差函数写出两刚度系数或两质量系数的互协方差。
12
最后，特征值问题的方程可表达为
[K?K(?)?PKGG?PKGC(?)?QKGD?QKGD(?)]x??[M?M(?)]x (7.5.6)
或 K*x=λM*x 。相应平均问题的方程为
?K?PKGC?QKGDx??Mx (7.5.7)
? 利用单元刚度、质量和几何刚度之间的协方差的表达式，等效刚度矩阵K*的协方差矩阵能用独立的a(x)、b(x)和p(x)来表示。解未扰动的特征值问题（平均值问题）可得特征值的均值。对于满足小扰动情况，可推导出求两个特征值之间的协方差公式。在建立了材料特性和几何尺寸有随机扰动的梁-柱的随机特征值问题的随机有限元公式后，利用局部平均理
[24]
论，随机过程由均值、方差和相关偏度来定义。上述公式可推广到二维或三维。 ● 线性结构系统的瞬态响应
设结构系统的质量矩阵M是确定性的，阻尼、刚度和外力的概率分布由广义协方差矩阵Cov(bi,bj)（i,j=1,2,?,q）表示，线性结构系统的运动方程为
?(b,t)?C(b)x?(b,t)?K(b)x(b,t)?F(b,t) (7.5.8) x M?式中，质量阵M、阻尼阵C(b)和刚度阵K(b)都是n阶方阵，位移向量x和荷载向量F为n阶列阵。用Taylor级数对随机向量b展开，并保留二阶项，则位移向量x(b,t)关于b的二
1阶摄动式为 x(b,t)?x(t)???xbi(t)?bi??2i?1q2i,j?1?xqbibj(t)?bi?bj (7.5.9)
其中，γ是小参数；b是b的均值；x(t)、xbi(t)和xbibi(t)分别表示位移向量的均值、位移向量相对于bi的一阶差商并在b处取值和位移向量相对于bi、bj的二阶差商并在b处取值；Δbi是bi与bi之差。同样可以对C(b)、K(b)和F(b,t)写出关于bi的二阶摄动展开式。将这些展开式代入运动方程（7.5.8）式中，归并γ，γ和γ项后，得到
0
2
?(t)?Cx?(t)?Kx(t)?F(t) 零阶方程：（7.5.10） M?x?bi(t)?Cx?bi(t)?Kxbi(t)?F1bi(x,t)一阶方程： M?x式中，F1bi(x,t)?Fbi(t)?Cbix(t)?Kbix(t)(i?1,2,?,q) (7.5.11)
??(i?1,2,?,q) (7.5.12)
?2(t)?Cx?2(t)?K~二阶方程：M?xx2(t)?F2(x,xbi，t) (7.5.13)
q??1?~?式中：F2(x,xbi,t)????Fbibi(t)?Cov(bi,bj)??i,j?1??2?~~~??1?1??(t)?Kbibix(t)?Cbix?bi(t)?Kbixbi(t)?Cov(bi,bj)?????Cbibix2?i,j?1??2?q (7.5.14)
1q~x2(t)??xbibi(t)Cov(bi,bj) (7.5.15)
2i,j?1
13
~求解递归微分方程组（7.5.10）～（7.5.15）式，可同时得到x(t)、xbi(t)、x2(t)。为节省
计算量，可利用特征正交性将协方差矩阵Cov(bi,bj)转化成不相关的方差对角阵Var(ci)。
在得到零阶方程的振型矩阵后，可使上述递归方程解耦。
假定阻尼与刚度成正比，对应于k阶模态的阻尼比是ξk，则各阶解耦的递归方程组为
2?k?2?k?kx?k??k零阶方程：?xxk?fkk?1,2,?，n (7.5.16)
kci?kci?2?k?kx?kci??kxkci?f一阶方程：?x式中，xkci是xci的第k个分量；fkci2k?1,2,?,n;i?1,2,?,r (7.5.17)
且fci??TF1ci(t)，i=1,2,?,r 是fci的第k个分量，
k22?k2?2?k?kx?k2??k二阶方程：?xxk2?fk?1,2,?,n (7.5.18)
式中，xk2是x2的第k个分量；fk2是f2的第k个分量，且f2??TF2(t)。
由递归方程组（7.5.16）～（7.5.18）式解得x、再变换成x(t)、xci和x2，xci(t)和x2(t)，然后得到位移、应变和应力的概率分布、均值、方差和互协方差。
● 随机有限元法的计算机实施
不确定结构系统动力响应分析的随机有限元法的计算机实施过程，包含对前面导出的一系列递归方程组的时间积分，可用Newmark-β法、Wilson-θ法的直接积分法计算。为了给出均值和方差，需要积分的方程数目共有(q+2)个。所有积分过程都用了相同的有效刚度矩阵，对运用并行计算方法十分有利，也说明随机有限元法在动力分析中的效用。
对于概率非线性系统，仍可采用隐式或显式的时域直接积分法。对均值方程求解，可用Newton-Raphson迭代进行。结构在具有时间随机性荷载作用下的动力响应分析，具有重要的现实意义，如受地震荷载、风荷载或车辆荷载作用，在某给定时刻荷载的一些物理参数还存在不确定性，因此有必要在n维状态空间建立结构动力行为的控制方程。假定施加于系统的激励在时间历程上是零均值的正态随机向量，则响应在时间历程上也是零均值的。利用随机有限元计算的结果，均值是二阶精度，方差是一阶精度。随机有限元法可以很方便地利用已有的有限元分析程序，因此应用前景十分广阔。
14