水翼非定?？栈鲃又型牧髂Ｐ脱芯?/h1>
2019-12-20 02:26:30李雨濛張亞太
火箭推進 2019年6期
李雨濛，陳 暉，項 樂，張亞太
(西安航天動力研究所 液體火箭發(fā)動機重點實驗室，陜西 西安 710100)
0 引言
近年來航天不斷發(fā)展，液體火箭發(fā)動機作為重要的動力來源也不斷更新?lián)Q代。新一代液體火箭發(fā)動機以液氫液氧或液氧煤油作為推進劑，比沖高，推力大，工作時間長，在大型運載火箭上得到廣泛應(yīng)用。渦輪泵作為液體火箭發(fā)動機的核心部件，在高轉(zhuǎn)速的環(huán)境下工作易導(dǎo)致渦輪泵內(nèi)出現(xiàn)空化，影響液體火箭發(fā)動機的可靠性[1-5]。
空化作為水力機械中經(jīng)常出現(xiàn)的一種現(xiàn)象，會造成水力機械性能明顯降低[6]，材料表面破壞，引起振動和噪聲等[7-8]?？栈瘜討B(tài)響應(yīng)特性的改變會使流動內(nèi)部出現(xiàn)不穩(wěn)定性，多個國家在研究液體火箭發(fā)動機渦輪泵時都碰到過空化不穩(wěn)定帶來的問題乃至事故。例如，在日本H-2火箭的第8次發(fā)射中，空化不穩(wěn)定誘發(fā)的脈動頻率與泵前導(dǎo)流葉片固有頻率相近，引起共振導(dǎo)致葉片疲勞斷裂，轉(zhuǎn)子失衡及摩擦，并最終導(dǎo)致發(fā)動機停機發(fā)射失敗。為了減小空化不穩(wěn)定帶來的危害，各國學(xué)者做了多方面研究[9]。由于水力機械幾何過于復(fù)雜，便從較為簡單的水翼非定?？栈胧?。時素果等[10]研究了熱力學(xué)效應(yīng)對空化水動力脈動特性的影響，得到了水翼升力在非定常空化階段的特征頻率。Wang等[11]研究了附著空化流動，對空化初生進行了探究。Leroux等[12]對水翼的云狀空化進行研究，發(fā)現(xiàn)翼型表面不同位置壓強隨時間變化不同，得到云狀空化兩種不同的動態(tài)特性。尹必行等[13]采用試驗研究與數(shù)值模擬相結(jié)合的方法研究繞水翼ys930的非定常空化流場結(jié)構(gòu)，得到片狀空化和云狀空泡兩個不同的階段。
實驗研究依靠實驗設(shè)備、實驗人員和實驗場地等條件，使得實驗研究內(nèi)容的多樣性受到了限制。除了實驗，數(shù)值仿真也是研究空化的重要手段。目前有3種方法可以計算湍流，雖然DNS和LES能提供更好的計算精度，但是會消耗大量計算資源，目前還無法廣泛實現(xiàn)工程應(yīng)用，現(xiàn)在工程上最常用的還是RANS方法。Launder和Spalding[14]提出了標準k-ε雙方程模型，在RANS方法中預(yù)測結(jié)果更具有穩(wěn)定性和可靠性[15]，適用于絕大多數(shù)的工程湍流模型?？栈鲃邮且环N復(fù)雜的流動現(xiàn)象，包含湍流、相變、多相流、可壓縮和非定常性等等。標準k-ε模型在處理時間相關(guān)流動上的缺陷使得對非定?？栈哪M存在一定誤差[16]。有學(xué)者從標準k-ε模型的渦黏項入手，提出了一些修正方案。Johansen等[17]在標準k-ε模型的基礎(chǔ)上增加濾波器提出FBM模型，并計算了方柱繞流的非定常流動，結(jié)果表明FBM模型部分提高了對非定常流動的預(yù)測能力。Huang等[18]將FBM模型和DCM模型[19]中的濾波系數(shù)加權(quán)，得到一種新的濾波模型FBDCM模型并用該模型計算了Clark-Y型水翼的繞流空化，結(jié)果顯示FBDCM模型基于網(wǎng)格分辨率和當?shù)亓鲃拥拿芏刃拚梢源龠M空化流動的非定常脫落。洪鋒等[20]通過FBM模型和DCM模型中濾波項的比較選擇新的濾波函數(shù)，發(fā)展出MFBM模型，同樣仿真Clark-Y型水翼的繞流空化取得了較好的結(jié)果。以上3種修正的湍流模型與標準k-ε模型相比，都從不同程度上增強了預(yù)測非定常流動的能力，但這3種修正湍流模型之間目前還沒有一個對比評價。
綜上，本文分別選用原湍流模型和上述3種修正湍流模型，對Clark-Y型水翼表面的非定常空化流動進行模擬，通過與實驗結(jié)果的對比，對不同模型修正方式的計算效果做出評價。
1 數(shù)值計算方法
1.1 控制方程
本文的數(shù)值方法采用均相流模型，混合流體被認為是各向同性的，相間無速度滑移，氣相和液相有相同的速度和壓力，則混合流體的密度表達式為
ρ=αvρv+(1-αv)ρl
(1)
式中：ρ為混合相密度；ρv，ρl分別為氣相和液相密度；αv為氣相體積分數(shù)。
在笛卡爾坐標系下，連續(xù)方程和動量方程分別為
(2)
(3)
混合流體中的氣相輸運方程為
(4)
式中m+和m-分別為氣相中的蒸發(fā)源項和凝結(jié)源項，控制氣體的產(chǎn)生與潰滅。
1.2 空化模型
Zwart等[21]結(jié)合簡化的Rayleigh-Plesset方程，推導(dǎo)出的蒸發(fā)源項和凝結(jié)源項如下
(5)
(6)
式中：Ce和Cc分別為蒸發(fā)和凝結(jié)常數(shù)；RB為空泡半徑；αnuc為不凝結(jié)氣相分數(shù)。
該模型重點考慮了空化初生和發(fā)展時空泡體積變化的影響，適于模擬空化的非定常特性。
1.3 湍流模型
1.3.1 標準k-ε模型
標準k-ε模型由Launder和Spalding[14]提出，是在代數(shù)渦黏模型基礎(chǔ)上發(fā)展得到的渦黏模型，與代數(shù)模型的主要區(qū)別是將渦黏系數(shù)與湍動能k及湍動能耗散ε聯(lián)系在一起，包含部分湍流統(tǒng)計量之間的歷史效應(yīng)，有一定的物理依據(jù)，其控制方程為
(7)
(8)
湍流黏性系數(shù)μt表示為k和ε的函數(shù)，即
(9)
式中Pt為湍動能生成項。模型常數(shù)分別為：Cε1=1.44，Cε2=1.92，σε=1.3，σk=1.0，Cμ=0.09。
1.3.2 FBM模型
RANS方法的計算精度除了與網(wǎng)格尺寸有關(guān)外，還與流場的渦黏系數(shù)即雷諾應(yīng)力有關(guān)。Johansen等人[13]認為標準k-ε模型對湍流預(yù)測能力不足是由于流場中渦黏性過大，為了解決這個問題提出FBM模型。FBM模型在原模型湍流黏性系數(shù)的基礎(chǔ)上增加濾波器，得到新的湍流黏性系數(shù)為
(10)
式中fFBM為濾波函數(shù)，由濾波器尺寸Δ和湍流結(jié)構(gòu)尺寸的比值決定，定義為
(11)
為了確保濾波與數(shù)值方法相容，所選的濾波器尺寸應(yīng)不小于計算網(wǎng)格區(qū)域的大小，即Δ>Δgrid，這里的網(wǎng)格大小定義為Δgrid=(Δx,Δy,Δz)1/3，Δx，Δy和Δz分別為網(wǎng)格在3個坐標方向的長度。
對于湍流尺度小于濾波器尺寸的流動，渦黏保持標準k-ε模型的黏度不變；對于湍流尺度大于濾波器尺寸的流動，由式(11)可知，其湍流黏度變?yōu)?/p>
(12)
1.3.3 FBDCM模型
Huang等[18]認為FBM模型沒有正確修正近壁的渦黏，在空穴周圍沒有相變的區(qū)域，F(xiàn)BM模型計算的渦黏仍然過大?？紤]到水翼繞流的動態(tài)特性，可以將空化區(qū)分為附著空穴和分離空穴兩個子域。在附著空穴區(qū)，蒸汽的相對體積分數(shù)較高，可以看成可壓的兩相流。在分離的空穴區(qū)內(nèi)是大尺度的氣液混合不穩(wěn)定運動，需要將當?shù)氐耐牧鹘Y(jié)構(gòu)尺寸考慮在內(nèi)。由此，Huang等在FBM模型的基礎(chǔ)上，結(jié)合Coutier-Delgosha等[19]提出的用流場的當?shù)孛芏却嫫骄芏鹊腄CM模型，提出了FBDCM模型。與FBM模型類似，只對渦黏項進行了修正
(13)
(14)
(15)
式中χρ/ρl用來加權(quán)兩種模型的濾波函數(shù)，定義為
(16)
式中：C1=4；C2=0.2。
1.3.4 MFBM模型
洪鋒等[20]認為，將FBM模型和DCM模型的濾波函數(shù)相互比較后得到的MFBM模型可以兼具兩種模型的優(yōu)點，則MFBM的湍流黏度系數(shù)為
(17)
fMFBM=min(fFBM,fDCM)
(18)
1.4 數(shù)值計算
本文以二維Clark-Y型水翼作為研究對象，翼型的弦長c=70 mm，攻角α=8°，空化數(shù)σ=0.8，計算域的幾何參數(shù)如圖1所示。邊界條件設(shè)置采用速度入口條件，流速V∞=7.8 m/s，湍流強度I=2%；出口為壓力條件，利用空化數(shù)可以推導(dǎo)出口的壓力值；上下壁面及水翼壁面均設(shè)置為無滑移壁面。計算介質(zhì)采用25 ℃的水和水蒸氣，飽和蒸汽壓為3 574 Pa。本文的實驗數(shù)據(jù)來自文獻[10]。數(shù)值計算平臺為ANSYS CFX 15.0。
圖1 計算域幾何示意圖Fig.1 Geometric diagram of computational fluid domain
首先利用標準k-ε模型計算繞二維水翼的定?？栈鲃觼眚炞C網(wǎng)格無關(guān)性。本文共驗證了5套網(wǎng)格，網(wǎng)格數(shù)分別為77 056(mesh1)、144 770(mesh2)、187 524(mesh3)、213 512(mesh4)和275 667(mesh5)，水翼附近采用O型網(wǎng)格劃分，并對靠近壁面的網(wǎng)格進行加密，結(jié)果顯示30利用這5套網(wǎng)格計算水翼在定?？栈癄顟B(tài)下的升力系數(shù)
(19)
式中：Clift為水翼的升力系數(shù);Flift為水翼表面的升力;V∞為無窮遠處的速度;S為參考面積。
在不同網(wǎng)格下計算的升力系數(shù)與實驗數(shù)據(jù)對比如圖2所示。
從圖2可以看出，mesh2的誤差是最小的，但只有在網(wǎng)格數(shù)量大于187 523(mesh3)后，計算結(jié)果不隨網(wǎng)格數(shù)改變而改變，說明mesh2的小誤差只是偶然的，不能作為最終的計算網(wǎng)格。mesh3計算的升力系數(shù)與實驗數(shù)據(jù)的誤差小于5%，滿足了網(wǎng)格無關(guān)性的要求。為了節(jié)省時間和計算資源，本文選用mesh3作為最終的計算網(wǎng)格，如圖3所示。
圖2 仿真計算的升力系數(shù)與實驗數(shù)據(jù)的對比Fig.2 Comparison of calculated lift coefficient with experiment data
圖3 計算網(wǎng)格Fig.3 Computational grids
2 結(jié)果與討論
2.1 仿真結(jié)果比較
2.1.1 翼型升力系數(shù)的對比
圖4是4種不同湍流模型計算水翼空化得到的基于升力系數(shù)的FFT變換。從4張圖中均能看出比較突出的主頻，表明這里的4種湍流模型都能夠捕捉該工況下的非定常流動特性。但標準k-ε模型的FFT變化與其他3種模型明顯不同，除了能夠分辨主頻，次頻也很明顯，主頻與次頻之間不存在其他雜頻。對于其他3種修正模型，雖然主頻也很明顯，但幾乎無法分辨次頻，主頻附近分布許多小雜頻，表明修正的湍流模型顯著增強了計算結(jié)果的非定常性。
將4種模型的計算結(jié)果與實驗數(shù)據(jù)進行對比，如表1所示，其中基于弦長的斯特勞哈爾數(shù)
(20)
式中f為升力系數(shù)變化的頻率。
圖4 關(guān)于升力系數(shù)的快速傅里葉變換Fig.4 Fast fourier transform of lift coefficient around hydrofoil
表1 不同湍流模型計算的St與實驗數(shù)據(jù)的對比
從表1可以看出，只有FBM模型的St要大于實驗值，即升力變化周期要??；而標準k-ε模型、FBDCM模型和MFBM模型的St要小于實驗值，即升力變化周期要長一些。4種湍流模型中標準k-ε模型的St是最小的，其他3種模型的St都大于0.151，說明修正模型的非定常性更加突出。在3種修正模型中，F(xiàn)BM模型的St過大，誤差高達36.9%，效果反而不如標準k-ε模型；FBDCM模型和MFBM模型的St與實驗數(shù)據(jù)較為吻合。MFBM模型的St是最接近實驗值的，與實驗數(shù)據(jù)的誤差不到1%，在4種湍流模型中能夠更好地捕捉流動的非定常特性。
2.1.2 空穴形態(tài)的對比
圖5給出了4種湍流模型計算的空穴形態(tài)與實驗數(shù)據(jù)的對比。從圖5中可以看出繞水翼空化由兩部分組成，頭部的附著空化和尾部的脫落空化。標準k-ε模型的空穴形態(tài)在長度和厚度以及脫落空化的面積上都與實驗存在較大的誤差，即使在空穴發(fā)展最充分的階段(t0+28 ms)，標準k-ε模型的部分空化長度剛剛超過弦長的一半，而實驗數(shù)據(jù)則顯示空穴覆蓋了整個吸力面，反映出標準k-ε模型在預(yù)測非定常空化流動時存在較大的缺陷。3種修正湍流模型采用相同的濾波尺寸進行計算，在空穴形態(tài)的預(yù)測上有較大的改進，模擬出的空化附著部分長度、厚度增加，脫落部分的空化也更為明顯。修正模型中，F(xiàn)BM模型的改進效果是最不明顯的，在空穴區(qū)變化周期的中間階段(t0+21 ms～t0+42 ms)，附著空化最長可以發(fā)展到水翼吸力面3/4的位置，但沒有產(chǎn)生超空化；FBDCM模型的空穴在3種修正湍流模型最長最厚，特別是在t0+28 ms和t0+35 ms兩個階段，比實驗數(shù)據(jù)顯示的空化區(qū)域還要大，是FBDCM模型的不足之處；MFBM模型仿真的結(jié)果是最為接近實驗數(shù)據(jù)的一組結(jié)果，空穴區(qū)的發(fā)展形態(tài)以及脫落后與實驗數(shù)據(jù)相比都能一一對應(yīng)。FBM模型以濾波尺寸為界限，沒有修正小于濾波尺寸網(wǎng)格內(nèi)的渦黏，而FBDCM模型和MFBM模型在此基礎(chǔ)上從不同角度修正了翼型近壁區(qū)的流場，是對全流場的修正，因而后兩種模型能更好地捕捉空化的非定常流動特性。從計算結(jié)果來看，MFBM模型的效果是最佳的。通過模型的控制方程分析，MFBM模型中，小于濾波器尺寸的近壁區(qū)網(wǎng)格用DCM模型計算，而近壁區(qū)的流動是非定常性較強的區(qū)域，伴隨著相變和密度的變化，DCM模型中用當?shù)孛芏却嫫骄芏冗M行計算是較為合理的。
圖5 不同湍流模型計算的氣相分數(shù)云圖與實驗數(shù)據(jù)的對比Fig.5 Comparison of vapor volume fraction contours calculated by four turbulent models with experimental data
2.1.3 流場渦黏系數(shù)的對比
渦黏系數(shù)跟流動非定常性直接相關(guān)，標準k-ε模型的渦黏系數(shù)與湍動能k與湍動能耗散ε直接相關(guān)，無法準確?；膭幽芘c湍動能耗散導(dǎo)致標準k-ε模型預(yù)測非定常流動效果不佳。為了分析不同湍流模型對非定常流動預(yù)測能力差異形成的原因，圖6給出了4種湍流模型計算的流場渦黏。從圖中可以看出標準k-ε模型水翼吸力面后半部分以及尾跡區(qū)的渦黏都非常大，這勢必會影響模型預(yù)測非定常流動的準確性，使得仿真結(jié)果出現(xiàn)較大的誤差。3種修正模型則大幅減小了流場的渦黏系數(shù)，這是修正模型能夠更好地預(yù)測非定?？栈鲃拥闹匾颉Ｓ捎诹鲃拥姆嵌ǔＬ匦缘靡酝癸@，用修正模型計算得到的結(jié)果會在真值附近出現(xiàn)較大的波動，所以在圖4的FFT變換中，修正模型的頻域圖中除了突出的主頻還伴隨著許多小雜頻。3種修正模型中，F(xiàn)BM模型計算的渦黏是最大的，因為沒有修正近壁處的流場；FBDCM模型和MFBM模型計算的渦黏非常相近且略小于FBM模型，實際上這兩種模型計算得到的基于升力變化的St也非常相似，與渦黏云圖的結(jié)果是一致的。
MFBM模型能夠更好地預(yù)測非定常空化流動特性，標準k-ε模型是在代數(shù)渦黏模型的基礎(chǔ)上發(fā)展起來的，它和代數(shù)渦黏模型的主要區(qū)別在于標準k-ε模型的渦黏系數(shù)包含部分歷史效應(yīng)[16]，MFBM模型作為標準k-ε模型的延伸，修正了渦黏過大對于流場仿真結(jié)果的影響，但并沒有改變湍動能和湍動能耗散的生成。因而，MFBM模型雖然可以更好地預(yù)測空化流動的非定常特性，但MFBM模型仍然無法準確預(yù)測湍流統(tǒng)計量關(guān)系之間的歷史效應(yīng)，導(dǎo)致仿真結(jié)果與實驗存在一定誤差。雖然MFBM模型結(jié)合了FBM模型和DCM模型來修正近壁區(qū)流場，但只是在計算中選擇使當?shù)販u黏更小的濾波系數(shù)進行模擬，物理依據(jù)上有所欠缺。如果能在FBM模型的基礎(chǔ)上更科學(xué)地修正小于濾波尺寸網(wǎng)格的區(qū)域，還可以更加準確地捕捉非定常流動的特性。
圖6 不同湍流模型計算的渦黏Fig.6 Eddy viscosity calculated by four turbulent models
2.2 水翼非定?？栈鲃臃治?/h3>
2.2.1 逆壓梯度的對比
圖7是一個周期內(nèi)水翼表面壓力系數(shù)與氣相體積分數(shù)分布圖。圖7(a)中水翼吸力面的低壓區(qū)的長度和位置與圖7(b)中空穴的長度和位置一一對應(yīng)，頭部的低壓區(qū)對應(yīng)附著空化，尾部的低壓區(qū)對應(yīng)脫落的空化?？栈跎?t0+7 ms)時，在吸力面的后半部分有一段長度為0.3c的空穴，是還沒有完全潰滅的脫落空化，在0到0.3c的區(qū)間存在逆壓梯度。隨著時間的發(fā)展，t0+14 ms時，尾部的脫落空化已經(jīng)消失，初生的空穴區(qū)后0.2c到0.7c的區(qū)間內(nèi)存在較強的逆壓梯度，說明伴隨著脫落空化的潰滅逆壓梯度被增強。t0+21 ms時，脫落空穴的影響基本消失，逆壓梯度強度減小，僅在吸力面0.75c的位置存在一個非常小的逆壓梯度，附著空化得以進一步發(fā)展，空穴長度達到0.7c。此時部分空化脫落較少，逆壓梯度對空穴生長的限制作用不明顯，吸力面空化會進一步生長成為超空化，從圖7(a)圖和圖7(b)中可以看出t0+28 ms和t0+35 ms時，翼型吸力面被空化覆蓋，表面也不存在壓力梯度。超空化后，空化開始大量脫落，不斷潰滅形成的局部高壓區(qū)再次加強逆壓梯度，逆壓梯度又會阻礙空穴生長，促進空化脫落，直到水翼表面的空穴消失，從t0+42 ms到t0+56 ms的3個時間段，逆壓梯度位置不斷向上游移動，強度也不斷增加，直至空穴完全脫落消失。以上的變化說明在一個周期內(nèi)，空穴的生長和脫落受逆壓梯度強度和位置的影響，而空穴的生長與脫落又影響著水翼表面的壓力分布。因而，翼型升力變化的周期與空穴的脫落周期是對應(yīng)的。
圖7 水翼表面升力系數(shù)與氣相分數(shù)分布Fig.7 Distribution of lift coefficient and vapor volume fraction around hydrofoil
圖8給出了空穴長度隨時間的變化。可見水翼表面的空化存在周期性脫落，脫落的周期為40 ms～50 ms，與升力變化周期56 ms[6]基本吻合，表明空化云周期性脫落是升力系數(shù)周期性變化的原因，這與前文中分析的結(jié)果是一致的。
圖8 空穴長度隨時間的變化Fig.8 Change of cavity length with time
2.2.2 反向射流的對比
圖9給出了一個周期內(nèi)翼型表面的速度矢量分布。可以看出翼型吸力面頭部速度最大，尾部速度最小并伴隨回流區(qū)。根據(jù)伯努利原理，水翼吸力頭部速度較大。隨著流體沿著翼型向下游運動，逆壓梯度會使流動出現(xiàn)分離，形成反向射流。t0+7 ms時，上一個周期的空化還未完全脫落，壓力分布上翼型表面存在較大的壓力梯度，速度分布上從圖9中可以看到在吸力面0.3c到尾部之間存在較強的回流區(qū)。t0+14 ms時，回流區(qū)移動到翼型尾部，即將脫離水翼，新周期的空化已經(jīng)初生。t0+21 ms時，隨著逆壓梯度的消失，回流區(qū)大部分消失在主流中，附著在翼型表面的回流區(qū)面積大大減小。t0+28 ms和t0+35 ms時，逆壓梯度作用不明顯，回流區(qū)面積基本沒有增加，使得部分空化發(fā)展成為超空化。在此之后，t0+42 ms時空穴充分發(fā)展，尾部的空化與回流區(qū)開始相互干涉?；亓髟鰪娏丝栈拿撀?，脫落空化潰滅瞬間的局部高壓造成的逆壓梯度又加強了回流，使得回流區(qū)的長度和厚度開始增加。另一方面，空化潰滅后反向射流會迅速填滿空泡破裂后的間隙，使得回流沿翼型向上游運動，進一步擴大回流區(qū)的范圍。t0+49 ms時，隨著回流厚度不斷增加，反向射流已經(jīng)運動到水翼的頭部，與空化區(qū)氣液交界面相互作用，造成空化云的大尺度脫落，完成空穴在一個周期內(nèi)的變化。對比圖7，反向射流的位置與逆壓梯度出現(xiàn)的位置是一一對應(yīng)的，說明反向射流和逆壓梯度共同影響空穴的生長與脫落。
圖9 水翼表面速度矢量分布Fig.9 Velocity vectors distribution around hydrofoil
3 結(jié)論
本文采用標準k-ε模型和3種修正湍流模型對二維Clark-Y水翼空化流動進行計算，通過仿真結(jié)果的相互比較及與實驗結(jié)果的對比，得出以下結(jié)論：
1)與標準k-ε模型相比，修正的湍流模型可以大幅減小流場中的渦黏，使得雷諾應(yīng)力的部分歷史效應(yīng)得以凸顯，從而增強模型對于非定?？栈鲃拥念A(yù)測能力。
2)MFBM模型在3種修正模型中，捕捉的非定?？栈鲃犹匦耘c實驗最為接近，效果最好，能夠比較有效地模擬二維水翼非定常空化流動。
猜你喜歡
空穴效應(yīng)下泡沫金屬復(fù)合相變材料熱性能數(shù)值模擬應(yīng)用數(shù)學(xué)和力學(xué)(2024年1期)2024-01-29 08:32:46 功率超聲作用下鋼液中空化泡尺寸的演變特性鋼鐵釩鈦(2023年5期)2023-11-17 08:48:34 噴油嘴內(nèi)部空穴流動試驗研究車用發(fā)動機(2021年5期)2021-10-31 05:48:38 波浪滑翔機橢圓形后緣水翼動力特性研究海洋技術(shù)學(xué)報(2021年3期)2021-08-19 11:07:04 基于MoOx選擇性接觸的SHJ太陽電池研究進展電源技術(shù)(2021年7期)2021-07-29 08:35:24 袖珍水翼突防潛艇的設(shè)計構(gòu)想及運用研究數(shù)字海洋與水下攻防(2020年5期)2021-01-04 01:40:20 三維扭曲水翼空化現(xiàn)象CFD模擬廈門理工學(xué)院學(xué)報(2016年1期)2016-12-01 04:50:53 不同運動形式下水物相互作用空化數(shù)值模擬水利科技與經(jīng)濟(2016年7期)2016-04-25 13:03:00 湍流進流誘發(fā)的二維水翼振動噪聲特性研究噪聲與振動控制(2015年4期)2015-01-01 07:08:04 負壓空化法提取銀杏葉中原花青素及其分離純化食品科學(xué)(2013年24期)2013-03-11 18:30:24
李雨濛，陳 暉，項 樂，張亞太
(西安航天動力研究所 液體火箭發(fā)動機重點實驗室，陜西 西安 710100)
0 引言
近年來航天不斷發(fā)展，液體火箭發(fā)動機作為重要的動力來源也不斷更新?lián)Q代。新一代液體火箭發(fā)動機以液氫液氧或液氧煤油作為推進劑，比沖高，推力大，工作時間長，在大型運載火箭上得到廣泛應(yīng)用。渦輪泵作為液體火箭發(fā)動機的核心部件，在高轉(zhuǎn)速的環(huán)境下工作易導(dǎo)致渦輪泵內(nèi)出現(xiàn)空化，影響液體火箭發(fā)動機的可靠性[1-5]。
空化作為水力機械中經(jīng)常出現(xiàn)的一種現(xiàn)象，會造成水力機械性能明顯降低[6]，材料表面破壞，引起振動和噪聲等[7-8]?？栈瘜討B(tài)響應(yīng)特性的改變會使流動內(nèi)部出現(xiàn)不穩(wěn)定性，多個國家在研究液體火箭發(fā)動機渦輪泵時都碰到過空化不穩(wěn)定帶來的問題乃至事故。例如，在日本H-2火箭的第8次發(fā)射中，空化不穩(wěn)定誘發(fā)的脈動頻率與泵前導(dǎo)流葉片固有頻率相近，引起共振導(dǎo)致葉片疲勞斷裂，轉(zhuǎn)子失衡及摩擦，并最終導(dǎo)致發(fā)動機停機發(fā)射失敗。為了減小空化不穩(wěn)定帶來的危害，各國學(xué)者做了多方面研究[9]。由于水力機械幾何過于復(fù)雜，便從較為簡單的水翼非定?？栈胧?。時素果等[10]研究了熱力學(xué)效應(yīng)對空化水動力脈動特性的影響，得到了水翼升力在非定常空化階段的特征頻率。Wang等[11]研究了附著空化流動，對空化初生進行了探究。Leroux等[12]對水翼的云狀空化進行研究，發(fā)現(xiàn)翼型表面不同位置壓強隨時間變化不同，得到云狀空化兩種不同的動態(tài)特性。尹必行等[13]采用試驗研究與數(shù)值模擬相結(jié)合的方法研究繞水翼ys930的非定常空化流場結(jié)構(gòu)，得到片狀空化和云狀空泡兩個不同的階段。
實驗研究依靠實驗設(shè)備、實驗人員和實驗場地等條件，使得實驗研究內(nèi)容的多樣性受到了限制。除了實驗，數(shù)值仿真也是研究空化的重要手段。目前有3種方法可以計算湍流，雖然DNS和LES能提供更好的計算精度，但是會消耗大量計算資源，目前還無法廣泛實現(xiàn)工程應(yīng)用，現(xiàn)在工程上最常用的還是RANS方法。Launder和Spalding[14]提出了標準k-ε雙方程模型，在RANS方法中預(yù)測結(jié)果更具有穩(wěn)定性和可靠性[15]，適用于絕大多數(shù)的工程湍流模型?？栈鲃邮且环N復(fù)雜的流動現(xiàn)象，包含湍流、相變、多相流、可壓縮和非定常性等等。標準k-ε模型在處理時間相關(guān)流動上的缺陷使得對非定?？栈哪M存在一定誤差[16]。有學(xué)者從標準k-ε模型的渦黏項入手，提出了一些修正方案。Johansen等[17]在標準k-ε模型的基礎(chǔ)上增加濾波器提出FBM模型，并計算了方柱繞流的非定常流動，結(jié)果表明FBM模型部分提高了對非定常流動的預(yù)測能力。Huang等[18]將FBM模型和DCM模型[19]中的濾波系數(shù)加權(quán)，得到一種新的濾波模型FBDCM模型并用該模型計算了Clark-Y型水翼的繞流空化，結(jié)果顯示FBDCM模型基于網(wǎng)格分辨率和當?shù)亓鲃拥拿芏刃拚梢源龠M空化流動的非定常脫落。洪鋒等[20]通過FBM模型和DCM模型中濾波項的比較選擇新的濾波函數(shù)，發(fā)展出MFBM模型，同樣仿真Clark-Y型水翼的繞流空化取得了較好的結(jié)果。以上3種修正的湍流模型與標準k-ε模型相比，都從不同程度上增強了預(yù)測非定常流動的能力，但這3種修正湍流模型之間目前還沒有一個對比評價。
綜上，本文分別選用原湍流模型和上述3種修正湍流模型，對Clark-Y型水翼表面的非定常空化流動進行模擬，通過與實驗結(jié)果的對比，對不同模型修正方式的計算效果做出評價。
1 數(shù)值計算方法
1.1 控制方程
本文的數(shù)值方法采用均相流模型，混合流體被認為是各向同性的，相間無速度滑移，氣相和液相有相同的速度和壓力，則混合流體的密度表達式為
ρ=αvρv+(1-αv)ρl
(1)
式中：ρ為混合相密度；ρv，ρl分別為氣相和液相密度；αv為氣相體積分數(shù)。
在笛卡爾坐標系下，連續(xù)方程和動量方程分別為
(2)
(3)
混合流體中的氣相輸運方程為
(4)
式中m+和m-分別為氣相中的蒸發(fā)源項和凝結(jié)源項，控制氣體的產(chǎn)生與潰滅。
1.2 空化模型
Zwart等[21]結(jié)合簡化的Rayleigh-Plesset方程，推導(dǎo)出的蒸發(fā)源項和凝結(jié)源項如下
(5)
(6)
式中：Ce和Cc分別為蒸發(fā)和凝結(jié)常數(shù)；RB為空泡半徑；αnuc為不凝結(jié)氣相分數(shù)。
該模型重點考慮了空化初生和發(fā)展時空泡體積變化的影響，適于模擬空化的非定常特性。
1.3 湍流模型
1.3.1 標準k-ε模型
標準k-ε模型由Launder和Spalding[14]提出，是在代數(shù)渦黏模型基礎(chǔ)上發(fā)展得到的渦黏模型，與代數(shù)模型的主要區(qū)別是將渦黏系數(shù)與湍動能k及湍動能耗散ε聯(lián)系在一起，包含部分湍流統(tǒng)計量之間的歷史效應(yīng)，有一定的物理依據(jù)，其控制方程為
(7)
(8)
湍流黏性系數(shù)μt表示為k和ε的函數(shù)，即
(9)
式中Pt為湍動能生成項。模型常數(shù)分別為：Cε1=1.44，Cε2=1.92，σε=1.3，σk=1.0，Cμ=0.09。
1.3.2 FBM模型
RANS方法的計算精度除了與網(wǎng)格尺寸有關(guān)外，還與流場的渦黏系數(shù)即雷諾應(yīng)力有關(guān)。Johansen等人[13]認為標準k-ε模型對湍流預(yù)測能力不足是由于流場中渦黏性過大，為了解決這個問題提出FBM模型。FBM模型在原模型湍流黏性系數(shù)的基礎(chǔ)上增加濾波器，得到新的湍流黏性系數(shù)為
(10)
式中fFBM為濾波函數(shù)，由濾波器尺寸Δ和湍流結(jié)構(gòu)尺寸的比值決定，定義為
(11)
為了確保濾波與數(shù)值方法相容，所選的濾波器尺寸應(yīng)不小于計算網(wǎng)格區(qū)域的大小，即Δ>Δgrid，這里的網(wǎng)格大小定義為Δgrid=(Δx,Δy,Δz)1/3，Δx，Δy和Δz分別為網(wǎng)格在3個坐標方向的長度。
對于湍流尺度小于濾波器尺寸的流動，渦黏保持標準k-ε模型的黏度不變；對于湍流尺度大于濾波器尺寸的流動，由式(11)可知，其湍流黏度變?yōu)?/p>
(12)
1.3.3 FBDCM模型
Huang等[18]認為FBM模型沒有正確修正近壁的渦黏，在空穴周圍沒有相變的區(qū)域，F(xiàn)BM模型計算的渦黏仍然過大?？紤]到水翼繞流的動態(tài)特性，可以將空化區(qū)分為附著空穴和分離空穴兩個子域。在附著空穴區(qū)，蒸汽的相對體積分數(shù)較高，可以看成可壓的兩相流。在分離的空穴區(qū)內(nèi)是大尺度的氣液混合不穩(wěn)定運動，需要將當?shù)氐耐牧鹘Y(jié)構(gòu)尺寸考慮在內(nèi)。由此，Huang等在FBM模型的基礎(chǔ)上，結(jié)合Coutier-Delgosha等[19]提出的用流場的當?shù)孛芏却嫫骄芏鹊腄CM模型，提出了FBDCM模型。與FBM模型類似，只對渦黏項進行了修正
(13)
(14)
(15)
式中χρ/ρl用來加權(quán)兩種模型的濾波函數(shù)，定義為
(16)
式中：C1=4；C2=0.2。
1.3.4 MFBM模型
洪鋒等[20]認為，將FBM模型和DCM模型的濾波函數(shù)相互比較后得到的MFBM模型可以兼具兩種模型的優(yōu)點，則MFBM的湍流黏度系數(shù)為
(17)
fMFBM=min(fFBM,fDCM)
(18)
1.4 數(shù)值計算
本文以二維Clark-Y型水翼作為研究對象，翼型的弦長c=70 mm，攻角α=8°，空化數(shù)σ=0.8，計算域的幾何參數(shù)如圖1所示。邊界條件設(shè)置采用速度入口條件，流速V∞=7.8 m/s，湍流強度I=2%；出口為壓力條件，利用空化數(shù)可以推導(dǎo)出口的壓力值；上下壁面及水翼壁面均設(shè)置為無滑移壁面。計算介質(zhì)采用25 ℃的水和水蒸氣，飽和蒸汽壓為3 574 Pa。本文的實驗數(shù)據(jù)來自文獻[10]。數(shù)值計算平臺為ANSYS CFX 15.0。
圖1 計算域幾何示意圖Fig.1 Geometric diagram of computational fluid domain
首先利用標準k-ε模型計算繞二維水翼的定?？栈鲃觼眚炞C網(wǎng)格無關(guān)性。本文共驗證了5套網(wǎng)格，網(wǎng)格數(shù)分別為77 056(mesh1)、144 770(mesh2)、187 524(mesh3)、213 512(mesh4)和275 667(mesh5)，水翼附近采用O型網(wǎng)格劃分，并對靠近壁面的網(wǎng)格進行加密，結(jié)果顯示30 利用這5套網(wǎng)格計算水翼在定?？栈癄顟B(tài)下的升力系數(shù) (19) 式中：Clift為水翼的升力系數(shù);Flift為水翼表面的升力;V∞為無窮遠處的速度;S為參考面積。 在不同網(wǎng)格下計算的升力系數(shù)與實驗數(shù)據(jù)對比如圖2所示。 從圖2可以看出，mesh2的誤差是最小的，但只有在網(wǎng)格數(shù)量大于187 523(mesh3)后，計算結(jié)果不隨網(wǎng)格數(shù)改變而改變，說明mesh2的小誤差只是偶然的，不能作為最終的計算網(wǎng)格。mesh3計算的升力系數(shù)與實驗數(shù)據(jù)的誤差小于5%，滿足了網(wǎng)格無關(guān)性的要求。為了節(jié)省時間和計算資源，本文選用mesh3作為最終的計算網(wǎng)格，如圖3所示。 圖2 仿真計算的升力系數(shù)與實驗數(shù)據(jù)的對比Fig.2 Comparison of calculated lift coefficient with experiment data 圖3 計算網(wǎng)格Fig.3 Computational grids 2.1.1 翼型升力系數(shù)的對比 圖4是4種不同湍流模型計算水翼空化得到的基于升力系數(shù)的FFT變換。從4張圖中均能看出比較突出的主頻，表明這里的4種湍流模型都能夠捕捉該工況下的非定常流動特性。但標準k-ε模型的FFT變化與其他3種模型明顯不同，除了能夠分辨主頻，次頻也很明顯，主頻與次頻之間不存在其他雜頻。對于其他3種修正模型，雖然主頻也很明顯，但幾乎無法分辨次頻，主頻附近分布許多小雜頻，表明修正的湍流模型顯著增強了計算結(jié)果的非定常性。 將4種模型的計算結(jié)果與實驗數(shù)據(jù)進行對比，如表1所示，其中基于弦長的斯特勞哈爾數(shù) (20) 式中f為升力系數(shù)變化的頻率。 圖4 關(guān)于升力系數(shù)的快速傅里葉變換Fig.4 Fast fourier transform of lift coefficient around hydrofoil 表1 不同湍流模型計算的St與實驗數(shù)據(jù)的對比 從表1可以看出，只有FBM模型的St要大于實驗值，即升力變化周期要??；而標準k-ε模型、FBDCM模型和MFBM模型的St要小于實驗值，即升力變化周期要長一些。4種湍流模型中標準k-ε模型的St是最小的，其他3種模型的St都大于0.151，說明修正模型的非定常性更加突出。在3種修正模型中，F(xiàn)BM模型的St過大，誤差高達36.9%，效果反而不如標準k-ε模型；FBDCM模型和MFBM模型的St與實驗數(shù)據(jù)較為吻合。MFBM模型的St是最接近實驗值的，與實驗數(shù)據(jù)的誤差不到1%，在4種湍流模型中能夠更好地捕捉流動的非定常特性。 2.1.2 空穴形態(tài)的對比 圖5給出了4種湍流模型計算的空穴形態(tài)與實驗數(shù)據(jù)的對比。從圖5中可以看出繞水翼空化由兩部分組成，頭部的附著空化和尾部的脫落空化。標準k-ε模型的空穴形態(tài)在長度和厚度以及脫落空化的面積上都與實驗存在較大的誤差，即使在空穴發(fā)展最充分的階段(t0+28 ms)，標準k-ε模型的部分空化長度剛剛超過弦長的一半，而實驗數(shù)據(jù)則顯示空穴覆蓋了整個吸力面，反映出標準k-ε模型在預(yù)測非定常空化流動時存在較大的缺陷。3種修正湍流模型采用相同的濾波尺寸進行計算，在空穴形態(tài)的預(yù)測上有較大的改進，模擬出的空化附著部分長度、厚度增加，脫落部分的空化也更為明顯。修正模型中，F(xiàn)BM模型的改進效果是最不明顯的，在空穴區(qū)變化周期的中間階段(t0+21 ms～t0+42 ms)，附著空化最長可以發(fā)展到水翼吸力面3/4的位置，但沒有產(chǎn)生超空化；FBDCM模型的空穴在3種修正湍流模型最長最厚，特別是在t0+28 ms和t0+35 ms兩個階段，比實驗數(shù)據(jù)顯示的空化區(qū)域還要大，是FBDCM模型的不足之處；MFBM模型仿真的結(jié)果是最為接近實驗數(shù)據(jù)的一組結(jié)果，空穴區(qū)的發(fā)展形態(tài)以及脫落后與實驗數(shù)據(jù)相比都能一一對應(yīng)。FBM模型以濾波尺寸為界限，沒有修正小于濾波尺寸網(wǎng)格內(nèi)的渦黏，而FBDCM模型和MFBM模型在此基礎(chǔ)上從不同角度修正了翼型近壁區(qū)的流場，是對全流場的修正，因而后兩種模型能更好地捕捉空化的非定常流動特性。從計算結(jié)果來看，MFBM模型的效果是最佳的。通過模型的控制方程分析，MFBM模型中，小于濾波器尺寸的近壁區(qū)網(wǎng)格用DCM模型計算，而近壁區(qū)的流動是非定常性較強的區(qū)域，伴隨著相變和密度的變化，DCM模型中用當?shù)孛芏却嫫骄芏冗M行計算是較為合理的。 圖5 不同湍流模型計算的氣相分數(shù)云圖與實驗數(shù)據(jù)的對比Fig.5 Comparison of vapor volume fraction contours calculated by four turbulent models with experimental data 2.1.3 流場渦黏系數(shù)的對比 渦黏系數(shù)跟流動非定常性直接相關(guān)，標準k-ε模型的渦黏系數(shù)與湍動能k與湍動能耗散ε直接相關(guān)，無法準確?；膭幽芘c湍動能耗散導(dǎo)致標準k-ε模型預(yù)測非定常流動效果不佳。為了分析不同湍流模型對非定常流動預(yù)測能力差異形成的原因，圖6給出了4種湍流模型計算的流場渦黏。從圖中可以看出標準k-ε模型水翼吸力面后半部分以及尾跡區(qū)的渦黏都非常大，這勢必會影響模型預(yù)測非定常流動的準確性，使得仿真結(jié)果出現(xiàn)較大的誤差。3種修正模型則大幅減小了流場的渦黏系數(shù)，這是修正模型能夠更好地預(yù)測非定?？栈鲃拥闹匾颉Ｓ捎诹鲃拥姆嵌ǔＬ匦缘靡酝癸@，用修正模型計算得到的結(jié)果會在真值附近出現(xiàn)較大的波動，所以在圖4的FFT變換中，修正模型的頻域圖中除了突出的主頻還伴隨著許多小雜頻。3種修正模型中，F(xiàn)BM模型計算的渦黏是最大的，因為沒有修正近壁處的流場；FBDCM模型和MFBM模型計算的渦黏非常相近且略小于FBM模型，實際上這兩種模型計算得到的基于升力變化的St也非常相似，與渦黏云圖的結(jié)果是一致的。 MFBM模型能夠更好地預(yù)測非定常空化流動特性，標準k-ε模型是在代數(shù)渦黏模型的基礎(chǔ)上發(fā)展起來的，它和代數(shù)渦黏模型的主要區(qū)別在于標準k-ε模型的渦黏系數(shù)包含部分歷史效應(yīng)[16]，MFBM模型作為標準k-ε模型的延伸，修正了渦黏過大對于流場仿真結(jié)果的影響，但并沒有改變湍動能和湍動能耗散的生成。因而，MFBM模型雖然可以更好地預(yù)測空化流動的非定常特性，但MFBM模型仍然無法準確預(yù)測湍流統(tǒng)計量關(guān)系之間的歷史效應(yīng)，導(dǎo)致仿真結(jié)果與實驗存在一定誤差。雖然MFBM模型結(jié)合了FBM模型和DCM模型來修正近壁區(qū)流場，但只是在計算中選擇使當?shù)販u黏更小的濾波系數(shù)進行模擬，物理依據(jù)上有所欠缺。如果能在FBM模型的基礎(chǔ)上更科學(xué)地修正小于濾波尺寸網(wǎng)格的區(qū)域，還可以更加準確地捕捉非定常流動的特性。 圖6 不同湍流模型計算的渦黏Fig.6 Eddy viscosity calculated by four turbulent models 2.2.1 逆壓梯度的對比 圖7是一個周期內(nèi)水翼表面壓力系數(shù)與氣相體積分數(shù)分布圖。圖7(a)中水翼吸力面的低壓區(qū)的長度和位置與圖7(b)中空穴的長度和位置一一對應(yīng)，頭部的低壓區(qū)對應(yīng)附著空化，尾部的低壓區(qū)對應(yīng)脫落的空化?？栈跎?t0+7 ms)時，在吸力面的后半部分有一段長度為0.3c的空穴，是還沒有完全潰滅的脫落空化，在0到0.3c的區(qū)間存在逆壓梯度。隨著時間的發(fā)展，t0+14 ms時，尾部的脫落空化已經(jīng)消失，初生的空穴區(qū)后0.2c到0.7c的區(qū)間內(nèi)存在較強的逆壓梯度，說明伴隨著脫落空化的潰滅逆壓梯度被增強。t0+21 ms時，脫落空穴的影響基本消失，逆壓梯度強度減小，僅在吸力面0.75c的位置存在一個非常小的逆壓梯度，附著空化得以進一步發(fā)展，空穴長度達到0.7c。此時部分空化脫落較少，逆壓梯度對空穴生長的限制作用不明顯，吸力面空化會進一步生長成為超空化，從圖7(a)圖和圖7(b)中可以看出t0+28 ms和t0+35 ms時，翼型吸力面被空化覆蓋，表面也不存在壓力梯度。超空化后，空化開始大量脫落，不斷潰滅形成的局部高壓區(qū)再次加強逆壓梯度，逆壓梯度又會阻礙空穴生長，促進空化脫落，直到水翼表面的空穴消失，從t0+42 ms到t0+56 ms的3個時間段，逆壓梯度位置不斷向上游移動，強度也不斷增加，直至空穴完全脫落消失。以上的變化說明在一個周期內(nèi)，空穴的生長和脫落受逆壓梯度強度和位置的影響，而空穴的生長與脫落又影響著水翼表面的壓力分布。因而，翼型升力變化的周期與空穴的脫落周期是對應(yīng)的。 圖7 水翼表面升力系數(shù)與氣相分數(shù)分布Fig.7 Distribution of lift coefficient and vapor volume fraction around hydrofoil 圖8給出了空穴長度隨時間的變化。可見水翼表面的空化存在周期性脫落，脫落的周期為40 ms～50 ms，與升力變化周期56 ms[6]基本吻合，表明空化云周期性脫落是升力系數(shù)周期性變化的原因，這與前文中分析的結(jié)果是一致的。 圖8 空穴長度隨時間的變化Fig.8 Change of cavity length with time 2.2.2 反向射流的對比 圖9給出了一個周期內(nèi)翼型表面的速度矢量分布。可以看出翼型吸力面頭部速度最大，尾部速度最小并伴隨回流區(qū)。根據(jù)伯努利原理，水翼吸力頭部速度較大。隨著流體沿著翼型向下游運動，逆壓梯度會使流動出現(xiàn)分離，形成反向射流。t0+7 ms時，上一個周期的空化還未完全脫落，壓力分布上翼型表面存在較大的壓力梯度，速度分布上從圖9中可以看到在吸力面0.3c到尾部之間存在較強的回流區(qū)。t0+14 ms時，回流區(qū)移動到翼型尾部，即將脫離水翼，新周期的空化已經(jīng)初生。t0+21 ms時，隨著逆壓梯度的消失，回流區(qū)大部分消失在主流中，附著在翼型表面的回流區(qū)面積大大減小。t0+28 ms和t0+35 ms時，逆壓梯度作用不明顯，回流區(qū)面積基本沒有增加，使得部分空化發(fā)展成為超空化。在此之后，t0+42 ms時空穴充分發(fā)展，尾部的空化與回流區(qū)開始相互干涉?；亓髟鰪娏丝栈拿撀?，脫落空化潰滅瞬間的局部高壓造成的逆壓梯度又加強了回流，使得回流區(qū)的長度和厚度開始增加。另一方面，空化潰滅后反向射流會迅速填滿空泡破裂后的間隙，使得回流沿翼型向上游運動，進一步擴大回流區(qū)的范圍。t0+49 ms時，隨著回流厚度不斷增加，反向射流已經(jīng)運動到水翼的頭部，與空化區(qū)氣液交界面相互作用，造成空化云的大尺度脫落，完成空穴在一個周期內(nèi)的變化。對比圖7，反向射流的位置與逆壓梯度出現(xiàn)的位置是一一對應(yīng)的，說明反向射流和逆壓梯度共同影響空穴的生長與脫落。 圖9 水翼表面速度矢量分布Fig.9 Velocity vectors distribution around hydrofoil 本文采用標準k-ε模型和3種修正湍流模型對二維Clark-Y水翼空化流動進行計算，通過仿真結(jié)果的相互比較及與實驗結(jié)果的對比，得出以下結(jié)論： 1)與標準k-ε模型相比，修正的湍流模型可以大幅減小流場中的渦黏，使得雷諾應(yīng)力的部分歷史效應(yīng)得以凸顯，從而增強模型對于非定?？栈鲃拥念A(yù)測能力。 2)MFBM模型在3種修正模型中，捕捉的非定?？栈鲃犹匦耘c實驗最為接近，效果最好，能夠比較有效地模擬二維水翼非定常空化流動。2 結(jié)果與討論
2.1 仿真結(jié)果比較
2.2 水翼非定?？栈鲃臃治?/h3>
3 結(jié)論
猜你喜歡
應(yīng)用數(shù)學(xué)和力學(xué)(2024年1期)2024-01-29 08:32:46
鋼鐵釩鈦(2023年5期)2023-11-17 08:48:34
車用發(fā)動機(2021年5期)2021-10-31 05:48:38
海洋技術(shù)學(xué)報(2021年3期)2021-08-19 11:07:04
電源技術(shù)(2021年7期)2021-07-29 08:35:24
數(shù)字海洋與水下攻防(2020年5期)2021-01-04 01:40:20
廈門理工學(xué)院學(xué)報(2016年1期)2016-12-01 04:50:53
水利科技與經(jīng)濟(2016年7期)2016-04-25 13:03:00
噪聲與振動控制(2015年4期)2015-01-01 07:08:04
食品科學(xué)(2013年24期)2013-03-11 18:30:24
