力学学报 2018 , 50 (3): 538-552https://doi.org/10.6052/0459-1879-17-424

流体力学

改进虚拟边界算法在超声速流动问题求解中的应用

张阳 ^* ^,,邹建锋 ^* ^,,郑耀

浙江大学航空航天学院, 杭州 310027

AN IMPROVED GHOST-CELL IMMERSED BOUNDARY METHOD FOR SOLVING SUPERSONIC FLOW PROBLEMS

Zhang Yang ^* ^,,Zou Jianfeng ^* ^,,Zheng Yao

College of Aeronautics and Astronautics , Zhejiang University , Hangzhou 310027, China

中图分类号:V211.3

文献标识码:A

通讯作者: 通讯作者:张阳, 博士研究生, 主要研究方向: 超声速流动不稳定性.E-mail: yangzhang@zju.edu.cn;邹建锋, 副教授, 主要研究方向: 流动不稳定性, 爆震燃烧技术.E-mail: zoujianfeng@zju.edu.cn 通讯作者:张阳, 博士研究生, 主要研究方向: 超声速流动不稳定性.E-mail: yangzhang@zju.edu.cn;邹建锋, 副教授, 主要研究方向: 流动不稳定性, 爆震燃烧技术.E-mail: zoujianfeng@zju.edu.cn

收稿日期:2017-12-25

接受日期:2018-03-15

网络出版日期:2018-06-10

基金资助: 国家自然科学基金资助项目(11372276, 11432013).

展开

摘要

提出了一种改进的虚拟单元浸没边界法, 并与一种高阶格式的有限差分算法相结合, 运用于求解超声速复杂几何绕流问题.该算法的核心思想在于在固体边界的内部和外部分别施加满足边界关系的作用点, 使得几何边界离散更加细化, 起到了壁面附近网格局部加密的作用.采用源空间内流体点作为反距离插值算法的重构点, 有效避免了插值点数目过少而与作用点相重合情况.通过对二维激波反射现象 (马赫数为 2.81) 和三维超声速球体绕流问题 (马赫数为 1.2) 的数值模拟, 与实验结果对比表明, 本文改进算法相对一般的虚拟边界法来说能显著提高数值精度, 减小计算误差.计算结果揭示了球体绕流中剪切层、压缩波系和尾迹的相互作用导致自由剪切层失稳的机理.剪切层厚度和湍流雷诺脉动经历了线性增长、大幅度震荡和小幅度波动三个阶段, 导致剪切层表面褶皱因子变化呈指数规律增长.其湍流结构表现出明显的各向异性, 具体在流向雷诺正应力在湍流脉动中占主导地位, 激波的压缩作用对不同方向雷诺正应力的影响存在空间迟滞效应.

关键词： 虚拟单元边界法 ; 剪切层 ; 超声速流动

Abstract

An improved ghost-cell immersed boundary method proposed in this paper, coupled with a high order finite difference solver, is applied to simulate the supersonic compressive flows around the complex obstacles.The main improvement of this algorithm is the treatment of the solid boundary that both ghost points inside the solid domain and forcing points inside the fluid domain due to the extension of the boundary are chosen to reconstruct the flow information considering the effect of solid wall on fluid.This brings refined boundary with discrete points and strengthens the wall conditions, which plays the role of local mesh refinement.The fluid points are limited in a certain source space as the interpolating points of the inverse distance algorithm, which effectively avoids the fact that the interpolating points are too few to possibly lead to coincide with the forcing points.Two problems of two dimensional shock reflection (Ma=2.81) and three dimensional flow around the smooth sphere (Ma=1.2) demonstrate the significant improvement of the numerical accuracy compared to the general ghost cell method. The results reveal the instability mechanism of the free shear layer as a result of the interaction between the shear layer, the compression wave system and the wake. The thickness and Reynolds fluctuation of the shear layer experience three regimes of linear growth, large amplitude oscillation and small amplitude fluctuation, resulting in an exponential growth of wrinkling factor.The turbulent structure near the shear layer shows obvious anisotropy because the streamwise Reynolds normal stress is dominant and a spatial hysteresis exists in the effect of the tail shock on Reynolds normal stresses in different directions.

Keywords： ghost cell immersed boundary method ; shear layer ; supersonic flow

PDF (19393KB)元数据多维度评价相关文章收藏文章

本文引用格式导出EndNote Ris Bibtex

张阳 , 邹建锋 , 郑耀 . 改进虚拟边界算法在超声速流动问题求解中的应用 [J]. 力学学报 , 2018, 50(3): 538-552https://doi.org/10.6052/0459-1879-17-424

Zhang Yang , Zou Jianfeng , Zheng Yao . AN IMPROVED GHOST-CELL IMMERSED BOUNDARY METHOD FOR SOLVING SUPERSONIC FLOW PROBLEMS [J]. Chinese Journal of Theoretical and Applied Mechanics , 2018, 50(3): 538-552https://doi.org/10.6052/0459-1879-17-424

引言

在计算流体动力学领域, 通常会遇到复杂几何绕流问题, 例如超声速飞行器气动力/热问题, 发动机内涡轮叶片流动换热研究, 燃烧室内旋流器稳流特性研究等. 面对任意复杂的几何外形, 高精度高效率的数值方法显得十分重要. 当前比较常用的处理方法是生成贴体的结构或非结构网格^[1,2,3], 网格的形状与几何边界变形保持一致, 这样很容易施加边界条件, 采用高阶计算格式, 并在边界附近进行局部网格加密, 可以获得高分辨率的流动特征. 但往往生成贴体网格的过程需要耗费很大的计算资源, 同时针对复杂的几何外形, 生成高质量的网格并不容易, 从而极大地影响了数值计算精度, 特别是当处理动边界问题时, 需要在每一迭代时刻重新生成新的网格体系, 并利用插值方法将上一时刻流场解插值到新的网格系统下, 这无疑带来了很大的计算量. 浸没边界法 (immersed boundary method, IBM) 作为一种可行的选择, 能较好地处理复杂几何绕流问题, 它的主要优势在于不需要生成复杂的贴体网格, 只需要在规则的笛卡尔网格体系下考虑浸没边界对流体的影响, 这样在处理复杂静止或运动的几何边界问题上具有较高的效率, 同时由于采用规则的正交网格, 相对非结构网格来说, 极大地节省了计算资源^[4], 同时在程序编写方面, 非常适用于并行处理, 这对于求解大规模高精度流场问题来说具有十分重要意义.

浸没边界法首次被 Peskin^[5]提出, 用于模拟笛卡尔网格下血液流过心脏瓣膜时的相互作用, 通过对动量方程施加附加力来代替瓣膜的存在, 并应用一种离散的 Dirac-delta 函数将拉格朗日作用力分散到欧拉网格点上. 后来 Goldstein 等^[6], Saiki 和 Biringen^[7]发展了 Peskin 的方法, 他们采用反馈力模型使得刚性静止壁面上流体速度为零; 但反馈力模型会带来虚假的数值震荡, 同时采用显式的时间步进格式, 严重限制了时间步长大小^[6,8]. 此外, Kolomenskiy 等^[9,10]采用惩罚模型对浸没边界进行处理, 但其惩罚常数依赖于所求解的具体问题, 影响了数值解的收敛性^[11]. 以上方法都是将附加力添加在连续 N-S 方程中, 与此不同的是 Mohd-Yusof^[12]首次提出了直接附加力 (direct forcing) 方法, 将边界条件直接施加在离散的控制方程中, 得到离散的附加力. 这样, 控制方程中没有引入额外的附加项, 因此不存在数值刚性问题^[11]. Fadlun 等^[13]将此方法与基于交错网格的有限差分方法相结合, 运用于求解三维复杂问题, 展示了此方法相比反馈力模型更加高效.

本文采用虚拟单元浸没边界法 (ghost-cell immersed boundary method) 对高速可压缩问题进行求解, 其基本思想与文献^[12]的方法相似, 将边界信息直接通过固体几何内虚拟点加以反映, 不需要在控制方程中加入附加项, 因而能直接与已有的流体动力学求解器相结合, 在处理复杂界面方面有着广泛应用. 虚拟单元边界法的关键在于如何重构虚拟点的流场信息, 使得几何边界信息的刻画更加准确和稳定. 一般来说, 虚拟点沿着边界法向镜像外插是最常见的处理方式^[14], 一方面诺依曼边界条件的施加更加方便, 另一方面避免了负值加权系数的产生^[15]. Tseng 和 Ferziger^[16]系统分析了一阶线性重构、二次多项式重构和反距离加权重构在处理不同边界条件上的优缺点, 虽然虚拟点高阶重构具有较高的数值精度, 但会带来数值不稳定问题^[17], 同时计算量也会增大. 除了重构格式外, 虚拟点层数和边界附近网格密度对计算精度和收敛性影响也较大. 对于边界附近的对流项处理, Ghias 等^[18]采用中心离散格式, 虚拟点的选取使其附近仅存在一个流体点. Chaudhuri 等^[19]采用 2 层虚拟单元, 使之与 5 阶加权基本无振荡 (weighted essentially non-oscillatory, WENO) 格式相匹配用于求解激波与障碍物之间相互作用. 从几何角度来看, 虚拟层数越多, 对应几何边界离散点越多, 有助于提高计算精度, 但虚拟点越远离边界, 重构解受边界条件影响越小. 通常情况下, 需要在边界附近进行网格加密, 使得边界附近的虚拟点增多, 从而边界条件的刻画越准确, 但无论网格的加密还是虚拟层数增加都会增加计算量.

此外, 超声速或高超声速绕流问题, 往往伴随着激波/障碍物、激波/边界层、激波/激波之间相互作用, 因此网格疏密、计算格式及其虚拟点的分布将对整个流场细节的模拟非常重要. 但目前虚拟单元边界法大多数侧重不可压缩问题的研究^[20,21], 对于高雷诺数可压缩黏性问题, 特别是高马赫数问题的运用很少. Qu 等^[17]改进了 Seo 和 Mittal^[22]提出的高阶虚拟点重构格式, 将其用于不同马赫数 (1.20~2.37)下超声速可压缩问题的研究, 并消除了该格式在处理移动边界时出现数值震荡问题. Haugen 和 Kragset^[23]将虚拟浸没边界法与一种高阶可压缩代码 (PENCIL) 相结合, 对横向流中颗粒撞击圆柱过程进行了直接数值模拟. Luo 等^[24,25]测试了虚拟单元边界法在处理可压缩流动传热问题上的优势, 其研究雷诺数和马赫数分别在 5000 和 0.5 以下. 由于虚拟单元浸没边界法在处理复杂几何问题上的优势, 再加上笛卡尔网格体系下非常容易构造高阶有限差分格式和大规模并行操作, 因而虚拟单元边界法在处理超声速或高超声速下可压缩复杂流场问题方面具有很大的潜力.

本文从固体几何边界细化的角度, 考虑到边界网格疏密对计算精度的影响, 对基于文献^[19]的虚拟单元边界法进行改进, 并与一种高阶精度的有限差分方法相结合, 运用于求解超声速可压缩流动问题. 本文组织结构如下: 在 2.1 节以三维空间激波/圆柱相互作用 (马赫数为 2.81) 为例对本文算法的精度进行对比验证, 在 2.2 节以三维空间超声速球体绕流 (马赫数为 1.2) 为例进一步说明本文算法对于高雷诺数下超声速湍流和激波结构的捕捉能力, 并对球体绕流的自由剪切层演化过程进行细致分析. 第 3 节对全文进行小结.

1 均匀数学模型及计算方法

1.1 三维非定常可压缩湍流控制方程

三维可压缩流动的基本方程是 Navier-Stokes (N-S) 方程, 本文采用直接数值模拟 (direct numerical simulation) 方法对其进行求解.在数值求解过程中通常将控制方程进行无量纲处理, 本文选取来流参数 $ρ_{\infty}, u_{\infty}, T_{\infty}, ρ_{\infty} u_{\infty}$ 分别对整个流场物理参量 $P, V = (u, v, w), T, p$ 进行无量纲化, 流场特征长度为 $L$ , 特征时间为 $L / u_{\infty}$ , 由此笛卡尔坐标系下三维非定常可压缩 N-S 方程^[26,27,28]如下

$\frac{\partial Q}{\partial t} + \frac{\partial F_{1}}{\partial x} + \frac{\partial F_{2}}{\partial y} + \frac{\partial F_{3}}{\partial z} - \frac{1}{Re} (\frac{\partial G_{1}}{\partial x} + \frac{\partial G_{2}}{\partial y} + \frac{\partial G_{3}}{\partial z}) = 0$ (1)

其中 $Q$ 是守恒变量, $F_{1}, F_{2}, F_{3}, G_{1}, G_{2}, G_{3}$ 分别为 $x, y, z$ 方向的无黏和黏性通量

${\pmb Q}=\begin{bmatrix}\rho\\ \rho u\\ \rho v\\ \rho w\\ E\\ \end{bmatrix}, \ {\pmb F}_{1}= \begin{bmatrix} \rho u\\ \rho u^{2}+p\\ \rho uv\\ \rho uw\\ (E+p)u\\ \end{bmatrix}, \ \ {\pmb F}_{2}= \begin{bmatrix} \rho v\\ \rho uv\\ \rho v^{2}+p\\ \rho vw\\ v(E+p)\\ \end{bmatrix} $

${\pmb F}_{3}= \begin{bmatrix} \rho w\\ \rho wu\\ \rho wv\\ \rho ww+p\\ w(E+p)\\ \end{bmatrix}, \ {\pmb G}_{1}= \begin{bmatrix} 0\\ \tau_{xx}\\ \tau_{xy}\\ \tau_{xz}\\ u\tau_{xx}+v\tau_{xy}+w\tau_{xz}+k\dfrac{\partial T}{\partial x}\\ \end{bmatrix} $

${\pmb G}_{2}= \begin{bmatrix} 0\\ \tau_{xy}\\ \tau_{yy}\\ \tau_{zy}\\ u\tau_{xy}+v\tau_{yy}+w\tau_{zy}+k\dfrac{\partial T}{\partial y}\\ \end{bmatrix} $

${\pmb G}_{3}= \begin{bmatrix} 0\\ \tau_{xz}\\ \tau_{yz}\\ \tau_{zz}\\ u\tau_{xz}+v\tau_{yz}+w\tau_{zz}+k\dfrac{\partial T}{\partial z}\\ \end{bmatrix}$

其中是 $E$ 单位质量流体总能量, 黏性应力项为

$\tau_{xx}=2\mu\dfrac{\partial u}{\partial x}-\dfrac{2}{3}\rm{\mu div}({\pmb V})$

$\tau_{yy}=2\mu\dfrac{\partial v}{\partial y}-\dfrac{2}{3}\rm{\mu div}({\pmb V})$

$\tau_{zz}=2\mu\dfrac{\partial w}{\partial z}-\dfrac{2}{3}\rm{\mu div}({\pmb V})$

$\tau_{xy}=\mu\left(\dfrac{\partial u}{\partial x}+\dfrac{\partial v}{\partial y}\right)$

$\tau_{yz}=\mu\left(\dfrac{\partial w}{\partial y}+\dfrac{\partial v}{\partial z}\right)$

$\tau_{zx}=\mu\left(\dfrac{\partial u}{\partial z}+\dfrac{\partial w}{\partial x}\right)$

流体动力学黏度系数 $μ$ 由 Sutherland 公式得出

$μ = \frac{T^{1.5} (1 + T_{s} / T_{\infty})}{T + T_{s} / T_{\infty}}$

其中, $T_{s} = 124 K$ 是 Sutherland 常数.热传导系数 $k = \frac{C_{p} μ}{\Pr}$ , 普朗特数 $\Pr = 0.72$ .

1.2 特征空间下无黏项流通矢量分裂 (flux vector splitting)

在实际的数值计算中, 考虑到流体黏性的影响, 往往需要在固体壁面附近进行局部网格加密.因而为了计算的方便, 需要将流动控制方程从笛卡尔坐标空间 $(x, y, z)$ 变换到一般的曲线坐标空间 $(ξ, η, ζ)$ , 得到

$\frac{\partial \tilde{Q}}{\partial t} + \frac{\partial {\tilde{F}}_{1}}{\partial ξ} + \frac{\partial {\tilde{F}}_{2}}{\partial η} + \frac{\partial {\tilde{F}}_{3}}{\partial ζ} - \frac{1}{Re} (\frac{\partial {\tilde{G}}_{1}}{\partial ξ} + \frac{\partial {\tilde{G}}_{2}}{\partial η} + \frac{\partial {\tilde{G}}_{3}}{\partial ζ}) = 0$ (2)

$J = |\frac{\partial (x, y, z)}{\partial (ξ, η, ζ)}| = |\begin{array}{l} x_{ξ} & x_{η} & x_{ζ} \\ y_{ξ} & y_{η} & y_{ζ} \\ z_{ξ} & z_{η} & z_{ζ} \end{array}|$ (3)

其中, $J$ 是坐标变换的雅可比系数矩阵.

根据扰动波传播的方向性, 对无黏通量 $F_{i} (i = 1,2, 3)$ 使用 Steger-Warming 分裂^[29], 将其分解成正通量 $F_{i}^{+}$ 和负通量 $F_{i}^{-}$ , 然后再分别对其空间导数进行差分离散, 这样保证了熵增的性质和格式的稳定性.

$A = \frac{\partial F_{i}}{\partial Q} (i = 1,2, 3)$ 为其系数雅可比矩阵, 其特征值为

$\left.\begin{array}{ll}\lambda_1=\lambda_2=\lambda_3=k_1u+k_2v+k_3w \\ \lambda_4=\lambda_1+ck,\lambda_5=\lambda_1-ck \end{array}\!\! \right \}$ (4)

其中, $c$ 是当地声速, $k = \sqrt[]{k_{1}^{2} + k_{2}^{2} + k_{3}^{2}}$ .将特征值分解成正和负两部分

$λ_{k} = λ_{k}^{+} + λ_{k}^{-} (k = 1,2, 3,4, 5)$ (5)

其中 $λ_{k}^{+} = (λ_{k} + \sqrt[]{λ_{k}^{2} + ε}) / 2, λ_{k}^{+} = (λ_{k} - \sqrt[]{λ_{k}^{2} + ε}) / 2, ε$ 是一个小量, 可以避免驻点和声速点附近数值振荡.将正、负特征值分别代入正、负流通矢量 $F_{i}^{+}, F_{i}^{-} (i = 1,2, 3)$ 的表达式, 即可得到

$F_{i} = F_{i}^{+} + F_{i}^{-}$ (6)

具体分解得到的正、负通量 $F_{i}^{+}$ 和 $F_{i}^{-}$ 表达式可以参考文献^[28].

1.3 时空离散处理

对于 N-S 方程中非线性无黏项, 本文采用 Jiang 和 Shu^[30]提出的 5 阶 WENO 有限差分格式, 使得在激波捕捉问题上具有较高的精度和数值稳定性, 其空间导数 $f_{j}^{'}$ 表示为

$f_{j}^{'} = \frac{F_{j + 1 / 2} - F_{j - 1 / 2}}{h}$ (7)

其中,h是网格间距,${\pmb F}_{j+1/2}=\sum\limits^{4}_{m=1}\omega_{m}{\pmb P}_{m}$

$P_{1} = (- 3 f_{j - 3} + 13 f_{j - 2} - 23 f_{j - 1} + 25 f_{j}) / 12$

$P_{2} = (f_{j - 2} - 5 f_{j - 1} + 13 f_{j} + 3 f_{j}) / 12$

$P_{3} = (- f_{j - 3} + 7 f_{j} + 7 f_{j + 1} - f_{j + 2}) / 12$

$P_{4} = (3 f_{j} + 13 f_{j + 1} - 5 f_{j + 2} + f_{j + 3}) / 12$

其中, $ω_{m}$ 是权值系数, 具体取值可以参考文献^[30].

对于黏性通量项 $G$ , 采用 6 阶精度的中心差分格式进行离散, 其 $x$ 方向的一阶导数表示为

${\frac{\partial G}{\partial x}|}_{j} = \frac{1}{60 Δ x} (- f_{j - 3} + 9 f_{j - 2} - 45 f_{j - 1} + 45 f_{i + 1} - 9 f_{j + 2} + f_{j + 3})$

对于时间项推进, 本文采用 3 阶 TVD 格式的 Runge-Kutta 方法

$\left.\begin{array}{ll}{\pmb L}=\dfrac{\partial {\pmb Q}}{\partial t} \\ {\pmb Q}^{(1)}={\pmb Q}^n+\Delta t{\pmb L}({\pmb Q}^n) \\ {\pmb Q}^{(2)}=\dfrac{3}{4}{\pmb Q}^n+\dfrac{1}{4}{\pmb Q}^{(1)}+\dfrac{1}{4}\Delta t{\pmb L}({\pmb Q}^{(1)}) \\ {\pmb Q}^{n+1}=\dfrac{1}{3}{\pmb Q}^n+\dfrac{2}{3}{\pmb Q}^{(2)}+\dfrac{2}{3}\Delta t{\pmb L}({\pmb Q}^{(2)}) \end{array}\!\! \right \}$ (8)

1.4 一种改进的虚拟单元浸没边界法

从求解N-S方程角度看, 为了在流体内部施加边界条件以区分流体与固体之间的界面, 需要在动量方程右端增加作用力项 $f$ , 代表固体对流体的作用.Mohd-Yusof^[11]最先提出了一种直接作用力方法, 将作用力项增加到离散的动量方程中

$\frac{u_{i}^{n + 1} - u_{i}^{n}}{Δ t} = RH S_{i} + f_{i}$ (9)

其中, $RH S_{i}$ 表示动量方程中 $i$ 方向的对流和黏性项, $u_{i}$ 是 $i$ 方向的速度分量, $Δ t$ 是时间步长.

假如在 $n + 1$ 时刻施加 $u_{i}^{n + 1} = u_{Ψ}$ , 其中 $u_{Ψ}$ 是所定义的边界条件, 因而得到作用力 $f_{i}$ 的表达式

$f_{i} = - RH S_{i} + \frac{u_{Ψ} - u_{i}^{n}}{Δ t}$ (10)

这样 $f_{i}$ 的作用使得在求解流场的每一时间步内均满足所定义的边界条件.这样处理的好处是所定义的速度边界条件直接施加在流、固界面上, 求解 N-S 方程不需计算新的增加项, 因而没有增加额外的 CPU 计算时间, 同时相比文献^[11,13], 处理方法也不受时间步长的限制.

通常对于复杂曲线几何来说, 网格节点位置并不一定正好落在几何表面, 这就需要对几何边界进行局部解重构.一种方法是采用虚拟单元, 在固体几何内部构造虚拟边界, 将物理边界上的信息通过内部流体节点, 外推至虚拟边界上, 从而保证流、固界面满足所需要的物理边界条件.下面以二维模型为例, 阐述虚拟单元浸没边界法的基本过程.

本文参考了文献^[19]的研究工作, 如图 1所示.在进行数值计算之前, 首先需要识别出虚拟边界, 并存储其节点上的有效信息.整个计算区域 $Ω$ 由固体区域 $Ω_{s olid}$ 和流体区域 $Ω_{f luid}$ 组成, 为了分离物理几何边界, 需要对所有网格节点进行标记

$\left.\begin{array}{ll}\Re_{\rm fluid}=\left\{\varphi_{i,j}=1,(x_i,y_i)\in\mathit{\Omega}_{\rm fluid}\right\} \\ \Re_{\rm solid}=\left\{\varphi_{i,j}=0,(x_i,y_i)\in\mathit{\Omega}_{\rm solid}\right\} \\ \quad\quad i=1,2,\cdots,N_{x},j=1,2,\cdots,N_{y} \end{array}\!\! \right \}$ (11)

由于本文用于超声速问题求解, 对流项和黏性项分别采用 5 阶 WENO 和 6 阶中心差分格式, 为了边界插值格式的需要, 虚拟边界选择 3 层虚拟点 (ghost point, GP), 虚拟点的标记可以通过如下判别

$\begin{equation}\Re_{GP}= \begin{cases} (x_{i},y_{j})\in\Re_{\rm solid},{\rm if}~\exists(x_{k},y_{l})\in\Re_{\rm fluid} \\ {\rm for}~k=i-3,\cdots,i+3,l=j-3,\cdots,j+3 \end{cases} \end{equation}$ (12)

虚拟节点上的流场信息 $ϕ_{G}$ 需要通过物理边界条件和流体点重构获得, 本文采用虚拟点 (GP) 相对物理边界的镜像点 (image point, IP), 通过 2 阶精度的线性逼近来确定任何 GP 点上的流场信息, 这种方法能有效避免产生负值重构加权系数

$ϕ_{G} = 2 ϕ_{B} - ϕ_{I}$ (13)

但在数值实践中发现, 当几何边界附近的网格不太密集的情况下, 采用虚拟点的镜像关系, 对于几何边界的细化来说精度往往达不到要求, 尽管可以通过增加虚拟点层数来进一步细化几何边界, 但越远离边界的虚拟点, 采用镜像关系的线性插值显然是不适合的.这就需要在几何边界附近进行局部网格加密, 但增加了求解计算量.为了解决这一问题, 本文在紧靠边界区域拓展了 1 层流体作用点 (forcing point, FP), 同样采用镜像关系的线性插值, 重构这一层流体作用点的解.如图 2对比了两种方法下边界曲线 AB 的细化情况.由于增加了 1 层流体作用点, AB 段离散增加至 14 段, 同时相较于增加虚拟点层数, 边界离散点 (boundary point, BP) 分布更趋于均匀

$ϕ_{F} = (ϕ_{B} + ϕ_{I}) / 2$ (14)

通常情况下, 镜像点 (IP) 并不与网格节点重合, 因而需要通过内部流体点对 IP 点物理信息进行插值重构.最常用的是线性插值, Luo 等^[24]根据 IP 点所在单元 4 个流体节点, 采用双线性插值求得 IP 点流场信息, Majumdar 等^[15]对比了双线性插值和二次曲线插值两种方法的区别, 发现两者的精度并没有明显差异.Chaudhuri 采用 Franke^[31]提出的一种反距离加权插值方法对镜像点流场信息进行重构, 保证了重构解在极值点附近的光滑性

$\left.\begin{array}{ll}\phi_{\rm I}=\displaystyle\sum\limits^4_{m=1}\omega_m\phi_m \\ \omega_m=\eta_m\left(\displaystyle\sum\limits^4_{k=1}\eta_k\right) \\ \eta_m=1/d^2_k \end{array}\!\! \right \}$ (15)

其中 $ϕ_{m} (m = 1,2, 3,4)$ 表示包含 IP 点单元的 4个相邻节点 NP 上的流场信息, $d_{k}$ 是第 $k$ 个 NP 节点到 IP 点的距离, 如图 1所示.

显示原图|下载原图ZIP|生成PPT

图 1 虚拟点镜像处理

Fig. 1 Immersed boundary treatment

由于扩展 1 层流体作用点后, 可能会出现很多插值节点 NP 与流体作用点 FP 重合的情况, 导致插值点 NP 数目减少到 3 个. 为了解决这一问题, 本文将 4 点反距离加权插值方法扩展到限定一定范围内多点插值, 如图 3所示, 以几何边界外的镜像点 IP 为源点, 限定半径为R的圆形区域内所有的流体点 (NP) 作为插值点.注意半径R选取越大, 则插值计算量越大

$\begin{equation}\phi_{\rm I}=\sum^{np}_{m=1}\omega_m\phi_m \end{equation}$(16)

其中 $np$ 是圆形区域内的 NP 点数.

显示原图|下载原图ZIP|生成PPT

图 2 浸没几何边界的细化

Fig. 2 Immersed boundary subdivision

对于上面提出的改进虚拟单元浸没边界法很容易扩展至 3 维空间, 只需将插值点限定在以 IP 点为球心,R为半径的球体内. 此外, 由于整个计算区域采用笛卡尔坐标系下的正交网格, 使得在信息传递接口并行编程计算方面比较容易实现, 一方面可以应用于大规模并行数值计算, 研究复杂的流动现象及其物理机理;另一方面, 由于网格系统并不随着内部几何模型的改变而需要重新划分, 这使得虚拟单元浸没边界法非常适用于优化设计.

显示原图|下载原图ZIP|生成PPT

图3 改进的浸没边界法

Fig. 3 Modified immersed boundary treatment

2 结果与讨论

2.1 激波反射现象

在超声速流动问题中, 激波反射现象在过去一直被大家广泛研究^{[32,33,34,35]}. 当一道入射激波在传播过程中遇到一个圆柱障碍物时, 首先会在圆柱表面发生常规激波反射, 形成交汇的二激波结构, 随着交汇点在圆柱表面逐渐向后移动, 常规反射会转捩成马赫反射, 形成三分叉点的三激波结构. 下面以二维空间入射激波与圆柱相互作用过程为例, 对本文的虚拟单元浸没边界算法进行验证. 平面入射激波马赫数为 2.81, 整个计算域长宽尺寸分别为 $L_{x} \times L_{y} = 240 m m \times 60 m m$ , 圆柱直径 $D$ 为 10 mm. 进、出口分别是超声速进、出口边界条件, 上、下壁面均为无滑移的绝热壁面. 整个计算区域采用均匀的正交网格, 采用 5 种不同的网格尺寸算例来检验网格的无关性, 如表 1所示.

表 1 不同疏密的网格尺寸

Table 1 Calculation mesh size in different cases

Case	$N_{x} \times N_{y}$	$Δ x = Δ y$
Mesh1	500×125	0.048D
Mesh2	667×167	0.036D
Mesh3	1000×250	0.024D
Mesh4	2000×500	0.012D
Mesh5	4000×1000	0.006D

新窗口打开

图 4是不同时刻密度等值面云图, 显示了激波绕射圆柱的演化过程, 包括马赫波的形成 (a)、马赫波沿着壁面传播 (b)、马赫波在尾迹处碰撞 (c) 和二次 3- 波点的出现 (d). 整个流场结构能清晰的识别出入射激波 (I.S.)、反射弓形激波 (R.S.)、马赫波 (M.S.1 和 M.S.2)、三波点 (T.P.1 和 T.P.2)、滑移线 (S.L.1 和 S.L.2) 等复杂波系结构, 同时还能观察到在滑移线后流场卷吸形成小涡结构 (V.). 特别是在圆柱后尾迹处, 存在上、下两道马赫激波相互碰撞, 并在尾迹中发展形成二次三波系结构. 与 Bryson 和 Gross^[36]实验可视化图比较, 计算结果较好的捕捉到了流场的特征.

显示原图|下载原图ZIP|生成PPT

图 4 激波反射流场纹影图对比 ((a)和(b): 数值摸拟结果, (c)和(d): 实验结果^{[

36
]})

Fig. 4 Schlieren pictures for the flow evolution of the shock reflection ((a)and(b): numerical simulation, (c)and(d): experimental results^{[

36
]})

图 5展示了整个激波绕射过程中圆柱表面 4 个位置的压力系数 $2 P / ρ_{0} u_{0}^{2}$ 变化曲线 (其中 $ρ_{0}$ 和 $u_{0}$ 分别是超声速进口密度和速度), 1 是圆柱前表面滞止点位置, 2, 3, 4 分别对应着顺时针方向 $6 0^{\circ}, 12 0^{\circ}, 18 0^{\circ}$ 位置. 从中看出, 随着网格尺寸加密, 数值解逐渐趋于收敛, 网格独立性在 Mesh4 和 Mesh5 达到流场解析要求. 特别是在 Mesh1 空间分辨率下, 前滞止点 1 位置出现了较小的数值震荡. 在圆柱后流场, 由于存在激波与尾迹之间、激波与激波之间等相互作用过程, 流场波系结构变化较为复杂, 这使得低的空间分辨率根本无法捕捉到详细的流场特征, 如图 6所示的两种网格尺度下流向速度等值图, 此外在圆柱壁面 3, 4 位置, 压力系数的波动也出现了较大的偏差.

显示原图|下载原图ZIP|生成PPT

图 5 圆柱表面压力系数变化

Fig. 5 Pressure coefficients on the cylindrical surface in different spatial resolutions

前文提到, 本文对虚拟单元浸没边界法的改进之处在于边界的细化, 通过考虑同时将边界条件关系施加到固体区域内部的虚拟点和固体区域附近的加力点, 可以对边界附近网格起到类似局部加密作用, 下面采用 3 种不同网格尺度对该方法的局部加密效果进行测试.

表 2展示了 6 种不同的计算工况, 其中虚拟点表示仅采用固体几何内虚拟点施加边界条件, 虚拟 + 加力点表示同时采用边界内外的虚拟和加力点施加边界条件(简称Ref.).

表 2 不同工况下网格尺寸和边界处理方法

Table 2 Calculation conditions with different implementation methods

Case	$N_{x} \times N_{y}$	Ghost or forcing points
C1	667×167	ghost+forcing points
C2	1000×250	ghost+forcing points
C3	2000×500	ghost+forcing points
C4	667×167	ghost points
C5	1000×250	ghost points
C6	2000×500	ghost points

新窗口打开

这里需要强调的是选用了 Mesh2, Mesh3, Mesh4 这 3 种中间尺度的网格系统用于测试, 是考虑到较大的网格尺度会造成流场解出现较大偏差, 边界细化所起的作用不大. 此外, 虽然得不到该算例的精确解, 但可以根据网格独立性检验, 选择 C3 工况作为基准参考解(简称Ref.).

图 7是两种虚拟单元边界法壁面压力系数的数据对比, 可以看出对于改进的虚拟单元方法, 其相对数值误差比原始的虚拟单元方法小. 对于网格数目最稀疏的 C1 和 C4 两种情况, 虽然圆柱前缘滞止点的压力系数与基准解相差不大, 但在圆柱后缘 3, 4 点位置, 两种方法均与基准解有一定的差距. 在 C2 和 C5 情况时, 改进的虚拟单元方法其计算结果与基准解大致吻合. 这说明施加基于加力点的边界条件不仅能细化几何边界起到强化边界信息的作用, 而且改善了数值解的计算精度.但当网格数目本身很密情况下, 这种效果不太明显.图 8给出 C2, C3, C5 三种情况下三波点（T.P.1 和 T.P.2）的轨迹图, 从中可以看出基准工况的计算结果与文献^[36,37]的实验结果吻合较好, 证明了本文采用的基于改进虚拟单元边界法的高阶有限差分求解器能准确的模拟激波与障碍物之间相互作用.

显示原图|下载原图ZIP|生成PPT

图 6 流场马赫数等值图对比

Fig. 6 Comparison of Mach number contours

显示原图|下载原图ZIP|生成PPT

图 7 改进虚拟边界法的圆柱表面压力系数对比

Fig. 7 Comparison of pressure coefficients with modified ghost-cell methods

显示原图|下载原图ZIP|生成PPT

图 8 三波点轨迹图

Fig. 8 Comparison of triple point trajectory

此外,图 9给出了 $t = 22.9 s$ 时刻 C3 工况下流场的流线图, 左边是对固体内部除去虚拟点外所有固体节点速度场归零处理, 右边是对固体内部不作处理. 从中对比可以看出, 两种处理方法对外部流场的计算结果是一致的, 说明了本文采用虚拟单元方法无需对固体内部进行处理, 也验证了文献^[11]的说法.

显示原图|下载原图ZIP|生成PPT

图 9 流线图对比 (上: 内部特殊处理; 下: 内部不做处理)

Fig. 9 Comparison of streamlines (top: internal treatment of body; bottom: no treatment)

2.2 超声速三维球体绕流

2.2.1 数值验证与计算精度分析

为了进一步验证本文所提出算法的稳定性, 对三维超声速圆球绕流算例进行研究. 过去几十年, 许多学者致力于亚临界和超临界雷诺数下圆球绕流流动分离现象的研究^{[38,39,40,41]}, 特别关注于剪切层不稳定发展和尾迹涡旋结构的演化规律探索^{[42,43,44,45,46]}. 本文超声速来流马赫数马赫数为 1.2, 雷诺数 $Re = 1.62 \times 1 0^{5}$ , 圆球直径D= 10 mm, 计算几何模型长度为 $L_{x} \times L_{y} \times L_{z} = 24 D \times 6 D \times 6 D$ .为减少网格计算量, 计算域采用分块加密网格, 球体附近区域的网格尺寸 $Δ δ$ 最小, 以提高边界流场的分辨率, 如图 10所示. 根据局部加密尺寸 $Δ δ$ 不同, 比较 8 个计算工况下改进虚拟单元边界法的计算精度.

显示原图|下载原图ZIP|生成PPT

图 10 光滑球体绕流计算域示意图

Fig. 10 Computational domain of the smooth sphere

图 11展示了 S4 工况下瞬态流场结构数值纹影图, 在球体正前方形成一道脱体激波, 自由剪切层从球体表面分离后向下游发展, 并与球体后尾迹结构相互作用, 导致剪切层内涡旋结构向外卷起或向内收缩, 在外部超声速流动作用下形成一系列向外传播的压缩波系 (A处), 压缩波系叠加产生一道强的尾激波结构 (B处). 在尾激波后剪切层出现了间断结构 (C处), 交替间断的结构逐渐演化成大尺度的涡结构. 在尾激波前一段距离内, 剪切层最先开始变形扭曲, 出现失稳现象. 这是因为一系列弱的压缩波向上游传播对剪切层进行了预压缩作用. 变形的剪切层、尾迹与激波之间的相互作用更加剧了剪切层结构的失稳, 导致流场结构的演化更加复杂.

显示原图|下载原图ZIP|生成PPT

图 11 中截面上瞬态流场密度梯度数值分布云图

Fig. 11 Transient density contours in the middle section

表 3 不同条件的计算工况

Table 3 Eight cases with different implementation methods

Case	$N_{x} \times N_{y} \times N_{z}$	$Δ δ$	Method
S1	409×109×109	0.048D	ghost+forcing points
S2	546×125×125	0.036D	ghost+forcing points
S3	606×156×156	0.024D	ghost+forcing points
S4	694×194×194	0.012D	ghost+forcing points
S5	409×109×109	0.048D	ghost points
S6	546×125×125	0.036D	ghost points
S7	606×156×156	0.024D	ghost points
S8	694×156×300	0.012D	ghost points

新窗口打开

图 12是球体中截面的壁面平均压力系数和摩擦系数与实验结果^[47]的对比, 可以看出, 在球体迎风侧, 各工况计算结果与实验值相差不大. 在背风侧, 随着球体附近网格依次加密, 数值解的误差逐渐减小, 说明流场解的精度随着网格分辨率的增大而提高. 当球体附近网格尺度达到 0.012D 即 S4 工况, 数值解与实验结果吻合较好, 因此可以将 S4 结果作为基准解, 对 S1~S8 工况整个计算域内的计算误差进行对比分析, 定义流场解误差的L2范数

$\varepsilon_{2}=\left[\dfrac{1}{N_{x}N_{y}N_{z}}\sum_{i=1}^{N_{x}}\sum_{j=1}^{N_{y}}\sum_{k=1}^{N_{z}}(\psi_{i,j,k}^{N_{x}\times N_{y}\times N_{z}}-\psi_{i,j,k}^{694\times 194\times 194})^2\right]^{1/2}$

显示原图|下载原图ZIP|生成PPT

图 12 网格独立性检查

Fig. 12 Grid independence check

其中 $ψ_{i, j, k}^{694 \times 194 \times 194}$ 表示 S4 工况 $(i, j, k)$ 网格节点处流场物理量. 选取流场的流向速度 $u$ 和球体表面压力系数 $C_{p}$ 两个物理量作误差分析, 如图 13所示. 在球体边界附近流场解达到了二阶精度. 对于改进虚拟单元方法 (S1, S2, S3), 由于同时对虚拟和加力点施加边界关系, 使得界面信息更加细化, 相对一般的虚拟单元方法 (S5, S6, S7) 来说, 计算精度得到了提高. 本文算法分别采用 5 阶 WENO 格式和 6 阶中心差分格式对空间进行离散, 因而在流场光滑区域, 数值解可以达到高阶精度, 但在激波和球体边界附近, 数值解的精度会降低. $L_{2} - u$ 表示所有三维网格点流向速度u的全局误差, 其精度在二阶以上. 这说明本文的算法能较准确的模拟高速可压缩问题, 在全局范围内能达到高阶计算精度.

显示原图|下载原图ZIP|生成PPT

图 13 流向速度u和球体表面压力系数Cp的误差L2范数

Fig. 13 Global errors inuandCpversus number of mesh points

2.2.2 剪切层湍流结构定性分析

为了对整个流场结构的演化过程进行细致分析, 特别关注自由剪切层失稳特性的观察,图 14给出了不同时刻下球体绕流后等马赫线分布. 根据流场马赫数特征, 可以把整个流场分为 4 个部分, 即脱体激波与声速线 (马赫数为1) 之间的亚声速区域, 剪切层外声速线与尾激波之间的超声速区域, 剪切层内球体后方的尾迹结构, 尾激波后流场尾迹涡旋结构.从图 14可以看出, 在尾激波前后区域, 剪切层结构出现较为明显的扭曲变形, 这种扭曲扰动逐渐向上游剪切层传播, 形成小的波浪状剪切层结构, 这是由于压缩波向上游传播的缘故. 尾激波后的剪切层结构出现了较大的弯曲变形, 并呈现周期性变化. 剪切层在发生扭曲变形的同时, 逐渐向内收缩, 将尾迹流场分成 2 部分, 一部分是球体正后方剪切层内未受扰动的尾迹结构;另一部分是向下游发展的尾迹结构, 这部分结构是上游尾迹、变形剪切层和压缩波系共同作用的结果, 且随着流场的演化, 这部分尾迹结构逐渐与上游剪切层脱落, 形成向下游发展的不同尺度脱落涡旋结构. 由此构成了一个完整的剪切层内涡旋脱落演化过程.为了进一步刻画剪切层结构的变化, 选取 $t = 826.4 s$ 时刻流向不同位置横截面的周向速度分布进行分析.

显示原图|下载原图ZIP|生成PPT

图 14 不同时刻流场马赫数分布等值图

Fig. 14 Mach number contours with different time

前面已经指出, 剪切层附近的压缩波系向上游传播, 使得剪切层出现小波浪形的弯曲, 触发了剪切层流动的不稳定, 随后小变形剪切层、尾激波和尾迹之间的相互作用加剧了剪切层的弯曲变形, 直至最终失稳破碎成间断的涡结构.图 15给出了球体下游不同位置的周向速度梯度分布云图, 可以看出剪切层在最初阶段为完整的环形, 但随着流场向下游发展, 局部出现个别卷起的涡 (如图 $x / D = 5.25$ 和 5.5);在 $x / D = 6.0$ 处由于受到尾激波的强烈压缩作用, 剪切层出现多个弯曲褶皱, 形成卷起的花瓣结构; 花瓣结构受到激波和尾迹共同作用, 由中心对称结构逐渐演化成不规则的带状结构; 尾激波后剪切层不再保持连续的带状结构, 而是分散成多段, 形成多个卷起的涡旋结构.图 16展示了三维流场涡结构的演化情况, 在尾激波前, 等值涡管表面出现了多个波浪状的凸起, 这时涡管的变形还不是很明显.但经过尾激波后, 涡管形态出现严重的扭曲和缠绕, 并在流场最前端形成脱落的涡旋结构, 向前逐步推进.

显示原图|下载原图ZIP|生成PPT

图 15 流向不同位置横截面的数值云图 (物理量： $∂ u / ∂ r$ )

Fig.15 Numerical schlieren pictures in different streamwise sections ( $∂ u / ∂ r$ )

显示原图|下载原图ZIP|生成PPT

图 16 三维流场的涡结构演化

Fig.16 Three-dimensional vortex structures

2.2.3 剪切层湍流结构定量分析

为了对剪切层的变形失稳过程进行定量分析, 考察剪切层厚度和褶皱因子的变化规律. 剪切层厚度 $δ_{t}$ 的定义为

$δ_{t} (x) = \frac{u_{1} - u_{2}}{(\partial u / \partial r)_{m ax}}$

其中 $u_{1}$ , $u_{2}$ 分别表示剪切层内外高速和低速区域流向速度.

图 17所示是剪切层厚度沿流向的变化情况, 可以看出剪切层厚度的变化呈现 3 个阶段：第 I 阶段为厚度的线性增长, 这是由于剪切层受内部尾迹回流影响, 再加上外面受到弱的压缩波扰动作用所导致的;第 II 阶段剪切层厚度呈现无规律的波动状态, 这时剪切层受到尾激波和尾迹的共同作用, 出现复杂的弯曲、缠绕和变形; 第 III 阶段剪切层厚度变化不大, 对应着尾迹涡旋向后逐步发展的过程. 需要注意的是, 虽然在剪切层的失稳过程中剪切层厚度变化处于上下波动状态, 但其变形扭曲却呈规律性的增长, 定义剪切层的褶皱因子 $Ω_{w}$ 为

$Ω_{w} = \frac{DL (x)}{S (x)}$

显示原图|下载原图ZIP|生成PPT

图 17 剪切层厚度的变化

Fig. 17 Change of the shear layer thickness

这里考虑流向不同横截面 $x = x_{0}$ 下二维平面褶皱因子,L和S分别代表环状剪切层的周长和面积.图 18是 2 个时刻流向剪切层褶皱因子的变化情况, $x / D$ 在6.0~7.5范围内, 剪切层由于受尾激波的作用, 虽然其褶皱因子呈现小幅度的波动, 但总体上以指数递增规律变化.通过对离散点的数据拟合, 得出了剪切层褶皱因子增长的曲线规律

显示原图|下载原图ZIP|生成PPT

图 18 剪切层褶皱因子的变化

Fig. 18 Wrinkling factor of the shear layer

$Ω_{w} (\frac{x}{D}) = 0.45 e^{0.39 \frac{x}{D}}$

褶皱因子的不断增长代表着流场剪切层与内部低速流体、外部高速流体的接触面积增大, 高低速流体的动量交换更加剧烈. 与此同时, 剪切层表面曲率变化越大, 容易产生向外卷起的涡结构, 在外部超声速流动下形成向外传播的压缩波, 进而对附近的剪切层产生干扰作用. 在 $x / D > 8$ 的尾迹区域, 由于远离尾激波的影响, 脱落的涡旋结构在向前推进过程中变化很小, 因而其褶皱因子基本保持不变.

图 19为剪切层速度分布云图.图 20和图 21给出了剪切层内湍流雷诺应力的演化规律, 剪切层的位置取为图 19中虚线, 在 D 处剪切层开始出现波浪状变形褶皱, S 处为尾激波所在位置.从湍流速度均方根分布可以看出, 3个方向的雷诺正应力呈现大致相同的特征趋势走向, 与剪切层厚度的演化规律一致, 第 I 阶段为雷诺正应力的线性增长阶段, 自由剪切层为稳定的平直状态; 第 II 阶段由于压缩波系的作用使得湍流结构更为复杂, 湍流脉动强度经历一个混沌的震荡过程, 并且在尾激波位置流向脉动强度达到最大值; 第 III 阶段为尾激波后涡旋的脱落向前推进过程. 还可以观察到, 超声速球体绕流的湍流雷诺正应力表现出明显的各项异性, 雷诺应力中流向正应力最大, 这说明在剪切层的演化过程中湍流结构的流向脉动占主导地位.此外,z方向的正应力最大值在尾激波后一段距离内取得, 说明了激波的压缩作用对不同方向雷诺正应力的影响存在空间迟滞效应, 进一步阐述了剪切层内湍流结构的各项异性.从湍流雷诺剪切应力分布可以看出, 最大雷诺剪切应力仍然出现在尾激波所在位置 (图中S), 不存在空间上的滞后, 但雷诺剪切应力的发展在激波附近存在一个先增加后减小的突变过程, 这一点与雷诺正应力的演化规律不同. 与此同时, 在 $x / D = 5.5$ 处 (图中D) 出现了雷诺剪切应力出现了局部最大值, 这是因为剪切层附近压缩波向前传播, 对 D 位置的剪切层产生了压缩作用, 正是第 II 阶段剪切层失稳过程的开始.

显示原图|下载原图ZIP|生成PPT

图 19 剪切层速度分布云图

Fig. 19 Streamwise velocity distribution of the shear layer

显示原图|下载原图ZIP|生成PPT

图 20 剪切层雷诺正应力分布

Fig. 20 Reynolds normal stress distribution of the shear layer

显示原图|下载原图ZIP|生成PPT

图 21 剪切层主切应力分布

Fig. 21 Primary shear stress distribution of the shear layer

3 结论

本文提出了一种改进的虚拟单元浸没边界法, 并与一种高阶格式的有限差分算法相结合, 运用于求解超声速复杂几何绕流问题. 通过对二维激波反射问题和三维超声速球体绕流问题的数值分析, 可以得到以下结论：

(1) 该改进算法的核心思想在于通过在固体边界的内部和外部同时施加作用点 (ghost point和forcing point), 使得几何边界离散更加细化, 强化了壁面边界关系, 起到了壁面附近局部网格加密的作用. 与此同时, 采用限定半径为R的源空间内流体点作为反距离插值算法的重构点, 有效地避免了插值点数目过少以至于与流体作用点相重合情况.

(2) 在相同网格计算量下, 本文提出算法相对一般的虚拟边界法可以达到提高计算精度, 减小数值误差的目的. 反距离插值算法对浸没边界重构点的处理, 使得流场在边界附近能达到 2 阶精度, 5 阶 WENO 格式和 6 阶中心差分格式对空间项的处理, 使得整个流场能达到高阶计算精度, 能较准确的捕捉到激波与固体壁面相互作用的流场特征, 刻画障碍物后湍流结构的流动细节.

(3) 揭示了三维球体超声速绕流自由剪切层失稳的机理, 剪切层附近的压缩波系向上游传播, 使得剪切层出现小波浪形的弯曲, 触发了剪切层流动的不稳定, 随后小变形剪切层、尾激波和尾迹之间的相互作用加剧了剪切层的弯曲变形, 直至最终失稳破碎成间断的涡结构.

(4) 超声速自由剪切层湍流特征的发展经历了 3 个阶段：第 I 阶段是剪切层厚度和湍流雷诺脉动线性增长过程; 第 II 阶段是剪切层失稳的大幅度震荡过程, 剪切层表面褶皱因子变化呈现指数规律增长; 第 III 阶段是剪切层破碎成涡旋结构导致小幅度的波动过程. 此外, 自由发展剪切层内湍流结构具有明显的各项异性特征, 具体表现在流向雷诺正应力在湍流脉动中占主导地位, 激波的压缩作用对不同方向雷诺正应力的影响存在空间迟滞效应.

The authors have declared that no competing interests exist.

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]	Zheng Y , Liou MS . A novel approach of three-dimensional hybrid grid methodology: Part 1. Grid generation . Computer Methods in Applied Mechanics and Engineering, 2003 , 192 ( 37 ): 4147 - 4171 [本文引用: 1]
[2]	Zhao D , Chen J , Zheng Y , et al . Fine-grained parallel algorithm for unstructured surface mesh generation . Computers& Structures, 2015 , 154 : 177 - 191 [本文引用: 1]
[3]	Kadoch B , Reimann T , Schneider K , et al . Comparison of a spectral method with volume penalization and a finite volume method with body fitted grids for turbulent flows . Computers& Fluids, 2016 , 133 : 140 - 150 [本文引用: 1]
[4]	Sotiropoulos F , Yang X . Immersed boundary methods for simulating fluid-structure interaction . Progress in Aerospace Sciences, 2014 , 65 : 1 - 21 [本文引用: 1]
[5]	Peskin CS . Flow patterns around heart valves: A numerical method . Journal of Computational Physics, 1972 , 10 ( 2 ): 252 - 271 [本文引用: 1]
[6]	Goldstein D , Handler R , Sirovich L . Modeling a no-slip flow boundary with an external force field . Journal of Computational Physics, 1993 , 105 ( 2 ): 354 - 366 [本文引用: 2]
[7]	Saiki E M , Biringen S . Numerical simulation of a cylinder in uniform flow: Application of a virtual boundary method . Journal of Computational Physics, 1996 , 123 ( 2 ): 450 - 465 [本文引用: 1]
[8]	Lee C . Stability characteristics of the virtual boundary method in three-dimensional applications . Journal of Computational Physics, 2003 , 184 ( 2 ): 559 - 591 [本文引用: 1]
[9]	Kolomenskiy D , Schneider K . A Fourier spectral method for the Navier-Stokes equations with volume penalization for moving solid obstacles . Journal of Computational Physics, 2009 , 228 ( 16 ): 5687 - 5709 [本文引用: 1]
[10]	Kadoch B , Kolomenskiy D , Angot P , et al . A volume penalization method for incompressible flows and scalar advection-diffusion with moving obstacles . Journal of Computational Physics, 2012 , 231 ( 12 ): 4365 - 4383 [本文引用: 1]
[11]	Iaccarino G , Verzicco R . Immersed boundary technique for turbulent flow simulations . Applied Mechanics Reviews, 2003 , 56 ( 3 ): 331 - 347 [本文引用: 5]
[12]	Mohd-Yosuf J . Combined immersed boundary/B-spline methods for simulations of flow in complex geometries, Annual Research Briefs, Center for Turbulence Research, NASA Ames Research Center/Stanford University 1997 , 317 - 328 [本文引用: 2]
[13]	Fadlun EA , Verzicco R , Orlandi P , et al . Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations . Journal of Computational Physics, 2000 , 161 ( 1 ): 35 - 60 [本文引用: 2]
[14]	Mittal R , Dong H , Bozkurttas M , et al . A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries . Journal of Computational Physics, 2008 , 227 ( 10 ): 4825 - 4852 [本文引用: 1]
[15]	Majumdar S , Iaccarino G , Durbin P , RANS solvers with adaptive structured boundary non-conforming grids, Annual Research Briefs, NASA Ames Research Center/Stanford University. Stanford, CA, 2001 , 353 - 366 [本文引用: 2]
[16]	Tseng YH , Ferziger JH . A ghost-cell immersed boundary method for flow in complex geometry . Journal of Computational Physics, 2003 , 192 ( 2 ): 593 - 623 [本文引用: 1]
[17]	Qu Y , Shi R , Batra RC . An immersed boundary formulation for simulating high-speed compressible viscous flows with moving solids . Journal of Computational Physics, 2018 , 354 : 672 - 691 [本文引用: 2]
[18]	Ghias R , Mittal R , Dong H . A sharp interface immersed boundary method for compressible viscous flows . Journal of Computational Physics, 2007 , 225 ( 1 ): 528 - 553 [本文引用: 1]
[19]	Chaudhuri A , Hadjadj A , Chinnayya A . On the use of immersed boundary methods for shock/obstacle interactions . Journal of Computational Physics, 2011 , 230 ( 5 ): 1731 - 1748 [本文引用: 3]
[20]	Bouchon F , Dubois T , James N . A second-order cut-cell method for the numerical simulation of 2D flows past obstacles . Computers& Fluids, 2012 , 65 : 80 - 91 [本文引用: 1]
[21]	Ikeno T , Kajishima T . Finite-difference immersed boundary method consistent with wall conditions for incompressible turbulent flow simulations . Journal of Computational Physics, 2007 , 226 ( 2 ): 1485 - 1508 [本文引用: 1]
[22]	Seo JH , Mittal R . A high-order immersed boundary method for acoustic wave scattering and low-Mach number flow-induced sound in complex geometries . Journal of Computational Physics, 2011 , 230 ( 4 ): 1000 - 1019 [本文引用: 1]
[23]	Haugen NEL , Kragset S . Particle impaction on a cylinder in a crossflow as function of Stokes and Reynolds numbers . Journal of Fluid Mechanics, 2010 , 661 : 239 - 261 [本文引用: 1]
[24]	Luo K , Zhuang Z , Fan J , et al . A ghost-cell immersed boundary method for simulations of heat transfer in compressible flows under different boundary conditions . International Journal of Heat and Mass Transfer, 2016 , 92 : 708 - 717 [本文引用: 2]
[25]	Luo K , Mao C , Zhuang Z , et al . A ghost-cell immersed boundary method for the simulations of heat transfer in compressible flows under different boundary conditions Part-II: Complex geometries . International Journal of Heat and Mass Transfer, 2017 , 104 : 98 - 111 [本文引用: 1]
[26]	Li XL , Fu DX , Ma YW , et al . Direct numerical simulation of shock/turbulent boundary layer interaction in a supersonic compression ramp . Science China Physics, Mechanics& Astronomy, 2010 , 53 ( 9 ): 1651 - 1658 [本文引用: 1]
[27]	Li X , Fu D , Ma Y . Direct numerical simulation of hypersonic boundary layer transition over a blunt cone with a small angle of attack . Physics of Fluids, 2010 , 22 ( 2 ): 025105 [本文引用: 1]
[28]	童福林 , 李新亮 , 唐志共 . 激波与转捩边界层干扰非定常特性数值分析 . 力学学报, 2017 , 49 ( 1 ): 93 - 104 [本文引用: 2] ( Tong Fulin , Li Xinliang , Tang Zhigong . Numerical analysis of unsteady motion in shock wave/transitional boundary layer interaction . Chinese Journal of Theoretical and Applied Mechanics, 2017 , 49 ( 1 ): 93 - 104 (in Chinese)) [本文引用: 2]
[29]	Steger JL , Warming RF . Flux vector splitting of the inviscid gasdynamic equations with application to finite-difference methods . Journal of Computational Physics, 1981 , 40 ( 2 ): 263 - 293 [本文引用: 1]
[30]	Jiang GS , Shu CW . Efficient implementation of weighted ENO schemes . Journal of Computational Physics, 1996 , 126 ( 1 ): 202 - 228 [本文引用: 2]
[31]	Franke R . Scattered data interpolation: Tests of some methods . Mathematics of Computation, 1982 , 38 ( 157 ): 181 - 200 [本文引用: 1]
[32]	Yang J , Liu Y , Lomax H . Computation of shock wave reflection by circular cylinders . AIAA Journal, 1987 , 25 ( 5 ): 683 - 689 [本文引用: 1]
[33]	Giepman RHM , Schrijer FFJ , Van Oudheusden BW . Flow control of an oblique shock wave reflection with micro-ramp vortex generators: Effects of location and size . Physics of Fluids, 2014 , 26 ( 6 ): 066101 [本文引用: 1]
[34]	Igra D , Takayama K , Igra O . Shock wave reflection over roughened wedges// 30th International Symposium on Shock Waves, Cham: Springer, 2017 : 621 - 625 [本文引用: 1]
[35]	Fang J , Yao Y , Zheltovodov AA , et al . Investigation of three dimensional shock wave turbulent boundary layer interaction initiated by a single fin . AIAA Journal, 2016 , 55 ( 2 ): 509 - 523 [本文引用: 1]
[36]	Bryson AE , Gross RWF . Diffraction of strong shocks by cones, cylinders, and spheres . Journal of Fluid Mechanics, 1961 , 10 ( 1 ): 1 - 16 [本文引用: 4]
[37]	Kaca J . An interferometric investigation of the diffraction of a planar shock wave over a semicircular cylinder . NASA STI/Recon Technical Report N, 1988 , 89 [本文引用: 1]
[38]	Constantinescu G , Squires K . Numerical investigations of flow over a sphere in the subcritical and supercritical regimes . Physics of Fluids, 2004 , 16 ( 5 ): 1449 - 1466 [本文引用: 1]
[39]	Rodriguez I , Borell R , Lehmkuhl O , et al . Direct numerical simulation of the flow over a sphere at Re= 3700. Journal of Fluid Mechanics, 2011 , 679 : 263 - 287 [本文引用: 1]
[40]	胡政敏 , 肖志坚 , 陶铸等 . 基于 LES 的亚临界和超临界圆球绕流数值研究 . 水动力学研究与进展, 2016 , 31 ( 4 ): 496 - 502 [本文引用: 1] ( Hu Zhengmin , Xiao Zhijian , Tao Zhu , et al . LES investigations of flow over a sphere in the subcritical and supercritical regimes . Chinese Journal of Hydrodynamics, 2016 , 31 ( 4 ): 496 - 502 (in Chinese)) [本文引用: 1]
[41]	Von Terzi DA , Sandberg RD , Fasel HF . Identification of large coherent structures in supersonic axisymmetric wakes . Computers& Fluids, 2009 , 38 ( 8 ): 1638 - 1650 [本文引用: 1]
[42]	Yun G , Kim D , Choi H . Vortical structures behind a sphere at subcritical Reynolds numbers . Physics of Fluids, 2006 , 18 ( 1 ): 015102 [本文引用: 1]
[43]	林孟达 , 崔桂香 , 张兆顺等 . 飞机尾涡演变及快速预测的大涡模拟研究 . 力学学报, 2017 , 49 ( 6 ): 1185 - 1200 [本文引用: 1] ( Lin Mengda , Cui Guixiang , Zhang Zhaoshun , et al . Large eddy simulation on the evolution and the fast-time prediction of aircraft wake vortices . Chinese Journal of Theoretical and Applied Mechanics, 2017 , 49 ( 6 ): 1185 - 1200 (in Chinese)) [本文引用: 1]
[44]	及春宁 , 花阳 , 许栋等 . 不同剪切率来流作用下柔性圆柱涡激振动数值模拟 . 力学学报, 2018 , 50 ( 1 ): 21 - 31 [本文引用: 1] ( Ji Chunning , Hua Yang , Xu Dong , et al . Numerical simulation of vortex-induced vibration of a flexible cylinder exposed to shear flow at different shear rates . Chinese Journal of Theoretical and Applied Mechanics, 2018 , 50 ( 1 ): 21 - 31 (in Chinese)) [本文引用: 1]
[45]	任安禄 , 李广望 , 邹建锋 . 中等雷诺数圆球绕流的数值研究 . 浙江大学学报(工学版), 2004 , 38 ( 5 ): 644 - 648 [本文引用: 1] ( Ren Anlu , Li Guangwang , Zou Jianfeng . Numerical study of uniform flow over sphere at intermediate Reynolds numbers . Journal of Zhejiang University( Engineering Science), 2004 , 38 ( 5 ): 644 - 648 (in Chinese)) [本文引用: 1]
[46]	邹建锋 , 任安禄 , 邓见 . 圆球绕流场的尾涡分析和升阻力研究 . 空气动力学报. 2004 , 22 ( 3 ): 303 - 308 [本文引用: 1] ( Zou Jianfeng , Ren Anlu , Deng Jian . Numerical investigations of wake and force for flow past a sphere . Acta Aerodynamica Sinica, 2004 , 22 ( 3 ): 303 - 308 (in Chinese)) [本文引用: 1]
[47]	Achenbach E . Experiments on the flow past spheres at very high Reynolds numbers . Journal of Fluid Mechanics, 1972 , 54 ( 3 ): 565 - 575 [本文引用: 1]

A novel approach of three-dimensional hybrid grid methodology: Part 1. Grid generation

2003

... 在计算流体动力学领域, 通常会遇到复杂几何绕流问题, 例如超声速飞行器气动力/热问题, 发动机内涡轮叶片流动换热研究, 燃烧室内旋流器稳流特性研究等. 面对任意复杂的几何外形, 高精度高效率的数值方法显得十分重要. 当前比较常用的处理方法是生成贴体的结构或非结构网格 ^{[

1
,

2
,

3
]}, 网格的形状与几何边界变形保持一致, 这样很容易施加边界条件, 采用高阶计算格式, 并在边界附近进行局部网格加密, 可以获得高分辨率的流动特征. 但往往生成贴体网格的过程需要耗费很大的计算资源, 同时针对复杂的几何外形, 生成高质量的网格并不容易, 从而极大地影响了数值计算精度, 特别是当处理动边界问题时, 需要在每一迭代时刻重新生成新的网格体系, 并利用插值方法将上一时刻流场解插值到新的网格系统下, 这无疑带来了很大的计算量. 浸没边界法 (immersed boundary method, IBM) 作为一种可行的选择, 能较好地处理复杂几何绕流问题, 它的主要优势在于不需要生成复杂的贴体网格, 只需要在规则的笛卡尔网格体系下考虑浸没边界对流体的影响, 这样在处理复杂静止或运动的几何边界问题上具有较高的效率, 同时由于采用规则的正交网格, 相对非结构网格来说, 极大地节省了计算资源 ^{[

4
]}, 同时在程序编写方面, 非常适用于并行处理, 这对于求解大规模高精度流场问题来说具有十分重要意义. ...

Fine-grained parallel algorithm for unstructured surface mesh generation

2015

Comparison of a spectral method with volume penalization and a finite volume method with body fitted grids for turbulent flows

2016

Immersed boundary methods for simulating fluid-structure interaction

2014

Flow patterns around heart valves: A numerical method

1972

... 浸没边界法首次被 Peskin ^{[

5
]}提出, 用于模拟笛卡尔网格下血液流过心脏瓣膜时的相互作用, 通过对动量方程施加附加力来代替瓣膜的存在, 并应用一种离散的 Dirac-delta 函数将拉格朗日作用力分散到欧拉网格点上. 后来 Goldstein 等 ^{[

6
]}, Saiki 和 Biringen ^{[

7
]}发展了 Peskin 的方法, 他们采用反馈力模型使得刚性静止壁面上流体速度为零; 但反馈力模型会带来虚假的数值震荡, 同时采用显式的时间步进格式, 严重限制了时间步长大小 ^{[

6
,

8
]}. 此外, Kolomenskiy 等 ^{[

9
,

10
]}采用惩罚模型对浸没边界进行处理, 但其惩罚常数依赖于所求解的具体问题, 影响了数值解的收敛性 ^{[

11
]}. 以上方法都是将附加力添加在连续 N-S 方程中, 与此不同的是 Mohd-Yusof ^{[

12
]}首次提出了直接附加力 (direct forcing) 方法, 将边界条件直接施加在离散的控制方程中, 得到离散的附加力. 这样, 控制方程中没有引入额外的附加项, 因此不存在数值刚性问题 ^{[

11
]}. Fadlun 等 ^{[

13
]}将此方法与基于交错网格的有限差分方法相结合, 运用于求解三维复杂问题, 展示了此方法相比反馈力模型更加高效. ...

Modeling a no-slip flow boundary with an external force field

1993

... [ 6 , 8 ]. 此外, Kolomenskiy 等 ^{[

9
,

10
]}采用惩罚模型对浸没边界进行处理, 但其惩罚常数依赖于所求解的具体问题, 影响了数值解的收敛性 ^{[

11
]}. 以上方法都是将附加力添加在连续 N-S 方程中, 与此不同的是 Mohd-Yusof ^{[

12
]}首次提出了直接附加力 (direct forcing) 方法, 将边界条件直接施加在离散的控制方程中, 得到离散的附加力. 这样, 控制方程中没有引入额外的附加项, 因此不存在数值刚性问题 ^{[

11
]}. Fadlun 等 ^{[

13
]}将此方法与基于交错网格的有限差分方法相结合, 运用于求解三维复杂问题, 展示了此方法相比反馈力模型更加高效. ...

Numerical simulation of a cylinder in uniform flow: Application of a virtual boundary method

1996

Stability characteristics of the virtual boundary method in three-dimensional applications

2003

A Fourier spectral method for the Navier-Stokes equations with volume penalization for moving solid obstacles

2009

A volume penalization method for incompressible flows and scalar advection-diffusion with moving obstacles

2012

Immersed boundary technique for turbulent flow simulations

2003

... [ 11 ]. Fadlun 等 ^{[

13
]}将此方法与基于交错网格的有限差分方法相结合, 运用于求解三维复杂问题, 展示了此方法相比反馈力模型更加高效. ...

... 从求解N-S方程角度看, 为了在流体内部施加边界条件以区分流体与固体之间的界面, 需要在动量方程右端增加作用力项

f

, 代表固体对流体的作用.Mohd-Yusof ^{[

11
]}最先提出了一种直接作用力方法, 将作用力项增加到离散的动量方程中 ...

... 这样

f_{i}

的作用使得在求解流场的每一时间步内均满足所定义的边界条件.这样处理的好处是所定义的速度边界条件直接施加在流、固界面上, 求解 N-S 方程不需计算新的增加项, 因而没有增加额外的 CPU 计算时间, 同时相比文献 ^{[

11
,

13
]}, 处理方法也不受时间步长的限制. ...

... 此外, 图 9 给出了

t = 22.9 s

时刻 C3 工况下流场的流线图, 左边是对固体内部除去虚拟点外所有固体节点速度场归零处理, 右边是对固体内部不作处理. 从中对比可以看出, 两种处理方法对外部流场的计算结果是一致的, 说明了本文采用虚拟单元方法无需对固体内部进行处理, 也验证了文献 ^{[

11
]}的说法. ...

Combined immersed boundary/B-spline methods for simulations of flow in complex geometries,

1997

... 本文采用虚拟单元浸没边界法 (ghost-cell immersed boundary method) 对高速可压缩问题进行求解, 其基本思想与文献 ^{[

12
]}的方法相似, 将边界信息直接通过固体几何内虚拟点加以反映, 不需要在控制方程中加入附加项, 因而能直接与已有的流体动力学求解器相结合, 在处理复杂界面方面有着广泛应用. 虚拟单元边界法的关键在于如何重构虚拟点的流场信息, 使得几何边界信息的刻画更加准确和稳定. 一般来说, 虚拟点沿着边界法向镜像外插是最常见的处理方式 ^{[

14
]}, 一方面诺依曼边界条件的施加更加方便, 另一方面避免了负值加权系数的产生 ^{[

15
]}. Tseng 和 Ferziger ^{[

16
]}系统分析了一阶线性重构、二次多项式重构和反距离加权重构在处理不同边界条件上的优缺点, 虽然虚拟点高阶重构具有较高的数值精度, 但会带来数值不稳定问题 ^{[

17
]}, 同时计算量也会增大. 除了重构格式外, 虚拟点层数和边界附近网格密度对计算精度和收敛性影响也较大. 对于边界附近的对流项处理, Ghias 等 ^{[

18
]}采用中心离散格式, 虚拟点的选取使其附近仅存在一个流体点. Chaudhuri 等 ^{[

19
]}采用 2 层虚拟单元, 使之与 5 阶加权基本无振荡 (weighted essentially non-oscillatory, WENO) 格式相匹配用于求解激波与障碍物之间相互作用. 从几何角度来看, 虚拟层数越多, 对应几何边界离散点越多, 有助于提高计算精度, 但虚拟点越远离边界, 重构解受边界条件影响越小. 通常情况下, 需要在边界附近进行网格加密, 使得边界附近的虚拟点增多, 从而边界条件的刻画越准确, 但无论网格的加密还是虚拟层数增加都会增加计算量. ...

Combined immersed-boundary finite-difference methods for three-dimensional complex flow simulations

2000

... 这样

f_{i}

A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries

2008

RANS solvers with adaptive structured boundary non-conforming grids,

2001

... 通常情况下, 镜像点 (IP) 并不与网格节点重合, 因而需要通过内部流体点对 IP 点物理信息进行插值重构.最常用的是线性插值, Luo 等 ^{[

24
]}根据 IP 点所在单元 4 个流体节点, 采用双线性插值求得 IP 点流场信息, Majumdar 等 ^{[

15
]}对比了双线性插值和二次曲线插值两种方法的区别, 发现两者的精度并没有明显差异.Chaudhuri 采用 Franke ^{[

31
]}提出的一种反距离加权插值方法对镜像点流场信息进行重构, 保证了重构解在极值点附近的光滑性 ...

A ghost-cell immersed boundary method for flow in complex geometry

2003

An immersed boundary formulation for simulating high-speed compressible viscous flows with moving solids

2018

... 此外, 超声速或高超声速绕流问题, 往往伴随着激波/障碍物、激波/边界层、激波/激波之间相互作用, 因此网格疏密、计算格式及其虚拟点的分布将对整个流场细节的模拟非常重要. 但目前虚拟单元边界法大多数侧重不可压缩问题的研究 ^{[

20
,

21
]}, 对于高雷诺数可压缩黏性问题, 特别是高马赫数问题的运用很少. Qu 等 ^{[

17
]}改进了 Seo 和 Mittal ^{[

22
]}提出的高阶虚拟点重构格式, 将其用于不同马赫数 (1.20~2.37)下超声速可压缩问题的研究, 并消除了该格式在处理移动边界时出现数值震荡问题. Haugen 和 Kragset ^{[

23
]}将虚拟浸没边界法与一种高阶可压缩代码 (PENCIL) 相结合, 对横向流中颗粒撞击圆柱过程进行了直接数值模拟. Luo 等 ^{[

24
,

25
]}测试了虚拟单元边界法在处理可压缩流动传热问题上的优势, 其研究雷诺数和马赫数分别在 5000 和 0.5 以下. 由于虚拟单元浸没边界法在处理复杂几何问题上的优势, 再加上笛卡尔网格体系下非常容易构造高阶有限差分格式和大规模并行操作, 因而虚拟单元边界法在处理超声速或高超声速下可压缩复杂流场问题方面具有很大的潜力. ...

A sharp interface immersed boundary method for compressible viscous flows

2007

On the use of immersed boundary methods for shock/obstacle interactions

2011

... 本文从固体几何边界细化的角度, 考虑到边界网格疏密对计算精度的影响, 对基于文献 ^{[

19
]}的虚拟单元边界法进行改进, 并与一种高阶精度的有限差分方法相结合, 运用于求解超声速可压缩流动问题. 本文组织结构如下: 在 2.1 节以三维空间激波/圆柱相互作用 (马赫数为 2.81) 为例对本文算法的精度进行对比验证, 在 2.2 节以三维空间超声速球体绕流 (马赫数为 1.2) 为例进一步说明本文算法对于高雷诺数下超声速湍流和激波结构的捕捉能力, 并对球体绕流的自由剪切层演化过程进行细致分析. 第 3 节对全文进行小结. ...

... 本文参考了文献 ^{[

19
]}的研究工作, 如图 1 所示.在进行数值计算之前, 首先需要识别出虚拟边界, 并存储其节点上的有效信息.整个计算区域

Ω

由固体区域

Ω_{s olid}

和流体区域

Ω_{f luid}

组成, 为了分离物理几何边界, 需要对所有网格节点进行标记 ...

A second-order cut-cell method for the numerical simulation of 2D flows past obstacles

2012

Finite-difference immersed boundary method consistent with wall conditions for incompressible turbulent flow simulations

2007

A high-order immersed boundary method for acoustic wave scattering and low-Mach number flow-induced sound in complex geometries

2011

Particle impaction on a cylinder in a crossflow as function of Stokes and Reynolds numbers

2010

A ghost-cell immersed boundary method for simulations of heat transfer in compressible flows under different boundary conditions

2016

A ghost-cell immersed boundary method for the simulations of heat transfer in compressible flows under different boundary conditions Part-II: Complex geometries

2017

Direct numerical simulation of shock/turbulent boundary layer interaction in a supersonic compression ramp

2010

... 三维可压缩流动的基本方程是 Navier-Stokes (N-S) 方程, 本文采用直接数值模拟 (direct numerical simulation) 方法对其进行求解.在数值求解过程中通常将控制方程进行无量纲处理, 本文选取来流参数

ρ_{\infty}, u_{\infty}, T_{\infty}, ρ_{\infty} u_{\infty}

分别对整个流场物理参量

P, V = (u, v, w), T, p

进行无量纲化, 流场特征长度为

L

, 特征时间为

L / u_{\infty}

, 由此笛卡尔坐标系下三维非定常可压缩 N-S 方程 ^{[

26
,

27
,

28
]}如下 ...

Direct numerical simulation of hypersonic boundary layer transition over a blunt cone with a small angle of attack

2010

ρ_{\infty}, u_{\infty}, T_{\infty}, ρ_{\infty} u_{\infty}

分别对整个流场物理参量

P, V = (u, v, w), T, p

进行无量纲化, 流场特征长度为

L

, 特征时间为

L / u_{\infty}

, 由此笛卡尔坐标系下三维非定常可压缩 N-S 方程 ^{[

26
,

27
,

28
]}如下 ...

激波与转捩边界层干扰非定常特性数值分析

2017

ρ_{\infty}, u_{\infty}, T_{\infty}, ρ_{\infty} u_{\infty}

分别对整个流场物理参量

P, V = (u, v, w), T, p

进行无量纲化, 流场特征长度为

L

, 特征时间为

L / u_{\infty}

, 由此笛卡尔坐标系下三维非定常可压缩 N-S 方程 ^{[

26
,

27
,

28
]}如下 ...

... 具体分解得到的正、负通量

F_{i}^{+}

和

F_{i}^{-}

表达式可以参考文献 ^{[

28
]}. ...

激波与转捩边界层干扰非定常特性数值分析

2017

ρ_{\infty}, u_{\infty}, T_{\infty}, ρ_{\infty} u_{\infty}

分别对整个流场物理参量

P, V = (u, v, w), T, p

进行无量纲化, 流场特征长度为

L

, 特征时间为

L / u_{\infty}

, 由此笛卡尔坐标系下三维非定常可压缩 N-S 方程 ^{[

26
,

27
,

28
]}如下 ...

... 具体分解得到的正、负通量

F_{i}^{+}

和

F_{i}^{-}

表达式可以参考文献 ^{[

28
]}. ...

Flux vector splitting of the inviscid gasdynamic equations with application to finite-difference methods

1981

... 根据扰动波传播的方向性, 对无黏通量

F_{i} (i = 1,2, 3)

使用 Steger-Warming 分裂 ^{[

29
]}, 将其分解成正通量

F_{i}^{+}

和负通量

F_{i}^{-}

, 然后再分别对其空间导数进行差分离散, 这样保证了熵增的性质和格式的稳定性. ...

Efficient implementation of weighted ENO schemes

1996

... 对于 N-S 方程中非线性无黏项, 本文采用 Jiang 和 Shu ^{[

30
]}提出的 5 阶 WENO 有限差分格式, 使得在激波捕捉问题上具有较高的精度和数值稳定性, 其空间导数

f_{j}^{'}

表示为 ...

... 其中,

ω_{m}

是权值系数, 具体取值可以参考文献 ^{[

30
]}. ...

Scattered data interpolation: Tests of some methods

1982

Computation of shock wave reflection by circular cylinders

1987

... 在超声速流动问题中, 激波反射现象在过去一直被大家广泛研究 ^{[

32
,

33
,

34
,

35
]}. 当一道入射激波在传播过程中遇到一个圆柱障碍物时, 首先会在圆柱表面发生常规激波反射, 形成交汇的二激波结构, 随着交汇点在圆柱表面逐渐向后移动, 常规反射会转捩成马赫反射, 形成三分叉点的三激波结构. 下面以二维空间入射激波与圆柱相互作用过程为例, 对本文的虚拟单元浸没边界算法进行验证. 平面入射激波马赫数为 2.81, 整个计算域长宽尺寸分别为

L_{x} \times L_{y} = 240 m m \times 60 m m

, 圆柱直径

D

为 10 mm. 进、出口分别是超声速进、出口边界条件, 上、下壁面均为无滑移的绝热壁面. 整个计算区域采用均匀的正交网格, 采用 5 种不同的网格尺寸算例来检验网格的无关性, 如表 1 所示. ...

Flow control of an oblique shock wave reflection with micro-ramp vortex generators: Effects of location and size

2014

L_{x} \times L_{y} = 240 m m \times 60 m m

, 圆柱直径

D

Shock wave reflection over roughened wedges//

2017

L_{x} \times L_{y} = 240 m m \times 60 m m

, 圆柱直径

D

Investigation of three dimensional shock wave turbulent boundary layer interaction initiated by a single fin

2016

L_{x} \times L_{y} = 240 m m \times 60 m m

, 圆柱直径

D

Diffraction of strong shocks by cones, cylinders, and spheres

1961

... 图 4是不同时刻密度等值面云图, 显示了激波绕射圆柱的演化过程, 包括马赫波的形成 (a)、马赫波沿着壁面传播 (b)、马赫波在尾迹处碰撞 (c) 和二次 3- 波点的出现 (d). 整个流场结构能清晰的识别出入射激波 (I.S.)、反射弓形激波 (R.S.)、马赫波 (M.S.1 和 M.S.2)、三波点 (T.P.1 和 T.P.2)、滑移线 (S.L.1 和 S.L.2) 等复杂波系结构, 同时还能观察到在滑移线后流场卷吸形成小涡结构 (V.). 特别是在圆柱后尾迹处, 存在上、下两道马赫激波相互碰撞, 并在尾迹中发展形成二次三波系结构. 与 Bryson 和 Gross ^{[

36
]}实验可视化图比较, 计算结果较好的捕捉到了流场的特征. ...

... 激波反射流场纹影图对比 ((a)和(b): 数值摸拟结果, (c)和(d): 实验结果 ^{[

36
]}) ...

... Schlieren pictures for the flow evolution of the shock reflection ((a)and(b): numerical simulation, (c)and(d): experimental results ^{[

36
]}) ...

... 图 7是两种虚拟单元边界法壁面压力系数的数据对比, 可以看出对于改进的虚拟单元方法, 其相对数值误差比原始的虚拟单元方法小. 对于网格数目最稀疏的 C1 和 C4 两种情况, 虽然圆柱前缘滞止点的压力系数与基准解相差不大, 但在圆柱后缘 3, 4 点位置, 两种方法均与基准解有一定的差距. 在 C2 和 C5 情况时, 改进的虚拟单元方法其计算结果与基准解大致吻合. 这说明施加基于加力点的边界条件不仅能细化几何边界起到强化边界信息的作用, 而且改善了数值解的计算精度.但当网格数目本身很密情况下, 这种效果不太明显. 图 8 给出 C2, C3, C5 三种情况下三波点（T.P.1 和 T.P.2）的轨迹图, 从中可以看出基准工况的计算结果与文献 ^{[

36
,

37
]}的实验结果吻合较好, 证明了本文采用的基于改进虚拟单元边界法的高阶有限差分求解器能准确的模拟激波与障碍物之间相互作用. ...

An interferometric investigation of the diffraction of a planar shock wave over a semicircular cylinder

1988

Numerical investigations of flow over a sphere in the subcritical and supercritical regimes

2004

... 为了进一步验证本文所提出算法的稳定性, 对三维超声速圆球绕流算例进行研究. 过去几十年, 许多学者致力于亚临界和超临界雷诺数下圆球绕流流动分离现象的研究 ^{[

38
,

39
,

40
,

41
]}, 特别关注于剪切层不稳定发展和尾迹涡旋结构的演化规律探索 ^{[

42
,

43
,

44
,

45
,

46
]}. 本文超声速来流马赫数马赫数为 1.2, 雷诺数

Re = 1.62 \times 1 0^{5}

, 圆球直径D= 10 mm, 计算几何模型长度为

L_{x} \times L_{y} \times L_{z} = 24 D \times 6 D \times 6 D

.为减少网格计算量, 计算域采用分块加密网格, 球体附近区域的网格尺寸

Δ δ

最小, 以提高边界流场的分辨率, 如图 10 所示. 根据局部加密尺寸

Δ δ

不同, 比较 8 个计算工况下改进虚拟单元边界法的计算精度. ...

Direct numerical simulation of the flow over a sphere at

2011

Re = 1.62 \times 1 0^{5}

, 圆球直径D= 10 mm, 计算几何模型长度为

L_{x} \times L_{y} \times L_{z} = 24 D \times 6 D \times 6 D

.为减少网格计算量, 计算域采用分块加密网格, 球体附近区域的网格尺寸

Δ δ

最小, 以提高边界流场的分辨率, 如图 10 所示. 根据局部加密尺寸

Δ δ

不同, 比较 8 个计算工况下改进虚拟单元边界法的计算精度. ...

LES 的亚临界和超临界圆球绕流数值研究

2016

Re = 1.62 \times 1 0^{5}

, 圆球直径D= 10 mm, 计算几何模型长度为

L_{x} \times L_{y} \times L_{z} = 24 D \times 6 D \times 6 D

.为减少网格计算量, 计算域采用分块加密网格, 球体附近区域的网格尺寸

Δ δ

最小, 以提高边界流场的分辨率, 如图 10 所示. 根据局部加密尺寸

Δ δ

不同, 比较 8 个计算工况下改进虚拟单元边界法的计算精度. ...

LES 的亚临界和超临界圆球绕流数值研究

2016

Re = 1.62 \times 1 0^{5}

, 圆球直径D= 10 mm, 计算几何模型长度为

L_{x} \times L_{y} \times L_{z} = 24 D \times 6 D \times 6 D

.为减少网格计算量, 计算域采用分块加密网格, 球体附近区域的网格尺寸

Δ δ

最小, 以提高边界流场的分辨率, 如图 10 所示. 根据局部加密尺寸

Δ δ

不同, 比较 8 个计算工况下改进虚拟单元边界法的计算精度. ...

Identification of large coherent structures in supersonic axisymmetric wakes

2009

Re = 1.62 \times 1 0^{5}

, 圆球直径D= 10 mm, 计算几何模型长度为

L_{x} \times L_{y} \times L_{z} = 24 D \times 6 D \times 6 D

.为减少网格计算量, 计算域采用分块加密网格, 球体附近区域的网格尺寸

Δ δ

最小, 以提高边界流场的分辨率, 如图 10 所示. 根据局部加密尺寸

Δ δ

不同, 比较 8 个计算工况下改进虚拟单元边界法的计算精度. ...

Vortical structures behind a sphere at subcritical Reynolds numbers

2006

Re = 1.62 \times 1 0^{5}

, 圆球直径D= 10 mm, 计算几何模型长度为

L_{x} \times L_{y} \times L_{z} = 24 D \times 6 D \times 6 D

.为减少网格计算量, 计算域采用分块加密网格, 球体附近区域的网格尺寸

Δ δ

最小, 以提高边界流场的分辨率, 如图 10 所示. 根据局部加密尺寸

Δ δ

不同, 比较 8 个计算工况下改进虚拟单元边界法的计算精度. ...

飞机尾涡演变及快速预测的大涡模拟研究

2017

Re = 1.62 \times 1 0^{5}

, 圆球直径D= 10 mm, 计算几何模型长度为

L_{x} \times L_{y} \times L_{z} = 24 D \times 6 D \times 6 D

.为减少网格计算量, 计算域采用分块加密网格, 球体附近区域的网格尺寸

Δ δ

最小, 以提高边界流场的分辨率, 如图 10 所示. 根据局部加密尺寸

Δ δ

不同, 比较 8 个计算工况下改进虚拟单元边界法的计算精度. ...

飞机尾涡演变及快速预测的大涡模拟研究

2017

Re = 1.62 \times 1 0^{5}

, 圆球直径D= 10 mm, 计算几何模型长度为

L_{x} \times L_{y} \times L_{z} = 24 D \times 6 D \times 6 D

.为减少网格计算量, 计算域采用分块加密网格, 球体附近区域的网格尺寸

Δ δ

最小, 以提高边界流场的分辨率, 如图 10 所示. 根据局部加密尺寸

Δ δ

不同, 比较 8 个计算工况下改进虚拟单元边界法的计算精度. ...

不同剪切率来流作用下柔性圆柱涡激振动数值模拟

2018

Re = 1.62 \times 1 0^{5}

, 圆球直径D= 10 mm, 计算几何模型长度为

L_{x} \times L_{y} \times L_{z} = 24 D \times 6 D \times 6 D

.为减少网格计算量, 计算域采用分块加密网格, 球体附近区域的网格尺寸

Δ δ

最小, 以提高边界流场的分辨率, 如图 10 所示. 根据局部加密尺寸

Δ δ

不同, 比较 8 个计算工况下改进虚拟单元边界法的计算精度. ...

不同剪切率来流作用下柔性圆柱涡激振动数值模拟

2018

Re = 1.62 \times 1 0^{5}

, 圆球直径D= 10 mm, 计算几何模型长度为

L_{x} \times L_{y} \times L_{z} = 24 D \times 6 D \times 6 D

.为减少网格计算量, 计算域采用分块加密网格, 球体附近区域的网格尺寸

Δ δ

最小, 以提高边界流场的分辨率, 如图 10 所示. 根据局部加密尺寸

Δ δ

不同, 比较 8 个计算工况下改进虚拟单元边界法的计算精度. ...

中等雷诺数圆球绕流的数值研究

2004

Re = 1.62 \times 1 0^{5}

, 圆球直径D= 10 mm, 计算几何模型长度为

L_{x} \times L_{y} \times L_{z} = 24 D \times 6 D \times 6 D

.为减少网格计算量, 计算域采用分块加密网格, 球体附近区域的网格尺寸

Δ δ

最小, 以提高边界流场的分辨率, 如图 10 所示. 根据局部加密尺寸

Δ δ

不同, 比较 8 个计算工况下改进虚拟单元边界法的计算精度. ...

中等雷诺数圆球绕流的数值研究

2004

Re = 1.62 \times 1 0^{5}

, 圆球直径D= 10 mm, 计算几何模型长度为

L_{x} \times L_{y} \times L_{z} = 24 D \times 6 D \times 6 D

.为减少网格计算量, 计算域采用分块加密网格, 球体附近区域的网格尺寸

Δ δ

最小, 以提高边界流场的分辨率, 如图 10 所示. 根据局部加密尺寸

Δ δ

不同, 比较 8 个计算工况下改进虚拟单元边界法的计算精度. ...

圆球绕流场的尾涡分析和升阻力研究

2004

Re = 1.62 \times 1 0^{5}

, 圆球直径D= 10 mm, 计算几何模型长度为

L_{x} \times L_{y} \times L_{z} = 24 D \times 6 D \times 6 D

.为减少网格计算量, 计算域采用分块加密网格, 球体附近区域的网格尺寸

Δ δ

最小, 以提高边界流场的分辨率, 如图 10 所示. 根据局部加密尺寸

Δ δ

不同, 比较 8 个计算工况下改进虚拟单元边界法的计算精度. ...

圆球绕流场的尾涡分析和升阻力研究

2004

Re = 1.62 \times 1 0^{5}

, 圆球直径D= 10 mm, 计算几何模型长度为

L_{x} \times L_{y} \times L_{z} = 24 D \times 6 D \times 6 D

.为减少网格计算量, 计算域采用分块加密网格, 球体附近区域的网格尺寸

Δ δ

最小, 以提高边界流场的分辨率, 如图 10 所示. 根据局部加密尺寸

Δ δ

不同, 比较 8 个计算工况下改进虚拟单元边界法的计算精度. ...

Experiments on the flow past spheres at very high Reynolds numbers

1972

... 图 12是球体中截面的壁面平均压力系数和摩擦系数与实验结果 ^{[

47
]}的对比, 可以看出, 在球体迎风侧, 各工况计算结果与实验值相差不大. 在背风侧, 随着球体附近网格依次加密, 数值解的误差逐渐减小, 说明流场解的精度随着网格分辨率的增大而提高. 当球体附近网格尺度达到 0.012D 即 S4 工况, 数值解与实验结果吻合较好, 因此可以将 S4 结果作为基准解, 对 S1~S8 工况整个计算域内的计算误差进行对比分析, 定义流场解误差的L2范数 ...

〈

〉

改进虚拟边界算法在超声速流动问题求解中的应用

AN IMPROVED GHOST-CELL IMMERSED BOUNDARY METHOD FOR SOLVING SUPERSONIC FLOW PROBLEMS

引言

1 均匀数学模型及计算方法

1.1 三维非定常可压缩湍流控制方程

1.2 特征空间下无黏项流通矢量分裂 (flux vector splitting)

1.3 时空离散处理

1.4 一种改进的虚拟单元浸没边界法

2 结果与讨论

2.1 激波反射现象

2.2 超声速三维球体绕流

3 结 论

参考文献 文献选项 原文顺序 文献年度倒序 文中引用次数倒序 被引期刊影响因子

3 结论

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子