張曉東

（海軍裝備部，北京 100841）

1 引 言

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

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

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

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

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

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

2 空泡模型

針對(duì)軸對(duì)稱(chēng)通氣空泡航行體（見(jiàn)圖1），Paryshev提出的空泡動(dòng)力學(xué)模型的DDEs方程如下：

圖1 Paryshev空泡動(dòng)力學(xué)模型的概念圖Fig.1 Paryshev’s dynamic model of cavity

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

式中H是初始水深。

3 DDEs的數(shù)值特點(diǎn)

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

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

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

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

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

4 局部擬合技術(shù)

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

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

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

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

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

經(jīng)驗(yàn)表明，一般在tk-τi前后各取3個(gè)時(shí)間點(diǎn)進(jìn)行擬合就足夠了，即取tk-Δt·（mki+3），tk-Δt·（mki+2），tk-Δt·（mki+1），tk-Δt·mki，tk-Δt·（mki-1），tk-Δt·（mki-2）6個(gè)點(diǎn)進(jìn)行擬合，擬合多項(xiàng)式可以采用前文所述的PCHIP。

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

5 計(jì)算結(jié)果

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

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

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

6 結(jié) 論

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

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