0 引　言

随着工程仿真领域大批国产化有限元（Finite Element Method，FEM）软件的设计与推广，软件核心算法和模型的有效性的验证与确认需求大大增加，数值仿真模型和算法经过完整规范的验证与确认之后再行应用推广，逐渐成为工程仿真领域业内的共识。

对舰船结构响应的仿真贯穿于舰船设计、建造、使用的全链路之中，舰船全生命周期生命力框架、舰船目标易损性分析都对舰船结构响应仿真模型的准确性提出了很高要求，因此，物理过程的数值仿真模型的实际应用必须经过一整套合理、严格、规范的验证与确认流程。

随着计算结构力学、计算流体力学建模和仿真验算法的高速发展，相关的验证与确认（Verification & Validation，V&V）理论框架和技术体系研究在国内外先后起步。美国ASME协会于2006年起成立了V&V分委员会，制定了计算流体力学V&V导则和计算结构力学V&V导则，对验证、确认、不确定度量化（Uncertainty Quantification，前三者常合称V&V&UQ）、可信度评估认证（Accreditation，与验证和确认常合称VV&A）等概念给出了定义^{[1 − 3]}。验证是确定数值计算模型是否准确反映数学模型和算法的过程，包含代码验证和计算验证。确认是确定在特定应用场景中仿真模拟对确认试验预测精确度的过程。不确定度量化对仿真结果和试验结果不确定度进行量化分析，从而给出仿真模型精确度的评价、量化建模仿真的预测能力。V&V分委员会结合受分布载荷悬臂梁响应的数值仿真案例，对计算结构力学数值仿真V&V导则、流程和方法进行了进一步的阐释。2001年，Stern等^[4]提出计算流体力学验证与确认方法和流程，并对RANS代码计算的船模阻力实施了验证与确认。国内自2003年研讨几种翼型标准模型数值模拟以来，尤其是2018年，中国空气动力研究与发展中心组织召开了第一届航空 CFD 可信度研讨会，推出了自主设计的单通道运输机标模CHN-T1；中航工业西安航空计算技术研究所开发了 CFD 验证与确认数据库，研制了 CFD 可信度分析平台 WiseCFD^[5]。

目前，常见的爆炸及毁伤力学数值仿真方法大多是对爆炸现象中物理过程的模拟，较少涉及或采用唯象模型间接涉及爆炸化学过程中物质相态状态等描述，求解的主体仍然是流体和结构的偏微分物理方程、多项式函数、复杂函数关系^[6]。同时，由于爆轰波、冲击波强间断及结构弹塑性力学行为的强非线性特性，爆炸毁伤力学对偏微分方程的具体解法和稳态的计算结构力学、计算流体力学有一定区别。

数十年间尤其是近年来，涌现出了许多针对舰船受水下爆炸冲击作用的结构毁伤的数学物理模型和仿真算法^[7-8]，例如结构有限元法^[9]、结构变分法^[10]、流体有限体积法^[11]、流体边界元法^[12]、流体有限元法^[13]、流体-结构欧拉-拉格朗日耦合^[14-15]、流体间断伽辽金法-结构有限元法耦合^[16-17]、光滑粒子法等。新算法的提出和仿真程序应用一般经过一定程度的试验验证。每一种算法都有其适用范围和工况限制，例如：流体-结构欧拉-拉格朗日耦合适用于计算舰船结构局部毁伤，受工况主要物理现象特性和仿真时效限制，其不适用于求解舰船总体毁伤；船体梁有限元适用于求解船体总纵强度毁伤但不适用于舰船局部毁伤求解。目前，舰船结构毁伤力学业内对于新算法的验证，一般验证一种算法（或仿真程序），实施一个试验，仍然缺乏系统性层层递进式的分层确认试验、标准模型试验等。综观整个爆炸力学和舰船结构毁伤力学仿真领域，与计算流体力学、计算结构力学领域相比，目前业内对仿真模型和算法的验证与确认仍然缺乏系统的认识、成熟的理论框架、专门的分析方法、标准的规范流程。

验证与确认贯穿于仿真模型和仿真程序的设计和编程整个过程中，例如飞机、汽车、高速列车、舰船等结构机械噪声的仿真程序往往从板壳受简谐激励的振动响应和辐射噪声仿真验证与确认起步，船舶螺旋桨空泡的仿真往往从敞水工况二维三维水翼的空化仿真验证与确认起步。参考国际上计算结构力学仿真验证与确认流程框架，提出舰船受水下爆炸气泡脉动压力冲击响应仿真模型验证与确认具体方法和流程。

1 仿真模型验证与分层确认方法

物理仿真建模过程主要分为概念模型建立、数学模型建立、计算模型建立等几个关键阶段^[2]。水下远场爆炸冲击舰船响应的概念模型为舰船主要产生船体总体响应，气泡脉动压力未达到船体总纵弯矩设计值时，船体主要产生弹性鞭状响应而非弹塑性响应^[18]。数学模型为船体材料力学行为主要在本构关系的线弹性段，可以采用二维船体梁有限元模型（多自由度结构系统）对船体建模。下文船体结构线性弹性响应模型主要包含模态展开法、逐步法2种关键算法（即计算模型），其中，模态展开法是求解多自由度结构系统振动响应的常用算法，逐步法特别针对非周期载荷、瞬态冲击载荷作用下结构的非线性响应。

根据图1的分层确认树状图对舰船受远场水下爆炸作用的结构响应数值仿真模型和算法代码进行验证与分层确认。单一问题层确认内容和试验如表1所示，在结构非线性冲击响应（基准过程层）中对2种算法进行跨层确认。基准过程层确认内容和确认试验如表2所示，系统层即舰船受水下远场爆炸气泡脉动压力冲击响应的确认内容为船体纵向应变、横向弯矩、鞭状响应行为、结构和附加阻尼等。

2 单一问题层的验证与确认 2.1 单自由度结构系统受一般动力载荷的响应模型

结构单自由度系统受一般动力载荷（非周期载荷、瞬态冲击载荷等）的响应有2种典型解法，叠加法和逐步法。

2.1.1 叠加法

叠加法即指杜哈梅（Duhamel）积分法。对于一单自由度有阻尼系统，其受迫运动方程为：

$ m\ddot{w}\left(t\right)+c\dot{w}\left(t\right)+kw\left(t\right)=p\left(t\right)。$

(1)

式中：w为位移；m为质量；c为阻尼；k为刚度；p为载荷；变量上标“·”为对时间t求导。

低临界阻尼系统受一般动力载荷的位移响应可由杜哈梅（Duhamel）积分求解：

$ w\left(t\right)=A\left(t\right){{\mathrm{sin}}}{\omega }_{\mathrm{d}}t-B\left(t\right){{\mathrm{cos}}}{\omega }_{\mathrm{d}}t。$

(2)

式中：$ {\omega }_{d}=\omega \sqrt{1-{\zeta }^{2}} $为有阻尼系统自振圆频率；$ \omega = \sqrt{{k}/{m}} $为无阻尼系统自振圆频率

$\begin{split}& A\left(t\right)=\frac{1}{m{\omega }_{{d}}}{\displaystyle\int }_{0}^{t}p\left(\tau \right){{\mathrm{exp}}}[-\zeta \omega (t-\tau )]{{\mathrm{cos}}}{\omega }_{{d}}\tau {\mathrm{d}}\tau，\\ & B\left(t\right)=\frac{1}{m{\omega }_{{d}}}{\displaystyle\int }_{0}^{t}p\left(\tau \right){{\mathrm{exp}}}[-\zeta \omega (t-\tau )]{{\mathrm{sin}}}{\omega }_{{d}}\tau {\mathrm{d}}\tau。\end{split} $

数值仿真中采用梯形法则数值积分求解$ A\left(t\right) $和$ B\left(t\right) $。

2.1.2 逐步法

本文采用的逐步法为Newmark-β法。非线性分析中运动方程的增量列式为：

$ \Delta {f}_{I}+\Delta {f}_{D}+\Delta {f}_{S}=\Delta p。$

(3)

式中：f_I为惯性力；f_D为阻尼力；f_S为弹簧力；p为载荷；0表示上一时刻，1表示下一时刻。

$ \left\{\begin{array}{l}\Delta {f}_{I}={f}_{I1}-{f}_{I0}=m\Delta \ddot{w}，\\ \Delta {f}_{D}={f}_{D1}-{f}_{D0}=c\left(t\right)\Delta \dot{w}，\\ \Delta {f}_{S}={f}_{S1}-{f}_{S0}=k\left(t\right)\Delta w，\\ \Delta p={p}_{1}-{p}_{0}。\end{array}\right. $

假设船体梁为线性系统，刚度和阻尼不随时间变化。时刻t的增量平衡方程为：

$ m\Delta\ddot{w}+c_0\Delta\dot{w}+k_0\Delta w=\Delta p。$

(4)

增量等效静力平衡方程为：

$ \stackrel{~}{k}\Delta w=\Delta \stackrel{~}{p}。$

(5)

采用线加速度，β=1/6，等效刚度为：

$ \stackrel{~}{k}={k}_{0}+\frac{3{c}_{0}}{\Delta t}+\frac{6m}{\Delta {t}^{2}}。$

(6)

等效载荷增量为：

$ \Delta \stackrel{~}{p}=\Delta p+{c}_{0}\left(3{\dot{w}}_{0}+\frac{\Delta t}{2}{\ddot{w}}_{0}\right)+m\left(\frac{6}{\Delta t}{\dot{w}}_{0}+3{\ddot{w}}_{0}\right)。$

(7)

速度增量为：

$ \Delta \dot{w}=\frac{3}{\Delta t}\Delta w-3{\dot{w}}_{0}-\frac{\Delta t}{2}{\ddot{w}}_{0}。$

(8)

Newmark-β法的每一时间步内计算步骤为：

步骤1　利用初始载荷、初始速度、初始位移、阻尼和刚度求出初始加速度。

$ \ddot{w}_0=\frac{1}{m}\left(p_0-f_{D0}-f_{S0}\right)。$

步骤2　计算等效刚度和等效载荷增量。

步骤3　计算位移增量和速度增量。

步骤4　求出下一时间步的速度和位移。

$ \left\{\begin{array}{c}{w}_{1}={w}_{0}+\Delta w，\\ {\dot{w}}_{1}={\dot{w}}_{0}+\Delta \dot{w}。\end{array}\right. $

2.2 多自由度结构系统受一般动力载荷的响应模型

多自由度结构系统受一般动力载荷激励的响应模型包含模态展开法和逐步法。

2.2.1 模态展开法

线性n自由度系统有阻尼强迫振动的运动方程为：

$ \boldsymbol{M}\ddot{\boldsymbol{w}}+\boldsymbol{C}\dot{\boldsymbol{w}}+\boldsymbol{K}\boldsymbol{w}=\boldsymbol{F}\left(t\right)。$

(9)

式中：位移为n行列向量，$ w={\left({w}_{1},{w}_{2},\cdots {w}_{n}\right)}^{\rm T} $，质量、刚度、阻尼为n行n列矩阵，载荷为n行列向量。与之对应的无阻尼自由振动运动方程为：

$ \boldsymbol{M}\ddot{\boldsymbol{w}}+\boldsymbol{K}\boldsymbol{w}=0。$

(10)

特征方程式为$ \left|\boldsymbol{K}-{\omega }^{2}\boldsymbol{M}\right|=0 $，ω为系统固有频率。

位移的振型展开解为$ \boldsymbol{w}\ =\ \boldsymbol{\varPhi }\boldsymbol{q} $，振型矩阵$ \boldsymbol{\varPhi }= \left[{\boldsymbol{\varphi }}_{1}{\boldsymbol{\varphi }}_{2}{\boldsymbol{\varphi }}_{3}\cdots {\boldsymbol{\varphi }}_{n}\right] $，振型坐标$ \boldsymbol{q}={\left({q}_{1},{q}_{2},\cdots {q}_{n}\right)}^{\rm T} $

将位移振型展开解代入式(10)，位移或振型坐标都是时间t的正弦函数，因此式(10)可化为：

$ ( - {{\boldsymbol{\omega}} ^2}{\boldsymbol{M}} + {\boldsymbol{K}}){\boldsymbol{\varPhi}} = 0。$

(11)

即特征值（固有频率的平方）与特征矢量（振型）一一对应。

对于有阻尼系统，将位移振型展开解代入运动方程，并每一项左乘Φ^T，得

$ \stackrel{~}{\boldsymbol{M}}\ddot{\boldsymbol{q}}+\stackrel{~}{\boldsymbol{C}}\dot{\boldsymbol{q}}+\stackrel{~}{\boldsymbol{K}}\boldsymbol{q}=\stackrel{~}{\boldsymbol{F}}。$

(12)

式中：${\boldsymbol{{M}}}={\boldsymbol{\varPhi }}^{\bf{T}}\boldsymbol{M}\boldsymbol{\varPhi } $；$ \stackrel{~}{\boldsymbol{C}}={\boldsymbol{\varPhi }}^{\bf{T}}\boldsymbol{C}\boldsymbol{\varPhi } $；$ \stackrel{~}{\boldsymbol{K}}={\boldsymbol{\varPhi }}^{\bf{T}}\boldsymbol{K}\boldsymbol{\varPhi } $；$\stackrel{~}{\boldsymbol{F}}={\boldsymbol{\varPhi }}^{\bf{T}} \boldsymbol{F}\left(t\right) $。

振型关于质量和刚度矩阵有如下正交条件：

$ {{\boldsymbol{\varphi }}_{i}}^{\rm T}\boldsymbol{M}{\boldsymbol{\varphi }}_{j}=\left\{\begin{array}{l}0,i\ne j，\\ {\stackrel{~}{m}}_{j},i=j。\end{array}\right. $

(13)

$ {{\boldsymbol{\varphi }}_{i}}^{\rm{T}}{\boldsymbol{K}}{\boldsymbol{\varphi }}_{j}=\left\{\begin{array}{c}0,i\ne j，\\ {\stackrel{~}{k}}_{j},i=j。\end{array}\right. $

(14)

假设系统的振型关于阻尼矩阵有正交关系：

$ {{\boldsymbol{\varphi }}_{i}}^{\rm T}\boldsymbol{C}{\boldsymbol{\varphi }}_{j}=\left\{\begin{array}{l}0,i\ne j，\\ {\stackrel{~}{c}}_{j}=2{\omega }_{j}{\zeta }_{j}{\stackrel{~}{m}}_{j},i=j。\end{array}\right. $

(15)

式中：ζ_j为第j阶振型阻尼比。

从而n阶耦合的运动方程可以解耦为n个相互独立的单自由度运动方程：

$ {\ddot{\boldsymbol{q}}}_{j}+2{\zeta }_{j}{\omega }_{j}{\dot{\boldsymbol{q}}}_{j}+{\omega }_{j}^{2}{\boldsymbol{q}}_{j}=\frac{{\stackrel{~}{f}}_{j}}{{\stackrel{~}{m}}_{j}}。$

(16)

当初始条件$ \boldsymbol{q}\left(0\right)=\dot{\boldsymbol{q}}\left(0\right)=0 $，式(16)解即振型坐标为：

$ {q}_{j}\left(t\right)=\frac{1}{{\stackrel{~}{m}}_{j}{\omega }_{{d}_{j}}}{\int }_{0}^{t}f\left(\tau \right){{\mathrm{exp}}}[-{\zeta }_{j}{\omega }_{j}(t-\tau )]\cdot {{\mathrm{sin}}}{\omega }_{{d}_{j}}(t-\tau ){\mathrm{d}}\tau。$

(17)

式中t ≥ 0，任意阶振型有阻尼自振频率$ {\omega }_{{d}_{j}}={\omega }_{j}\sqrt{1-{\zeta }_{j}^{2}} $。

2.2.2 逐步法

多自由度系统的Newmark-β法定义与公式与单自由度系统的基本一致，只是其位移、速度、加速度、载荷是矢量形式，并且，结构系统刚度、阻尼、质量为矩阵形式。

2.3 单自由度钢架受冲击响应逐步法模型的验证与确认

单自由度钢架如图2所示，钢架质量m=0.1 kips.s²/m、刚度k=5 kips/m、阻尼c=0.2 kips.s/m，钢架受到横向冲击载荷，外载荷时程见图3(a)。根据上述结构单自由度系统的逐步法编写仿真程序代码，求解钢架响应并与数值基准解（Numerical Benchmark Solution）进行对比。由图3可见，模型仿真程序代码计算结果与数值基准解吻合，该数值仿真程序代码正确性、数值计算精确度得以验证。

2.4 多自由度钢架受冲击响应模态展开逐步法模型的验证与确认

一个3层钢架如图4所示，各层所受横向正弦爆炸冲击载荷为p₁(t)=500cos(π/0.02t)，p₂(t)=1 000cos(π/0.02t)，p₃(t)=1000cos(π/0.02t)，单位kips。采用结构三自由度系统对钢架建模，质量M=[1 0 0;0 1.5 0;0 0 2.0] kips·s²/m、刚度K=600×[1 −1 0;−1 3 −2;0 −2 5] kips/m，系统无阻尼。采用多自由度结构系统的模态展开逐步法求解系统受横向正弦冲击的响应，并且与模态展开叠加法计算结果（属于数值基准解）进行对比，如图5所示。可见各节点逐步法计算值与叠加法计算值吻合较好，幅值相对误差平均值−4.39%，多自由度结构系统的模态展开逐步法程序代码正确性、数值计算精确性得以确认。

3 基准过程层的确认 3.1 船体梁自由弯曲振动有限元模型

船体梁等大多数梁是非均匀梁，难以求得振型的精确表达式，这里采用有限元法引入瑞利-里兹法的基本思想，以近似方法求响应。考虑一种梁单元体，节点广义坐标分别为左端节点横向位移、左端节点转角、右端节点横向位移、右端节点转角。任意梁单元体内各点的位移采用三次多项式展开，$ w(x,t)={\alpha }_{1}\left(t\right)+ {\alpha }_{2}\left(t\right)x+{\alpha }_{3}\left(t\right){x}^{2}+{\alpha }_{4}\left(t\right){x}^{3} $。基函数$ \boldsymbol{H}\left(x\right)=\left[1\;\; x\;\; {x}^{2}\;\; {x}^{3}\right] $，广义坐标$ \boldsymbol{\alpha }\left(t\right)=[{\alpha }_{1}(t)\;\;\; {\alpha }_{2}(t)\;\;\; {\alpha }_{3}(t)\;\;\; {\alpha }_{4}(t)]^{\rm T} $。对于欧拉-伯努利梁（纯弯曲梁），梁单元节点位移坐标为$ {\boldsymbol{w}}_{N}\left(t\right)= [w(0,t)\;\;\;\; w'(0,t)\;\;\;\; w(l,t)\;\;\;\; w'(l,t)]^{\rm T} $。节点坐标和广义坐标的关系为$ {\boldsymbol{w}}_{N}={\boldsymbol{B}}^{-1}\boldsymbol{\alpha } $，其中，

$ {\boldsymbol{B}}^{-1}=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 1& l& {l}^{2}& {l}^{3}\\ 0& 1& 2l& 3{l}^{2}\end{array}\right]。$

梁单元体内各点位移可采用节点位移坐标表示：

$ w(x,t)=\boldsymbol{H}\left(x\right)\boldsymbol{\alpha }\left(t\right)=\boldsymbol{H}\left(x\right)\boldsymbol{B}{\boldsymbol{w}}_{N}=\boldsymbol{{\mathrm{\mathit\Psi}} }{\boldsymbol{w}}_{N}。$

(18)

式中：$ {{\mathrm{\mathit{\boldsymbol \Psi}}} }\left(x\right) $为梁单元体形状函数。

梁单元为均匀直梁，其自由运动方程为：

$ EIw^{\text{IV}}+m\ddot{w}=0。$

(19)

梁单元刚度矩阵为：

$ \begin{split} {K}_{{e}} = {\int }_{0}^{l}{\boldsymbol{{\mathrm{\mathit\Psi}} }}''^{\rm T}EI{\boldsymbol{{\mathrm{\mathit\Psi}} }}^{''} {\rm d}x = \text{}\frac{EI}{{l}^{3}} \left[\begin{array}{cccc} 12& 6l& -12& 6l \\ 6l& 4{l}^{2}& -6l& 2{l}^{2} \\ -12& -6l& 12& -6l \\ 6l& 2{l}^{2}& -6l& 4{l}^{2} \end{array}\right]。\end{split}$

(20)

梁单元质量矩阵为：

$ \begin{split} {\boldsymbol{M}}_{{e}} = {\int }_{0}^{l}{\boldsymbol{{\mathrm{\mathit\Psi}} }}^{\rm T}m\boldsymbol{{\mathrm{\mathit\Psi}} }{\rm d}x = \text{}\frac{ml}{420} \left[\begin{array}{cccc} 156& 22l& 54& -13l \\ 22l& 4{l}^{2}& 13l& -3{l}^{2} \\ 54& 13l& 156& -22l \\ -13l& -3{l}^{2}& -22l& 4{l}^{2} \end{array}\right]。\end{split}$

(21)

设船体梁包含n个梁单元，组装系统刚度矩阵、系统质量矩阵（一致质量法）和系统外力矢量：

$ \left\{\begin{array}{c}K={\displaystyle\sum }_{r=1}^{n}{\boldsymbol{A}}^{\left(r\right)T}{\boldsymbol{K}}_{\mathrm{e}}^{\left(r\right)}{\boldsymbol{A}}^{\left(r\right)}，\\ M={\displaystyle\sum }_{r=1}^{n}{\boldsymbol{A}}^{\left(r\right)T}{\boldsymbol{M}}_{\mathrm{e}}^{\left(r\right)}{\boldsymbol{A}}^{\left(r\right)}。\end{array}\right. $

(22)

式中：$ {\boldsymbol{A}}^{\left(r\right)} $为梁系统节点坐标U至第r个梁单元节点坐标w_N^(r)的转化矩阵，$ {{\boldsymbol{w}}_{N}}^{\left(r\right)}={\boldsymbol{A}}^{\left(r\right)}\boldsymbol{U} $。

$ {\boldsymbol{A}}^{\left(r\right)} $矩阵维度为4行2(n+1)列，除下列元素为1以外，转化矩阵的其余元素为0。

$ \mathit{A}^{\left(r\right)}\left(i,2\left(r-1\right)+i\right)=1,(i=\mathrm{1,2},\mathrm{3,4})。$

(23)

船体梁有限元无阻尼自由弯曲振动方程为：

$ \mathit{M}\ddot{\mathit{U}}+\mathit{K}\mathit{U}=0。$

(24)

船体梁除计入自身质量外，还需要计入附连水质量。船体梁垂向振动的单位长度附连水质量可以根据下式和经验图谱求解。

$ {m}_{\text{ad}}=\frac{1}{2}{\alpha }_{V}{K}_{i}{C}_{V}\rho {\text π} {b}^{2}。$

(25)

式中：α_V为浅水修正系数；K_i为三维流动修正系数；C_V为垂向振动附连水质量系数；b为剖面的水线半宽。

船体材料弹性模量、材料密度或梁单元质量、船体梁各站横截面惯性矩、附连水质量等共同构成了船体梁的总体刚度和总体质量，根据船体梁无阻尼横向自由运动方程，得出其对应的特征方程，求解该特征方程即得到船体梁横向自由振动固有频率和振型。由于细长船体的鞭状响应以弹性响应为主，故此处滤除刚体运动模态及其对应零频率。

3.2 水下爆炸气泡脉动压力模型

船体梁所受水下爆炸气泡脉动压力载荷由Geers-Hunter模型^[19]数值计算得出，数值计算方法为龙格-库塔法。Geers-Hunter模型可以解得气泡半径、气泡迁移量、气泡迁移速度、气泡径向速度时间历程，气泡半径为a，球状气泡体积为V。船体鞭状响应主要由气泡脉动压力诱导产生，爆炸初始的冲击波压力不计入船体外载荷。气泡脉动阶段脉动气泡对远场空间任一点的压力为：

$ P=\frac{\rho }{4{\text π} R}\ddot{V}。$

(26)

考虑入射角（船体梁节点与爆心连线与船体纵轴夹角）为θ，并考虑壁面反射，则任一船体梁微元（半宽b，长度dl）所受气泡脉动压力垂向载荷为

$ p=\kappa \frac{\rho }{4{\text π} R}\ddot{V}\cdot {{\mathrm{sin}}}\theta \cdot 2b\cdot {\mathrm{d}}l。$

(27)

式中：κ为壁面反射系数，取1.8。外载荷单位是1 kg·m/s²即1 N，量纲为L·m·t⁻²。

3.3 船体梁自由弯曲振动有限元模型的确认

对比船体梁横向自由振动固有振型和频率的仿真计算值与试验值并统计误差，以验证与确认船体梁有限元模型的正确性。

案例1船模^[20]，长2800 mm，宽度300 mm，船高80 mm，舱壁数量6，舱壁间距400 mm，板厚3 mm。船舶附连水质量采用一阶附连水质量。采用二维有限单元梁建模，船模划分成12个梁单元，横向弯曲一阶固有频率计算值20.9 Hz，两端自由梁横向弯曲固有频率理论公式值20.5 Hz，试验值23.3 Hz，代码计算值与试验值相对误差10.3%。前3阶固有模态如图6所示。

案例2船模为中国船舶科学研究中心推出的受水下远场爆炸船体鞭状响应的标准模型，长7700 mm，型宽890 mm，型深400 mm，吃水280 mm，重量约900 kg，共21站，各站半宽、吃水、质量、截面垂向惯性矩、截面垂向剪切面积、水线面面积、浸水面积参见文献[21]。采用二维有限单元梁建模，船模划分成21个梁单元，横向弯曲湿模态一阶固有频率计算值10.4 Hz，试验值10.6 Hz，代码计算值与试验值相对误差1.9%。前3阶固有模态如图7所示。

2个案例中，船体梁有限元模型固有频率计算值与试验值误差较小，验证与确认了船体梁有限元模型也即船体梁刚度、质量、附加质量代码的正确性和计算的准确性。

3.4 水下爆炸气泡运动和脉动压力模型的确认

对于水下远场（远离自由表面、固壁，有重力）爆炸气泡，确认内容包括第一次脉动周期、气泡半径、气泡第一次脉动压力波形和峰值。

案例1工况为5 g球形TNT装药水下1.04 m深度爆炸，爆距1 m。气泡脉动阶段气泡半径如图8、压力时程曲线如图9所示，Geers-Hunter模型龙格-库塔法程序代码计算结果张阿漫气泡统一理论中水下爆炸气泡数值仿真程序计算值^[22]进行了对比，代码计算值与Cole经验公式值误差见表1，三者的气泡最大半径R_max、气泡第一次脉动周期T吻合程度较高。气泡脉动压力峰值计算值与Zamyshlyayev公式值、张阿漫仿真软件计算值较为吻合。

案例2工况为90 g球形TNT装药水下1.5355 m深度爆炸，爆距2.44 m。气泡脉动阶段气泡半径见图10、压力时程曲线见图11。Geers-Hunter模型龙格-库塔法程序计算结果和张阿漫气泡水下爆炸气泡数值仿真程序计算值进行了对比，代码计算值与Cole经验公式值误差见表3，三者的气泡最大半径、气泡第一次脉动周期吻合程度较高。气泡脉动压力峰值计算值与Zamyshlyayev公式值、张阿漫仿真软件计算值较为吻合。

4 系统层的确认 4.1 船体梁受迫弯曲运动有限元模型

假设船体梁中任意梁单元只有左端节点受到集中外力：

$ \boldsymbol{F}_e=\left[\begin{array}{cccc}f_1\left(t\right) & 0 & 0 & 0\end{array}\right]^{\rm{T}}。$

(28)

n个梁单元构成的船体梁系统外力矢量为$ \boldsymbol{F}= {\sum }_{r=1}^{n}{\boldsymbol{A}}^{\left(r\right){\rm T}}{\boldsymbol{F}}_{\boldsymbol{e}}^{\left(r\right)} $。

船体梁有限元受迫弯曲运动方程为：

$ \boldsymbol{M}\ddot{\boldsymbol{U}}+\boldsymbol{C}\dot{\boldsymbol{U}}+\boldsymbol{K}\boldsymbol{U}=\boldsymbol{F}。$

(29)

n个梁单元的船体梁有限元系统具有2(n+1)个自由度，采用模态展开法将该多自由度系统解耦成m个相互独立的单自由度系统，再采用模态展开法和逐步法求解船体运动响应。

采用瑞利阻尼，为计入流体的影响，振型阻尼比采用文献[23]经验公式计算：

$ \zeta=1\ 085\left(\frac{KSF}{f_1\cdot L^{1.5}}\right)^{1.31}。$

(30)

式中：${KSF}={\left[1+\sin \left(\alpha \right)\right]}/{2}\cdot {\sqrt{W}}/{R} $为龙骨冲击因子；α为装药迎角；W为装药质量；R为爆距；f₁为船体振动一阶频率；L为船长。船体结构阻尼还可根据复变指数法^[24]求解。

船体梁纵向（x向）和横向（y向）应变的求解。对于欧拉梁假设，船体梁各处弯矩和纵向应变为，

$ M=EI\frac{{\rm d}^{2}w}{{\rm{d}}{x}^{2}} \text{，} {\text{ε}}_{x}=\frac{M}{EI}y。$

(31)

式中：y为梁横截面一点至中性轴的距离。

4.2 水下远场爆炸气泡脉动压力冲击舰船响应模型的确认

案例1炸药位于船舯正下方，船模中部舱底处x向应变（微应变με）的仿真计算值与试验值见图12。采用评估舰船抗远场水下爆炸冲击数值仿真准确性的Russell误差因子^[25]（包含幅值误差因子R_M、相位误差因子R_P、综合误差因子R_C），仿真值与试验值（统计时间为0～0.1 s）的R_M=0.13，R_P=0.23，R_C=0.23，综合误差因子R_C在(0.15,0.28]范围内，仿真结果可接受。船模鞭状响应幅值的仿真与试验一致性、周期的仿真与试验一致性均较好。

案例2爆心对船模攻角39°（药心与船中观测点连线与水平面的夹角）。采用模态展开逐步法程序代码对船模鞭状响应进行了仿真，船模中部弯矩仿真计算值、试验值对比见图13。

船体梁中部弯矩计算值与试验值（统计时间为0～0.2 s）的R_M=0.001，R_P=0.12，R_C=0.11。综合误差因子R_C在(0,0.15]范围内，仿真结果较好。船模鞭状响应幅值的仿真与试验一致性、周期的仿真与试验一致性均较好。

5 结　语

本文借鉴计算结构力学、计算流体力学数值仿真的验证与确认关键技术，以水下远场爆炸气泡脉动冲击舰船响应数值仿真程序的验证与确认为例，提出舰船抗水下远场爆炸冲击结构毁伤数值仿真的验证与确认方法、策略与流程。

模型验证包括代码验证和数值计算验证，分别检验仿真程序的错误纰漏、数值计算的精确性。模型确认主要对仿真结果和试验结果进行一致性分析，从而评估仿真模型在特定工况域的可信度。

水下远场爆炸气泡脉动冲击舰船响应仿真模型的确认策略为：单一问题层确认结构响应的模态展开法和逐步法；基准过程层确认结构非线性冲击响应、水下爆炸气泡脉动压力和船体梁自由弯曲振动；系统层确认水下远场爆炸气泡脉动压力作用下船体鞭状响应。

水下近场爆炸作用下流体与结构耦合作用较强，确认策略为：单一问题层确认材料属性和弹塑性本构或状态方程、爆轰产物状态方程、水的状态方程、自由场水下爆炸气泡脉动过程；基准过程层确认自由场水下爆炸冲击波传播、流场空化现象、近壁面近自由面气泡脉动过程；子系统确认水下爆炸对板壳结构毁伤过程；系统层确认水下爆炸对典型板架、舱段、缩比船模毁伤过程。