二维拉氏辐射流体力学人为解构造方法

刘学哲; 林忠; 王瑞利; 余云龙

doi:10.11883/bzycj-2017-0199

二维拉氏辐射流体力学人为解构造方法

doi: 10.11883/bzycj-2017-0199

北京应用物理与计算数学研究所, 北京 100094

基金项目:

国家自然科学基金 11201035

国家自然科学基金 91630312

国防基础科研计划 C1520110002

中国工程物理研究院科学发展基金 2015B0202045

科学挑战专题 TZ2018001

详细信息

作者简介:
刘学哲(1977-), 男, 博士, 副研究员

通讯作者:
王瑞利, wang_ruili@iapcm.ac.cn

中图分类号: O368;O241
计量
- 文章访问数: 4883
- HTML全文浏览量: 1635
- PDF下载量: 34
- 被引次数: 1
出版历程
- 收稿日期: 2017-06-10
- 修回日期: 2017-08-09
- 刊出日期: 2019-01-05

A constructed method of manufactured solutions and code verification for 2D Lagrangian radiation hydrodynamic equations

Institute of Applied Physics and Computational Mathematics, Beijing 100094, China

摘要

摘要: 人为构造解方法是复杂多物理过程耦合程序正确性验证的重要方法之一，适用于二维拉氏大变形网格的流体、辐射耦合人为解模型较为少见。针对拉氏辐射流体力学程序正确性验证的需要，从二维拉氏辐射流体力学方程组出发，基于坐标变换技术，给出了拉氏空间到欧氏空间的物理变量导数关系式，开展了辐射流体耦合的人为解构造方法研究，构造了一类质量方程无源项的二维人为解模型，并应用于非结构拉氏程序LAD2D辐射流体力学计算的正确性考核，为流体运动网格上的辐射扩散计算提供了一种有效手段。数值结果显示观测到的数值模拟收敛阶与理论分析一致。
- 二维拉氏辐射流体力学 /
- 人为解构造方法 /
- 欧氏空间 /
- 拉氏空间 /
- 坐标变换 /
- 程序验证
Abstract: The method of manufactured solutions (MMS) is a fundamental part of code verification for the coupling codes with complex multi-physics process, the MMS for multi-dimension Lagrangian radiation hydrodynamics with large deformation meshes is very sparse yet. In this paper, a new constructed method of manufactured solutions for 2D Lagrangian radiation hydrodynamic equations was proposed based on the derivative relations of physical variables between the Lagrangian space and the Eulerian space. The manufactured solution models without additional source terms in the mass equation could be used to solve diffusion problems on fluid moving meshes and were applied to the verification of the 2D Lagrangian radiation hydrodynamic codes. The numerical results show that the observed order of accuracy matches the formal order of accuracy.
- 2D Lagrangian radiation hydrodynamic equations /
- method of manufactured solutions (MMS) /
- Eulerian space /
- Lagrangian space /
- coordinate transformation /
- code verification

HTML全文

辐射流体力学在武器物理、惯性约束聚变、磁约束聚变、地质学和天体物理等诸多领域有广泛的应用，描述多物理过程耦合的辐射流体力学方程组具有强耦合、强间断、强非线性等特征，很难得到解析解，利用计算机进行数值求解是重要的研究手段。

数值模拟是否正确求解了方程，模拟结果能否再现真实的物理过程，都需要对数值模拟程序进行验证与确认^[1-3]。人为构造解方法^[4]是程序验证的重要方法之一，对某些特殊需求(正确性验证、收敛性考察，等)具有不可替代性。在求解双曲型流体力学方程组程序人为解验证方面已有很好的工作^[5-14]，但大多基于欧氏方程，对抛物型辐射扩散方程解析解及人为解验证方面也有很多工作^[15-16]，而适用于求解辐射与流体耦合方程组及拉氏程序验证的人为解模型尚不多见。

针对拉氏辐射流体力学程序正确性验证的需要，文献[17-18]中构造了一类一维拉氏流体力学和辐射流体力学人为解模型，但一维模型的网格运动仅改变网格点的疏密，不涉及到网格的变形，而网格随流体运动而变形是拉氏计算的特点，因此，构造适应流体大变形的二维拉氏辐射流体力学人为解模型对实际应用程序的正确性验证很有意义。本文中基于二维坐标变换关系式，研究了拉氏辐射流体力学人为解方法，构造了适用于辐射流体力学程序验证的二维人为解模型，并应用于非结构拉氏应用程序LAD2D^[19]辐射流体力学计算的正确性考核，取得了很好的效果。

1. 二维拉氏辐射流体力学方程组

考虑如下拉氏形式的二维辐射流体力学方程组：

$\frac{{{\rm{d}}\rho }}{{{\rm{d}}t}} = - \rho \left( {\frac{{\partial u}}{{\partial x}} + \frac{{\partial v}}{{\partial y}}} \right)$

(1)

$\frac{{{\rm{d}}u}}{{{\rm{d}}t}} = - \frac{1}{\rho }\frac{{\partial p}}{{\partial x}} + {s_u}$

(2)

$\frac{{{\rm{d}}v}}{{{\rm{d}}t}} = - \frac{1}{\rho }\frac{{\partial p}}{{\partial y}} + {s_v}$

(3)

$\frac{{{\rm{d}}e}}{{{\rm{d}}t}} = - \frac{p}{\rho }\left( {\frac{{\partial u}}{{\partial x}} + \frac{{\partial v}}{{\partial y}}} \right) + \frac{1}{\rho }\nabla \cdot \left( {\kappa \nabla T} \right) + {s_e}$

(4)

式中： $\frac{\text{d}}{\text{d}t}=\frac{\partial }{\partial t}+u\frac{\partial }{\partial x}+v\frac{\partial }{\partial y}$ 为拉氏随体导数，(x, y, t)为欧氏时空坐标，ρ为密度，u、v分别为x、y方向的流体速度，s_u、s_v为动量方程源项，s_e为能量方程源项，e为比内能，p为压力，T为温度，κ为导热系数。给定状态方程：

$\begin{array}{*{20}{c}} {p = p\left( {\rho ,T} \right),}&{e = e\left( {\rho ,T} \right),}&{\kappa = \kappa \left( {\rho ,T} \right)} \end{array}$

(5)

和适当的初边值条件，则计算模型封闭。

2. 二维拉氏辐射流体力学人为解构造

2.1 二维坐标变换

引入拉氏坐标x₀=(x₀, y₀), 给定拉氏空间(x₀, y₀, τ)到欧氏空间(x, y, t)的坐标变换关系式：

$x = x\left( {{x_0},{y_0},\tau } \right),\;\;\;\;y = y\left( {{x_0},{y_0},\tau } \right),t = \tau$

(6)

其微分关系为：

${\rm{d}}x = \frac{{\partial x}}{{\partial {x_0}}}{\rm{d}}{x_0} + \frac{{\partial x}}{{\partial {y_0}}}{\rm{d}}{y_0} + \frac{{\partial x}}{{\partial \tau }}{\rm{d}}\tau ,\\{\rm{d}}y = \frac{{\partial y}}{{\partial {x_0}}}{\rm{d}}{x_0} + \frac{{\partial y}}{{\partial {y_0}}}{\rm{d}}{y_0} + \frac{{\partial y}}{{\partial \tau }}{\rm{d}}\tau ,\;\;\;\;\;{\rm{d}}t = {\rm{d}}\tau$

(7)

记J(x₀, y₀, τ)为坐标变换的Jacobi矩阵：

$\mathit{\boldsymbol{J}}\left( {{x_0},{y_0},\tau } \right) = \left( {\begin{array}{*{20}{c}} {{J_{11}}}&{{J_{12}}}\\ {{J_{21}}}&{{J_{22}}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} {\frac{{\partial x\left( {{x_0},{y_0},\tau } \right)}}{{\partial {x_0}}}}&{\frac{{\partial x\left( {{x_0},{y_0},\tau } \right)}}{{\partial {y_0}}}}\\ {\frac{{\partial y\left( {{x_0},{y_0},\tau } \right)}}{{\partial {x_0}}}}&{\frac{{\partial y\left( {{x_0},{y_0},\tau } \right)}}{{\partial {y_0}}}} \end{array}} \right)$

(8)

其行列式记为：

$J = \left| {\begin{array}{*{20}{c}} {\frac{{\partial x}}{{\partial {x_0}}}}&{\frac{{\partial x}}{{\partial {y_0}}}}\\ {\frac{{\partial y}}{{\partial {x_0}}}}&{\frac{{\partial y}}{{\partial {y_0}}}} \end{array}} \right| = \frac{{\partial x}}{{\partial {x_0}}}\frac{{\partial y}}{{\partial {y_0}}} - \frac{{\partial x}}{{\partial {y_0}}}\frac{{\partial y}}{{\partial {x_0}}}$

(9)

$\begin{array}{*{20}{c}} {u\left( {{x_0},{y_0},\tau } \right) = \frac{{\partial x\left( {{x_0},{y_0},\tau } \right)}}{{\partial \tau }},}&{v\left( {{x_0},{y_0},\tau } \right) = \frac{{\partial y\left( {{x_0},{y_0},\tau } \right)}}{{\partial \tau }}} \end{array}$

(10)

为流体速度。

函数u(x₀, y₀, τ)、v(x₀, y₀, τ)与坐标变换的Jacobi矩阵J(x₀, y₀, τ)应满足Cauchy相容关系：

$\begin{array}{*{20}{c}} {\frac{{\partial u}}{{\partial {x_0}}} = \frac{{\partial {J_{11}}}}{{\partial \tau }},}&{\frac{{\partial u}}{{\partial {y_0}}} = \frac{{\partial {J_{12}}}}{{\partial \tau }};}&{\frac{{\partial v}}{{\partial {x_0}}} = \frac{{\partial {J_{21}}}}{{\partial \tau }},}&{\frac{{\partial v}}{{\partial {y_0}}} = \frac{{\partial {J_{22}}}}{{\partial \tau }}} \end{array}$

(11)

下面给出欧氏空间中的空间导数和拉氏随体导数在拉氏空间的表达式。

设物理变量z在欧氏空间关于坐标(x, y, t)的函数关系为f，在拉氏空间关于坐标(x₀, y₀, τ)的函数关系为g，即：

$\begin{array}{*{20}{c}} {z = f\left( {x,y,t} \right) = f\left[ {x\left( {{x_0},{y_0},\tau } \right),y\left( {{x_0},{y_0},\tau } \right),\tau } \right] \equiv }\\ {g\left( {{x_0},{y_0},\tau } \right) = g\left[ {x\left( {{x_0},{y_0},\tau } \right),y\left( {{x_0},{y_0},\tau } \right),t} \right] \equiv f\left( {x,y,t} \right)} \end{array}$

(12)

则有：

$\begin{array}{*{20}{c}} {\frac{{\partial z}}{{\partial {x_0}}} = \frac{{\partial g\left( {{x_0},{y_0},\tau } \right)}}{{\partial {x_0}}} = \frac{{\partial f\left[ {x\left( {{x_0},{y_0},\tau } \right),y\left( {{x_0},{y_0},\tau } \right),\tau } \right]}}{{\partial {x_0}}} = }\\ {\frac{{\partial f\left( {x,y,t} \right)}}{{\partial x}}\frac{{\partial x\left( {{x_0},{y_0},\tau } \right)}}{{\partial {x_0}}} + \frac{{\partial f\left( {x,y,t} \right)}}{{\partial y}}\frac{{\partial y\left( {{x_0},{y_0},\tau } \right)}}{{\partial {x_0}}} = \frac{{\partial z}}{{\partial x}}\frac{{\partial x}}{{\partial {x_0}}} + \frac{{\partial z}}{{\partial y}}\frac{{\partial y}}{{\partial {x_0}}}} \end{array}$

(13)

$\begin{array}{*{20}{c}} {\frac{{\partial z}}{{\partial {y_0}}} = \frac{{\partial g\left( {{x_0},{y_0},\tau } \right)}}{{\partial {y_0}}} = \frac{{\partial f\left[ {x\left( {{x_0},{y_0},\tau } \right),y\left( {{x_0},{y_0},\tau } \right),\tau } \right]}}{{\partial {y_0}}} = }\\ {\frac{{\partial f\left( {x,y,t} \right)}}{{\partial x}}\frac{{\partial x\left( {{x_0},{y_0},\tau } \right)}}{{\partial {y_0}}} + \frac{{\partial f\left( {x,y,t} \right)}}{{\partial y}}\frac{{\partial y\left( {{x_0},{y_0},\tau } \right)}}{{\partial {y_0}}} = \frac{{\partial z}}{{\partial x}}\frac{{\partial x}}{{\partial {y_0}}} + \frac{{\partial z}}{{\partial y}}\frac{{\partial y}}{{\partial {y_0}}}} \end{array}$

(14)

假定坐标变换Jacobi行列式 $J=\frac{\partial x}{\partial {{x}_{0}}}\frac{\partial y}{\partial {{y}_{0}}}-\frac{\partial x}{\partial {{y}_{0}}}\frac{\partial y}{\partial {{x}_{0}}}\ne 0$ ，由式(13)~(14)可以解出欧氏空间导数：

$\begin{array}{*{20}{c}} {\frac{{\partial z}}{{\partial x}} = \frac{1}{J}\left( {\frac{{\partial z}}{{\partial {x_0}}}\frac{{\partial y}}{{\partial {y_0}}} - \frac{{\partial z}}{{\partial {y_0}}}\frac{{\partial y}}{{\partial {x_0}}}} \right),}&{\frac{{\partial z}}{{\partial y}} = \frac{1}{J}\left( {\frac{{\partial z}}{{\partial {y_0}}}\frac{{\partial x}}{{\partial {x_0}}} - \frac{{\partial z}}{{\partial {x_0}}}\frac{{\partial x}}{{\partial {y_0}}}} \right)} \end{array}$

(15)

而:

$\begin{array}{*{20}{c}} {\frac{{\partial z}}{{\partial \tau }} = \frac{{\partial g\left( {{x_0},{y_0},\tau } \right)}}{{\partial \tau }} = \frac{{\partial f\left[ {x\left( {{x_0},{y_0},\tau } \right),y\left( {{x_0},{y_0},\tau } \right),\tau } \right]}}{{\partial \tau }} = }\\ {\frac{{\partial f\left( {x,y,t} \right)}}{{\partial t}}\frac{{\partial f\left( {x,y,t} \right)}}{{\partial x}} + \frac{{\partial x\left( {{x_0},{y_0},\tau } \right)}}{{\partial \tau }} + \frac{{\partial f\left( {x,y,t} \right)}}{{\partial y}}\frac{{\partial y\left( {{x_0},{y_0},\tau } \right)}}{{\partial \tau }} = \frac{{{\rm{d}}f}}{{{\rm{d}}t}} = \frac{{{\rm{d}}z}}{{{\rm{d}}t}}} \end{array}$

即欧氏空间的拉氏随体导数就等于拉氏空间的时间导数：

$\frac{{{\rm{d}}z\left( {x,y,t} \right)}}{{{\rm{d}}t}} = \frac{{\partial z\left( {{x_0},{y_0},\tau } \right)}}{{\partial \tau }}$

(16)

2.2 二维拉氏辐射流体力学人为解构造方法

从拉氏计算中使用的二维拉氏辐射流体力学方程(1)~(5)出发，给定位移解函数和温度解函数，利用二维坐标变换关系式，由无源项的质量方程解出密度解函数，再将已知解函数代入动量方程和能量方程，残余项视为方程源项。具体构造过程如下:

(1) 给定位移解函数x=x(x₀, y₀, τ)，y=y(x₀, y₀, τ)和温度解函数T=T(x₀, y₀, τ)。

(2) 根据式(9)和(10)，由位移解函数x=x(x₀, y₀, τ), y=y(x₀, y₀, τ), 可计算出坐标变换Jacobi行列式J=J(x₀, y₀, τ)和流体速度u=u(x₀, y₀, τ), v=v(x₀, y₀, τ)。

(3) 解连续性方程求出密度解函数。

利用空间导数关系式(15)和Cauchy相容关系式(11)求出速度散度：

$\begin{array}{*{20}{c}} {\nabla \cdot \mathit{\boldsymbol{v}} = \frac{{\partial u}}{{\partial x}} + \frac{{\partial v}}{{\partial y}} = \frac{1}{J}\left( {\frac{{\partial u}}{{\partial {x_0}}}\frac{{\partial y}}{{\partial {y_0}}} - \frac{{\partial u}}{{\partial {y_0}}}\frac{{\partial y}}{{\partial {x_0}}}} \right) + \frac{1}{J}\left( {\frac{{\partial v}}{{\partial {y_0}}}\frac{{\partial x}}{{\partial {x_0}}} - \frac{{\partial v}}{{\partial {x_0}}}\frac{{\partial x}}{{\partial {y_0}}}} \right) = }\\ {\frac{1}{J}\left[ {\frac{\partial }{{\partial \tau }}\left( {\frac{{\partial x}}{{\partial {x_0}}}} \right)\frac{{\partial y}}{{\partial {y_0}}} - \frac{\partial }{{\partial \tau }}\left( {\frac{{\partial x}}{{\partial {y_0}}}} \right)\frac{{\partial y}}{{\partial {x_0}}}} \right] + \frac{1}{J}\left[ {\frac{\partial }{{\partial \tau }}\left( {\frac{{\partial y}}{{\partial {y_0}}}} \right)\frac{{\partial x}}{{\partial {x_0}}} - \frac{\partial }{{\partial \tau }}\left( {\frac{{\partial y}}{{\partial {x_0}}}} \right)\frac{{\partial x}}{{\partial {y_0}}}} \right] = }\\ {\frac{1}{J}\left[ {\frac{\partial }{{\partial \tau }}\left( {\frac{{\partial x}}{{\partial {x_0}}}} \right)\frac{{\partial y}}{{\partial {y_0}}} + \frac{\partial }{{\partial \tau }}\left( {\frac{{\partial y}}{{\partial {y_0}}}} \right)\frac{{\partial x}}{{\partial {x_0}}} - \frac{\partial }{{\partial \tau }}\left( {\frac{{\partial x}}{{\partial {y_0}}}} \right)\frac{{\partial y}}{{\partial {x_0}}} - \frac{\partial }{{\partial \tau }}\left( {\frac{{\partial y}}{{\partial {x_0}}}} \right)\frac{{\partial x}}{{\partial {y_0}}}} \right] = }\\ {\frac{1}{J}\left[ {\frac{\partial }{{\partial \tau }}\left( {\frac{{\partial x}}{{\partial {x_0}}}\frac{{\partial y}}{{\partial {y_0}}} - \frac{{\partial x}}{{\partial {y_0}}}\frac{{\partial y}}{{\partial {x_0}}}} \right)} \right] = \frac{1}{J}\frac{{\partial J}}{{\partial \tau }}} \end{array}$

(17)

代入连续性方程(1)，并利用欧氏空间的拉氏随体导数与拉氏空间的时间导数之间的关系式(16)，得：

$\frac{{\partial \rho \left( {{x_0},{y_0},\tau } \right)}}{{\partial \tau }} = - \rho \left( {{x_0},{y_0},\tau } \right)\frac{1}{{J\left( {{x_0},{y_0},\tau } \right)}}\frac{{\partial J\left( {{x_0},{y_0},\tau } \right)}}{{\partial \tau }}$

即 $J\frac{\partial \rho }{\partial \tau }+\rho \frac{\partial J}{\partial \tau }=0$ ，进一步改写为 $\frac{\partial \left(\rho J \right)}{\partial \tau }=0$ ，可得出密度解函数的精确表达式：

$\rho \left( {{x_0},{y_0},\tau } \right) = \frac{{{\rho _0}\left( {{x_0},{y_0}} \right){J_0}\left( {{x_0},{y_0}} \right)}}{{J\left( {{x_0},{y_0},\tau } \right)}}$

(18)

式中：ρ₀(x₀, y₀)=ρ(x₀, y₀, 0)，J₀(x₀, y₀)=J(x₀, y₀, 0)。

(4) 根据温度解函数确定动量方程源项和内能方程源项。

由动量方程(2)~(3)导出源项：

$\begin{array}{*{20}{c}} {{s_u} = \frac{{{\rm{d}}u}}{{{\rm{d}}t}} + \frac{1}{\rho }\frac{{\partial p}}{{\partial x}},}&{{s_v} = \frac{{{\rm{d}}v}}{{{\rm{d}}t}} + \frac{1}{\rho }\frac{{\partial p}}{{\partial y}}} \end{array}$

(19)

根据时间导数关系式(16)和空间导数关系式(15)，可得：

$\begin{array}{*{20}{c}} {{s_u} = \frac{{\partial u\left( {{x_0},{y_0},\tau } \right)}}{{\partial \tau }} + \frac{1}{{\rho \left( {{x_0},{y_0},\tau } \right)J\left( {{x_0},{y_0},\tau } \right)}} \cdot }\\ {\left[ {\frac{{\partial p\left( {{x_0},{y_0},\tau } \right)}}{{\partial {x_0}}}\frac{{\partial y\left( {{x_0},{y_0},\tau } \right)}}{{\partial {y_0}}} - \frac{{\partial p\left( {{x_0},{y_0},\tau } \right)}}{{\partial {y_0}}}\frac{{\partial y\left( {{x_0},{y_0},\tau } \right)}}{{\partial {x_0}}}} \right]} \end{array}$

(20)

$\begin{array}{*{20}{c}} {{s_v} = \frac{{\partial v\left( {{x_0},{y_0},\tau } \right)}}{{\partial \tau }} + \frac{1}{{\rho \left( {{x_0},{y_0},\tau } \right)J\left( {{x_0},{y_0},\tau } \right)}} \cdot }\\ {\left[ {\frac{{\partial p\left( {{x_0},{y_0},\tau } \right)}}{{\partial {y_0}}}\frac{{\partial x\left( {{x_0},{y_0},\tau } \right)}}{{\partial {x_0}}} - \frac{{\partial p\left( {{x_0},{y_0},\tau } \right)}}{{\partial {x_0}}}\frac{{\partial x\left( {{x_0},{y_0},\tau } \right)}}{{\partial {y_0}}}} \right]} \end{array}$

(21)

由(18)式可得：

$\begin{array}{*{20}{c}} {{s_u} = \frac{{\partial u\left( {{x_0},{y_0},\tau } \right)}}{{\partial \tau }} + \frac{1}{{{\rho _0}\left( {{x_0},{y_0}} \right){J_0}\left( {{x_0},{y_0}} \right)}} \cdot }\\ {\left[ {\frac{{\partial p\left( {{x_0},{y_0},\tau } \right)}}{{\partial {x_0}}}\frac{{\partial y\left( {{x_0},{y_0},\tau } \right)}}{{\partial {y_0}}} - \frac{{\partial p\left( {{x_0},{y_0},\tau } \right)}}{{\partial {y_0}}}\frac{{\partial y\left( {{x_0},{y_0},\tau } \right)}}{{\partial {x_0}}}} \right]} \end{array}$

(22)

$\begin{array}{*{20}{c}} {{s_v} = \frac{{\partial v\left( {{x_0},{y_0},\tau } \right)}}{{\partial \tau }} + \frac{1}{{{\rho _0}\left( {{x_0},{y_0}} \right){J_0}\left( {{x_0},{y_0}} \right)}} \cdot }\\ {\left[ {\frac{{\partial p\left( {{x_0},{y_0},\tau } \right)}}{{\partial {y_0}}}\frac{{\partial x\left( {{x_0},{y_0},\tau } \right)}}{{\partial {x_0}}} - \frac{{\partial p\left( {{x_0},{y_0},\tau } \right)}}{{\partial {x_0}}}\frac{{\partial x\left( {{x_0},{y_0},\tau } \right)}}{{\partial {y_0}}}} \right]} \end{array}$

(23)

其中用到了状态方程关系式(5)：

$p\left( {{x_0},{y_0},\tau } \right) = p\left( {\rho \left( {{x_0},{y_0},\tau } \right),T\left( {{x_0},{y_0},\tau } \right)} \right)$

由内能方程(4)，导出源项：

$\begin{array}{*{20}{c}} {{s_e} = \frac{{{\rm{d}}e}}{{{\rm{d}}t}} + \frac{p}{\rho }\nabla \cdot \mathit{\boldsymbol{v}} - \frac{1}{\rho }\nabla \cdot \left( {\kappa \nabla T} \right) = \frac{{{\rm{d}}e}}{{{\rm{d}}t}} + \frac{p}{\rho }\nabla \cdot \mathit{\boldsymbol{v}} \\- \frac{1}{\rho }\left[ {\frac{\partial }{{\partial x}}\left( {\kappa \frac{{\partial T}}{{\partial x}}} \right) + \frac{\partial }{{\partial y}}\left( {\kappa \frac{{\partial T}}{{\partial y}}} \right)} \right] = }\\ {\frac{{{\rm{d}}e}}{{{\rm{d}}t}} + \frac{p}{\rho }\nabla \cdot \mathit{\boldsymbol{v}} \\- \frac{1}{\rho }\left( {\frac{{\partial \kappa }}{{\partial x}}\frac{{\partial T}}{{\partial x}} + \kappa \frac{{{\partial ^2}T}}{{\partial {x^2}}} + \frac{{\partial \kappa }}{{\partial y}}\frac{{\partial T}}{{\partial y}} + \kappa \frac{{{\partial ^2}T}}{{\partial {y^2}}}} \right)} \end{array}$

(24)

利用时间导数关系式(16)、速度散度表达式(17)和空间导数关系式(15)，得：

$\begin{array}{*{20}{c}} {{s_e} = \frac{{\partial e}}{{\partial \tau }} + \frac{p}{{\rho J}}\frac{{\partial J}}{{\partial \tau }} - \frac{1}{{\rho J}}\left[ {\left( {\frac{{\partial \kappa }}{{\partial {y_0}}}\frac{{\partial x}}{{\partial {x_0}}} - \frac{{\partial \kappa }}{{\partial {x_0}}}\frac{{\partial x}}{{\partial {y_0}}}} \right)\left( {\frac{{\partial T}}{{\partial {x_0}}}\frac{{\partial y}}{{\partial {y_0}}} - \frac{{\partial T}}{{\partial {y_0}}}\frac{{\partial y}}{{\partial {x_0}}}} \right)} \right] - }\\ {\frac{\kappa }{{\rho J}}\left\{ {\frac{\partial }{{\partial {x_0}}}\left[ {\frac{1}{J}\left( {\frac{{\partial T}}{{\partial {x_0}}}\frac{{\partial y}}{{\partial {y_0}}} - \frac{{\partial T}}{{\partial {y_0}}}\frac{{\partial y}}{{\partial {x_0}}}} \right)} \right]\frac{{\partial y}}{{\partial {y_0}}} - \frac{\partial }{{\partial {y_0}}}\left[ {\frac{1}{J}\left( {\frac{{\partial T}}{{\partial {x_0}}}\frac{{\partial y}}{{\partial {y_0}}} - \frac{{\partial T}}{{\partial {y_0}}}\frac{{\partial y}}{{\partial {x_0}}}} \right)} \right]\frac{{\partial y}}{{\partial {x_0}}}} \right\} - }\\ {\frac{1}{{\rho J}}\left[ {\left( {\frac{{\partial \kappa }}{{\partial {y_0}}}\frac{{\partial x}}{{\partial {x_0}}} - \frac{{\partial \kappa }}{{\partial {x_0}}}\frac{{\partial x}}{{\partial {y_0}}}} \right)\left( {\frac{{\partial T}}{{\partial {x_0}}}\frac{{\partial y}}{{\partial {y_0}}} - \frac{{\partial T}}{{\partial {y_0}}}\frac{{\partial y}}{{\partial {x_0}}}} \right)} \right] - }\\ {\frac{\kappa }{{\rho J}}\left\{ {\frac{\partial }{{\partial {x_0}}}\left[ {\frac{1}{J}\left( {\frac{{\partial T}}{{\partial {x_0}}}\frac{{\partial y}}{{\partial {y_0}}} - \frac{{\partial T}}{{\partial {y_0}}}\frac{{\partial y}}{{\partial {x_0}}}} \right)} \right]\frac{{\partial y}}{{\partial {y_0}}} - \frac{\partial }{{\partial {y_0}}}\left[ {\frac{1}{J}\left( {\frac{{\partial T}}{{\partial {x_0}}}\frac{{\partial y}}{{\partial {y_0}}} - \frac{{\partial T}}{{\partial {y_0}}}\frac{{\partial y}}{{\partial {x_0}}}} \right)} \right]\frac{{\partial y}}{{\partial {x_0}}}} \right\}} \end{array}$

(25)

其中用到了状态方程关系式(5)：

$\begin{array}{*{20}{c}} {e\left( {{x_0},{y_0},\tau } \right) = e\left( {\rho \left( {{x_0},{y_0},\tau } \right),T\left( {{x_0},{y_0},\tau } \right)} \right),}&\\{\kappa \left( {{x_0},{y_0},\tau } \right) = \kappa \left( {\rho \left( {{x_0},{y_0},\tau } \right),T\left( {{x_0},{y_0},\tau } \right)} \right)} \end{array}$

3. 二维拉氏辐射流体力学程序验证

3.1 拉氏辐射流体二维人为解模型

为方便计，在不致引起混淆的情况下，时间变量统一用t表示。

由上节的人为解方法可知，选择不同的位移解函数和温度解函数，可以构造出不同的人为解模型。

给定计算区域:

$\begin{array}{*{20}{c}} {0 \le {x_0},{y_0} \le 1;}&{0 \le x,y \le 1;}&{0 \le t \le 1} \end{array}$

和状态方程:

$\begin{array}{*{20}{c}} {p = 0,\;\;e = {c_V}T,\;\;\kappa = {k_0}}&&{{c_V} = 1,\;\;{\kappa _0} = 1} \end{array}$

选择位移解函数：

$x\left( {{x_0},{y_0},\tau } \right) = {x_0} + {x_0}\left( {1 - {x_0}} \right)b\sin \left( {2{\rm{ \mathsf{ π} }}t} \right)\cos \left( {{\rm{ \mathsf{ π} }}{y_0}} \right),$

$y\left( {{x_0},{y_0},\tau } \right) = {y_0} + {y_0}\left( {1 - {y_0}} \right)b\sin \left( {2{\rm{ \mathsf{ π} }}t} \right)\cos \left( {{\rm{ \mathsf{ π} }}{x_0}} \right)\;\;\;\;\;\;\;\;b = 0.8$

和温度解函数:

$T = T\left( {x,y,t} \right) = {T_1} + \sin \left( {2{\rm{ \mathsf{ π} }}x} \right)\cos \left( {2{\rm{ \mathsf{ π} }}y} \right)\;\;\;\;\;\;{T_1} = 2$

初始密度:

${\rho _0}\left( {{x_0},{y_0}} \right) = {\rho _0} = 1$

则可推导出以下。

(1) 坐标函数初、边值：

$0 \le {x_0} \le 1,\;\;\;\;0 \le {y_0} \le 1$

$\begin{array}{*{20}{c}} {{x_{{\rm{left}}}}\left( {{y_0},t} \right) = 0,\;\;\;\;{x_{{\rm{right}}}}\left( {{y_0},t} \right) = 1,\;\;\;\;\\{y_{{\rm{top}}}}\left( {{x_0},t} \right) = 1,\;\;\;\;\;{y_{{\rm{bottom}}}}\left( {{x_0},t} \right) = 0}&{0 \le t \le 1} \end{array}$

且:

${x_{{\rm{left}}}}\left( {{y_0},{t_0}} \right) = {x_{0{\rm{left}}}} = 0,{x_{{\rm{right}}}}\left( {{y_0},{t_0}} \right) = {x_{{\rm{0right}}}} = 1,\\{y_{{\rm{top}}}}\left( {{x_0},{t_0}} \right) = {y_{{\rm{0top}}}} = 1,{y_{{\rm{bottom}}}}\left( {{x_0},{t_0}} \right) = {y_{{\rm{0bottom}}}} = 0$

式中：x_left、x_right、y_top、y_bottom分别为计算区域的左、右、上、下边界。

(2) 坐标变换Jacobi行列式：

$\begin{array}{*{20}{c}} {J\left( {{x_0},{y_0},t} \right) = \frac{{\partial x}}{{\partial {x_0}}}\frac{{\partial y}}{{\partial {y_0}}} - \frac{{\partial x}}{{\partial {y_0}}}\frac{{\partial y}}{{\partial {x_0}}} = }\\ {\left[ {1 + \left( {1 - 2{x_0}} \right)b\sin \left( {2{\rm{ \mathsf{ π} }}t} \right)\cos \left( {{\rm{ \mathsf{ π} }}{y_0}} \right)} \right]\left[ {1 + \left( {1 - 2{y_0}} \right)b\sin \left( {2{\rm{ \mathsf{ π} }}t} \right)\cos \left( {{\rm{ \mathsf{ π} }}{x_0}} \right)} \right] - }\\ {\left[ {{\rm{ \mathsf{ π} }}{x_0}\left( {1 - {x_0}} \right)b\sin \left( {2{\rm{ \mathsf{ π} }}t} \right)\sin \left( {{\rm{ \mathsf{ π} }}{y_0}} \right)} \right]\left[ {{\rm{ \mathsf{ π} }}{y_0}\left( {1 - {y_0}} \right)b\sin \left( {2{\rm{ \mathsf{ π} }}t} \right)\sin \left( {{\rm{ \mathsf{ π} }}{x_0}} \right)} \right]} \end{array}$

不难验证，当0≤x₀≤1，0≤y₀≤1，0 < b < 1时，J(x₀, y₀, t)>0;当t=0时，J(x₀, y₀, t)=1。

(3) 速度解函数：

$u\left( {{x_0},{y_0},t} \right) = 2{\rm{ \mathsf{ π} }}{x_0}\left( {1 - {x_0}} \right)b\cos \left( {2{\rm{ \mathsf{ π} }}t} \right)\cos \left( {{\rm{ \mathsf{ π} }}{y_0}} \right)$

$v\left( {{x_0},{y_0},t} \right) = 2{\rm{ \mathsf{ π} }}{y_0}\left( {1 - {y_0}} \right)b\cos \left( {2{\rm{ \mathsf{ π} }}t} \right)\cos \left( {{\rm{ \mathsf{ π} }}{x_0}} \right)$

(4) 密度解函数：

$\rho \left( {{x_0},{y_0},t} \right) = \frac{1}{{J\left( {{x_0},{y_0},t} \right)}}$

(5) 压力解函数：

$p\left( {{x_0},{y_0},t} \right) = 0$

(6) 比内能解函数：

$e\left( {{x_0},{y_0},t} \right) = T$

(7) 内能方程源项：

$\begin{array}{*{20}{c}} {{s_e} = \frac{{{\rm{d}}e}}{{{\rm{d}}t}} + \frac{p}{\rho }\nabla \cdot \mathit{\boldsymbol{v}} - \frac{1}{\rho }\nabla \cdot \left( {\kappa \nabla T} \right) = 4{{\rm{ \mathsf{ π} }}^2}{c_V}b\cos \left( {2{\rm{ \mathsf{ π} }}t} \right) \cdot }\\ {\left[ {\cos \left( {2{\rm{ \mathsf{ π} }}x} \right)\cos \left( {2{\rm{ \mathsf{ π} }}y} \right){x_0}\left( {1 - {x_0}} \right)\cos \left( {{\rm{ \mathsf{ π} }}{y_0}} \right) \\- \sin \left( {2{\rm{ \mathsf{ π} }}x} \right)\sin \left( {2{\rm{ \mathsf{ π} }}y} \right){y_0}\left( {1 - {y_0}} \right)\cos \left( {{\rm{ \mathsf{ π} }}{x_0}} \right)} \right] + }\\ {8{{\rm{ \mathsf{ π} }}^2}{k_0}\sin \left( {2{\rm{ \mathsf{ π} }}x} \right)\cos \left( {2{\rm{ \mathsf{ π} }}y} \right)J} \end{array}$

(8) 动量方程源项：

${s_u} = \frac{{{\rm{d}}u}}{{{\rm{d}}t}} + \frac{1}{\rho }\frac{{\partial p}}{{\partial x}} = \frac{{\partial u\left( {{x_0},{y_0},\tau } \right)}}{{\partial t}} = - 4{{\rm{ \mathsf{ π} }}^2}{x_0}\left( {1 - {x_0}} \right)b\sin \left( {2{\rm{ \mathsf{ π} }}t} \right)\cos \left( {{\rm{ \mathsf{ π} }}{y_0}} \right)$

${s_v} = \frac{{{\rm{d}}v}}{{{\rm{d}}t}} + \frac{1}{\rho }\frac{{\partial p}}{{\partial y}} = \frac{{\partial v\left( {{x_0},{y_0},\tau } \right)}}{{\partial t}} = - 4{{\rm{ \mathsf{ π} }}^2}{y_0}\left( {1 - {y_0}} \right)b\sin \left( {2{\rm{ \mathsf{ π} }}t} \right)\cos \left( {{\rm{ \mathsf{ π} }}{x_0}} \right)$

(9) 边界条件。

流体为固壁边界条件：

$\begin{array}{*{20}{c}} {{u_{{\rm{left}}}}\left( {{y_0},t} \right) = 0,}&{{u_{{\rm{right}}}}\left( {{y_0},t} \right) = 0,}&{{v_{{\rm{top}}}}\left( {{x_0},t} \right) = 0,}&{{v_{{\rm{bottom}}}}\left( {{x_0},t} \right) = 0} \end{array}$

辐射左右边界为第一类边界条件，上下边界为第二类边界条件：

$\begin{array}{*{20}{c}} {{T_{{\rm{left}}}}\left( {{y_0},t} \right) = {T_1},}&{{T_{{\rm{right}}}}\left( {{y_0},t} \right) = {T_1},}&{{{\frac{{\partial T}}{{\partial y}}}_{{\rm{top}}}}\left( {{x_0},t} \right) = 0,}&{{{\frac{{\partial T}}{{\partial y}}}_{{\rm{bottom}}}}\left( {{x_0},t} \right) = 0} \end{array}$

不难看出，构造的人为解满足保正的物理条件，即对任意的:

$\left( {{x_0},{y_0},t} \right) \in \left( {{x_{{\rm{left}}}},{x_{{\rm{right}}}}} \right) \times \left( {{y_{{\rm{top}}}},{y_{{\rm{bottom}}}}} \right) \times {{\bf{R}}^ + }$

都有:

$J\left( {{x_0},{y_0},t} \right) > 0,\;\;\;\;\rho \left( {{x_0},{y_0},t} \right) > 0,\;\;\;\;T\left( {{x_0},{y_0},t} \right) > 0,\;\;\;\;e\left( {{x_0},{y_0},t} \right) > 0$

3.2 二维拉氏辐射流体力学LAD2D程序正确性验证

将构造的人为解模型应用于二维非结构网格拉氏辐射流体力学应用程序LAD2D^[19]的正确性验证。

初始时，计算区域均匀划分为32×32的方形网格，人为解中参数b=0.8，取固定时间步长Δt=7.812 5×10^－4，计算到t=1.0时刻。

由构造的人为解位移解函数可知，它为关于时间t的周期函数，这里选取了t=0，0.125，0.25等3个典型时刻，图 1给出了不同时刻流体运动网格图。

图 1 不同时刻流体运动网格

Figure 1. Fluid meshes at different times

下载: 全尺寸图片幻灯片

图 2给出了t=0.125时刻流体网格上的密度、温度、x方向速度、y方向速度的数值解和人为精确解的等值线，红线为数值解，绿线为人为精确解，两者基本一致，验证了程序辐射流体求解的正确性。

图 2 t=0.125时刻流体运动网格上密度、温度、x方向速度、y方向速度等值线(红：数值解绿：精确解)

Figure 2. Contours of density, temperature, x and y components of velocity(red: numerical solution; green: exact solution)

下载: 全尺寸图片幻灯片

通常的扩散方程求解，无论是正交网格还是大变形扭曲网格，求解过程中网格始终保持不变，在二维运动网格上考察扩散问题目前尚未见文献报道。由本文人为解构造可知，通过设定特殊形式的状态方程(p=0)，可以在流体运动网格上求解扩散模型。图 3给出了不同时刻流体运动网格上温度数值解和人为精确解的比较，红线为数值解，绿线为人为精确解，从图中可以看出，从初始时刻开始，网格在不断运动，网格变形程度逐渐增大，而我们构造的温度解函数在欧氏空间来看不随时间变化，数值解较好地保持了精确解的这一性质。

图 3 不同时刻流体运动网格上温度等值线

(红：数值解绿：精确解)

Figure 3. Contours of temperature on fluid meshes at different time

(red: numerical solution; green: exact solution)

下载: 全尺寸图片幻灯片

进一步考察人为解的数值收敛性，在不同的网格规模下，网格比Δt/Δx²=0.8保持不变，计算到t=1.0时刻，表 1给出了温度的误差收敛性分析。数值结果显示温度L2模和最大模误差二阶收敛，验证了计算结果与理论分析的一致性。

表 1 温度收敛误差和收敛阶

Table 1. Convergence errors and orders for temperature

网格数	L2模误差	L2模误差收敛阶	最大模误差	最大模误差收敛阶
8×8	2.82×10^-2		1.04×10^-1
16×16	7.13×10^-3	1.98	3.06×10^-2	1.76
32×32	1.79×10^-3	1.99	7.97×10^-3	1.94
64×64	4.49×10^-4	2.00	2.01×10^-3	1.99
128×128	1.12×10^-4	2.00	5.03×10^-4	2.00

下载: 导出CSV

| 显示表格

4. 结论

从拉氏计算中使用的辐射流体力学方程组出发，基于二维坐标变换技术，研究了一类质量方程无源项，动量方程和能量方程包含源项的人为解构造方法，构造了一类适用于辐射流体力学程序验证的二维人为解模型，并应用于非结构拉氏程序LAD2D辐射流体力学计算的正确性考核，为流体运动网格上的辐射扩散计算提供了新的途径。

图 1 不同时刻流体运动网格

Figure 1. Fluid meshes at different times

下载: 全尺寸图片幻灯片

图 2 t=0.125时刻流体运动网格上密度、温度、x方向速度、y方向速度等值线(红：数值解绿：精确解)

Figure 2. Contours of density, temperature, x and y components of velocity(red: numerical solution; green: exact solution)

下载: 全尺寸图片幻灯片

图 3 不同时刻流体运动网格上温度等值线

(红：数值解绿：精确解)

Figure 3. Contours of temperature on fluid meshes at different time

(red: numerical solution; green: exact solution)

下载: 全尺寸图片幻灯片

表 1 温度收敛误差和收敛阶

Table 1. Convergence errors and orders for temperature

网格数	L2模误差	L2模误差收敛阶	最大模误差	最大模误差收敛阶
8×8	2.82×10^-2		1.04×10^-1
16×16	7.13×10^-3	1.98	3.06×10^-2	1.76
32×32	1.79×10^-3	1.99	7.97×10^-3	1.94
64×64	4.49×10^-4	2.00	2.01×10^-3	1.99
128×128	1.12×10^-4	2.00	5.03×10^-4	2.00

下载: 导出CSV

参考文献(19)

[1]	OBERKAMPF W L, ROY C J. Verification and validation in scientific computing[M]. New York:Cambridge University Press, 2010.
[2]	王瑞利, 江松.多物理耦合非线性偏微分方程与数值解不确定度量化数学方法[J].中国科学:数学, 2015, 45(6):723-738. DOI: 10.1360/N012014-00115. WANG Ruili, JIANG Song. Mathematical methods for uncertainty quantification in nonlinear multi-physics systems and their numerical simulations[J]. Scientia Sinica Mathematica, 2015, 45(6):723-738. DOI: 10.1360/N012014-00115.
[3]	邓小刚, 宗文刚, 张来平, 等.计算流体力学中的验证与确认[J].力学进展, 2007, 37(2):279-288. doi: 10.3321/j.issn:1000-0992.2007.02.011 DENG Xiaogang, ZONG Wengang, ZHANG Laiping, et al. Verification and validation in computational fluid dynamics[J]. Advances in Mechanics, 2007, 37(2):279-288. doi: 10.3321/j.issn:1000-0992.2007.02.011
[4]	SALARI K, KNUPP P. Code verification by the method of manufactured solutions: SAND200-1444[R]. Sandia National Laboratories, 2000.
[5]	ROY C J, NELSON C C, SMITH T M, et al. Verification of Euler/Navier-Stokes codes using the method of manufactured solutions[J]. International Journal for Numerical Methods in Fluids, 2004, 44:599-620. DOI: 10.1002/fld.660.
[6]	王瑞利, 刘全, 刘希强, 等.人为解方法及其在流体力学程序验证中的应用[J].计算机应用与软件, 2012, 29(11):4-7. DOI: 10.3969/j.issn.1000-386x.2012.11.002. WANG Ruili, LIU Quan, LIU Xiqiang, et al. Artificial solution and its applications in verification of hydrodynamics program[J]. Computer Applications and Software, 2012, 29(11):4-7. DOI: 10.3969/j.issn.1000-386x.2012.11.002.
[7]	刘全, 王瑞利, 刘希强, 等.流体力学方程组一类人为解构造方法[J].数学的实践与认识, 2013, 43(8):176-182. doi: 10.3969/j.issn.1000-0984.2013.08.024 LIU Quan, WANG Ruili, LIU Xiqiang, et al. One method of artificial solution to 2D hydrodynamics euler equation[J]. Mathematics in Practice and Theory, 2013, 43(8):176-182. doi: 10.3969/j.issn.1000-0984.2013.08.024
[8]	刘希强, 王瑞利, 刘全.流体力学方程组人为解的构造及其应用[J].聊城大学学报(自然科学版), 2013, 26(2):1-3. doi: 10.3969/j.issn.1672-6634.2013.02.001 LIU Xiqiang, WANG Ruili, LIU Quan. Manufactured solution of fluid mechanics equation of construction and its application[J]. Journal of Liaocheng University (Natural Science Edition), 2013, 26(2):1-3. doi: 10.3969/j.issn.1672-6634.2013.02.001
[9]	杨振虎.CFD程序验证的虚构解方法及其边界精度匹配问题[J].航空计算技术, 2007, 37(6):5-9. doi: 10.3969/j.issn.1671-654X.2007.06.002 YANG Zhenhu. CFD code verification via the method of manufactured solution and its boundary accuracy match problem[J]. Aeronautical Computing Technique, 2007, 37(6):5-9. doi: 10.3969/j.issn.1671-654X.2007.06.002
[10]	刘全, 王瑞利, 林忠, 等.流体力学拉氏程序收敛性及数值计算不确定度初探[J].计算物理, 2013, 30(3):346-352. doi: 10.3969/j.issn.1001-246X.2013.03.004 LIU Quan, WANG Ruili, LIN Zhong, et al. Asymptotic convergence analysis and quantification of uncertainty in Lagrangian computations[J]. Chinese Journal of Computational Physics, 2013, 30(3):346-352. doi: 10.3969/j.issn.1001-246X.2013.03.004
[11]	刘全, 王瑞利, 林忠, 等.爆轰计算JWL状态方程参数不确定度研究[J].爆炸与冲击, 2013, 33(6):647-654. doi: 10.3969/j.issn.1001-1455.2013.06.014 LIU Quan, WANG Ruili, LIN Zhong, et al. Uncertainty quantification for JWL EOS parameters in explosive numerical simulation[J]. Explosion and Shock Waves, 2013, 33(6):647-654. doi: 10.3969/j.issn.1001-1455.2013.06.014
[12]	王瑞利, 刘全, 温万治.非嵌入式多项式混沌法在爆轰产物JWL参数评估中的应用[J].爆炸与冲击, 2015, 35(1):9-15. http://www.bzycj.cn/CN/abstract/abstract9419.shtml WANG Ruili, LIU Quan, WEN Wanzhi. Non-intrusive polynomial chaos methods and its application in the parameters assessment of explosion product JWL[J]. Explosion and Shock Waves, 2015, 35(1):9-15. http://www.bzycj.cn/CN/abstract/abstract9419.shtml
[13]	梁霄, 王瑞利.混合不确定度量化方法及其在计算流体动力学迎风格式中的应用[J].爆炸与冲击, 2016, 36(4):509-515. http://www.bzycj.cn/CN/abstract/abstract9625.shtml LIANG Xiao, WANG Ruili. Mixed uncertainty quantification and its application in upwind scheme for computational fluid dynamics(CFD)[J]. Explosion and Shock Waves, 2016, 36(4):509-515. http://www.bzycj.cn/CN/abstract/abstract9625.shtml
[14]	王瑞利, 梁霄.基于误差马尾图量化爆轰数值模拟结果的置信度[J].爆炸与冲击, 2017, 37(6):893-900. http://www.bzycj.cn/CN/abstract/abstract9811.shtml WANG Ruili, LIANG Xiao. Confidence level of numerical simulation of detonation through quantifying the horsetail of errors[J]. Explosion and Shock Waves, 2017, 37(6):893-900. http://www.bzycj.cn/CN/abstract/abstract9811.shtml
[15]	杨容, 杭旭登, 李敬宏.二维柱对称辐射输运基准模型及程序考核[J].计算物理, 2010, 27(4):533-540. doi: 10.3969/j.issn.1001-246X.2010.04.007 YANG Rong, HANG Xudeng, LI Jinghong. A two-dimensional cylindric symmetric radiative transfer benchmark model and code tests[J]. Chinese Journal of Computational Physics, 2010, 27(4):533-540. doi: 10.3969/j.issn.1001-246X.2010.04.007
[16]	刘学哲, 余云龙, 王瑞利, 等.非结构任意多边形网格辐射扩散方程有限体积格式[J].数值计算与计算机应用, 2010, 31(4):259-270. http://d.old.wanfangdata.com.cn/Periodical/szjsyjsjyy201004003 LIU Xuezhe, YU Yunlong, WANG Ruili, et al. A cell-centered finite volume scheme for discretizing diffusion equation[J]. Journal on Numerical Methods and Computer Applications, 2010, 31(4):259-270. http://d.old.wanfangdata.com.cn/Periodical/szjsyjsjyy201004003
[17]	刘全, 林忠, 王瑞利, 等.流体力学拉氏方程组一类人为解构造方法[J].河北师范大学学报(自然科学版), 2014, 38(5):451-455. DOI: 10.11826/j.issn.1000-5854.2014.05.004. LIU Quan, LIN Zhong, WANG Ruili, et al. One method of manufactured solution to 2D hydrodynamics Euler equations[J]. Journal of Hebei Normal University (Natural Science Edition), 2014, 38(5):451-455. DOI: 10.11826/j.issn.1000-5854.2014.05.004.
[18]	余云龙, 林忠, 王瑞利, 等.辐射流体力学Lagrange方程组一类人为解构造方法[J].应用数学与力学, 2015, 36(1):110-118. DOI: 10.3879/j.issn.1000-0887.2015.01.010. YU Yunlong, LIN Zhong, WANG Ruili, et al. A method of manufactured solution for verification of Lagrangian radiation hydrodynamic codes[J]. Applied Mathematics and Mechanics, 2015, 36(1):110-118. DOI: 10.3879/j.issn.1000-0887.2015.01.010.
[19]	王瑞利, 林忠, 温万治, 等.自主多介质拉氏自适应流体动力学软件LAD2D研制及其应用[J].计算机辅助工程, 2014, 23(2):1-7. DOI:10.13340/j.cae. 2014.02.001. WANG Ruili, LIN Zhong, WEN Wanzhi, et al. Development and application of adaptive multi-media Lagrangian fluid dynamics software LAD2D[J]. Computer Aided Engineering, 2014, 23(2):1-7. DOI:10.13340/j.cae. 2014.02.001.