基于Isight-MSC.Patran/Nastran联合仿真的弧形闸门动力模型修正

吴思远，王正中，王 岳，徐 超，刘计良，吴守军

(1.同济大学 土木工程防灾国家重点实验室，上海 200092；2.西北农林科技大学 旱区寒区水工程安全研究中心，陕西 杨凌 712100；3.西安理工大学 水利水电学院，西安 710048)

弧形闸门作为水库大坝的调蓄建筑物，在完成预定复杂工况(静、动水各开度启闭)的同时，要保证结构自身稳定，以满足整体枢纽全生命周期安全运行的要求。通常弧形闸门可简化为主梁支臂框架体系，其属于中柔度结构形式[1]，加之启闭杆及周边止水对其结构刚度约束有限[2]，在运行过程中，结构失稳破坏问题普遍存在[3]。以往弧形闸门按规范要求，首先应满足静力稳定，在此基础上增加动力安全系数[4]，此类设计方法，存在局部强度雍余，不满足节能降耗的发展要求，同时忽略了实际工程中所存在的复杂流激振动现象，如强迫振动、自激振动和参数振动等耦合工况[5]。众多调查结果表明不良振动是闸门破坏的主要因素[6]。因此为了提高弧形闸门设计水平，迫切需要从动力学角度进行结构优化及减振控制研究。而结构动力分析的前提是建立一个合理高效的分析模型。

传统弧形闸门有限元分析基于原型数据直接建模，而实际结构中，阻尼的不确定性，虚杆、虚板的影响及长期运行造成的锈蚀缺陷等问题，凸显了有限元建模过程的理想化，而与实际运行工况存在一定差距，导致实践中危险估计不足，威胁闸门乃至大坝整体运行[7]。而工程人员在进行有限元模型误差分析时往往仅通过经验判断进而改良，缺乏理论依据。同时文献[8]从弧形闸门建模方式的角度充分说明了除实际结构中所存在的随机误差因素外，建模方式的不同同样会导致结构力学特性的差异。由此可知弧形闸门实测数据与有限元模拟结果难以统一，给进一步结构振动特性分析及振动控制的研究带来困难，有必要对模型进行优化修正。

现有动力模型修正理论主要包括矩阵型修正法、参数型修正法和基于神经网络的模型修正方法[9]，在航空及精密制造领域，以往所使用的矩阵型修正法需投入巨大精力，通过精确测量来获取全模态信息，而在神经网络的模型修正方法应用过程中，又需反复试算选取网络模型及隐含神经元，上述两种方法难以与大型复杂工程实践相结合。随着有限元软件的发展，基于有限元软件优化模块的参数型辅助修正方法被广泛接受。胡仔溪[10]、苏忠亭等[11]和李志刚等[12]分别基于Nastran和ANSYS有限元软件，以结构参数(几何尺寸、弹性模量、密度)为优化变量，对模型进行修正，以此减小结构动力特性实测与仿真间的误差，梁忠仔等[13]在此基础上集成多个有限元软件，发挥各自优势，提高了优化效率。上述文献均从有限元模型整体优化的角度分析，以期能提高动力模型精度，但优化变量的选取不合理，其结果改变了真实结构的几何特性和材料属性，存在修正过度导致模型失真的问题。同时优化过程局限于某一软件自带功能，没能充分发挥各软件特色。针对优化变量选取问题，杨世浩等[14]提出利用节点集中质量修正，得到了弧形闸门简化动力力学模型，但文中并未提及修正方法，不具有广泛性。王轲基于节点质量重分布，提出了一套详尽的静动力模型转换方法[15]，并给出了相应的模型评估方法[16]，但任存在应用过程中简单通用、快速优化收敛的问题。

为解决上述问题，本文提出一种基于Isight优化平台的MSC.Patran/Nastran联合仿真节点质量特征灵敏度参数优化快速建模方法。Isight是一个将设计过程集成化、自动化的数字化平台，整合各类仿真软件及设计流程，并通过内嵌优化算法，实现简易高效优化设计。

1 特征灵敏度优化分析

特征灵敏度优化是一种参数型优化修正理论，基本思路是以修正模型与实验模型之间在同一激励下动力特性的误差为目标函数，选择一定的修正量使该误差满足最小化，以此来达到优化目的。灵敏度从数学角度表述为：若一函数F(x)可导，对任意变量xJ的一阶灵敏度可表示为

(1)

在弧形闸门振动中，变量xJ可理解为动力特性参数(自振频率和振型相对应的特征值λ和特征向量φ)对结构参数Pl(质量、刚度、阻尼)的改变率，即所谓的特征值灵敏度∂λ/Pl和特征向量灵敏度∂φ/Pl。对于弧形闸门有限元模型具有N个自由度的黏性阻尼系统，其自由振动运动方程为

(2)

式(2)中的模态质量特征值灵敏度可展开为

(3)

式中:φik为第k阶模态中第i个节点振型值；φjk为第k阶模态中第j个节点振型值。

特征向量灵敏度表示为特性向量线性组合的形式

(4)

φs表示第s阶模态任意节点处振型值。

在上述特征值灵敏度分析的基础上，通过泰勒级数展开可得到得到自振频率与振型修正量为

(5)

(6)

同理，可得到刚度和阻尼特征灵敏度修正量。值得注意的是，特征灵敏度优化方法在不同有限元软件中所体现的功能具有一定的差别，在ANSYS中，通过概率设计系统(PDS)模块获取灵敏度[17]，其实质是Spearman秩相关系数[18]，仅是从概率角度分析不同物理参数(密度、模量、尺寸等)对动力响应特性的影响。而MSC.Nastran提供了设计灵敏度及优化模块(Design Sensitivity Analysis and Optimization)[19]，以此得到由设计变量改变引起的结构刚度和质量矩阵的偏导，即特征值灵敏度和特征向量灵敏度相对应的设计灵敏度系数，可用于定量分析。

2 弧形闸门优化动力模型设计

本文将结构固有频率与动力响应作为约束条件，以各节点质量修正量作为优化变量，保证了结构质量、质心、转动惯量和刚度与实际结构的一致性。在此基础上运用特征灵敏度优化设计方法，通过Isight-MSC.Patran/Nastran联合仿真确定设计变量的取值，经过多次迭代使模型与实际结构的动力特性差异趋于最小。

2.1 优化模型

2.1.1 设计变量

以往文献[10-12]中模型修正基于体型尺寸优化，选取众多变量，尽管能达到有限元模拟的动力特性与实测结果相符，但模型与实测结构的几何属性大相径庭，结构刚度也随之改变，造成结构承载力变化，影响结构稳定分析。其次众多改变量的选取在实质上弱化了修正概念，优化后模型丧失了原有的对称性及原模型各单元间连接信息。针对上述问题本文仅以弧形闸门节点处质量的修正量Δm作为设计变量X，可表示为

X=(Δm1,Δm2,…,Δmn)

(7)

式中：n为有限元模型节点数。

2.1.2 目标函数

为满足模型结构各节点相关信息偏差最小的目的，以各节点质量修正量平方和最小为目标函数

(8)

2.1.3 等式约束条件

修正后结构重要几何特征质量和质心应严格满足静力学约束条件。

(1) 修正后各节点集中质量之和等于总质量m

(9)

式中：mAi为初步建立的有限元质量矩阵MA中i节点集中质量。

(10)

式中：xi，yi，zi为有限元节点坐标。

2.1.4 不等式约束条件

为保证自变量约束范围内存在有效解，适当放宽约束条件，加入容许误差ri。

(1) 各节点集中质量的转动惯量和与结构的总转动惯量(Ix，Iy，Iz)约束条件如下

(11)

(2) 自振频率约束条件为

(12)

(3) 在有限元软件分析中，模态振型是一个无量纲相对坐标值，本文中取正则振型，振型约束条件中第k阶模态振型可表述为

(13)

2.2 优化方案

2.2.1 优化算法

上述优化模型属于非线性规划，其中目标函数是非线性，而等式约束和不等式约束是线性，可以转化为二次规划(QP)问题，标准形式为

(14)

H为n×n阶对称矩阵；Gi(X)为向量函数等式约束与不等式约束在X处的值；Vlb和Vub分别为X的上下限。若优化变量存在全局最优解，式(14)应满足K-T条件，即要求H为正定矩阵，此时目标函数为凸函数，优化问题为严格二次规划。

根据本文优化目标，转化成标准形式

(15)

H正定，故本文优化结果存在全局最优解。本文基于Isight采用序列联系规划法(DONLP)，此方法是对拉式函数H矩阵的Pantoja-Mayne更新，其对armiji-type算法改进，对变量上下界采用类梯度法处理。

2.2.2 Isight-MSC.Patran/Nastran联合仿真

首先在Patran中建立弧形闸门有限元几何模型，赋予材料属性得到初始模型，获取模型几何特征，即总质量、质心和惯矩，并采用集中质量法进行模态分析提取质量矩阵(每个元素作为优化变量对应模型各节点质量)。然后把获取的质量矩阵离散至不含密度属性的有限元模型上，得到集中质量分布模型，再通过模态分析获取集中质量分布模型的自振频率和正则振型，利用灵敏度分析模块获取各节点质量对应的频率和振型响应的灵敏度。利用Isight提取的相关信息进行优化计算，从而得到质量修正量Δmi，通过修改bdf文件，导入Nastran中进行动力分析，获取自振频率并预设收敛准则，迭代误差依据实测与初始集中质量有限元模型自振频率相对误差大小决定。循环上述过程即可得到准确的弧形闸门动力有限元模型，并获取自振频率与振型数据，通过模态置信度准则MAC和频率相对误差进行误差判断，控制流程图如图1所示。

3 方法验证

某高水头弧形闸门的主框架有限元模型如图2(a)所示(水电站机电设计手册编写组1988)，图中圆圈的编号为节点号，有限元模型共有16个单元，17个节点，每个节点有3个自由度。构件截面尺寸如图2(b)所示，主梁采用工字形截面，弯曲平面内的惯性矩Ix1=1.235×10-2m4，单位长度质量m=430 kg/m；支臂采用箱形截面，弯曲平面内的惯性矩Ix2=2.319×10-3m4，单位长度质量m=386 kg/m；弹性模量E=210 GPa[20]。以上述模型为无误差状态，在此基础上，将第3、7、12号节点集中质量+5%～+40%的摄动，4、8、14号节点集中质量-5～-40%的摄动，以模拟建模误差，同时保证总质量不变，具体参数如表1所示。

图1 Isight-MSC.Patran/Nastran联合仿真流程图Fig.1 The joint simulation flowchart of Isight and MSC.Patran/Nastran

(a) 主框架有限元模型(m)

(b) 构件截面尺寸(mm)

表1 节点摄动误差Tab.1 Node modeling error

利用本文方法对初始有限元模型进行修正，可得到相应的优化模型与目标频率相对误差迭代收敛曲线。

由图3计算结果可知，本文方法可识别并修正摄动误差，经五次迭代优化后有限元模型与目标频率最大误差小于1%。以此说明本文方法能够快速建立满足目标状态下振动特性的结构动力有限元模型，进一步验证了本文方法实际应用的可能性。

图3 误差迭代收敛曲线Fig.3 Identified error convergence curve

在上述验证基础上与文献[10-12]有限元多参数变量优化方法进行对比，优化结果及力学特性分析如表2。

表2 不同优化方法结果分析Tab.2 Calculation results of different optimization

通过与现有文献中多参数变量优化方法对比表明，在保证低阶频率修正相对误差小于0.5%的情况下，本文方法仅对节点质量进行修正，其优化模型不改变原有承载力，振型相关性较好，随机振动位移响应误差较小。而文献[10-12]通过灵敏度相关分析选取众多变量(节点质量、截面尺寸、弹性模量)，尽管满足了模态频率的一致性，但极大的改变了原有结构形态，造成结构承载力改变，同时文献中忽略振型约束条件，导致优化模型与原有结构振型相关性较差，造成了随机振动过程中，位移响应存在较大误差。

4 弧形闸门动力模型优化实例

以图4某水电站深孔闸门1∶20的“全水弹性相似”模型实验为例，实际最大挡水水头为93 m，孔口尺寸为4 m×5 m(宽×高)。钢结构密度为7 850 kg/m3，弹性模量为210 GPa，泊松比为0.25。实验仪表为：航天工业总公司701所研制的带压电式力传感器的力锤，扬州无线电二厂生产的CA-YD-141型压电式三向加速度传感器，2692型电荷放大器及东方振动与噪声研究所研制的INV206D(M)E型智能信息采集系统。

有限元模型均采用壳单元，划分单元数为10 661，节点数为9 216，模型如图5所示。同时考虑到闸门振动属流体诱发振动，根据Westergaard公式加入流固耦合引起的“附加质量”。本文迭代收敛准则为前五次迭代相对误差不超过10%，此后不超过5%，优化结果及评价如表3、图6所示。

图4 弧形闸门完全水弹性相似模型Fig.4 Physical model of radial gate

图5 弧形闸门有限元模型Fig.5 Finite element model of radial gate

表3 有限元模型计算频率与实测结果对比Tab.3 Frequencies comparison between the FEM and testing

图6 相关置信因子图Fig.6 Figure of MAC

从表3可以看出，初始模型的计算结果与实测频率存在较大误差，通过特征值灵敏度优化迭代后，误差满足收敛预设的误差范围。通常MAC对角元素大于0.8，非对角元素小于0.2可以说明振型相关性较好，从图6可以看出，实测振型与优化后弧形闸门动力模型相关度较好。弧形闸门流激振动过程可由振型叠加法得到，以此判断本文弧形闸门动力模型可合理模拟流激振动响应。

5 结 论

本文针对弧形闸门有限元模型动力分析过程中存在偏差的问题，提出一种Isight-MSC.Patran/Nastran联合仿真与特征灵敏度优化相结合的节点质量重分布快速建模方法。通过对弧形闸门主框架有限元模型的摄动检测，证明了本文方法的有效性，进而与以往参数优化方法进行对比说明了本文在尽可能保证与原有结构几何、材料性质一致情况下，能够获得更为精确的动力有限元分析模型。在此基础上根据弧形闸门完全水弹性振动实验数据，利用本文方法建立考虑流固耦合情况下的有限元模型，结果表明，该方法可以较好地用于修正模型动力仿真中存在的误差，为提高弧形闸门动力设计水平，确定其减振控制策略提供了更精确的动力分析模型。