﻿ 气体动理学统一算法中相容性条件不满足引起的数值误差及其影响研究
«上一篇
 文章快速检索 高级检索

 力学学报  2015, Vol. 47 Issue (1): 163-168  DOI: 10.6052/0459-1879-14-083 0

引用本文 [复制中英文]

[复制中文]
Jiang Dingwu, Mao Meiliang, Li Jin, Deng Xiaogang. STUDY ON THE NUMERICAL ERROR INTRODUCED BY DISSATISFYING THE CONSERVATION CONSTRAINT IN UGKS AND ITS EFFECTS[J]. Chinese Journal of Ship Research, 2015, 47(1): 163-168. DOI: 10.6052/0459-1879-14-083.
[复制英文]

文章历史

2014-06-09收稿
2014-09-12录用
2014-09-17网络版发表

1. 中国空气动力研究与发展中心计算空气动力研究所, 绵阳 621000;
2. 国防科技大学, 长沙 410073

1 UGKS方法介绍 1.1 控制方程及求解过程

 $\begin{gathered} {\text{f}} = \bar f/\left( {{\rho _\infty }/U_\infty ^3} \right),\;{{\text{f}}^ + } = {{\bar f}^ + }/\left( {{\rho _\infty }/U_\infty ^3} \right) \hfill \\ {g_M} = {{\bar g}_M}/\left( {{\rho _\infty }/U_\infty ^3} \right),\;tau = \bar \tau /(L/{U_\infty }) \hfill \\ p = \bar p/({\rho _\infty }U_\infty ^2),\;rho = \bar \rho /{\rho _\infty },\;\mu = \bar \mu /{\mu _\infty } \hfill \\ t = \bar t/(L/{U_\infty }),\lambda = \bar \lambda /\left( {1/U_\infty ^2} \right) \hfill \\ \left( {x,y,z} \right) = \left( {\bar x,\bar y,\bar z} \right)/L,\;(u,v,w) = \left( {\bar u,\bar v,\bar w} \right)/{U_\infty } \hfill \\ \end{gathered}$

 $f_t + uf_x + vf_y = \dfrac{f^ + - f}{\tau }$ (1)

 ${g_M} = \rho {\left( {\frac{\lambda }{\pi }} \right)^{\tfrac{{K + 2}}{2}}}\exp \left( { - \lambda \left( {{{\left( {u - U} \right)}^2} + {{\left( {v - V} \right)}^2} + {\xi ^2}} \right)} \right)$

 ${g^ + } = {g_{\text{M}}}(1 - Pr)c \cdot q({c^2}/RT - 5)/(5pRT)$
${c} = {u} - {U}$是随机速度，$T$为温度，${q}$是热流，$Pr$为普朗特数.

 ${\left( {\rho ,\rho U,\rho V,E} \right)^{\text{T}}} = {\text{ }}\int {{\psi ^{\text{T}}}f{\text{d}}} {\text{ }}{\rm E}$ (2)

 $\begin{gathered} \frac{{\partial Q}}{{\partial t}} + \frac{{\partial F}}{{\partial x}} + \frac{{\partial G}}{{\partial y}} = 0 \hfill \\ F = \int {uf{\psi _\alpha }d{\rm E}} ,\;\;G = \int {uf{\psi _\alpha }d{\rm E}} \hfill \\ \end{gathered}$ (3)

 $\int {(f - {f^ + }){\psi _\alpha }d{\rm E}} {\text{ }} = 0,\;\;\alpha = 1,2,3,4$ (4)

 $\begin{gathered} \Delta Q = - {V^{ - 1}}\int {_{{t^n}}^{{t^{n + 1}}}} [{(J \cdot S)_{i + 1/2,j}} - {(J \cdot S)_{i - 1/2,j}} + \hfill \\ \qquad {(J \cdot S)_{i,j + 1/2}} - {(J \cdot S)_{i,j - 1/2}}]{\text{d}}t \hfill \\ J = Fi + Gj \hfill \\ \end{gathered}$ (5)

 $\begin{gathered} \iiint {{u_n}}{f_p}{(1\;\;u\;\;v\;\;\frac{1}{2}({u^2} + {v^2}))^{\text{T}}}dudvd\xi + \hfill \\ \qquad \iiint {{u_n}}{f_p}{(0\;\;0\;\;0\;\;\frac{1}{2}{\xi ^2})^{\text{T}}}du{\text{d}}vd\xi = \hfill \\ \qquad \iint {\hat g}{(1\;\;u\;\;v\;\;\frac{1}{2}({u^2} + {v^2}))^{\text{T}}}{u_n}dudv + \hfill \\ \qquad \frac{1}{2}\iint (0\;\;0\;\;0\;\;\hat h{)^{\text{T}}}{u_n}dudv \hfill \\ \end{gathered}$

1.2 边界条件

(1) 壁面边界：本文采用动理学壁面边界，详见文献[17].

(2) 超声速来流边界：守恒变量为自由来流值，分布函数按来流变量取平衡态麦克斯韦分布.

(3) 超声速出流边界：守恒变量、分布函数都由内场外推得到.

(4) 对称边界：关于$x$轴对称的情况，$\rho$采用对称条件，$v$采取反对称条件.

2 相容性条件不满足引起的数值误差

 $\int {(f - {f^ + }){\psi _\alpha }{\text{d}}{\rm E}} = Err(N - C)$ (6)
${E r r}$为采用数值积分公式带来的误差. 该误差在速度空间点数增加到一定程度之后将不再减小，这是由数值积分方法本身决定的[18]. 因而与原始方程中隐含的碰撞过程中质量、动量和能量守恒是不相容的.

 $SS = \frac{1}{{\Delta t}}\int {_{{t^n}}^{{t^{n + 1}}}} \left[{\frac{1}{\tau }Err(N - C)} \right]{\text{ d}}t \approx {\left[{\frac{1}{\tau }Err(N - C)} \right]^{n + 1}}$ (7)

 $\tau = \frac{\mu }{{pR{e_\infty }}} \propto \frac{{K{n_\infty }}}{{{M_\infty }}}\frac{\mu }{p}$ (8)

3 守恒型离散速度坐标法

 $\begin{gathered} \iiint {\frac{{{f^ + } - f}}{\tau }}\psi _1^{\text{T}}dudvd\xi = \hfill \\ \qquad \frac{1}{\tau }{\left( {0,0,0,0,- 2/3{q_x},- 2/3{q_y}} \right)^{\text{T}}} \hfill \\ \end{gathered}$ (9)

 $\begin{gathered} \psi _1^{\text{T}} = (1,u,v,\frac{1}{2}\left( {{u^2} + {v^2} + {\xi ^2}} \right),dfrac12\left( {u - U} \right){c^2},\hfill \\ \qquad \frac{1}{2}\left( {v - V} \right){c^2}{)^{\text{T}}} \hfill \\ {c^2} = {\left( {u - U} \right)^2} + {\left( {v - V} \right)^2} + {\xi ^2} \hfill \\ \end{gathered}$

 $\iiint f\psi _1^{\text{T}}{\text{ d}}u{\text{d}}v{\text{ d}}\xi = {\left( {\rho ,\rho U,\rho V,{\text{E}},{q_x},{q_y}} \right)^{\text{T}}}$

 $\begin{gathered} \Sigma {f^ + }\psi _2^{\text{T}} - {\left( {\rho ,\rho U,\rho V,\rho E,{q_x},{q_y}} \right)^{\text{T}}} = \hfill \\ \qquad {\left( {0,0,0,0,- 2/3{q_x},- 2/3{q_y}} \right)^{\text{T}}} \hfill \\ \end{gathered}$ (10)

 $\begin{gathered} \psi _2^{\text{T}} = (1,u,v,\frac{1}{2}\left( {{u^2} + {v^2} + {\xi ^2}} \right),\frac{1}{2}\left( {u - U'} \right){{c'}^2},\hfill \\ \qquad \frac{1}{2}\left( {v - V'} \right){{c'}^2}{)^{\text{T}}} \hfill \\ {{c'}^2} = {\left( {u - U'} \right)^2} + {\left( {v - V'} \right)^2} + {\xi ^2} \hfill \\ \end{gathered}$

4 数值模拟结果与分析

$M_{\infty }=1.8$，$Kn_{\infty }=0.01$,0.40以及$M_{\infty }=10$，$Kn_{\infty }=0.002$,0.050.

 图 1 紧挨壁面第一层网格单元|SS(1)| 分布 Fig.1 $\left| {{SS} (1)} \right|$ on the solid surface

 图 2 阻力中压力项贡献随速度空间网格点数的变化 Fig.2 Drag coefficient vs points in u(v) direction

 图 3 阻力随速度空间网格点数的变化 Fig.3 Drag coefficient vs points in u(v) direction

5 结论

(1)在传统的离散速度坐标法下，相容性条件得不到严格满足，来流马赫数越高，来流努森数越低，数值误差越大，对计算结果影响越大.

(2)采用CDOM之后，相容性条件得到严格满足，尽管单步计算量有所增加，但是能够在较稀的速度空间网格上得到可信的结果. 因而总耗时却能够减少(具体减少多少与来流条件有关)，有助于提高计算效率.

 [1] 费飞, 樊菁. 高雷诺数方腔流动的分子模拟.力学学报, 2013, 45 (5): 653-659 (Fei Fei, Fan Jing. Molecular simulation of driven cavity flows with high Reynolds number. Chinese Journal of Theoretical and Applied Mechanics, 2013, 45(5): 653-659 (in Chinese)) [2] Xu K, Huang JC. A unified gas-kinetic scheme for continuum and rarefied flows. Journal of Computational Physics, 2010, 229(20): 7747-7764 [3] Huang JC, Xu K, Yu PB. A unified gas-kinetic scheme for continuum and rarefied flows Ⅱ: multi-dimensional cases. Communications in Computational Physics, 2012, 12(3): 662-690 [4] Chen SZ, Xu K, Lee C, et al.. A unified gas kinetic scheme with moving mesh and velocity space adaptation. Journal of Computational Physics, 2012, 231(20): 6643-6664 [5] Huang JC, Xu K, Yu PB. A unified gas-kinetic scheme for continuum and rarefied flows Ⅲ: Microflow simulations. Communications in Computational Physics, 2013, 14(5): 1147-1173 [6] 李启兵, 徐昆. 气体动理学格式研究进展. 力学进展, 2012, 42(5): 522-537(Li Qibing, Xu Kun. Progress in gas kinetic scheme. Advances in Mechanics, 2012, 42(5):522-537 (in Chinese)) [7] 江定武, 毛枚良, 邓小刚 等. 二维缝隙跨流域流动数值模拟. 见:计算流体力学研究进展, 第十五届全国计算流体力学会议, 山东烟台, 2012 (Jiang Dingwu, Mao Meiliang, Deng Xiaogang, et al. Numerical simulation of the flow through a slit covering various regimes. In: Progress in Computational Fluid Dynamics, The 15th Computational Fluid Dynamics Conference, Yantai, Shangdong, 2012 (in Chinese)) [8] 李志辉, 张涵信. 稀薄流到连续流的气体运动论统一数值算法初步研究. 空气动力学学报, 2000, 18(3): 255-263 (Li Zhihui, Zhang Hanxin. Study on gas kinetic algorithm for flows from rarefied transition to continuum. Acta Aerodynamica Sinica, 2000, 18(3): 255-263 (in Chinese)) [9] 李志辉, 张涵信. 稀薄流到连续流的气体运动论统一算法研究. 空气动力学学报, 2003, 21(3): 255-266(Li Zhihui, Zhang Hanxin. Study on gas kinetic unified algorithm for flows from rarefied transition to continuum. Acta Aerodynamica Sinica, 2003, 21(3): 255-266((in Chinese)) [10] 李志辉, 吴俊林, 蒋新宇 等. 含转动非平衡效应Boltzmann模型方程统一算法与跨流域绕流问题模拟研究. 空气动力学学报, 2014, 32(2): 137-145 (Li Zhihui, Wu Junlin, Jiang Xinyu, et al.. The unified algorithm for various flow regimes solving boltzmann model equation in rotational non-equilibrium. Acta Aerodynamica Sinica, 2014, 32(2): 137-145 (in Chinese)) [11] 李志辉, 蒋新宇, 吴俊林 等. 转动非平衡玻尔兹曼模型方程统一算法与全流域绕流计算应用. 力学学报, 2014, 46(3): 336-351 (Li Zhihui, Jiang Xinyu, Wu Junlin, et al. Gas-Kinetic unified algorithm for boltzmann model equation in rotational nonequilibrium and its application to the whole range flow regimes. Chinese Journal of Theoretical and Applied Mechanics, 2014, 46(3): 336-351 (in Chinese)) [12] Mieussens L. Discrete velocity models and numerical schemes for the Boltzmann-BGK equation in plane and axisymmetric geometries. Journal of Computational Physics, 2000, 162(2): 429-466 [13] Arslanbekov RR, Kolobov VI, Frolova AA, Kinetic solvers with adaptive mesh in phase space. In:Abe T ed. Rarefied Gas Dynamics, Melville: Amer Inst Physics, 2013: 294-301 [14] Chigullapalli S, Alexeenko A. Unsteady 3D rarefied flow solver based on Boltzmann-ESBGK model kinetic equations. AIAA paper2011-3993, 2011 [15] Huang JC. A conservative discrete ordinate method for model boltzmann equations. Computers & Fluids, 2011, 45(1): 261-267 [16] Shakhov EM. Generalization of the krook kinetic equation. Fluid Dynamics, 1968, 3(5): 95-96 [17] Li QB, Fu S, Xu K. Application of gas kinetic scheme with kinetic boundary conditions in hypersonic flow. AIAA Journal, 2005, 43(10): 2170-2176 [18] 奚梅成, 刘儒勋. 数值分析方法. 合肥:中国科技大学出版社, 2003 (Xi Meicheng, Liu Ruxun. Numerical Analysis Methods. Hefei:University of Science and Technology of China Press, 2003 (in Chinese)) [19] Titarev VA. Conservative numerical methods for model kinetic equations. Computers & Fluids, 2007, 36(9): 1446-1459
STUDY ON THE NUMERICAL ERROR INTRODUCED BY DISSATISFYING THE CONSERVATION CONSTRAINT IN UGKS AND ITS EFFECTS
Jiang Dingwu, Mao Meiliang, Li Jin, Deng Xiaogang
1. Computational Aerodynamics Institute, China Aerodynamics Research and Development Center, Mianyang 621000, China;
2. National University of Defense Technology, Changsha 410073, China
Fund: The project was supported by the National Natural Science Foundation of China (11402287).
Abstract: The Unified Gas Kinetic Scheme (UGKS) was developed by Xu Kun especially for simulating rarefied flows. In UGKS, the conservation constraint can not be satisfied with conventional Discrete Ordinate Method (DOM) and numerical error would be introduced. We found theoretically and numerical experimentally that the error was proportional to the freestream Mach number and inversely proportional to the freestream Knudsen number. Conservative discrete ordinate method (CDOM) was introduced with which the conservation constraint was satisfied at discrete level. The supersonic and hypersonic flows around a cylinder at different Knudsen numbers were simulated and results showed that whether the conservative constraint was satisfied would have a big effect on the variable distributions and the drag coefficients when the Mach number was high or the Knudsen number was low. With CDOM, the grid-free solution could be obtained with a relative coarse velocity grid, resulting in a maximum cost reduction of 2/3.
Key words: Boltzmann model equation    unified scheme    conservation constraint    conservative discrete ordinate method