水下爆炸作用下舰船整体弹塑性运动响应预测方法及系统

摘要

本发明公开了一种水下爆炸作用下舰船整体弹塑性运动响应预测方法及系统，方法包括步骤：S1、将舰船简化为等截面自由船体梁；S2、将水下爆炸作用下自由船体梁载荷压力划分为五个时间阶段；S3、求解自由船体梁在水下爆炸五个时间阶段内的运动位移；S4、求解船体梁的弯矩；S5、求解船体梁正向塑性运动参数；S6、求解船体梁反向弹性卸载运动参数；S7、求解船体梁反向塑性运动参数；S8、求解船体梁正向弹性运动参数。本发明还提供了实现上述方法的系统。本发明在综合考虑冲击波和气泡载荷对舰船总体的冲击作用的基础上分析舰船塑性运动过程，可以准确且简便地实现对水下爆炸作用下舰船整体弹塑性运动变形的工程预报。

说明书

技术领域

本发明属于舰船总体结构抗水下爆炸冲击变形技术领域，更具体地， 涉及一种水下爆炸作用下舰船整体弹塑性运动响应预测方法及系统。

背景技术

现代海战中，舰船遭受的水下爆炸攻击多以近距作用为主，其中水下 近距非接触爆炸形成的冲击作用遍及整船，往往造成大范围设备破坏甚至 造成舰船整体折断。2010年3月26日发生的韩国“天安号”护卫舰遭水下 爆炸攻击发生瞬间折断沉没事件，就可看作是水下近距非接触爆炸造成舰 船折断破坏的典型例子。

对于水下近距非接触爆炸作用下的舰船而言，其最先受到冲击波的作 用，随后爆炸气泡脉动及辐射压力对结构产生后续影响。考虑到冲击波和 气泡载荷发生的先后关系，已有研究往往将两者分开研究或者直接忽视冲 击波对结构整体变形的影响，但从物理过程而言，冲击波造成的结构响应 应当作为计算后续气泡造成的结构响应的初始输入条件，其影响不应忽视。 另外，水下近距爆炸作用会造成舰船发生明显的整体变形，根据爆炸工况 的不同，可能发生中垂或中拱弯曲破坏，这点已被实船爆炸或模型试验研 究所证实。目前，国内尚缺乏一套完整考虑爆炸冲击波和气泡联合作用， 计算舰船发生反复中拱或中垂塑性变形的理论计算方法。另外，基于所建 立的舰船整体响应计算方法，进一步研究舰船整体响应特点及损伤模式， 对于提高舰船抗水下爆炸冲击防护设计水平，特别是对指导水中兵器攻击 效能优化设计及攻击方式选取等均具有重要意义。

发明内容

针对现有技术的缺陷和技术需求，本发明提供了一种水下爆炸作用下 舰船整体弹塑性运动响应预测方法及系统，本发明综合考虑冲击波和气泡 载荷对舰船总体的冲击作用，能合理反映舰船反复加载、卸载的响应过 程，可以较为简便且准确地实现对近距爆炸作用下舰船弹塑性运动变形的 工程预报。

为实现上述目的，按照本发明，提供了一种水下爆炸作用下舰船整体 弹塑性运动响应预测方法，所述方法包括步骤：

S1、将舰船等效为等截面自由船体梁，由此将对舰船运动响应的预测 等效为对自由船体梁运动响应的预测；

S2、将水下爆炸作用下自由船体梁载荷压力过程划分为五个时间阶段， 采集冲击波压力峰值P_m、气泡脉动阶段负压峰值P_b、气泡第一次脉动压力 峰值P_s；所述五个时间阶段为：0≤t＜t₁、t₁≤t＜t₂、t₂≤t＜t₃、t₃≤t＜t₄、t₄≤t＜t₅； 其中t₁＝θ，θ为冲击波衰减常数；P₀为炸药处静水压力，P_atm为大气压，c为水中 声速，r₀为装药半径，R为爆距；m_e为装药当量，ρ_w为水的密度，g为重力加速度； $t_5=3290\frac{r_0}{P_0^{0.71}}+t_4;$

S3、求解自由船体梁在水下爆炸五个时间阶段内的运动位移w(x,t)：

其中，x为自由船体梁任意点处的横坐标值，k₂＝t₂-t₁， 为自由船体梁第一阶固有振形， 其中ζ₁为振形幅值，μ₁为梁运动 频率参数，l为自由船体梁长度；为五个时间阶段的积分常数，ω₁为自由梁的一阶振形固有频率；φ₁为积分常 数，m为考虑附连水质量的单位梁长度质量，p(x)为自由 船体梁不同时间阶段压力分布特征函数，对于0≤t＜t₁、t₁≤t＜t₂时间阶段， 对于第t₂≤t＜t₃、t₃≤t＜t₄、t₄≤t＜t₅时间阶段， $p\left(x\right)=(1-\frac{2x-l}{l})·\exp [-8{\left(\frac{2x-l}{l}\right)}^2+4{\left(\frac{2x-1}{l}\right)}^3]+{\left(\frac{2l-2x}{l}\right)}^{1.5}·\left(\frac{2x-l}{l}\right);$

S4、根据船体梁在五个时间阶段的弹性运动位移w(x,t)，求解船体梁在 不同阶段的速度和弯矩值；

S5、当所述船体梁中点弯矩绝对值超过船体梁塑性极限弯矩M_s绝对值 时，船体梁进入正向塑性运动，求解船体梁正向塑性运动时的运动幅值H'(t) 和相对转动的转角α(t)：

$\left(\begin{array}{c}{\stackrel{··}{H}}^{\prime }\left(t\right)=\frac{3M_s-(l{\xi }_3+3{\xi }_4)p\left(t\right)}{(l{\xi }_1+3{\xi }_2)m}\\ \stackrel{··}{\alpha }\left(t\right)=\frac{24({\xi }_1{\xi }_4-{\xi }_2{\xi }_3)p\left(t\right)-24M_s{\xi }_1}{(l{\xi }_1+3{\xi }_2){\mathrm{ml}}^2}\end{array}\right)$

ξ₁＝λ₁(l)-λ₁(l/2)，ξ₂＝λ₂(l)-λ₂(l/2)-λ₁(l)·l/2，

ξ₃＝η₁(l/2)-η₁(l)，ξ₄＝η₂(l/2)-η₂(l)+η₁(l)·l/2；

S6、当所述船体梁正向塑性变形达到最大即α(t)达到最大值时，其正向 塑性变形阶段结束，进入反向弹性卸载，根据公式 $\left(\begin{array}{c}\alpha \left(t\right)(x-l/2)+{\phi }_1\left(x\right)H^{\prime }\left(t\right)={\phi }_1\left(x\right)H\left(t\right)\\ \stackrel{·}{\alpha }\left(t\right)(x-l/2)+{\phi }_1\left(x\right){\stackrel{·}{H}}^{\prime }\left(t\right)={\phi }_1\left(x\right)\stackrel{·}{H}\left(t\right)\end{array}\right)$求解船体梁反向弹性卸载过程的运动幅值 H(t)；

S7、在船体梁反向弹性变形过程中，当梁中点弯矩绝对值超过塑性极 限弯矩M_s绝对值时，船体梁进入反向塑性运动，求解梁进入反向塑性变形 后的运动幅值H₁'(t)和相对转动的转角α₁(t)：

$\left(\begin{array}{c}{\stackrel{··}{H}}_1^{\prime }\left(t\right)=\frac{3(-M_s)-(l{\xi }_3+3{\xi }_4)p\left(t\right)}{(l{\xi }_1+3{\xi }_2)m}\\ {\stackrel{··}{\alpha }}_1\left(t\right)=\frac{24({\xi }_1{\xi }_4-{\xi }_2{\xi }_3)p\left(t\right)-24\left({-M}_s\right){\xi }_1}{(l{\xi }_1+3{\xi }_2){\mathrm{ml}}^2}\end{array}\right);$

S8、当所述船体梁反向塑性变形达到最大即α₁(t)达到最大值时，其反 向塑性变形阶段结束，进入正向弹性卸载，根据 $\left(\begin{array}{c}{\alpha }_1\left(t\right)(x-l/2)+{\phi }_1\left(x\right)H_1^{\prime }\left(t\right)={\phi }_1\left(x\right)H_1\left(t\right)\\ {\stackrel{·}{\alpha }}_1\left(t\right)(x-l/2)+{\phi }_1\left(x\right){\stackrel{·}{H}}_1^{\prime }\left(t\right)={\phi }_1\left(x\right){\stackrel{·}{H}}_1\left(t\right)\end{array}\right)$求解船体梁正向弹性卸载阶段的运动幅值 H₁(t)。

相应地，本发明还提供一种水下爆炸作用下舰船整体弹塑性运动响应 预测系统，所述系统包括：

第一模块，用于将舰船等效为等截面自由船体梁，由此将对舰船运动 响应的预测等效为对自由船体梁运动响应的预测；

第二模块，用于将水下爆炸作用下自由船体梁载荷压力过程划分为五 个时间阶段，采集冲击波压力峰值P_m、气泡脉动阶段负压峰值P_b、气泡第 一次脉动压力峰值P_s；所述五个时间阶段为：0≤t＜t₁、t₁≤t＜t₂、t₂≤t＜t₃、 t₃≤t＜t₄、t₄≤t＜t₅；其中t₁＝θ，θ为冲击波衰减常数；P₀为炸药处静水压力，P_atm为大气压， c为水中声速，r₀为装药半径，R为爆距；m_e为装药当量，ρ_w为水的密度，g为重力加速度； $t_5=3290\frac{r_0}{P_0^{0.71}}+t_4;$

第三模块，用于求解自由船体梁在水下爆炸五个时间阶段内的运动位 移w(x,t)：

其中，x为自由船体梁任意点处的横坐标值，k₂＝t₂-t₁， 为自由船体梁第一阶固有振形， 其中ζ₁为振形幅值，μ₁为梁运动 频率参数，l为自由船体梁长度；为五个时间阶段的积分常数，ω₁为自由梁的一阶振形固有频率；φ₁为积分 常数，m为考虑附连水质量的单位梁长度质量，p(x)为自 由船体梁不同时间阶段压力分布特征函数，对于0≤t＜t₁、t₁≤t＜t₂时间阶段， 对于第t₂≤t＜t₃、t₃≤t＜t₄、t₄≤t＜t₅时间阶段， $p\left(x\right)=(1-\frac{2x-l}{l})·\exp [-8{\left(\frac{2x-l}{l}\right)}^2+4{\left(\frac{2x-1}{l}\right)}^3]+{\left(\frac{2l-2x}{l}\right)}^{1.5}·\left(\frac{2x-l}{l}\right);$

第四模块，用于根据船体梁在五个时间阶段的弹性运动位移w(x,t)，求 解船体梁在不同阶段的速度和弯矩值；

第五模块，用于当所述船体梁中点弯矩绝对值超过船体梁塑性极限弯 矩M_s绝对值时，判断船体梁进入正向塑性运动，求解船体梁正向塑性运动 时的运动幅值H'(t)和相对转动的转角α(t)：

$\left(\begin{array}{c}{\stackrel{··}{H}}^{\prime }\left(t\right)=\frac{3M_s-(l{\xi }_3+3{\xi }_4)p\left(t\right)}{(l{\xi }_1+3{\xi }_2)m}\\ \stackrel{··}{\alpha }\left(t\right)=\frac{24({\xi }_1{\xi }_4-{\xi }_2{\xi }_3)p\left(t\right)-24M_s{\xi }_1}{(l{\xi }_1+3{\xi }_2){\mathrm{ml}}^2}\end{array}\right),$

ξ₁＝λ₁(l)-λ₁(l/2)，ξ₂＝λ₂(l)-λ₂(l/2)-λ₁(l)·l/2，

ξ₃＝η₁(l/2)-η₁(l)，ξ₄＝η₂(l/2)-η₂(l)+η₁(l)·l/2；

第六模块，用于当所述船体梁正向塑性变形达到最大即α(t)达到最大值 时，判断船体梁正向塑性变形阶段结束，进入反向弹性卸载，根据公式 $\left(\begin{array}{c}\alpha \left(t\right)(x-l/2)+{\phi }_1\left(x\right)H^{\prime }\left(t\right)={\phi }_1\left(x\right)H\left(t\right)\\ \stackrel{·}{\alpha }\left(t\right)(x-l/2)+{\phi }_1\left(x\right){\stackrel{·}{H}}^{\prime }\left(t\right)={\phi }_1\left(x\right)\stackrel{·}{H}\left(t\right)\end{array}\right)$求解船体梁反向弹性卸载过程的运动幅值 H(t)；

第七模块，用于在船体梁反向弹性变形过程中，当梁中点弯矩绝对值 超过塑性极限弯矩M_s绝对值时，判断船体梁进入反向塑性运动，求解梁进 入反向塑性变形后的运动幅值H₁'(t)和相对转动的转角α₁(t)：

$\left(\begin{array}{c}{\stackrel{··}{H}}_1^{\prime }\left(t\right)=\frac{3(-M_s)-(l{\xi }_3+3{\xi }_4)p\left(t\right)}{(l{\xi }_1+3{\xi }_2)m}\\ {\stackrel{··}{\alpha }}_1\left(t\right)=\frac{24({\xi }_1{\xi }_4-{\xi }_2{\xi }_3)p\left(t\right)-24\left({-M}_s\right){\xi }_1}{(l{\xi }_1+3{\xi }_2){\mathrm{ml}}^2}\end{array}\right);$

第八模块，用于当所述船体梁反向塑性变形达到最大即α₁(t)达到最大 值时，判断船体梁反向塑性变形阶段结束，进入正向弹性卸载，根据 $\left(\begin{array}{c}{\alpha }_1\left(t\right)(x-l/2)+{\phi }_1\left(x\right)H_1^{\prime }\left(t\right)={\phi }_1\left(x\right)H_1\left(t\right)\\ {\stackrel{·}{\alpha }}_1\left(t\right)(x-l/2)+{\phi }_1\left(x\right){\stackrel{·}{H}}_1^{\prime }\left(t\right)={\phi }_1\left(x\right){\stackrel{·}{H}}_1\left(t\right)\end{array}\right)$求解船体梁正向弹性卸载阶段的运动幅值 H₁(t)。

总体而言，通过本发明所构思的以上技术方案与现有技术相比，主要 具备以下的技术优点：本发明综合考虑冲击波和气泡载荷对舰船总体的冲 击作用，将水下爆炸载荷压力简化为五个阶段，通过建立相应的数学模 型，计算出舰船在不同阶段的运动位移，从而得到舰船速度、加速度以及 弯矩，并进一步根据上述参数分析舰船正向塑性变形、反向弹性卸载、反 向塑性变形和正向弹性卸载的运动过程。相比于现有技术，本发明可以较 为准确且简便地实现对水下近距爆炸作用下舰船弹塑性运动变形的工程预 报，对于提高舰船抗水下爆炸冲击防护设计水平、优化水中兵器攻击效能 及攻击方式等均具有指导借鉴意义。

附图说明

图1为本发明水下爆炸作用下舰船整体弹塑性运动响应预测方法流程 图；

图2为本发明船体梁爆炸工况示意图；

图3为本发明水下爆炸载荷五个阶段示意图；

图4为本发明理想弹塑性梁反复加载、卸载过程中的应力应变关系示 意图；

图5为本发明一个实施例中梁模型结构形式及尺寸示意图；

图6为本发明一个实施例中全弹性运动条件下船体梁中点弯矩时程曲 线示意图；

图7为本发明一个实施例中船体梁发生弹塑性运动时的中点位移时程 曲线示意图；

图8为本发明一个实施例中不同时刻船体梁长度方向变形情况比较示 意图；

图9为本发明一个实施例中不同时刻船体梁塑性运动过程中转角α变 化时程曲线示意图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图 及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体 实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的 本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可 以相互组合。

本发明提供了一种水下爆炸作用下舰船整体弹塑性运动响应预测方 法，如图1所示，方法首先将舰船运动等效为自由船体梁运动，综合考虑 冲击波和气泡两者对舰船的影响，将水下爆炸作用下自由船体梁载荷压力 划分为五个阶段，并建立相应的数学模型，进而求解出不同阶段内舰船的 弹性运动位移，根据位移求解各阶段内船体梁运动位移、速度、加速度及 弯矩，并根据上述参数进一步求解分析梁正向塑性运动过程、梁反向弹性 卸载过程、梁反向塑性变形过程和梁正向弹性卸载过程的运动响应。其具 体实现思路为：

步骤一，按照总纵强度等效及相似原则将舰船简化为等截面自由船体 梁；

步骤二，将水下爆炸载荷压力曲线简化为五个阶段，其中，第Ⅰ、Ⅱ 阶段压力满足线性变化规律，第Ⅲ阶段满足正弦函数变化关系，第Ⅳ、Ⅴ 阶段满足线性变化规律；

步骤三，确定船体梁的固有振形函数；

步骤四，根据冲击波和气泡脉动压力载荷变化特点，确定各阶段内船 体梁长度方向上的压力分布函数；

步骤五，根据不同压力阶段船体梁长度方向上的压力分布函数，确定 船体梁运动变形的主坐标函数；

步骤六，根据不同压力阶段船体梁的振形函数和主坐标函数，确定该 阶段内梁的运动位移函数；

步骤七，根据初始运动条件及各阶段连续运动条件，求解各阶段内船 体梁运动位移、速度、加速度及弯矩；

步骤八，根据弹性运动阶段梁内弯矩分布情况，确定梁最初进入塑性 变形的时刻，合理选取梁进入塑性变形后的位移函数，并求解正向塑性运 动过程；

步骤九，根据梁正向塑性运动过程，确定其发生反向弹性卸载的时刻， 结合弹性和塑性阶段的运动方程以及连续条件，获得弹性卸载阶段梁的运 动过程；

步骤十，根据梁反向弹性卸载过程的弯矩分布，确定梁进入反向塑性 变形的时刻，利用连续边界条件和弹塑性运动控制方程，确定反向塑性运 动过程；

步骤十一，根据梁反向塑性变形最大值确定其进入正向弹性卸载的时 刻，利用连续边界条件和弹塑性运动控制方程，确定正向塑性运动过程。

以下对上述步骤的具体实现做详细说明。

当炸药在舰船中部正下方爆炸时，其运动响应相对明显。本发明一个 具体实施例中，以该典型工况为研究对象。如图2所示，在此工况条件下， 舰船在水下近距爆炸冲击波和气泡联合作用下的整体弹塑性运动响应预测 方法包括如下步骤：

步骤一，进行舰船结构的船体梁近似。

将舰船简化为等截面自由船体梁，相比于实尺度船体梁和缩比船体模 型而言，其简化过程满足以下两个原则：

(a)保持原型和模型总纵惯性矩几何相似；

(b)舰船结构和炸药按相同尺度几何缩比后，保持模型一阶湿频率与 缩比药量的爆炸气泡脉动频率吻合。

步骤二，将水下爆炸载荷压力曲线简化为五个阶段

将水下爆炸载荷压力曲线简化为如图3所示的Ⅰ～Ⅴ共五个阶段。其中， 第I、II阶段为冲击波载荷衰减阶段，其满足线性变化规律；第III阶段为 边界条件下气泡膨胀收缩运动形成的流场负压变化阶段，其满足正弦函数 关系；第Ⅳ、Ⅴ阶段为气泡收缩产生二次脉动压力的上升、下降阶段，其 满足线性变化规律，则列出各阶段压力载荷计算公式如下：

(1)第Ⅰ阶段：p(t)＝P_m·(1-t/k₁)，        0≤t＜t₁

(2)第Ⅱ阶段：$p\left(t\right)=P_m/e·(1-\frac{t-t_1}{k_2}),$       t₁≤t＜t₂

(3)第Ⅲ阶段：p(t)＝P_b·sinβ(t-t₂)          t₂≤t＜t₃

(4)第Ⅳ阶段：$p\left(t\right)=P_s·\left(\frac{t-t_3}{k_4}\right)$            t₃≤t＜t₄

(5)第Ⅴ阶段：$p\left(t\right)=P_s·(1-\frac{t-t_4}{k_5})$         t₄≤t＜t₅

式中：$\left(\begin{array}{cc}{P}_m=K_1·{\left(\frac{{m_e}^{1/3}}{R}\right)}^{A_1}& k_1=\frac{\mathrm{e\theta }}{e-1},\end{array}\right)$t₁＝θ，$\theta =K_2·{m_e}^{1/3}·{\left(\frac{{m_e}^{1/3}}{R}\right)}^{A_2},$k₂＝t₂-t₁， $t_2=(\frac{850}{{\bar{P}}_0^{0.18}}-\frac{20}{{\bar{P}}_0^{1/3}}+n)·\frac{r_0}{c},\overline{P_0}=\frac{P_0}{P_{\mathrm{atm}}},$P₀＝P_atm+ρ_wgH₀，$n=11.4-\frac{10.6}{{\bar{r}}^{0.13}}+\frac{1.51}{{\bar{r}}^{1.26}},$$\bar{r}=\frac{R}{r_0},\beta =\frac{\pi }{k_3},$k₃＝t₃-t₂，$t_4=T=2.11\frac{{m_e}^{1/3}}{{(P_0/{\rho }_{w}g)}^{5/6}},P_s=\frac{39\times {10}^6+24P_0}{\bar{r}},$$P_b=\frac{-3.1\times {10}^4{m_e}^{1/3}{(H_0+10)}^{2/3}}{R},$k₄＝t₄-t₃，k₅＝t₅-t₄，$k_4=k_5=3290\frac{r_0}{P_0^{0.71}}.$

其中，m_e为TNT装药当量，R为爆距，P_m为冲击波压力峰值，K₁、 K₂、A₁、A₂为冲击波常数，k₁、k₂、k₃、k₄、k₅分别为与五个压力阶段的持 续时间相关的参数，P₀为炸药处静水压力，为无量纲压力参数，P_b为气泡 脉动过程中的负压峰值，P_s为脉动压力峰值，为表征爆炸距离的无量纲参 数，θ为冲击波衰减常数，H₀为装药深度，r₀为装药半径，ρ_w为水的密度， c为水中声速，P_atm为大气压，g为重力加速度，T为气泡脉动周期，β为气 泡负压阶段压力函数角频率值，各参数物理量均采用国际单位制。

步骤三，确定水下爆炸作用下船体梁振动的固有振型函数

当水下爆炸载荷较小时，自由船体梁在爆炸压力P(x,t)作用下作弹性运 动，其运动控制方程为

$\mathrm{EI}\frac{{\partial }^{4}w}{\partial x^4}+m\frac{{\partial }^{2}w}{\partial t^2}=P(x,t)---\left(1\right)$

其一般解为

其中，x为自由船体梁任意点处的横坐标值，w为运动位移函数，E为 弹性模量，I为梁横截面惯性矩，m为考虑附连水的单位梁长度质量。

其中，是梁第i阶固有振形，H_i(t)是梁第i阶振形对应的主坐标。 为了确定梁的运动位移函数，首先应分别确定固有振形及对应的主坐标。

梁固有振形的一般表达式为

其中，μ为与梁固有频率相关的变量，l为梁长，ζ₁、ζ₂、ζ₃、ζ₄为常 数。

对于自由边界梁，两端弯矩和剪力为零，因此在x＝0和x＝l时满足如下 边界条件

将式(3)代入上述边界条件，推导可得有ζ₁＝ζ₃,ζ₂＝ζ₄，且

$\left(\begin{array}{c}{\zeta }_2(\mathrm{ch\mu }-\cos \mu )+{\zeta }_1(\mathrm{sh\mu }-\sin \mu )=0\\ {\zeta }_2(\mathrm{sh\mu }+\sin \mu )+{\zeta }_1\left(\text{chμ-cosμ}\right)=0\end{array}\right)---\left(5\right)$

对于常数ζ₁和ζ₂来说，上述方程组的一个特解要求矩阵的行列式系数 等于零，即

(chμ-cosμ)²-(sh²μ-sin²μ)＝0          (6)

上述方程的解确定了自由梁的固有频率。另外，梁的固有频率还满足

$\omega ={\left(\frac{\mu }{l}\right)}^2\sqrt{\frac{\mathrm{EI}}{m}}---\left(7\right)$

当μ取第一个非零解μ₁＝4.730时，梁的一阶振形固有频率 此时ζ₂/ζ₁＝(sinμ₁-shμ₁)/(chμ₁-cosμ₁)≈-1，则梁的一阶振形 函数可近似表示为

步骤四，确定水下爆炸载荷在船体梁长度方向的压力分布函数

对于水下爆炸压力函数P(x,t)而言，其满足P(x,t)＝p(t)·p(x)，其中p(t)为 梁中点处压力时程衰减曲线，p(x)为相对于梁中点的爆炸载荷压力分布特征 函数。

对于冲击波载荷衰减阶段(第I、II阶段)，梁上压力分布特征函数p_s(x) 可表征为

$p_s\left(x\right)=\frac{R^{1.13}}{{\sqrt{R^2+{(x-l/2)}^2}}^{1.13}}---\left(9\right)$

对于气泡脉动阶段(第Ⅲ-V阶段)，梁上压力分布特征函数p_b(x)可表 征为

$p_b\left(x\right)=(1-\frac{2x-l}{l})·\exp [-8{\left(\frac{2x-l}{l}\right)}^2+4{\left(\frac{2x-1}{l}\right)}^3]+{\left(\frac{2l-2x}{l}\right)}^{1.5}·\left(\frac{2x-l}{l}\right)---\left(10\right)$

分别确定冲击波和气泡阶段船体梁上的压力分布特征函数p_s(x)和p_b(x) 后，结合五个阶段内的压力时程变化公式p(t)，即可得压力函数P(x,t)。

步骤五，确定不同压力阶段内船体梁运动的主坐标函数

对于强迫振动条件下的主坐标H_i(t)，可以表示为

$H_i\left(t\right)=a_i\cos {\omega }_{i}t+b_i\sin {\omega }_{i}t+\frac{1}{{\omega }_i}{\int }_0^t{f}_i\left(\tau \right)\sin {\omega }_i(t-\tau )\mathrm{d\tau }---\left(11\right)$

式中为与广义质量对应的广义激振力，a_i、b_i为积 分常数，由初始条件(位移和速度均为零)和运动连续条件决定。

明确水下爆炸各阶段的压力分布特征函数p_s(x)后，为简化后续公式推 导的参数表达形式，特引入参数

步骤六，确定不同压力阶段内梁的运动位移函数

考虑到炸药在梁中部下方爆炸时主要激起梁的低阶运动响应，为简化 问题，假设此时船体梁主要呈现一阶运动模态。根据式(8)、(11)分别确 定振形函数(x)和其对应的主坐标函数H₁(t)后，可得到梁的近似位移函数 为

由于水下爆炸载荷分为5个不同阶段，每个阶段的振形函数保持 一致，而对应的强迫振动主坐标函数H₁(t)不同，按照不同阶段的压力计算 公式，可得其对应的运动位移函数为(依次分别为式(13)至式(17))

第Ⅰ阶段：

第Ⅱ阶段：

第Ⅲ阶段：

第Ⅳ阶段：

第Ⅴ阶段：

步骤七，确定不同压力阶段内梁的运动参数，根据运动参数预测船体 梁鞭状运动响应情况。

式(13)～(17)中时间t以各压力阶段的起始时刻为零点，系数由初始运动条件(位移和速度均为零)及各阶段连续运动条件确定， 参数φ₁根据不同阶段的压力分布特征函数确定。

各阶段的位移函数w(x,t)确定后，梁的速度、加速度以及弯矩可以通过 位移函数对时间t或变量x求偏导数得到。

步骤八，梁由弹性变形进入塑性变形后的正向塑性运动过程

随着水下爆炸强度的增大，船体梁将全弹性运动响应向弹塑性运动响 应转变，包含如下变形过程：弹性运动→正向塑性变形→反向弹性卸载→ 反向塑性变形→正向弹性卸载。对于理想弹塑性船体梁而言，其运动变形 过程以及应力应变关系如图4所示。

在水下爆炸冲击载荷的某一阶段，若船体梁中部弯矩超过塑性极限弯 矩时，假设其中部形成一个固定塑性铰，梁继续保持一阶运动变形的同时， 整个梁还绕着塑性铰发生相对转动。考虑船体梁的对称性，以右半船体梁 为研究对象，列出其位移函数如下

式中，l/2<x≤l，为一阶振形函数，H'(t)为与之对应的运动幅值 函数，α(t)为发生相对转动时的转角函数。

在梁初始弹性运动过程中，其中部的弯矩最大，当其最大弯矩值超过 塑性极限弯矩M_s的绝对值时，中部出现固定塑性铰，梁将由弹性阶段进入 塑性变形阶段，梁中点弯矩值保持为M_s不变。此时，梁的运动控制方程满 足

$\frac{{\partial }^{2}M}{\partial x^2}+m\frac{{\partial }^{2}w}{\partial t^2}=P(x,t)---\left(19\right)$

将式(18)代入上式，并对其进行一、二阶积分，结合利用边界条件： x＝l/2时，Q＝0，M＝M_s，Q为梁内的剪力，M为梁的弯矩；x＝l时，Q＝0， M＝0，可以分别得到函数H'(t)和α(t)的计算方程：

$\left(\begin{array}{c}{\stackrel{··}{H}}^{\prime }\left(t\right)=\frac{3M_s-(l{\xi }_3+3{\xi }_4)p\left(t\right)}{(l{\xi }_1+3{\xi }_2)m}\\ \stackrel{··}{\alpha }\left(t\right)=\frac{24({\xi }_1{\xi }_4-{\xi }_2{\xi }_3)p\left(t\right)-24M_s{\xi }_1}{(l{\xi }_1+3{\xi }_2){\mathrm{ml}}^2}\end{array}\right)---\left(20\right)$

其中

ξ₁＝λ₁(l)-λ₁(l/2)，

ξ₂＝λ₂(l)-λ₂(l/2)-λ₁(l)·l/2，

ξ₃＝η₁(l/2)-η₁(l)，

ξ₄＝η₂(l/2)-η₂(l)+η₁(l)·l/2。

上述公式中引入了如下积分函数

${\int }_0^{x}p\left(x\right)\mathrm{dx}={\eta }_1\left(x\right),$

${\int }_0^x{\int }_0^{x}p\left(x\right)\mathrm{dxdx}={\eta }_2\left(x\right),$

梁在发生弹塑性运动转变过程中，满足动量和能量平衡条件，其中动 量平衡条件为

$2{\int }_{l/2}^{l}m[{\stackrel{·}{\alpha }}_0(x-l/2)+{\stackrel{·}{H}}_0^{\prime }{\phi }_1\left(x\right)]\mathrm{dx}={\int }_0^{l}m{\phi }_1\left(x\right)\stackrel{·}{H}\left(t_s\right)\mathrm{dx}---\left(21\right)$

能量平衡条件为

$2{\int }_{l/2}^{l}0.5m{[{\stackrel{·}{\alpha }}_0(x-l/2)+{\stackrel{·}{H}}_0^{\prime }{\phi }_1\left(x\right)]}^2\mathrm{dx}={\int }_0^{l}0.5m{\phi }_1^2\left(x\right){\stackrel{·}{H}}^2\left(t_s\right)\mathrm{dx}---\left(22\right)$

式中t_s为塑性运动发生的初始时刻，为塑性运动发生时位移函 数的初始参数值。

将方程组(20)对时间t进行一、二次积分，结合梁在弹塑性运动转换 时的动量和能量平衡条件，可以求出任何压力时间段内梁塑性变形阶段的 位移及相对转角。

步骤九，梁反向弹性卸载过程

船体梁进入塑性运动后，当其中部塑性变形达到最大值时，梁两端绕 着固定塑性铰发生相对转动的转角α(t)也将达到最大值，此后梁开始反向运 动，结构应力释放，梁中点带动两端进入弹性卸载。

假设弹性卸载过程中，梁位移函数仍满足式(12)要求，根据塑性、 弹性转换时刻的位移、速度连续条件，求解船体梁反向弹性卸载过程的运 动幅值H(t)，方程如下：

$\left(\begin{array}{c}\alpha \left(t\right)(x-l/2)+{\phi }_1\left(x\right)H^{\prime }\left(t\right)={\phi }_1\left(x\right)H\left(t\right)\\ \stackrel{·}{\alpha }\left(t\right)(x-l/2)+{\phi }_1\left(x\right){\stackrel{·}{H}}^{\prime }\left(t\right)={\phi }_1\left(x\right)\stackrel{·}{H}\left(t\right)\end{array}\right)---\left(23\right)$

以梁中点为考察点，结合弹性和塑性阶段得到的运动方程，同时利用 以上连续条件，可以得到弹性卸载阶段梁的运动方程。需要注意的是，除 梁中点外，其它部分的运动位移需要在弹性变形基础上，叠加一个相对转 动位移。

步骤十，梁反向塑性运动过程

在弹性卸载及反向运动过程中，当梁中点弯矩绝对值超过塑性极限弯 矩M_s绝对值时，梁将进入反向塑性运动过程，中部出现固定塑性铰，梁继 续保持一阶运动变形的同时，其两端绕着塑性铰发生相对转动。假设反向 塑性运动过程中，梁的位移函数及运动控制方程仍满足式(18)、(19)的 要求，保持梁中点弯矩值为-M_s不变，利用连续边界条件，可以得到梁反向 塑性变形时的运动幅值H₁'(t)及相对转角α₁(t)：

$\left(\begin{array}{c}{\stackrel{··}{H}}_1^{\prime }\left(t\right)=\frac{3(-M_s)-(l{\xi }_3+3{\xi }_4)p\left(t\right)}{(l{\xi }_1+3{\xi }_2)m}\\ {\stackrel{··}{\alpha }}_1\left(t\right)=\frac{24({\xi }_1{\xi }_4-{\xi }_2{\xi }_3)p\left(t\right)-24\left({-M}_s\right){\xi }_1}{(l{\xi }_1+3{\xi }_2){\mathrm{ml}}^2}\end{array}\right)---\left(24\right)$

此处需要注意反向塑性变形阶段内，其压力载荷p(x,t)的变化情况，即 判断p(x,t)是否跨越了多个载荷阶段。根据爆炸载荷p(x,t)的实际情况，分阶 段确定运动控制方程。

步骤十一，梁正向弹性卸载过程

当船体梁中部反向塑性变形达到最大值时，梁两端绕着固定塑性铰发 生相对转动的转角α₁(t)也将达到最大值，此后梁开始正向运动，结构应力 释放，梁中点带动其它部分再次弹性卸载。

假设梁正向弹性卸载过程中，梁位移函数仍满足式(12)要求，根据 塑性、弹性转换时刻的位移、速度连续条件，仍可得形如式(23)的计算 方程，即船体梁正向弹性卸载阶段的运动幅值H₁(t)满足如下公式：

$\left(\begin{array}{c}{\alpha }_1\left(t\right)(x-l/2)+{\phi }_1\left(x\right)H_1^{\prime }\left(t\right)={\phi }_1\left(x\right)H_1\left(t\right)\\ {\stackrel{·}{\alpha }}_1\left(t\right)(x-l/2)+{\phi }_1\left(x\right){\stackrel{·}{H}}_1^{\prime }\left(t\right)={\phi }_1\left(x\right){\stackrel{·}{H}}_1\left(t\right)\end{array}\right)---\left(25\right).$

以梁中点为考察点，结合弹性和塑性阶段得到的运动方程，同时利用 以上连续条件，可以得到正向弹性卸载阶段梁的运动方程。需要注意的是， 除梁中点外，其它部分的运动位移需要在弹性变形基础上，叠加一个相对 转动位移。

在本发明一个具体实施例中，选取某梁模型为分析对象，其相关尺寸 参数为：梁长2.8m，宽0.3m，高0.08m，板厚1mm，梁塑性极限弯矩1.8e4 Nm，具体结构形式如图5所示。选取TNT药量5g、爆距0.7m的爆炸工况 来分析梁的整体运动响应过程。

图6给出了全弹性运动假设条件下船体梁中点弯矩时程曲线。可以看 出，在气泡膨胀过程中，梁中点弯矩值已超过塑性极限弯矩M_s，表明梁在 气泡膨胀初期就已进入塑性运动过程。梁进入塑性变形后，其后续运动响 应需要根据弹性卸载、反向加载情况，进一步明确是否存在反向塑性变形 和卸载过程。

图7给出了船体梁发生弹塑性运动时的中点位移时程曲线。可以看出， 该工况条件下船体梁仍呈现明显的升沉运动，且气泡溃灭阶段梁的变形快 速增大，中拱变形明显，超过初次中拱变形(3.8cm)和中垂弯曲变形值 (5.1cm)，呈现鞭状运动；相比于1m爆距工况，0.7m爆距下气泡运动造 成的梁中垂变形相对明显，这与低压流场分布及变化特性相关；梁的运动 响应经历了中拱塑性变形→弹性卸载及反向加载→中垂塑性变形→再次弹 性卸载及中拱塑性变形这一复杂过程。与图6对比可知，塑性变形导致梁 的往复运动周期增大，一阶响应频率降低。

图8给出了典型时刻船体梁长度方向变形对比情况。可以看出，梁总 体呈现一阶变形，但塑性变形的存在导致梁并不严格绕着固定的驻点发生 正反向运动；在37.1ms时刻，梁发生最大中垂塑性变形，中部可见比较明 显的塑性铰；对于某一特定时刻，梁两端的运动位移超过中点位移值。

以塑性相对转角α为研究对象，进一步分析梁的塑性变形过程。图9 给出了转角α变化时程曲线。可以看出，初次中拱变形阶段，梁塑性铰转 动速度相对平缓，最大转角为-1.7e-3；气泡膨胀后期至收缩阶段，梁中垂塑 性变形快速增大，其塑性铰由-1.7e-3转变为1e-2；气泡溃灭后，脉动压力 造成梁发生快速中拱塑性变形，至计算结束时，其塑性转角为-8e-3。

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