基于流固耦合理論的關(guān)機水擊特性

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

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

0 引言

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

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

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

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

1 物理模型和計算方法

1.1 流體控制方程

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

1)連續(xù)方程

(1)

2)動量方程

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

(2)

3)湍流模型

采用標(biāo)準(zhǔn)k-ε湍流模型來使上述方程封閉。

k方程：

(3)

ε方程：

(4)

1.2 固體控制方程

管路結(jié)構(gòu)動力學(xué)控制方程如下：

(5)

σ=FST

(6)

(7)

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

(8)

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

1.3 流固耦合方程

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

(9)

式中：下標(biāo)f為流體；s為固體。

2 算例驗證

2.1 算例驗證模型

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

圖1 模型及球閥結(jié)構(gòu)Fig.1 Model and structure of ball valve

圖2 管橫截面網(wǎng)格Fig.2 Mesh of profile of pipe

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

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

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

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

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

結(jié)構(gòu)計算中，設(shè)置管內(nèi)壁為流固耦合面。system coupling設(shè)置中，設(shè)置流體與固體的耦合面實現(xiàn)力與位移的數(shù)據(jù)傳遞。

2.3 計算結(jié)果分析

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

峰值壓力誤差：

(10)

衰減速率誤差：

(11)

水擊頻率誤差：

(12)

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

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

表1 仿真值與試驗值誤差

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

3 水擊特性

3.1 模型及邊界條件設(shè)置

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

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

表2 不同工況工作條件

3.2 水擊特性分析

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

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

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

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

圖8為工況1條件下，考慮流固耦合的水擊仿真和不考慮流固耦合的CFD水擊仿真水擊壓強變化曲線。從圖中可以看出，考慮流固耦合仿真結(jié)果的峰值壓力比不考慮流固耦合的水擊峰值壓力小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

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

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

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

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

(13)

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

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

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

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

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

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

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

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

4 結(jié)論

利用基于流固耦合的數(shù)值模擬方法，對算例進(jìn)行仿真。驗證了使用Workbench15.0 平臺構(gòu)建雙向流固耦合分析系統(tǒng)仿真的可行性。基于軌控發(fā)動機高室壓、大流量、快響應(yīng)的發(fā)展趨勢，設(shè)計8組計算工況，分析了背壓、流量、關(guān)閥時間對水擊特性的影響，得到了以下結(jié)論：

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

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

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

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

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