朱国树, 邵亚会

(合肥工业大学 土木与水利工程学院,安徽 合肥 230009)

随着桥梁设计施工理论的进步和新型材料的应用,现代桥梁的跨度不断增大,结构刚度也随之降低,使得桥梁风振问题日益突出。在桥梁的各种风振形式中,颤振是一种危险的自激发散风致振动,一旦发生就可能会造成整座桥梁的坍塌,引起灾难性的后果。近几十年来,研究人员将风洞试验、计算流体动力学和理论分析等方法应用到颤振研究中,取得了一系列成果[1-5],使得颤振分析理论,尤其是频域分析理论日益成熟,发展出了多模态和全模态2类方法,被广泛用于桥梁三维颤振分析中。实际工程中,由于设计、施工和风洞试验过程中不可控因素的影响,桥梁的自身特性(如刚度、质量和气动参数等)不可避免地存在不确定性,而且桥址处来流风速也是随机的,从而导致桥梁颤振失稳成为一个随机事件。因此,采用概率方法合理地考虑这些不确定性因素的影响,进而求解桥梁颤振失效概率尤为重要。

桥梁抗风设计中,颤振临界风速代表着结构的抗力,求解过程涉及大量的试验和数值计算,是失稳分析中最为关键的部分。早期的颤振概率性评价方法通常是基于经验或少量的试验来确定颤振临界风速的统计参数[6-7],并使用显式的功能函数,这类方法因过于简略而难以准确描述结构的失效模式。为了解决这一问题,基于有限元和颤振理论建立结构极限状态方程的颤振可靠度分析方法陆续被提出,例如,基于多模态分析理论,文献[8-9]分别提出颤振可靠度分析的随机有限元法和响应面法。影响颤振临界风速的诸多参数中,颤振导数被认为是最为重要的因素之一,其微小扰动即有可能对结果产生显著影响[10];文献[11]认为不同折减风速下的颤振导数均应被视为单独的随机变量,这一观点的合理性在后续的试验和理论研究中得到了验证,因此被广泛采用[12]。然而,这也不可避免地导致变量的数目大大增加,使得颤振失稳成为高维可靠度问题。文献[13]在研究颤振约束下的悬索桥可靠性优化设计时考虑了颤振导数之间相关性的影响,发现忽略颤振导数的相关性通常会导致相对保守的结果;文献[14]通过开展大量风洞试验,研究半开口分离双箱截面颤振导数的统计分布规律,结果表明颤振导数的相关性在分析中是不可忽略的。

综上所述,桥梁颤振失稳作为一类具有高维非线性和非显性功能函数的可靠度问题,使用传统的一次阶法、二次矩法和响应面法求解的精度难以让人满意[15],直接抽样的蒙特卡洛(Monte Carlo,MC)方法虽然有很高的可信度,但庞大的计算量导致其不再适用。为了解决这些问题,提高颤振可靠度分析的计算效率和精度,本文提出一种基于径向基函数(radial basis function,RBF)神经网络代理模型和自适应随机采样策略的桥梁颤振可靠度分析方法。首先,基于全模态理论建立参数化颤振临界风速分析有限元模型,将试验设计生成的样本代入有限元模型求解得到训练集,构建逼近特征量与设计参数间复杂映射的初代RBF代理模型;然后,通过定义的自适应采样策略,在迭代过程中不断从MC模拟生成的大量样本中筛选出失效面附近的样本点添加至训练集,重新训练RBF模型直至满足收敛要求,并使用该RBF模型识别处于失效域的MC样本点,从而直接计算出失效概率;最后通过2个数值算例对该方法进行了验证,结果表明,本文方法能够在保证计算精度的前提下,大幅度降低计算复杂度,有效提高计算效率。

1 全模态颤振分析理论

全模态方法属于桥梁颤振频域研究方法,它将风-桥梁作为耦合系统,直接在物理坐标中分析耦合系统的复模态特性来判断系统的气动稳定性,避免了多模态方法叠加模态选择和精度难以验证的问题[16-17]。桥梁断面气动自激力如图1所示。

图1 桥梁断面气动自激力

按照Scanlan的理论,主梁受到的气动升力Lse、气动扭矩Mse和阻力Dse与其竖向位移h、扭转位移α、侧向位移p及其一阶导数成线性关系,使用18个颤振导数表示的单位长度主梁气动自激力模型如下：

(1)

使用有限元方法进行颤振分析时,主梁被离散为若干个梁单元。每个梁单元所受的气动自激力可由(1)式求得,使用简单的堆积原则将沿主梁分布的自激力转换成等效节点荷载,即

(2)

因此,均匀流中风-桥耦合系统(简称系统)的振动方程可以表示为:

(3)

将(3)式移项,可得:

(4)

当给定风速和振动频率后,使用复模态分析方法[2]求解(4)式可得系统的n对共轭特征值λj=μj±iωj,其中n为系统的自由度数。如果所有的特征值实部μ都小于0,那么系统就是动力稳定的;若存在任意一对特征值的实部大于0,则系统是发散失稳的。因此,全模态方法是一种双参数搜索方法,即搜索使系统失稳的最小风速(颤振临界风速)及其对应的颤振临界频率ω。

基于上述全模态分析理论,本文采用MATLAB和ANSYS APDL混合编程技术建立参数化颤振临界风速分析有限元模型。其中,MATLAB负责生成样本点和后台调用ANSYS,ANSYS则用于实现颤振临界风速搜索过程[16]。ANSYS的自定义单元Matrix 27能够很好地模拟气动自激力,同桥梁有限元模型相结合,即可建立颤振分析系统。由于使用了混合编程技术,MATLAB和ANSYS能够实现数据交换,并按照事先设定的逻辑交替运行,进而实现整个分析过程的参数化和自动化。

2 颤振失效概率分析理论基础

2.1 失效概率

结构的失效概率被定义为:

(5)

其中:f(X)为影响结构性能的多维随机变量X的联合概率密度函数;G(X)为结构功能函数;极限状态面G(X)=0将X的设计空间分为失效域G(X)≤0和安全域G(X)>0。

(6)

(7)

2.2 极限状态函数

参考已有的颤振可靠度分析模型[6],本文将桥梁颤振失稳发生的极限状态定义为桥址处极值风速达到颤振临界风速,相对应的结构功能函数为:

G(X)=Uf(x1,x2,…,xn)-Us

(8)

其中:Uf为桥梁颤振临界风速;Us为桥面高度年平均最大风速,即Uf和Us分别代表结构的抗力和荷载效应;xi(i=1,2,…,n)为影响Uf的各种随机变量;Uf(x1,x2,…,xn)表示求解桥梁颤振临界风速Uf的过程。

2.3 随机变量统计分布

年平均最大风速Us的概率分布可以由桥址处附近的气象站资料统计分析得出,也可在桥面处安装风速仪采集数据,从而得到更加精确的统计参数。一般来说,从历史风速数据中提取的极值风速样本的累积分布函数能够很好地被极值分布所描述。我国大多数气象站资料统计分析表明,Us的最优分布为极值Ⅰ型(Gumbel)分布。

颤振临界风速Uf是结构刚度、质量、阻尼和桥梁断面气动特性的隐式函数,其不确定性取决于上述因素的随机分布。基于此,本文将主梁的弹性模量、密度、模态阻尼比和颤振导数均视为随机变量,分别用于表征结构刚度、质量、阻尼和气动特性的不确定性。根据一般经验,弹性模量和阻尼比服从对数正态分布,密度服从正态分布。考虑到现代桥梁的跨度较大,可将全桥主梁划分为若干段,每段主梁的弹性模量和密度均由独立抽样确定。

颤振导数一般由节段模型风洞试验确定,在试验过程中逐级增加风速,并使用模态识别方法对采集的振动信号进行处理,识别出系统阻尼和振动频率[18-19],进而确定当前折减风速下的颤振导数。然而,由于风洞试验模拟的风环境特性(湍流强度、风攻角和风速不均匀度等)的不确定性、试验初始条件和测试误差的随机性、静风引起的附加质量和附加阻尼等因素的共同作用,使结构颤振导数表现出复杂的变异性。参考已有的研究[11-14],本文将每个折减风速下的颤振导数均视为服从正态分布的随机变量,同时也可考虑颤振导数之间的相关性。原则上,不同类型的桥梁断面应通过专门的重复试验来确定颤振导数的统计参数。当缺乏资料时,也可按下述方法来近似地考虑颤振导数的不确定性[13]:最大折减风速下颤振导数变异系数取较大的值,折减风速为0时变异系数取为0,其他折减风速下的变异系数由线性插值来确定,如图2所示。

图2 颤振导数随机分布

3 颤振失效概率高效算法

MC方法的缺点在于计算量大,尤其是对于单个样本计算时间较长的情况,该方法将不再适用。针对这一问题,通常将方差缩减技术和代理模型技术用于减少MC方法的计算成本[15]。方差缩减技术通过减少失效概率估计值的方差来加速收敛,但是仍然需要数以万计的样本。代理模型技术则是使用人工神经网络、克里金模型和支持向量机等计算复杂度低的数学模型来逼近随机变量与响应之间的映射关系,从而缩减单个样本的计算时间。在众多的代理模型中,RBF神经网络以学习效率高、非线性映射能力强和无局部最小问题而闻名,因此本文选择RBF神经网络构建颤振失效预测模型。

3.1 RBF神经网络

RBF神经网络由输入层、隐含层和输出层构成,具有以任意精度逼近任意连续函数能力的前馈神经网络,其结构如图3所示,输入层接收的外界信息X通过隐含层的径向基函数映射到隐空间H,并最终通过线性变换得到输出Y。径向基函数通常是一个非线性径向对称函数,能够对输入信号在局部产生响应,且当输入信息靠近径向基函数的中央范围时会产生较大的输出,从而使RBF神经网络具备局部逼近能力。最常用径向基函数为高斯基函数,如图3所示。经试验设计确定训练集后,可通过自组织选取中心法和正交最小二乘法等学习策略来确定模型参数。

图3 RBF神经网络结构

3.2 自适应随机采样方法

使用代理模型求解失效概率之前,需要预先获取样本点用以训练模型,使其满足一定的精度要求。因此,训练样本的试验设计方法成为决定代理模型预测精度的主要因素。按照抽样流程,试验设计方法可分为直接抽样法和序贯抽样法[20]。直接抽样法是通过拉丁超立方设计、均匀设计和正交设计等基于空间填充思想的采样方法一次确定全部训练样本。因样本点遍布设计空间,该方法通常能够较好地保证模型的全局精度。然而随着问题维度上升,直接抽样法所需的样本点数目会呈指数上升,因而难以在复杂工程中应用。序贯抽样法也称自适应采样方法[21],其基本思想是首先使用少量的样本构造初始代理模型,然后利用已有样本和当前代理模型的信息,基于一定的自适应采样策略,在迭代过程中选取对模型预测精度有较大影响的高价值样本点添加至训练集,从而实现用较少的样本点构造高精度的代理模型。

由(5)式可知,可靠度问题中,代理模型的预测值大小并不重要,其正负才是关注的重点。换而言之,代理模型只需要在极限状态面附近,尤其是最可能失效点附近的极限状态曲面有足够精确的逼近即可,因此可靠度分析中的高价值样本点应具有以下特征[15,22]:

(1) 接近极限状态面,即|G(X)|≈0。

(2) 失效模式发生可能性较大,即样本点概率密度函数值f(X)较大。

鉴于此,本文提出一种基于加权抽样和动态距离约束的自适应策略。首先,使用Rosenblatt变换,得到MC方法生成的标准正态空间样本Ui(i=1,2,…,nMC)及其对应的非标准正态变量样本Xi(i=1,2,…,nMC),并使用当前代理模型预测相应的功能函数值;然后,选取前m个符合约束条件的最靠近极限状态面的样本点作为候选点,约束条件为:

(9)

其中:dU*为样本U*与已有ndoe个训练样本之间的欧式距离最小值;D为距离阈值,取已有训练样本之间的欧式距离最小值中的最大值[23]的1/2,即

D=max{Di|Di=0.5min{‖Xi-Xj‖},

i≠j;i,j=1,2,…,ndoe}

(10)

最后进行加权抽样,在(0.5,1)区间生成m个随机数ri,并计算各个候选点重要性系数:

ci=lnri×‖Ui‖

(11)

选择m个候选样本中ci值最大的样本作为最优样本添加至训练集。

距离约束能够保证选取的样本点均匀分布在极限状态面左右,避免出现样本点“聚集”的现象,从而减少冗余信息[21]。加权采样则是为了更多地选择概率密度值较大的样本点,同时也能够兼顾设计空间边缘的样本点。

3.3 计算流程

本文所提方法的具体步骤为:

(1) 建立颤振分析有限元模型。分析影响桥梁颤振失稳的随机变量X=(x1,x2,…,xn),基于全模态分析理论,建立参数化桥梁颤振临界风速有限元模型。给定随机变量取值后,有限元模型即可自动求解颤振临界风速。

(2) 生成初始训练集。在标准正态空间内,采用拉丁超立方抽样方法生成ndeo个样本点,并将样本点转换成相应非标准正态变量Xi(i=1,2,…,ndeo),代入有限元模型和(8)式,得到初始训练集:

Sdeo={(Ui,G(Xi))|i=1,2,…,ndeo}

(12)

初始训练集的数目通常根据问题的维度和非线性程度确定,也可进行估算[21],具体公式如下:

ndeo=max{2n+1,0.5nmax}

(13)

其中:n为随机变量维度;nmax为允许调用有限元分析的次数。

(3) 生成MC样本。使用步骤(2)抽样方法,在标准正态空间内生成nMC个样本点。此时,这些样本点无需调用有限元模型进行计算。

(4) 创建RBF模型。基于初始训练集Sdeo,使用交叉验证方法确定最优参数,训练得到初始RBF神经网络模型。此处所用工具是MATLAB Neural Network工具箱。

(5) 更新训练集。将MC样本代入RBF模型求解,按照3.2节定义的自适应采样策略选取最优样本点,代入真实功能函数求解,不断更新训练集Sdeo,建立新的RBF模型。

(7) 收敛准则判定。若迭代过程中,2次计算的失效概率相对变化率小于容许范围ε(通常取0.005),且添加的新训练样本点数目大于最小允许数目,则进入下一步,否则重复步骤(5)～步骤(6)。

由上述步骤可以看出,本文建立的颤振可靠度方法并没有对功能函数和随机变量的性质做特定假设。因此,只需要将步骤(1)所建立的颤振分析模型更换为其他相应的数值模型,即可将本文方法推广到其他可靠度分析问题中。

4 数值算例

本文选用一个具有非线性功能函数的6维可靠度分析问题和一个具有理想平板断面的简支梁颤振失稳可靠度问题作为数值算例。其中,以算例1为例,将本文所提基于加权抽样和动态距离约束的自适应采样策略和现有方法进行比较,分析该策略的精度和计算效率;算例2则用于验证本文所提颤振可靠度分析方法的可行性。

4.1 算例1

本算例为非线性振荡器动力学性能的可靠度分析问题,非线性动力系统如图4所示,功能函数如下:

(14)

图4 非线性动力系统

各变量的统计参数见表1所列。

表1 算例1变量的统计参数

为对比分析本文所提自适应采样策略的精度、计算效率及初始样本点的数目对其影响,本算例以105次MC模拟结果作为精确解,并使用同样的MC样本,考虑初始样本点ndeo为5、10、15的3种情况,将本文方法分别独立执行10次得到平均失效概率,具体计算结果见表2所列。由表2可知,本文所提自适应采样策略与已有方法在精度方面大致相当,但在计算效率上具有一定优势。

不同初始样本点数目的收敛情况如图5所示,结合表2发现,随着初始样本点数目的增加,计算精度无明显变化,计算过程收敛速度则逐渐加快, 3种情况下总的计算次数无明显变化,总体均在40次以下。

表2 算例1不同方法计算结果的比较

图5 失效概率收敛情况

4.2 算例2

本算例为具有理想平板截面的简支梁桥。该简支梁桥长为300 m,截面宽为40 m,截面面积为8 m2,竖向抗弯惯性矩为10 m4,横向抗弯惯性矩为85.714 m4,每延米质量惯性矩为4.5×106kg·m2/m,不考虑阻尼的影响。

表3 算例2变量统计信息

该简支梁桥有限元模型如图6所示,系统由30个梁单元(Beam4),29个质量单元(Mass25)和29对气动单元(Matrix27)构成。将上述随机变量取均值,由ANSYS有限元分析软件计算得到颤振系统临界风速确定性分析结果为139.64 m/s,与理论解139.90 m/s十分接近[24],证明了有限元模型的正确性和有效性,可用于后续计算。

图6 简支梁有限元模型

为验证本文所提颤振可靠度分析方法的可行性,首先采用拉丁方超立方抽样方法生成105个样本点,代入上述有限元模型求解颤振临界风速,得到失效概率为8.88×10-3,并以此作为该简支梁颤振失效概率的精确解;然后按照本文所提计算流程,采用同样的MC样本,将初始训练集数目设定为100,同样将本文方法独立执行10次;最终由本文方法得到的平均失效概率为8.81×10-3,具体计算结果见表4所列。由表4可以看出,本文方法计算次数仅为369次,较大程度上减少了MC方法的计算成本。

表4 算例2计算结果

计算过程中本文方法与MC方法样本点分布情况如图7所示。

图7 本文方法与MC方法样本点示意图

从图7可以看出,新添加的样本点分布在极限状态面附近,这表明本文所提出的基于加权抽样和动态距离约束的自适应采样策略能够准确识别靠近系统极限状态面的样本点。

5 结论

针对因颤振失稳功能函数高维非线性和非显性导致颤振失效概率难以求解的问题,本文采用RBF神经网络替代复杂耗时的全模态颤振临界风速搜索过程,并引入自适应随机采样方法以提高神经网络建模效率,提出了一种直接计算失效概率的颤振可靠度分析方法。本文得出的主要结论如下:

(1) 本文提出了一种基于加权抽样和动态距离约束的自适应采样策略。数值算例表明该采样策略能够准确识别靠近极限状态面的样本点,相较于已有的自适应采样方法,在精度相当的前提下,该采样策略调用功能函数计算的次数更少。

(2) 基于定义的自适应采样策略,本文提出的颤振可靠度分析方法将使用RBF神经网络在变量设计空间全域拟合颤振失稳功能函数的问题转换为拟合极限状态面的问题,在保证计算精度的前提下,有效地降低了计算复杂度和计算成本。

(3) 考虑到获取MC精确解的计算成本问题,本文使用具有理想平板断面简支梁桥,而非实桥,作为验证本文方法有效性和准确性的算例。实际上,实桥的颤振分析原理和过程与本文算例是一致的,只是模型复杂程度和计算耗费更大。

(4) 本文方法对于结构功能函数和随机变量性质未做特定假设。因此,理论上通过修改该方法的数值模型,可将本文方法应用于工程中其他的可靠度问题。