徐嘉啟,梅志遠,劉志華
(海軍工程大學艦船與海洋學院,430033,武漢)
聲吶自噪聲對聲吶工作有不利影響,并在一定程度上影響了艦艇的對敵探測能力。當艦艇航速達10 kn以上時,水動力噪聲將成為自噪聲的主要分量[1]。實船測試表明,低航速下聲吶自噪聲主要源自于機械噪聲和螺旋槳噪聲,高航速下聲吶自噪聲的主要源頭之一則是結(jié)構(gòu)流激振動噪聲[1]。邊界層壁面脈動壓力載荷是聲吶導流罩流激振動的源頭,因此對邊界層壁面脈動壓力的研究是聲吶導流罩結(jié)構(gòu)流激振動的基礎。在縮尺模型載荷及結(jié)構(gòu)振動試驗中,往往因為縮尺模型與實船雷諾數(shù)的差異較大,而需要借助人工轉(zhuǎn)捩技術(shù)以實現(xiàn)壁面脈動壓力載荷相似。
本文的研究對象是絆線上下游壁面脈動壓力,其中,絆線上下游邊界層屬于非平衡邊界層,壁面脈動壓力研究方法屬于大渦模擬方法,絆線則是激起湍流的人工轉(zhuǎn)捩裝置。
對于邊界層壁面脈動壓力的研究始于平板平衡湍流邊界層。Blake、Farabee、Keith等針對平板平衡湍流邊界層壁面脈動壓力功率譜進行了大量研究,對其特征有了系統(tǒng)、全面的認識[2-4]。近年來,研究方向逐漸向含壓力梯度的三維邊界層發(fā)展。McGrath、Cipolla等研究了非平衡邊界層中壓力梯度對湍流邊界層壁面脈動壓力功率譜的影響[5-6]。Magionesi、Salze、Hu、龐業(yè)珍等分別對艦首聲吶導流罩、坡度板、可變攻角翼型下方平板等非平衡邊界層壁面脈動壓力功率譜進行了試驗研究,結(jié)果表明,縱向逆壓力梯度與零壓力梯度情況相比,湍流邊界層壁面脈動壓力自功率譜平臺區(qū)較高[7-10]。
壁面脈動壓力的研究方法一般包括試驗法、直接數(shù)值模擬法和大渦模擬法,這三種方法借助功率譜估計得到脈動壓力功率譜。此外,通過功率譜半經(jīng)驗模型也可直接得到脈動壓力功率譜。
Abraham在水洞中采用陣列壓力傳感器研究了平板湍流邊界層壁面脈動壓力的一維波數(shù)-頻率譜[11]。Salze、Arguillat等在風洞中采用可旋轉(zhuǎn)的陣列壓力傳感器研究了湍流脈動壓力的二維波數(shù)-頻率譜[8,12]。
經(jīng)典的功率譜半經(jīng)驗模型有Corcos模型[13]和Chase模型[14]。后來,又發(fā)展出適用于水中低速湍流邊界層的波數(shù)-頻率譜改進型Corcos模型[15]以及自功率譜Goody模型[16]等。
Wang、張曉龍等采用大渦模擬研究了湍流邊界層壁面脈動壓力功率譜[17-18]。張曉龍等給出的平板湍流壁面脈動壓力自功率譜大渦模擬計算值[18]與Abraham的試驗值[11]誤差在6 dB以內(nèi)。
功率譜半經(jīng)驗模型大多針對平板平衡湍流邊界層壁面脈動壓力,試驗研究受儀器性能限制較大、周期較長。綜上,本文采用大渦模擬對絆線上下游邊界層進行數(shù)值模擬。
邊界層人工轉(zhuǎn)捩裝置有絆線、粗糙元或金剛砂等[19]。田永強等在某后掠機翼人工轉(zhuǎn)捩研究中指出,在低雷諾數(shù)工況下,采用粗糙帶轉(zhuǎn)捩可近似模擬高雷諾數(shù)工況下的邊界層自由轉(zhuǎn)捩[20]。王猛等在風洞中測量了帶鋸齒形絆線的翼型壁面脈動壓力,以脈動壓力均方根為轉(zhuǎn)捩判據(jù)[21]。Rengasamy等對多種絆線的轉(zhuǎn)捩效果進行了風洞試驗研究,以近壁面速度均方根沿壁面法向的分布規(guī)律為轉(zhuǎn)捩判據(jù)[22]。近年來,對含人工轉(zhuǎn)捩裝置的非平衡邊界層壁面脈動壓力研究相對較少。此外,前后臺階流動與絆線上下游流動有一定的相似性,臺階流的試驗及數(shù)值模擬研究[23]有一定的借鑒意義。
本文通過大渦模擬方法實現(xiàn)了平板以及人工轉(zhuǎn)捩邊界層流場的數(shù)值模擬,以壁面脈動壓力功率譜等統(tǒng)計學特性描述了邊界層人工轉(zhuǎn)捩現(xiàn)象,并以時均速度、壁面脈動壓力自功率譜、波數(shù)-頻率譜特征為判據(jù),對比了不同絆線的人工轉(zhuǎn)捩效果,驗證了絆線實現(xiàn)邊界層人工轉(zhuǎn)捩的可行性。
(1)
在統(tǒng)計學中,定常、均勻的湍流邊界層壓力場為寬平穩(wěn)、各態(tài)遍歷隨機過程,可研究均方值、功率譜等統(tǒng)計特性。根據(jù)能量分布域的不同,功率譜可分為自功率譜、空間-頻率譜、波數(shù)-頻率譜等。
(2)
式中:f為頻率,ω為圓頻率。由式(2)以及ω=2πf知,Φ(ω)=Φ(f)/2π。瞬時壓力試驗值、大渦模擬計算值均為離散時間序列。功率譜估計前預先對瞬時壓力序列做去均值運算,得到脈動壓力序列。采用Welch周期圖法,由脈動壓力的N點時間序列估計其單邊自功率譜,公式為
(3)
對于平板平衡湍流邊界層,歸一化的壁面脈動壓力自功率譜Goody模型公式為
(4)
統(tǒng)計定常、均勻的湍流,壁面脈動壓力自功率譜不隨空間位置變化,空間-頻率譜(互功率譜)與空間間隔有關(guān)[11]。若ξ為流向間距,η為展向間距,則空間-頻率譜Φ、相干函數(shù)Γ、相位角Θ、遷移速度Uc的定義分別為
Φ(ξ,η,ω)=Φ(ω)Γ(ξ,η,ω)eiθ
(5)
(6)
(7)
Uc(ξ,η,ω)=ωξ/Θ(ξ,η,ω)
(8)
二維波矢-頻率譜是空間-頻率譜的二維傅里葉變換,其公式為
(9)
式中:kx為流向波數(shù);kz為展向波數(shù)。
(10)
一維波矢-頻率譜簡稱波數(shù)-頻率譜,減去瞬時壓力離散時空序列的時均值,可得脈動壓力離散時空序列p′(xm,tn)。采用Welch周期圖法估計脈動壓力流向波數(shù)-頻率譜(雙邊),公式為
p′(xm,tn)e-i[(2πf)tn-kxxm]|2〉(2πMNCW2)-1
(11)
采用計算流體動力學軟件Star-CCM+進行數(shù)值模擬。采用RANS(ReynoldsaveragedNavier-Stokesequation)方法求解不可壓縮流體Navier-Stokes方程,從而計算流場速度及壓力的時均項,采用大渦模擬求解速度及壓力的脈動項。
采用大渦模擬對流場變量進行空間濾波,空間濾波后的速度場為
(12)
式中:G(x,x′)為空間濾波函數(shù),x=(x,y,z),x′=(x′,y′,z′);ui為任意坐標方向的瞬時速度,變量上方加“~”代表空間濾波后的結(jié)果。
空間濾波后的動量方程為
(13)
對平板平衡湍流邊界層進行大渦模擬,并進行網(wǎng)格無關(guān)性驗證,以此作為大渦模擬方法的驗證。
RLS算法引入0<λ> <1的遺忘因子,可使得離當前時刻較遠的采樣數(shù)據(jù)對統(tǒng)計量估計的作用減小,由此使得時變信道被有效跟蹤,此外,要使算法適用于不同的信道時變速率環(huán)境,可以通過調(diào)節(jié)的大小的方法來解決。λ
為進行網(wǎng)格無關(guān)性驗證,建立了稀、中、密3套網(wǎng)格,分別稱為網(wǎng)格方案a、b、c。3套網(wǎng)格方案對應的網(wǎng)格數(shù)、測試區(qū)網(wǎng)格尺寸及量綱一網(wǎng)格尺寸見表1,其中Δx、Δz分別為流向和展向網(wǎng)格尺寸,測試區(qū)域內(nèi)流向網(wǎng)格尺寸達最小值且均勻分布,對應著數(shù)值計算的最小采樣空間間隔,Δx+、Δz+分別為Δx和Δz量綱一形式,Δx+=uτΔx/ν,Δz+=uτΔz/ν。由于需要求解包括黏性底層在內(nèi)的整個湍流邊界層,因此近壁面網(wǎng)格必須足夠精細,3套網(wǎng)格的壁面首層網(wǎng)格厚度均保證y+<1,,y+為量綱一壁面距離,y+=uτy/ν,y為網(wǎng)格中心點與壁面間距,近壁面處網(wǎng)格厚度增長比為1.05。
表1 不同網(wǎng)格方案的網(wǎng)格尺寸及量綱一網(wǎng)格尺寸
流場入口為均勻速度入口,入口對應坐標原點,出口為壓力出口,其余為無滑移壁面邊界。先采用RANS方法及SSTk-ω湍流模型[26]生成初始的穩(wěn)態(tài)流場,并計算湍流邊界層時均參數(shù),再采用大渦模擬進行非定常計算。
對平板湍流邊界層進行模擬,設置自由流流速為4.6m/s,流體密度為999.04kg/m3,流體動力黏度為1.14×10-3Pa·s,流體運動黏度為1.14×10-6m2/s。由于時間步長為10-4s,設置最小流向網(wǎng)格尺寸為0.001 06m,因此CFL(CourantFriedrichLewy)數(shù)小于1。根據(jù)文獻[11],布置48個壓力測點,測點1距離速度入口1.63m,采樣空間間隔即測點流向間距Δxp=4.22mm,采樣波數(shù)ks=2π/Δxp=1 490rad/m,分析截止波數(shù)為745rad/m。本文的數(shù)值模擬采樣頻率為10kHz,分析截止頻率為5kHz;文獻[11]的采樣頻率為2kHz,分析截止頻率為1kHz。數(shù)值計算網(wǎng)格如圖1所示,網(wǎng)格較密的測試區(qū)對應計算域中部。
圖1 數(shù)值計算網(wǎng)格
圖2 測點1處的時均速度剖面
從圖2可得出結(jié)論:黏性底層(y+<10)、過渡層(10
平板湍流邊界層的邊界層厚度δ、位移厚度δ*、動量損失厚度θ以及壁面摩擦速度uτ的計算值和試驗值對比見表2,可見計算值與試驗值相對誤差較小。綜上,本文在平板平衡湍流邊界層數(shù)值模擬中的近壁面處理是合理的。
表2 平板湍流邊界層時均參數(shù)的計算值和試驗值
首先,研究網(wǎng)格精度對壁面脈動壓力自功率譜計算值的影響。3套網(wǎng)格方案的脈動壓力自功率譜計算值(測點48)、試驗值[11]對比如圖3所示,自功率譜采用混合變量尺度律歸一化,以壁面切應力τw作為壓力尺度因子,以邊界層外邊界速度Ue作為時間尺度因子。
圖3 不同網(wǎng)格方案的脈動壓力自功率譜計算值與試驗值
由圖3得出結(jié)論:隨著網(wǎng)格加密,計算值與試驗值愈發(fā)匯聚,中高頻段譜線變化尤為明顯,使得平臺區(qū)對應的最高頻率增大;對比網(wǎng)格方案b和c,可以看出網(wǎng)格加密造成的計算值差異不大,表明平板湍流大渦模擬時,方案c的網(wǎng)格尺寸是合理的。
其次,采用網(wǎng)格方案c,對壁面脈動壓力自功率譜大渦模擬計算值(測點1、24、48)、試驗值[11]、Goody模型計算值進行對比,自功率譜采用混合變量尺度律歸一化,結(jié)果如圖4所示。
圖4 脈動壓力自功率譜的大渦模擬計算值、試驗值、Goody模型計算值
從圖4得出結(jié)論:大渦模擬計算值與試驗值吻合良好;不同測點的自功率譜經(jīng)尺度律歸一化后,在中低頻、重疊區(qū)匯聚,在頻率稍高時有所發(fā)散,與Abraham的試驗研究結(jié)論一致[11]。以上結(jié)論表明本文方法合理有效。
圖5展示了網(wǎng)格精度對流向波數(shù)-頻率譜計算值的影響,功率譜幅值以聲壓級表示,參考值為1Pa2/(Hz·rad·m-1),為與文獻[11]對應,將波數(shù)-頻率譜的負波數(shù)區(qū)平移至正波數(shù)區(qū),且只選取1kHz以下的部分。
(a)網(wǎng)格方案a
(b)網(wǎng)格方案b
(c)網(wǎng)格方案c圖5 不同網(wǎng)格方案的流向波數(shù)-頻率譜計算值
從圖5得出結(jié)論:對于不同精度的網(wǎng)格,遷移脊峰值均在-20dB左右;隨著網(wǎng)格精度降低,高頻區(qū)功率譜降低,與不同網(wǎng)格精度下脈動壓力自功率譜計算值的趨勢相仿,其原因可能是隨著網(wǎng)格精度降低,較大尺度渦的模擬由亞格子渦模型完成,對渦能量耗散的建模愈發(fā)不準確,計算精度下降。
圖6是流向波數(shù)-頻率譜試驗值的等高線圖[11],功率譜幅值以聲壓級表示,參考值為1Pa2/(Hz·rad·m-1),與圖5c進行對比可知:大渦模擬計算值與試驗值的遷移脊斜率即遷移速度吻合良好,遷移脊峰值均在-20dB左右,誤差小于5dB。以上結(jié)論表明本文方法合理有效。
圖6 流向波數(shù)-頻率譜試驗值[11]
參考文獻[21-22]風洞試驗中方形、鋸齒形絆線的設計,以吉賓斯雷諾數(shù)為依據(jù),設計高度d=1mm的方形、鋸齒形絆線。在平板上距入口0.2m處,沿展向布置方形或鋸齒形絆線,方形絆線寬4mm,鋸齒形絆線如圖7所示。
圖7 鋸齒形絆線
測點1~11(p1~p11)位于絆線上游,測點1距入口0.15m,測點12~20(p12~p20)位于絆線下游,測點12距入口0.208m,測點流向間距為0.004m,測點與絆線的示意圖見圖8。對鋸齒形絆線,編號加“c”的測點位于“前凸”下游,不加“c”的測點位于“后凹”下游。
(a)方形絆線
(b)鋸齒形絆線圖8 測點與絆線的示意圖
設置3種絆線方案,方案1不使用絆線,方案2使用方形絆線,方案3使用鋸齒形絆線。近壁面首層網(wǎng)格厚度保證y+<1,網(wǎng)格厚度增長比為1.05。以平板湍流邊界層數(shù)值模擬的網(wǎng)格精度為基礎,參考臺階流動數(shù)值模擬研究[23]中的臺階下游網(wǎng)格劃分,劃分絆線方案的數(shù)值計算網(wǎng)格,各絆線方案對應的網(wǎng)格數(shù)、測試區(qū)網(wǎng)格尺寸及量綱一網(wǎng)格尺寸見表3。絆線附近網(wǎng)格見圖9。
設置自由流流速U0=1.5m/s,流向網(wǎng)格尺寸Δx=0.5mm,時間步長Δt=0.2ms,易得CFL數(shù)小于1,采樣頻率fs=50kHz,分析截止頻率為2.5
表3 不同絆線方案的網(wǎng)格數(shù)、網(wǎng)格尺寸和量綱一網(wǎng)格尺寸
(a)方形絆線 (b)鋸齒形絆線圖9 絆線附近的網(wǎng)格
kHz。計算流向波數(shù)-頻率譜時,設置絆線上下游壓力測點流向間距Δxp=2mm,采樣波數(shù)ks=2π/Δxp=3 140rad/m,分析截止波數(shù)為1 570rad/m。
當距離速度入口的流向距離l=0.2m時,當?shù)乩字Z數(shù)Rel=Uel/ν≈3×105,根據(jù)邊界層理論,此時邊界層開始轉(zhuǎn)捩,但未達到充分發(fā)展湍流。以邊界層外緣流速Ue為特征速度,絆線高度為特征長度的雷諾數(shù),Red=Ued/ν≈1 300,大于吉賓斯雷諾數(shù)826,此時絆線可立即引起轉(zhuǎn)捩[19],并且絆線下游邊界層動量損失厚度雷諾數(shù)Reθ=Ueθ/ν≈4 000>2 000,根據(jù)文獻[22],絆線產(chǎn)生的上游歷史記憶效應對下游基本無影響。
絆線上游邊界層時均速度剖面由RANS方法獲得,絆線下游時均速度剖面由大渦模擬計算瞬時值再取平均獲得。對無絆線、方形絆線上下游、使用鋸齒形絆線時“前凸”上下游、使用鋸齒形絆線時“后凹”上下游的時均速度剖面進行對比,并以黏性底層壁面律u+=y+和對數(shù)律層壁面律u+=2.44·ln(y+)+5.0作參照,結(jié)果如圖10所示。
由圖10得出以下結(jié)論。
(1)在方形絆線緊后方,自壁面至絆線高度(y+≈66)區(qū)域內(nèi)有回流區(qū),且此時外層尾跡區(qū)較平板平衡邊界層有較大抬升;下游是再附著區(qū)以及松弛區(qū),此時的時均速度剖面接近平板平衡邊界層的情況。
(2)與平板平衡邊界層時均速度剖面相比,鋸齒形絆線“前凸”下游的時均速度剖面更接近平板平衡邊界層的情況;鋸齒形絆線“后凹”的下游絆線高度處流速稍低,且外層尾跡區(qū)有較大抬升。
方形絆線、鋸齒形絆線上下游的邊界層縱向壓力梯度分布如表4所示。
(a)無絆線
(b)方形絆線
(c)鋸齒形絆線“前凸”上下游
(d)鋸齒形絆線“后凹”上下游圖10 時均速度剖面
測點l/m縱向壓力梯度/Pa·m-1方形絆線鋸齒形絆線30.158223.73261.1270.174888.581 109.07110.1906 444.485 650.16120.20825 099.536 127.27140.2166 115.1297.92160.224-558.05-727.96180.232-559.94-467.88200.240-319.82-317.56
由表4得出結(jié)論:方形絆線、鋸齒形絆線上游以及絆線緊后方為縱向逆壓力梯度,逆壓力梯度較大時,時均速度剖面外層尾跡區(qū)較平板平衡邊界層有明顯抬升;絆線下游逐漸變?yōu)榭v向順壓力梯度,順壓力梯度值較小,時均速度剖面接近平板平衡邊界層的情況;相較于同等高度的方形絆線,鋸齒形絆線緊后方的縱向逆壓力梯度值較小,且較快地轉(zhuǎn)變?yōu)榭v向順壓力梯度。
(14)
式中Reτ為壁面摩擦速度雷諾數(shù),Reτ=uτδ/ν。
由圖11得出以下結(jié)論:
(1)絆線下游的壁面脈動壓力均方根明顯高于絆線上游,并且高于Farabee公式經(jīng)驗值,表明絆線下游的湍流脈動能量較高;
(a)方形絆線
(2)在與文獻[5]工況的邊界層壓力梯度相近的情況下,雖然絆線下游邊界層的壁面摩擦速度雷諾數(shù)較低,但脈動壓力均方根仍然達到甚至超過了試驗值。
(b)鋸齒形絆線圖11 絆線上下游的壁面脈動壓力均方根
圖12展示了無絆線、方形絆線和鋸齒形絆線上下游壓力測點的壁面脈動壓力的自功率譜,功率譜以聲壓級表示,參考值為1Pa2/Hz。試驗測量時自由流流速為1.5m/s,測點距速度入口1.63m[11]。為了保證流速相同時,絆線上游邊界層仍處于轉(zhuǎn)捩區(qū),絆線布置的較為靠前,距入口0.2m,與試驗工況[11]相比,當?shù)乩字Z數(shù)較小。
從圖12得出以下結(jié)論。
(1)有絆線時,l=0.2m處上游的脈動壓力自功率譜在中低頻段有一突躍,這與轉(zhuǎn)捩區(qū)T-S波有關(guān),表明邊界層已開始轉(zhuǎn)捩;無絆線時,l=0.2m處下游的脈動壓力自功率譜中低頻段也有這種突躍,表明邊界層仍處在轉(zhuǎn)捩區(qū)。
(2)兩種絆線下游的脈動壓力自功率譜均包含明顯的平臺區(qū),具備充分發(fā)展湍流脈動壓力自功率譜的一般特征,初步證明了兩種絆線的人工轉(zhuǎn)捩效果。
(3)雖然數(shù)值計算時當?shù)乩字Z數(shù)低于文獻[11]試驗時的當?shù)乩字Z數(shù),但是絆線下游壁面脈動壓力自功率譜平臺區(qū)的計算值卻高于文獻[11]的試驗值。
(a)無絆線
(b)方形絆線上下游
(c)鋸齒形絆線上下游圖12 壁面脈動壓力自功率譜
(4)絆線下游離絆線越遠,脈動壓力自功率譜幅值越低,測點16與20流向間距為0.016m,其脈動壓力自功率譜幅值差達10dB。文獻[22]指出,絆線下游邊界層流態(tài)可能從充分發(fā)展湍流退化為層流。因此,絆線實現(xiàn)轉(zhuǎn)捩的有效區(qū)域仍有待進一步研究。
圖13展示了無絆線、方形絆線和鋸齒形絆線上下游壁面脈動壓力流向波數(shù)-頻率譜,將波數(shù)-頻率譜的負波數(shù)區(qū)平移至正波數(shù)區(qū),并選取1.2kHz以下的部分,功率譜以聲壓級表示,參考值為1Pa2/(Hz·rad·m-1)。
由圖13得出以下結(jié)論。
(1)絆線下游與平板湍流邊界層的壁面脈動壓力流向波數(shù)-頻率譜分布規(guī)律一致,包含遷移脊且遷移脊斜率dω/dk=Uc≈0.65U0,U0為自由流流速。這與絆線下游壁面脈動壓力均方根、自功率譜特征一同證明了絆線的轉(zhuǎn)捩效果。
(a) 無絆線x=0.2 m處 (b) 無絆線x=0.2 m處 的上游 的下游
(c)方形絆線上游 (d)方形絆線下游
(2)無絆線的壁面脈動壓力流向波數(shù)-頻率譜均不具有明顯的遷移脊,表明自然轉(zhuǎn)捩時,測試區(qū)邊界層流態(tài)仍未達充分發(fā)展湍流。
(e)鋸齒形絆線上游 (f)鋸齒形絆線下游圖13 壁面脈動壓力流向波數(shù)-頻率譜
采用大渦模擬以及功率譜估計,獲得了平板平衡湍流邊界層的壁面脈動壓力自功率譜、波數(shù)-頻率譜,進行了網(wǎng)格無關(guān)性驗證,并與試驗值對比驗證了數(shù)值計算方法的合理性與有效性。
研究了絆線下游的非平衡湍流邊界層壁面脈動壓力自功率譜、波數(shù)-頻率譜,通過與平衡湍流邊界層脈動壓力自功率譜、波數(shù)-頻率譜的對比,得出以下結(jié)論:
(1)設計了方形絆線、鋸齒形絆線,當絆線高度雷諾數(shù)達到吉賓斯雷諾數(shù)時,可以實現(xiàn)邊界層轉(zhuǎn)捩;
(2)與同流速平板湍流邊界層相比,雖然當?shù)乩字Z數(shù)較小,絆線下游壁面脈動壓力自功率譜平臺區(qū)仍然高出10dB左右,絆線下游稍遠處,平臺區(qū)數(shù)值有所降低;
(3)從絆線下游壁面脈動壓力自功率譜、波數(shù)-頻率譜特征來看,本文設計的等高的方形絆線、鋸齒形絆線的轉(zhuǎn)捩效果基本相當;
(4)鋸齒形絆線緊后方的回流稍弱,不存在太過劇烈的逆壓力梯度,為避開絆線緊后方的回流,在進行湍流邊界層流激結(jié)構(gòu)振動模型試驗時,宜將絆線置于測試區(qū)上游的一定位置。
未來可增加對照組進一步歸納湍流壁面脈動壓力功率譜級與絆線特征參數(shù)之間的函數(shù)關(guān)系,從而更加精確地預估絆線下游的邊界層壁面脈動壓力功率譜級。