CN103064999A - 一种用于抽水蓄能电站地下厂房结构的模型修正方法 - Google Patents

Abstract

本发明涉及一种用于抽水蓄能电站地下厂房结构的模型修正方法。其技术方案是：首先根据抽水蓄能电站地下厂房结构的先验信息，确定该地下厂房结构需要修正的设计变量x，通过实验测试得到特征响应f_cs；其次，在先验信息的指导下进行试验设计，得到k个样本点，再利用有限元软件进行计算得到特征响应f_FEM；再次，利用最小二乘法确定二次响应面函数g(x)的系数，构造响应面并进行自身特性检验和精度检验；最后，对满足条件的二次响应面函数g(x)进行优化，更新加权系数至满足收敛条件和模型修正的评价指标。本发明具有计算效率高和修正效果好的特点，本发明还适用于其它大型厂房结构的模型修正。

Description

一种用于抽水蓄能电站地下厂房结构的模型修正方法

技术领域

本发明属于模型修正方法技术领域。具体涉及一种用于抽水蓄能电站地下厂房结构的模型修正方法。

背景技术

数值仿真和模型试验是进行结构动力分析设计的两种主要方法。对于大型结构，如航天器、抽水蓄能电站地下厂房等，由于其结构的复杂性（荷载和边界条件的不确定性等）而造成的模型误差，使得有限元模型预示的结果与实验测试结果之间存在明显差异，需要对其结构模型进行修正以正确预示结构动力行为。目前的共识是，经过修正后的结构模型是进行荷载识别、结构损伤识别、结构健康监测、结构性能评价和优化设计的基础。

模型修正属于反问题的一种，反问题中的不确定性、非线性、计算量大等特性在模型修正问题中都是存在的。从最早出现的矩阵修正技术到近年发展的参数修正技术，模型修正法在不断完善，但是对于大型复杂结构的模型修正方法，主要还存在以下三个难点：第一，有限元模型是数值形式，不是函数形式，难以找到符合设计要求的功能函数和影响因素自变量之间的关系，即使存在这种关系，对有限元中各个矩阵的修正也会导致其稀疏性消逝，从而难以得到实际应用；第二，有限元模型误差来源与实测结果的误差来源完全不同，误差性质缺乏可比性，直接利用实测结果对有限元模型进行修正具有盲目性；第三，结构模型修正是一个渐进的优化过程，需要经历一系列的修正和检验的步骤，计算效率有待提高。因此，对于大型复杂结构，传统的模型修正技术难以真正实施。

相对于常规水电站而言，抽水蓄能电站具有水头高、容量大、双向运行、工况转换频繁等特点，因此发电厂房的振动问题较常规电站更为严重，有些甚至因振动造成结构开裂，对水电站安全运行造成威胁。但是，目前国内缺少专门的关于抽水蓄能电站的动力设计规范，采用常规水电站的设计规范又过于笼统和一般化，难以从根本上解决地下厂房的振动问题。

发明内容

本发明旨在克服现有技术缺陷，目的是提出一种计算效率高和修正效果好的用于抽水蓄能电站地下厂房结构的模型修正方法。

为实现上述目的，本发明采用的技术方案的具体步骤是：

步骤1：根据抽水蓄能电站地下厂房结构的先验信息，确定该地下厂房结构需要修正的设计变量，通过实验测试得到特征响应f_cs。

步骤2：根据抽水蓄能电站地下厂房结构的特点，在先验信息的指导下进行试验设计，得到k个样本点，再采用ANSYS有限元软件或ABAQUS有限元软件对厂房结构进行计算得到特征响应f_FEM。

步骤3：对k个样本点进行归一化处理，建立二次响应面函数g(x)

$g\left(x\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{[\mathrm{\γ}{f}_{\mathrm{FEM}}\left(i\right)-f_{\mathrm{cs}}\left(i\right)]}^2={\mathrm{\β}}_0+\underset{i=1}{\overset{k^{\′}}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_i+\underset{i=1}{\overset{k^{\′}}{\mathrm{\Σ}}}{\mathrm{\β}}_{\mathrm{ii}}{x}_i^2---\left(1\right)$

式（1）中：

x为厂房结构需要修正的设计变量；

f_FEM为利用有限元软件计算得到的特征响应；

f_cs为实验测试得到的特征响应；

γ为二次响应面函数g(x)的加权系数，为采用有限元软件计算得到的信息与实验测试信息的融合程度；

k为样本点的个数；

k′为设计变量x的个数；

β₀为采用最小二乘法所确定的二次响应面函数g(x)中的常数项；

β_i为采用最小二乘法所确定的二次响应面函数g(x)中的一次项系数；

β_ii为采用最小二乘法所确定的二次响应面函数g(x)中的二次项系数。

步骤4：根据式（1），设定初始加权系数γ＝1.0，得到初始二次响应面函数g₀(x)；对初始二次响应面函数g₀(x)进行自身特性检验和精度检验，若不能满足全部条件，则返回步骤2，重新定义k个样本点至满足自身特性检验和精度检验的全部条件；若能满足自身特性检验和精度检验的全部条件，则进入步骤5。

步骤5：对满足要求的二次响应面函数g(x)进行优化计算，得到二次响应面函数g(x)最小值ming(x)所对应的一组设计变量x^*，再利用ANSYS有限元软件或ABAQUS有限元软件进行计算，得到对应的二次响应面函数g(x^*)。

迭代过程中二次响应面函数g(x)迭代第j步的加权系数γ_j

${\mathrm{\γ}}_j=\frac{\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{f}_{\mathrm{FEM}}\left(i\right)\·f_{\mathrm{cs}}\left(i\right)}{{\left[f_{\mathrm{FEM}}\left(i\right)\right]}^2}---\left(2\right)$

式（2）中：

k为样本点的个数；

f_FEM为利用有限元软件计算得到的特征响应；

f_cs为实验测试得到的特征响应。

若二次响应面函数g(x)在迭代过程中满足以下收敛条件：

$\left\{\begin{array}{c}|g\left(x^{*}\right)-\min g\left(x\right)|\≤{\mathrm{\ϵ}}_1\\ |{\mathrm{\γ}}_{j-1}-{\mathrm{\γ}}_j|\≤{\mathrm{\ϵ}}_2\end{array}\right.---\left(3\right)$

式（3）中：

g(x^*)是以x^*为设计变量代入有限元软件计算得到的二次响应面函数g(x)的值；

ming(x)为优化二次响应面函数g(x)后得到的最小值；

ε₁和ε₂为收敛阈值。

则设计变量x^*收敛，停止计算，进入步骤7，否则进入步骤6。

步骤6：按照式（2）更新加权系数γ_j，再返回步骤3，构造新的二次响应面函数g′(x)，重复步骤4和步骤5至满足式（3）的收敛条件。

步骤7：经过j步迭代，反归一化处理后，输出优化后的二次响应面函数g(x)所得到的设计变量

步骤8：建立模型修正后的评价指标

$\mathrm{ER}(f_{\mathrm{cs}},f_{\mathrm{mod}})=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}\left|\frac{f_{\mathrm{cs}}\left(i\right)-f_{\mathrm{mod}}\left(i\right)}{f_{\mathrm{cs}}\left(i\right)}\right|\×100\%---\left(4\right)$

式（4）中：

f_cs为实验测试得到的特征响应；

f_mod为经过模型修正后的特征响应；

k为样本点的个数。

若ER(f_cs，f_mod)≤10%时，则表明该模型修正方法的结果接近于真实情况；若ER(f_cs，f_mod)＞10%时，则返回步骤2，直至ER(f_cs，f_mod)≤10%。

所述的设计变量x是利用先验信息得到的设计变量，包括上下游刚度值、左侧粘结力值、约束分布密度值、质量点大小值和质量点坐标值。

所述的先验信息是：

①刚度悬殊：主体结构如蜗壳、机墩和风罩等刚度较大，次要结构如板梁柱等刚度较小；

②结构类型固定：厂房结构为3～4层，为厚板结构或为板梁结构；

③楼板开孔较大：楼板上开有吊物孔和楼梯孔，振型往往围绕开孔处发生；

④约束条件：上下游围岩作用以及左右侧厂房之间的相互作用力对结构的特征频率影响十分明显；

⑤设计规范要求：在激励荷载频率一定的条件下，地下厂房结构的特征频率要满足

所述的特征响应是地下厂房结构的特征频率或模态位移。

所述的试验设计为均匀设计和正交设计中的一种。

所述的自身特性检验是指同时满足如下三个条件是：

①可旋转性：响应量的方差仅是从设计中心到设计点距离的函数，与方向无关；

②斜率可旋转性：与设计中心的等距离点上的响应面的偏导数的方差是一个常数；

③稳健性：在个别数据有误或误差非正态分布的条件下，k个样本点中所含误差的样本点与不含误差的样本点对响应面结果引起的偏差在10%以内。

所述的精度检验是指同时满足如下两个条件：

①误差平方和RMSE→0

$\mathrm{RMSE}=\sqrt{\frac{\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{[f_{\mathrm{FEM}}\left(i\right)-f_{\mathrm{mod}}\left(i\right)]}^2}{k-k^{\′}}}---\left(5\right)$

式（5）中：

k为样本点的个数；

k′为设计变量的个数；

f_FEM为利用有限元软件计算得到的特征响应；

f_mod为经过模型修正后的特征响应。

②响应面精度R²→1

$R^2=\frac{\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{[f_{\mathrm{mod}}\left(i\right)-{\stackrel{-}{f}}_{\mathrm{FEM}}\left(i\right)]}^2}{\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{[f_{\mathrm{FEM}}\left(i\right)-{\stackrel{-}{f}}_{\mathrm{FEM}}\left(i\right)]}^2}---\left(6\right)$

式（6）中：

为利用有限元软件计算得到的特征响应的平均值。

由于采用上述技术方案，本发明与现有技术相比具有如下优点：

第一、传统的模型修正技术是直接对有限元矩阵元素进行修正，这样会导致矩阵元素的对称性、正定性和稀疏性的消逝而难以得到实际应用。本发明采用概率方法，在修正过程中把实测数据和有限元计算数据都视为概率统计中的样本，为模型修正这个反问题的求解提供了一个新的模式。

第二、数值仿真过程中会因结构简化和参数误差等因素引起误差，实验测试也会因外界干扰和设备仪器故障等因素带来误差，两者的误差来源完全不相同，因此直接利用实验测试结果对有限元模型进行修正带有盲目性，而本发明的响应面中采用加权系数用以考虑实测信息和有限元计算信息之间的信息融合，通过响应面自身特性检验和精度检验能保证计算精度，使模型修正更接近真实，修正效果好。

第三、采用响应面模型代替有限元模型进行模型修正迭代计算，能显著提高计算效率，尤其是根据模型修正流程通过二次开发，形成基于计算机语言和有限元APDL命令的模型修正软件，进一步提高计算效率。

因此，本发明具有计算效率高和修正效果好的特点。本发明还适用于其它大型厂房结构的模型修正。

附图说明

图1为本发明的一种有限元模型结构示意图；

图2为图1的剖面示意图；

图3为图1的模型修正过程中的一种响应面示意图；

图4为图1的模型修正过程中的另一种响应面示意图。

具体实施方式

下面结合附图和具体实施方式对本发明作进一步的描述，并非对其保护范围的限制。

实施例1

一种用于抽水蓄能电站地下厂房结构的模型修正方法。其具体步骤是：

步骤1、抽水蓄能电站地下厂房结构的先验信息是：

①刚度悬殊：主体结构如蜗壳、机墩和风罩的刚度较大，次要结构如板梁柱的刚度较小；

②结构类型固定：厂房结构为3～4层，为厚板结构或为板梁结构；

③楼板开孔较大：楼板上开有吊物孔和楼梯孔，振型往往围绕开孔处发生；

④约束条件：上下游围岩作用以及左右侧厂房之间的相互作用力对结构的特征频率影响十分明显；

⑤设计规范要求：在激励荷载频率一定的条件下，地下厂房结构的特征频率要满足

根据先验信息，确定该地下厂房结构需要修正的设计变量分别是：

上下游刚度值x₁（N/m）；左侧粘结力值x₂（N/m）；约束分布密度值x₃（无量纲）；质量点大小值x₄（Kg）；质量点X坐标值x₅（m）；质量点Y坐标值x₆（m）。

本实施例实验测试得到的前五阶特征频率分别为f₁＝12.343Hz、f₂＝14.234Hz、f₃＝16.237Hz、f₄＝18.451Hz和f₅＝20.452Hz。

步骤2、根据抽水蓄能电站地下厂房结构的特点，在先验信息的指导下，采用试验设计得到25个样本点，每个样本点6个因素，每个因素5个水平，采用有限元计算建立如图1和图2所示的地下厂房结构有限元模型，建立范围是沿厂房长度方向取一个机组段，自上而下从发电机层至尾水管底板，包括蜗壳外围混凝土、机墩、风罩和楼板，上下游边界为岩石，左右边界为结构缝，底面按固定约束，顶面为自由面，再采用有限元软件计算出厂房结构的特征频率f_FEM。

本实施例的试验设计采用均匀试验设计，有限元软件为ANSYS有限元软件。

步骤3、对k＝25个样本点进行归一化处理，建立二次响应面函数g(x)

$g\left(x\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{[\mathrm{\γ}{f}_{\mathrm{FEM}}\left(i\right)-f_{\mathrm{cs}}\left(i\right)]}^2={\mathrm{\β}}_0+\underset{i=1}{\overset{k^{\′}}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_i+\underset{i=1}{\overset{k^{\′}}{\mathrm{\Σ}}}{\mathrm{\β}}_{\mathrm{ii}}{x}_i^2---\left(1\right)$

式（1）中：

x为厂房结构需要修正的设计变量，x为步骤1中的x₁～x₆；

f_FEM为步骤1所述利用有限元软件计算得到的前五阶特征频率；

f_cs为步骤1所述实验测试得到的前五阶特征频率f₁～f₅；

γ为二次响应面函数g(x)的加权系数，为采用有限元软件计算得到的信息与实验测试信息的融合程度；

k′为设计变量x的个数，k′=6；

β₀为采用最小二乘法所确定的二次响应面函数g(x)中的常数项；

β_i为采用最小二乘法所确定的二次响应面函数g(x)中的一次项系数；

β_ii为采用最小二乘法所确定的二次响应面函数g(x)中的二次项系数。

步骤4、根据式（1），设定初始加权系数γ＝1.0，得到初始二次响应面函数 $g_0\left(x\right)=0.5x_1^2+1.75x_2^2+1.20x_3^2-0.82x_4^2+2.4x_5^2+3.4x_6^2+0.7x_1+1.34x_2+3.1x_3-4x_4+5.5x_5-{1.8x}_6+12.5.$ 二次响应面函数和其中两个设计变量之间的响应曲面如图3和图4所示。

本实施例对初始二次响应面函数g₀(x)进行自身特性检验，其结果满足可旋转性、斜率可旋转性和稳健性的三个条件：

①可旋转性：响应量的方差仅是从设计中心到设计点距离的函数，与方向无关；

②斜率可旋转性：与设计中心的等距离点上的响应面的偏导数的方差是一个常数；

③稳健性：在个别数据有误或误差非正态分布的条件下，k个样本点中所含误差的样本点与不含误差的样本点对响应面结果引起的偏差在10%以内。

本实施例同时对初始二次响应面函数g₀(x)进行精度检验：

$\mathrm{RMSE}=\sqrt{\frac{\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{[f_{\mathrm{FEM}}\left(i\right)-f_{\mathrm{mod}}\left(i\right)]}^2}{k-k^{\′}}}---\left(2\right)$

式（2）中：

k为样本点的个数，本实施例中k＝25；

k′为设计变量的个数，本实施例中k′=6；

f_FEM为步骤1所利用有限元软件计算得到的前五阶特征频率；

f_mod为经过模型修正后的前五阶特征频率；

经计算，误差平方和RMSE＝0.05，即误差平方和RMSE→0。

$R^2=\frac{\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{[f_{\mathrm{mod}}\left(i\right)-{\stackrel{-}{f}}_{\mathrm{FEM}}\left(i\right)]}^2}{\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{[f_{\mathrm{FEM}}\left(i\right)-{\stackrel{-}{f}}_{\mathrm{FEM}}\left(i\right)]}^2}---\left(3\right)$

式（3）中：

为利用有限元软件计算得到的前五阶特征频率的平均值。

计算得到响应面精度R²＝0.98，即R²→1。

即本实施例同时满足：①误差平方和RMSE→0和②响应面精度R²→1。进入步骤5。

步骤5、对满足要求的二次响应面函数g(x)进行优化计算，得到二次响应面函数g(x)最小值ming(x)＝0.3所对应的设计变量x^*，再利用ANSYS有限元软件进行计算，得到对应的二次响应面函数g(x^*)＝0.6。

根据

${\mathrm{\γ}}_j=\frac{\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{f}_{\mathrm{FEM}}\left(i\right)\·f_{\mathrm{cs}}\left(i\right)}{{\left[f_{\mathrm{FEM}}\left(i\right)\right]}^2}---\left(4\right)$

式（4）中：

k为样本点的个数，本实施例中k＝25；

f_FEM为步骤1所利用有限元软件计算得到的前五阶特征频率；

f_cs为步骤1实验测试得到的前五阶特征频率。

得到迭代过程中二次响应面函数g(x)迭代第2步的加权系数γ₂＝0.981。

根据

$\left\{\begin{array}{c}|g\left(x^{*}\right)-\min g\left(x\right)|\≤{\mathrm{\ϵ}}_1\\ |{\mathrm{\γ}}_{j-1}-{\mathrm{\γ}}_j|\≤{\mathrm{\ϵ}}_2\end{array}\right.---\left(5\right)$

式（5）中：

g(x^*)是以x^*为设计变量代入有限元软件计算得到的二次响应面函数g(x)的值；

ming(x)为优化二次响应面函数g(x)后得到的最小值；

ε₁和ε₂为收敛阈值，本实施例中，阈值取ε₁＝ε₂＝0.1。

将g(x^*)＝0.6、ming(x)＝0.3和|γ_j-1-γ_j|=|0.981-1.0|代入式（3），得到

$\left\{\begin{array}{c}|0.6-0.3|>0.1\\ |0.981-1.0|\≤0.1\end{array}\right.$

即不满足式（5）的收敛条件，则进入步骤6。

步骤6、按照式（4）更新加权系数γ₂，再返回步骤3，构造新的二次响应面函数g′(x)，重复步骤4和步骤5至满足式（5）的收敛条件。

本实施例中，迭代第3步的加权系数γ₃＝0.951，g(x^*)＝0.4，ming(x)＝0.2，不满足式（5）的收敛条件；迭代第4步的加权系数γ₄＝0.960，g(x^*)＝0.18，ming(x)＝0.11，二次响应面函数g(x)在第4次迭代过程中满足式（5）的收敛条件，即满足

$\left\{\begin{array}{c}|0.18-0.11|\≤0.1\\ |0.951-0.960|\≤0.1\end{array}\right.$

步骤7、经过4步迭代，反归一化处理后，输出优化后的二次响应面函数g(x)所得到的设计变量x^*分别为： $x_1^{*}=4.91e7(N/m);$ $x_2^{*}=4.99e7(N/m);$

（无量纲）； $x_4^{*}=1655$ $\left(\mathrm{Kg}\right);$ $x_5^{*}=16.69\left(m\right);$ $x_6^{*}=22.71\left(m\right).$

步骤8、建立模型修正后的评价指标

$\mathrm{ER}(f_{\mathrm{cs}},f_{\mathrm{mod}})=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}\left|\frac{f_{\mathrm{cs}}\left(i\right)-f_{\mathrm{mod}}\left(i\right)}{f_{\mathrm{cs}}\left(i\right)}\right|\×100\%---\left(6\right)$

式（6）中：f_cs为步骤1实验测试得到的前五阶特征频率；

f_mod为经过模型修正后的前五阶特征频率，修正后的特征频率为f₁＝12.211Hz，f₂＝14.114Hz，f₃＝16.263Hz，f₄＝18.272Hz，f₅＝20.431Hz；

k为样本点的个数，本实施例k＝25。

计算得到修正后的ER值2.15%小于修正前的ER值22.54%，修正后的频率更接近实测频率，为工程设计采用。

实施例2

一种用于抽水蓄能电站地下厂房结构的模型修正方法。其具体步骤是：

步骤1、除下述特征响应外，其余同实施例1中的步骤1。

本实施例所选择的特征响应为厂房结构的模态位移，实验测试得到模态位移为f＝5.34mm。

步骤2、本实施例除特征频率f_FEM为模态位移外，其余同实施例1的步骤2。

步骤3、对k＝25个样本点进行归一化处理，建立本实施例的二次响应面函数g(x)

$g\left(x\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{[\mathrm{\γ}{f}_{\mathrm{FEM}}\left(i\right)-f_{\mathrm{cs}}\left(i\right)]}^2={\mathrm{\β}}_0+\underset{i=1}{\overset{k^{\′}}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_i+\underset{i=1}{\overset{k^{\′}}{\mathrm{\Σ}}}{\mathrm{\β}}_{\mathrm{ii}}{x}_i^2---\left(1\right)$

式（1）中：除f_FEM为利用有限元软件计算得到的模态位移和f_cs为实验测试得到的模态位移外，其余同实施例的式（1）。

步骤4、根据式（1），设定初始加权系数γ＝1.0，得到初始二次响应面函数 $g_0\left(x\right)=0.2x_1^2+2.12x_2^2+4.50x_3^2-1.25x_4^2+1.3x_5^2+2.0x_6^2+3.2x_1+2.34x_2+5.4x_3-2x_4+9.1x_5-{3.4x}_6+10.1.$

本实施例对初始二次响应面函数g₀(x)进行自身特性检验，其结果满足旋转性、斜率可旋转性和稳健性的三个条件：

①可旋转性：响应量的方差仅是从设计中心到设计点距离的函数，与方向无关；

②斜率可旋转性：与设计中心的等距离点上的响应面的偏导数的方差是一个常数；

③稳健性：在个别数据有误或误差非正态分布的条件下，k个样本点中所含误差的样本点与不含误差的样本点对响应面结果引起的偏差在10%以内。

本实施例同时对初始二次响应面函数g₀(x)进行精度检验：

$\mathrm{RMSE}=\sqrt{\frac{\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{[f_{\mathrm{FEM}}\left(i\right)-f_{\mathrm{mod}}\left(i\right)]}^2}{k-k^{\′}}}---\left(2\right)$

式（2）中：

k为样本点的个数，本实施例中k＝25；

k′为设计变量的个数，本实施例中k′=6；

f_FEM为步骤1所利用有限元软件计算得到的模态位移；

f_mod为经过模型修正后的模态位移。

计算得到误差平方和RMSE＝0.03，即误差平方和RMSE→0。

$R^2=\frac{\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{[f_{\mathrm{mod}}\left(i\right)-{\stackrel{-}{f}}_{\mathrm{FEM}}\left(i\right)]}^2}{\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{[f_{\mathrm{FEM}}\left(i\right)-{\stackrel{-}{f}}_{\mathrm{FEM}}\left(i\right)]}^2}---\left(3\right)$

式（3）中：

为利用有限元软件计算得到的模态位移的平均值。

计算得到响应面精度R²＝0.95，即R²→1。

即本实施例同时满足：①误差平方和RMSE→0和②响应面精度R²→1。进入步骤5。

步骤5、对满足要求的二次响应面函数g(x)进行优化计算，得到二次响应面函数g(x)最小值ming(x)＝0.4所对应的设计变量x^*，再利用ANSYS有限元软件进行计算，得到对应的二次响应面函数g(x^*)＝0.5。

根据

${\mathrm{\γ}}_j=\frac{\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{f}_{\mathrm{FEM}}\left(i\right)\·f_{\mathrm{cs}}\left(i\right)}{{\left[f_{\mathrm{FEM}}\left(i\right)\right]}^2}---\left(4\right)$

式（4）中：

k为样本点的个数，本实施例中k＝25；

f_FEM为步骤1所利用有限元软件计算得到的模态位移；

f_cs为步骤1实验测试得到的模态位移。

得到迭代过程中二次响应面函数g(x)迭代第2步的加权系数γ₂＝0.979。

根据

$\left\{\begin{array}{c}|g\left(x^{*}\right)-\min g\left(x\right)|\≤{\mathrm{\ϵ}}_1\\ |{\mathrm{\γ}}_{j-1}-{\mathrm{\γ}}_j|\≤{\mathrm{\ϵ}}_2\end{array}\right.---\left(5\right)$

式（5）中：

g(x^*)是以x^*为设计变量代入有限元软件计算得到的二次响应面函数g(x)的值；

ming(x)为优化二次响应面函数g(x)后得到的最小值；

ε₁和ε₂为收敛阈值，本实施例中，阈值取ε₁＝ε₂＝0.1。

将g(x^*)＝0.5、ming(x)＝0.4、|γ_j-1-γ_j|=|0.979-1.0|代入式（5），得到

$\left\{\begin{array}{c}|0.5-0.4|\≤0.1\\ |0.979-1.0|\≤0.1\end{array}\right.$

即满足式（5）的收敛条件，停止计算。进入步骤7。

步骤7、经过2步迭代，反归一化处理后，输出优化后的二次响应面函数g(x)所得到的设计变量x^*分别为： $x_1^{*}=4.97e7(N/m);$ $x_2^{*}=4.74e7(N/m);$

（无量纲）； $x_4^{*}=5000$ $\left(\mathrm{Kg}\right);$ $x_5^{*}=12.85\left(m\right);$ $x_6^{*}=3.56\left(m\right).$

步骤8、建立模型修正后的评价指标

$\mathrm{ER}(f_{\mathrm{cs}},f_{\mathrm{mod}})=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}\left|\frac{f_{\mathrm{cs}}\left(i\right)-f_{\mathrm{mod}}\left(i\right)}{f_{\mathrm{cs}}\left(i\right)}\right|\×100\%---\left(6\right)$

式（6）中：f_cs为实验测试得到的模态位移；

f_mod为经过模型修正后的模态位移；修正后的模态位移为f＝5.21mm；

k为样本点的个数，本实施例k＝25。

计算得到修正后的ER值1.34%小于修正前的ER值15.21%，修正后的模态位移更接近实测模态位移，可为工程设计采用。

实施例3

一种用于抽水蓄能电站地下厂房结构的模型修正方法。其具体步骤是：

步骤1、同实施例1中的步骤1。

步骤2、本实施例除试验设计为正交试验设计和有限元软件采用ABAQUS有限元软件外，其余同实施例1的步骤2。

步骤3、同实施例1中的步骤3。

步骤4、根据步骤3的式（1），设定初始加权系数γ＝1.0，得到初始二次响应面函数 $g_0\left(x\right)=10.1x_1^2+9.56x_2^2-20.1x_3^2+3.26x_4^2+2.33x_5^2+8.54x_6^2-2.1x_1+5.12x_2+2.9x_3-12x_4+23.1x_5+{11.2x}_6-8.7.$

本实施例对初始二次响应面函数g₀(x)进行自身特性检验，其结果满足旋转性、斜率可旋转性和稳健性的三个条件：

①可旋转性：响应量的方差仅是从设计中心到设计点距离的函数，与方向无关；

②斜率可旋转性：与设计中心的等距离点上的响应面的偏导数的方差是一个常数；

③稳健性：在个别数据有误或误差非正态分布的条件下，k个样本点中所含误差的样本点与不含误差的样本点对响应面结果引起的偏差在10%以内。

本实施例同时对初始二次响应面函数g₀(x)进行精度检验：

根据

$\mathrm{RMSE}=\sqrt{\frac{\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{[f_{\mathrm{FEM}}\left(i\right)-f_{\mathrm{mod}}\left(i\right)]}^2}{k-k^{\′}}}---\left(2\right)$

式（2）中：

k为样本点的个数，本实施例中k＝25；

k′为设计变量的个数，本实施例中k′=6；

f_FEM为步骤1所利用有限元软件计算得到的前五阶特征频率；

f_mod为经过模型修正后的前五阶特征频率；

经计算，误差平方和RMSE＝0.51，并不趋近于0。

根据

$R^2=\frac{\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{[f_{\mathrm{mod}}\left(i\right)-{\stackrel{-}{f}}_{\mathrm{FEM}}\left(i\right)]}^2}{\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{[f_{\mathrm{FEM}}\left(i\right)-{\stackrel{-}{f}}_{\mathrm{FEM}}\left(i\right)]}^2}---\left(3\right)$

式（3）中：

为利用有限元软件计算得到的前五阶特征频率的平均值。

计算得到响应面精度R²＝0.65，并不趋近于1。

即本实施例不满足：①误差平方和RMSE→0和②响应面精度R²→1。返回步骤2。

步骤4.2、重新定义25个样本点，其余同实施例1的步骤2。

步骤4.3、同实施例1的步骤3。

步骤4.4、按照实施例1的步骤4，得到的结果满足旋转性、斜率可旋转性和稳健性的三个条件，同时满足①误差平方和RMSE＝0.05→0和②响应面精度R²＝0.96→1。进入步骤5。

步骤5、对满足要求的二次响应面函数g(x)进行优化计算，得到二次响应面函数g(x)最小值ming(x)＝0.2所对应的设计变量x^*，再利用ANSYS有限元软件进行计算，得到对应的二次响应面函数g(x^*)＝0.3。

根据

${\mathrm{\γ}}_j=\frac{\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{f}_{\mathrm{FEM}}\left(i\right)\·f_{\mathrm{cs}}\left(i\right)}{{\left[f_{\mathrm{FEM}}\left(i\right)\right]}^2}---\left(4\right)$

式（4）中：

k为样本点的个数，本实施例中k＝25；

f_FEM为利用有限元软件计算得到的前五阶特征频率；

f_cs为实验测试得到的前五阶特征频率。

得到迭代过程中二次响应面函数g(x)迭代第2步的加权系数γ₂＝0.982。

根据

$\left\{\begin{array}{c}|g\left(x^{*}\right)-\min g\left(x\right)|\≤{\mathrm{\ϵ}}_1\\ |{\mathrm{\γ}}_{j-1}-{\mathrm{\γ}}_j|\≤{\mathrm{\ϵ}}_2\end{array}\right.---\left(5\right)$

式（5）中：

g(x^*)是以x^*为设计变量代入有限元软件计算得到的二次响应面函数g(x)的值；

ming(x)为优化二次响应面函数g(x)后得到的最小值；

ε₁和ε₂为收敛阈值，本实施例中，阈值取ε₁＝ε₂＝0.1。

将g(x^*)＝0.3、ming(x)＝0.2、|γ_j-1-γ_j|=|0.982-1.0|代入式（5），得到

$\left\{\begin{array}{c}|0.3-0.2|\≤0.1\\ |0.982-1.0|\≤0.1\end{array}\right.$

即满足式（5）的收敛条件，停止计算，进入步骤7。

步骤7、经过2步迭代，反归一化处理后，输出优化后的二次响应面函数g(x)所得到的设计变量x^*分别为： $x_1^{*}=4.96e7(N/m);$ $x_2^{*}=2.53e7(N/m);$

（无量纲）； $x_4^{*}=3009$ $\left(\mathrm{Kg}\right);$ $x_5^{*}=13.16\left(m\right);$ $x_6^{*}=2.26\left(m\right).$

步骤8、建立模型修正后的评价指标

$\mathrm{ER}(f_{\mathrm{cs}},f_{\mathrm{mod}})=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}\left|\frac{f_{\mathrm{cs}}\left(i\right)-f_{\mathrm{mod}}\left(i\right)}{f_{\mathrm{cs}}\left(i\right)}\right|\×100\%---\left(6\right)$

式（6）中：f_cs为实验测试得到的前五阶特征频率；

f_mod为经过模型修正后的前五阶特征频率，修正后的前五阶特征频率为f₁＝13.123Hz，f₂＝14.560Hz，f₃＝17.453Hz，f₄＝19.678Hz，f₅＝21.564Hz；

k为样本点的个数，本实施例k＝25。

计算得到修正后的ER值12.56%大于10%，，则表明该模型修正方法的结果与真实情况相差较大，则返回步骤2。

步骤8.2、再重新定义25个样本点，其余同实施例1的步骤2。

步骤8.3、同实施例1的步骤3。

步骤8.4、按照实施例1的步骤4，得到的结果满足旋转性、斜率可旋转性和稳健性的三个条件，同时满足①误差平方和RMSE＝0.03→0和②响应面精度R²＝0.95→1。进入步骤5。

步骤8.5、对满足要求的二次响应面函数g(x)进行优化计算，得到二次响应面函数g(x)最小值ming(x)＝0.4所对应的设计变量x^*，再利用ANSYS有限元软件进行计算，得到对应的二次响应面函数g(x^*)＝0.48。

根据

${\mathrm{\γ}}_j=\frac{\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{f}_{\mathrm{FEM}}\left(i\right)\·f_{\mathrm{cs}}\left(i\right)}{{\left[f_{\mathrm{FEM}}\left(i\right)\right]}^2}---\left(4\right)$

式（4）中：

k为样本点的个数，本实施例中k＝25；

f_FEM为步骤1所利用有限元软件计算得到的前五阶特征频率；

f_cs为步骤1实验测试得到的前五阶特征频率。

得到迭代过程中二次响应面函数g(x)迭代第2步的加权系数γ₂＝0.976。

根据

$\left\{\begin{array}{c}|g\left(x^{*}\right)-\min g\left(x\right)|\≤{\mathrm{\ϵ}}_1\\ |{\mathrm{\γ}}_{j-1}-{\mathrm{\γ}}_j|\≤{\mathrm{\ϵ}}_2\end{array}\right.---\left(5\right)$

式（5）中：

g(x^*)是以x^*为设计变量代入有限元软件计算得到的二次响应面函数g(x)的值；

ming(x)为优化二次响应面函数g(x)后得到的最小值；

ε₁和ε₂为收敛阈值，本实施例中，阈值取ε₁＝ε₂＝0.1。

将g(x^*)＝0.48、ming(x)＝0.4、|γ_j-1-γ_j|=|0.976-1.0|代入式（5），得到

$\left\{\begin{array}{c}|0.48-0.4|\≤0.1\\ |0.976-1.0|\≤0.1\end{array}\right.$

即本实施例满足式（5）的收敛条件，停止计算。进入步骤8.7。

步骤8.7、经过2步迭代，反归一化处理后，输出优化后的二次响应面函数g(x)所得到的设计变量x^*分别为： $x_1^{*}=4.50e7(N/m);$ $x_2^{*}=2.21e7(N/m);$

（无量纲）； $x_4^{*}=3236$ $\left(\mathrm{Kg}\right);$ $x_5^{*}=10.12\left(m\right);$ $x_6^{*}=5.34\left(m\right).$

步骤8.8、建立模型修正后的评价指标

$\mathrm{ER}(f_{\mathrm{cs}},f_{\mathrm{mod}})=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}\left|\frac{f_{\mathrm{cs}}\left(i\right)-f_{\mathrm{mod}}\left(i\right)}{f_{\mathrm{cs}}\left(i\right)}\right|\×100\%---\left(6\right)$

式（6）中：f_cs为实验测试得到的前五阶特征频率；

f_mod为经过模型修正后的前五阶特征频率，修正后的前五阶特征频率为f₁＝12.432Hz，f₂＝14.123Hz，f₃＝16.546Hz，f₄＝18.213Hz，f₅＝20.213Hz；

k为样本点的个数，本实施例k＝25。

计算得到修正后的ER值2.11%小于修正前的ER值25.65%，修正后的频率更接近实测频率，为工程设计采用。

本具体实施方式与现有技术相比有如下优点：

第一、传统的模型修正技术是直接对有限元矩阵元素进行修正，这样会导致矩阵元素的对称性、正定性和稀疏性的消逝而难以得到实际应用。本具体实施方式采用概率方法，在修正过程中把实测数据和有限元计算数据都视为概率统计中的样本，本具体实施方式为模型修正这个反问题的求解提供了一个新的模式。

第二、数值仿真过程中会因结构简化和参数误差等因素引起误差，实验测试也会因外界干扰和设备仪器故障等因素带来误差，两者的误差来源完全不相同，因此直接利用实验测试结果对有限元模型进行修正带有盲目性，而本具体实施方式的响应面中采用加权系数用以考虑实测信息和有限元计算信息之间的信息融合，通过响应面自身特性检验和精度检验能保证计算精度，使模型修正更接近真实，修正效果好。

第三、采用响应面模型代替有限元模型进行模型修正迭代计算，能显著提高计算效率，尤其是根据模型修正流程通过二次开发，形成基于计算机语言和有限元APDL命令的模型修正软件，进一步提高计算效率。

因此，本具体实施方式具有计算效率高和修正效果好的特点。本发明还适用于其它大型厂房结构的模型修正。

Claims (7)

1.一种用于抽水蓄能电站地下厂房结构的模型修正方法，其特征在于：该模型修正方法的具体步骤是：

步骤1：根据抽水蓄能电站地下厂房结构的先验信息，确定该地下厂房结构需要修正的设计变量x，通过实验测试得到特征响应f_cs；

步骤2：根据抽水蓄能电站地下厂房结构的特点，在先验信息的指导下进行试验设计，得到k个样本点，再采用ANSYS有限元软件或ABAQUS有限元软件对厂房结构进行计算得到特征响应f_FEM；

步骤3：对k个样本点进行归一化处理，建立二次响应面函数g(x)

$g\left(x\right)=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{[\mathrm{\γ}{f}_{\mathrm{FEM}}\left(i\right)-f_{\mathrm{cs}}\left(i\right)]}^2={\mathrm{\β}}_0+\underset{i=1}{\overset{k^{\′}}{\mathrm{\Σ}}}{\mathrm{\β}}_i{x}_i+\underset{i=1}{\overset{k^{\′}}{\mathrm{\Σ}}}{\mathrm{\β}}_{\mathrm{ii}}{x}_i^2---\left(1\right)$

式（1）中：

x为厂房结构需要修正的设计变量，

f_FEM为利用有限元软件计算得到的特征响应，

f_cs为实验测试得到的特征响应，

γ为二次响应面函数g(x)的加权系数，为采用有限元软件计算得到的信息与实验测试信息的融合程度，

k为样本点的个数，

k′为设计变量x的个数，

β₀为采用最小二乘法所确定的二次响应面函数g(x)中的常数项，

β_i为采用最小二乘法所确定的二次响应面函数g(x)中的一次项系数，

β_ii为采用最小二乘法所确定的二次响应面函数g(x)中的二次项系数；

步骤4：根据式（1），设定初始加权系数γ＝1.0，得到初始二次响应面函数g₀(x)；对初始二次响应面函数g₀(x)进行自身特性检验和精度检验，若不能满足全部条件，则返回步骤2，重新定义k个样本点至满足自身特性检验和精度检验的全部条件；若能满足自身特性检验和精度检验的全部条件，则进入步骤5；

步骤5：对满足要求的二次响应面函数g(x)进行优化计算，得到二次响应面函数g(x)最小值ming(x)所对应的一组设计变量x^*，再利用ANSYS有限元软件或ABAQUS有限元软件进行计算，得到对应的二次响应面函数g(x^*)；

迭代过程中二次响应面函数g(x)迭代第j步的加权系数γ_j

${\mathrm{\γ}}_j=\frac{\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{f}_{\mathrm{FEM}}\left(i\right)\·f_{\mathrm{cs}}\left(i\right)}{{\left[f_{\mathrm{FEM}}\left(i\right)\right]}^2}---\left(2\right)$

式（2）中：

k为样本点的个数，

f_FEM为利用有限元软件计算得到的特征响应，

f_cs为实验测试得到的特征响应；

若二次响应面函数g(x)在迭代过程中满足以下收敛条件：

$\left\{\begin{array}{c}|g\left(x^{*}\right)-\min g\left(x\right)|\≤{\mathrm{\ϵ}}_1\\ |{\mathrm{\γ}}_{j-1}-{\mathrm{\γ}}_j|\≤{\mathrm{\ϵ}}_2\end{array}\right.---\left(3\right)$

式（3）中：

g(x^*)是以x^*为设计变量代入有限元软件计算得到的二次响应面函数g(x)的值，

ming(x)为优化二次响应面函数g(x)后得到的最小值，

ε₁和ε₂为收敛阈值；

则设计变量x^*收敛，停止计算，进入步骤7，否则进入步骤6；

步骤6：按照式（2）更新加权系数γ_j，再返回步骤3，构造新的二次响应面函数g′(x)，重复步骤4和步骤5至满足式（3）的收敛条件；

步骤7：经过j步迭代，反归一化处理后，输出优化后的二次响应面函数g(x)所得到的设计变量

步骤8：建立模型修正后的评价指标

$\mathrm{ER}(f_{\mathrm{cs}},f_{\mathrm{mod}})=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}\left|\frac{f_{\mathrm{cs}}\left(i\right)-f_{\mathrm{mod}}\left(i\right)}{f_{\mathrm{cs}}\left(i\right)}\right|\×100\%---\left(4\right)$

式（4）中：f_cs为实验测试得到的特征响应，

f_mod为经过模型修正后的特征响应，

k为样本点的个数；

若ER(f_cs，f_mod)≤10%时，则表明该模型修正方法的结果接近于真实情况；若ER(f_cs，f_mod)＞10%时，则返回步骤2，直至ER(f_cs，f_mod)≤10%。

2.根据权利要求1所述的用于抽水蓄能电站地下厂房结构的模型修正方法，其特征在于所述的设计变量x是利用先验信息得到的设计变量，包括上下游刚度值、左侧粘结力值、约束分布密度值、质量点大小值和质量点坐标值。

3.根据权利要求1或2所述的用于抽水蓄能电站地下厂房结构的模型修正方法，其特征在于所述的先验信息是：

①刚度悬殊：主体结构如蜗壳、机墩和风罩等刚度较大，次要结构如板梁柱等刚度较小；

②结构类型固定：厂房结构为3～4层，为厚板结构或为板梁结构；

③楼板开孔较大：楼板上开有吊物孔和楼梯孔，振型往往围绕开孔处发生；

④约束条件：上下游围岩作用以及左右侧厂房之间的相互作用力对结构的特征频率影响十分明显；

⑤设计规范要求：在激励荷载频率一定的条件下，地下厂房结构的特征频率要满足

4.根据权利要求1所述的用于抽水蓄能电站地下厂房结构的模型修正方法，其特征在于所述的特征响应是地下厂房结构的特征频率或模态位移。

5.根据权利要求1所述的用于抽水蓄能电站地下厂房结构的模型修正方法，其特征在于所述的试验设计为均匀设计和正交设计中的一种。

6.根据权利要求1所述的用于抽水蓄能电站地下厂房结构的模型修正方法，其特征在于所述的自身特性检验是指同时满足如下三个条件是：

①可旋转性：响应量的方差仅是从设计中心到设计点距离的函数，与方向无关；

②斜率可旋转性：与设计中心的等距离点上的响应面的偏导数的方差是一个常数；

③稳健性：在个别数据有误或误差非正态分布的条件下，k个样本点中所含误差的样本点与不含误差的样本点对响应面结果引起的偏差在10%以内。

7.根据权利要求1所述的用于抽水蓄能电站地下厂房结构的模型修正方法，其特征在于所述的精度检验是指同时满足如下两个条件：

①误差平方和RMSE→0

$\mathrm{RMSE}=\sqrt{\frac{\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{[f_{\mathrm{FEM}}\left(i\right)-f_{\mathrm{mod}}\left(i\right)]}^2}{k-k^{\′}}}---\left(5\right)$

式（5）中：

k为样本点的个数，

k′为设计变量的个数，

f_FEM为利用有限元软件计算得到的特征响应，

f_mod为经过模型修正后的特征响应；

②响应面精度R²→1

$R^2=\frac{\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{[f_{\mathrm{mod}}\left(i\right)-{\stackrel{-}{f}}_{\mathrm{FEM}}\left(i\right)]}^2}{\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{[f_{\mathrm{FEM}}\left(i\right)-{\stackrel{-}{f}}_{\mathrm{FEM}}\left(i\right)]}^2}---\left(6\right)$

式（6）中：

为利用有限元软件计算得到的特征响应的平均值。

Images