張凱宏，江 欣，肖明杰，梁樹強 ，胡 偉，冀 鵬

(1.西安航天動力研究所，液體火箭發(fā)動機重點實驗室，陜西 西安 710100;2.航天推進技術研究院，陜西 西安 710100)

0 引言

隨著航天技術的發(fā)展，航天器對其推進系統(tǒng)的技術指標要求也越來越高。特別是對于起軌道機動功能的軌控動力系統(tǒng)，其推力越來越大、室壓越來越高、響應時間越來越短[1-3]，導致管路系統(tǒng)流量大、壓力高，關機時形成的管路水擊壓力急劇增加，嚴重時引起管路的泄漏和破壞。因此，了解管路系統(tǒng)水擊特性，采取措施減小水擊峰值壓力成為軌控發(fā)動機管路系統(tǒng)設計中需解決的突出問題。

水擊[4-6]也稱液壓沖擊或“水錘”現(xiàn)象，是由于液體局部加速度過大所致的管內壓強急劇變化的一種動態(tài)現(xiàn)象，這種現(xiàn)象經常出現(xiàn)在液體火箭發(fā)動機的起動、關機以及轉工況過程中。尤其在關機過程中出現(xiàn)的水擊峰值壓力可以達到管路正常工作壓力的幾倍甚至幾十倍，這種突然的壓力大幅升高會導致管路系統(tǒng)破裂、推進劑泄漏，進而導致推進系統(tǒng)工作異常。水擊的數值計算方法主要包括特征線法和有限體積法等。特征線法[7-8]理論比較完善，且不受管道系統(tǒng)復雜程度的限制，但不利于解決多特征插值與非積分項誤差的問題。有限體積法[9-10]能較好地反映水擊壓強，同時克服了特征線法計算的水擊壓強波動幅度較大、衰減較慢的問題。

水擊波在沿管道運動的過程中，會誘發(fā)管路結構振動，這一振動又影響壓力的波動過程和流體的流動過程，形成以流體的瞬時流動、流體波動以及結構振動為一體的振動系統(tǒng)。3種運動形式間的相互作用稱為流固耦合[11-13]。流固耦合效應會影響水擊的峰值壓力和衰減速率，基于流固耦合對水擊特性進行分析能得到更準確的結果。對于水擊壓力比較高的環(huán)境，忽略結構變形對水擊的影響，將會對水擊仿真帶來較大誤差。近年來，國內外學者開始利用流固耦合對水擊等問題進行仿真。2014，Cesana等[14]對兩根同軸管路的水錘響應進行了流固耦合分析，并驗證了基于流固耦合模擬水擊壓強的準確性。2015年，朱雨[15]將入口設置為脈沖流體速度入口邊界條件，研究了管路水擊壓力脈動誘導結構振動的流固耦合相互作用機理和結構動力學響應規(guī)律。

為了驗證考慮流固耦合的水擊仿真適用性，本文基于Workbench協(xié)同仿真環(huán)境[16-17]，分別在Fluent仿真軟件和Transient Structural軟件內設置流體和固體相對應的作用面，并在System Coupling內實現(xiàn)數據的相互傳遞。

1 物理模型和計算方法

1.1 流體控制方程

閥門關閉過程是一個動態(tài)過程，管道中流體為非穩(wěn)態(tài)湍流流動，控制方程如下：

1)連續(xù)方程

(1)

2)動量方程

根據運動守恒定律，管路流體的運動方程的基本形式可寫成：

(2)

3)湍流模型

采用標準k-ε湍流模型來使上述方程封閉。

k方程：

(3)

ε方程：

(4)

1.2 固體控制方程

管路結構動力學控制方程如下：

(5)

σ=FST

(6)

(7)

S=λsTr(E)I+2μsE

(8)

式中:σ為應力張量；F為變形梯度；S為二階Piola應力；E為格林-拉格朗日應變；I為單位矩陣；λs和μs為Lame常數；ρs為管路材料的密度，kg/m3；fs為體積力。

1.3 流固耦合方程

流固耦合遵循力和位移的守恒原則，所以在流固耦合交界面處，應該保證流體與固體結構的應力、位移等基本參數數值相等。即滿足如下2個方程：

(9)

式中：下標f為流體；s為固體。

2 算例驗證

2.1 算例驗證模型

本文基于文獻[18]中的管路在閥門關閉時所產生的水擊過程進行基于流固耦合的數值研究。試驗管道內徑10 cm，壁δ=7 mm，長度L=24.47 m,沿程損失系數由試驗測得λ=0.033 1，文獻中沒有說明閥門的具體形式，在此選擇比較常見的球閥。如圖1所示取球閥閥芯及前后直管段中全部流體為計算區(qū)域。靠近管壁處的網格用邊界層網格加密，如圖2所示。

圖1 模型及球閥結構Fig.1 Model and structure of ball valve

圖2 管橫截面網格Fig.2 Mesh of profile of pipe

2.2 仿真環(huán)境及設置

為了實現(xiàn)數值模擬管路水擊雙向流固耦合作用過程，本文在研究中利用 Workbench15.0 平臺，通過聯(lián)立 Fluent + Transient Structural + System Coupling三部分分析模塊構建雙向流固耦合分析系統(tǒng)。流固耦合仿真平臺如圖3所示。

圖3 流固耦合仿真平臺Fig.3 The simulation platform of FSI

閥門的開關過程采用Fluent中的滑移網格技術實現(xiàn)，將球閥與左右管路流體部分的接觸面定義為interface交界面。在閥芯的轉動過程中，兩個交界面互相接觸的部分為interior內部面，流體可以通過該面進行流動，不接觸的部分設置為不可流通的固體壁面(wall)，沒有流體通過。當兩個交界面完全脫離時，全部交界面成為壁面，以此來模擬閥門完全關閉時的狀態(tài)。閥門的轉動規(guī)律使用UDF控制，閥門從0.2 s開始作動，以1.75 rad/s的角速度繞中心軸勻速旋轉0.77 s，從而實現(xiàn)閥門關閉過程。

根據試驗數據，管道進口截面壓力保持恒定，定義為壓力入口，總壓16.4 kPa。出口截面定義為壓力出口，靜壓0 Pa。管壁粗糙度設置為0.000 6，流體為可壓縮液體水，設置管壁為動網格，動網格類型為system coupling。

結構計算中，設置管內壁為流固耦合面。system coupling設置中，設置流體與固體的耦合面實現(xiàn)力與位移的數據傳遞。

2.3 計算結果分析

水擊過程的試驗與仿真的壓強變化如圖4所示，為了評估仿真結果與試驗結果的吻合程度，利用峰值壓力、衰減速率、水擊頻率的誤差來進行評估。

峰值壓力誤差：

(10)

衰減速率誤差：

(11)

水擊頻率誤差：

(12)

式中：下標1為試驗值；下標2為仿真值；pm為峰值壓力，定義為水擊最大壓力與水擊前初始壓力之差；t為衰減時間，由于穩(wěn)態(tài)壓力相對于峰值壓力較小，認為當壓力在穩(wěn)態(tài)壓力的±20%之間達到穩(wěn)態(tài)，衰減時間為閥門開始作動到達到穩(wěn)態(tài)所用的時間；f為水擊頻率，用峰值壓力到第二個波谷的時間(1.5個周期)來計算水擊頻率。

從圖4中可以看出，基于流固耦合的水擊仿真可以較好的模擬水擊壓力第一次上升和下降的過程，與試驗曲線吻合度較高。表1為仿真結果與試驗的誤差值。從表1中看出，仿真值與試驗值的峰值壓力誤差和衰減速率誤差分別為7.47%，5.56%?；诹鞴恬詈系乃畵舴抡婵梢暂^好地模擬水擊的峰值壓力和衰減過程。仿真值與試驗值的水擊頻率誤差為27.47%，雖然誤差超過20%，但依然可以看出水擊過程中的振蕩趨勢。因此，基于流固耦合的水擊仿真，可以得出較為一致的水擊特性。

表1 仿真值與試驗值誤差

圖4 仿真與試驗水擊壓強變化曲線Fig.4 Water hammer pressure of simulation and experiment

3 水擊特性

3.1 模型及邊界條件設置

根據軌控發(fā)動機管路短、半徑小的結構特點與壓力高、流量大、響應時間短工作特點，對研究管路的尺寸和工作條件作如下更改：減小管路直徑、壁厚和管路長度，直徑為0.028 m，壁厚為0.001 m，管路長度為0.5 m；增加出口背壓和流量，背壓設置為1 MPa，4 MPa，6 MPa，流量設置為7.1±0.1 kg/s，7.55±0.1 kg/s，8.0±0.1 kg/s，8.45±0.1 kg/s；縮短閥門關閉時間，閥門關閉時間為0.02 s，0.04 s，0.06 s。研究共設置8個工況，具體設置如表2所示。

仿真環(huán)境設置分為結構組件設置、流體組件設置、耦合組件設置，時間步長均為0.001 s，結束時間0.2 s。結構組件設置：選擇workbench中Transient Structural計算結構，設置管路內壁面為流固耦合面。流體組件設置：選擇workbench中Fluid Flow(Fluent)計算流體，type選擇Pressure-Based，time選擇Transient，model選擇Standardk-ε湍流模型，Materials設置為water-liquid (compressible-liquid)，閥門關閉使用udf控制閥芯運動，閥芯從0.05 s開始作動；進口為pressure-inlet，總壓根據流量要求設置，出口為pressure-outlet，靜壓根據要求設置，溫度均設置為298 K；使用interface實現(xiàn)不同域之間數據的連通，將管壁設置為system-coupling動網格；Solution Method設置為Coupled Second Order；時間步長0.001 s，計算時間步數200，每時間步計算200步。耦合組件設置：將結構域的內壁和流體域的管壁創(chuàng)建數據傳遞，實現(xiàn)位移與壓力的數據傳遞。

表2 不同工況工作條件

3.2 水擊特性分析

圖5為工況1不同時刻壓力和流線分布圖。0.05 s時，閥芯開始轉動，此時壓力分布均勻，流動趨于穩(wěn)定。0.06 s時，閥芯轉動25%，此時閥門處壓力升高，壓力波向入口處傳播。0.07 s時，閥芯轉動50%，壓力急劇上升，閥芯前管路出現(xiàn)渦，但渦流較小，壓力波已經傳播到入口處，入口處壓力升高。0.08 s時，閥門轉動75%，閥芯前管路渦增大，壓力波繼續(xù)向入口處傳遞，壓力繼續(xù)升高，此時貼近壁面處未出現(xiàn)逆流。0.085 s時，閥芯轉動87.5%，壓力波開始向閥門處傳遞，壓力開始降低，閥芯前管路渦增大，在管路中間區(qū)域，在壁面位置存在逆流，在壁面位置首次出現(xiàn)細長渦流。0.09 s時，閥門完全關閉，壓力波繼續(xù)向閥門處傳遞，壓力繼續(xù)降低，壁面位置渦繼續(xù)增加。通過壓力分布的變化可以看出壓力波的傳播過程，先從閥芯處傳遞至入口處，此時壓力升高；之后壓力波從入口處傳遞至閥芯處，此時壓力降低。從圖中可以看出，在水擊發(fā)生過程中，管路中心位置為正向流動，貼近壁面的流動為逆向流動，因此形成渦流。

圖5 不同時刻壓力和流線分布圖Fig.5 Pressure and streamline at different time

圖6為0.05 s時入口壓力分布圖，此時流動速度較高。圖7為閥門關閉后入口壓力分布圖。0.085 s時，中心處速度約為0.2 m/s，邊界層為逆流，與流線分布圖吻合。0.09 s時，入口處全部為逆流，最大速度達到-0.4 m/s，且邊界層處逆流速度更大。0.1 s時，入口中心處流動變?yōu)檎蛄鲃樱吔鐚右廊粸槟媪鳌?.115 s時，中心處流動為逆向流動。0.125 s時，中心處流動為正向流動。0.2 s時，入口處流動速度趨于0。分析可知，中心處流動出現(xiàn)正向流動和逆向流動交替出現(xiàn)的現(xiàn)象，這和壓力波的傳遞有關。中心處速度的絕對值在逐漸減小，到0.2 s時，流動速度趨于0。從0.085 s開始，邊界層流動一直為逆向流動，流動速度與中心處流動速度和方向有關，中心處正向流動時，邊界層流速較小，中心處速度越大，邊界層流速越??；中心處逆向流動時，邊界層流速較大，中心處逆流速度越大，邊界層流速越大。

圖6 0.05 s時入口速度分布圖Fig.6 Velocity distribution of inlet at 0.05 s

圖8為工況1條件下，考慮流固耦合的水擊仿真和不考慮流固耦合的CFD水擊仿真水擊壓強變化曲線。從圖中可以看出，考慮流固耦合仿真結果的峰值壓力比不考慮流固耦合的水擊峰值壓力小3.1%。這是因為在水擊發(fā)生過程中，由于流體壓力的增加，管路截面積變大，管壁變形吸收了一部分水擊能量，使得水擊峰值壓力減小。

圖7 不同時刻入口速度分布圖Fig.7 Velocity distribution of inlet at different time

圖8 CFD和FSI計算方法水擊壓強變化曲線Fig.8 Water hammer pressure at CFD and FSI

綜合上述分析，水擊的能量耗散主要有以下幾種：通過動能和壓力能的方式從入口逆流出管路；流體的粘性損失；通過結構的變形吸收能量。由圖5、圖7中可知，入口處壓力一直處于較低水平，且入口處逆流速度很小，因此只有很少一部分水擊能量通過動能和壓力能的方式從入口逆流出管路。而結構吸收的能量使水擊壓力減小了3.1%，因此可以認為大部分水擊能量的耗散是由于流體的粘性損失。從圖5、圖7可以看出：在壓力下降過程中，管內出現(xiàn)了渦流，且壓力下降到最低時，渦流最大；管壁附近一直為逆流，且逆流速度相對中心處較高。因此，通過渦流和壁面處摩擦導致的黏性損失較大，是水擊能量耗散的主要原因。

圖9為不同流量下水擊壓強變化曲線，從0.05 s閥芯作動開始，壓強開始急劇增加，在0.075 s時壓力達到峰值壓力，0.096 s時壓力達到最低值，經過兩次的振蕩達到穩(wěn)態(tài)。不同流量的水擊變化趨勢一致，衰減速率和頻率一致，水擊的峰值壓力隨流量的增大而增大。

圖9 不同流量下水擊壓強變化曲線Fig.9 Water hammer pressure at different mass flow rate

水擊的峰值壓力與流量的關系圖如圖10所示，峰值壓力的變化和流量呈一次函數關系。水擊峰值壓力計算公式為：

(13)

式中：ρ為液體密度；a為聲速；u為流動速度；T為水擊周期；t為閥門關閉時間。

圖10 峰值壓力與流量關系圖Fig.10 Relationship between pressure and mass flow rate

水擊峰值壓力與流速成正比，與仿真結果一致。

閥門關閉時間是影響水擊特性的關鍵因素，圖11為不同閥門關閉時間下水擊壓強變化曲線。在閥門關閉時間為0.02 s條件下，0.05 s時閥芯開始作動，0.065 s時壓力達到峰值壓力，之后閥門壓力急劇降低，0.078 s時壓力達到最小值，經過2.5個周期后壓力趨于穩(wěn)態(tài)值。在閥門關閉時間為0.04 s條件下， 0.075 s時壓力達到峰值壓力，之后閥門壓力急劇降低，0.098 s時壓力達到最小值，經過2個周期后壓力趨于穩(wěn)態(tài)值。在閥門關閉時間為0.06 s條件下， 0.085 s時壓力達到峰值壓力，之后閥門壓力急劇降低，0.122 s時壓力達到最小值，經過1個周期后壓力趨于穩(wěn)態(tài)值。閥門關閉時間的縮短增加了峰值壓力和水擊頻率，增加了衰減時間，對管路造成惡劣的影響。閥門關閉時間短，會使水擊峰值壓力急劇增加，水擊頻率增大，衰減速率變慢。圖12表示峰值壓力p與關閥時間的倒數1/t的關系，從圖中可以看出，p與1/t呈一次函數關系，與式(13)相吻合。

圖11 不同關閥時間下水擊壓強變化曲線Fig.11 Water hammer pressure at different valve closure time

圖12 峰值壓力與關閥時間關系圖Fig.12 Relationship between pressure and valve closure time

圖13為不同背壓對水擊特性的影響，在流量、關閥時間等其他條件相同、背壓不同的情況下，峰值壓力、水擊頻率、衰減速率誤差低于2%。因此，在進行軌控發(fā)動機高室壓水力試驗時，可以在保證流量和關閥時間情況下，減小出口背壓，可以得出與高背壓一致的水擊壓力變化曲線。

圖13 不同背壓下水擊壓強變化曲線Fig.13 Water hammer pressure at different backpressure

4 結論

利用基于流固耦合的數值模擬方法，對算例進行仿真。驗證了使用Workbench15.0 平臺構建雙向流固耦合分析系統(tǒng)仿真的可行性?；谲壙匕l(fā)動機高室壓、大流量、快響應的發(fā)展趨勢，設計8組計算工況，分析了背壓、流量、關閥時間對水擊特性的影響，得到了以下結論：

1)Workbench15.0平臺構建雙向流固耦合分析系統(tǒng)不僅可以模擬水擊發(fā)生過程中壓力變化，還可以通過壓力和流線的分布圖分析壓力波的傳播過程和能量的耗散過程。

2)水擊的能量只有很少一部分通過動能和壓力能的方式從入口逆流出管路和通過結構的變形吸收能量這兩種方式耗散。大部分水擊能量的耗散是由于流體的粘性損失。粘性損失主要有渦流和壁面處摩擦兩種方式。

3)流量是影響水擊峰值壓力關鍵因素，在管路截面積不變的情況下，流量越大，水擊峰值壓力越大，且流量與水擊峰值壓力呈一次函數關系。流量對水擊頻率和衰減速率沒有影響。閥門關閉時間的縮短增加了峰值壓力和水擊頻率，增加了衰減時間，對管路造成惡劣的影響。峰值壓力p與關閥時間的倒數1/t呈一次函數關系。

4)在流量、關閥時間等其他條件相同、背壓不同的情況下，峰值壓力、水擊頻率、衰減速率誤差低于2%。因此，在進行軌控發(fā)動機高室壓水力試驗時，在保證流量和關閥時間情況下，減小出口背壓，可以得出與高背壓較為一致的水擊壓力變化曲線。

對于更為復雜的軌控發(fā)動機管路系統(tǒng)，基于流固耦合的三維水擊數值仿真不僅可以得出水擊壓力的變化，還可以得出水擊發(fā)生過程中流動過程和水擊波的傳遞過程。同時，三維仿真不需要對分支管路、彎管、變截面管路、異型截面管路進行簡化，更加保證了仿真的準確性。本文作為前期工作，為基于流固耦合的水擊仿真方法在復雜管路系統(tǒng)中的應用提供了參考。