水击问题的Fourier谱方法计算

陈宏玉，刘红军，刘 上

（西安航天动力研究所，陕西 西安710100）

0 引言

液体火箭发动机关机时，推进剂供应系统由于阀门关闭或涡轮泵减速等因素影响，不可避免地会在管路中形成水击和压力振荡。这种水击压力波动会造成诸如推进剂泄漏、结构破坏和推力输出精度降低等一系列问题，在发动机设计中需要对其进行深入分析。

水击问题可由一维瞬变流动基本方程描述，它们是一组非线性双曲型偏微分方程，一般情况下不存在解析解，只能采用数值计算方法求解。特征线法由于其具有较高的计算精度，数值方法简便，运算稳定性较好和易于编程等特点广泛应用在水击问题的分析中[1-4]。但由于推进剂供应管路系统是包括液路、气路和容腔的复杂管网系统，因此采用特征线法在时间步长协调上比较困难，边界条件的处理也比较复杂，而且非线性迭代解算收敛较慢[5]。常用求解水击问题的方法还有有限元法[6-8]、有限体积法[9]和无网格法[10]等。

Quarteroni,Canuto和郭本瑜等[11-12]提出并发展的谱方法 (Spectral Method)提供了除有限差分和有限元法之外的第三类数值计算方法,开辟了数值计算的全新领域。谱方法的最大优点是所谓“无穷阶收敛性”，即如果原问题的解充分光滑，那么用适当的谱方法所求得的近似解将以基函数个数的任意次幂速度收敛于精确解。本文采用谱方法中的Fourier谱方法就推进剂供应系统管路中阀门关闭形成的水击问题进行数值计算，给出了数学模型，求解原理、求解方法及仿真结果，并将该结果与已发表的采用特性线法和有限元法的求解结果进行了比较。

1 一维瞬变流动控制方程

假设管道内推进剂流动是一维绝热有摩擦的液体瞬变流，摩擦损失均按准稳态处理，则描述管道内液体瞬变流的基本方程包括连续方程和运动方程，是一组拟线性双曲型偏微分方程，具体形式如下[4]：

式中：p为x断面的流体平均压力；u为x断面的流体平均流速；ρ为流体密度；a为压力波传播速度；f为管壁摩擦损失系数；g为质量场加速度；θ为管流轴线方向与加速度方向的夹角。

式 (1)和 (2)即为管内非定常流的基本方程，又称波动方程。在一般流体管道中，有a＞＞u，故可略去式 (1)和 (2)中第二项，对于水平直管还可忽略重力场影响。在此基础上，令

那么，可得无量纲波动方程：

式中：p0为管道内液体初始压力。

2 Fourier谱方法求解

2.1 Fourier谱方法

谱方法源于经典的Ritz-Galerkin方法，以整体无限光滑的函数系 (例如：Fourier级数，Chebyshev多项式，Legendre多项式等)作为基函数，对定义在计算域中的函数进行近似，然后结合加权余量法对偏微分方程进行求解。推进剂供应系统管路中流体瞬变特性具有时间周期性，因此可通过选用特殊的周期函数 (Fourier级数)作为基函数，将未知函数分解展开成空间变量基与时间函数系数的分离变量形式的解式。这里将液体无量纲压力φ(y,τ)作Fourier正弦级数展开，液体无量纲流速U(y,τ)作Fourier余弦级数展开。

式 (5)是通过φy(y,τ)的Fourier余弦级数展开，然后逐项积分所得。式 (5)和 (6)中，与时间相关的系数 φk(τ)和 Uk(τ)为：

式中，k∈［1,2,… ］N 。这样，非线性无量纲控制方程 (3)和 (4)求解就转换为确定系数函数φk(τ)和Uk(τ)的形式。Fourier谱展开的项数N选取取决于计算精度。

将液体无量纲压力φ (y,τ)、流速U(y,τ)的Fourier谱展开式 (5)和 (6)代入方程 (3)和 (4)中，得：

值得注意的是，非线性摩擦积分项出现在方程(10)中，令

对摩擦项积分，W的Fourier余弦系数可表达成离散的Fourier变换形式

式中：yp＝p/K。

式 (9)和 (10)未考虑边界问题，在具体应用时，带入边界条件和初值，可得展开系数φk(t)，Uk(t)(k=0,1,2,…N)。再将其带入式 (5)和 (6)中，即可求得管路内液体瞬变流的近似解。为了适应谱方法的高精度并保证时程积分的稳定性，采用具有四阶精度的Runge-Kutta积分法。

2.2 非物理性振荡问题

应用谱方法求解管道内流体瞬变流动时，常会遇到带有大梯度甚至间断的问题 (如水击问题)。在它们附近数值解会出现大幅度的Gibbs类型的振荡。为了克服这种困难，学者们提出很多改进方法，其中包括超粘性方法、滤波法和谱粘性消去法 (SVV)等[13-14]。其中滤波法主要思想是用一个允许低频波通过的滤波因子使得导致数值不稳定的高频波大大减弱，那么谱方法在处理水击和其它间断处的振荡基本上被抑制。

文中采用滤波法消除因间断导致的非物理振荡，对无量纲压力和无量纲流量的展开系数φk(τ)和 Uk(τ)乘以滤波因子 σk。σk为

式中：系数α和η的取值参考文献 [14]。

3 算例及结果分析

为了验证Fourier谱方法计算水击问题的有效性，对文献 [15]中的单管阀门关断问题 (见图1)进行了仿真计算。问题描述：流体介质为水，原始参数：pT=15.0 MPa，pE=14.8 MPa，水平圆管长度L=4.8 m，内径D=0.01 m，ρ=1 000 kg/m3，α=1 200 m/s，CVAV=0.7 A，A为管道流通面积，摩擦系数f=0.018。

考虑阀门两种关闭方式：

(a)管道阀门突然关闭；

图1 简单管路系统Fig.1 A simple pipeline system with valve and tank

按Fourier谱方法 (SFM)的计算结果与按特征线方法 (MOC)[1]和分段集中参数有限单元法(FEM)[8]的计算结果比较见图2和图3。计算时，Fourier展开项数取8，为保证特征线方法和有限单元的精度，特征线方法分段数取为40，分段集中参数有限单元法取分段数为50。由图可见：采用Fourier谱方法计算单管水击问题无论是最大水击压力还是水击压力波的衰减速率都能与特征线方法和分段集中参数有限单元法保持一致。

滤波法的作用是避免因大梯度参数变化或间断问题引起计算时的非物理振荡。从图4可以看出，如果不采用滤波技术，在管内流体实际振荡曲线上叠加一个高频振荡，而采用滤波技术能较好的抑制这一高频振荡。

图4谱滤波后计算结果对比Fig.4 Contrast of calculated results after filtering by spectral-Fourier method

图5 给出了谱展开项数不同，利用无滤波的Fourier谱方法计算单管水击问题时，阀门上游压力变化。可以看出，当Fourier展开项数由8项增加到16和32项时，由数值计算引起的高频振荡未能得到有效抑制，因此采用Fourier谱方法计算大梯度或间断问题时，增加展开项数并不能有效抑制非物理振荡。

图5 Fourier展开项数对数值振荡的影响Fig.5 Influence of different Fourier terms on numerical oscillation

4 结论

推进剂供应系统管路水击问题的Fourier谱方法求解结果表明：

1）将Fourier谱方法应用于推进剂供应系统管路水击这一强非线性问题计算时，求得结果与特征线法和分段集中参数有限单元法结果一致，表明方法和仿真结果可信。

2）采用滤波技术能够较好地抑制数值解在间断处产生的非物理振荡。增加Fourier展开项数并不能有效抑制非物理振荡。

3）采用Fourier谱方法求解一维瞬变流动模型时，边界处理简单，可以方便地连接集中参数元件，适用于复杂流体网络的建模与仿真，工程实际应用有广阔前景。

[1]WYLIE E B,STREETER V L.Fluid transient in systems[M].Endlewood Cliffs,USA:Prentice Hall,1993.

[2]RUTH E K,AIM H,BAKER R L,et al.Advanced liquid rocket engine transient model,AIAA1990-2299[R].USA:AIAA,1990.

[3]SAGNICK H D,KRFILL G.Numerical simulation of transients in feed systems for cryogenic rocket engines,AIAA 95-2967[R].USA:AIAA,1995.

[4]林景松,王平阳,高红,等.液体火箭发动机关机水击的数值模拟[J].上海航天,2008,25(3):53-57.

[5]张育林,刘昆,程谋森.液体火箭发动机动力学理论与应用[M].北京:科学出版社,2005.

[6]KOLCIO K,HELMICKI A J,JAWEED S.Propusion system modeling for condition monitoring and control:partⅠtheoretical foundation,AIAA 94-3227[R].USA:AIAA,1994.

[7]KOLCIO K,HELMICKI A J,JAWEED S.Propusion system modeling for condition monitoring and control:partⅡapplication to the SSME AIAA 94-3228[R].USA:AIAA,1994.

[8]刘昆,张育林.一维可压缩流的有限元状态空间模型[J].推进技术,1999,19(5):62-66.

[9]MAJUMDAR A k.Numerical modeling of fluid transient by a finite volume procedure for rocket propulsion systems[C]//Proceedings of4th ASME/JSME JointFluids Engineering Conference.Honolulu,Hawaii,USA:ASME/JSME,2003:2967-2974.

[10]TIJSSELING A S,BERGANT A.Meshless computation of water hammer[C]//Proceedings of 2nd IAHR International Meeting of the Workgroup on Cavitation and Dynamic Problems in Hydraulic Machinery and Systems.Timisoara,Romania:IAHR,2007:1-34.

[11]CANUTO C,QUARTERONI A.Approximation results for orthogonal polynomials in Sobolev spaces[J].Math Compu,1982,38(1):67-86.

[12]GUO B Y.Spectral methods and their applications[M].Singapore:World Scientific,1998.

[13]MA H P.Chebyschev-Legendre super spectral viscosity method for nonlinear conservation laws[J].SIAM J.Numer.Anal.,1998,35(3):893-908.

[14]SHEN Jie,TANG Tao.Spectral and high-order methods with applications[M].Beijing:Science Press,2006.

[15]刘昆,张育林.推进剂供应管道的集中参数近似模型研究-模态近似模型[J].推进技术,1998,19(4):41-45.