孫鵬飛,徐鴻鵬,李濤,李大海,李寶童,洪軍

(1.西安交通大學機械工程學院,710049,西安;2.西安交通大學現(xiàn)代設計及轉子軸承系統(tǒng)教育部重點實驗室,710049,西安;3.西安航天動力試驗技術研究所,710100,西安)

渦輪流量計作為大推力液體火箭發(fā)動機燃料供給系統(tǒng)的關鍵精密測量組件,其測量精度對研制發(fā)動機、獲取發(fā)動機的準確特性具有非常重要的現(xiàn)實意義。作為典型的速度式流量計,渦輪流量計由流體推動渦輪旋轉,根據(jù)渦輪旋轉角速度與流體流速的比例關系,實現(xiàn)對流量的監(jiān)測。此外,渦輪流量計因具有重復性好、量程范圍寬、適應性強、精度高、體積小等特點,被廣泛應用于多種領域[1-2]。目前,渦輪流量計的測量性能通常采用水或柴油進行鑒定,而在大推力液體火箭發(fā)動機燃料供給系統(tǒng)中,渦輪流量計的工作環(huán)境為低溫液氧,介質和溫度的轉變會導致渦輪流量計的測量偏差。

為提升渦輪流量計的測量精度,常見的方法為優(yōu)化渦輪流量計的幾何構型,例如導流架的類型[3-4]、導頭形狀[5-6]、凹槽寬度[7]、導程長度[8]以及渦輪的構型[9-10]。其中,導流架通過對流體介質的流動方向進行導向,維持了流體介質的穩(wěn)定運動方向。對具有不同導頭的導流架研究結果表明,在大流量條件下,相比于傳統(tǒng)球形導流架,流線型導流架因具有更緩的漸變梯度導頭,有效降低了渦輪流量計運行管道內的速度梯度,提升了渦輪流量計的測量精度[5]。此外,合理減小導流架的凹槽寬度還能有效降低渦流引起較大的壓力損失[7],從而保證了流量計的測量精度。然而,在實際應用中,發(fā)動機的供液系統(tǒng)通過測量流量計的渦輪轉速,實現(xiàn)對流體介質流量的監(jiān)測。與導流架相比,優(yōu)化渦輪構型能進一步有效地提升其測量精度。在渦輪的優(yōu)化中,通常將渦輪的葉片數(shù)、輪轂半徑、導程等參數(shù)作為優(yōu)化變量[7,11]。例如,對于小口徑渦輪流量計,通過增大渦輪的輪轂半徑、輪轂長度和葉片導程,使其儀表系數(shù)線性度誤差由5.23%降低到4.69%,實現(xiàn)了渦輪流量計測量精度的提升[10]。可以渦輪流量計的渦輪導程、葉片數(shù)和前后導流架的輪轂長度為設計變量,在高黏度流體介質下進行響應面和正交試驗,得到了具有儀表系數(shù)線性度誤差最小化的優(yōu)化變量組合,使儀表系數(shù)線性度誤差下降了0.12%[12]。綜上所述,現(xiàn)有研究對渦輪流量計的優(yōu)化設計側重于小口徑流量計渦輪的整體構型和導流構件,且關于大口徑渦輪流量計的葉片構型對測量精度影響的研究較少。因此,有必要對大口徑渦輪流量計展開研究,提升其應用價值。

本文以大推力液體火箭發(fā)動機供氧系統(tǒng)中的試驗用大口徑DN600渦輪流量計為研究對象,采用六自由度(6DOF)模型[13-14]結合流體仿真方法,研究了渦輪流量計在水介質和液氧中的儀表系數(shù)等效關系,預估了渦輪流量計在液氧中的測量精度,并以儀表系數(shù)線性度誤差最小化為優(yōu)化目標,對渦輪葉片進行優(yōu)化設計。

1 研究對象及方法

1.1 研究對象

圖1為大推力液體火箭發(fā)動機供氧系統(tǒng)。其中,Dn為管道的內徑尺寸。本文以大推力液體火箭發(fā)動機中的DN600渦輪流量計為研究對象,對渦輪葉片構型進行優(yōu)化,實現(xiàn)對其測量精度的提升。根據(jù)渦輪流量計的平面參數(shù)化圖形,提取其葉片的基本結構參數(shù),如表1所示,并建立渦輪葉片三維結構模型,如圖2所示。此外,為便于流量計優(yōu)化中模型重構參數(shù)的選定,將其幾何參數(shù)保留一位有效位數(shù)。

表1 DN600渦輪流量計葉片基本結構參數(shù)

葉片的厚度、中徑來流角和寬度分別如圖3所示,中徑來流角為渦輪葉片中間段與流體介質流動方向的夾角。葉片寬度為葉片在x方向的跨度。葉片重合度即n個葉片在zoy平面上的投影面積與渦輪工作平面面積的比例系數(shù)。

圖1 大推力液體火箭發(fā)動機供氧系統(tǒng)Fig.1 The oxygen feeding system of high-thrust liquid rocket engine

圖2 渦輪葉片三維結構模型Fig.2 The three-dimensional model of turbine blades

在葉片優(yōu)化過程中,為便于獲得具有不同構型的葉片,根據(jù)其幾何參數(shù)的耦合關系,建立葉片數(shù)n、葉片重合度k、葉根偏角α1、葉頂偏角α2、葉根圓心角θ1、葉頂圓心角θ2、葉根半徑Rl、葉頂半徑Rt、中徑來流角α3和葉片寬度h的函數(shù)表征關系,如式(1)～式(4)所示。

葉根偏角α1、葉根圓心角θ1、葉根半徑Rl與葉片寬度h的關系為

θ1=2arcsin(htanα1/2Rl)

(1)

葉頂偏角α2、葉頂圓心角θ2、葉頂半徑Rt與葉片寬度h的關系為

θ2=2arcsin(htanα2/2Rt)

(2)

(a)葉片厚度 (b)中徑來流角 (c)葉片寬度

(d)葉根偏角 (e)葉頂偏角圖3 渦輪葉片俯視圖Fig.3 The top view of turbine blades

葉根偏角α1、葉頂偏角α2與中徑來流角α3的關系為

α3=(α1+α2)/2

(3)

葉根圓心角θ1、葉頂圓心角θ2及重合度k的關系為

θ1=θ2=360k/n

(4)

渦輪流量計運行機制為前后導流架固定,渦輪繞中心軸旋轉運動,故對仿真流域進行三段式建模,分別為進口上游區(qū)域、渦輪區(qū)域和出口下游區(qū)域,渦輪流量計的局部管道模型如圖4所示。DN600渦輪流量計沿x方向的總跨度為811.6 mm,葉頂直徑為588.0 mm。

圖4 渦輪流量計的局部管道模型Fig.4 The local pipeline model of the turbine flowmeter

1.2 數(shù)值模擬方法與矯正

渦輪流量計受到流體沖擊時,渦輪在流體驅動力矩作用下被動旋轉。因此,為得到不同流量下渦輪流量計的轉速,分析其儀表系數(shù)線性度誤差,本文采用ANSYS/Fluent商業(yè)軟件,分別對渦輪流量計模型展開定常與非定常分析。由渦輪的運行機制可知,其初速度為0,在定常分析中僅需設定入口流速和出口壓強進行仿真計算,在本文中出口均定義為自由出口。在非定常分析時,通過6DOF模型將流量計渦輪區(qū)域的運動形式定義為被動,自由度定義為繞x軸的單一轉動自由度[14],并將定常分析結果作為初始值進行求解[15-17]。

1.2.1 網格無關性驗證

在有限元分析中,不同數(shù)量級的網格對物理場結果的影響不可忽視。隨著有限元網格數(shù)增加,分析結果越精確,但是求解過程也越復雜、計算效率越低。因此,需綜合考慮仿真計算精度、計算機的運算能力和工時等多方面因素的影響,選取合適的網格數(shù)以保證仿真方法的可靠性和經濟性。在水介質條件下,本文分別采用121萬、160萬、200萬、250萬和334萬5種數(shù)量級的網格進行有限元仿真分析。水介質仿真網格無關性驗證如圖5所示。

圖5 水介質仿真網格無關性驗證Fig.5 The mesh independence validation

由圖5可知,出口壓力隨網格數(shù)的增加而逐漸上升,當網格數(shù)達到250萬以后,出口壓力趨于穩(wěn)定,且不斷增加網格數(shù)對于仿真結果沒有實質性的影響,可認為當網格數(shù)達到250萬以后,仿真結果與網格無關。因此,在水介質條件下,均采用250萬數(shù)量級的網格劃分策略,管道模型的進口上游區(qū)域和出口下游區(qū)域采用四面體單元的體網格,網格大小為35 mm;渦輪區(qū)域采用四面體單元的面網格,網格大小為5 mm;前后導流架的壁面采用四面體單元的面網格,網格大小為7 mm。渦輪流量計的流體域有限元模型如圖6所示。

圖6 渦輪流量計的流體域有限元模型Fig.6 The fluid-field finite element model of the turbine flowmeter

1.2.2 數(shù)值模擬方法的矯正

在運行過程中,渦輪流量計的渦輪因流體驅動力矩Td、流體黏性阻力矩Trf、軸承摩擦機械阻力矩Tb和電磁阻力矩Tm的作用而發(fā)生旋轉運動,其轉動的微分方程[18]如下式

(5)

式中:J為渦輪的轉動慣量;ω為渦輪穩(wěn)定狀態(tài)時的轉速。

為分析渦輪流量計的力矩,從葉片任意半徑r處將渦輪展開成直列葉珊,葉片入口和出口的速度三角形如圖7所示。

圖7 葉片入口和出口的速度三角形Fig.7 The velocity triangle of blades in inlet and outlet

葉片沿高度方向受到的驅動力矩微分方程為

dTd=rdF

(6)

式中:dF為葉片半徑r處的驅動力,由動量定理可得知

dF=(vtanθi-rω)2πρvrdr

(7)

式中:v為流體流速;ρ為流體介質密度。tanθi的表達式如下

(8)

由式(8)可得

(9)

從渦輪葉根半徑到葉頂半徑對式(9)積分得到總驅動力矩,如下式

(10)

(11)

由式(11)可得,渦輪流量計所受到的驅動力矩為

(12)

式中:A為進出口通流面積;Q為流體流量。

本文流體黏性阻力矩分別為渦輪輪轂表面的黏性阻力矩Th、葉片表面黏性阻力矩Ts以及葉頂與管壁間的黏性阻力矩Tt。

根據(jù)順流放置平板表面的阻力研究結果,渦輪輪轂表面的黏性阻力矩[19]為

(13)

Ah=2πRhLh-Ntbhch

(14)

(15)

(16)

式中:Vzh為輪轂位置流體的軸向速度;tbh為輪轂位置葉片厚度;U∞ch為輪轂位置出入口流體的平均相對速度;Rh為輪轂半徑;Lh為渦輪的導程;ch為輪轂位置的弦長。本文中Re均表示雷諾數(shù)。

葉片表面黏性阻力矩的表達式為

(17)

At=2Nctr

(18)

(19)

(20)

式中ct為葉頂弦長。

根據(jù)同軸圓筒壁面間的摩擦力矩,葉頂與管壁間的黏性阻力矩如下

(21)

(22)

式中tbt為葉頂?shù)暮穸?本文中即為葉片厚度h。

在實際條件下,轉軸與軸承間的油膜不充分導致潤滑不足,從而存在機械摩擦阻力矩。在渦輪開始旋轉或逐漸停止旋轉時,軸與軸承間存在黏性與機械摩擦共存的混合狀態(tài),且電磁阻力矩相對較小,對渦輪流量計的影響很小,難以精確獲取機械摩擦阻力矩和電磁力矩。因此,本文在6DOF模型中施加預設阻力矩以等效機械阻力矩和電磁阻力矩,對有限元仿真模型進行矯正。

為獲得DN600渦輪流量計在不同流量下的儀表系數(shù)真實值,在水介質條件下,對不同流量下的DN600渦輪流量計儀表系數(shù)分別進行多次試驗測量,渦輪流量計試驗測量平臺如圖8所示。通過數(shù)值平均得到流量計儀表系數(shù),如表2所示。根據(jù)表2中的結果,通過下式分別計算仿真模型在不同流量下的入口流速

(23)

式中R為流體域入口截面半徑。

表2 渦輪流量計試驗測量數(shù)據(jù)

根據(jù)下式分別計算不同流速的雷諾數(shù)[20-21],分析流場的狀態(tài)

(24)

式中:D為水力直徑,對于圓形管道,水力直徑即為圓形管道橫截面的直徑;μ為流體介質的動力黏度。

圖8 渦輪流量計試驗測量平臺Fig.8 The experimental measurement platform of the turbine flowmeter

供氧管道的入口流速和雷諾數(shù)如表3所示。由表3可知,所得到的雷諾數(shù)的數(shù)量級均為106,流場狀態(tài)可視為湍流狀態(tài)。此外,根據(jù)渦輪流量計的運行機制分析,計算模型可選用標準的k-ε模型[22]。

表3 供氧管道的入口流速和雷諾數(shù)

根據(jù)渦輪穩(wěn)定狀態(tài)的轉速,可計算得到渦輪旋轉頻率f,如下式

(25)

基于渦輪旋轉頻率與流量Q的關系,可得到渦輪流量計儀表系數(shù)K,如下式[23]

(26)

根據(jù)表3的入口水介質流速條件,對渦輪流量計進行仿真分析,提取渦輪的穩(wěn)定轉速并計算相應的儀表系數(shù)。渦輪流量計仿真與試驗結果對比如表4所示。

表4 渦輪流量計仿真與試驗結果對比

由表4可知,本文所采用的方法與實測數(shù)據(jù)間的誤差均小于0.9%。誤差是因為6DOF模型中未考慮摩擦阻力矩和電磁阻力矩的作用。然而,根據(jù)實際工程需求,儀表系數(shù)誤差應不超過0.5%。因此,基于上述數(shù)值結果,在6DOF模型中施加預載荷力矩對仿真模型進行誤差矯正。通過在不同流量條件下施加多組預載荷力矩進行對比分析,得到如表5所示的矯正結果。

表5 不同流量點下的矯正預載荷力矩

2 水介質和液氧的等效分析

在實際應用中,液體火箭發(fā)動機的真實工作環(huán)境難以模擬,且采用液氧進行試驗以獲得渦輪流量計儀表系數(shù)的成本高。因此,本文通過數(shù)值模擬方法,獲得渦輪流量計在水介質和液氧下的儀表系數(shù)等效關系,實現(xiàn)對液氧條件下渦輪流量計測量精度的預估。液氧的物性參數(shù)如表6所示。

表6 液氧的物性參數(shù)

2.1 渦輪流量計熱變形分析

由于液氧溫度為-183 ℃,渦輪流量計在液氧中會發(fā)生冷縮變形,影響其測量精度,需對渦輪流量計進行熱變形分析。因此,本文采用熱固耦合的仿真方法揭示渦輪流量計的熱變形情況。渦輪流量計有限元模型如圖9(a)所示,網格類型為6 mm四面體單元的面網格。渦輪流量計的溫度載荷邊界和約束邊界分別如圖9(b)和9(c)所示,在渦輪流量計外表面施加-183 ℃的溫度載荷,并固定其中心軸。

渦輪流量計的熱應變云圖如圖10所示。由圖10可知,渦輪流量計的最大熱變形主要集中在導流架的葉片頂部。

2.2 液氧條件下的仿真分析

(a)渦輪流量計有限元模型

(b)溫度載荷邊界

(c)約束邊界圖9 渦輪流量計有限元模型與邊界條件Fig.9 The finite element model and boundary conditions of the turbine flowmeter

圖10 渦輪流量計的熱應變云圖Fig.10 The total thermal deformation of the turbine flowmeter

首先提取熱變形模型的幾何信息,對其進行模型重構。然后,在液氧環(huán)境下,對重構的渦輪流量計模型進行流場有限元分析,獲得渦輪流量計的儀表系數(shù)。為保證仿真結果的正確性,分別采用了116萬、193萬、252萬、338萬和400萬5種數(shù)量級的網格進行有限元仿真,分析出口壓力的變化以進行液氧仿真網格無關性驗證,如圖11所示。

圖11 液氧仿真網格無關性驗證Fig.11 The mesh independent validation

由圖11可知,當網格數(shù)達到338萬時,出口壓力趨于穩(wěn)定,且不斷增加網格數(shù)對于仿真結果沒有實質性的影響,可認為當網格數(shù)達到338萬后,仿真結果與網格無關,滿足網格無關性的要求。因此,采用338萬數(shù)量級的網格劃分策略,分別在流量為3 599、3 208、2 901和2 507 m3/h的液氧環(huán)境下進行仿真分析,獲取相應的儀表系數(shù)。在338萬的網格劃分策略中,管道模型的入口上游區(qū)域和出口下游區(qū)域采用27 mm四面體單元的體網格;渦輪區(qū)域采用5 mm四面體單元的面網格;前后導流架的壁面采用5 mm四面體單元的面網格。渦輪流量計在液氧與水介質條件下的儀表系數(shù)如圖12所示。由圖12可知,在液氧和水介質環(huán)境下,DN600渦輪流量計的儀表系數(shù)成正相關,比例系數(shù)近似為0.999 5。

圖12 水介質和液氧介質仿真儀表系數(shù)Fig.12 The meter coefficient in water and liquid oxygen simulations

3 渦輪流量計葉片優(yōu)化

為提高渦輪流量計的測量精度,對渦輪葉片的中徑來流角、葉片寬度和葉片重合度進行優(yōu)化。根據(jù)實際工程經驗,優(yōu)化參數(shù)取值的最大變化范圍為30%,且相較于原始構型,優(yōu)化的葉片構型變化不宜過大。因此,本文分別取中徑來流角為45°、50°和55°,葉片寬度為45 mm、55 mm和65 mm,重合度為0.9、1.0和1.1。通過正交試驗,分別使用α3、h、k表示中徑來流角、葉片寬度和重合度3個因素,如表7所示[24]。每個因素有3個水平,分別為1、2、3。選用L9(34)正交表安排試驗,如表8所示。

表7 渦輪葉片參數(shù)的因素水平表

由圖1可知,在渦輪流量計的真實管道模型中,存在的偏心流會影響流量計的測量精度。因此,為避免因偏心流對優(yōu)化結果的影響,在仿真過程中采用直管道進行有限元仿真分析,將圖4中的入口上游區(qū)域定義為10.5Dn,出口下游區(qū)域定義為10.5Dn。

表8 渦輪流量計正交試驗表

本文采用儀表系數(shù)線性度誤差以評估渦輪流量計的測量精度,如下式[25]

(27)

式中:Kmax,i為渦輪流量計在不同流量下得到的儀表系數(shù)最大值;Kmin,i為渦輪流量計在不同流量下得到的儀表系數(shù)最小值。

渦輪流量計正交試驗結果與均值如表9和表10所示。由表10可知,α3對應的第3水平均值最小,表明α3取第3水平上的值時線性度誤差最小,同理可以得到h和k的取值均為第2水平。根據(jù)渦輪流量計在液氧和水介質中的儀表系數(shù)比例關系可知,理論上在液氧環(huán)境中的最優(yōu)水平組合為(α3)3h2k2,即中徑來流角、葉片寬度和重合度的優(yōu)化組合為55°、55 mm、1。

表9 渦輪流量計正交試驗結果

表10 渦輪流量計的線性度誤差均值

為對比優(yōu)化前后渦輪流量計的測量精度,分別對優(yōu)化所得的渦輪流量計進行有限元分析,得到不同流量下的儀表系數(shù),如圖13所示。根據(jù)圖13中結果計算相應的儀表系數(shù)線性度誤差,結果表明,優(yōu)化的渦輪流量計儀表系數(shù)線性度誤差為0.277 0%,相比于原始結構的0.381 5%,降低了0.104 5%,顯著提升了渦輪流量計的測量精度。

圖13 優(yōu)化前后不同流量下的儀表系數(shù)對比Fig.13 The comparison of the initial and optimized meter coefficients

(a)優(yōu)化渦輪

(b)原始渦輪圖14 渦輪速度分布云圖Fig.14 The velocity distribution of the turbine

為進一步揭示優(yōu)化后渦輪流量計的測量精度,優(yōu)化前后渦輪流量計在3 599 m3/h流量下的渦輪速度分布和流場截面速度分布云圖分別如圖14和15所示。由圖14可知,渦輪的最大轉速在葉片頂部,相對于原始渦輪流量計,優(yōu)化的渦輪流量計的轉速更大,運轉更加順暢。由圖15可知,流體在流經優(yōu)化的渦輪葉片時,流場速度分布較原始流量計更加均勻。綜上所述,優(yōu)化的流量計渦輪結構在流場中的旋轉穩(wěn)定性更好,從而使測量精度得到提高。

(a)優(yōu)化結構

(b)原始結構圖15 流場截面速度分布云圖Fig.15 The velocity distribution of the flow field

4 結 論

本文以大推力液體火箭發(fā)動機供氧系統(tǒng)中的試驗用DN600渦輪流量計為研究對象,通過數(shù)值模擬獲得了渦輪流量計在水介質和液氧環(huán)境中儀表系數(shù)的等效關系。以儀表系數(shù)線性度誤差最小為優(yōu)化目標,通過正交仿真試驗優(yōu)化了渦輪流量計葉片,主要結論如下。

(1)數(shù)值仿真結果表明,渦輪流量計在水介質和液氧環(huán)境下的儀表系數(shù)呈正相關,為液氧流量計的儀表系數(shù)修正提供了數(shù)據(jù)參考。

(2)渦輪流量計中徑來流角、葉片寬度和重合度的優(yōu)化組合為55°、55 mm、1;相比于原始構型,優(yōu)化的DN600渦輪流量計的線性度誤差降低了0.104 5%,測量精度顯著提升。