工程结构在使用过程中会受到多种因素影响，其结构参数在空间上往往存在变异性，称这些参数为随机参数，对具有随机参数的结构进行分析即为随机分析．目前，学者们已经将随机分析应用于多种工程结构，文献[1]在考虑空间变异性的条件下研究了随机损伤演化引起的结构倒塌．文献[2]将随机参数引入到温度场中，分析了寒冷地区宽路基高速公路的热积累效应．文献[3-4]认为钢轨的不规则性是由多个随机结构参数引起的，并将这种不确定性引入到车桥耦合振动系统的动力分析中．文献[5-6]考虑土体的空间变异性，对边坡进行了可靠性分析．从以上研究可以看出：随机分析方法主要包含随机场的离散和随机变量的计算方法两部分．若要提出一种新的随机分析算法，就要从以上两个方向切入．Karhunen-Loève (K-L)级数展开法是当前最常用的随机场离散方法，具有高效、准确、不与有限元离散耦合等诸多优点，但求解第二类Fredholm积分方程的步骤较为繁琐．求解该方程的方法主要分为解析法与数值法：解析法由于其数学格式严谨，只能求解具有特定类型协方差函数的积分方程，如指数型、三角型等；而数值法则无此限制，其中Wavelet-Galerkin法[7]将信号处理中常用的小波变换技术与传统的Galerkin法结合，实现了高效、准确地求解第二类Fredholm积分方程，且不受协方差函数类型限制，被广泛应用．伴随着随机场理论的发展，随机有限元法也被越来越多的学者研究．最初，人们将蒙特卡罗法与有限元法进行结合，通过大量的抓取样本反复计算获得结构响应的统计参数，但计算效率过低，往往被用来验证其他方法的准确性．与其相对的，是基于数学和力学方法进行公式推导的摄动随机有限元法，但这种方法须取系统矩阵对随机变量的偏导，分析复杂问题时，其应用往往受限，且分析大变异结构时精度较低．针对以上缺点，文献[8-10]将摄动随机有限元法加以改进并推广到桁架结构、金属材料等领域．文献[11-12]提出了改进的摄动随机有限元法(MPSFEM)，跳过了求偏导的步骤，在精度上有提升，并且针对不同类型的随机问题，计算格式统一，无须重复推导，提升了计算和编程的效率．QR法[13]是一种非常有效的框架结构计算方法，该方法与有限元法计算格式相通，组装刚度矩阵的计算过程更为简洁，因此本研究引入QR法来替换MPSFEM的有限元计算格式，再将其与K-L法进行结合，提出了一种框架结构随机分析的新算法，即基于K-L级数的改进摄动随机QR法(简称为KLSMPSQRM)．1 K-L级数展开法K-L法将随机场分为确定部分和随机部分，其中随机部分用一组不相关的标准正态分布随机变量与特征值和特征函数的乘积来描述．假设X(x,θ)为一个平稳均匀的一维高斯随机场，则其K-L展开为X(x,θ)=X¯(x,θ)+∑i=1Mξi(θ)λifi(x)，式中：X¯(x,θ)为随机场的均值；x为空间坐标；θ为概率空间坐标；M为截取项数；ξi(θ)为一组不相关的标准正态分布随机变量，均值为0，标准差为1；λi和fi(x)分别为特征值与特征函数．λi和fi(x)可以通过求解下面的第二类Fredholm积分方程得到∫ΩC(x1,x2)fi(x2)dx2=λifi(x1)，式中：Ω为随机场区域；C(x1,x2)为协方差函数；x1和x2为随机场中两个不同位置的空间坐标．采用Wavelet-Galerkin法来求解积分方程，截取项数采用比能因子[14]判断，即η=∑i=1Mλi/(Ωσ2)，式中σ为随机场标准差，可由均值乘以变异系数得到．当η≥95%时，截取前M项特征值与特征函数．2 QR法与改进的摄动随机有限元法2.1　QR法如图1所示的框架结构，在水平方向上对其进行样条结点离散，以步长h划分为N份，在高度方向上采用正交多项式描述，因此该结构的位移场为u=∑m=1R∑i=0Nφi(x)aimXm(y)=∑m=1Rφ(x)amXm(y) ;v=∑m=1R∑i=0Nφi(x)bimYm(y)=∑m=1Rφ(x)bmYm(y) ;θ=∑m=1R∑i=0Nφi(x)cimHm(y)=∑m=1Rφ(x)cmHm(y),式中：u为水平位移；v为竖向位移；θ为转角位移；φi(x)为x方向的样条函数；R为多项式级数，一般与框架结构的最大层数相同；Xm(y)，Ym(y)和Hm(y)为反映高度方向变形规律的正交多项式，当柱底为固定支座时10.13245/j.hust.210704.F001图1框架结构的样条离散Xm(y)=Ym(y)=Hm(y)=∑n=1m(-1)n-1(m+n)!(m-n)!(n+1)!(n-1)!yHn .假定转角位移以逆时针为正，则上式可用矩阵形式表示为u=Nδ，(1)式中：δ为QR法的基本未知量，是与样条结点、级数项有关的广义位移参数，δ={δ1T，δ2T，⋯，δRT}，其中δm={δ0mT,δ1mT,⋯,δNmT}，δim={ai，bi,ci}mT(m=1,2,⋯,R;i=0,1,⋯,N)；N为形函数矩阵，N=[N1,N2,⋯,NR]，Nm= [N0m,N1m,⋯,NNm]，Nim=diag{Xm(y)φi(x)，Ym(y)φi(x)，Hm(y)φi(x)}，取φi(x)=103φ3xh-i-43φ3xh-i+12-43φ3xh-i-12+16φ3xh-i+1+16φ3xh-i-1  ,其中φ3(x)为三次B样条函数φ3(x)=16(x+2)3    (x∈(-2,-1)),(x+2)3-4(x+1)3    (x∈[-1,0)),(2-x)3-4(1-x)3    (x∈[0,1)),(2-x)3    (x∈[1,2)),0    (|x|≥2).QR法的支配方程为Kδ=F，(2)式中K=∑eK¯e，F=∑eF¯e，K¯e与F¯e分别为QR法的单元刚度矩阵与单元荷载向量，K¯e=NeTk¯eNe；F¯e=NeTf¯e，其中，k¯e为用整体坐标描述的单元刚度矩阵，f¯e为用整体坐标描述的单元等效荷载向量，Ne为单元的形函数矩阵，一个单元包含两个结点Ne=[Ni,Nj]T，Ni与Nj为单元两端的形函数矩阵，将单元结点坐标(xi,yi)和(xj,yj)代入形函数矩阵表达式中即可得到．求解式(2)可以得到QR法的结构广义位移，代入Ve=Neδ=[Ni,Nj]Tδ可以获得单元两端结点的位移．在QR法中，整体刚度矩阵K与整体荷载向量F可由QR变换后的单元刚度矩阵与单元荷载向量直接叠加得到，无须装配，这是QR法对比传统有限元法的重要优势．2.2　改进的摄动随机有限元法在传统的摄动随机有限元法中，位移和刚度矩阵在均值处的二阶泰勒展开式为U=U0+∑i=1qUiIαi+∑i=1q∑j=iqUijIIαiαj+O(αi3)；K=K0+∑i=1qKiIαi+∑i=1q∑j=iqKijIIαiαj+O(αi3)，式中：U0和K0分别为位移向量和刚度矩阵的均值；αi为随机变量；UiI，UijII，KiI和KijII分别为位移与刚度矩阵对随机变量的一阶和二阶偏导，通过有限元法支配方程可以由KiI和KijII得到UiI和UijII，从而得到位移的均值与方差；O( )为截断误差．在一些复杂的随机问题中，获取KiI和KijII较为困难，制约了该方法的应用．改进的摄动随机有限元法，根据随机变量之间的相关性可以分为三种情况，由于K-L法得到的随机变量是一组不相关的标准正态分布随机变量，因此选择随机变量不相关且具有对称联合概率密度函数的情况．在MPSFEM中，结构的随机向量α被替换为确定性向量bi，即bi=0,0,⋯,0︸i-1,ρiiiiσi,0,0,⋯,0︸q-i，式中：ρiiii为随机变量αi，αi，αi和αi均值与其标准差之比；q为随机变量总数，因此替换要总共进行q次；σi为第i个随机变量的标准差．当随机变量服从高斯分布时，ρiiii=3，则bi可改写为bi=0,0,⋯,0︸i-1,3,0,0,⋯,0︸q-i．在MPSFEM中，位移在均值处的四阶泰勒展开式为：U(bi)=U0+UiIρiiii12σi+UiiIIρiiiiσi2+UiiiIIIρiiii32σi3+UiiiiIVρiiii2σi4+O(||σ||∞5) ; (3)U(-bi)=U0-UiIρiiii12σi+UiiIIρiiiiσi2-UiiiIIIρiiii32σi3+UiiiiIVρiiii2σi4+O(||σ||∞5)  , (4)式中UiiiIII和UiiiiIV分别为位移对随机变量αi，αi，αi与αi，αi，αi，αi的三阶和四阶偏导．将式(3)与式(4)相加，再用式(4)减去式(3)可得zi2ρiiii+O(||σ||∞6)=UiiIIσi2+UiiiiIVρiiiiσi4；wi2ρiiii+O(||σ||∞5)=UiIσi+UiiiIIIρiiiiσi3，式中：zi=U(bi)+U(-bi)-2U0；wi=U(bi)-U(-bi)，U(±bi)可以从下式中求得K(±bi)U(±bi)=F，其中，F为荷载向量，K(±bi)和U(±bi)为随机向量α被替换为确定性向量bi与-bi后得到的包含bi与-bi的刚度矩阵K(bi)，K(-bi)和位移向量U(bi)，U(-bi)．将zi和wi代入MPSFEM的位移均值与方差矩阵表达式中，可以得到：E(U)=U0+12∑i=1qziρiiii+O(||σ||∞4)；Cov(U,U)=14∑i=1qwiwiTρiiii+ziziTρiiii2(ρiiii-1)+O(||σ||∞4).显然，MPSFEM通过随机向量的替换巧妙地跳过了传统摄动随机有限元法中刚度矩阵对随机变量求偏导的步骤，节省了大量的字符运算时间，编程效率也高，是一种高效的随机分析手段．3 改进摄动随机QR法MPSFEM是基于有限元法进行开发的，其支配方程为KU=F，而QR法同样是基于有限元法开发的，其支配方程(式(2))与有限元法支配方程形式相同，且QR法与有限元法实际上是通用的，因此将MPSFEM中的有限元法计算格式替换为QR法的计算格式是可行的．将有限元法的位移向量U(±bi)替换为QR法的广义位移向量δ(±bi)进行泰勒展开，经过同样的公式推导可以得到：E(δ)=δ0+12∑i=1qziρiiii+O(||σ||∞4)；Cov(δ,δ)=14∑i=1qwiwiTρiiii+ziziTρiiii2ρiiii-1+O(||σ||∞4),式中δ(±bi)可由下式求出K(±bi)δ(±bi)=F，(5)K(±bi)与F则为QR法的整体刚度矩阵与整体荷载向量．这样，K-L法、QR法与MPSFEM就结合到了一起，形成基于K-L级数的改进摄动随机QR法．4 算例如图2所示的双跨四层平面框架结构，柱截面面积Ac=0.012 m2，惯性矩Ic=2.05×10-4 m4，梁截面面积为Ab=0.011 m2，惯性矩Ib=1.7×10-4 m4，弹性模型为2.1×1011 Pa，框架左侧每层都有作用于柱顶的集中荷载P=20 kN，每层梁上作用均布荷载Q=100 kN/m．10.13245/j.hust.210704.F002图2双跨四层平面框架结构及其单元划分分别使用有限元法与QR法对结构进行计算，取结点7，8，9，10侧向位移及这些点在⑤，⑥，⑦，⑧单元中的弯矩列于表1中进行对比．在本算例中，竖直方向的多项式级数与层数相同，即R=4，水平方向上以步长h=4 m将框架离散为2份，即N=2．由表1可知：QR法精度很高，计算结果与有限元法几乎无异，因此使用QR法替换有限元法是可行的．10.13245/j.hust.210704.T001表1各结点位移与弯矩结点编号有限元法QR法位移/mm弯矩/(kN•m)位移/mm弯矩/(kN•m)72.662 137.844 02.662 137.844 086.067 943.509 96.067 943.509 998.569 031.233 08.568 931.233 0109.978 618.634 79.978 618.634 5将图2平面框架结构的弹性模量设置为随机场，均值为2.1×1011 Pa，协方差函数为高斯型，各梁和柱随机场相互独立，梁随机场相关长度为2 m，柱随机场相关长度为1 m，计算不同变异系数下结点15的水平位移均值与标准差．在本算例中，小波阶数为9阶，根据比能因子，K-L级数截取3项．除提出的KLSMPSQRM外，基于上述思想，分别将蒙特卡罗有限元法和传统的摄动随机有限元法与K-L法和QR法进行了结合，简称为MCM[15]和PSM，并使用Matlab编写了相应的程序．为保证精度，蒙特卡罗法采用拉丁超立方抽样技术抓取了5×104个样本，以此作为算例的精确解，来检验其他算法的精度．对于这类含有多个相互独立的随机场的问题，采用的离散方法是独立离散，然后将离散得到的随机变量组合成一个随机向量，依次代入计算即可．本算例中包含20个相互独立的随机场，根据截取项数得到60个随机变量，单元可分为ξ1(θ)={ξ1(θ),ξ2(θ),ξ3(θ)};ξ2(θ)={ξ4(θ),ξ5(θ),ξ6(θ)};⋮ξ20(θ)={ξ58(θ),ξ59(θ),ξ60(θ)},则这20个随机向量可以组合为一个随机向量ξ(θ)={ξ1(θ),ξ2(θ),⋯,ξ20(θ)}={ξ1(θ),ξ2(θ),⋯,ξ60(θ)},之后根据随机变量编号代入本文算法中计算即可．为方便对比结果，各随机场变异系数采取统一数值．结点15水平位移均值和标准差如表2和表3所示．可知：当变异系数为0.05~0.15时，3种方法精度几乎没有差别，而当变异系数为0.20和0.25时，KLSMPSQRM的精度要明显高于PSM，这是因为PSM通常只有一至二阶泰勒展开，而KLSMPSQRM通过高阶的泰勒展开与随机向量替换使算法精度提升，拓展了应用范围．从效率上看，蒙特卡罗法须计算5×104次，而KLSMPSQRM只须计算一次即可得到最终结果，节省了大量计算时间．10.13245/j.hust.210704.T002表2结点15水平位移均值变异系数MCMKLSMPSQRMPSM0.059.869 79.869 79.869 70.109.892 39.892 29.892 10.159.930 79.929 89.929 20.209.986 09.982 99.981 30.2510.059 310.051 910.048 2mm10.13245/j.hust.210704.T003表3结点15水平位移标准差变异系数MCMKLSMPSQRMPSM0.051.052 21.050 91.048 60.102.124 62.113 92.097 20.153.245 33.201 73.145 70.204.461 14.328 94.194 30.255.699 35.512 75.242 910-4 m5 结语使用QR法代替了改进的摄动随机有限元法中有限元计算格式，结合Karhunen-Loève级数展开法，提出了基于K-L级数的改进摄动随机QR法，并与其他方法进行了对比计算．算例分析得到：在精度方面，KLSMPSQRM可以在变异系数0~0.25范围内保持良好的精度，适用范围比PSM更广；在效率方面，KLSMPSQRM的计算速度要远快于MCM，节省了大量的计算时间．本文方法也存在一定的局限性，框架结构中每根杆件都视为一个独立的随机场进行离散，这样处理会产生大量的随机变量，降低了计算效率，因此可在框架结构的空间变异性上做进一步的研究，寻找一种快速便捷的随机场布置与离散方式．