水下航行体非定常垂直空泡长度的计算

张晓东

（海军装备部，北京 100841）

1 引 言

水下航行体在水下高速航行时，通常会在航行体后面产生空泡，而空泡的出现会极大地改变航行体所受到的流体动力。确定带空泡或超空泡航行体的流体动力是一个复杂而又重要的流体力学问题。显然，空泡流体动力必然与空泡的特性有关，因此对空泡动力学特性的研究是空泡流体动力研究的基础。

在实际问题中，如何计算在压力变化流场中的空泡长度是工程上最关心的问题之一，本文研究的水下航行体垂直空泡长度的计算就是这类问题。当航行体从水下朝水面高速运动时，受重力影响，不同水深位置的流体静压不同，导致不同水深的空泡所受到的外界流场压力不断发生变化，从而使得空泡长度也受到重力场的影响。影响空泡长度的另一关键物理特性就是空泡的记忆效应，它会导致超空泡自身的自激振动[1-2]，或带空泡航行体的复杂动力学行为[4-6]。因此，当航行体从水下朝水面运动时，空泡必然受到重力场和记忆效应的相互耦合作用，导致水下航行体空泡长度的变化相当复杂。

建立空泡动力学的数学模型是研究空泡的关键。带有记忆效应的系统通常可以用一个或一组延迟微分方程（Delay Differential Equations，简写为 DDEs）来描述。Paryshev［1-2］在 Logvinovich［3］的独立膨胀原理上建立了通气空泡的DDEs数学模型。Kirschner［4］利用自己提出的算法通过求解Paryshev的模型获得了航行体水平运动时的超空泡形态变化特征，但是利用Paryshev模型研究水下航行体在有限水深重力场中的垂直空泡还未见有文献报导。显然，航行体水平运动时，空泡轴线方向上并不存在重力梯度，而当航行体垂直运动时，沿空泡轴线方向的重力梯度及流体静压的变化增加了空泡变化的复杂性。

空泡模型建立后，对空泡DDEs模型的数值求解是另一个研究难点，求解DDEs方程并不能简单地采用通常的常微分方程（ODEs）数值方法。因为对于DDEs方程，初始值必须严格地在一个时间间隔内全部给出，而不是一个单点，这使得DDEs所描述的延迟系统对初始状态非常敏感，并呈现出特殊的现象。Bellen和Zennaro[7]总结了DDEs在数学上的一些令人感兴趣的特殊特征，例如：（1）DDEs的初始值必须在一个时间区间内都指定，而不是只在一个点上指定。（2）DDEs的解不是光滑地与初始函数连接。（3）不同的初始值有可能导致同一个解。（4）对于依赖于状态的延迟，初始函数的不规则会引起解的唯一性的丧失，或在有限次迭代后终止求解。5）延迟时间的介入，会强烈地改变ODEs算法的稳定性。因此，Bellen和Zennaro指出了这样的事实：即将许多常用的ODEs求解方法运用到DDEs中，将导致解的精度下降或失去稳定性。

前文所述的Kirschner［4］提出的算法具有较好的鲁棒性。这算法有两个关键技术，（1）在计算过程中，采用分段三次Hermite插值多项式（piecewise-cubic Hermite interpolating polynomial，PCHIP）来描述解的时间历程，PCHIP具有保形和消除伪振荡的特点，在离散解的连续化重构方面是一个很有吸引力的方法。（2）为了避免在每一时间步进行解的重构以及相应的函数赋值，算法采用PCHIP重叠方式（PCHIPs）来构造解的时间历程。虽然重叠拟合方法在求解DDEs方程中显示了它的优点，但是它还是需要对所有历史解进行拟合，这无疑会浪费计算时间。

针对上述问题，本文做了两方面的工作，一是在Paryshev空泡模型中考虑了重力的影响；二是为了能快速而又准确地求解强非线性的复杂Paryshev空泡模型，对离散解的拟合重构技术进行了改进，提出了一种效率更高的在延迟时间点进行局部拟合的方法。最后利用新方法求解了Paryshev空泡模型，计算了重力场中垂直空泡的长度变化。

2 空泡模型

针对轴对称通气空泡航行体（见图1），Paryshev提出的空泡动力学模型的DDEs方程如下：

图1 Paryshev空泡动力学模型的概念图Fig.1 Paryshev’s dynamic model of cavity

其中，h=x+z，x是空化器到水面的距离。p0是大气压，ρ是水密度。再补充空化器到水面的距离x的方程

式中H是初始水深。

3 DDEs的数值特点

正如引言所述，求解空泡DDEs方程并不能简单地采用常微分方程数值方法，DDEs方程有其特点。为了对比，这里先给出常微分方程（ODEs）的一般形式

其中，y（t）是t时刻的系统状态，y（ t0）是初始时刻t0的初始状态。

延迟微分方程（DDEs）的一般形式则为

其中，τi（i=1,2,…n）就是描述记忆效应的一组非负延迟变量。根据延迟变量τi的性质，DDEs可分为三类：（1） 常数类 τi=const；（2） 依赖时间 t的类 τi=τi（t），（3） 依赖历史状态 y（t）的类 τi=τi（t,y （t））。

在利用数值方法求解DDEs方程（2）时，假设在第k次迭代得到新的 τki（i=1,2,…n）值，然后需要求出对应的状态量y（ tk-τki）（i=1,2,…n），因为 （tk-τki）不一定就是以前迭代的时间点tk（k=1,2,…），所以y（ tk-τki）需要通过对已求出的离散解（tk, y （ tk））k=1,2,…进行拟合得到。 这个拟合过程不仅给计算带来了很多消耗，而且拟合带来的误差会传播到后续的迭代过程中，因此拟合方法的好坏非常重要，采用何种拟合方法一般需要根据具体问题来定。

4 局部拟合技术

前文所述，Kirschner提出了一个利用PCHIP的重叠方式（PCHIPs）来构造解的时间历程，方法的原理图见图2，每条PCHIP拟合曲线最多只覆盖固定数量的解，相邻的两条PCHIP拟合曲线之间在最后两个区间上有重叠。

图2 对DDE离散解进行连续化的PCHIP覆盖方法Fig.2 The PCHIP covering method for continuous extension of the discrete DDE solution

PCHIPs的覆盖拟合技术虽然能降低计算时间，但是PCHIPs还是需要对全部历史区间进行覆盖或拟合，实际上这是不必要的。仔细分析可以看出，在第k次迭代得到一个新的τi后，仅需要对tk-τi附近局部区间内的几个解进行拟合就可以求出状态变量y（ tk- τi）（见图 3）。

图3 局部区间拟合方法示意图Fig.3 The skech map of local fit

图3显示tk-τi落在离散时间的一个区间内，然后选择这个区间的相邻几个区间的解进行拟合就可以计算出y（ tk- τki），因此第一步就是需要判断 （tk- τki）落在哪个离散时间区间内。

经验表明，一般在tk-τi前后各取3个时间点进行拟合就足够了，即取tk-Δt·（mki+3），tk-Δt·（mki+2），tk-Δt·（mki+1），tk-Δt·mki，tk-Δt·（mki-1），tk-Δt·（mki-2）6个点进行拟合，拟合多项式可以采用前文所述的PCHIP。

可以看出，这种拟合方法仅仅对需要计算的区间进行局部拟合，而对其他无关区间不做处理，因此这种局部区间拟合方法可以大幅度简化计算过程，提高计算效率。

5 计算结果

由于在工程问题中最关心的是通气空泡的长度变化，因此这里只计算通气空泡长度l（t）的变化。为了简化，这里认为空泡压力pc（t）和航行体速度V（t）可以由试验测量得到，这样方程中与通气率in和泄气率Q˙相关的方程和计算都可以去掉。

从图4可看出，航行体运动到不同水深位置的计算空泡长度与试验测量结果一致，表明本文提出的计算非定常垂直空泡长度的方法是有效的。图中=x/R=2位置上的尖峰是由于通气装置在通气瞬间产生的压力跳跃现象造成的。

图4 局部拟合方法的计算结果与试验对比Fig.4 Comparison of the length of unsteady cavity

6 结 论

航行体从水下到水面的高速运动过程中，沿空泡轴线的自由流体压力不仅具有梯度，而且在不断变化，这种重力导致的压力梯度与空泡自身的记忆效应耦合，使得水下航行体垂直空泡呈现出特殊的非定常特点，计算这种非定常垂直空泡长度是工程上最为关键的问题之一。本文在Paryshev空泡模型的基础上，增加了重力的影响项，提出了局部拟合方法，改进了求解空泡模型DDEs方程的数值计算技术，建立起一种可快速计算水下航行体非定常垂直空泡长度变化规律的方法，在工程上具有良好的应用前景。

[1]Paryshev E V.A system of nonlinear differential equations with a time delay,Describing the Dynamics of Non-Stationary,Axially Symmetric Cavities[M].Trudy,TsAGI,No.1907,Moscow,Russia;translated from the Russian,1978.

[2]Paryshev E V.Approximate mathematical models in high-speed hydrodynamics[J].Journal of Engineering Mathematics,1978,55:41-64.

[3]Logvinovich G V.Hydrodynamics of Free-Boundary Flows[M].Naukova Dumka,Kiev,Ukraine,1969.translated from the Russian by the Israel Program for Scientific Translations,Jerusalem,1972.

[4]Kirschner I N.A robust solver for delay differential equations[C]//Proceedings of the 2006 Conference on High-Speed Hydrodynamics&Numerical Simulation(HSH&NS2006).Kemerovo State University,Kemerovo,Russia,2006.

[5]Kirschner I N,Kring D C,Stokes A W,Fine N E,Uhlman J S.Control strategies for supercavitating vehicles[J].J Vibration and Control,2002,82.(Also presented at the Eighth International Symposium on Nonlinear Dynamical Systems,Blacksburg,VA.)

[6]Kirschner I N,Rosenthal B J,Uhlman J S.Simplified dynamical systems analysis of supercavitating high-speed bodies[C]//Cav03-OS-7-005,Proceedings of the Fifth International Symposium on Cavitation(CAV2003).Osaka,Japan,2003.

[7]Kirschner I N,Arzoumanian S H.Implementation and extension of Paryshev’s model of cavity dynamics[C]//SuperFAST’2008.Russion,2008.