CN108334670A - 一种基于摄动法的复杂模型不确定性有限元模型修正方法 - Google Patents

Abstract

本发明提出一种基于摄动法的复杂模型不确定性有限元模型修正方法，包括：推导复杂模型待修正参数的均值和协方差矩阵的迭代形式；确定复杂模型修正目标均值和协方差矩阵；确定复杂模型待修正参数初始估计值；通过有限元分析计算复杂模型待修正参数的灵敏度矩阵及模态频率；读取复杂模型有限元计算结果，计算中间变量转换矩阵及协方差矩阵；分别对复杂模型待修正参数的均值和方差进行修正，得到更新后的参数；循环迭代进行参数修正。本发明基于商用数学软件MATLAB和有限元仿真软件NASTRAN，建立了文件读写和数据传输的通道；适用于解决复杂模型的系统参数与试验数据均存在不确定性的模型修正问题，修正后参数的均值精度较高。

Description

一种基于摄动法的复杂模型不确定性有限元模型修正方法

技术领域

本发明涉及不确定性有限元模型修正方法，特别涉及一种基于摄动法的复杂模型不确定性有限元模型修正方法。

背景技术

准确的结构动力学模型是动响应预测、力学环境预示、损伤检测以及动态优化设计的前提。近几十年来，模型修正技术作为建立精确动力学模型的有效方法得到迅速发展，并在工程领域广泛应用于改善有限元模型与实验数据的相关性、提高有限元模型的分析精度和可靠性。

传统的有限元模型修正仅采用某次试验的数据修正初始有限元模型，目前已经有较多的方法。然而，实际问题中一般存在两类不确定性：(1)结构参数的不确定性，例如新型复合材料的弹性参数、螺栓和铆钉等连接方式以及边界条件难以准确模拟，工程结构在试验或运行过程中动力学特性发生变化等；(2)试验数据的不确定性，在试验过程中由于试验噪声、试验设备及试验条件影响所产生的不确定性，例如，热模态试验目前在试验方法和试验手段上并不成熟，条件难以控制，导致试验结果离散性较大。由于众多不确定性的存在，针对某次试验数据的模型修正结果只能模拟结构特定情况下的动力学特性，而无法对结构的动力学特性进行完整的描述；因此，在获得多次试验数据样本情况下，考虑不确定性的模型修正方法能够得到具有统计特征的动力学模型。

模型修正分为对有限元模型的矩阵进行修正的矩阵型修正方法和直接对结构的物理、几何参数进行修正的参数型修正方法。由于矩阵型修正方法物理意义不明确，人们在工程应用中广泛采用参数型修正方法。其基本思路为，建立理论模型后采用优化方法通过对结构参数进行修正使理论模型与结构模型的动力特性参数误差最小。而参数性模型修正通常是基于灵敏度分析展开的，需要将工程问题抽象为低阶的数学模型，并利用数值方法求解近似解，而摄动法就是一种重要的近似方法。摄动法的常规做法为：对微分方程的变量加以微小的变化，之后对方程进行推导，得到一个原微分方程的近似表达式，用以来代替原方程的解析解。采用摄动法进行模型修正时，若待修正参数不确定性较小并且服从高斯分布，能高效地得到修正结果。

目前，基于摄动法的不确定性模型修正仅停留在理论分析及简单的三自由度质量-弹簧系统层面上，而对于复杂模型未有涉及，原因有以下几点：

a、对于复杂模型的模态分析及灵敏度求解，传统的数值分析方法很难实现，通常需要借助有限元软件来完成；

b、模型修正问题是一个反复迭代的过程，每一个迭代步都需要将新的结构参数输入到有限元软件中进行模态分析及灵敏度分析等；

c、每一个迭代步都需要将有限元分析结果读取出来，用于优化算法中，修正得到新的结构参数；

d、目前有限元分析软件和优化算法之间并没有一套通用的数据接口。

发明内容

本发明的目的在于克服现有技术之不足，提出一种基于摄动法的复杂模型不确定性有限元模型修正方法，基于商用数学软件MATLAB和有限元仿真软件NASTRAN，建立了文件读写和数据传输的通道；适用于解决复杂模型的系统参数与试验数据均存在不确定性的模型修正问题，修正后参数的均值精度较高。

本发明解决其技术问题所采用的技术方案是：

一种基于摄动法的复杂模型不确定性有限元模型修正方法，包括：

步骤101，推导复杂模型待修正参数的均值和协方差矩阵的迭代形式；

步骤102，确定复杂模型修正目标均值和协方差矩阵；

步骤103，确定复杂模型待修正参数初始估计值；

步骤104，通过有限元分析计算复杂模型待修正参数的灵敏度矩阵及模态频率；

步骤105，读取复杂模型有限元计算结果，计算中间变量转换矩阵及协方差矩阵；

步骤106，分别对复杂模型待修正参数的均值和方差进行修正，得到更新后的参数；

步骤107，循环迭代进行参数修正。

优选的，所述步骤101包括：

根据确定性模型修正理论，将待修正参数的估计值θ_j+1用先验值θ_j进行迭代更新，如下

θ_j+1＝θ_j+T_j(z_m-z_j) (1)

其中，z_j是在先验值θ_j下的系统输出值，z_m是试验测得的系统输出值，T_j为转换矩阵，j为大于等于1的正整数；

现考虑待修正参数及系统输出值的不确定性，以结构的模态参数为系统输出值并将其作为修正目标，对基于摄动法的不确定性模型修正迭代方程进行推导，过程如下。

输出的模态参数表示如下：

其中，表示z_m的均值，表示z_j的均值，△z_m表示z_m的不确定性，△z_j表示z_j的不确定性；

第j次迭代过程中的输入值表示如下：

其中，表示θ_j的均值，△θ_j表示θ_j的不确定性；

考虑不确定性的模型修正问题表述如下：

其中，表示T_j的均值，△T_j表示T_j的不确定性；

应用摄动理论分离关于不确定性因素的零阶项和一阶项得到如下表达：

O(△⁰):

O(△¹):

其中，O(△⁰)表示不确定性因素的零阶项；O(△¹)表示不确定性因素的一阶项；

记△z_m＝[△z_m1,△z_m2,...,△z_mn]^H，其中n为向量△z_m的维度，表示结构的前n阶模态参数，△z_mn为试验测得的第n阶模态参数的不确定性，上角标H为矩阵的转置运算符号；则转换矩阵的不确定性可表示为：

其中，为偏导数运算符，k为正整数；则：

记则由(7)式、(9)式可知：

其中，表示矩阵A_j的转置矩阵；表示矩阵的转置矩阵；

假设测量值和仿真的输入输出值是不相关的，即Cov(△z_m,△z_j)＝0、Cov(△z_m,△θ_j)＝0，推导出参数协方差的迭代形式如下：

优选的，所述步骤102包括：

以模态频率作为修正目标，进行若干次次模态试验，测得结构的前n阶固有频率z_m；将若干组测量结果取均值作为修正目标均值并求出若干次试验固有频率值的协方差矩阵

优选的，所述步骤103包括：

假设待修正参数初始估计值为θ₁服从正态分布N(μ,σ²)，其中μ、σ²分别为其均值和方差，则得到Cov(△θ₁,△θ₁)＝Cov(θ₁,θ₁)＝σ²。

若待修正参数个数为p，则其中μ_p、分别为第p个参数的均值和方差。

优选的，所述步骤104包括：

在Matlab中将步骤103中的待修正参数初始值写入到有限元模型中，提交有限元软件Nastran计算，通过模态分析得到结构模态数据再通过灵敏度分析得到结构各阶模态频率对各个待修正参数的灵敏度矩阵

优选的，所述步骤105包括：

在Matlab中，将步骤104中得到的有限元分析结果文件读取到Matlab工作空间，对灵敏度矩阵求广义逆得到转换矩阵均值如下

其中，表示矩阵的转置矩阵；上角标-1为矩阵的求逆运算符；

若目标模态频率对待修正参数的灵敏度矩阵奇异，则采用求解不适定问题的正则化方法计算转换矩阵，如下：

其中，E为单位矩阵；λ为正则化参数，通过以为x轴，以为y轴作L-curve曲线求得λ的值。

同时求解迭代中间量

优选的，所述步骤106包括：

将步骤102至步骤105中的数据代入公式(6)与公式(9)中，分别对待修正参数的均值和方差进行修正，得到更新后的参数。

优选的，所述步骤107包括：

将更新后的参数通过Matlab写入到有限元软件中进行有限元分析，重复步骤103至步骤106，直到修正后的有限元模态频率与试验模态频率的误差满足设定的精度值，退出迭代循环，模型参数修正完成。

本发明具有如下有益效果：

(1)本发明方法实现了计算软件Matlab与有限元分析软件Nastran之间的数据传输互通，并基于此实现了摄动法在任意复杂模型的不确定性修正问题当中的应用；采用非对称H型梁模型对该方法的可行性进行了验证，实验结果表明，本发明方法修正精度较高且收敛较快；

(2)本发明方法首先推导了待修正参数的均值和协方差矩阵的迭代形式；然后建立了MATLAB与NASTRAN的双向数据传输接口；再利用有限元分析结果对设计参数进行多次迭代修正，直至满足精度要求。

以下结合附图及实施例对本发明作进一步详细说明，但本发明的一种基于摄动法的复杂模型不确定性有限元模型修正方法不局限于实施例。

附图说明

图1为本发明方法主流程图；

图2为本发明实施非对称H型梁模型及其尺寸示意图；

图3为本发明实施例修正目标频率对待修正参数的灵敏度矩阵图示；

图4为本发明实施例参数均值修正收敛图；

图5为本发明实施例参数变异系数修正收敛图。

具体实施方式

参见图1所示，本发明一种基于摄动法的复杂模型不确定性有限元模型修正方法，包括：

步骤101，推导复杂模型待修正参数的均值和协方差矩阵的迭代形式；

步骤102，确定复杂模型修正目标均值和协方差矩阵；

步骤103，确定复杂模型待修正参数初始估计值；

步骤104，通过有限元分析计算复杂模型待修正参数的灵敏度矩阵及模态频率；

步骤105，读取复杂模型有限元计算结果，计算中间变量转换矩阵及协方差矩阵；

步骤106，分别对复杂模型待修正参数的均值和方差进行修正，得到更新后的参数；

步骤107，循环迭代进行参数修正。

本实施例中，以非对称H型梁模型为例进行说明。具体的，运用本发明中的方法，以材料的杨氏模量为待修正参数，以模态频率为修正目标，实施基于摄动法的复杂模型不确定性有限元模型修正。梁结构模型参见图2所示，模型被划分为12个单元，具体实施过程包括：

1、灵敏度分析，确定待修正参数及修正目标。

计算模型前8阶模态频率对12个单元的材料杨氏模量的灵敏度矩阵，取灵敏度较高的材料参数和模态频率作为待修正参数和修正目标。本实施例中取单元5、6、9、10的杨氏模量为待修正参数，第1、2、4、5、6、8阶模态频率为修正目标，进行不确定性模型修正，其灵敏度矩阵参见图3所示。

2、确定待修正参数初始值及修正目标。

取待修正参数初始值为：

E₅～N(6.2,1.24²)×10¹⁰N·m^-2，E₆～N(6.2,1.24²)×10¹⁰N·m^-2

E₉～N(8.5,1.7²)×10¹⁰N·m^-2，E₁₀～N(8.5,1.7²)×10¹⁰N·m^-2

参数的初始变异系数(均值与方差之比)均为0.2。

参数修正目标值为：

E₅～N(7.2,0.36²)×10¹⁰N·m^-2，E₆～N(7.2,0.36²)×10¹⁰N·m^-2

E₉～N(7.2,0.36²)×10¹⁰N·m^-2，E₁₀～N(7.2,0.36²)×10¹⁰N·m^-2

参数的目标变异系数(均值与方差之比)均为0.05。

依据参数目标值对参数进行抽样，以10000次抽样为例，可得到10000组参数样本，将其分别带入有限元模型中可获得10000组模态频率样本，如表1所示，计算模态频率样本的均值和协方差，作为修正目标。

表1

3、实施本发明中的复杂模型不确定性修正方法，得出修正结果

非对称H型梁模型的修正结果参见图4和图5所示，分别为均值和变异系数(均值与方差之比)的15次迭代收敛过程，可以看出，在经过6-7次迭代修正后，参数的均值和方差均以较高的精度收敛到修正目标值，验证了本发明方法的实用性及可靠性。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种基于摄动法的复杂模型不确定性有限元模型修正方法，其特点在于，包括：

步骤101，推导复杂模型待修正参数的均值和协方差矩阵的迭代形式；

步骤102，确定复杂模型修正目标均值和协方差矩阵；

步骤103，确定复杂模型待修正参数初始估计值；

步骤104，通过有限元分析计算复杂模型待修正参数的灵敏度矩阵及模态频率；

步骤105，读取复杂模型有限元计算结果，计算中间变量转换矩阵及协方差矩阵；

步骤106，分别对复杂模型待修正参数的均值和方差进行修正，得到更新后的参数；

步骤107，循环迭代进行参数修正。

2.根据权利要求1所述的基于摄动法的复杂模型不确定性有限元模型修正方法，其特征在于，所述步骤101包括：

根据确定性模型修正理论，将待修正参数的估计值θ_j+1用先验值θ_j进行迭代更新，如下θ_j+1＝θ_j+T_j(z_m-z_j) (1)

其中，z_j是在先验值θ_j下的系统输出值，z_m是试验测得的系统输出值，T_j为转换矩阵，j为大于等于1的正整数；

现考虑待修正参数及系统输出值的不确定性，以结构的模态参数为系统输出值并将其作为修正目标，对基于摄动法的不确定性模型修正迭代方程进行推导，过程如下。

输出的模态参数表示如下：

其中，表示z_m的均值，表示z_j的均值，△z_m表示z_m的不确定性，△z_j表示z_j的不确定性；

第j次迭代过程中的输入值表示如下：

其中，表示θ_j的均值，△θ_j表示θ_j的不确定性；

考虑不确定性的模型修正问题表述如下：

其中，表示T_j的均值，△T_j表示T_j的不确定性；

应用摄动理论分离关于不确定性因素的零阶项和一阶项得到如下表达：

其中，O(△⁰)表示不确定性因素的零阶项；O(△¹)表示不确定性因素的一阶项；

记△z_m＝[△z_m1,△z_m2,...,△z_mn]^H，其中n为向量△z_m的维度，表示结构的前n阶模态参数，△z_mn为试验测得的第n阶模态参数的不确定性，上角标H为矩阵的转置运算符号；则转换矩阵的不确定性可表示为：

其中，为偏导数运算符，k为正整数；则：

记则由(7)式、(9)式可知：

其中，表示矩阵A_j的转置矩阵；表示矩阵的转置矩阵；

假设测量值和仿真的输入输出值是不相关的，即Cov(△z_m,△z_j)＝0、Cov(△z_m,△θ_j)＝0，推导出参数协方差的迭代形式如下：

3.根据权利要求2所述的基于摄动法的复杂模型不确定性有限元模型修正方法，其特征在于，所述步骤102包括：

以模态频率作为修正目标，进行若干次次模态试验，测得结构的前n阶固有频率z_m；将若干组测量结果取均值作为修正目标均值并求出若干次试验固有频率值的协方差矩阵

4.根据权利要求3所述的基于摄动法的复杂模型不确定性有限元模型修正方法，其特征在于，所述步骤103包括：

假设待修正参数初始估计值为θ₁服从正态分布N(μ,σ²)，其中μ、σ²分别为其均值和方差，则得到Cov(△θ₁,△θ₁)＝Cov(θ₁,θ₁)＝σ²。

若待修正参数个数为p，则其中μ_p、分别为第p个参数的均值和方差。

5.根据权利要求4所述的基于摄动法的复杂模型不确定性有限元模型修正方法，其特征在于，所述步骤104包括：

在Matlab中将步骤103中的待修正参数初始值写入到有限元模型中，提交有限元软件Nastran计算，通过模态分析得到结构模态数据再通过灵敏度分析得到结构各阶模态频率对各个待修正参数的灵敏度矩阵

6.根据权利要求5所述的基于摄动法的复杂模型不确定性有限元模型修正方法，其特征在于，所述步骤105包括：

在Matlab中，将步骤104中得到的有限元分析结果文件读取到Matlab工作空间，对灵敏度矩阵求广义逆得到转换矩阵均值如下

其中，表示矩阵的转置矩阵；上角标-1为矩阵的求逆运算符；

若目标模态频率对待修正参数的灵敏度矩阵奇异，则采用求解不适定问题的正则化方法计算转换矩阵，如下：

其中，E为单位矩阵；λ为正则化参数，通过以为x轴，以为y轴作L-curve曲线求得λ的值。

同时求解迭代中间量

7.根据权利要求6所述的基于摄动法的复杂模型不确定性有限元模型修正方法，其特征在于，所述步骤106包括：

将步骤102至步骤105中的数据代入公式(6)与公式(9)中，分别对待修正参数的均值和方差进行修正，得到更新后的参数。

8.根据权利要求7所述的基于摄动法的复杂模型不确定性有限元模型修正方法，其特征在于，所述步骤107包括：

将更新后的参数通过Matlab写入到有限元软件中进行有限元分析，重复步骤103至步骤106，直到修正后的有限元模态频率与试验模态频率的误差满足设定的精度值，退出迭代循环，模型参数修正完成。