超空泡射弹结构响应的流固耦合仿真研究

张劲生，张嘉钟，王 聪，魏英杰

（哈尔滨工业大学 航天学院，哈尔滨 150001）

1 引 言

超空泡减阻是目前实现水下航行体超高速运动的新原理、新途径。利用超空泡减阻效应，水下航行体可以达到非常高的速度。超空泡射弹发射时会受到扰动，导致射弹在前进的同时，还绕着其头部来回转动，使射弹的尾部与空泡壁面发生连续的碰撞和反弹（尾拍）[1]。极大的轴向冲击载荷和剧烈的尾拍使得超空泡射弹的结构响应与传统的航行体有很大的不同。过大的动应力会导致射弹结构失效，过大的弯曲变形可以使射弹运动时发生翻滚，导致运动失稳[2]。因此研究超空泡射弹的结构响应具有重要意义，并为射弹的结构优化设计提供前提条件。

由于射弹空间体积狭小，同时极强的瞬态冲击载荷对传感器等测量器件又有潜在的破坏作用，因此用试验测试方法，研究射弹的结构响应的难度比较大[3]。运用理论分析和数值计算，文献[4-6]对存在尾拍时，航行体的动力学特性进行了研究，将航行体视为刚体，并对尾拍力的形式进行了简化处理；文献[7]将匀速运动的航行体作为弹性体，研究其结构响应，但对尾拍载荷的处理上与文献[4-5]是一样的。已有的对超空泡射弹结构响应的研究，都先对尾拍力采用高度简化的数学关系式，忽略弹体与水的耦合作用，将弹体视为刚体或者让弹体作匀速运动。

本文将弹体作为弹性体，并利用弹体与流体间的流固耦合作用来处理尾拍力，研究了依靠惯性作自由飞行的射弹的结构响应。先在Logvinovich独立膨胀原理基础上对流体域给出了一种简化计算模型，然后采用有限元法对射弹的结构响应进行仿真计算。并分析了射弹的长细比对其速度响应、动应力响应和弯曲变形的影响。

2 射弹模型与受力分析

射弹模型前部为平头圆锥体，平头直径为0.002 m，其后部为圆柱体。采用体积为3.4×10-6m3，但长细比不同的系列模型。材料为钢，弹性模量为1.93×1011Pa，泊松比0.3，密度7750 kg/m3，且只考虑弹性变形。

本文不考虑空泡内气体阻力和弹体重力的影响；弹体载荷包括前端面上的流体阻力FR和尾部的尾拍冲击力FI；在惯性作用下，弹体作减速自由飞行，同时发生尾拍。如图1所示。

FR方向沿模型轴向，其大小为[8]：FR=ρAcDv2/2，其中ρ，cD分别为流体的密度和阻力系数，A，v分别为模型前端面积和前进速度。模型的加速度a和惯性力FT由下式确定：FT=－FR=－ma。

由于运动的相对性，弹体模型的前端面中心点o看作静止不动[4]，在20 rad/s的初始扰动角速度的作用下，模型绕该中心在铅垂面内做上下摆动，尾部与空泡壁面不断发生碰撞。

3 流体域的计算模型

3.1 流体域的简化计算模型

由于多相流CFD计算以及空泡流与弹体非线性振动间的耦合计算耗时巨大，本文使用了流体域的简化计算模型，可以显著降低计算的时间耗费。

尾拍过程使流体域的内边界（空泡壁面）不断发生变化，使该处的边界条件很难设定。所以本文将整个流体域设定为多物质ALE算法，只需要对流体域的外边界进行边界条件的设定。将空泡的内部空间设为空白物质单元。如图2（a）所示。

根据Logvinovich的独立膨胀原理[9]，超空泡的每个横截面在其所在的平面内由中心几乎独立地向外扩张，而在射弹飞行方向（纵向）上几乎没有速度，这已经为试验所证实[10]。又考虑到尾拍的瞬时性，所以尾拍过程可以看作是运动着的弹体与静止的水面之间的拍击过程。也可以看作是，弹体模型做定轴转动时其尾部与高速前进的水面间的撞击。所以，通过圆筒形的运动流体来代替射弹周围的水流、用圆筒的内表面来代替空泡壁面（见图2中的a图），就完全可以模拟模型的尾拍过程。

3.2 流体域尺寸与边界条件

流体域长度0.74 m，内径0.022 m，外径0.2 m，其内径由超空泡半径的经验公式来确定[11]。可以算得本文中射弹尾部的空泡半径变化很小，故采用0.022 m作为流体域内径。

由于对称性，流体域只取一半来计算。对称面为对称边界条件，其他边界设定见图2（b）所示。边界速度的变化区间取为300～1400 m/s，其大小随时间变化，并按下式确定[4]：v（t）=v0/（1+ αv0t ），其中，v0，α分别为射弹的初始速度和系统常数。

4 数值方法

射弹模型采用拉格朗日算法，流体域采用ALE有限元算法，可以实现流固耦合的动态分析。ALE有限元矩阵方程[12]：

（1）质量守恒方程

式中：Mρ，Lρ和 Kρ分别为容量,转换和散度矩阵。

（2）动量守恒方程

式中：M和L分别为广义质量和传递矩阵，vr为对应于在参考构形描述下的速度；fint和fext分别为内力和外力向量。

本文采用LS-DYNA软件进行仿真计算。弹体采用全积分S/R六面体单元，算法为拉格朗日算法。流体域采用中心单点积分的带空白材料的单物质ALE单元，将空泡内部空间设为空白材料，使用欧拉算法。弹体与流体间的耦合采用欧拉—拉格朗日耦合算法。单元网格采用六面体结构化网格，固体单元与流体单元之间采用无侵蚀的罚函数耦合。为避免发生渗透现象对弹体尾部处的流体网格进行了局部加细。

5 计算结果与分析

射弹的长细比β分别取为：8.50、10.48、12.57、15.22、17.54和20。对弹体模型的速度响应、动应力响应和横向弯曲变形进行了数值计算，并对结果进行了比较分析。

5.1 射弹的速度响应

对弹体中部一节点的速度响应进行快速傅立叶变换，得到其频谱图，将频谱图中幅值最大的频率记为K。然后对不同长细比模型的速度响应的K值进行比较，可得到图3所示的模型长细比β与频率K的关系图。

从图中可以看出，随着β的增加，K逐渐降低，即弹体振动的频率变小了。β的增加使弹体变得更细长，抗弯刚度降低，使结构振动响应的周期也随之变长，从而使K值变小。

5.2 射弹的动应力响应

对射弹模型的第一主应力σ1、第三主应力σ3、最大剪切应力σs和Von Mises等效应力σv-m进行计算。结果显示，在每个单元内，这些应力的大小都近似于同步变化，并同时达到各自的最大值。对于弹体内的应力分布，由于弹体中部弯曲变形最显著，导致各应力的最大值也都出现在弹体中部位置，如图4所示。

该位置位于弹体中部的外表面上，通过应力分析可知，该处的应力状态可近似为单向应力状态，横截面近似为应力主平面。当该处处于拉伸变形时，σ1＞0和σ3=0；当处于压缩时，σ1=0和σ3＜0。仿真得到该位置的σ1随时间的变化如图5所示。从图中可以看出，在弹体上下摆动的一个周期内，有近一半的时间σ1＞0，另一半时间σ1=0，这与弹体的动态弯曲变形过程是一致的。并且，在运动的初期，第一主应力σ1达到最大值，然后迅速减小。在该位置上，应力、σs和σv-m的变化情况也与之类似。

对于长细比不同的模型，随着β的增加，弹体内σ1、、σs和σv-m的最大值的变化较复杂，忽高忽低。如图6所示。这主要是由于，长细比的改变不仅仅使弹体的变形能力发生了变化，而且对尾拍载荷也产生了很大的影响，这些影响因素的共同作用使最大应力与长细比β的关系变得复杂、不明确。但从图中可以看出，如果将β分段来看，在β不同的取值范围内，长细比与最大应力的关系还是可以确定的。

4.3 射弹的横向弯曲变形

建立射弹模型的局部坐标系，坐标原点位于前端面中心点O，坐标系随模型一起转动，X轴始终由O指向弹体尾端中心点，Z轴与对称面垂直。如图所示。则弹体的横向弯曲变形大小可用局部坐标系内Y方向的位移量Dy来衡量。

结果显示，弹体横向变形最大的位置在弹体的中部，如图4所示的n6-n8截面处。以β=10.48的弹体为例，当t=0.006 4 s时，Dy达到最大值，在该时刻，弹体轴线的弯曲变形情况如图7。对于其他长细比的模型，当出现最大弯曲时，弹体变形情况也是如此。

n6-n8截面处的Dy随时间的变化如图8所示。从图中可见，弹体的横向弯曲变形在运动开始不久就达到极值，然后很快衰减。这与弹体所受外载荷的变化情况相符合。

弹体内Dy的最大值与弹体长细比的关系如图9所示。可见，随着长细比增加，弹体变得更细长，从而使其弯曲变形有明显增大的趋势。可见，较细长的弹体飞行时的抖动更严重，这将使射弹的运动更加不稳定。

6 结 论

本文通过对射弹-流体的流固耦合仿真分析，得到如下结论：

（1）弹体的长细比变大将使其速度响应的频率变小，弹体结构振动的频率降低。

（2）弹体的第一主应力、第三主应力、最大剪切应力和Von Mises等效应力同时达到最大值，最大值出现在弹体中部。在运动的初期，各应力达到极值，然后迅速衰减。弹体长细比的变化对应力最大值的影响较复杂，应将长细比分段考虑。

（3）弹体横向变形最大的位置在弹体的中部。射弹运动开始不久横向变形达到最大值。随着长细比增加，弯曲变形程度明显增加。

[1]曹 伟,王 聪,魏英杰等．自然超空泡形态特性的射弹试验研究[J].工程力学,2006,23(12):175-179.

[2]Hrubes J D.High-speed imaging of supercavitating underwater projectiles[J].Experiments in Fluids,2001,30:57-64.

[3]吴三灵,温 波,于永强等.火炮动力学实验[M].北京:国防工业出版社,2004.

[4]Rand R,Pratap R,Ramani D.Impact Dynamics of a supercavitating underwater projectile[C]//Proceedings of the 1997 ASME Design Engineering Technical Conferences.Sacramento:ASME,1997.

[5]Kulkarni S S,Pratap R.Studies on the dynamics of a supercavitating projectile[J].Applied Mathematical Modelling,2000,24:113-129.

[6]孟庆昌,张志宏,顾建农等.超空泡射弹尾拍分析与计算[J].爆炸与冲击,2009,29(1):57-60.

[7]Ruzzene M,Soranna F.Impact dynamics of elastic supercavitating underwater vehicles[C]//9th AIAA/ISSMO Symposium on Multidisciplinary Analysis and Optimization.Atlanta:AIAA,2002.

[8]Kirschner I N,Kring D C,Stokes A W,et al.Control strategies for supercavitating vehicles[J].Journal of Vibration and Control,2002,8:219-242.

[9]Vasin A D.The Principle of independence of the cavity sections expansion(Logvinovich’s principle)as the basis for investigation on cavitation flows[C]//VKI Special Courses on Supercavitating Flows.Brussels:RTO-AVT and VKI,2001,RTO-EN-010(8).

[10]Semenenko V N.Artificial supercavitation:Physics and calculation[C]//VKI Special Courses on Supercavitating Flows.Brussels:RTO-AVT and VKI,2001,RTO-EN-010(11).

[11]Savchenko Y N.Experimental investigation of supercavitating motion of bodies[C]//VKI Special Courses on Supercavitating Flows.Brussels:RTO-AVT and VKI,2001,RTO-EN-010(4).

[12]李 政,金先龙,亓文果.流体-结构耦合问题的有限元并行计算研究[J].计算力学学报,2007,24(6):727-732.