張曉東

（海軍裝備部，北京 100841）

1 引 言

水下航行體在水下高速航行時，通常會在航行體后面產(chǎn)生空泡，而空泡的出現(xiàn)會極大地改變航行體所受到的流體動力。確定帶空泡或超空泡航行體的流體動力是一個復雜而又重要的流體力學問題。顯然，空泡流體動力必然與空泡的特性有關(guān)，因此對空泡動力學特性的研究是空泡流體動力研究的基礎。

在實際問題中，如何計算在壓力變化流場中的空泡長度是工程上最關(guān)心的問題之一，本文研究的水下航行體垂直空泡長度的計算就是這類問題。當航行體從水下朝水面高速運動時，受重力影響，不同水深位置的流體靜壓不同，導致不同水深的空泡所受到的外界流場壓力不斷發(fā)生變化，從而使得空泡長度也受到重力場的影響。影響空泡長度的另一關(guān)鍵物理特性就是空泡的記憶效應，它會導致超空泡自身的自激振動[1-2]，或帶空泡航行體的復雜動力學行為[4-6]。因此，當航行體從水下朝水面運動時，空泡必然受到重力場和記憶效應的相互耦合作用，導致水下航行體空泡長度的變化相當復雜。

建立空泡動力學的數(shù)學模型是研究空泡的關(guān)鍵。帶有記憶效應的系統(tǒng)通?？梢杂靡粋€或一組延遲微分方程（Delay Differential Equations，簡寫為 DDEs）來描述。Paryshev［1-2］在 Logvinovich［3］的獨立膨脹原理上建立了通氣空泡的DDEs數(shù)學模型。Kirschner［4］利用自己提出的算法通過求解Paryshev的模型獲得了航行體水平運動時的超空泡形態(tài)變化特征，但是利用Paryshev模型研究水下航行體在有限水深重力場中的垂直空泡還未見有文獻報導。顯然，航行體水平運動時，空泡軸線方向上并不存在重力梯度，而當航行體垂直運動時，沿空泡軸線方向的重力梯度及流體靜壓的變化增加了空泡變化的復雜性。

空泡模型建立后，對空泡DDEs模型的數(shù)值求解是另一個研究難點，求解DDEs方程并不能簡單地采用通常的常微分方程（ODEs）數(shù)值方法。因為對于DDEs方程，初始值必須嚴格地在一個時間間隔內(nèi)全部給出，而不是一個單點，這使得DDEs所描述的延遲系統(tǒng)對初始狀態(tài)非常敏感，并呈現(xiàn)出特殊的現(xiàn)象。Bellen和Zennaro[7]總結(jié)了DDEs在數(shù)學上的一些令人感興趣的特殊特征，例如：（1）DDEs的初始值必須在一個時間區(qū)間內(nèi)都指定，而不是只在一個點上指定。（2）DDEs的解不是光滑地與初始函數(shù)連接。（3）不同的初始值有可能導致同一個解。（4）對于依賴于狀態(tài)的延遲，初始函數(shù)的不規(guī)則會引起解的唯一性的喪失，或在有限次迭代后終止求解。5）延遲時間的介入，會強烈地改變ODEs算法的穩(wěn)定性。因此，Bellen和Zennaro指出了這樣的事實：即將許多常用的ODEs求解方法運用到DDEs中，將導致解的精度下降或失去穩(wěn)定性。

前文所述的Kirschner［4］提出的算法具有較好的魯棒性。這算法有兩個關(guān)鍵技術(shù)，（1）在計算過程中，采用分段三次Hermite插值多項式（piecewise-cubic Hermite interpolating polynomial，PCHIP）來描述解的時間歷程，PCHIP具有保形和消除偽振蕩的特點，在離散解的連續(xù)化重構(gòu)方面是一個很有吸引力的方法。（2）為了避免在每一時間步進行解的重構(gòu)以及相應的函數(shù)賦值，算法采用PCHIP重疊方式（PCHIPs）來構(gòu)造解的時間歷程。雖然重疊擬合方法在求解DDEs方程中顯示了它的優(yōu)點，但是它還是需要對所有歷史解進行擬合，這無疑會浪費計算時間。

針對上述問題，本文做了兩方面的工作，一是在Paryshev空泡模型中考慮了重力的影響；二是為了能快速而又準確地求解強非線性的復雜Paryshev空泡模型，對離散解的擬合重構(gòu)技術(shù)進行了改進，提出了一種效率更高的在延遲時間點進行局部擬合的方法。最后利用新方法求解了Paryshev空泡模型，計算了重力場中垂直空泡的長度變化。

2 空泡模型

針對軸對稱通氣空泡航行體（見圖1），Paryshev提出的空泡動力學模型的DDEs方程如下：

圖1 Paryshev空泡動力學模型的概念圖Fig.1 Paryshev’s dynamic model of cavity

其中，h=x+z，x是空化器到水面的距離。p0是大氣壓，ρ是水密度。再補充空化器到水面的距離x的方程

式中H是初始水深。

3 DDEs的數(shù)值特點

正如引言所述，求解空泡DDEs方程并不能簡單地采用常微分方程數(shù)值方法，DDEs方程有其特點。為了對比，這里先給出常微分方程（ODEs）的一般形式

其中，y（t）是t時刻的系統(tǒng)狀態(tài)，y（ t0）是初始時刻t0的初始狀態(tài)。

延遲微分方程（DDEs）的一般形式則為

其中，τi（i=1,2,…n）就是描述記憶效應的一組非負延遲變量。根據(jù)延遲變量τi的性質(zhì)，DDEs可分為三類：（1） 常數(shù)類 τi=const；（2） 依賴時間 t的類 τi=τi（t），（3） 依賴歷史狀態(tài) y（t）的類 τi=τi（t,y （t））。

在利用數(shù)值方法求解DDEs方程（2）時，假設在第k次迭代得到新的 τki（i=1,2,…n）值，然后需要求出對應的狀態(tài)量y（ tk-τki）（i=1,2,…n），因為 （tk-τki）不一定就是以前迭代的時間點tk（k=1,2,…），所以y（ tk-τki）需要通過對已求出的離散解（tk, y （ tk））k=1,2,…進行擬合得到。 這個擬合過程不僅給計算帶來了很多消耗，而且擬合帶來的誤差會傳播到后續(xù)的迭代過程中，因此擬合方法的好壞非常重要，采用何種擬合方法一般需要根據(jù)具體問題來定。

4 局部擬合技術(shù)

前文所述，Kirschner提出了一個利用PCHIP的重疊方式（PCHIPs）來構(gòu)造解的時間歷程，方法的原理圖見圖2，每條PCHIP擬合曲線最多只覆蓋固定數(shù)量的解，相鄰的兩條PCHIP擬合曲線之間在最后兩個區(qū)間上有重疊。

圖2 對DDE離散解進行連續(xù)化的PCHIP覆蓋方法Fig.2 The PCHIP covering method for continuous extension of the discrete DDE solution

PCHIPs的覆蓋擬合技術(shù)雖然能降低計算時間，但是PCHIPs還是需要對全部歷史區(qū)間進行覆蓋或擬合，實際上這是不必要的。仔細分析可以看出，在第k次迭代得到一個新的τi后，僅需要對tk-τi附近局部區(qū)間內(nèi)的幾個解進行擬合就可以求出狀態(tài)變量y（ tk- τi）（見圖 3）。

圖3 局部區(qū)間擬合方法示意圖Fig.3 The skech map of local fit

圖3顯示tk-τi落在離散時間的一個區(qū)間內(nèi)，然后選擇這個區(qū)間的相鄰幾個區(qū)間的解進行擬合就可以計算出y（ tk- τki），因此第一步就是需要判斷 （tk- τki）落在哪個離散時間區(qū)間內(nèi)。

經(jīng)驗表明，一般在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。

可以看出，這種擬合方法僅僅對需要計算的區(qū)間進行局部擬合，而對其他無關(guān)區(qū)間不做處理，因此這種局部區(qū)間擬合方法可以大幅度簡化計算過程，提高計算效率。

5 計算結(jié)果

由于在工程問題中最關(guān)心的是通氣空泡的長度變化，因此這里只計算通氣空泡長度l（t）的變化。為了簡化，這里認為空泡壓力pc（t）和航行體速度V（t）可以由試驗測量得到，這樣方程中與通氣率in和泄氣率Q˙相關(guān)的方程和計算都可以去掉。

從圖4可看出，航行體運動到不同水深位置的計算空泡長度與試驗測量結(jié)果一致，表明本文提出的計算非定常垂直空泡長度的方法是有效的。圖中=x/R=2位置上的尖峰是由于通氣裝置在通氣瞬間產(chǎn)生的壓力跳躍現(xiàn)象造成的。

圖4 局部擬合方法的計算結(jié)果與試驗對比Fig.4 Comparison of the length of unsteady cavity

6 結(jié) 論

航行體從水下到水面的高速運動過程中，沿空泡軸線的自由流體壓力不僅具有梯度，而且在不斷變化，這種重力導致的壓力梯度與空泡自身的記憶效應耦合，使得水下航行體垂直空泡呈現(xiàn)出特殊的非定常特點，計算這種非定常垂直空泡長度是工程上最為關(guān)鍵的問題之一。本文在Paryshev空泡模型的基礎上，增加了重力的影響項，提出了局部擬合方法，改進了求解空泡模型DDEs方程的數(shù)值計算技術(shù)，建立起一種可快速計算水下航行體非定常垂直空泡長度變化規(guī)律的方法，在工程上具有良好的應用前景。

[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.