CN103093026B - 液体附加质量振动反演方法 - Google Patents

Abstract

本发明提出一种流体附加质量振动反演算法，该方法不需将流体假定为理想流体，能够真实反映结构物周围流体的自身及流动特性，本发明利用实际结构的振动模态信息作为流体附加质量计算的输入数据，同时将结构的环境状态分为两种，即空气中及流体中，通过两种状态的振动测试，既可以考虑实际结构与有限元模型之间的模型误差，又可以在此基础上以矩阵的形式获得每一个处于流体中单元的附加质量，所计算的流体附加质量以矩阵的形式存在，便于工程计算、分析。

Description

液体附加质量振动反演方法

技术领域

本发明涉及一种利用结构振动加速度数据进行结构物周围流体附加质量计算的新方法，特别是针对海洋平台的流体附加质量算法。

背景技术

在船舶与海洋工程领域，外部海洋动力环境因素对深海工程结构作用，可以通过流固耦合分析确定深海平台系统的运动响应，并判断平台系统的稳定性，这是进行海洋油气开采所必需的一项关键性的科学问题。在深海环境中，非线性作用将大大加强，平台结构及其附属的柔性结构系统都将呈现出明显的非线性运动响应特征，平台整体会发生大尺度振荡和摇摆，立管、海底管道等会出现剧烈的涡激振动。为保证深海平台系统在一般状况下的稳定作业和极端状况下的安全避险，非线性流固耦合过程的准确模拟和分析具有极其重要的意义。

准确的流固耦合计算离不开符合实际情况的水动力参数，以及因结构振动而引起的附加水动力参数，如当流体及结构二者的密度在一个数量级时，附加质量的模拟准确与否将严重影响流固耦合分析结果。虽然API规范中也给出了常规截面的附加质量系数，但经验性因素占据相当大的成分。因为附加质量是一个比较复杂的问题，关于其定义以及求解，一直存在着争议。正如T.Sarpkaya所言：附加质量是流体动力学著名的但难于理解、易于混淆的特征量之一。和物体的质量一样，仅当物体做加速度运动的时候才会显示出附加质量的存在。它决定于细长海洋结构物或其周围的流体及尾流的运动形式、与其它物体的接近程度、自由液面和时间。Vikestad（2000）通过试验研究了涡激振动中附加质量的变化，试验结果表明附加质量会随着来流速度的变化出现较大的变化，甚至会出现负的附加质量。王艺(2006)从涡激振动结构的运动平衡方程入手，推导出结构出现涡激振动锁频阶段的附加质量的表达式；并根据试验结果分析附加质量的变化情况，进一步推导出附加质量的估算公式。周期平均的附加质量系数可能是负值。在论及流动相对光滑和粗糙柱体做正弦运动的问题时，Keulegan & Carpenter（1958）和Sarpkaya（1976）都曾讨论过类似的“负的”附加质量。之后，Vandiver（1993），Vikestad（2000）等也研究过。Vandiver（1993）曾说过，“负的附加质量不过反映了作用在柱体上的流体力的符号与加速度同相位”。考虑到附加质量系数是质量输运的总和在周期内（加速段和减速段）的平均值，Sarpkaya（2004）指出“附加质量为负值意味着减速过程的漂移质量要比加速过程中的大”。在计算流体与结构相互作用方面，Conca(1997)指出对于浸没在不可压缩流体中运动的物体，其附加质量并不依赖于粘滞系数，并且流体的粘滞性可用时间卷积阻尼项来模拟。Causin（2005）讨论了应用弱耦合时间推移算法计算流固耦合时的数值稳定性问题。其成果主要应用于人体动脉中血液的流固耦合问题，同时也可以应用到流体中弹性体的流固耦合计算，他们的算法中考虑了流体作用在结构上的附加质量效应。Miguel（2006）对具有强附加质量效应的流固耦合问题应用半隐式方法进行了研究。其算法在提高计算效率方面具有一定的优势。Degroote（2010）研究了流固耦合的分离模拟方法，并且发展了具有雅克比求逆的近似解准牛顿耦合算法。

常规附加质量的求解，一般是先通过理论求解结构所受到的流体作用力，通过将其中与加速度有关的惯性分量分离出，以获得流体的附加质量。并且，附加质量求解过程中往往假定流体为理想流体，即流体满足无旋、不可压缩等条件。而这些条件又与做大尺度振荡和摇摆的平台环境条件差别较大，导致所计算结果与实际情况差别较大，从而造成结构设计时计算误差。

发明内容

本发明提出一种流体附加质量振动反演的算法，该方法不需将流体假定为理想流体，能够真实反映结构物周围流体的自身及流动特性，其采用的技术方案如下：

一种液体附加质量振动反演算法，其特征在于包括如下步骤：

Ａ、建立结构有限元数值模型，获取该有限元模型的质量矩阵与刚度矩阵，分别标记为M、K，并存于专用存储器中；

Ｂ、在结构上布置加速度/速度/位移传感器；

Ｃ、将结构置于空气中获得未放流体状态下结构的模态参数并存于专用存储器中，其中j为模态阶次；

Ｄ、将结构置于流体中获得流体中结构的模态参数并存于专用存储器中；

Ｅ、流体附加质量计算：

e1、从专用存储器中读取上述步骤A中存储的数据，即M、K；

e2、结构模态参数提取

利用模态参数识别方法提取模态频率f、模态振型Φ，其中 $f_i=\sqrt{{\mathrm{\λ}}_i}\×180/\left(2\mathrm{\π}\right);$

e3、建立未放流体状态结构自由振动方程；

e4、建立流体中结构自由振动方程；

e5、将步骤e3、e4中结构自由振动方程转换到模态领域后求得液体附加质量。

进一步地，所述步骤e3中，未放流体状态结构自由振动方程为其中M′、K′分别为实际结构的质量、刚度矩阵；x为结构位移，为加速度，其中未放流体状态结构与步骤A结构中M、K的关系为：

$M^{\′}=M+\underset{\mathrm{nv}=1}{\overset{\mathrm{nv}=\mathrm{Nmv}}{\mathrm{\Σ}}}{\mathrm{\η}}_{\mathrm{nv}}{M}_n+\underset{\mathrm{nw}=1}{\overset{\mathrm{nw}=\mathrm{Nmw}}{\mathrm{\Σ}}}{\mathrm{\η}}_{\mathrm{nw}}{M}_n$

$K^{\′}=K+\underset{n=1}{\overset{n=\mathrm{Nk}}{\mathrm{\Σ}}}{\mathrm{\κ}}_n{K}_n$

其中M_n、K_n为有限元模型中第n个单元质量、刚度矩阵在整体坐标系下的形式；κ_n为单元总数为Nk时第n个单元的修正系数；η_nv为位于流体中的单元总数为Nmv时第nv个单元的修正系数；η_nw为位于流体中的单元总数为Nmw时第nw个单元的修正系数。

进一步地，所述步骤e4中，流体中结构自由振动方程为其中M″、K″分别为实际结构并考虑流体附加质量的质量、刚度矩阵，以M_a代表流体附加质量，流体中结构与空气中结构质量矩阵、刚度矩阵关系为：

$M^{\′\′}=M^{\′}+M_a=M+\underset{\mathrm{nv}=1}{\overset{\mathrm{nv}=\mathrm{Nmv}}{\mathrm{\Σ}}}{\mathrm{\η}}_{\mathrm{nv}}{M}_n+\underset{\mathrm{nw}=1}{\overset{\mathrm{nw}=\mathrm{Nmw}}{\mathrm{\Σ}}}{\mathrm{\η}}_{\mathrm{nw}}{M}_n+\underset{\mathrm{nw}=1}{\overset{\mathrm{nw}=\mathrm{Nmw}}{\mathrm{\Σ}}}{\mathrm{\β}}_{\mathrm{nw}}{M}_n$

K″＝K′

其中β_nw为单元总数为Nmw时位于流体中的第n个单元的质量附加系数。

进一步地，用有限元模型的第i阶模态振型Φ_i分别前乘步骤e5中得到的两种状态下模态领域振动方程，

${\mathrm{\Φ}}_i^T{K}^{\′}{\mathrm{\Φ}}_j^{\′}+{\mathrm{\λ}}_j^{\′2}{\mathrm{\Φ}}_i^T{M}^{\′}{\mathrm{\Φ}}_j^{\′}=0$ 与

${\mathrm{\Φ}}_i^T{K}^{\′\′}{\mathrm{\Φ}}_l^{\′\′}+{\mathrm{\λ}}_j^{\′\′2}{\mathrm{\Φ}}_i^T{M}^{\′\′}{\mathrm{\Φ}}_j^{\′\′}=0$

上述两式相加，得

其中j为结构处于空气中时测得的模态阶次，l为结构处于流体中时测得的模态阶次；

将步骤e3、e4带入上述第3式中得

${\mathrm{\Φ}}_i^T(K+\underset{n=1}{\overset{n=\mathrm{Nk}}{\mathrm{\Σ}}}{\mathrm{\κ}}_n{K}_n){\mathrm{\Φ}}_j^{\′}+{\mathrm{\λ}}_j^{\′2}{\mathrm{\Φ}}_i^T(m+\underset{\mathrm{nv}=1}{\overset{\mathrm{nv}=\mathrm{Nmv}}{\mathrm{\Σ}}}{\mathrm{\η}}_{\mathrm{nv}}{M}_n+\underset{\mathrm{nw}=1}{\overset{\mathrm{nw}=\mathrm{Nmw}}{\mathrm{\Σ}}}{\mathrm{\η}}_{\mathrm{nw}}{M}_n){\mathrm{\Φ}}_j^{\′}+$

${\mathrm{\Φ}}_i^T(k+\underset{n=1}{\overset{n=\mathrm{Nk}}{\mathrm{\Σ}}}{\mathrm{\κ}}_n{K}_n){\mathrm{\Φ}}_l^{\′\′}+{\mathrm{\λ}}_l^{\′\′2}{\mathrm{\Φ}}_i^T(M+\underset{\mathrm{nv}=1}{\overset{\mathrm{nv}=\mathrm{Nmv}}{\mathrm{\Σ}}}{\mathrm{\η}}_{\mathrm{nv}}{M}_n+\underset{\mathrm{nw}=1}{\overset{\mathrm{nw}=\mathrm{Nmw}}{\mathrm{\Σ}}}{\mathrm{\η}}_{\mathrm{nw}}{M}_n+\underset{\mathrm{nw}=1}{\overset{\mathrm{nw}=\mathrm{Nmw}}{\mathrm{\Σ}}}{\mathrm{\β}}_n{M}_n){\mathrm{\Φ}}_l^{\′\′}=0$

进一步化简为

$\underset{n=1}{\overset{n=\mathrm{Nk}}{\mathrm{\Σ}}}{\mathrm{\κ}}_n(A_1+A_2)+\underset{\mathrm{nv}=1}{\overset{\mathrm{Nmv}}{\mathrm{\Σ}}}{\mathrm{\η}}_{\mathrm{nv}}(B_1+B_2)+\underset{\mathrm{nw}=1}{\overset{n=\mathrm{Nmw}}{\mathrm{\Σ}}}{\mathrm{\η}}_{\mathrm{nw}}(C_1+C_2)+\underset{\mathrm{nw}=1}{\overset{N_{\mathrm{mw}}}{\mathrm{\Σ}}}{\mathrm{\β}}_{\mathrm{nw}}D=F$

其中 $A_1={\mathrm{\Φ}}_i^T{K}_n{\mathrm{\Φ}}_j^{\′},$ $A_2={\mathrm{\Φ}}_i^T{K}_n{\mathrm{\Φ}}_l^{\′\′},$ $B_1={\mathrm{\λ}}_j^{\′2}{\mathrm{\Φ}}_i^T{M}_n{\mathrm{\Φ}}_j^{\′},$ $B_2={\mathrm{\λ}}_l^{\′\′2}{\mathrm{\Φ}}_i^T{M}_n{\mathrm{\Φ}}_l^{\′\′},$ $C_1={\mathrm{\λ}}_j^{\′2}{\mathrm{\Φ}}_i^T{M}_n{\mathrm{\Φ}}_j^{\′},$ $C_2={\mathrm{\λ}}_l^{\′\′2}{\mathrm{\Φ}}_i^T{M}_n{\mathrm{\Φ}}_l^{\′\′},$ $D={\mathrm{\β}}_{\mathrm{nw}}{\mathrm{\λ}}_l^{\′\′2}{\mathrm{\Φ}}_i^T{M}_n{\mathrm{\Φ}}^{\′\′},$ $F=-({\mathrm{\Φ}}_i^{T}K{\mathrm{\Φ}}_j^{\′}+{\mathrm{\λ}}_j^{\′2}{\mathrm{\Φ}}_i^T{\mathrm{M\Φ}}_j^{\′}+{\mathrm{\Φ}}_i^{T}K{\mathrm{\Φ}}_l^{\′\′}+{\mathrm{\λ}}_l^{\′\′2}{\mathrm{\Φ}}_i^{T}M{\mathrm{\Φ}}_l^{\′\′})$

写成矩阵的形式

[A₁ A₂ B₁ B₂ C₁ C₂ D]·Δ＝F

进一步简化为

·Δ＝F

其中□＝[A₁ A₂ B₁ B₂ C₁ C₂ D]， $\mathrm{\Δ}=\left\{\begin{array}{c}\mathrm{\κ}\\ {\mathrm{\η}}_v\\ {\mathrm{\η}}_w\\ \mathrm{\β}\end{array}\right\}.$

应用线性规划方法求解

Δ＝lsqlin(□F，A，b，Aeq，Beq，Lb，Ub）

其中lsqlin为软件Matlab命令，并且满足A□Δ≤b，Aeq□Δ＝Beq，且Lb≤Δ≤Ub，A,b,Aeq,beq,Lb,Ub均为施加的约束条件，

流体附加质量提取，并记为：

${\tilde{M}}_a=\underset{\mathrm{nw}=1}{\overset{\mathrm{nw}=\mathrm{Nmw}}{\mathrm{\Σ}}}{\mathrm{\β}}_{\mathrm{nw}}{M}_n.$

进一步地，所述动力响应信号是加速度和/或是速度和/或是位移。

与现有技术相比，本发明的优点和积极效果是：

1、本发明利用实际结构的振动模态信息作为流体附加质量计算的输入数据，较现有方法而言，该算法既不需要对流体自身特性做任何假定，同时所计算流体附加质量亦能够充分反映结构周围流体的真实特性，计算结果也更符合工程实际；

2、本发明将结构的环境状态分为两种，即空气中及流体中。通过两种状态的振动测试，既可以考虑实际结构与有限元模型之间的模型误差，又可以在此基础上以矩阵的形式获得每一个处于流体中单元的附加质量矩阵。而传统方法以简单的流体附加系数考虑该部分流体附加质量，故而本发明计算结果更便于工程结构的设计、分析；

3、本发明的输入数据为结构的振动模态信息，该算法的优势在于本算法仅需要少数低阶实测模态。对于大型结构，借助动力测试数据获得该结构的所有模态是不可能实现的，尤其对于环境激励下结构的动力响应测试问题，因此，本发明可以进行环境激励下的海洋平台结构流体附加质量的计算问题，具有实际应用价值。

结合附图阅读本发明的具体实施方式后，本发明的其他特点和优点将变得更加清楚。

附图说明

图1为本发明的导管架式海洋平台结构有限元模型图；

图2为导管架式海洋平台刚度参数建模误差；

图3为导管架式海洋平台处于流体中单元的质量参数建模误差；

图4为导管架式海洋平台处于空气中单元的质量参数建模误差；

图5为导管架式海洋平台处于流体中单元流体附加系数；

图6为导管架式海洋平台处于流体中单元的计算质量系数与计算值对比；

图7为导管架式海洋平台处于空气中单元的计算质量系数与计算值对比；

图8为导管架式海洋平台处于流体中单元的流体附加系数与计算值对比。

具体实施方式

本发明的思想如下：浸没于流体中的结构之所以与其在空气中的频率、振型存在差异，其原因在很大程度上是因为结构周围的流体因其流速等变化对结构自身产生了附加水动力项，即附加质量及附加阻尼。考虑到无论该附加水动力项怎样随时间变化，其最终将在结构自身的振动响应信号中有所体现的特点，捕捉该信息并加以处理，即可反演流体附加质量的分布情况。基于该思想，本发明提出一种流体附加质量振动反演的算法。

本发明采用以下技术方案予以实现：

A建立结构有限元数值模型，获取该有限元模型的质量矩阵与刚度矩阵，分别标记为M、K并存储入专用存储器中；

b）在结构上布置加速度/速度/位移传感器，该类传感器要求具有防水性能；

c）将结构置于空气中（未放流体状态），获取结构动力响应信号，其或是加速度，和/或是速度，和/或是位移。

d）利用模态参数识别技术得到未放流体状态下结构的模态参数其中j为模态阶次，并将其存储入所述专用存储器中；

e）将结构置于流体中，获取其动力响应信号，其或是加速度，和/或是速度，和/或是位移。

f）利用模态参数识别技术得到流体中结构的模态参数并将其存储入所述专用存储器中；

g）流体附加质量计算：

①从专用存储器中读取上述步骤A中存储的数据，即M、K；

②结构模态参数提取

利用模态参数识别方法提取模态频率f、模态振型Φ，其中 $f_i=\sqrt{{\mathrm{\λ}}_i}\×180/\left(2\mathrm{\π}\right).$

③建立未放流体状态结构自由振动方程：其中M′、K′分别为实际结构的质量、刚度矩阵；x为结构位移，为加速度。其中未放流体状态结构与步骤A结构中M、K的关系为：

$M^{\′}=M+\underset{\mathrm{nv}=1}{\overset{\mathrm{nv}=\mathrm{Nmv}}{\mathrm{\Σ}}}{\mathrm{\η}}_{\mathrm{nv}}{M}_n+\underset{\mathrm{nw}=1}{\overset{\mathrm{nw}=\mathrm{Nmw}}{\mathrm{\Σ}}}{\mathrm{\η}}_{\mathrm{nw}}{M}_n$

$K^{\′}=K+\underset{n=1}{\overset{n=\mathrm{Nk}}{\mathrm{\Σ}}}{\mathrm{\κ}}_n{K}_n$

其中M_n、K_n为有限元模型中第n个单元质量、刚度矩阵在整体坐标系下的形式；κ_n为单元总数为Nk时第n个单元的修正系数；η_nv为位于流体中的单元总数为Nmv时第nv个单元的修正系数；同理，η_nw为位于流体中的单元总数为Nmw时第nw个单元的修正系数。

本步骤在计算流体附加质量之前，可以考虑有限元模型的建模误差，并且可以将流体中单元与空气中单元分别处理，从而提高流体附加质量的计算精度。该步骤的工程意义在于本发明可以区分引起模态信息变化的因素是有限元建模误差还是流体附加质量，便于工程设计时量化有限元建模误差以及流体附加质量的大小，从而更高精度的进行结构流固耦合计算。

④建立流体中结构自由振动方程，其中M″、K″分别为实际结构并考虑流体附加质量的质量、刚度矩阵；以M_a代表流体附加质量，流体中结构与空气中结构质量矩阵、刚度矩阵关系为：

$M^{\′\′}=M^{\′}+M_a=M+\underset{\mathrm{nv}=1}{\overset{\mathrm{nv}=\mathrm{Nmv}}{\mathrm{\Σ}}}{\mathrm{\η}}_{\mathrm{nv}}{M}_n+\underset{\mathrm{nw}=1}{\overset{\mathrm{nw}=\mathrm{Nmw}}{\mathrm{\Σ}}}{\mathrm{\η}}_{\mathrm{nw}}{M}_n+\underset{\mathrm{nw}=1}{\overset{\mathrm{nw}=\mathrm{Nmw}}{\mathrm{\Σ}}}{\mathrm{\β}}_{\mathrm{nw}}{M}_n$

K″＝K′

其中β_nw为单元总数为Nmw时位于流体中的第n个单元的质量附加系数。

本步骤流体附加质量计算时，所涉及的单元只有处于流体中的部分，尤其是流体附加质量计算所采用的信息源为处于流体环境中的结构振动信息，并不过多涉及流体的特性，从而既避开了流体模拟的难题，又可以真实反映实际结构所处的流体自身及流体特点。该步骤的意义在于可以仅利用输出信息判断流体附加质量，而且输出信息仅为低阶模态，具有工程可操作性。

⑤将结构自由振动方程转换到模态领域，即

${\mathrm{\Φ}}_i^T{K}^{\′}{\mathrm{\Φ}}_j^{\′}+{\mathrm{\λ}}_j^{\′2}{\mathrm{\Φ}}_i^T{M}^{\′}{\mathrm{\Φ}}_j^{\′}=0$

${\mathrm{\Φ}}_i^T{K}^{\′\′}{\mathrm{\Φ}}_l^{\′\′}+{\mathrm{\λ}}_j^{\′\′2}{\mathrm{\Φ}}_i^T{M}^{\′\′}{\mathrm{\Φ}}_l^{\′\′}=0$

上述两式相加，得

${\mathrm{\Φ}}_i^T{K}^{\′}{\mathrm{\Φ}}_j^{\′}+{\mathrm{\λ}}_j^{\′2}{\mathrm{\Φ}}_i^T{M}^{\′}{\mathrm{\Φ}}_j^{\′}+{\mathrm{\Φ}}_i^T{K}^{\′\′}{\mathrm{\Φ}}_l^{\′\′}+{\mathrm{\λ}}_l^{\′\′2}{\mathrm{\Φ}}_i^T{M}^{\′\′}{\mathrm{\Φ}}_l^{\′\′}=0$

其中Φ_i为有限元模型的第i阶模态振型，j为结构处于空气中时测得的模态阶次，l为结构处于流体中时测得的模态阶次。

⑥将步骤③④带入步骤⑤中第3式，得

${\mathrm{\Φ}}_i^T(K+\underset{n=1}{\overset{n=\mathrm{Nk}}{\mathrm{\Σ}}}{\mathrm{\κ}}_n{K}_n){\mathrm{\Φ}}_j^{\′}+{\mathrm{\λ}}_j^{\′2}{\mathrm{\Φ}}_i^T(m+\underset{\mathrm{nv}=1}{\overset{\mathrm{nv}=\mathrm{Nmv}}{\mathrm{\Σ}}}{\mathrm{\η}}_{\mathrm{nv}}{M}_n+\underset{\mathrm{nw}=1}{\overset{\mathrm{nw}=\mathrm{Nmw}}{\mathrm{\Σ}}}{\mathrm{\η}}_{\mathrm{nw}}{M}_n){\mathrm{\Φ}}_j^{\′}+$

${\mathrm{\Φ}}_i^T(k+\underset{n=1}{\overset{n=\mathrm{Nk}}{\mathrm{\Σ}}}{\mathrm{\κ}}_n{K}_n){\mathrm{\Φ}}_l^{\′\′}+{\mathrm{\λ}}_l^{\′\′2}{\mathrm{\Φ}}_i^T(M+\underset{\mathrm{nv}=1}{\overset{\mathrm{nv}=\mathrm{Nmv}}{\mathrm{\Σ}}}{\mathrm{\η}}_{\mathrm{nv}}{M}_n+\underset{\mathrm{nw}=1}{\overset{\mathrm{nw}=\mathrm{Nmw}}{\mathrm{\Σ}}}{\mathrm{\η}}_{\mathrm{nw}}{M}_n+\underset{\mathrm{nw}=1}{\overset{\mathrm{nw}=\mathrm{Nmw}}{\mathrm{\Σ}}}{\mathrm{\β}}_n{M}_n){\mathrm{\Φ}}_l^{\′\′}=0$

⑦步骤⑥进一步化简为

$\underset{n=1}{\overset{n=\mathrm{Nk}}{\mathrm{\Σ}}}{\mathrm{\κ}}_n(A_1+A_2)+\underset{\mathrm{nv}=1}{\overset{\mathrm{Nmv}}{\mathrm{\Σ}}}{\mathrm{\η}}_{\mathrm{nv}}(B_1+B_2)+\underset{\mathrm{nw}=1}{\overset{n=\mathrm{Nmw}}{\mathrm{\Σ}}}{\mathrm{\η}}_{\mathrm{nw}}(C_1+C_2)+\underset{\mathrm{nw}=1}{\overset{N_{\mathrm{mw}}}{\mathrm{\Σ}}}{\mathrm{\β}}_{\mathrm{nw}}D=F$

其中 $A_1={\mathrm{\Φ}}_i^T{K}_n{\mathrm{\Φ}}_j^{\′},$ $A_2={\mathrm{\Φ}}_i^T{K}_n{\mathrm{\Φ}}_l^{\′\′},$ $B_1={\mathrm{\λ}}_j^{\′2}{\mathrm{\Φ}}_i^T{M}_n{\mathrm{\Φ}}_j^{\′},$ $B_2={\mathrm{\λ}}_l^{\′2}{\mathrm{\Φ}}_i^T{M}_n{\mathrm{\Φ}}_l^{\′\′},$ $C_1={\mathrm{\λ}}_j^{\′2}{\mathrm{\Φ}}_i^T{M}_n{\mathrm{\Φ}}_j^{\′},$ $C_2={\mathrm{\λ}}_l^{\′\′2}{\mathrm{\Φ}}_i^T{M}_n{\mathrm{\Φ}}_l^{\′\′},$ $D={\mathrm{\β}}_{\mathrm{nw}}{\mathrm{\λ}}_l^{\′\′2}{\mathrm{\Φ}}_i^T{M}_n{\mathrm{\Φ}}^{\′\′},$

$F=-({\mathrm{\Φ}}_i^{T}K{\mathrm{\Φ}}_j^{\′}+{\mathrm{\λ}}_j^{\′2}{\mathrm{\Φ}}_i^T{\mathrm{M\Φ}}_j^{\′}+{\mathrm{\Φ}}_i^{T}K{\mathrm{\Φ}}_l^{\′\′}+{\mathrm{\λ}}_l^{\′\′2}{\mathrm{\Φ}}_i^{T}M{\mathrm{\Φ}}_l^{\′\′})$

⑧步骤⑦写成矩阵的形式

[A₁ A₂ B₁ B₂ C₁ C₂ D]·Δ＝F

进一步简化为

·Δ＝F

其中□＝[A₁ A₂ B₁ B₂ C₁ C₂ D]， $\mathrm{\Δ}=\left\{\begin{array}{c}\mathrm{\κ}\\ {\mathrm{\η}}_v\\ {\mathrm{\η}}_w\\ \mathrm{\β}\end{array}\right\}.$

⑨应用线性规划方法求解步骤⑧

Δ＝lsqlin(□,F,A,b,Aeq,Beq,Lb,Ub)

其中lsqlin为软件Matlab命令，并且满足A□Δ＝b，Aeq□Δ＝Beq，且Lb≤Δ≤Ub。A,b,Aeq,beq,Lb,Ub均为施加的约束条件。

⑩流体附加质量提取，并记为

${\tilde{M}}_a=\underset{\mathrm{nw}=1}{\overset{\mathrm{nw}=\mathrm{Nmw}}{\mathrm{\Σ}}}{\mathrm{\β}}_{\mathrm{nw}}{M}_n$

本发明的具体应用实例如下，以四腿导管架平台模型进行研究。

1、三维海洋平台有限元数值模型建立

如图1所示，该模型共有135个管单元组成，其中前50个单元位于流体中，其余单元位于空气中，所用材料的杨氏模量为2.1×10¹¹Pa，泊松比为0.3，密度为7860Kg/m³。

因有限元建模必然与实际结构之间存在误差，故而假定实际结构与有限元模型质量、刚度矩阵均存在一定差异，该差异通过一系列随机数考虑，即步骤③中κ_n、η_nv及η_nw，其值如图2，如3及图4所示。流体附加质量亦是借助随机数来考虑，即β_nw，具体参考图5。表1为有限元模型、实际结构在空气中，以及实际结构在流体中的前10阶频率及阻尼系数。

表1导管架式海洋平台结构有限元模型、置于空气中模型，以及置于流体中模型的前10阶频率，单位：Hz。

	有限元模型	置于空气中	置于流体中
				1	3.12	3.01	2.95
2	3.16	3.05	2.99
				3	3.67	3.56	3.49
4	5.53	5.23	4.81
				5	5.60	5.40	4.97
6	6.51	6.32	5.86
7	10.01	9.17	7.84
				8	10.19	9.27	8.19
9	11.19	10.76	9.34
				10	11.64	10.88	9.61

2、流体附加质量计算

假定置于空气中结构只有前2阶模态能够测得，置于流体中的结构也有1阶模态测得，并且有限元模型的模态全部取出。应用本发明，可以得到建模质量、刚度误差与真实值（图2-图3）的对比情况如图6、图7，附加质量各单元附加系数如图8。图5-图8表明本发明在只有低阶模态信息时可以准确计算有限元的建模误差以及流体附加质量。

总之，本实施例充分证明：1）、本发明在只有少数低阶实测模态时，既可以考虑有限元模型与实际结构的建模误差，又能准确计算流体的附加质量，并且能够以矩阵的形式表示；2）、流体附加质量计算过程中无需对结构周围流体特性进行任何假定；3）、流体附加质量可以归结到每一个结构单元上，并可以矩阵的形式存在，并且附加质量系数可以为负，符合现有物理模型试验的规律。

以上所述，仅是本发明的较佳实施例而已，并非是对本发明作其它形式的限制，任何熟悉本专业的技术人员可能利用上述揭示的技术内容加以变更或改型为等同变化的等效实施例应用于其它领域有类似要求的液体、气体的控制，但是凡是未脱离本发明技术方案内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与改型，仍属于本发明技术方案的保护范围。

Claims (1)

1.一种液体附加质量振动反演方法，其特征在于包括如下步骤：

A、建立结构有限元数值模型，获取该有限元模型的质量矩阵与刚度矩阵，分别标记为M、K，并存于专用存储器中；

B、在结构上布置加速度/速度/位移传感器；

C、将结构置于空气中获得未放液体状态下结构的模态参数λ′_j、Φ′_j，并将其存储入所述专用存储器中，其中j为模态阶次；

D、将结构置于液体中获得液体中结构的模态参数λ″_j、Φ″_j，并将其存储入所述专用存储器中；

E、液体附加质量计算：

e1、从专用存储器中读取上述步骤A中存储的数据，即M、K；

e2、结构模态参数提取

利用模态参数识别方法提取模态频率f、模态振型Φ；

e3、建立未放液体状态结构自由振动方程；所述步骤e3中，未放液体状态结构自由振动方程为其中M′、K′分别为实际结构的质量、刚度矩阵；x为结构位移，为加速度，其中未放液体状态结构与步骤A结构中M、K的关系为：

$M^{\′}=M+\underset{\mathrm{nv}=1}{\overset{\mathrm{nv}=\mathrm{Nmv}}{\mathrm{\Σ}}}{\mathrm{\η}}_{\mathrm{nv}}{M}_n+\underset{\mathrm{nw}=1}{\overset{\mathrm{nw}=\mathrm{Nmw}}{\mathrm{\Σ}}}{\mathrm{\η}}_{\mathrm{nw}}{M}_n$

$K^{\′}=K+\underset{n=1}{\overset{n=\mathrm{Nk}}{\mathrm{\Σ}}}{\mathrm{\κ}}_n{K}_n$

其中M_n、K_n为有限元模型中第n个单元质量、刚度矩阵在整体坐标系下的形式；κ_n为单元总数为Nk时第n个单元的修正系数；η_nv为位于液体中的单元总数为Nmv时第nv个单元的修正系数；η_nw为位于液体中的单元总数为Nmw时第nw个单元的修正系数；

e4、建立液体中结构自由振动方程；所述步骤e4中，液体中结构自由振动方程为其中M″、K″分别为实际结构并考虑液体附加质量的质量、刚度矩阵，以M_a代表液体附加质量，液体中结构与空气中结构质量矩阵、刚度矩阵关系为：

$M^{\′\′}=M^{\′}+M_a=M+\underset{\mathrm{nv}=1}{\overset{\mathrm{nv}=\mathrm{Nmv}}{\mathrm{\Σ}}}{\mathrm{\η}}_{\mathrm{nv}}{M}_n+\underset{\mathrm{nw}=1}{\overset{\mathrm{nw}=\mathrm{Nmw}}{\mathrm{\Σ}}}{\mathrm{\η}}_{\mathrm{nw}}{M}_n+\underset{\mathrm{nw}=1}{\overset{\mathrm{nw}=\mathrm{Nmw}}{\mathrm{\Σ}}}{\mathrm{\β}}_{\mathrm{nw}}{M}_n$

K″＝K′

其中β_nw为单元总数为Nmw时位于液体中的第n个单元的质量附加系数；

e5、利用步骤e3、e4中结构自由振动方程求得液体附加质量；用有限元模型的第i阶模态振型Φ_i分别前乘步骤e5中得到的两种状态下模态领域振动方程，

(1)、 ${\mathrm{\Φ}}_i^T{K}^{\′}{\mathrm{\Φ}}_j^{\′}+{\mathrm{\λ}}_j^{\′2}{\mathrm{\Φ}}_i^T{M}^{\′}{\mathrm{\Φ}}_j^{\′}=0$

与

(2)、 ${\mathrm{\Φ}}_i^T{K}^{\′\′}{\mathrm{\Φ}}_j^{\′\′}+{\mathrm{\λ}}_j^{\′\′2}{\mathrm{\Φ}}_i^T{M}^{\′\′}{\mathrm{\Φ}}_j^{\′\′}=0$

上述两式相加，得

(3)、 ${\mathrm{\Φ}}_i^T{K}^{\′}{\mathrm{\Φ}}_j^{\′}+{\mathrm{\λ}}_j^{\′2}{\mathrm{\Φ}}_i^T{M}^{\′}{\mathrm{\Φ}}_j^{\′}+{\mathrm{\Φ}}_i^T{K}^{\′\′}{\mathrm{\Φ}}_l^{\′\′}+{\mathrm{\λ}}_l^{\′\′2}{\mathrm{\Φ}}_i^T{M}^{\′\′}{\mathrm{\Φ}}_l^{\′\′}=0$

其中j为结构处于空气中时测得的模态阶次，l为结构处于液体中时测得的模态阶次；

将步骤e3、e4带入上述第(3)式中得

$\begin{array}{c}{\mathrm{\Φ}}_i^T(K+\underset{n=1}{\overset{n=\mathrm{Nk}}{\mathrm{\Σ}}}{\mathrm{\κ}}_n{K}_n){\mathrm{\Φ}}_j^{\′}+{\mathrm{\λ}}_j^{\′2}{{\mathrm{\Φ}}_i}^T(M+\underset{\mathrm{nv}=1}{\overset{\mathrm{nv}=\mathrm{Nmv}}{\mathrm{\Σ}}}{\mathrm{\η}}_{\mathrm{nv}}{M}_n+\underset{\mathrm{nw}=1}{\overset{\mathrm{nw}=\mathrm{Nmw}}{\mathrm{\Σ}}}{\mathrm{\η}}_{\mathrm{nw}}{M}_n){\mathrm{\Φ}}_j^{\′}+\\ {\mathrm{\Φ}}_i^T(K+\underset{n=1}{\overset{n=\mathrm{Nk}}{\mathrm{\Σ}}}{\mathrm{\κ}}_n{K}_n){\mathrm{\Φ}}_j^{\′\′}+{\mathrm{\λ}}_j^{\′\′2}{{\mathrm{\Φ}}_i}^T(M+\underset{\mathrm{nv}=1}{\overset{\mathrm{nv}=\mathrm{Nmv}}{\mathrm{\Σ}}}{\mathrm{\η}}_{\mathrm{nv}}{M}_n+\underset{\mathrm{nw}=1}{\overset{\mathrm{nw}=\mathrm{Nmw}}{\mathrm{\Σ}}}{\mathrm{\η}}_{\mathrm{nw}}{M}_n+\underset{\mathrm{nw}}{\overset{\mathrm{nw}=\mathrm{Nmw}}{\mathrm{\Σ}}}{\mathrm{\β}}_{\mathrm{nw}}{M}_n){\mathrm{\Φ}}_l^{\′\′}=0\end{array}$

进一步化简为

$\underset{n=1}{\overset{n=\mathrm{Nk}}{\mathrm{\Σ}}}{\mathrm{\κ}}_n(A_1+A_2)+\underset{\mathrm{nv}=1}{\overset{\mathrm{Nmv}}{\mathrm{\Σ}}}{\mathrm{\η}}_{\mathrm{nv}}(B_1+B_2)+\underset{\mathrm{nw}=1}{\overset{n=\mathrm{Nmw}}{\mathrm{\Σ}}}{\mathrm{\η}}_{\mathrm{nw}}(C_1+C_2)+\underset{\mathrm{nw}=1}{\overset{N_{\mathrm{mw}}}{\mathrm{\Σ}}}{\mathrm{\β}}_{\mathrm{nw}}D=F$

其中 $A_1={\mathrm{\Φ}}_i^T{K}_n{\mathrm{\Φ}}_j^{\′},A_2={\mathrm{\Φ}}_i^T{K}_n{\mathrm{\Φ}}_l^{\′\′},B_1={\mathrm{\λ}}_j^{\′2}{\mathrm{\Φ}}_i^T{M}_n{\mathrm{\Φ}}_j^{\′},B_2={\mathrm{\λ}}_l^{\′\′2}{\mathrm{\Φ}}_i^T{M}_n{\mathrm{\Φ}}_l^{\′\′},$ $C_1={\mathrm{\λ}}_j^{\′2}{\mathrm{\Φ}}_i^T{M}_n{\mathrm{\Φ}}_j^{\′},C_2={\mathrm{\λ}}_l^{\′\′2}{\mathrm{\Φ}}_i^T{M}_n{\mathrm{\Φ}}_l^{\′\′},D={\mathrm{\β}}_{\mathrm{nw}}{\mathrm{\λ}}_l^{\′\′2}{\mathrm{\Φ}}_i^T{M}_n{\mathrm{\Φ}}^{\′\′},$ $F=-({\mathrm{\Φ}}_i^T{\mathrm{K\Φ}}_j^{\′}+{\mathrm{\λ}}_j^{\′2}{\mathrm{\Φ}}_i^{T}M{\mathrm{\Φ}}_j^{\′}+{\mathrm{\Φ}}_i^{T}K{\mathrm{\Φ}}_l^{\′\′}+{\mathrm{\λ}}_l^{\′\′2}{\mathrm{\Φ}}_i^{T}M{\mathrm{\Φ}}_l^{\′\′})$

写成矩阵的形式

[A₁ A₂ B₁ B₂ C₁ C₂ D]·Δ＝F

进一步简化为

□·Δ＝F

其中□＝[A₁ A₂ B₁ B₂ C₁ C₂ D]， $\mathrm{\Δ}=\left\{\begin{array}{c}\mathrm{\κ}\\ {\mathrm{\η}}_v\\ {\mathrm{\η}}_w\\ \mathrm{\β}\end{array}\right\},$

应用线性规划方法求解

Δ＝lsqlin(□,F,A,b,Aeq,Beq,Lb,Ub)

其中lsqlin为软件Matlab命令，并且满足A□Δ≤b，Aeq□Δ＝Beq，且Lb≤Δ≤Ub，A,b,Aeq,beq,Lb,Ub均为施加的约束条件，

流体附加质量提取，并记为：

${\tilde{M}}_a=\underset{\mathrm{nw}=1}{\overset{\mathrm{nw}=\mathrm{Hmw}}{\mathrm{\Σ}}}{\mathrm{\β}}_{\mathrm{nw}}{M}_n.$