于向陽(yáng)，姚凌虹，孟慶昌，劉巨斌，張志宏

(1. 海軍航空大學(xué)青島校區(qū)，山東 青島 266000；2. 海軍工程大學(xué)，湖北 武漢 430033)

0 引 言

潛艇水下操縱性能研究，與潛艇的安全航行密切相關(guān)，是其進(jìn)行戰(zhàn)術(shù)機(jī)動(dòng)的重要保證，同時(shí)有著直接的作戰(zhàn)應(yīng)用背景。實(shí)際航行中的潛艇是一種復(fù)雜運(yùn)動(dòng)的組合體。一方面，其高速、大攻角及強(qiáng)機(jī)動(dòng)的運(yùn)動(dòng)特征，使其呈現(xiàn)出非定常性與強(qiáng)非線性；另一方面，潛艇作有攻角機(jī)動(dòng)時(shí)，即使攻角保持不變，由于邊界層分離產(chǎn)生漩渦，流動(dòng)仍然會(huì)呈現(xiàn)出非定常特征。這就給理論和試驗(yàn)研究帶來(lái)了很大的困難，采用數(shù)值模擬是切實(shí)可行的研究方法[1–3]。

本文以DARPA2潛艇模型為研究對(duì)象，采用基于非結(jié)構(gòu)化網(wǎng)格的有限體積法，數(shù)值求解非定常、不可壓縮的RANS方程，湍流模型采用SA模式，在SA模式中，速度-壓力耦合SIMPLE方法處理，對(duì)流項(xiàng)采用2階迎風(fēng)格式，非定常項(xiàng)采用2階精度離散，隱式迭代求解方程組，非穩(wěn)態(tài)迭代時(shí)間步長(zhǎng)=2.322 4×10–3s（無(wú)量綱時(shí)間步長(zhǎng)=3.73×10–2），步進(jìn)200歩，非穩(wěn)態(tài)流場(chǎng)速度由UDF給定，得到非定常粘性流場(chǎng)和水動(dòng)力數(shù)值計(jì)算結(jié)果，并與文獻(xiàn)[1]的實(shí)驗(yàn)結(jié)果進(jìn)行對(duì)比。

1 數(shù)值離散方法

采用CFD軟件數(shù)值求解DARPA2潛艇模型粘性流場(chǎng)和水動(dòng)力，采用基于非結(jié)構(gòu)化網(wǎng)格的有限體積法，數(shù)值求解非定常、不可壓縮的RANS方程。

湍流模型采用SA模式，在SA模式中，生成項(xiàng)所采用的形式為：

速度-壓力耦合采用SIMPLE方法處理，對(duì)流項(xiàng)采用2階迎風(fēng)格式，非定常項(xiàng)采用2階精度離散，隱式迭代求解方程組。計(jì)算雷諾數(shù)Re=5.52×106，初值收斂準(zhǔn)則為所有求解變量的無(wú)量綱殘差下降4個(gè)量級(jí)，非穩(wěn)態(tài)迭代時(shí)間步長(zhǎng)=2.322 4×10–3s（無(wú)量綱時(shí)間步長(zhǎng)=3.73×10–2），步進(jìn)200歩，模型的運(yùn)動(dòng)狀態(tài)由UDF定義。

1.1 計(jì)算模型

計(jì)算對(duì)象為DARPA2潛艇模型（見(jiàn)圖1），模型主要參數(shù)如表1所示，包括艇主體、指揮臺(tái)圍殼、尾附體。數(shù)值計(jì)算在慣性坐標(biāo)系中進(jìn)行，采用直角坐標(biāo)系，原點(diǎn)取在模型頭部頂端，x軸與艇體的對(duì)稱(chēng)軸重合，方向與無(wú)攻角時(shí)來(lái)流方向一致，z軸豎直向上，y軸水平，計(jì)算時(shí)直角坐標(biāo)系保持為模型處于零攻角時(shí)的位置，當(dāng)模型作操縱運(yùn)動(dòng)時(shí)，坐標(biāo)不變，而模型與坐標(biāo)系產(chǎn)生相對(duì)運(yùn)動(dòng)；在所定義的坐標(biāo)系下， 規(guī)定來(lái)流的z方向速度分量正時(shí)為正攻角，初始攻角為0°。

圖 1 全附體DARPA2模型Fig. 1 DARPA2 submarine model with sail and stern appendages mounted

表 1 DARPA2潛艇模型主要參數(shù)Tab. 1 DARPA2 model scale

考慮到模型本身的對(duì)稱(chēng)性，同時(shí)操縱運(yùn)動(dòng)僅為xz面內(nèi)的非定?；剞D(zhuǎn)運(yùn)動(dòng)，在不影響數(shù)值計(jì)算的準(zhǔn)確性前提下，僅對(duì)半艇體流動(dòng)進(jìn)行數(shù)值計(jì)算，以減少計(jì)算量；計(jì)算域?yàn)榘雸A柱體區(qū)域，范圍：–2.0≤x/L≤5.2，–2≤z/L≤2，0≤y/L≤2，其中 L為艇體長(zhǎng)度，L=2.236 m，計(jì)算域如圖2所示。

圖 2 模型的邊界條件設(shè)置Fig. 2 Boundary setting of the computational domain

1.2 邊界條件

邊界條件設(shè)置可分為初始流場(chǎng)狀態(tài)和非穩(wěn)態(tài)運(yùn)動(dòng)狀態(tài)2個(gè)階段。

初始邊界條件（見(jiàn)圖2）：半圓柱體的左半圓面AEB（x/L=–2.0處）定義為速度入口，給定速度大小和方向；半圓柱體的右半圓面DFC（x/L=5.2處）定義為壓力出口；半圓柱面的上、下兩部分（EBCF/AEFD），及主體、指揮臺(tái)圍殼、尾附體，定義為無(wú)滑移壁面，速度分量和湍動(dòng)能為0，以加速收斂，獲得較好的初值；對(duì)稱(chēng)面上，法向速度分量為0，平行于對(duì)稱(chēng)面的速度分量的法向?qū)?shù)和所有標(biāo)量的法向?qū)?shù)為0；初始攻角為0。

非穩(wěn)態(tài)時(shí)模型作操縱運(yùn)動(dòng)，攻角發(fā)生改變，變化范圍由0°～–25°，須調(diào)整邊界條件，以適應(yīng)運(yùn)動(dòng)狀態(tài)的變化，此時(shí)定義半圓柱面的上半部分EBCF（z＞0）為速度入口；半圓柱面的下半部分AEFD（z≤0）為壓力出口，考慮到在迭代過(guò)程中，有可能存在數(shù)值上的反向流動(dòng)，設(shè)置出口回流方向的方式為垂直于出口面；由于采用慣性坐標(biāo)系，體網(wǎng)格是運(yùn)動(dòng)的，定義體網(wǎng)格和物面旋轉(zhuǎn)中心（0.659 62，0，0）和旋轉(zhuǎn)軸y軸，角速度由UDF定義。角速度曲線，具體描述如如圖3及表2所示（實(shí)驗(yàn)數(shù)據(jù)來(lái)自文獻(xiàn)[1]）。

1.3 網(wǎng)格設(shè)置

體網(wǎng)格采用基于四面體的非結(jié)構(gòu)化網(wǎng)格形式，邊界層網(wǎng)格進(jìn)行加密處理，網(wǎng)格密度由模型表面至遠(yuǎn)場(chǎng)由密到疏分布。

圖 3 角速度ω隨時(shí)間變化曲線Fig. 3 ω vs. time

表 2 不同模型的角速度函數(shù)Tab. 2 The function of ω for different model

表 3 網(wǎng)格主要參數(shù)Tab. 3 The grids for simulation

考慮到對(duì)模型進(jìn)行非穩(wěn)態(tài)運(yùn)動(dòng)模擬時(shí)，網(wǎng)格動(dòng)態(tài)變化，網(wǎng)格需要均勻銜接，以避免網(wǎng)格的奇異變化導(dǎo)致迭代不收斂。

2 動(dòng)態(tài)邊界處理

本文中應(yīng)用動(dòng)網(wǎng)格技術(shù)處理含有動(dòng)態(tài)邊界問(wèn)題的非穩(wěn)態(tài)運(yùn)動(dòng)狀態(tài)，邊界運(yùn)動(dòng)由UDF定義。

2.1 UDF編寫(xiě)

用戶自定義函數(shù)UDF（User-Defined Function）是CFD軟件提供的一個(gè)用戶接口，使用者可以通過(guò)它與CFD軟件的內(nèi)部數(shù)據(jù)進(jìn)行交流，方便用戶解決一些標(biāo)準(zhǔn)的CFD軟件模塊不能解決的問(wèn)題[4–6]。

圖 4 攻角–25°物面（model-1）Fig. 4 for the wall at –25° （model-1）

圖 5 攻角–25°物面 （model-2）Fig. 5 for the wall at –25° （model-2）

本文考慮到非穩(wěn)態(tài)問(wèn)題的復(fù)雜性，以及計(jì)算機(jī)本身的運(yùn)行環(huán)境（已安裝Visual Studio C++ 開(kāi)發(fā)軟件），采用編譯的UDF定義其運(yùn)動(dòng)（非穩(wěn)態(tài)場(chǎng)的角速度函數(shù)），以提高執(zhí)行速度[7]。

2.2 網(wǎng)格更新

動(dòng)網(wǎng)格模型，即在每一個(gè)時(shí)間歩迭代之前，根據(jù)邊界或物體的運(yùn)動(dòng)和變形來(lái)更新和重新構(gòu)建計(jì)算域網(wǎng)格，從而達(dá)到求解非定常問(wèn)題的目的[8–10]。而邊界的形變和運(yùn)動(dòng)由UDF加以定義。計(jì)算中的網(wǎng)格更新情況如圖6所示。

圖 6 網(wǎng)格更新示意圖（model-1）Fig. 6 The view of mesh motion （model-1）

由圖6可知，艇體和附體的物面邊界層都參與了網(wǎng)格的重構(gòu)過(guò)程，網(wǎng)格質(zhì)量沒(méi)有受到影響，網(wǎng)格更新過(guò)程中沒(méi)有負(fù)體積和大長(zhǎng)寬比網(wǎng)格出現(xiàn)，從而保證了數(shù)值迭代的延續(xù)。

3 初始流場(chǎng)計(jì)算

為了得到較好的非定常結(jié)果，需要對(duì)初始流場(chǎng)進(jìn)行討論。本文中非定常因素，一方面，大攻角時(shí)，分離流動(dòng)導(dǎo)致伴隨渦的出現(xiàn)；另一方面，給定的模型運(yùn)動(dòng)狀態(tài)，即角速度的定義本身就是隨時(shí)間變化的函數(shù)。圖7給出不同計(jì)算參數(shù)和邊界條件下的初始流場(chǎng)變量收斂曲線。最初，分析初始時(shí)刻時(shí)，沒(méi)有過(guò)多地考慮零攻角的流場(chǎng)特征，將邊界條件半圓柱面的上、下兩部分（EBCF/AEFD）分別定義為壓力出口和速度出口，在迭代過(guò)程中流場(chǎng)變量波動(dòng)起伏，近似周期性變化，但收斂精度未達(dá)到要求，如圖7（a）所示；分析上述收斂精度持續(xù)不下的情況，可能是松弛因子過(guò)大的原因，將部分松弛因子降低30%，流場(chǎng)變量的波動(dòng)情況消失，但收斂精度未出現(xiàn)改觀，如圖7（b）所示；觀察到速度分量殘差未發(fā)生變化，追溯根源，可能是邊界條件有待優(yōu)化，在圖7（a）所定義的邊界條件的基礎(chǔ)上，將半圓柱面的上、下兩部分（EBCF/AEFD）分別定義為無(wú)滑移壁面，速度分量殘差迅速收斂，如圖7（c）所示；圖7（d）再次調(diào)整部分松弛因子，流場(chǎng)變量在925次迭代后達(dá)到收斂精度。

4 時(shí)間步長(zhǎng)選取

圖 7 初始流場(chǎng)變量收斂歷史Fig. 7 Convergence history of flowfied variables

初始流場(chǎng)由零攻角定常結(jié)果給出，進(jìn)行非定常計(jì)算時(shí)需給定時(shí)間步長(zhǎng)。依據(jù)文獻(xiàn)[1]中的模型實(shí)驗(yàn)，給定角速度隨無(wú)量綱時(shí)間t'的變化關(guān)系如圖8所示，無(wú)量綱時(shí)間定義為，無(wú)量綱時(shí)間取值范圍0～7.463 7時(shí)，對(duì)應(yīng)的有量綱時(shí)間為0～0.464 48，計(jì)算時(shí)采用有量綱時(shí)間進(jìn)行計(jì)算，非定常攻角范圍0～–25°，考慮到時(shí)間步長(zhǎng)對(duì)數(shù)值計(jì)算結(jié)果的影響，認(rèn)為10–3s數(shù)量級(jí)已滿足計(jì)算精度的要求。

圖 8 角速度ω（rad/s）隨無(wú)量綱時(shí)間變化曲線Fig. 8 ω（rad/s） vs. non-dimensional time

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

應(yīng)用表2中的3種模型曲線，對(duì)模型運(yùn)動(dòng)進(jìn)行定義。初始化流場(chǎng)后，進(jìn)行非定常數(shù)值模擬。

為追求曲線擬合的精確性，采用最小二乘法高階多項(xiàng)式擬合實(shí)驗(yàn)數(shù)據(jù)，如圖3（model-3）所示，當(dāng)n=18時(shí)，除波動(dòng)較復(fù)雜的拐點(diǎn)處略有差別外，整條曲線基本重合；但進(jìn)行非定常數(shù)值計(jì)算時(shí)，角速度給定值持續(xù)增加，在不到100歩的迭代過(guò)程中，迭代時(shí)間不到半數(shù)，殘差竟迅速增至102數(shù)量級(jí)，超出預(yù)定值近100倍，最后導(dǎo)致計(jì)算結(jié)果發(fā)散。分析原因，認(rèn)為是時(shí)間精度不夠的原因，將時(shí)間步長(zhǎng)下降一個(gè)量級(jí)后，結(jié)果仍舊如前所述。嘗試時(shí)間步長(zhǎng)的一味降低，勢(shì)必導(dǎo)致計(jì)算量的增大。

重新分析計(jì)算結(jié)果發(fā)散的原因，非穩(wěn)態(tài)問(wèn)題本身需要進(jìn)行空間和時(shí)間的離散，變量之間的耦合關(guān)系將導(dǎo)致非線性增強(qiáng)，同時(shí)高階曲線擬合容易造成較大的截?cái)嗾`差，綜合上述原因，應(yīng)用model-2，降低多項(xiàng)式階次擬合實(shí)驗(yàn)數(shù)據(jù)，以驗(yàn)證是否是截?cái)嗾`差過(guò)大導(dǎo)致結(jié)果發(fā)散；非定常結(jié)果穩(wěn)定，且角速度按給定值正常迭代。

但model-2對(duì)實(shí)驗(yàn)數(shù)據(jù)擬合效果不夠理想，需要尋找其他函數(shù)，進(jìn)行數(shù)值計(jì)算。本例采用model-1對(duì)實(shí)驗(yàn)數(shù)據(jù)進(jìn)行擬合，同時(shí)略去角速度變化較大的尾端，擬合效果較好，數(shù)值計(jì)算結(jié)果穩(wěn)定。

5.1 水動(dòng)力分析

圖9給出了模型進(jìn)行非定常數(shù)值計(jì)算時(shí)，升力系數(shù)CL隨時(shí)間的變化情況。圖中曲線unsteady，Quais-Steady+Time Lag和Quais-Steady為文獻(xiàn)[1]中的實(shí)驗(yàn)結(jié)果，曲線unsteady為某一次非定常試驗(yàn)結(jié)果，Quais-Steady+Time Lag為20次非定常試驗(yàn)結(jié)果的平均值，Quais-Steady為準(zhǔn)定常狀態(tài)下的試驗(yàn)結(jié)果。

圖 9 升力系數(shù)與時(shí)間的關(guān)系Fig. 9 Normal force development against time

從圖中可以看到，定義相同的運(yùn)動(dòng)狀態(tài)，unsteady的升力系數(shù)卻表現(xiàn)出較強(qiáng)的不確定性，Quais-Steady+Time Lag曲線在描述unsteady曲線時(shí)符合程度也存在較大的偏差。分析上述問(wèn)題產(chǎn)生的原因（參考圖3），本文在對(duì)實(shí)驗(yàn)所定義的運(yùn)動(dòng)曲線的模擬上也存在局限性；同時(shí)，不得不指出非定常試驗(yàn)本身所具有的不確定性因素，如風(fēng)洞中溫度、風(fēng)速的變化、及氣流對(duì)操縱運(yùn)動(dòng)裝置的影響（易使模型發(fā)生振動(dòng)）等，隨著時(shí)間的變化，流場(chǎng)中物理量的變化容易受到這些不確定因素的影響。

比較2種模型的數(shù)值計(jì)算結(jié)果，升力系數(shù)隨時(shí)間變化平穩(wěn)，在前0.25 s與實(shí)驗(yàn)值符合較好；模型model-1對(duì)升力系數(shù)的計(jì)算值在中間段（0.15～0.25 s）有較明顯的減小過(guò)程（參考圖3），這與model-1所定義的角速度曲線在中間段有一個(gè)負(fù)向加速過(guò)程相對(duì)應(yīng)；模型model-2在初始階段，升力系數(shù)由正轉(zhuǎn)負(fù)，是與其所定義的角速度值在初始階段大于0相一致的；要取得與實(shí)驗(yàn)一致的模擬效果，由實(shí)驗(yàn)中的不確定因素產(chǎn)生的誤差需要加以考慮。

5.2 粘性流場(chǎng)分析

圖10～圖13給出了分別采用2種模型得到攻角–25°時(shí)物面壓強(qiáng)系數(shù)和物面切應(yīng)力系數(shù)的數(shù)值計(jì)算結(jié)果，兩者在物面上的連貫性和均勻性均較好，在迎風(fēng)處出現(xiàn)極值，這與理論上此處對(duì)應(yīng)速度的最小值是相符合的。與定常攻角–25°（見(jiàn)圖14和圖15）時(shí)的結(jié)果相比，在變化趨勢(shì)和數(shù)值分布上一致。

圖16給出了攻角–25°時(shí)x/L=0.863截面速度向量圖（model-1），在背風(fēng)面可以清晰地看到有分離渦出現(xiàn)，反映了流場(chǎng)的粘性分離特性。

6 結(jié) 語(yǔ)

圖 10 攻角–25°物面壓強(qiáng)系數(shù)（mode-l）Fig. 10 Surface pressure cofficient at –25° （mode-2）

圖 11 攻角–25°物面壓強(qiáng)系數(shù)（mode-2）Fig. 11 Surface pressure cofficient at –25° （mode-1）

圖 12 攻角–25°物面剪切力系數(shù)（mode-l）Fig. 12 Surface friction cofficient at –25° （mode-l）

圖 13 攻角–25°物面剪切力系數(shù)（mode-2）Fig. 13 Surface friction cofficient at –25° （mode-2）

圖 14 攻角–25°物面壓強(qiáng)系數(shù)Fig. 14 Surface pressure cofficient at –25°

圖 15 攻角–25°物面切應(yīng)力系數(shù)Fig. 15 Surface friction cofficient at –25°

采用CFD軟件數(shù)值求解非定常的DARPA2潛艇模型粘性流場(chǎng)和水動(dòng)力，采用基于非結(jié)構(gòu)化網(wǎng)格的有限體積法，數(shù)值求解非定常、不可壓縮的RANS方程，湍流模型采用SA模式，在SA模式中，速度-壓力耦合采用SIMPLE方法處理，對(duì)流項(xiàng)采用2階迎風(fēng)格式，非定常項(xiàng)采用2階精度離散，隱式迭代求解方程組，非穩(wěn)態(tài)迭代時(shí)間步長(zhǎng)=2.322 4×10–3s（無(wú)量綱時(shí)間步長(zhǎng)=3.73×10–2），步進(jìn)200歩，非穩(wěn)態(tài)流場(chǎng)速度由UDF給定。采用編譯的UDF定義其運(yùn)動(dòng)（非穩(wěn)態(tài)場(chǎng)的角速度函數(shù)）。

為了得到較好的非定常結(jié)果，需要對(duì)初始流場(chǎng)進(jìn)行討論。通過(guò)調(diào)整邊界條件和部分松弛因子，初始流場(chǎng)變量在925次迭代后達(dá)到收斂精度。

圖 16 攻角–25°時(shí)x/L=0.863截面速度向量圖Fig. 16 Velocity vector diagram at x/L=0.863 section and –25°angle of attack against time

在初始流場(chǎng)穩(wěn)定的基礎(chǔ)上，應(yīng)用表2中的3種模型曲線，分別定義運(yùn)動(dòng)狀態(tài)，進(jìn)行非定常數(shù)值模擬。過(guò)分追求曲線擬合的精確性，采用高階多項(xiàng)式擬合實(shí)驗(yàn)數(shù)據(jù)導(dǎo)致將計(jì)算結(jié)果發(fā)散；降低多項(xiàng)式階次，得到穩(wěn)定的非定常結(jié)果；采用model-1對(duì)實(shí)驗(yàn)數(shù)據(jù)進(jìn)行擬合，擬合效果較好，數(shù)值計(jì)算結(jié)果穩(wěn)定。

要取得與實(shí)驗(yàn)一致的模擬效果，由實(shí)驗(yàn)中的不確定因素產(chǎn)生的誤差需要加以考慮。攻角–25°時(shí)物面壓強(qiáng)系數(shù)和物面切應(yīng)力系數(shù)的數(shù)值計(jì)算結(jié)果，在壁面上的連貫性和均勻性均較好，與定常攻角–25°時(shí)的結(jié)果在變化趨勢(shì)和數(shù)值分布上一致。