, 邹建锋, 郑耀浙江大学航空航天学院, 杭州 310027
AN IMPROVED GHOST-CELL IMMERSED BOUNDARY METHOD FOR SOLVING SUPERSONIC FLOW PROBLEMS
ZhangYang, ZouJianfeng, ZhengYao中图分类号:V211.3
文献标识码:A
通讯作者:
收稿日期:2017-12-25
接受日期:2018-03-15
网络出版日期:2018-06-10
版权声明:2018《力学学报》编辑部《力学学报》编辑部 所有
基金资助:
展开
摘要
关键词：
Abstract
Keywords：
-->0
PDF (19393KB)元数据多维度评价相关文章收藏文章
本文引用格式导出EndNoteRisBibtex收藏本文-->
引言
在计算流体动力学领域, 通常会遇到复杂几何绕流问题, 例如超声速飞行器气动力/热问题, 发动机内涡轮叶片流动换热研究, 燃烧室内旋流器稳流特性研究等. 面对任意复杂的几何外形, 高精度高效率的数值方法显得十分重要. 当前比较常用的处理方法是生成贴体的结构或非结构网格[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) 方法对其进行求解.在数值求解过程中通常将控制方程进行无量纲处理, 本文选取来流参数其中
${\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}$
其中是
$\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)$
流体动力学黏度系数
其中,
1.2 特征空间下无黏项流通矢量分裂 (flux vector splitting)
在实际的数值计算中, 考虑到流体黏性的影响, 往往需要在固体壁面附近进行局部网格加密.因而为了计算的方便, 需要将流动控制方程从笛卡尔坐标空间其中,
根据扰动波传播的方向性, 对无黏通量
$\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)
其中,
其中
具体分解得到的正、负通量
1.3 时空离散处理
对于 N-S 方程中非线性无黏项, 本文采用 Jiang 和 Shu[30]提出的 5 阶 WENO 有限差分格式, 使得在激波捕捉问题上具有较高的精度和数值稳定性, 其空间导数其中,h是网格间距,${\pmb F}_{j+1/2}=\sum\limits^{4}_{m=1}\omega_{m}{\pmb P}_{m}$
其中,
对于黏性通量项
对于时间项推进, 本文采用 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方程角度看, 为了在流体内部施加边界条件以区分流体与固体之间的界面, 需要在动 量方程右端增加作用力项其中,
假如在
这样
通常对于复杂曲线几何来说, 网格节点位置并不一定正好落在几何表面, 这就需要对几何边界进行局部解重构.一种方法是采用虚拟单元, 在固体几何内部构造虚拟边界, 将物理边界上的信息通过内部流体节点, 外推至虚拟边界上, 从而保证流、固界面满足所需要的物理边界条件.下面以二维模型为例, 阐述虚拟单元浸没边界法的基本过程.
本文参考了文献 [19] 的研究工作, 如图 1 所示.在进行数值计算之前, 首先需要识别出虚拟边界, 并存储其节点上的有效信息.整个计算区域
$\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)
虚拟节点上的流场信息
但在数值实践中发现, 当几何边界附近的网格不太密集的情况下, 采用虚拟点的镜像关系, 对于几何边界的细化来说精度往往达不到要求, 尽管可以通过增加虚拟点层数来进一步细化几何边界, 但越远离边界的虚拟点, 采用镜像关系的线性插值显然是不适合的.这就需要在几何边界附近进行局部网格加密, 但增加了求解计算量.为了解决这一问题, 本文在紧靠边界区域拓展了 1 层流体作用点 (forcing point, FP), 同样采用镜像关系的线性插值, 重构这一层流体作用点的解.如图 2 对比了两种方法下边界曲线 AB 的细化情况.由于增加了 1 层流体作用点, AB 段离散增加至 14 段, 同时相较于增加虚拟点层数, 边界离散点 (boundary point, BP) 分布更趋于均匀
通常情况下, 镜像点 (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)
其中
显示原图|下载原图ZIP|生成PPT图 1虚拟点镜像处理
-->Fig. 1Immersed 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)
其中
显示原图|下载原图ZIP|生成PPT图 2浸没几何边界的细化
-->Fig. 2Immersed boundary subdivision
-->
对于上面提出的改进虚拟单元浸没边界法很容易扩展至 3 维空间, 只需将插值点限定在以 IP 点为球心,R为半径的球体内. 此外, 由于整个计算区域采用笛卡尔坐标系下的正交网格, 使得在信息传递接口并行编程计算方面比较容易实现, 一方面可以应用于大规模并行数值计算, 研究复杂的流动现象及其物理机理;另一方面, 由于网格系统并不随着内部几何模型的改变而需要重新划分, 这使得虚拟单元浸没边界法非常适用于优化设计.
显示原图|下载原图ZIP|生成PPT图3改进的浸没边界法
-->Fig. 3Modified immersed boundary treatment
-->
2 结果与讨论
2.1 激波反射现象
在超声速流动问题中, 激波反射现象在过去一直被大家广泛研究[32,33,34,35]. 当一道入射激波在传播过程中遇到一个圆柱障碍物时, 首先会在圆柱表面发生常规激波反射, 形成交汇的二激波结构, 随着交汇点在圆柱表面逐渐向后移动, 常规反射会转捩成马赫反射, 形成三分叉点的三激波结构. 下面以二维空间入射激波与圆柱相互作用过程为例, 对本文的虚拟单元浸没边界算法进行验证. 平面入射激波马赫数为 2.81, 整个计算域长宽尺寸分别为Table 1
表 1
表 1不同疏密的网格尺寸
Table 1Calculation mesh size in different cases
| Case | ||
|---|---|---|
| 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): 实验结果[
-->Fig. 4Schlieren pictures for the flow evolution of the shock reflection ((a)and(b): numerical simulation, (c)and(d): experimental results[
-->
图 5展示了整个激波绕射过程中圆柱表面 4 个位置的压力系数
显示原图|下载原图ZIP|生成PPT图 5圆柱表面压力系数变化
-->Fig. 5Pressure coefficients on the cylindrical surface in different spatial resolutions
-->
前文提到, 本文对虚拟单元浸没边界法的改进之处在于边界的细化, 通过考虑同时将边界条件关系施加到固体区域内部的虚拟点和固体区域附近的加力点, 可以对边界附近网格起到类似局部加密作用, 下面采用 3 种不同网格尺度对该方法的局部加密效果进行测试.
表 2展示了 6 种不同的计算工况, 其中虚拟点表示仅采用固体几何内虚拟点施加边界条件, 虚拟 + 加力点表示同时采用边界内外的虚拟和加力点施加边界条件(简称Ref.).
Table 2
表 2
表 2不同工况下网格尺寸和边界处理方法
Table 2Calculation conditions with different implementation methods
| Case | 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. 6Comparison of Mach number contours
-->
显示原图|下载原图ZIP|生成PPT图 7改进虚拟边界法的圆柱表面压力系数对比
-->Fig. 7Comparison of pressure coefficients with modified ghost-cell methods
-->
显示原图|下载原图ZIP|生成PPT图 8三波点轨迹图
-->Fig. 8Comparison of triple point trajectory
-->
此外, 图 9 给出了
显示原图|下载原图ZIP|生成PPT图 9流线图对比 (上: 内部特殊处理; 下: 内部不做处理)
-->Fig. 9Comparison 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, 雷诺数
显示原图|下载原图ZIP|生成PPT图 10光滑球体绕流计算域示意图
-->Fig. 10Computational domain of the smooth sphere
-->
图 11展示了 S4 工况下瞬态流场结构数值纹影图, 在球体正前方形成一道脱体激波, 自由剪切层从球体表面分离后向下游发展, 并与球体后尾迹结构相互作用, 导致剪切层内涡旋结构向外卷起或向内收缩, 在外部超声速流动作用下形成一系列向外传播的压缩波系 (A处), 压缩波系叠加产生一道强的尾激波结构 (B处). 在尾激波后剪切层出现了间断结构 (C处), 交替间断的结构逐渐演化成大尺度的涡结构. 在尾激波前一段距离内, 剪切层最先开始变形扭曲, 出现失稳现象. 这是因为一系列弱的压缩波向上游传播对剪切层进行了预压缩作用. 变形的剪切层、尾迹与激波之间的相互作用更加剧了剪切层结构的失稳, 导致流场结构的演化更加复杂.
显示原图|下载原图ZIP|生成PPT图 11中截面上瞬态流场密度梯度数值分布云图
-->Fig. 11Transient density contours in the middle section
-->
Table 3
表 3
表 3不同条件的计算工况
Table 3Eight cases with different implementation methods
| Case | 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. 12Grid independence check
-->
其中
显示原图|下载原图ZIP|生成PPT图 13流向速度u和球体表面压力系数Cp的误差L2范数
-->Fig. 13Global errors in u and Cp versus number of mesh points
-->
2.2.2 剪切层湍流结构定性分析
为了对整个流场结构的演化过程进行细致分析, 特别关注自由剪切层失稳特性的观察, 图 14 给出了不同时刻下球体绕流后等马赫线分布. 根据流场马赫数特征, 可以把整个流场分为 4 个部分, 即脱体激波与声速线 (马赫数为1) 之间的亚声速区域, 剪切层外声速线与尾激波之间的超声速区域, 剪切层内球体后方的尾迹结构, 尾激波后流场尾迹涡旋结构.从图 14 可以看出, 在尾激波前后区域, 剪切层结构出现较为明显的扭曲变形, 这种扭曲扰动逐渐向上游剪切层传播, 形成小的波浪状剪切层结构, 这是由于压缩波向上游传播的缘故. 尾激波后的剪切层结构出现了较大的弯曲变形, 并呈现周期性变化. 剪切层在发生扭曲变形的同时, 逐渐向内收缩, 将尾迹流场分成 2 部分, 一部分是球体正后方剪切层内未受扰动的尾迹结构;另一部分是向下游发展的尾迹结构, 这部分结构是上游尾迹、变形剪切层和压缩波系共同作用的结果, 且随着流场的演化, 这部分尾迹结构逐渐与上游剪切层脱落, 形成向下游发展的不同尺度脱落涡旋结构. 由此构成了一个完整的剪切层内涡旋脱落演化过程.为了进一步刻画剪切层结构的变化, 选取
显示原图|下载原图ZIP|生成PPT图 14不同时刻流场马赫数分布等值图
-->Fig. 14Mach number contours with different time
-->
前面已经指出, 剪切层附近的压缩波系向上游传播, 使得剪切层出现小波浪形的弯曲, 触发了剪切层流动的不稳定, 随后小变形剪切层、尾激波和尾迹之间的相互作用加剧了剪切层的弯曲变形, 直至最终失稳破碎成间断的涡结构. 图 15 给出了球体下游不同位置的周向速度梯度分布云图, 可以看出剪切层在最初阶段为完整的环形, 但随着流场向下游发展, 局部出现个别卷起的涡 (如图
显示原图|下载原图ZIP|生成PPT图 15流向不同位置横截面的数值云图 (物理量：
-->Fig.15Numerical schlieren pictures in different streamwise sections (
-->
显示原图|下载原图ZIP|生成PPT图 16三维流场的涡结构演化
-->Fig.16Three-dimensional vortex structures
-->
2.2.3 剪切层湍流结构定量分析
为了对剪切层的变形失稳过程进行定量分析, 考察剪切层厚度和褶皱因子的变化规律. 剪切层厚度
其中
图 17所示是剪切层厚度沿流向的变化情况, 可以看出剪切层厚度的变化呈现 3 个阶段：第 I 阶段为厚度的线性增长, 这是由于剪切层受内部尾迹回流影响, 再加上外面受到弱的压缩波扰动作用所导致的;第 II 阶段剪切层厚度呈现无规律的波动状态, 这时剪切层受到尾激波和尾迹的共同作用, 出现复杂的弯曲、缠绕和变形; 第 III 阶段剪切层厚度变化不大, 对应着尾迹涡旋向后逐步发展的过程. 需要注意的是, 虽然在剪切层的失稳过程中剪切层厚度变化处于上下波动状态, 但其变形扭曲却呈规律性的增长, 定义剪切层的褶皱因子
显示原图|下载原图ZIP|生成PPT图 17剪切层厚度的变化
-->Fig. 17Change of the shear layer thickness
-->
这里考虑流向不同横截面
显示原图|下载原图ZIP|生成PPT图 18剪切层褶皱因子的变化
-->Fig. 18Wrinkling factor of the shear layer
-->
褶皱因子的不断增长代表着流场剪切层与内部低速流体、外部高速流体的接触面积增大, 高低速流体的动量交换更加剧烈. 与此同时, 剪切层表面曲率变化越大, 容易产生向外卷起的涡结构, 在外部超声速流动下形成向外传播的压缩波, 进而对附近的剪切层产生干扰作用. 在
图 19为剪切层速度分布云图.图 20 和图 21 给出了剪切层内湍流雷诺应力的演化规律, 剪切层的位置取为图 19 中虚线, 在 D 处剪切层开始出现波浪状变形褶皱, S 处为尾激波所在位置.从湍流速度均方根分布可以看出, 3个方向的雷诺正应力呈现大致相同的特征趋势走向, 与剪切层厚度的演化规律一致, 第 I 阶段为雷诺正应力的线性增长阶段, 自由剪切层为稳定的平直状态; 第 II 阶段由于压缩波系的作用使得湍流结构更为复杂, 湍流脉动强度经历一个混沌的震荡过程, 并且在尾激波位置流向脉动强度达到最大值; 第 III 阶段为尾激波后涡旋的脱落向前推进过程. 还可以观察到, 超声速球体绕流的湍流雷诺正应力表现出明显的各项异性, 雷诺应力中流向正应力最大, 这说明在剪切层的演化过程中湍流结构的流向脉动占主导地位.此外, z 方向的正应力最大值在尾激波后一段距离内取得, 说明了激波的压缩作用对不同方向雷诺正应力的影响存在空间迟滞效应, 进一步阐述了剪切层内湍流结构的各项异性.从湍流雷诺剪切应力分布可以看出, 最大雷诺剪切应力仍然出现在尾激波所在位置 (图中S), 不存在空间上的滞后, 但雷诺剪切应力的发展在激波附近存在一个先增加后减小的突变过程, 这一点与雷诺正应力的演化规律不同. 与此同时, 在
显示原图|下载原图ZIP|生成PPT图 19剪切层速度分布云图
-->Fig. 19Streamwise velocity distribution of the shear layer
-->
显示原图|下载原图ZIP|生成PPT图 20剪切层雷诺正应力分布
-->Fig. 20Reynolds normal stress distribution of the shear layer
-->
显示原图|下载原图ZIP|生成PPT图 21剪切层主切应力分布
-->Fig. 21Primary 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] | . |
| [2] | . |
| [3] | . |
| [4] | . |
| [5] | . |
| [6] | . |
| [7] | . |
| [8] | . |
| [9] | . |
| [10] | . |
| [11] | . |
| [12] | |
| [13] | . |
| [14] | . |
| [15] | |
| [16] | . |
| [17] | . |
| [18] | . |
| [19] | . |
| [20] | . |
| [21] | . |
| [22] | . |
| [23] | . |
| [24] | . |
| [25] | . |
| [26] | . |
| [27] | . |
| [28] | . . |
| [29] | . |
| [30] | . |
| [31] | . |
| [32] | . |
| [33] | . |
| [34] | |
| [35] | . |
| [36] | . |
| [37] | . |
| [38] | . |
| [39] | |
| [40] | . . |
| [41] | . |
| [42] | . |
| [43] | . . |
| [44] | . . |
| [45] | . . |
| [46] | . . |
| [47] | . |
