基于Kriging代理模型的中频结构振动可靠性分析

方琛 秦朝红 任方 张亚辉

(1 大连理工大学 工程力学系 工业装备结构分析国家重点实验室 大连 116024;2 北京强度环境研究所 可靠性与环境工程技术重点实验室 北京 100076)

0 引言

工程结构在制造、运输、装配及服役过程不可避免地存在各种不确定性因素，所以在动态激励作用下单纯地计算结构的动力学响应或者响应的极值会忽视不确定性的影响。随着计算能力和精细化建模需求的提升，结合概率论准确描述不确定性因素对动力学响应的影响，开展结构的动力学可靠度分析受到了广泛关注。低频问题在采用诸如有限元（Finite Element Method，FEM）等数值方法实现模型化后，通过概率模型结合抽样或矩方法等即可得到结构的失效概率[1]。随着频率的升高，继续采用FEM建模需要将网格划分得足够细，这会导致庞大的自由度数。在计算量难以负担的同时，还会出现参数敏感性问题。当频率足够高时，充分考虑不确定性的集总参数的统计能量分析[2]（Statistical Energy Analysis，SEA）应运而生。可以通过响应的方差及近似概率分布描述非参数随机的SEA模型的离散性[2-3]。

然而，实际工程结构在宽频激励下，不同组件的振动波长、模态密度等动力学特性差异较大，就会表现出中频振动行为。传统的低频或高频建模方法都无法直接应用。目前主要有两种解决方案，一是传统方法向中频段的扩展[4-5]，二是基于子系统思想分别采用不同建模方法得到的中频混合模型。其中Shorter和Langley基于扩散场互易关系[6]提出的混合有限元-统计能量分析方法（Hybrid Finite Element-Statistical Energy Analysis method，简称混合FE-SEA方法）[7]理论最为完善，在诸如航空航天等领域应用广泛[8-9]。后续学者尝试采用不同的方法替换FEM和SEA建立了不同的混合模型[10],[11]。为了实现中频结构振动的可靠性分析，本文将采用混合FE-SEA方法实现模型化。

需要强调的是，混合FE-SEA方法将短波组件定义为统计性子系统，采用了集总参数的SEA建模，计算得到响应是已经考虑了短波组件非参数不确定性的。而长波组件是完全确定性的，与实际不符。随着现在计算能力的发展和对不确定性研究的需求，安全系数法已经不能满足现有结构的设计分析需求。Cicirello基于混合FE-SEA方法首次开展了中频振动的结构可靠性分析[12]。Cicirello分别采用蒙特卡洛模拟（Monte-Carlo Simulation，MCS）和渐进近似拉普拉斯方法计算得到了弹簧-板组合系统的中频结构可靠度。

混合FE-SEA方法的计算受制于复数总动刚度矩阵的求逆运算及谱分析形式，同等自由度规模下计算负担远大于FEM。那么对混合模型直接采用MCS求解可靠度的将进一步加剧计算负担。此外基于混合FE-SEA方法建立的中频问题的功能函数也格外复杂，采用传统的可靠性分析方法求解比较困难。本文采用Kriging代理模型[13]-[15]直接建立随机参数与响应的关系，从而避免复杂的原物理模型重分析。Kriging模型的精度依赖于样本点的选取，对于混合方法建立的复杂的中频振动求解模型，需要选取大量的样本点，这样会产生较大的计算量。为了在减少样本数量的前提下保证计算精度，可以采用拉丁超立方试验设计[16]（Latin Hypercube Sampling，LHS）得到布满样本空间的试验组合，实现Kriging代理模型的全满空间设计。故本文主要工作是基于混合FE-SEA方法的模型化，进一步考虑随机参数构建了可靠性分析的功能函数，采用LHS结合Kriging模型构建代理模型实现中频振动的结构可靠度分析。

1 中频振动的混合FE-SEA模型

要实现中频振动的结构可靠性分析，首先需要建立中频振动的响应分析模型。本文采用经典的混合FE-SEA方法实现中频振动的模型化。混合FE-SEA方法将复杂系统的各个组件区分为长波变形子系统和短波变形子系统，然后分别采用FEM和SEA建模，最后通过扩散场互易关系实现系统响应的耦合求解。本节将对中频结构的混合FE-SEA方法的基本理论进行介绍。

1.1 总运动方程

如图1所示，混合FE-SEA方法将短波组件定义为统计性子系统，认为其是“直接场”与“混响场”的叠加。然后将混响场对直接场确定性边界的作用定义为“混响场荷载”[6]。

图1 混合FE-SEA模型的耦合示意图 Fig. 1 Coupling of the hybrid FE-SEA method

为开展可靠性分析，本文进一步通过概率模型考虑长波组件的参数随机特性，并将长波组件定义为随机子系统。首先建立系统的总运动方程，其描述的是直接场效应和长波组件，写作

其中()H表示共轭转置，E[]表示集合平均，nj为第j个统计性子系统的模态密度，ω为激励频率。

为了联立式(2)所示的扩散场互易关系，求总运动方程(1)建立的自由度q的互谱，表达式为

1.2 功率流平衡方程

根据功率流平衡关系，第j个统计性子系统的功率流平衡方程表示为[7]

hij表示第j个统计性子系统混响场的单位模态密度输入到第i个统计性子系统直接场的功率，表达式为

综合式(5)～(7)，将m个统计性子系统的功率流平衡方程联立为矩阵形式，表示为

E是由m个统计性子系统的能量组成的未知向量，是随机子系统对各统计性子系统的输入功率，N是各统计性子系统的模态密度组成的对角矩阵。L为系数矩阵，该系数矩阵的表达式为

求解式(8)即可得到所有统计性子系统的能量，回代到式(3)即可求解得到的随机子系统的响应Sqq。从系数矩阵L的表达式可以发现，混合模型求解需要首先计算总动刚度矩阵的逆，然后进行大量的求和运算，计算负担很大。

2 中频振动的结构可靠性分析

上一节通过混合FE-SEA方法实现了中频振动的模型化，随后即可开展结构可靠性分析。本节将建立可靠性分析的功能函数，并通过构造Kriging代理模型以避免复杂且计算量庞大的原物理模型的求解。

2.1 功能函数列式

基于混合FE-SEA方法进行结构中频振动的模型化时，统计性子系统采用集总参数的SEA方法建模时充分考虑了非参数随机性。在忽略了短波组件离散性基础上，本文假设可靠性分析的n个随机参数x=[x1,…xn]来源于定义为随机子系统的长波组件的相关参数。将通过混合FE-SEA中频随机模型求解的响应表示为Q(x)。进一步可以将结构的极限状态通过功能函数来表示

G(x)由结构的集合平均响应Q(x)和响应界限值B来决定。Z＞0、Z＜0和Z=0就表示结构处于安全、失效和极限状态。通过采样方法计算可靠度表达式为

即可靠度为安全状态的样本数np(Z＞0)占总样本数nsum的百分比。

2.2 Kriging代理模型

中频问题的功能函数(10)由于响应Q(x)与随机参数x=[x1,…xn]之间运算关系格外复杂，难以采用传统的动力学可靠度方法求解。同时采用混合FE-SEA方法建立1.2节的功率流平衡方程时，必须直接对总动刚度进行求逆，然后进行一系列求和运算。如果随机子系统的自由度规模较大，则计算负担将难以承受。故采用Kriging方法建立等效响应求解的代理模型，结合满样本空间的LHS控制样本点的选取，就可以建立式(10)所需的响应求解的等效方程。

Kriging代理模型是一个基于线性回归分析建立的改进拟合模型。该模型包括回归项和随机项两部分，对于一个函数y(x0)，其可以表示为

x0表示待测点；fT(x0)表示给定的多项式回归模型；β表示待定的回归系数；z(x0)表示随机过程。z(x0)的均值为0，标准差为σz，协方差表示为

R(x1,x2)表示样本点中任意两点x1、x2之间与距离dj有关的变异函数。R(x1,x2)由核函数计算，有诸多相关方程[17]，广泛使用的是高斯相关方程，写作

θj代表各向异性参数；p∈[0,2]表示光滑程度参数。

通过样本数据X和对应响应Y，可以得到相关函数矩阵R。考虑待测点x0与N个样本点X =[x1,…xN]的相关向量r (x0)，采用最大似然估计得到θj优化值。进一步采用最小二乘法可得最大似然估计因子

z(x)的方差估计表达式为

则Kriging代理模型在待测点x0处的预测值为

至此建立得到Kriging代理模型。

建立的Kriging代理模型的计算精度与效率依赖于样本点的选取，为了有效的选择充满空间的少量样本点，本文选择比较主流的LHS[16]采样方法。通过LHS建立满空间的nk个样本点，然后通过物理模型计算得到相应响应Q(x)。最后结合式(17)就可以建立Kriging代理模型，即一个近似的响应分析模型。

首先，建立代理模型所需的LHS样本点要远少于直接针对物理模型进行MCS的样本点数。其次代理模型的计算效率要远高于原物理模型，所以针对代理模型进行MCS即可迅速得到响应可靠度。本文采用成熟的DACE工具箱[18]建立Kriging模型，该工具箱简便实用，被工程师和研究者们广泛采用。

3 数值算例

本节将对如图2所示的一个框架梁栓接板组合系统开展中频结构动力学可靠度分析。该结构由一个边长0.7m为正方体框架梁和四个长1.1m、宽0.6m的矩形板组成。梁截面为矩形，各板通过四个点栓接在框架梁上，各栓接点距离框架梁顶点0.1m。在原点o作用一个沿y轴正向的200Hz的单频点荷载。框架梁的材料为钢，阻尼损耗因子为0.01，板为铝，阻尼损耗因子为0.02。

图2 框架梁栓接板结构 Fig. 2 A framework-panel system

该模型其余材料属性、几何参数、荷载幅值等均为随机参数，并考虑环境影响的不确定性。相关随机参数的具体统计信息及概率分布见表1。

表1 随机变量的概率分布参数 Table 1 Probability distribution parameters of random variables

为开展可靠性分析，首先需要基于混合FE-SEA方法实现模型化。而混合模型的建模需要首先根据振动波长进行子系统划分，200Hz激励下梁和板的波长见表2，计算采用名义参数。如图3 a)，将长波变形的梁和板面内振动定义为随机子系统，采用FEM建模，网格尺寸为0.05m。如图3 b)，采用SEA建模，将板短波变形的弯曲振动定义为统计性子系统。

表2 振动波长 Table 2 Vibrational wavelength

图3 混合FE-SEA模型的划分 Fig. 3 Subsystems of hybrid FE-SEA model

通过混合FE-SEA方法实现模型化后，即可开展可靠性分析。对于随机子系统，选择梁上的测点A(0.5, 0.7, 0.7)的速度幅值互谱Svv=1.66 ×10-5m2/s2为临界变形。对于统计性子系统，选择z=0.7m平面上的板的能量E4=5.42 ×10-4J 作为临界能量。

对混合FE-SEA方法建立的原物理模型进行nm=105个样本点的MCS[12]，将其计算结果作为可靠度参考解。对物理模型选取数量nk=500,600，1000的LHS样本（随机参数置信区间不小于99%）建立Kriging代理模型，可靠度计算结果如表3所示。通过nk=500个样本点建立的代理模型计算的可靠度已经接近MCS的结果，误差均小于2%。随着LHS样本点的增多，误差逐渐减小到低于0.5%。值得注意的是，该混合模型单次采样计算都非常耗时，其中大部分时间耗费在总动刚度矩阵的求逆及相关的求和运算中。而采用代理模型计算近似响应的时间可以忽略不计。所以采用LHS采样建立Kriging代理模型，求解可靠度的方法通过较小的计算量实现了中频振动的结构可靠度求解。即使是LHS采样数nk=1000，代理模型的计算时间只是nm=105次MCS的1%而已。而LHS采样数为nk=500,600次时间，计算效率可以在MCS基础上分别提升200、167倍。

表3 结构可靠度分析结果(%) Table 3 Reliability analysis results (%)

进一步验证更高可靠度下的结果。将测点A(0.5, 0.7, 0.7)的速度幅值互谱和z=0.7m平面上的板能量对应的临界值重新给定为Svv=1.94 ×10-5m2/s2和E4=5.71 ×10-4J 。计算结果如表4所示，此时MCS得到的可靠度均为99%。采用nk= 1000个LHS样本建立了Kriging代理模型，计算得到的可靠度误差均小于1%。

表4 结构可靠度分析结果(%) Table 4 Reliability analysis results(%)

4 结论

本文开展了中频振动的结构可靠性分析。混合FE-SEA方法将短波组件采用SEA建模，本文在忽略短波组件的离散性后通过概率模型描述长波变形组件的随机性来开展可靠性分析。通过满空间的LHS选取样本点即可建立Kriging代理模型。针对一个框架梁栓接板组合结构计算可靠度。发现较之对物理模型通过大量采样的直接MCS的方法，通过少量的LHS样本点建立的代理模型可以迅速计算得到可靠度，计算效率很高。随着LHS样本点增多，误差会逐渐减小。