对于大多数实际工程问题，都会存在不可避免的制造误差、工艺不稳定性和结构损伤等，表示结构可变性的不确定性无处不在[1-2]，这些不确定性将会影响结构的最终性能．在不确定性分析中，概率方法和非概率方法常用于量化不确定性对结构力学行为的影响．然而，概率方法在构造精确的概率分布函数时，需要大量的信息．但对于许多工程问题来说，收集关于不确定性的足够信息往往非常困难且成本高昂．相比之下，非概率方法中的区间分析非常适合处理边界明确但缺少充分的概率分布函数信息的不确定问题[3]．目前，有许多改进的区间分析方法[4]，比较常用的可分为顶点法[5-6]、区间摄动法[7-9]和配置法[10-11]．对于求解线性区间方程的精确边界，顶点法是一种较好的选择[5]．但该方法的应用仅限于单调情况，且计算成本会随着区间变量的增加呈指数增长．结构越复杂，所含的区间参数越多．为了提高复杂结构不确定性问题的计算效率，文献[7]基于一阶泰勒展开和参数摄动，求解了具有区间参数的结构静力位移问题．上述区间摄动法只考虑了一阶摄动近似而忽略了高阶项．面对非线性问题和不确定性较大的情况，该方法往往会给出粗略的结果．为了解决这一问题，子区间法被引入到区间摄动法中[12-13]．但随着区间参数和被划分子区间数的增加，其计算量将大大增加．当问题是非线性且区间参数的数量较少时，配置法被引入到区间分析中．配置法通过在不确定参数空间中设置一定的样本点来构造合理的替代模型，从而得到结构响应的边界[14]．其中，Chebyshev多项式是较为常用的一种替代模型．文献[10]基于Chebyshev多项式，提出了一种新的动态非线性系统的不确定性分析方法．文献[11]为了求解非概率区间过程模型下结构的响应边界，提出了一种基于Chebyshev替代模型的数值方法．但在这些问题中，不确定参数的数量相对较少，随着不确定参数数量的增加，这些配置方法的计算复杂度将呈指数级增加，陷入维度灾难．为了解决配置法在构造替代模型时计算成本昂贵的问题，本研究提出了一种基于降维算法和Chebyshev多项式(dimension reduction method and Chebyshev polynomial，DC)的区间不确定性分析方法．首先基于降维算法将高维响应函数分解为多个一维函数的和，再利用Chebyshev多项式来拟合一维响应函数，由此建立了一个关于区间参数的实际结构响应的替代模型．利用辅助方法来计算代理模型的极值，就可得到结构响应的上下界．1 问题的提出在有限元法中，静态位移和外部载荷的平衡方程可表示为Ku=F，(1)式中：K∈Rl×l为结构的刚度矩阵，l为结构的自由度；u∈Rl×1为结构的位移向量；F∈Rl×1为外载荷向量．令αI∈Rn为区间参数，上标I表示区间．那么，考虑区间参数的外部载荷和静位移的平衡方程可改写为K(αI)u(αI)=F(αI)，(2)式中αI=[α1I,α2I,…,αnI]且αiI=[α̲i,α¯i]．区间参数αiI的不确定水平可表示为λi=α¯i-α̲iα̲i+α¯i×100%，(3)式(2)的解可表示为u={u(α)|K(α)u(α)=F(α),α∈αI}．(4)区间有限元法的目的是确定最小的超立方体，来包围式(4)所定义的几何形状．超立方体可以用其边界来表示，其边界为uI=[u̲,u¯]，(5)式中：u̲=min{u=K-1(α)F(α)|α∈αI};u¯=max{u=K-1(α)F(α)α∈αI}, (6)求解式(5)和(6)的最直接的方法是蒙特卡罗模拟(Monte Carlo simulation，MCS)．它会在不确定域内生成样本点，并通过有限元法来计算这些样本点的结构静态响应．所得结果的极值可以视为响应区间的界限．由于在应用MCS时会进行大量的有限元分析，因此MCS常被用作参考解．2 Chebyshev替代模型利用Chebyshev多项式来近似式(4)中含不确定参数的位移响应u(α)．首先，将不确定参数α中的任意元素αi (i=1,2,…,n)通过变量转换，线性的变化为xi∈[-1,1]．变量转换表达式如下xi=2αi-(α̲i+α¯i)α¯i-α̲i∈[-1,1]. (7)然后，令Chebyshev多项式函数uC(x)在不确定域内对u(α)进行近似，其中x=[x1,x2,…,xn]，可得Chebyshev替代模型(Chebyshev surrogate model，CM)为u(α)≈uC(x)=∑i1=0k∑i2=0k⋯∑in=0k(1/2)pui1i2…inCi1i2…in(x), (8)式中：p为下标i1i2⋯ik中出现零的个数；ui1i2…in为n维Chebyshev多项式系数；Ci1i2…in(x)为n维Chebyshev多项式．若用Ck来表示阶数为k的Chebyshev多项式，则可以得到Ck(xi)=cos kθi,θi=arccos(xi)∈[0,π]．(9)对于n维问题，Ci1i2…in(xi)是每个一维多项式Ck(xi)的张量积，那么Ci1i2…in(x)=Ci1(x1)Ci2(x2)…Cin(xn)．(10)式(8)中ui1i2…in可使用Gauss-Chebyshev插值积分公式计算ui1i2…in=(2/π)n∙∫-11∫-11⋯∫-11u(x)Ci1i2…in(x)1-x121-x22⋯1-xn2dx1dx2…dxn≈(2/m)n∑j1=1m∑j2=1m⋯∑jn=1mu(xj1,xj2,…,xjn)∙Ci1i2…in(xj1,xj2,…,xjn), (11)式中：m为插值点的个数(本研究m=k+1)；xjn为插值点；u(xj1,xj2,…,xjn)为插值点处的位移值．每一维的插值点都是k+1阶Chebyshev多项式的零点．插值点由下式确定xjn=cos2jn-1k+1π2    (jn=1,2,…,k+1) ．(12)可以发现：对于一维问题，所需插值点的数量为k+1个．若为n维问题，则需要(k+1)n个插值点．显然，插值点的数量随着不确定参数的增加呈指数增长，而每一个插值点都须进行一次有限元计算．因此对于含很多不确定参数的位移响应问题，若直接使用Chebyshev替代模型则会产生巨大的计算成本．3 DC近似模型的构建为了减少利用配置法求解响应边界的计算量，首先利用降维算法将n维响应函数分解为n个一维函数的和，再基于Chebyshev多项式拟合一维响应．3.1　降维算法降维算法的实质是将多元函数进行维数分解，得到输入变量随维数增加的简单分量函数的有限和．u(α)的降维表达式为[15]u(α)=u0+∑i=1nui(αi)+∑1≤i1i2≤nui1i2(αi1,αi2)+⋯+∑1≤i1⋯is≤nui1i2⋯is(αi1,αi2,⋯,αis)+∑1≤i1⋯is+1≤nui1i2⋯is+1(αi1,αi2,⋯,αis+1)+⋯+u12⋯n(α1,α2,⋯,αn), (13)式中：u0为常数项；ui1i2⋯is(αi1,αi2,⋯,αis)为量化s个输入变量αi1,αi2,⋯,αis对u的协同影响的s元分量函数；u12⋯n为u的残余项．忽略式(13)的残余项可以得到u(α)≈us(α)=u0+∑i=1nui(αi)︸u1(α)+∑1≤i1i2≤nui1i2(αi1,αi2)+⋯+∑1≤i1⋯is≤nui1i2⋯is(αi1,αi2,⋯,αis), (14)式中us(α)为u(α)的s维近似函数．其中，单变量近似函数u1(α)利用了式(14)中有限分解的前两项．同样地，通过适当地选择s值，可以导出二元和其他更高阶的近似函数．当s无限趋近于n时，us(α)收敛于u(α)．选取不确定参数空间α中的任一点集c=(c1,c2,…,ci,…,cn)，将其作为参考点．那么，式(14)中的ui(αi)可表示为ui(αi)=u(c1,c2,⋯,ci-1,αi,ci+1,⋯,cn)-u0. (15)只保留式(14)的前两项，则u(α)的降维表达为u(α)=∑i=1nu(c1,c2,⋯,ci-1,αi,ci+1,⋯,cn)-n-1u0. (16)c一般选定为不确定参数的中点αC[16]，所以u(α)的单变量降维表达式为u(α)=∑i=1nu(α1C,α2C,⋯,αi-1C,αi,αi+1C,⋯,αnC)-(n-1)u0=∑i=1nu(αi)-(n-1)u0．(17)由式(17)可以看出，u(α)中求和的每一项都是仅含一个不确定参数的函数．3.2　DC近似模型对于式(17)中n个一维函数，本研究利用式(7)中的Chebyshev近似模型进行拟合．在n个一维不确定性域内，将式(17)中的每个单变量函数近似为u(αi)≈uC(xi)=∑j=0k(1/2)pui,jCj(xi)，(18)式中：uC(xi)为只含xi的一维Chebyshev多项式，参数αi和xi的转换关系见式(7)；ui,j和Cj(xi)分别为只含第i个不确定参数的系数和Chebyshev多项式，根据式(11)，ui,j可表示为ui,j=2m∑l=1mu(xi,l)Cj(xi,l)，(19)其中xi,l为一维函数uC(xi)的第l个插值点．然后，将式(9)和(19)代入到式(18)中，可得每个一维函数的Chebyshev多项式表达．再将式(18)代入到式(17)中，得到DC近似模型u(α)≈uDC(x)=∑i=1nuC(xi)-n-1u0=∑i=1n∑j=0k(1/2)pui,jCj(xi)-n-1u0．(20)利用式(20)拟合原始响应函数的优势在于：首先，利用单变量降维表达式将原始n维响应的求解问题转换为一维．然后，相比于式(8)需要(k+1)n个样本点，式(20)的样本数为n×k+1(k为偶数)或n×(k+1)+1(k为奇数)个，避免了不确定参数较多时发生维度灾难．3.3　DC近似模型的验证为了证明本研究所提模型的有效性，采用的测试函数[17]为f=x1sin(5πx1)-x2sin(5πx2+π)+1，(21)式中x1和x2为区间参数，取值范围均为[-0.2，0.2]．引入指标最大误差E1来描述DC区间模型的精度，E1=max(|f-fC|)/max(|f |)，(22)式中fC为DC近似模型的值．测试函数在不确定性域内的变化曲线如图1所示．10.13245/j.hust.210803.F001图1测试函数图2给出了DC和CM的最大误差(1g E1)随多项式阶数(N)的变化情况．应该注意，图2中所示的误差是对数误差，因此y轴的值越小表示误差越小．从图2中可以看出，DC和CM均有着较高的拟合精度．10.13245/j.hust.210803.F002图2DC和CM相对测试函数的误差图3和图4分别给出了六阶DC和CM近似模型及其相应的样本分布，图中黑色虚线对应的黑色点为样本点．可以看出：在拟合精度都非常高的前提10.13245/j.hust.210803.F003图3六阶DC近似模型及其相应的样本分布10.13245/j.hust.210803.F004图4六阶CM近似模型及其相应的样本分布下，CM的样本数为49个，而DC仅为13个．可见，DC有着更少的样本点，且对于高维问题计算成本的优势更加明显．4 基于DC近似模型的区间分析在原始位移响应u(α)中考虑区间不确定性．由于原始函数是通过近似模型uDC(x)逼近，因此求解u(αI)的区间问题就是计算含区间参数uDC(xI)的边界问题．DC区间模型为uDC(xI)=∑i=1n∑j=0k(1/2)pui,jCj(xiI)-(n-1)u0．(23)在得到DC区间模型后，下一步是计算其边界．首先，将式(23)等号右边的常数项提出，DC区间模型可改写为uDC(xI)=(∑i=1nui,0/2-(n-1)u0)+∑i=1n∑j=1kui,jCj(xiI)=C0+∑i=1n∑j=1kC˜j|ui,j|, (24)式中：C0=∑i=1nui,0/2-(n-1)u0；C˜j为Cj(xiI)在xiI=[-1,1]内的界限．由式(9)可知，当xiI=[-1,1]时，C˜j∈[-1,1]．因此，基于区间运算，DC区间模型的下界为u̲DC(xI)=C0-∑i=1n∑j=1k|ui,j|．(25)同时，DC区间模型的上界为u¯DC(xI)=C0+∑i=1n∑j=1k|ui,j|．(26)由式(25)和(26)可以看出：在求解过程中，主要的任务是根据式(19)计算Chebyshev多项式系数ui,j．5 数值算例图5所示的四杆桁架由4根杆件组成，长度L为500 mm．其中，节点2承受水平作用力P1，节点3承受垂直作用力P2．由于制造和测量误差，因此将4根杆件的截面积、弹性模量和外载荷视为区间参数．P1和P2的区间范围分别为[19，21] kN和[24，26] kN，其余8个区间参数的信息见表1．10.13245/j.hust.210803.F005图5四杆桁架10.13245/j.hust.210803.T001表1四杆桁架中区间参数的信息参数杆件号(1)(2)(3)(4)截面积/mm2[99，101][79，81][59，61][39，41]弹性模量/MPa[205，215][265，275][315，325][385，395]将桁架节点2的水平位移u2和节点3的竖直位移v3作为目标函数．利用MCS，Chebyshev区间分析方法和本研究方法分别计算u2和v3的边界值．其中，MCS的计算方法见文献[18]．表2和表3分别给出了小不确定水平(λ=5%)和大不确定水平(λ=10%)下的位移边界值和误差．为方便对照，DC和CM的Chebyshev多项式的阶数均为2．在MCS中，样本数量为1×105，其结果作为参考值．由表中的计算结果可以得出：10.13245/j.hust.210803.T002表2小不确定水平下节点位移的边界值和误差(λ=5%)方法边界u2v3样本数取值误差/%取值误差/%MCS下界0.328 3-0.338 91×105上界0.443 2-0.251 0CM下界0.320 82.28-0.338 60.0959 049上界0.443 00.04-0.244 82.47DC下界0.324 71.09-0.335 60.9721上界0.439 10.92-0.248 11.1610.13245/j.hust.210803.T003表3大不确定水平下节点位移的边界值和误差(λ=10%)方法边界u2v3样本数取值误差/%取值误差/%MCS下界0.283 4-0.395 61×105上界0.517 3-0.216 7CM下界0.253 710.45-0.394 90.1859 049上界0.515 90.28-0.192 011.40DC下界0.269 94.74-0.381 13.6721上界0.499 73.42-0.205 85.03a. 在小不确定水平下，CM的最大边界误差为2.47%，但是需要59 049个样本点．相比于CM，DC的结果与MCS更接近，最大边界误差为1.16%．同时，DC仅需21个样本点，即可获得可接受的准确度．b. 当不确定水平的较大时，CM的最大边界误差为11.40%，而DC的最大边界误差为5.03%．可见，随着不确定水平的增大，CM和DC的误差也在增大．但在样本数远小于CM的前提下，DC的误差增长幅度还较小．在传统的区间分析中，总是通过减少不确定性参数的数量来简化计算模型．例如，只考虑对结构影响较大参数的不确定性，忽略影响较小的参数．如果忽略弹性模量的不确定性[19-20]，那么本例中的区间参数从10个减少到6个．图6为利用本研究方法计算的含6个和10个区间参数的位移边界．10.13245/j.hust.210803.F006图6含6个和10个区间参数的位移边界从图6中可以看出：当不确定参数减小到6个时，所预测位移区间值的下边界变大，上边界变小．同时，随着不确定水平的增加，所预测上边界和下边界的误差也越来越大．这意味着当不确定性参数不足以描述结构中的不确定性时，无法预测最坏的情况．随着区间参数范围的增大，上述这种现象将更加明显．6 结论为了求解具有大规模不确定参数的结构响应问题，提出了一种基于DC区间模型的区间分析方法．与现有的配置方法不同，该方法在保证精度的前提下，可有效地减少计算成本．a. 利用降维算法和Chebyshev多项式展开所得到的区间模型，将n维响应函数的边界求解问题转换为求解n个一维函数的边界问题，且有较好的拟合精度．b. 数值算例结果表明：对于含大规模且不确定水平较大参数的结构响应，本研究方法相比于CM有着更少的计算样本和更高的计算精度．c. 随着结构的复杂化，不确定参数的数量急剧增加．结果表明：采用较少不确定参数的简化模型不能保证预测到结构的最坏情况，应充分考虑不确定参数．