6.3 Methodologies for Unsteady Flows 非定常流动的求解方法[cfd-6-3]

在许多工程学科中ï¼非定常流动现象的模拟正变得日益重要。例子包括:涡轮机械中静止部件与旋转部件的相互作用、活塞式发动机、流固耦合、直升机空气动力学、气动声学、湍流的DNS或LES、爆轰ï¼等等。显然ï¼模拟必须高效地进行ï¼并且精度要与所求解的问题相称。en

The simulation of unsteady flow phenomena is becoming increasingly important in many engineering disciplines. Examples are the interaction between stationary and rotating parts in turbomachinery, piston engines, fluid-structure interaction, helicopter aerodynamics, aeroacoustics, DNS or LES of turbulent flows, detonations, etc. Clearly, the simulation has to be conducted efficiently and with accuracy adequate for the problem being solved.

当时间尺度与“空间尺度除以特征值”相当时ï¼即物理决定的CFL数为1的量级时ï¼显式格式对某些非定常应用是最好的选择。气动声学、DNS和LES即属此类。对此类应用ï¼显式Runge-Kutta格式相当流行。由于这些应用中全局物理现象的演化远慢于解的局部变化ï¼必须在很长的物理时间内积分。为准确做到这一点ï¼显式格式的时间分辨率须达三阶或更高。这就需要使用与6.1节所介绍不同的Runge-Kutta方法(参见[95]、[96])。en

Explicit schemes represent the best choice for certain unsteady applications when the time scales are comparable to the spatial scales over the eigenvalue, i.e., when the CFL number dictated by the physics is of the order of unity. This is for example the case in aeroacoustics, DNS, and LES. For such applications, explicit Runge-Kutta schemes are quite popular. Since the global physical phenomena evolve much slower than the solution changes locally in these applications, it is necessary to integrate over a long period of physical time. In order to do this accurately, the temporal resolution of the explicit scheme has to be of 3rd or higher order. This requires the use of Runge-Kutta methods different to that we presented in Section 6.1 (see, e.g., [95], [96]).

在另一些情形ï¼物理时间尺度相比“空间尺度除以特征值”要大(例如颤振、转子-静子干涉等),CFL数可以取几百甚至几千的量级而不损害模拟精度。显然ï¼这类情形下隐式格式更为合适。下面我们将讨论一种被称为双时间步进(dual time-stepping)的特定技术ï¼它经常用于非定常流动。en

In other cases, where the physical time scales are large in comparison to the spatial scales divided by the eigenvalue (e.g., flutter, rotor-stator interaction, etc.), the CFL number can be chosen in the order of several hundreds or even thousands without impairing the accuracy of the simulation. Obviously, in such cases an implicit scheme is more appropriate. In the following, we shall discuss a particular technique known as the dual time-stepping approach, which is very often employed for unsteady flows.

双时间步进方法基于式(6.2)基本非线性格式的二阶时间精度版本。为此ï¼在式(6.2)中取 $\beta = 1$ 、 $\omega = 1/2$ 。于是得到en

The dual time-stepping approach is based on the second-order time accurate version of the basic non-linear scheme in Eq. (6.2). For this purpose we set $\beta = 1$ and $\omega = 1/2$ in Eq. (6.2). Hence, we obtain

\frac{3\left(\Omega\bar{M}\right)^{n+1}_{I}\vec{W}^{n+1}_{I} - 4\left(\Omega\bar{M}\right)^{n}_{I}\vec{W}^{n}_{I} + \left(\Omega\bar{M}\right)^{n-1}_{I}\vec{W}^{n-1}_{I}}{2\Delta t} = -\vec{R}^{n+1}_{I}, \tag{6.67}

其中 $\Delta t$ 表示全局物理时间步长, $\bar{M}$ 表示质量矩阵。方程(6.67)构成对式(6.1)中时间导数的三点后向差分(后向Euler)近似。为了求解式(6.67)给出的非线性方程组ï¼可以使用Newton法或时间推进方法。后者可写为en

where $\Delta t$ denotes the global physical time step and $\bar{M}$ the mass matrix, respectively. Equation (6.67) constitutes a 3-point backward-difference (backward Euler) approximation of the time derivative in Eq. (6.1). In order to solve the system of non-linear equations given by Eq. (6.67), we can use either Newton's method or a time-stepping methodology. The latter can be written as

\frac{\partial}{\partial t^{*}}\left(\Omega^{n+1}_{I}\vec{W}^{*}_{I}\right) = -\vec{R}^{*}_{I}\left(\vec{W}^{*}\right), \tag{6.68}

其中 $\vec{W}^{*}$ 是 $\vec{W}^{n+1}$ 的近似, $t^{*}$ 表示伪时间变量。注意ï¼时间导数中没有质量矩阵。非定常残差定义为en

where $\vec{W}^{*}$ is the approximation to $\vec{W}^{n+1}$ and $t^{*}$ denotes a pseudo-time variable. Note that there is no mass matrix in the time derivative. The unsteady residual is defined as

\vec{R}^{*}_{I}\left(\vec{W}^{*}\right) = \vec{R}_{I}\left(\vec{W}^{*}\right) + \frac{3}{2\Delta t}\left(\Omega\bar{M}\right)^{n+1}_{I}\vec{W}^{*}_{I} - \vec{Q}^{*}_{I}. \tag{6.69}

式(6.68)中在时间推进期间保持不变的所有项都归入一个源项ï¼即en

All terms which are constant during the time-stepping in Eq. (6.68) are gathered in a source term, i.e.,

\vec{Q}^{*}_{I} = \frac{2}{\Delta t}\left(\Omega\bar{M}\right)^{n}_{I}\vec{W}^{n}_{I} - \frac{1}{2\Delta t}\left(\Omega\bar{M}\right)^{n-1}_{I}\vec{W}^{n-1}_{I}. \tag{6.70}

在网格运动和/或变形的情形ï¼控制体的新尺寸ï¼即式(6.68)中的 $\Omega^{n+1}$ ï¼必须满足几何守恒定律(详见附录A.5及其所引文献)。en

In the case of moving and/or deforming grids, the new size of the control volume, i.e., $\Omega^{n+1}$ in Eq. (6.68) has to satisfy the Geometry Conservation Law (see Appendix A.5 for details and references).

式(6.68)的定常解对应于新时间层上的流动变量ï¼即 $\vec{W}^{*} = \vec{W}^{n+1}$ 。由于伪时间内定态时 $\vec{R}^{*}_{I} = 0$ ï¼方程(6.67)便得到满足。此前介绍过的任何显式或隐式时间推进格式都可以用来在伪时间内求解方程组(6.68)。下面我们将讨论双时间步进方法在显式多级格式与隐式格式中的实现。en

The stationary solution of Eq. (6.68) corresponds to the flow variables at the new time level, i.e., $\vec{W}^{*} = \vec{W}^{n+1}$ . Since $\vec{R}^{*}_{I} = 0$ at steady state in pseudo time, Equation (6.67) is fulfilled. Any of the previously presented explicit or implicit time-marching schemes can be employed for the solution of the system of equations (6.68) in the pseudo time. In the following, we shall discuss the implementation of the dual time-stepping approach for explicit multistage and implicit schemes.

6.3.1 Dual Time-Stepping for Explicit Multistage Schemes 显式多级格式的双时间步进[cfd-6-3-1]

Jameson[97]首先用由局部时间步进与多重网格加速的显式多级格式实现了双时间方法。该方法的重要优点是物理时间步长不像显式方法通常那样受限制ï¼可以完全依据流动物理来选取。另一方面ï¼只需为源项 $\vec{Q}^{*}$ 增加额外存储ï¼这使得该方法非常有吸引力。求解伪时间问题(6.68)的 $$m$$ 级显式格式为en

Jameson [97] first implemented the dual-time methodology using an explicit multistage scheme accelerated by local time-stepping and multigrid. The significant advantage of this approach is that the physical time step is not restricted as usual in explicit methods. It can be chosen based solely on the flow physics. On the other hand, additional storage is needed only for the source term $\vec{Q}^{*}$ , which makes the approach very attractive. An $$m$$ -stage explicit scheme for the solution of the pseudo-time problem (6.68) reads

\begin{aligned} \vec{W}^{(0)}_{I} &= \left(\vec{W}^{*}_{I}\right)^{l}\\ \vec{W}^{(1)}_{I} &= \vec{W}^{(0)}_{I} - \frac{\alpha_{1}\Delta t^{*}_{I}}{\Omega^{n+1}_{I}}\vec{R}^{*}_{I}\left(\vec{W}^{(0)}\right)\\ \vec{W}^{(2)}_{I} &= \vec{W}^{(0)}_{I} - \frac{\alpha_{2}\Delta t^{*}_{I}}{\Omega^{n+1}_{I}}\vec{R}^{*}_{I}\left(\vec{W}^{(1)}\right)\\ &\ \ \vdots\\ \left(\vec{W}^{*}_{I}\right)^{l+1} &= \vec{W}^{(0)}_{I} - \frac{\alpha_{m}\Delta t^{*}_{I}}{\Omega^{n+1}_{I}}\vec{R}^{*}_{I}\left(\vec{W}^{(m-1)}\right), \end{aligned} \tag{6.71}

其中 $$l$$ 表示当前伪时间层, $$(l+1)$$ 表示新伪时间层。时间推进过程或者以 $\left(\vec{W}^{*}\right)^{l} = \vec{W}^{n}$ 启动ï¼或者以由先前物理时间步外推的值启动ï¼例如[98]en

where $$l$$ denotes the actual and $$(l+1)$$ the new pseudo-time level, respectively. The time-marching process is started either with $\left(\vec{W}^{*}\right)^{l} = \vec{W}^{n}$ or with a value extrapolated from previous physical time steps, e.g., [98]

\left(\vec{W}^{*}_{I}\right)^{l} = \vec{W}^{n} + \frac{3\vec{W}^{n} - 4\vec{W}^{n-1} + \vec{W}^{n-2}}{2}. \tag{6.72}

推进一直持续到 $\left(\vec{W}^{*}_{I}\right)^{l+1}$ 以足够的精度逼近 $\vec{W}^{n+1}_{I}$ (通常当残差 $\vec{R}^{*}_{I}$ 降低了两个或三个数量级时)。此后进行下一个物理时间步。伪时间步长 $\Delta t^{*}$ 按6.1.4小节中同样方式计算。en

It is continued until $\left(\vec{W}^{*}_{I}\right)^{l+1}$ approximates $\vec{W}^{n+1}_{I}$ with sufficient accuracy (usually when the residual $\vec{R}^{*}_{I}$ was reduced by two or three orders of magnitude). After that, the next physical time step is conducted. The pseudo time step $\Delta t^{*}$ is computed in the same way that we saw in Subsection 6.1.4.

Arnone等人[98]指出ï¼当物理时间步长 $\Delta t$ 与伪时间步长 $\Delta t^{*}$ 同量级或更小时ï¼多级格式(6.71)会变得不稳定。Melson等人[99]证明ï¼不稳定性由如下项引起en

Arnone et al. [98] pointed out that the multistage scheme (6.71) becomes unstable when the physical time step $\Delta t$ is of the order of the pseudo time step $\Delta t^{*}$ or smaller. Melson et al. [99] demonstrated that the instability is caused by the term

\frac{3}{2\Delta t}\left(\Omega\bar{M}\right)^{n+1}_{I}\vec{W}^{*}_{I} \tag{3}

它出现在式(6.69)中ï¼在 $\Delta t$ 较小时变得显著。他们建议对这一项作隐式处理。于是ï¼必须修改式(6.71)的多级格式ï¼使第 $$k$$ 级变为[99]en

in Eq. (6.69), which becomes significant for small $\Delta t$ . They suggested an implicit treatment of this term. Thus, we have to modify the multistage scheme in Eq. (6.71) such that the $$k$$ -th stage becomes [99]

\begin{aligned} \vec{W}^{(k)}_{I} &= \vec{W}^{(0)}_{I} - \frac{\alpha_{k}\Delta t^{*}_{I}}{\Omega^{n+1}_{I}}\left[\bar{I} + \frac{3}{2\Delta t}\alpha_{k}\Delta t^{*}_{I}\bar{M}^{n+1}\right]^{-1}\\ &\quad\cdot\left[\vec{R}_{I}\left(\vec{W}^{(k-1)}\right) + \frac{3}{2\Delta t}\left(\Omega\bar{M}\right)^{n+1}_{I}\vec{W}^{(0)}_{I} - \vec{Q}^{*}_{I}\right]. \end{aligned} \tag{6.73}

同样的方法也适用于混合多级格式(见6.1.2小节)。上述形式(6.73)对任何物理时间步长 $\Delta t$ 都是稳定的[99]。en

The same methodology can also be applied to a hybrid multistage scheme (see Subsection 6.1.2). The above formulation (6.73) is stable for any physical time step $\Delta t$ [99].

对单元中心格式ï¼式(6.73)中的质量矩阵 $\bar{M}^{n+1}$ 可以集中化(用单位矩阵代替)而不降低解的精度。这样ï¼式中的项en

In the case of cell-centred schemes, the mass matrix $\bar{M}^{n+1}$ in Eq. (6.73) can be lumped (substituted by the identity matrix) without reducing the solution accuracy. In this way, the term

\left[\bar{I} + \frac{3}{2\Delta t}\alpha_{k}\Delta t^{*}_{I}\bar{M}^{n+1}\right]^{-1} \tag{5}

便变成一个标量值。然而ï¼对单元顶点空间离散格式ï¼则必须考虑质量矩阵ï¼否则多级格式在小物理时间步长下会不稳定。为了避免对 $\bar{M}^{n+1}$ 的昂贵求逆,Venkatakrishnan[100]以及Venkatakrishnan与Mavriplis[101]对式(6.73)提出了如下修改en

in Eq. (6.73) is turned into a scalar value. However, we have to account for the mass matrix in the case of a cell-vertex spatial discretisation scheme. Otherwise, the multistage scheme will be unstable for small physical time steps. In order to circumvent the expensive inversion of $\bar{M}^{n+1}$ , Venkatakrishnan [100] and Venkatakrishnan and Mavriplis [101] suggested the following modification to Eq. (6.73)

\begin{aligned} \vec{W}^{(k)}_{I} &= \vec{W}^{(0)}_{I} - \frac{\alpha_{k}\Delta t^{*}_{I}}{\Omega^{n+1}_{I}}\left[1 + \frac{3}{2\Delta t}\alpha_{k}\Delta t^{*}_{I}\beta\right]^{-1}\\ &\quad\cdot\left[\vec{R}^{*}_{I}\left(\vec{W}^{(k-1)}\right) - \frac{3}{2\Delta t}\Omega^{n+1}_{I}\beta\,\vec{W}^{(k-1)}_{I}\right]. \end{aligned} \tag{6.74}

参数 $\beta$ 此时可用来稳定时间推进格式。实践中发现取 $\beta = 2$ 即已足够[100]、[101]。en

The parameter $\beta$ can now be utilised to stabilise the time-stepping scheme. In practice, choosing $\beta = 2$ was found sufficient [100], [101].

用显式多级格式获得伪时间内解的双时间步进方法应用广泛。当多级格式由局部时间步进(在 $t^{*}$ 内)与多重网格加速时ï¼计算效率最高。原因是多重网格格式能快速(在少数几个循环内)收敛到式(6.68)的定常解。结构网格上的应用实例可见文献[97]-[99]以及[102]-[104];该方法在非结构网格上的实现例如见于[100]、[101]与[105]。en

The dual time-stepping approach, where the solution in pseudo time is obtained by an explicit multistage scheme, is widely used. The highest computational efficiency results when the multistage scheme is accelerated by local time-stepping (in $t^{*}$ ) and multigrid. The reason is that the multigrid scheme converges quickly (in few cycles) to the stationary solution of Eq. (6.68). Examples of applications on structured grids can be found in Refs. [97]-[99] as well as in [102]-[104]. Implementations of the methodology on unstructured grids were described, e.g., in [100], [101] and [105].

6.3.2 Dual Time-Stepping for Implicit Schemes 隐式格式的双时间步进[cfd-6-3-2]

用隐式格式在伪时间 $t^{*}$ 内求解式(6.68)ï¼其实现方式与6.2节所述相同。首先ï¼把式(6.68)写成非线性隐式格式ï¼即en

The implementation of an implicit scheme for the solution of Eq. (6.68) in pseudo time $t^{*}$ proceeds in the same way as outlined in Section 6.2. First of all, we formulate Eq. (6.68) as an nonlinear implicit scheme, i.e.,

\frac{\partial}{\partial t^{*}}\left(\Omega^{n+1}_{I}\vec{W}^{*}_{I}\right) = -\left(\vec{R}^{*}_{I}\right)^{l+1} \tag{6.75}

其中 $$(l+1)$$ 为新伪时间层。再次注意ï¼时间导数中没有 $\bar{M}$ 。式(6.69)定义的非定常残差可在伪时间内线性化如下en

with $$(l+1)$$ being the new pseudo-time level. Note again the absence of $\bar{M}$ in the time derivative. The unsteady residual, which is defined in Eq. (6.69), can be linearised in pseudo time as follows

\left(\vec{R}^{*}\right)^{l+1} \approx \left(\vec{R}^{*}\right)^{l} + \frac{\partial\vec{R}^{*}}{\partial\vec{W}^{*}}\Delta\vec{W}^{*}, \tag{6.76}

其中 $\Delta\vec{W}^{*} = \left(\vec{W}^{*}\right)^{l+1} - \left(\vec{W}^{*}\right)^{l}$ ï¼通量雅可比定义为en

where $\Delta\vec{W}^{*} = \left(\vec{W}^{*}\right)^{l+1} - \left(\vec{W}^{*}\right)^{l}$ and the flux Jacobian is defined as

\frac{\partial\vec{R}^{*}}{\partial\vec{W}^{*}} = \frac{\partial\vec{R}}{\partial\vec{W}} + \frac{3}{2\Delta t}\left(\Omega\bar{M}\right)^{n+1}. \tag{6.77}

把上述线性化代入式(6.75)ï¼便得到未分解的隐式格式[106]en

If we insert the above linearisation into Eq. (6.75), we obtain the unfactored implicit scheme [106]

\left[\left(\frac{1}{\Delta t^{*}_{I}} + \frac{3}{2\Delta t}\right)\left(\Omega\bar{M}\right)^{n+1}_{I} + \left(\frac{\partial\vec{R}}{\partial\vec{W}}\right)_{I}\right]\Delta\vec{W}^{*} = -\left(\vec{R}^{*}_{I}\right)^{l}. \tag{6.78}

6.2节介绍的任何方法都可用于求解方程组(6.78)。关于时间精确隐式方法的详细讨论可参见[107];近年实现的例子可参见例如[106]与[108]。en

Any of the methodologies presented in Section 6.2 can be employed for the solution of the system (6.78). A detailed discussion of time-accurate implicit methods can be found in [107]. For recent examples of implementations, the reader is referred to, e.g., [106] and [108].