三维粘性流的数值模拟

1、2001.2.8三维粘性流的数值模拟1.控制方程三维非定常可压缩NS方程在直角坐标下的形式为（1）这里式中 （2）为总能量，粘性系数 由Sutherland公式给出 （3）这里的 ，。而 为热传导系数， 是Prandtl数。粘性应力张量另外还要给出状态方程 （5）2.无量纲化取来流的密度 、速度 、温度 、黏性系数 和飞行器的特征长度 为特征量，定义无量纲量利用这些无量纲量，可将方程组（1）无量纲化，得到（6）这里，在省略无量纲量上面的“”后，无粘通量 、，总能量 和粘性应力张量等表达式的形式不变，定压比热从而粘性通量无量纲化后（省略无量纲量上面的“”）成为Sutherland公式本身就是无量

2、纲形式的，而状态方程的无量纲形式为省略无量纲量上面的“”，就是 （7）以上的推导用以到下列关系式：Mayer公式 ，音速 ，来流马赫数 ，以及来流的状态方程 。现在来改写粘性通量。为简化下面的推导，记 ， 定义向量和矩阵则方程组（6）中的粘性通量可写成（8）3.坐标变换通过坐标变换（9）可将方程组（6）变换到计算域，形式为（10）式中这里 （11）而是坐标变换的Jacobi行列式。方程组（10）的推导过程如下。用坐标变换的Jacobi行列式 乘以方程组（6）的两边，得（12）利用求导的链式规则，对上式左边的无粘通量，有但是由（11）式同理就有同样地，对（12）式右边的粘性通量，也有注意到坐标变

3、换（9）与时间无关， 可直接进入（12）式左边的时间导数项，就得到方程组（10）。4. 粘性通量的简化记 ， 则由（8）式，粘性通量又可写成 （13）定义矩阵 ，则 ，可用消去法求得。先求得再求逆矩阵矩阵中的最后一行是由（2）式（已无量纲化）得到的。即其中 。所以（13）式中的下面开始对粘性通量（13）式进行一系列的简化。首先，类似于薄层NS方程，在（13）式对 求和的三项中只保留 这一项；其次，对（13）式括号中的双重求和，只保留 的三项。记（14）则粘性通量就简化成为注意到简化过程中保留下来的矩阵 和矩阵 都是下三角矩阵，因而（14）式中导数前面的整个系数矩阵仍是下三角矩阵，其特征值就是其

4、对角线元素，即设 ，记 ，求得该矩阵的谱半径为再用该谱半径近似代替（14）式中的系数矩阵，最终将粘性通量简化成（15）以上对粘性通量的简化仅用于差分格式的隐式部分。5. 无粘通量的Jacobi矩阵方程组（10）中无粘通量的Jacobi矩阵定义为记向量和矩阵则无粘通量而无粘通量的Jacobi矩阵就是若记则 于是向量 可写成用 的各分量可将 的各分量表成这里用到了从（2）式导出的表达式通过计算，得其中 是焓, 是音速。下面是计算过程。为方便起见，改记 ，则 ， ，首先，对 ， ，所以而对 ，借助Dirac函数 ，有最后， ，所以而所以 ，而6. 三维无粘气体力学方程组的非守恒形式我们不直接计算矩阵

5、 的特征值和特征向量，转而考虑非守恒形式的三维气体力学方程组（已无量纲化）（16）利用向量和矩阵则方程组（16）也可以写成向量形式（17）而方程组（6）在暂时略去粘性项之后写成即或与（17）式比较，可知定义矩阵则 （18）此式表明矩阵 与矩阵 相似。因而，要计算矩阵 的特征值，只需矩阵 的特征值。而要计算矩阵 的特征向量，可先计算矩阵 的特征向量。下面计算矩阵 的特征值。考虑行列式方程该行列式方程的根即为矩阵 的，也就是矩阵 的，特征值，这里 。设 、 是与矩阵 的各特征值 、 、 分别对应的特征向量，则总可以使它们成为一组线性无关的向量。以它们为列构成矩阵 ，则 为可逆矩阵。同理，以矩阵 的

6、一组线性无关的特征向量 、 为列也构成一个可逆矩阵 。记对角矩阵 ，则有 ， ，即 ，从而由（18）式，有此式表明矩阵 与 之间有关系式 ， 因此，只需求出矩阵 及其逆矩阵 ，并利用矩阵 及其逆矩阵 ，即可求得矩阵 及其逆矩阵 。欲求矩阵 的特征向量，需求方程组 的非零解。而欲求逆矩阵 ，可求 的转置矩阵 的特征向量，也就是方程组 的非零解。设 、 是矩阵 的特征向量，则也可以使它们成为一组线性无关的向量，并满足正交条件如果 不是重特征值，上述正交条件自动满足。再对 、 进行归一化，即用 除向量 ，还可以有 。于是，以这组特征向量为行，即构成逆矩阵 。对特征值 ，方程组 成为即或这里 。取 ，

7、即得到非零解 ，也就是分别对应于特征值 和 的特征向量。对特征值 ，方程组 成为即或按归一化的要求，取 ，即得到分别对应于特征值 和 的特征向量 。对特征值 ，方程组 成为即 任意，在有些文献中取下面三个向量作为其非零解但这组解中，当 时前两个成为线性相关的，当 或 时甚至有一个成为零向量。这样一来，矩阵 就不是可逆矩阵了。因此这一做法是有问题的。这里取因 、 不能全为零，这组向量总是线性无关的非零向量。于是对特征值 ，方程组 成为即所求之特征向量 除需满足这一方程组外，还需满足正交条件由这些关系式，有即 ，也就是 。此时上述正交条件成为按归一化的要求，取 ，即得到特征向量类似地，可求得最终得

8、到对任一向量 ,记7.隐式格式及LU分解对方程组（10）使用隐式格式求解时，差分格式的形式为（19）这里， 是流动参量的时间增量。将（19）式改写成（20）此式的左边，即格式的隐式部分，其无粘通量可进行线性化处理。因 ，故 ，从而有得到具有二阶时间精度的格式（21）至于粘性通量，因 ，同理， ，所以由（13）式，有而由（15）式 因而如果只要求格式具有一阶时间精度，则上述两个同样具有一阶时间精度的量可以互相代替，从而将格式（21）进一步简化为（22）取 ，其中 表示矩阵的谱半径，并记则显然有 ，于是有根据迎风原则对上式左边各项分别进行向后和向前差分，就得到式中，下标 被省略， 是 在 处的值，

9、 是 在 处的值， 是 在 处的值。整理，得（23）记就有 （24）由近似LU分解即得到隐式格式的LU分解算法 （25）它可用下述两步格式求解即 （26）注意到 ，则是对角矩阵，因而在按对角线扫描时，算法（26）可完全避免矩阵求逆。至于右端项 的计算，对其中的粘性通量可简单地使用中心差分离散，但需使用未简化的形式。例如这里 ， ，余类推。右端项中的无粘通量的计算，可使用MUSCL处理。即这里 是数值通量。以 为例，表达式为（27）这里 ， 是Jacobi矩阵 的特征向量矩阵，向量 的分量定义为（28）式中 为矩阵 的特征值， 是向量 的分量。下标为 的量，可由 、 通过ROE平均或算术平均计算。、 的定义为（29）式中 ，取 。函数 （30）对含激波的流场，为了提高激波分辨率，