郜冶,譚大力,李海旭,王金玲,劉長猛
(1.哈爾濱工程大學航天與建筑工程學院,黑龍江哈爾濱150001;2.海軍裝備研究院,北京100192;3.中國船舶工業系統工程研究院,北京100094)
大型水面艦船的機動性對軍事戰術的實施尤其對艦載機的起飛降落過程起著決定性作用,這就要求艦船在航行過程中必須確定相對于優勢風的正確航向和特定的前進速度,確保艦載機安全著艦[1]。目前大型艦船主要依靠安裝在船橋中央桅桿上的測風儀器測量風速,風速儀本身的測量精度符合要求,但由于受艦船甲板及其艦上島型建筑周圍渦流的影響安裝在渦流區的風速儀測得的風場與實際真風存在很大的誤差。在國外,M.L.Thiebaux利用風洞試驗研究了艦載風速儀的誤差校正問題[2];P.K.Taylor和 E.C.Kent等針對商用船,利用實驗和CFD方法研究了受船體影響其周圍風場產生的扭曲和變形對風速儀測量精度的影響,并與衛星測報數據進行了對比[3-4];Philippe L.Nacass采用 CFD計算船載風速儀周圍的氣流失真現象[5]。在國內,劉連吉等用實驗方法研究了船上風速儀測量風速與浮標測量風速的誤差,并提出了測量風速的改進意見[6-7];劉慧等將船測風場資料與QuikSCAT衛星遙感風場資料進行對比分析研究了船測風場偏差統計特征[8]。近年來郜冶等致力于艦船流場的CFD計算,文獻[9]中運用CFD計算空氣正向流過滑躍式艦船產生的甲板氣流場,討論分析不同網格劃分形式對計算結果的影響,通過對全尺寸的滑躍式甲板艦船的數值計算分析得到艦船甲板氣流場渦旋結構特征,并將不同湍流模型應用于流場計算,分析了艦船流場的特點。文獻[10]利用FLUENT軟件UDF接口將LK和MMK模型引入FLUENT進行流場計算改進了計算結果。
本文采用數值模擬方法對不同工況下風速儀周圍的氣流場進行計算分析,得出風速儀更加合適的安裝位置。如計算條件允許,理論上可以研究流體在任何條件下的運動,再現風速儀周圍的風場。
為便于計算,將風速儀安裝位置處物理模型(包括梯形柱臺、雷達、天線和其他附屬結構)進行合理簡化,忽略掉不影響流場特征的細微結構如毫米級的天線、柱臺上的小孔等,簡化后的模型如圖1所示。

圖1 計算模型Fig.1 The calculation model
計算區域X軸表示船長方向,由船頭指向船尾為正方向,Y軸表示船寬方向,由左舷指向右舷為正方向,Z軸表示豎直方向,由下方指向上方為正方向,模型關于中心位置左右對稱。坐標原點位于梯形柱體底部中心處,現有方案風速儀安裝在如圖2中所示監測點位置,其坐標分別為D(0.0 m,-3.68 m,13.23 m)、D'(0.0 m,3.68 m,13.23 m)。

圖2 風速儀監測點位置俯視圖Fig.2 The location of the anemometers in the top viewport
艦船及其風速儀周圍的流場屬于三維復雜湍流場,為簡化計算設空氣為不可壓縮流體且符合Boussinesq假設,流動為穩態紊流,忽略固體壁面間的熱輻射。采用標準k-ε兩方程模型對氣流場進行數值模擬。建立流場穩定后的湍流控制通用方程[11-12]:

式中:div(ρVφ)表示對流項,div(Γgrad φ)表示擴散頂,S表示源項。控制方程包括連續方程如式(2)、動量方程如式(3)、能量方程如式(4)、k方程如式(5)和ε方程如式(6)。

式中:Gk是由平均速度梯度引起的湍動能k的產生項,μ是湍動粘度,可表示t成湍動能k和耗散率 ε 的函數,即 μt=ρCμk2/ε;C1ε、C2ε為經驗常數;σk、σε分別是湍動能k和耗散率ε對應的Prandtl數。在數值模擬過程中,各系數取為Cμ=0.09,C1ε=1.44,C2ε=1.92,σk=1.0,σε=1.3,σT=0.9。
流場計算區域取風速儀前為10倍柱臺底面邊長(柱臺底面邊長為5.75 m),風速儀后為20倍柱臺底面邊長,左右均為5倍柱臺底面邊長,垂向為5倍柱臺底面邊長,基于柱臺底面邊長為特征尺度的雷諾數量級為1.0×107。由于計算模型比較復雜采用非結構網格,對梯形柱臺、雷達、天線及附屬結構區域進行網格加密。全部網格由ANSYS ICEM生成,網格總數為8.77×106,計算區域及全流場網格如圖3所示,局部網格加密如圖4所示。

圖3 計算區域及全流場整體網格Fig.3 The calculation region and grids for the whole flow field

圖4 局部網格加密Fig.4 The refinement of local grid
計算中入口設為速度入口(具體速度方向大小隨不同計算工況改變),梯形柱體、雷達、天線及附屬結構所有表面均設定為無滑移壁面,其余邊界均設為自由滑移壁面,出口設為壓力出口。
采用文獻[11]中的方法驗證CFD計算方法的準確性,將相同計算方法及網格劃分方法的美國LHA級艦船的CFD計算結果與文獻[15]中的風洞實驗結果進行對比分析。LHA模型為1/48縮比模型,全長為 5.2 m,甲板寬為 0.75 m,網格數量約為 2.0×106,U0=6.858 m/s,使用標準k-ε 模型進行計算。驗證計算中,著艦點位置如圖5中所示。

圖5LHA著艦點的位置(m)Fig.5 The touchdown points of LHA(m)
圖6為CFD計算結果與文獻[15]中風洞試驗V-22旋翼通過艦點2的傾斜軌道上垂向速度對比。橫坐標為坐標值與甲板寬度d的比值,縱坐標為Z向速度W與入口速度U0的比值。圖6中數據對比顯示CFD計算結果與風洞試驗數據吻合較好,誤差不超過7%。證明論文中運用的網格劃分方法及標準k-ε模型可以用于預測船體繞流場的特征。

圖6 著艦點2傾斜軌道上垂向速度對比Fig.6 The vertical velocity of inclined orbit comparison at touchdown point 2
為討論不同船速、風速及不同風向角對風速儀測量位置處流場的影響,對現有安裝方案設置3種工況進行對比,當船速3 m/s,風速3 m/s(低航速低風速)時記為工況A;當船速14 m/s,風速3 m/s(高航速低風速)時記為工況 B;當船速14 m/s,風速14 m/s(高航速高風速)時記為工況C。并對每種工況設置4種不同的風向角進行計算,圖3中給出了風向角示意表示,從船體正前方吹風時設為是0°風,從船體左舷吹風時設是左舷風,由于模型關于中心位置左右對稱,因此只計算左舷風。計算時風向角分別選取 0°、10°、20°、30°。將每種工況左右監測點D、D'處數據進行提取,分離出此處風向角α'及風速υ',并與無窮遠處真風風向角α及風速υ進行對比,得到風向角差值 Δα和風速差值的模,其計算方法如圖7所示。現定義風向角差值Δα方向:與實際真風相比,監測點真風逆時針偏轉為正,順時針偏轉為負。

圖7 風向角差值和風速差值計算示意圖Fig.7 The calculation schematic diagram for wind angle and vector difference
圖8和圖9分別為工況A、B、C在不同風向角α時α'及e(υ')計算結果,其中e(υ')表示風速值誤差。風速儀在現有安裝位置得到的α'及e(υ')均較大,其中工況A與工況C誤差相差不大,而工況B中α'及e(υ')遠遠高于工況A和工況C,即低風速高船速時風速儀的測量誤差急劇升高;隨α值增大,各工況近風側D測點(左舷測點)α'誤差逐漸減小,遠風側D'測點(右舷測點)α'誤差逐漸增大,且遠風側誤差明顯小于近風側誤差;但隨風向角α增大,工況B及工況C中風速值誤差e(υ')減小,工況A在α=30°時風速值誤差e(υ')明顯增大,且遠風側風速值誤差始終略大于近風側誤差。

圖8 3種工況不同位置處α'計算結果Fig.8 The calculation result α'for the 3 conditions

圖9 3種工況不同位置風速值誤差計算結果Fig.9 The wind speed error for the 3 conditions
鑒于風速儀在現有安裝位置風向角及風速值測量誤差均較大,特設置風速儀不同安裝位置進行計算,試圖找到更合適的安裝位置以減小測量誤差。
不改變風速儀安裝高度,分別令x為0 m、-4.6 m 時,y為±3、±4、±5、±6、±7,±8 m,即將風速儀安裝位置在x=0 m與z=13.23 m及x=-4.6 m與z=13.23 m交線上分別由內向外進行平移。風速儀監測點分布線如圖2中所示。
3種工況都分別對α為0°和30°時,風速儀安裝在不同測點時的流場風速進行計算。圖10和11分別為3種工況在不同風向角α值時不同測點的α'計算結果,圖12和13分別為3種工況在不同風向角α值時不同測點的風速值誤差e(υ')。

圖10 α為0°時3種工況不同位置處α'計算結果Fig.10 α'calculation result for the 3 conditions in 0°

圖11 α為30°時3種工況不同位置處α'計算結果Fig.11 α'calculation result for the 3 conditions in 30°

圖12 α為0°時3種工況不同位置風速值誤差計算結果Fig.12 Wind speed error for the 3 conditions in 0°

圖13 α為30°時3種工況不同位置風速值誤差計算結果Fig.13 Wind speed error for the 3 conditions in 30°
計算結果顯示隨著風速儀安裝位置由內向外移動,風向角α'及風速值誤差e(υ')均隨y值絕對值(風速儀監測點距中心的距離)增大而減小,即將風速儀安裝位置向外移動有助于減小測量誤差。在遠風側,3種工況x=-4.6 m處的α'誤差普遍略大于x=0處,且隨α增大2位置處的α'差距加大。在近風側,工況A和C在x=-4.6 m處的α'誤差小于x=0處,且隨風向角α增大兩位置處的α'差距加大;而工況B在x=-4.6 m處的風向角α'誤差略大于x=0 m處,且隨風向角α增大兩位置處的α'差距減小。對于風速值,3種工況在x=-4.6m處誤差均明顯小于x=0處,且隨風向角增大兩位置處風速值誤差差距加大。即將風速儀安裝位置向前移到天線位置處可以明顯減小風速值測量誤差且風向角測量誤差不會顯著增大。
以上計算中僅考慮風速儀安裝位置的物理模型,忽略了船體結構的影響。為討論船體結構對風速儀測量位置處流場結構的影響,將船體模型和風速儀模型合并,如圖14所示。當船速為0,風速為5 m/s時為工況D,當船速為0,風速為20 m/s時為工況E,且分別對風速儀模型及合并模型進行計算。工況D對風向角為0°,工況E對風向角為0°和30°時,風速儀不同測點的流場進行計算。

圖14 船體和風速儀合并模型Fig.14 The merge model of ship and anemoscope
圖15~20分別為工況D、E風向角及風速值誤差計算結果。計算結果顯示2種工況風速儀單獨模型與合并模型測量誤差隨y值的變化趨勢是相同的,隨風速儀安裝位置向外延伸兩模型的風向角及風速值測量誤差差距均減小。

圖15 α為0°時工況D風向角差值α'計算結果Fig.15 α'calculation results for condition D in 0°

圖16 α為0°時工況E風向角差值α'計算結果Fig.16 α'calculation results for condition E in 0°

圖17 α為30°時工況E風向角差值α'計算結果Fig.17 α'calculation results for condition E in 30°

圖18 α為0°時工況D風速值誤差計算結果Fig.18 e(υ')calculation results for condition D in 0°

圖19 α為0°時工況E風速值誤差計算結果Fig.19 e(υ')calculation results for condition E in 0°

圖20 α為30°時工況E風速值誤差計算結果Fig.20 e(υ')calculation results for condition E in 30°
α=0°時,距離中心位置3 m以外,在左右舷對稱位置2種模型測得的風向角α'數值基本相等但方向相反,且2種模型測量結果的差值基本維持在2°以內;對于風速值測量誤差,在距離中心位置1~4 m,風速儀單獨模型測量結果明顯大于合并模型的測量結果,在距離中心位置4 m以外,合并模型測量的風速值誤差略微大于風速儀單獨模型的測量結果但兩模型的測量誤差在3%以內。
α=30°時,遠風側的風向角及風速值誤差測量結果均大于近風側的測量結果,近風側在4 m以外兩模型風向角測量值的差值在5°以內,風速值誤差的差值在2%以內,而在遠風側,2種模型的測量誤差差距相對較大。
1)風速儀在現有安裝位置測量誤差較大,其中低風速高船速時風速儀的測量誤差急劇升高。
2)隨吹風角度增大,近風側風向角測量誤差逐漸減小,遠風側風向角測量誤差逐漸增大,且遠風側風向角測量誤差明顯小于近風側風向角測量誤差。
3)隨吹風角度增大,風速值測量誤差變化不大,且遠風側風速值測量誤差始終略大于近風側風速值測量誤差。
4)將風速儀安裝位置向外平移有利于減小風向角和風速值測量誤差;將風速儀安裝位置向前移到天線位置處可以明顯減小風速值測量誤差且風向角測量誤差不會顯著增大。
5)船體本身會對風速儀附近的流場產生較大的影響。在此區域內,合并模型風速儀測量精度明顯較風速儀單獨模型低,隨風速儀安裝位置外移,2種模型的測量差距減小,尤其在近風側2種模型的測結果基本相同。此時可以使用風速儀單獨模型代替合并模型進行定性分析,但是定量分析時應該使用合并模型。
[1]宋劍,何建忠.航空母艦在隨機海況下的運動統計特性[J].艦船電子工程,2011,31(3):67-69.SONG Jian,HE Jianzhong.Statistical properties of the movement of aircraft carrier in the random sea conditions[J].Ship Electronic Engineering,2011,31(3):67-69.
[2]THIEBAUX M L.Wind tunnel experiments to determine correction functions for shipborne anemometers[J].Hydrog and Ocean Sci,1990,36(2):57-61.
[3]TAYLOR P K,KENT E C,YELLAND M J,et al.The accuracy of wind observations from ships[J].Advances in the Applications of Marine Climatology,1995,2(2):132-155.
[4]TAYLOR P K,KENT E C.The accuracy of meteorological observations from voluntary observing ships-present status and future requirements[C]//WMO,First Session of the CMM Subgroup on Voluntary Observing Ships.Athen,Geneva,1999:12.
[5]PHILIPPE L N.Shipborne wind measurements corrected for airflow distortion by computational fluid dynamics[R]//Brétigny-sur-Orge:Centred'Aviation Météorologique,2001:228-233.
[6]劉連吉,趙永平.船上氣溫、濕度和風速觀測的代表性[J].山東海洋學院學報,1981,11(3):24-31.LIU Lianji,ZHAO Yongping.The reliability of the wind speed air temperature and humidity observed on ship[J].Journal of Shandong College of Oceanology,1981,11(3):24-31.
[7]劉連吉.對現用海洋水文氣象要素調查規范的改進意見[J].海洋技術,1982,1(2):70-73.LIU Lianji.Suggestions on the instruction manual for obtaining oceanographic and meteorological data[J].Acta Oceanologica Sinica,1982,1(2):70-73.
[8]劉慧,胡松,鄒曉榮.船測資料與智利外海QuikSCAT風場比較分析[J].遙感技術與應用,2012,27(5):763-769.LIU Hui,HU Song,ZOU Xiaorong.Comparison between the fishing vessels and QuikSCAT scatterometer wind data of the offshore Chile[J].Remote Sensing Technology and Application,2012,27(5):763-769.
[9]劉長猛,郜冶.滑躍式甲板氣流場數值模擬[J].華中科技大學學報,2012,40(10):68-71.LIU Changmeng,GAO Ye.Numerical simulation of airwake on ski jump deck[J].Journal of Huazhong University of Science and Technology,2012,40(10):68-71.
[10]郜冶,劉長猛.護衛艦氣流場數值計算研究[J].哈爾濱工程大學學報,2013,34(1):1-5.GAO Ye,LIU Changmeng.Numerical simulation of frigate ship airwake[J].Journal of Harbin Engineering University,2013,34(1):1-5.
[11]LAUNDER B E,SPALDING D B.The numerical computation of turbulent flows[J].Comput Methods Appl Mech Eng,1974,80(3):269-289.
[12]COMINI G,GIUDICE S A.(k-ε )model of turbulent flow[J].Numer Heat transfer,1985,90(8):299-316.
[13]RAJAGOPALAN G,SCHALLER D,WADCOSK A J,et al.Experimental and computational simulation of a model ship in a wind tunnel[C]//43rd AIAA Aerospace Sciences Meeting.[S.l.],2005.