0 引言

直接数值模拟是研究流动问题的有效方法。由于计算域总是小于实际流动区域, 选择恰当的边界条件是保证计算可靠的必要条件。对可压缩流, 较常用的是基于特征关系的边界条件, 例如见文献[1-4]。但在非定常流的计算中, 这种方法有时并不能给出令人满意的结果, 特别是要求无反射的出流边界条件。其主要原因是, 从特征关系角度看, 只要是亚声速区, 就总有入射特征波。除非能事先给出此波之值, 否则总会有误差。而事先给出其值一般来说是不可能的。

在实现无反射出流条件的方法中, 有一种所谓的嵌边法(fringe method), 也称缓冲区法、海绵区法等(见文献[5])。其主要思想是:在出流处增加一段特殊区域——嵌边区, 在该区域内, 经过修正的Navier-Stokes方程可将出流扰动衰减至足够小以尽可能地减弱出流反射波的强度。这一方法对不可压缩流的有效性已在实际应用中被证实, 本文将通过对若干可压缩流的试用, 看其在可压缩流中的有效性, 还研究了有关参数的选取问题。

1 嵌边法原理

在原计算域要应用此方法的边界外增加一段计算域, 在此域内对N-S方程增加一项特殊的源项。经过修正的守恒型N-S方程可以表示为:

$ \frac{\partial U}{\partial t}+\frac{\partial E}{\partial x}+\frac{\partial F}{\partial y}+\frac{\partial G}{\partial z}=\frac{\partial E_{v}}{\partial x}+\frac{\partial F_{v}}{\partial y}+\frac{\partial G_{v}}{\partial z}+H $

(1)

$ \begin{array}{*{20}{c}} {H = \lambda (x)\left( {{U_e} - U} \right)}\\ {U = \left( {\rho , {\rho _u}, {\rho _v}, {\rho _w}, \rho {e_s}} \right)}\\ {{e_s} = \frac{{{a^2}}}{{\gamma (\gamma - 1)}} + \frac{1}{2}\left( {{u^2} + {v^2} + {w^2}} \right)} \end{array} $

式中, U为守恒型通量, $\frac{\partial E}{\partial x}、\frac{\partial F}{\partial y} $和$\frac{\partial G}{\partial z}$为对流项, $\frac{\partial E_{v}}{\partial x}、\frac{\partial F_{v}}{\partial y}$和$\frac{\partial G_{v}}{\partial z}$为粘性项, H为嵌边法对N-S方程的修正项。H中, U是要求解的通量, U_e是预先指定的嵌边区出流处U应取的值。λ(x)在嵌边区内总为非负数。

嵌边函数λ(x)的选取要保证未知通量U在嵌边区内平缓地趋近U_e, 其解析表达式为:

$\lambda(x)=\lambda_{\max }\left[S\left(\frac{x-x_{\text {start }}}{\Delta_{\text {rise }}}, 1\right)-S\left(\frac{x-x_{\text {end }}}{\Delta_{\text {fall }}}+1, -1\right)\right]$

(2)

式中形函数S(x, ϕ)为:

$ S(x, \phi ) = \left\{ {\begin{array}{*{20}{c}} 0\\ {\frac{1}{2}\left[ {\tanh \left( {7x - \frac{{7 + \phi }}{2}} \right) + 1} \right]}\\ 1 \end{array}} \right\}\begin{array}{*{20}{c}} {x \le 0}\\ {0 < x < 1}\\ {x \ge 1} \end{array} $

(3)

λ_max, x_start, x_end, Δ_rise, Δ_fall的含义如图 1所示。

2 数值计算方法

本文采用的符号及方程均与文献[6]相同, 为节省篇幅, 不再重复。在对守恒型N-S方程进行数值模拟时将对流项分裂, 采用五阶精度的弱迎风紧致格式, 粘性项则采用六阶中心型紧致差分格式求两次一阶导数, 时间导数采用二阶Runge-Kutta方法的中点格式(见文献[7])。

3 算例 3.1 亚、超声速平板边界层扰动波的演化

当计算域较平板短时, 就有扰动波是否能顺利从计算域下游出口传出去的问题。

设来流马赫数为0.7, 以无穷远处的速度、密度、粘性系数及计算域入口处边界层位移厚度作为无量纲化参数, 则雷诺数为9000。有效计算域长100, 另加长为30的嵌边区。

首先要计算基本流。为此先求相似性解, 以作为初始流场。然后保持入口处不变, 求解N-S方程, 直至定常为止。

在计算域入口处引入T-S波。在下游的嵌边区内采用修正的N-S方程进行计算。计算足够多整周期, 在扰动波已穿过嵌边区后, 即停止计算。

为了考察嵌边法的有效性, 又计算了一个参考算例与之比较:取足够长的计算域, 使得计算若干整周期后, 扰动波已穿过原计算域出口(嵌边区入口), 但下游还有足够长的计算域, 以致对计算结果不会有影响。比较两个算例在原计算域(不含嵌边区)内两者的结果, 即可断定嵌边法是否行之有效。

在入口处所加T-S波的振幅为0.01, 频率为0.154。嵌边法算例中计算域(包含嵌边区)为130×20, 网格数为130×200, 法向采用变网格, 越靠近壁面网格越密。嵌边函数如图 2所示。有关参数为:L_f=30, λ_max=1.0, Δ_rise=5/16L_f, Δ_fall=3/16L_f, 嵌边区出口用预先指定的出流条件给定。

图 3给出距壁面1.85处扰动速度v的分布。图中reference表示参考算例结果, fringe表示嵌边法算例结果(其它图同)。可以看出, 二者在原计算域中吻合得很好, 而在嵌边法算例的嵌边区内, 扰动波有明显的衰减, 通量逐渐改变至预先指定的出流通量值。

我们试用文献[1]中的无反射出流条件计算该问题, 反射系数在0与0.3之间调整, 但总不能得到令人满意的结果。实际上, 该条件是利用出口压力与环境压力差来确定入射特征波值, 但在非定常情况中无法保证其值与真实情况一致。

还计算了超声速平板边界层中扰动波的演化问题。来流马赫数为4.5, 雷诺数为90000, 可以得到与亚音速相同的结论。把扰动幅值增大至0.05, 或把扰动波换成振幅、频率各不相同的2~3个T-S波, 都能得到令人满意的结果。见图 4和图 5 (图中30~36为嵌边区)。实际上, 如出口处全是超声速, 则不会有入射特征波(不考虑粘性时), 因而出口边界条件很容易处理。而超声速边界层近壁面处总是亚声速, 但所占区域很小, 其影响随着马赫数提高而减小。因而超声速边界层出口边界条件的处理比亚声速边界层要容易一些。

3.2 超声速流中涡的传播

文献[1]中有一个涡向下游传播穿过出流边界的算例。初始的涡如图 6所示。

图 7、图 8的涡量和流向速度的等值线图中负值用虚线表示, 正值用实线表示。流向速度u₁输出绘图时用(u₁-u₀) /u₀表示, 因此虚线对应于小于u₀的速度值, 实线对应于大于u₀的速度值。u₀是初始时刻入口处流向速度。

用嵌边法处理出流边界, 嵌边区长度L_f=30, 其他参数选取同例3.1。

图 7、图 8为嵌边法得到的结果。从涡量角度看, 可以认为已无反射通过出口。但从速度来看, 仍有小量反射。与文献[1]中结果比, 要稍好一些。

上面曾提到, 不可能事先知道出口处入射特征波之值是无法做到完全无反射的根本原因。因此比较出口处从计算结果反推的入射波和实际应有的入射波是鉴别方法好坏的更有效指标。而实际应有的入射波则可以通过加长计算域, 使得在同一计算时间内扰动尚未到达出口处, 再反推原出口处的负通量而得到。下面将从这一角度比较嵌边法和特征关系边界条件。

图 9中给出了ct/l=2时刻实际出流位置处三种算法, 即用加长计算域(参考算例)、嵌边法及文献[1]中方法(结果用NSCBC表示)压力反射系数取为0.01时算得的x方向对流项的负通量沿y的变化。可以看出, 在涡穿过边界的过程中, 嵌边法比用特征关系边界条件更接近于参考算例。改变特征关系边界条件的压力反射系数, 增大至0.0525, 结果变化不大, 其差别在10^-6之内。

3.3 亚声速平板边界层尾缘扰动的演化

在气动声学的研究中, 人们对边界层中的T-S波在平板(机翼)尾缘处的绕射问题感兴趣。若声波不能顺利地从计算域传出, 将影响计算结果。以下我们研究一无限薄的二维平板的上边界层中的T-S波脱离边界层后能否顺利传出计算域边界的问题, 看嵌边法是否能让扰动顺利传出。

研究问题模型如图 10所示。虚线所示的矩形区域为计算域, 坐标原点取在平板尾缘。为计算基本流场, 初始流场取为:平板上下均由层流边界层的相似性解给出, 而尾缘后的尾流速度剖面, 用式(4)来近似。

$u=u^{*}\left(y^{*}\right)\left(1-x / x_{l}\right)+\left(x / x_{l}\right)\left(1-a e^{-b y^{2}}\right)$

(4)

其中, u^* (y^*)为入口处由Blasius相似性解得到的层流边界层流向速度, x、x_l如图所示, 本文计算中a=0.7, b=0.5。

与例3.1相同, 先由初始流场出发, 保持入口处不变, 计算至流场定常为止, 得到基本流场。然后在平板上侧的基本流入口处引入T-S波, 分别用嵌边法和文献[1] (简称NSCBC)中的出流条件, 算出两个结果。还计算了一个参考算例以与之比较。

平板上下两侧来流相同, 马赫数为0.7, 雷诺数为9000。

引入的T-S波振幅为0.05, 频率为0.41。嵌边法算例中-27~30是有效计算域, 后加20个长度的嵌边区。

图 11是第99个无量纲时间距壁面4.04个无量纲长度处参考算例和嵌边区算例流向速度u的扰动演化。从图中可以看出, 在实际物理出流处, 即嵌边区的入口处(本例中为x=30)之前的计算域中, 两者的结果吻合得很好, 在嵌边区域内, 波有明显的衰减。从通量守恒的角度比较三个算例在实际物理出流处, 即x=30处x方向负通量沿y的变化, 见图 12。显然, 嵌边法比NSCBC处理出流边界更令人满意, 事实上, 嵌边法与参考算例的负通量在图中几乎没有区别。

4 嵌边法参数的选取

嵌边法的使用无疑会增加计算量。显然, 其长度越短越好, 但它受几个因素制约。

首先, 嵌边区内扰动波会逐渐衰减到零, 衰减的快慢与λ_max、Δ_rise有关。λ_max越大, Δ_rise越小, 衰减能力越强, 即可保证在嵌边区内(U-U_e)尽快地减小至零。其次, λ(x)的下降段Δ_fall不是必须的, 但上升段Δ_rise是不可少的, 而且必须平滑地趋于λ_max, 不能太快, 否则在该处会引起扰动反传入计算域。第三, 给定λ_max和Δ_rise, 需要足够长度的嵌边区使得(U-U_e)衰减至零, 但是继续加长嵌边区或改变Δ_fall的值不会给计算结果带来明显变化。第四, λ_max并非越大越好, Δ_rise也并非越小越好, 因为λ_max太大, Δ_rise太小都无法保证在嵌边区入口端的平滑过渡。因此参数的选择实际是各方面平衡的结果。为了证实上述论断, 下面给出一些试算结果。

对于超声速平板边界层, 计算三个嵌边法算例, 计算域为36×7, 其中30~36为嵌边区, 嵌边区函数分别为λ₁、λ₂和λ₃, 如图 13 (a) (d)中所示。λ₁代表快速上升的函数, λ₃代表快速下降的函数。(b)、(c)和(e)、(f)分别为它们和参考算例横向速度的演化, 及在实际出流位置处通量的比较。从中可以看出嵌边函数上升太快会引起较大误差, 而若(U-U_e)在嵌边区出口处已衰减至接近于零时嵌边函数下降快慢对结果影响不大。

5 结论

本文以嵌边法作为出口边界条件, 通过对若干可压缩流场中的扰动演化问题的试算, 证实了嵌边法在可压缩流计算中, 比一般基于特征波分析的方法, 可以取得更接近于真正无反射边界条件的结果。特别是在亚音速情况更是如此, 因为在亚音速情况, 实际上总存在入射波。还分析了嵌边函数的参数选取对于计算结果的影响。虽然增加嵌边区会增加计算工作量, 但一般的实际问题, 其计算域都比较大, 增加嵌边区所增加的计算量所占比重并不大。例如, 有实际意义的边界层转捩问题的计算, 其流向计算域至少要有几百个边界层厚度, 或不小于100个T-S波的波长。而通常取1~2个T-S波波长的嵌边区就可以得到令人满意的结果。