基于改进的同伦随机有限元法的随机参数结构弹性稳定性分析

张 衡，项 煦，刘宇浩，黄 斌，曾 磊

(1.长江大学城市建设学院，荆州 434023；2.武汉理工大学土木工程与建筑学院，武汉 430070)

结构稳定性是关系到结构安全的主要问题之一[1−2]。传统的结构稳定性分析，如复合板材[3]、混凝土灌注桩[4−5]、壳体结构[6− 8]和管道[9 − 10]等的弹性屈曲荷载求解，均是在确定性结构系统的框架下进行。然而现实工程中，材料性能、几何尺寸以及边界约束总是存在一定程度的不确定性，一些研究已注意到这些不确定性将对结构的稳定性产生不可忽视的影响[4,11−12]。因此，在结构建模中考虑这些不确定性将更加符合工程实际，对结构的安全评估具有重要意义。

为研究材料不确定性(如弹性模量的随机性)对结构屈曲荷载的影响，学者提出了不同的分析方法，主要包括蒙特卡洛模拟法[13−15]和基于摄动技术的随机有限元法[16−23]。

SHINOZUKA 和ASTILL [13]利用蒙特卡洛模拟法研究了支撑于弹性地基上梁柱构件屈曲荷载的均值和方差；SONG 等[14]采用逆迭代法结合蒙特卡洛模拟实现了高斯随机场下第一阶屈曲特征值的求解，并指出当第一阶、二阶屈曲特征值为密特征值时，由于屈曲模态在结构参数变化后可能互换，需将两个屈曲模态同时作为候选向量才能保证结果的准确性。蒙特卡洛模拟法的应用不受结构规模、随机变量变异性大小以及变量维度的影响，但其求解过程十分耗时。一些学者在摄动法的基础上提出了计算效率优异的随机有限元法，应用于结构的弹性稳定性分析。

对于考虑材料不确定性的轴心受压柱，JEONG[16]用摄动随机有限元法得到了屈曲荷载的均值和方差；ZHANG 和ELLINGWOOD [17]用低阶摄动法计算弹性地基上梁屈曲荷载的均值和方差，并指出随着随机场变异性的增大，二阶摄动求得屈曲荷载统计特征值精度将降低；RAMU 和GANESAN[18]考虑弹性模量和质量密度的随机性，研究了置于Winkler 地基模型上梁结构的稳定性；GRAHAM和 SIRAGY [19]对加筋板进行了均方差为0.1 的小变异参数下的弹性稳定性分析；ALTUS 和TOTRY[20]利用函数式摄动法得到了考虑材料不确定时形式较简单结构屈曲荷载表达式；近期，GUPTA 和ARUN[21]将钢材的弹性模量考虑为高斯分布或对数分布随机场，结合无网格伽辽金法和摄动技术求解了随机参数变异性较小时轴心受压柱随机屈曲荷载的均值和方差。考虑到低阶摄动法的计算精度只在随机参数变异性较小时有保证，KAMIŃSKI和ŚWITA[22]提出用广义随机有限元法(GSFEM)分析弹性稳定性问题，该方法最高采用了十阶摄动技术求解临界应力的前四阶统计矩。但该方法中，为便于推导随机响应级数表达式的高阶项，要求随机变量服从高斯分布且为单一变量。

除随机有限元法外，学者们也提出了不确定弹性稳定性分析的其它方法。如WU 等 [23 − 24]将信息不完备的不确定性参数描述为区间量，给出了屈曲荷载的区间值，但各区间参数的波动范围不超过 ±24%；LI 等[25]指出对数正态分布比高斯分布更适合于描述弹性模量的随机性，并提出用等几何谱随机有限元法来研究考虑材料不确定时板的弹性稳定性，但随机变量的变异系数小于0.1；在响应面法的基础上，ALIBRANDI 等[26]对随机框架结构的屈曲极限状态进行了可靠性分析；徐亚洲和白国良 [27]采用概率密度演化法研究了大型冷却塔在自重、内吸力及外部风压综合作用下的屈曲承载能力，并进行可靠性分析；SU 等[28]利用格林函数法得到了材料和几何参数小变异时结构的屈曲荷载。总结来说，这些方法都是试图对随机屈曲荷载和屈曲模态进行准确而高效的求解。遗憾的是，除概率密度演化法之外，这些方法在求解随机屈曲荷载时通常要求结构参数变异性较小且为高斯分布，对于屈曲荷载的分析也局限于低阶统计矩(均值和方差)。众所周知，现实工程中结构参数常为非高斯分布且具有较大变异性；同时，仅参考低阶统计矩对于精确描述随机响应的概率分布特征是不够的。特别是当结构失稳为小概率事件时，屈曲荷载概率分布的尾部能否准确拟合就显得尤为重要。而利用高阶统计矩(如偏度和峰度)可进一步保证对随机响应概率分布(特别是尾部)拟合的准确性。因此，建立一种不局限于随机变量分布类型、能处理较大变异性并获得高精度高阶统计矩的随机屈曲荷载高效求解方法具有重要意义。

近期文献[29 − 31]基于同伦分析法提出了同伦随机有限元法(HSFEM)来处理线性和几何非线性静力响应问题[29−30]及随机特征值问题[31]。该方法继承了同伦分析法[32−33]不受小参数限制的特点，当结构参数有较大变异性时仍具有良好的求解精度，且对随机变量的分布类型没有限制。与现有常用随机有限元法如摄动随机有限元法和正交多项式展开法(PC 法)相比，在展开阶数相同的情况下，HSFEM 比摄动随机有限元法具有更大的收敛域和更优的求解精度，且有效避免了高阶摄动法可能出现的发散问题。同时，其计算工作量和求解难度显著小于PC 法：设结构自由度为N，含n个随机变量，级数取r阶展开，对于随机特征值问题，利用伽辽金法确定PC 法的投影系数，将把原N×N维特征方程扩阶为初始值未知的N(n+r)!/n!r!维非线性方程组迭代求解问题，不易保证所求解的稳定性。相比之下，HSFEM 中的系数则是通过1+rn(rn−r−n+5)/4 次N维矩阵的加、减和乘法运算求得。

在同伦随机有限元法中，结构的随机响应首先以同伦级数表示，然后，通过在随机样本空间内选取一定数量的样本来确定同伦级数的最优形式，因此，该方法的计算精度会受所选样本点的影响，且对不同问题样本点的选择需依赖一定的人工经验。故寻求一种最优同伦级数的自适应确定方法有利于大大提高该方法的稳定性。

本文基于随机残差最小化思想改进同伦随机有限元法，并将其应用于随机结构的弹性稳定性分析。首先，推导了弹性结构随机屈曲荷载和屈曲模态的同伦级数表达式，并给出了同伦级数中任意阶系数的显式递推表达。然后，定义了弹性稳定方程解的随机残余误差，通过最小化该随机残余误差实现了同伦级数的优化求解。本文方法解决了目前同伦随机有限元法计算精度会受所选样本点影响的问题，具有很好的稳定性。同时，本文方法不受参数分布类型限制，适用于结构参数变异性较大的情形，有效避免了高阶摄动法计算大变异随机参数结构屈曲荷载时容易出现发散的现象，并能获得屈曲荷载及模态的高精度高阶统计矩。通过函数算例以及变截面轴心受压柱和平面框架结构的弹性稳定性分析，验证了本文方法的有效性。

1 结构弹性稳定性分析的有限元方程

不失一般性，以承受外荷载P的确定性梁结构为例，设每个梁单元的内力为Fe，则每个有限单元的势能可写为：

式中：E、I和le分别为弹性模量、截面惯性矩以及单元长度；we为单元位移场，是单元节点位移De与形状函数Ne的乘积，可表示为we=NeDe，将它代入式(1)可得：

式中，λ 和D分别为屈曲特征值和特征向量。其中，λ 也称为荷载比例因子，λ 中的最小特征值λcr对应无缺陷结构的屈曲荷载Pcr=λcrP。

2 求解结构屈曲特征值的改进同伦随机有限元法

2.1 随机屈曲特征对的多变量同伦级数表达

本文考虑材料弹性模量的不确定性，将其描述为随机场或用相互独立的随机变量来表达，这样式(3)中的弹性刚度矩阵K可用式(4)表示：

式中：K0 为弹性模量取均值时的整体刚度矩阵；ΔK 为弹性刚度矩阵的随机部分，它包含n 个相互

独 立 的 随 机 变 量，其 向 量 形 式 为ξ={ξ1, ξ2, …,ξn}；Ki 为对应第i 个随机变量的确定性系数矩阵，可通过K-L 展开、中心点法、局部平均法等随机场离散方法得到。那么，屈曲特征值λ 和屈曲特征向量D 将是随机向量ξ 的函数，此时随机结构弹性稳定性方程为：

通常采用摄动法和正交多项式展开法对形如式(5)的含随机变量的广义特征方程进行求解。摄动法一般采用低阶摄动技术，因此，在随机变量变异性较大时计算精度难以保证；正交多项式展开法采用低阶展开表达式时的计算精度明显优于摄动法，但是在随机变量维度和结构自由度增加时，方程的求解工作量会急剧增大，同时，因方程为迭代初始值未知的高维非线性方程组，求解十分困难[34]。

廖世俊[32]最早在博士论文中提出了同伦分析法用以处理强非线性问题，该方法较传统的摄动法有着显著的优越性，已在流体力学、固体力学、应用数学和金融等领域得到广泛的应用[35]。本文基于同伦分析方法的思想，对式(5)进行求解。下面给出具体求解过程。

对于式(5)，设随机变量取均值时，对应的解为λ0和D0，同伦分析方法通过重新构造式(5)，建立起真实解λ(ξ)和D(ξ)与均值系统下λ0和D0之间的关系。首先将式(5)重新构造如下：

式 中：p∈[0, 1]、h≠0 为 辅 助 参 数；Γ(ξ;h,p)和Θ(ξ;h,p)为随机向量ξ 和辅助参数的函数。为简化表达，引入运算符号L和N，将式(6)写为：

观察式(7)可知，当p=0 时，有：

当p=1 时，有：

那么，伴随着参数p从0 逐渐增加至1，Γ(ξ;h,p) 和 Θ(ξ;h,p)由λ0和D0变化到式(5)的真实解λ(ξ)和D(ξ)，这样就建立了关于Γ (ξ;h,p) 和Θ(ξ;h,p)的同伦。将 Γ(ξ;h,p) 和 Θ(ξ;h,p)做关于参数p的泰勒展开，将有：

关于式(10)的收敛性，文献[31]讨论了一个变量的情况，并说明只要适当选择非线性方程中的线性算子，式(10)总会在p=1 处收敛；对于多个变量的情形，这个性质仍然成立，具体证明可参考文献[31]。那么，当p=1 时有：

式中，λm和Dm(m≥1)可通过以下步骤得到：

1) 将式(7)对参数p求导m次，然后令p=0可得：

具体如下：

3) 将式(14)代入式(13)可得系数Dm的显式递推表达式：

观察式(14)和式(15)可知，λm和Dm的递推表达式是由已知矩阵K0、ΔK、Kg、λ0和D0等简单的乘法与求和得到，且只涉及一次矩阵求逆(即η的计算)，故而可方便地通过MAPLE 等符号运算软件推出。将λm和Dm代入式(11)整理后，便可得到随机屈曲特征值λ(ξ)和屈曲模态D(ξ)的表达式如下：

式中：Z(ξ;h) 为 屈曲特征值或屈曲模态；Zi1,Zi1i2,···等为确定性系数，其与λm和Dm之间有如下对应关系：

式中：m=0, 1, …, ∞为式(16)的展开阶数；βm,l(h)(l=1, 2, …,m)为只含辅助参数h的函数，称为趋近函数，其中h∈(−2, 0)。

式(16)为同伦级数，该级数随参数h∈(−2, 0)取值的不同而变化，相应有不同的收敛范围。可以证明，传统的泰勒级数只是同伦级数的特解(即h=−1 时)。同时，注意到推导过程无需指定随机变量的分布类型，同伦级数可适用于任意分布类型的随机变量。

2.2 同伦级数的降维与降阶

理论上，当同伦级数式(16)的展开阶数m趋近无穷时，就是式(5)的精确解。在实际应用中，为减少计算工作量，可只取有限阶展开式；另外，考虑到当变量增多时，高维交叉项对级数精度的贡献将减小[29−30]，可忽略同伦级数中的部分交叉项以进一步提高计算效率。

如果同伦级数式(16)中保留至多含q个不同随机变量的交叉项，则称其为q维同伦级数展开，通常一维或二维同伦级数展开即可获得较满意的计算精度[30]。表1 列出了一维、二维和完整同伦级数在不同随机变量个数、不同展开阶数时，需要计算系数的个数。

表1 同伦级数中需要计算系数的个数Table 1 Number of terms in homotopy series

从表1 可见，当随机变量个数或级数展开阶数增加时，完整的同伦级数中，所需计算的系数个数将急剧增加；而采用一维或者二维同伦级数时，相应的系数计算工作量将大幅减少。具体地，当随机变量数为20 时，3 阶展开同伦级数中，一维系数个数约为完整级数的1/30，二维系数个数约为完整级数的1/3；若取4 阶展开，则一维和二维系数个数将分别降为完整级数展开的1/130 和1/9。

当通过级数截断和降维后，设同伦级数的近似解为Zˆ(ξ,h) ，该近似解与精确解Ze(ξ)之间的相对误差为：

式中，|·|为绝对值运算符。可以观察到，该误差会随着h值的调整而变化。

2.3 确定辅助参数h 值的随机残差最小化法

为使近似解Zˆ(ξ,h)的误差尽可能小，需要确定参数h的合适取值。在确定性问题中，同伦分析方法给出了三种确定h值的方法：h曲线法、残差法和比值法[33,35]，这些方法均不便于直接用于随机问题的求解。已有的HSFEM 方法采用抽样技术来确定h值[29−31]，但该方法计算精度依赖于样本值的选取，且对不同问题需要有一定的选点经验。本文提出利用随机残差最小化法来实现辅助参数h值的自动求解。

由于精确解Ze(ξ)未知，无法直接通过优化式(19)来确定h值。注意到随机弹性稳定性方程式(5)是已知的，屈曲特征值和屈曲模态的精确解λe(ξ)与De(ξ)对随机变量样本空间内任意一点，有式(20)恒成立：

式中，N为结构自由度。第j(j=1, 2, …,N)阶屈曲特征值与特征向量和在整个样本区间( ξ ∈Ω)上满足：

式中，gξ(ξ)为联合概率密度函数。将经过2.2 节降维与降阶后的屈曲特征值和特征向量同伦级数代入式(21)，则有：

式中：fξi(ξi) (i=1, 2, …,n)为 ξi的概率密度函数；E{•}为求期望运算符；wi(h)(i=1, 2, …,N)为关于参数h的标量函数。

式(22)中的W(h)描述了在随机变量全域上第j阶屈曲特征对近似解基于弹性稳定性方程产生的随机残差，是关于参数h的向量函数。注意到W(h)越小意味着近似效果越好，于是有如下优化问题：

这里推荐最小化W(h)的2 范数来确定h值，即：

因式(24)仅含一个未知量h，且任一wi(h)均为多项式形式，故式(24)的求解工作量不大。

通过以上随机残差最小化法得到参数h的优化值h*后，便可确定随机屈曲特征值和屈曲模态的优化同伦级数展开，称此方法为改进的同伦随机有限元法。同时约定，经由2.2 节降维的同伦级数中，采用只含单个随机变量时的方法称为本文方法-1，至多含两个随机变量交叉项时的方法称为本文方法-2。

3 算例分析

算例分析中将采用不同方法拟合非线性函数，求解轴心受压柱和平面框架结构的屈曲特征值及屈曲模态，以说明本文方法(包括本文方法-1 和本文方法-2)的有效性。其它方法包括泰勒级数展开法、正交多项式展开法(PC 法)，文献[17]提出的低阶摄动随机有限元法(以下简记为文献[17]方法)，文献[22]提出的广义随机有限元法(以下简记为文献[22]方法)以及基于样本点的同伦随机有限元法HSFEM，并以蒙特卡洛模拟法的结果作为参考值。所有方法均采用MATLAB 编程计算。

3.1 算例1.强非线性函数

函数Y如式(25)所示，该函数具有强非线性。分别采用泰勒级数展开法、正交多项式展开法(PC 法)和同伦级数展开法逼近。各方法均取4 阶展开，得到的表达式见式(26)～式(28)：

分别采用HSFEM 方法和本文方法确定同伦级数展开中参数h的值。在HSFEM 方法中的样本点坐标分别取为1σ、2σ 和3σ(σ=1 为随机变量的均方差)，所得h值分别为−0.96、−0.53 和−0.90，相应计算结果以HSFEM(X)表示，其中X为样本点坐标；本文方法利用随机残差最小化法得到h值为−0.86。另以蒙特卡洛模拟法20 万样本的计算结果作为标准参考值。不同方法计算Y的前4 阶统计矩结果见表2。

表2 不同方法计算Y 的前4 阶统计矩Table 2 First four order statistical moments of Y calculated by different methods

从表2 中可以看出，对于均值和均方差，四种方法均能得到满意的计算结果；对于偏度和峰度，泰勒级数结果已经发散，PC 法精度较好；HSFEM 方法的结果较泰勒级数展开法有明显改进，但是精度受所选样本点影响较大；本文方法的结果精度最高，与蒙特卡洛模拟法相比峰度误差仅为8%。

为进一步说明本文方法较HSFEM 方法的优势，绘制了随机残余误差式(24)随h值变化的关系曲线如图1 所示。图1 中，h=−1 处即为泰勒级数展开对应的随机残余误差。可观察到，HSFEM方法和本文方法的误差较泰勒级数展开法有明显减小，但HSFEM 方法所选样本点并不能保证找到h值的最优点，且样本点坐标的选取依赖于经验，本文方法则能自动找到全域上的误差最小值。

图1 随机残余误差关于辅助参数h 值变化曲线Fig.1 Stochastic residual error with respect to the auxiliary parameter h

3.2 算例2.变截面轴心受压柱

如图2 所示变截面轴心受压柱，AB和BC两段长均为3 m，在柱的端部有一集中力荷载P=1 kN。柱的有限元模型被划分为12 个单元，每个单元节点含横向侧移和转角两个自由度。下面分两种工况讨论AB段和BC段抗弯刚度随机性对该柱屈曲特征值的影响。

图2 变截面柱 /mFig.2 Variable cross-section column

工况1：假定AB段抗弯刚度为确定值EIAB=16 kN·m2，考虑BC段抗弯刚度的随机性并以式(29)表达：

图3 变截面柱屈曲特征值前4 阶统计矩(工况1)Fig.3 First four order statistical moments of the buckling eigenvalue for Case 1

从图3 可观察到，由于本工况中抗弯刚度变异性较大，随着展开阶数的增加，文献[22]方法只是对于均值存在收敛的趋势，对于均方差、偏度和峰度，随着展开阶数的增加计算结果反而会发散；相比之下，基于同伦级数展开的HSFEM(7σ)和本文方法得到的各阶统计矩值在4 阶展开后均趋于稳定；对于高阶统计矩，如峰度，HSFEM(7σ)在8 阶展开后可收敛，而本文方法只需要4 阶展开即可，说明本文方法在计算随机弹性屈曲特征值时具有良好的收敛性。

工况2：考虑多维随机变量，AB段的抗弯刚度为确定值EIAB=16 kN·m2。BC段每个单元抗弯刚度为对数正态分布并以式(30)表示：

图4 变截面柱屈曲特征值前4 阶统计矩(工况2)Fig.4 First four order statistical moments of the buckling eigenvalue in Case 2

图4 表明：与蒙特卡洛模拟法的结果相比，文献[17]方法在计算均值和均方差时有较好的精度，但是对于偏度和峰度在 δ大于0.3 后会有显著的误差；如果采用高阶摄动(4 阶摄动)，其求解精度不仅没有改进，反而在 δ大于0.5 后迅速发散。基于同伦级数展开的HSFEM(7σ)、本文方法-1 和本文方法-2 均能得到更加准确的结果，特别是对偏度和峰度这些高阶矩信息的求解没有出现发散现象。同时，本文方法计算的偏度值比HSFEM(7σ)更好。

本文方法中h值和趋近函数 βm,l(h) 随 δ变化的关系如图5 所示。从图5(a)可见，随着随机变量变异性的增加，h值由接近−1 逐渐变化至−0.6，HSFEM(7σ)中确定的h值在随机残差法得到h值附近波动，说明对于变异性不同的随机变量，需采用不同的样本点才能取得最优h值；图5(b)中，βm,l(h)、sigβm,l(h)和dou βm,l(h) (m=4,l=1, 2, 3, 4)分别为HSFEM(7σ)、本文方法-1 和本文方法-2 同伦级数中的趋近函数系数。可见，随着随机变量变异性的增大，3 阶、4 阶趋近函数值逐渐接近于0，可以显著避免高阶展开系数带来的发散现象。

图5 辅助参数h 值和趋近函数随δ 变化关系图(工况2)Fig.5 Values of the auxiliary parameter h and the approaching function in the homotopy series expansion for Case 2

3.3 算例3.平面框架结构

如图6 所示7 层平面框架结构，每根框架柱受P=150 kN 轴向荷载。该框架有限元建模中，每隔1 m 将梁柱划分为一个梁单元，每单元3 个自由度，模型共有383 个单元，自由度为1059 个。首先进行内力分析计算图示荷载作用下各单元应力，进而确定各单元几何刚度矩阵。

图6 7 层框架结构 /mFig.6 A 7-story frame structure

1)考虑6 根框架柱弹性模量的随机性，形成6个相互独立的对数分布随机场，每个随机场表示为：

式中，ξi为第i个标准高斯分布随机变量，指数部分是协方差核满足式(32)的高斯随机场：

式中：σ为高斯随机场的均方差；y1和y2为沿Y方向的坐标值；相关长度l=10 m。6 根柱子弹性模量随机场中的μc均为 5.327。为将该对数分布随机场转化为级数展开的形式，可用正交多项式展开法(PC 法)来逼近式(31)，具体表达式如下：

式 中：ξ={ξ1,ξ2,···,ξ4}，gk(ξ)为 关 于 随 机 变 量ξi(i=1,2,···,4 )的正交基；lk(x)为各正交基相应的确定性系数，其求解方法可参见文献[36]。

框架柱随机场的变异系数 δEc分别假设为 0.30和0.60。表3 和表4 分别列出了用摄动法、HSFEM和本文方法采用2 阶～4 阶展开时求出的屈曲特征值前4 阶统计矩值，其中HSFEM 选用样本点坐标为7 σ，并以蒙特卡洛模拟法(10 万样本)计算结果作为参考值。

表3 框架结构屈曲特征值的前4 阶统计矩( δEc=0.3)Table 3 First four order statistical moments of the buckling eigenvalue ( δEc=0.3)

表4 框架结构屈曲特征值的前4 阶统计矩( δEc=0.6)Table 4 First four order statistical moments of the buckling eigenvalue ( δEc=0.6)

从表3 可看出，在随机场变异性较小的情况下，摄动法、HSFEM 和本文方法求得的4 阶矩值都与蒙特卡洛模拟方法很吻合，随着展开阶数的提高，计算精度也越高，且HSFEM 和本文方法计算的偏度和峰度略优于摄动法。

当随机场变异性增大时，参见表4，可发现2 阶展开的摄动法计算的峰度和偏度相对误差超过20%。若提高展开阶数，摄动法的计算精度并未提高，反而出现了发散现象，说明高阶摄动并不一定适用于变异性较大时的非高斯分布随机变量。HSFEM 的结果较摄动法虽然有了明显的改进，但偏度和峰度的结果仍有较为显著的误差(特别是3 阶展开时)，原因是本算例中样本点坐标在随机变量变异性不同时、或展开阶数不同时均统一选取为7 σ。换言之，HSFEM 需针对不同变异性大小和展开阶数提出更为具体的选点原则才能保证该方法计算精度的稳定性。同时，可发现在随机变量变异性较大时，HSFEM 精度对所选样本将更加敏感。

而对于本文提出的两种降维方法，随着展开阶数的增加，仍能保证计算精度的提高，说明本文方法具有良好的收敛性和很好的稳定性。

图7 给出了框架柱随机场变异系数0.6 时，由4 阶展开本文方法-1 计算的框架结构第1 阶屈曲模态的均值、均方差、偏度和峰度。从图7 中可观察到，对屈曲模态的前4 阶统计矩，本文方法的计算精度都很高。

图7 一阶屈曲模态前4 阶统计矩Fig.7 First four order statistical moments of the 1st-order buckling mode

表5 是摄动法和本文方法采用不同展开阶数时计算屈曲特征值所需要的时间。本算例中共考虑了24 个随机变量，由于本文方法采取降维策略舍弃了交叉项，在高阶展开时其计算量将少于同阶展开的摄动随机有限元法，同时远少于蒙特卡洛模拟法的计算量。而HSFEM 方法需要额外进行24 次确定性有限元分析以确定h值，此部分工作量需耗时约170 s。

表5 各方法计算时间 /sTable 5 Computation time for different methods

2)除考虑6 根框架柱弹性模量的随机性外，同时考虑横梁弹性模量的随机性。设该7 层框架每层横梁的弹性模量可表示为式(34)的形式：

表6 框架结构屈曲特征值的前4 阶统计矩( δEc=0.6, δb=0.3)Table 6 First four order statistical moments of the buckling eigenvalue ( δEc=0.6, δb=0.3)

4 结论

本文基于同伦分析法推导了屈曲特征值和屈曲模态的多变量同伦级数表达式，并给出级数中任意阶系数的显式递推关系式；借助关于弹性稳定性方程的随机残差最小化法，确定了随机屈曲特征值和屈曲模态同伦级数解的最优形式。

函数数值算例表明：本文方法能够实现随机响应同伦级数表达式中辅助参数h值的自动寻优，解决了已有同伦随机有限元法HSFEM 结果易受所选样本的影响和不稳定的问题。并且本文方法结果比同阶泰勒级数展开结果精度更高。

变截面轴心受压杆和框架结构的稳定性分析数值算例表明：相比于现有的各类低阶或高阶摄动法，本文方法求解的随机屈曲特征值和屈曲模态均展现出优异的计算精度，相比于蒙特卡洛模拟法计算效率大大提高。同时发现，当随机参数变异性较大时，利用高阶摄动法得到的高阶统计矩(如偏度和峰度)可能比低阶摄动更差，而本文方法十分有效地避免了高阶展开级数解可能出现的发散现象，体现了本文方法的稳定性。此外，对于平面框架结构，当同时考虑柱和横梁弹性模量的随机性时，本文方法中交叉项对屈曲特征值计算精度的贡献已不可忽略，建议采用至少含两个随机变量交叉项的同伦级数求解屈曲特征值的高阶统计矩。

本文方法的提出为随机结构的弹性稳定性分析提供了一种新的高效求解途径。考虑到现有同伦随机有限元法已被应用于几何非线性问题的求解，本文方法也可方便拓展到处理随机结构非线性静力及几何非线性稳定问题。