1 引 言

非定常流动问题广泛存在于航空航天领域, 比如直升机旋翼 (Choi et al.2008)、叶轮机 (Guo et al. 2012)、气动弹性问题 (Farhat & Lesoinne 2000)等. 根据脉动响应的周期性特征,非定常流动可以分为周期性流动、准周期流动和非周期流动.比如直升机定直前飞时旋翼的绕流就是典型的周期性流动,机翼的气弹响应绕流为准周期流动, 湍流槽道流动属于非周期流动.对于周期性或准周期流动,根据非定常减缩频率大小分为高减缩频率流动、低减缩频率流动以及未知周期的流动.比如扑翼流动属于高减缩频率流动,刚性飞行器动导数计算中强迫振荡属于低减缩频率流动,静止非定常圆柱绕流属于未知周期的流动.

非定常流动的研究手段主要有理论模型、数值模拟和风洞试验.理论模型虽然计算量小, 但多基于线化势流理论,精度不足且无法用于跨声速或大迎角等非线性流动. 风洞试验虽然可靠性较高,但除了难度大、周期长、成本高外, 还存在洞壁及支撑干扰等复杂非定常修正问题.从国内外发展趋势来看, 通过计算流体力学方法 (computational fluid dynamic,CFD) 已经逐渐成为非定常流动研究的主要手段.

相对于定常计算, 非定常数值模拟存在计算量和存储量大的问题,很难在工程中广泛应用. 即使是在计算机配置日新月异的今天,计算机的计算能力也远远跟不上非定常流动模拟对其的需求.解决这种矛盾亟需发展高效高精度的流场求解算法. 因此,在数值计算中如何能够实现非定常流场的高效高精度求解便成了当前的研究热点.但以往多数计算流体力学研究者将关注点放在空间离散与迭代算法的改进上,忽略了时间离散方法. 近期的很多研究发现, 采用合适的时间离散方法也可使CFD计算的精度和效率有非常明显的提高. 因此,近些年尤其是进入21世纪以来时间离散方法得到了迅速发展,各种时间离散方法如雨后春笋般涌现.

当前针对非定常流场的时间离散格式主要是时域向后差分法 (backward difference format, BDF) (Butcher 2016). 其主要思想是先求解定常流场,然后通过时域推进顺次求解不同时刻的流场.最常见的就是二阶时域向后差分离散(BDF2).但这类方法在时间上只能够达到二阶精度,并且对于周期性流动需要若干周期才能达到稳定, 没有利用流动的周期性特征,计算效率不理想.

由于在实际工程应用中很多涉及到的非定常流动均具有周期性特征,于是针对周期性非定常流动, 近年来出现了线化频域法 (linear frequency domain solver, LFD) (Hall & Crawley 1989)、非线性频域方法(non-linear frequency domain method, NLFD) (Ning & He 1997,Mcmullen et al. 2001)、谐波平衡法 (harmonic balance method, HB)(Hall et al. 2002)、时间谱方法 (time spectral method, TS)(Gopinath & Jameson 2005, Van der Weide et al.2005)、时间有限差分法 (finite difference method in time, FDMT)(Leffell et al. 2016) 等时间离散格式. 对于准周期问题,发展了混合向后差分时间谱方法 (BDF/TS) (Yang & Mavriplis 2010) 以及切比雪夫伪谱法 (Dinu et al. 2006). 对于纯非周期问题,发展了时间谱元法 (spectral element method in time, SEMT) (Kurdi & Beran 2008).这些新发展的时间离散格式虽然精度和效率相比于传统的BDF有较为明显的提升,但适用范围有限. 如果时间离散格式应用不当,其计算精度和效率会反不如BDF, 甚至失效. 比如Jameson & Shankaran (2009) 将时间谱方法应用到减缩频率大的扑翼中,经过研究表明时间谱方法计算扑翼无论是精度还是收敛性均不如传统的BDF.因此, 研究并确定每一种时间离散格式的适用范围至关重要.

对于这些新发展的时间离散格式还需要较全面的横向比较. Ronch等 (2010, 2013)

对比了时域推进法、线化频域法和谐波平衡法在数值预测动导数精度、计算效率与消耗内存方面的能力.

Leffell等 (2016)

对于周期性非定常流动中的时域向后差分法、时间谱方法以及时间有限差分法做了比较详细的精度与效率对比,

并对时空并行算法进行了阐述. 但这对于所有的时间离散方法来说只是冰山一角,

有必要对近些年来蓬勃发展的多样化时间离散格式做一个系统的总结与对比.

本文总结目前非定常流场计算中主流的时间离散方法,根据每种方法时间离散形式的不同分为4大类：时域推进法、频域谐波法、时域配点法和混合方法.分别介绍各类方法的具体原理并比较各方法的计算精度、效率以及适用范围. 另外,对于新发展的时间离散格式在工程问题中的应用及效果也做了相应的归纳. 最后,对未来时间离散方法的发展方向进行了展望.

2 时域推进法

时域推进法的步骤是先算初始定常流动,之后沿物理时间向前推进顺次求出往后各时刻的流场. 主要代表就是时域向后差分法,这也是最初始且目前应用最广的时间离散方法.

时域向后差分法, 顾名思义,即对于时间导数进行向后差分离散之后进行实时间步的推进.二阶BDF关于时间导数的离散形式如下

$$\dfrac{{\rm d}\pmb U _i^n }{{\rm d}t} = \dfrac{3\pmb U _i^n - 4\pmb U _i^{n - 1} + \pmb U _i^{n - 2} }{2\Delta t}(1)$$

式中, $\pmb U $表示守恒量, $\Delta t$表示时间间隔.离散后的非定常控制方程为

$$V_i \dfrac{3\pmb U _i^n - 4\pmb U _i^{n - 1} + \pmb U _i^{n - 2} }{2\Delta t} + \pmb R \left( {\pmb U _i^n } \right)=\pmb 0(2)$$ 式中, $V_i $表示网格体积, $\pmb R $表示残差.Briley 和Mcdonald (1975) 与Beam 和Warming (1976)等通过线化非线性残差来隐式求解式(2).但用直接法求解如此庞大的系数矩阵的逆矩阵, 计算量难以接受,因此需要迭代求解, 即所谓的``单时间步法''. 之后, Jamson (1991)提出了"虚拟时间亚迭代法", 该方法又称"双时间步法".即在冻结的物理时间点上巧妙地引入虚拟时间迭代过程,可将定常流动计算中的加速收敛措施应用到非定常流动的计算中,思想简单、实现容易. 因此, 双时间步BDF很快得到了广泛的应用 (Briley & Mcdonald 1975). 其非定常控制方程如下所示

$$ V_i \dfrac{{\rm d}\pmb U _i^n }{{\rm d}\tau} + V_i \dfrac{3\pmb U_i^n - 4\pmb U _i^{n - 1} + \pmb U _i^{n - 2} }{2\Delta t} + \pmb R \left( {\pmb U _i^n } \right) =\pmb 0(3)$$

式中$\tau$表示伪时间. BDF的优点在于其可以适用于任意的非定常流动,方法容易理解且在定常流动计算程序上的改动工作量很小.但从计算精度上而言, 三阶或以上的BDF会导致产生非物理解甚至发散(Hull et al. 1972).因此BDF的精度很难达到高阶且需要较多的时刻点才能保证足够的时间方向准确度.

从计算效率上而言, 对于工程中常见的周期性流动,BDF需要若干周期才能达到稳定, 计算时间之长是在工程应用中难以忍受的.也正因为BDF有这些缺陷, CFD工作者需要寻找更优的时间离散方法,促使了近期各种时间离散方法的诞生.

3 频域谐波法

由于在工程中多数的非定常流动均具有周期性特征,通过利用流动的周期性特征将原非定常流场的控制方程做傅里叶变换,使之转化为频域的一组定常方程, 进而对该频域方程组进行求解获得各频域谐波分量,最后通过傅里叶反变换获得时域响应. 其计算主要在频域内进行,因而称之为频域谐波法.

3.1 线性频域谐波方法

早期频域方法在非定常气动力模型的构建 (Zeiler 2000) 中使用的,比如著名的西奥道生气动力模型 (Brunton & Rowley 2013).近期Hall等 (1989)在用欧拉方程数值模拟涡轮机的实践中发展了频域的非定常流场数值求解方法------线性谐波法,也由此拉开了时间离散方法变革的序幕.

Hall等 (1989) 假设流场原始变量由时均量以及弱扰动量叠加而成.如式(4)所示

$$\left. \begin{array}{l} \rho \left( {x,y,t} \right) = \bar {\rho }\left( {x,y} \right) + \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {\rho }} \left( {x,y,t} \right) \\ u\left( {x,y,t} \right) = \bar {u}\left( {x,y} \right) + \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {u}} \left( {x,y,t} \right) \\ v\left( {x,y,t} \right) = \bar {v}\left( {x,y} \right) + \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {v}} \left( {x,y,t} \right) \\ p\left( {x,y,t} \right) = \bar {p}\left( {x,y} \right) + \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {p}} \left( {x,y,t} \right) \\ \end{array}\right\}(4)$$

将式(4)代入欧拉方程中, 把扰动部分当成小量来处理并忽略高阶项,将原欧拉方程分为与时间无关的定常欧拉方程以及与时间有关的线化扰动方程.求解定常的时均流动方程如下

$$\left. \begin{array}{l} \dfrac{\partial \pmb F }{\partial x} + \dfrac{\partial \pmb G }{\partial y} =\pmb 0 \\ \pmb F = \left[ {{\begin{array}{*{20}c} {\bar {\rho }\bar {u}} \\ {\bar {\rho }\bar {u}^2 + \bar {p}} \\ {\bar {\rho }\bar {u}\bar {v}} \\ {\dfrac{\gamma }{\gamma - 1}\bar {p}\bar {u} + \dfrac{1}{2}\bar {\rho }\bar {u}\left( {\bar {u}^2 + \bar {v}^2} \right)} \\ \end{array} }} \right] \\ \pmb G = \left[ {{\begin{array}{*{20}c} {\bar {\rho }\bar {v}} \\ {\bar {\rho }\bar {u}\bar {v}} \\ {\bar {\rho }\bar {v}^2 + \bar {p}} \\ {\dfrac{\gamma }{\gamma - 1}\bar {p}\bar {v} + \dfrac{1}{2}\bar {\rho }\bar {v}\left( {\bar {u}^2 + \bar {v}^2} \right)} \\ \end{array} }} \right] \\ \end{array}\right\}(5)$$

对应弱扰动量的非定常欧拉方程如下所示

$$\left. \begin{array}{l} \dfrac{\partial }{\partial t}\pmb B _1 \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {\pmb W}} + \dfrac{\partial }{\partial x}\pmb B _2 \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {\pmb W}} + \dfrac{\partial }{\partial y}\pmb B _3 \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {\pmb W}} =\pmb 0 \\ \pmb B _1 = \left[ {{\begin{array}{*{20}c} 1 \quad & 0 \quad & 0 \quad & 0 \\ \bar {u} \quad & \bar {\rho } \quad & 0 \quad & 0 \\ \bar {v} \quad & 0 \quad & \bar {\rho } \quad & 0 \\ {\dfrac{1}{2}\left( {\bar {u}^2 + \bar {v}^2} \right)} \quad & {\bar {\rho }\bar {u}} \quad & {\bar {\rho }\bar {v}} \quad & {\dfrac{1}{\gamma - 1}} \\ \end{array} }} \right] \\ \pmb B _2 = \left[ {{\begin{array}{*{20}c} \bar {u} \quad & \bar {\rho } \quad & 0 \quad & 0 \\ {\bar {u}^2} \quad & {2\bar {\rho }\bar {u}} \quad & 0 \quad & 1 \\ {\bar {u}\bar {v}} \quad & {\bar {\rho }\bar {v}} \quad & {\bar {\rho }\bar {u}} \quad & 0 \\ {\dfrac{1}{2}\bar {u}\left( {\bar {u}^2 + \bar {v}^2} \right)} \quad & {\dfrac{\gamma }{\gamma - 1}\bar {p} + \dfrac{3}{2}\bar {\rho }\bar {u}^2 + \dfrac{1}{2}\bar {\rho }\bar {v}^2} \quad & {\bar {\rho }\bar {u}\bar {v}} \quad & {\dfrac{\gamma \bar {u}}{\gamma - 1}} \\ \end{array} }} \right] \\ \pmb B _3 = \left[ {{\begin{array}{*{20}c} \bar {v} \quad & 0 \quad & \bar {\rho } \quad & 0 \\ {\bar {u}\bar {v}} \quad & {\bar {\rho }\bar {v}} \quad & {\bar {\rho }\bar {u}} \quad & 0 \\ {\bar {v}^2} \quad & 0 \quad & {2\bar {\rho }\bar {v}} \quad & 1 \\ {\dfrac{1}{2}\bar {v}\left( {\bar {u}^2 + \bar {v}^2} \right)} \quad & {\bar {\rho }\bar {u}V} \quad & {\dfrac{\gamma }{\gamma - 1}\bar {p} + \dfrac{1}{2}\bar {\rho }\bar {u}^2 + \dfrac{3}{2}\bar {\rho }\bar {v}^2} \quad & {\dfrac{\gamma \bar {v}}{\gamma - 1}} \\ \end{array} }} \right] \\ \end{array}\right\}(6)$$

其中$\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {\pmb W}} $为原始变量的扰动分量, 即$\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {\pmb W}} = \left[ {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {\rho }} ,\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {u}} ,\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {v}} ,\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {p}} } \right]^{\rm T}$. 从式(6)不难发现,由于时均量在时均流动方程中均已求出, 此扰动方程为线化方程.假设流动的弱扰动量为谐振形式, 即可表达为如下形式

$$\left.\begin{array}{l} \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {\rho}} \left( {x,y,t} \right) = \tilde {\rho }\left( {x,y} \right){\rm e}^{{\rm j}wt} \\ \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {u}} \left( {x,y,t} \right) = \tilde {u}\left( {x,y} \right){\rm e}^{{\rm j}wt} \\ \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {v}} \left( {x,y,t} \right) = \tilde {v}\left( {x,y} \right){\rm e}^{{\rm j}wt} \\ \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {p}} \left( {x,y,t} \right) = \tilde {p}\left( {x,y} \right){\rm e}^{{\rm j}wt} \\ \end{array}\right\}(7)$$

式中$\omega $为角频率. 代入式(6)中整理得到

$$jwB_1 \tilde {\pmb W} + \dfrac{\partial}{\partial x}B_2 \tilde {\pmb W} + \dfrac{\partial }{\partial y}B_3 \tilde {\pmb W} =\pmb 0(8)$$

其中$\tilde {\pmb W} = \left[ {\tilde {\rho }\left( {x,y} \right),\tilde {u}\left( {x,y} \right),\tilde {v}\left( {x,y} \right),\tilde {p}\left( {x,y} \right)} \right]^{\rm T}$. 通过求解式(8)得到$\tilde {\pmb W}$进而得到弱扰动量$\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {\pmb W}}$,最后进行时均量与弱扰动量的叠加得到非定常流场的时域响应.

为了验证此求解方法的精度和效率, Hall和 Crawley (1989)采用亚声速翼型俯仰算例,其压力对比如 图1所示. 从 图1可以看出, 在满足假设条件前提下,当前频域计算方法与传统时域求解方法吻合很好, 最大误差不超过0.01,而效率提高了几十倍.

显示原图| 下载原图ZIP| 生成PPT

图 1   翼型表面压力系数$C_p$ 对比 (Hall & Crawley 1989)

   

Hall和Clark (1991) 针对守恒量进行谐波展开, Hall和 Crawley (1994)针对跨声速流动运用激波捕获法并套用Lax-Wendroff格式推导出相应的公式如下

$$\left( {1 - {\rm e}^{{\rm j}w\Delta t}} \right)u_i - \dfrac{\Delta t}{2\Delta x}\left[ {\left. {\dfrac{\partial F}{\partial U}} \right|_{i + 1} u_{i + 1} - \left. {\dfrac{\partial F}{\partial U}} \right|_{i - 1} u_{i - 1} } \right]+ $$ \n \quad $$ \dfrac{\Delta t^2}{2\Delta x^2}\left[ {\left. {\dfrac{\partial F}{\partial U}} \right|_{i + 1 / 2} \left( {\left. {\dfrac{\partial F}{\partial U}} \right|_{i + 1} u_{i + 1} - \left. {\dfrac{\partial F}{\partial U}} \right|_i u_i } \right) - \left. {\dfrac{\partial F}{\partial U}} \right|_{i - 1 / 2} \left( {\left. {\dfrac{\partial F}{\partial U}} \right|_i u_i - \left. {\dfrac{\partial F}{\partial U}} \right|_{i - 1} u_{i - 1} } \right)} \right]+$$ \n \quad $$ \dfrac{\Delta t^2}{2\Delta x^2}\left[ {\left. {\dfrac{\partial ^2F}{\partial U^2}} \right|_{i + 1 / 2} u_{i + 1 / 2} \left( {F_{i + 1} - F_i } \right) - \left. {\dfrac{\partial ^2F}{\partial U^2}} \right|_{i - 1 / 2} u_{i - 1 / 2} \left( {F_i - F_{i - 1} } \right)} \right] = 0(9)$$

Pechloff和Laschka (2006)也发展了一种基于小扰动理论的线性谐波法来计算非定常RANS方程.在此之后Widhalm等 (2010)对于守恒量进行时均值与扰动分量的叠加并通过自动微分法求解扰动方程中的雅可比矩阵,并取名为线化频域法 (linearized frequency domain solver, LFD).

线性频域谐波方法虽然很大程度上减少了计算量,但是仅仅局限于小扰动弱非线性的周期性流动, 局限性较强从而限制了其应用,而对于一些复杂的周期性流动, 就必须考虑流动的非线性. 因此,就有了非线性频域谐波方法的发展.

3.2 非线性频域谐波方法

为了能够将频域方法应用于非线性较强的周期性流动问题, Ning 和He(1997) 基于准三维非定常欧拉方程提出了非线性谐波法,在线性谐波法时均方程中引入额外的非定常应力项,其时均方程如式(10)所示

$$\left.\begin{array}{l} \dfrac{\partial \bar{\pmb F}}{\partial x} + \dfrac{\partial \bar{\pmb G}}{\partial y} =\pmb 0 \\ \bar{\pmb F} = \left[ {{\begin{array}{*{20}c} {\overline {\rho U} - \overline {\rho U_g } } \\ {\bar {U}\left( {\overline {\rho U} - \overline {\rho U_g } } \right) + \bar {P} + \overline {\left( {\widehat{\rho U}} \right)\widehat{U}} - \overline {\left( {\widehat{\rho U_g }} \right)\widehat{U}} } \\ {\bar {V}\left( {\overline {\rho U} - \overline {\rho U_g } } \right) + \overline {\left( {\widehat{\rho U}} \right)\widehat{V}} - \overline {\left( {\widehat{\rho U_g }} \right)\widehat{V}} } \\ {\overline H \left( {\overline {\rho U} - \overline {\rho U_g } } \right) + \bar {P}\overline U _g + \overline {\widehat{H}\left( {\widehat{\rho U}} \right)} + \overline {\widehat{P}\widehat{U_g }} - \overline {\widehat{H}\left( {\widehat{\rho U_g }} \right)} } \\ \end{array} }} \right] \\ \bar{\pmb G} = \left[ {{\begin{array}{*{20}c} {\overline {\rho V} - \overline {\rho V_g } } \\ {\bar {U}\left( {\overline {\rho V} - \overline {\rho V_g } } \right) + \bar {P} + \overline {\left( {\widehat{\rho V}} \right)\widehat{U}} - \overline {\left( {\widehat{\rho V_g }} \right)\widehat{U}} } \\ {\bar {V}\left( {\overline {\rho V} - \overline {\rho V_g } } \right) + \overline {\left( {\widehat{\rho V}} \right)\widehat{V}} - \overline {\left( {\widehat{\rho V_g }} \right)\widehat{V}} } \\ {\overline H \left( {\overline {\rho V} - \overline {\rho V_g } } \right) + \bar {P}\overline V _g + \overline {\widehat{H}\left( {\widehat{\rho V}} \right)} + \overline {\widehat{P}\widehat{V_g }} - \overline {\widehat{H}\left( {\widehat{\rho V_g }} \right)} } \\ \end{array} }} \right] \\ \end{array}\right\}(10)$$

式中, $U_g $与$V_g $表示网格运动速度.采用四阶龙格库塔将时均方程与扰动方程进行耦合求解, 求解流程如 图2所示.

显示原图| 下载原图ZIP| 生成PPT

图 2   非线性谐波法求解流程(Ning & He 1997)

   

通过在时均方程中引入非定常应力项, 在单倍谐波情况下非线性谐波法相比于线性谐波法可以把扰动幅值对流动特征的影响更好地刻画出来, 如 图3所示. 同时非线性谐波法可以通过控制扰动源数目和傅里叶阶数来控制求解精度, 然而求解阶数越高, 计算耗费的时间也就越多. 相比于线性谐波法效率会有所折扣 (Vilmin et al. 2006), 但对于非线性较强的周期性问题非线性谐波法则有更高的精度.

显示原图| 下载原图ZIP| 生成PPT

图 3   跨声速槽道流中非定常压力扰动项(Ning & He 1997)

   

He 和 Ning (1998) 引入伪时间项,并将非线性谐波法扩展到二维黏性流动, 应用到叶轮机数值模拟中.之后Chen等 (2000) 将其应用到三维湍流问题,通过对德国宇航院的单级跨声速压气机研究认为二阶谐波是在计算成本和计算精度间最佳折衷.

然而, 随着谐波数增加, 非线性谐波法中通量的频域求解就会变得复杂且计算效率会严重降低. 之后基于非线性谐波法基础上, 2001年McMullen提出了非线性频域法(Non-linear frequency domain method, NLFD) (Mcmullen et al. 2001, 2002, 2006, Mcmullen 2003). 在一个周期内取$N$个采样点, 之后对守恒量和通量均进行傅里叶变换, 并引入伪时间项得到一组定常频域方程. 其频域方程如下

$$ V\dfrac{{\rm d} \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {\pmb U}} }{{\rm d}\tau } + j\omega V \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {\pmb U}} + \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {\pmb R}}_k =\pmb 0(11)$$

式中, $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {\pmb U}} $表示守恒量在频域内的分量, $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {\pmb R}} $为残差在频域内分量. 其计算流程如{\bf 图4}所示,对于初始守恒量的谐波分量$\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {\pmb U}}_k^0 $通过傅里叶逆变换得到守恒量时域值$\pmb U ( t )$,计算对应的时域残差$\pmb R ( t )$并通过傅里叶变换得到其频域分量$\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {\pmb R}} _k^0 $, 进而通过式(11)求出守恒量的谐波分量$\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\frown$}}\over {\pmb U}} _k^1 $,再进行下一步迭代, 直至收敛.

显示原图| 下载原图ZIP| 生成PPT

图 4   非线性频域法求解流程(Mcmullen & Jameson et al. 2001)

   

此算法相比于非线性谐波法的优点在于其在时域内计算残差,可以有效地计算扰动大的周期性问题以及复杂流动问题且计算精度和效率更高.McMullen 和 Jameson (2006)对于非线性频域法的计算效率及其关于减缩频率和迎角等敏感度的研究,结果表明NLFD方法计算效率相比于时域推进法普遍提高了3$\sim $9倍.

4 时域配点法

进入21世纪后, 时域配点方法开始兴起. 最早是Hall等 (2002)提出的时域谐波平衡法,在此之后又出现了时间谱方法、时间谱元法等一系列时间离散格式.这类方法的特点在于通过将非定常流场的时间响应在给定的一组正交多项式系上投影,然后在非定常流场中选$N$个时刻,将非定常流动问题转化为$N$个时刻点相互耦合的定常问题. 这类方法叫做时域配点法.

4.1 时间谱方法

2002年, Hall等 (2002) 提出了应用于流体数值求解的谐波平衡法,将守恒量在频域内进行时间离散后代回时域方程进行求解.将周期性的非定常问题转化为$N$个时刻定常方程的耦合.利用傅里叶变换和逆变换求出谱矩阵, 最终的迭代计算仍然在时域内进行.但谐波平衡法这一概念在较早的计算力学中已经提出.其原理为假设近似解为傅里叶级数的形式代入方程通过各自谐波自相平衡来求得傅里叶系数进而求得原方程的解,是典型的频域方法. 因此, 为了区分, 后续研究将流体中谐波平衡法的概念进行推广,谐波平衡法作为一大类包括频域的谐波平衡法和时域的谐波平衡法.2002年Hall提出的谐波平衡法又被称为时域谐波平衡法或高维谐波平衡法.

2005年Gopinath等 (Gopinath & Jameson 2005, Van der Weide et al. 2005)进一步提出了时间谱方法,并针对Hall提出的谐波平衡法推出了更为简洁的时间谱矩阵的表达式.谱方法的概念是建立一个满足周期性边界条件的完备函数族所构成的谱空间,将方程的解在给定的谱空间上投影后再带回原方程进行求解,离散是在谱空间完成的. 对于足够光滑的物理问题,谱方法能达到很高的精度和效率.但对于某些完备函数族不满足周期性边界条件,则这类完备函数族所构成的谱空间称为伪谱空间,对应的求解方法称为伪谱法. 下面对时间谱方法的公式进行推导.

如果流动具有周期性, 则在每一个计算单元内的守恒量$\pmb U _i$随时间均呈周期性变化. 在一个周期$T$内设置$N$个等时间间距的采样点,若$N$为奇数, 对守恒量$\pmb U _i $进行正向傅里叶变换得

$$\widehat{\pmb U }_i^k \approx \dfrac{1}{N}\sum\limits_{n = 0}^{N - 1} {\pmb U _i^n {\rm e}^{ - {\rm j}k\omega n\Delta t}}(12)$$

对应的傅里叶反变换为

$$ \pmb U _i^n \approx \sum\limits_{k = - (N - 1) / 2}^{(N - 1) / 2} {\widehat{\pmb U }_i^k } {\rm e}^{{\rm j}k\omega n\Delta t}(13)$$

将式(13)对时间求导可得

$$\dfrac{{\rm d}}{{\rm d}t}\pmb U _i^n \approx \sum\limits_{k = - (N - 1) / 2}^{(N - 1) / 2} {jk\omega \widehat{\pmb U }_i^k } {\rm e}^{{\rm j}k\omega n\Delta t}(14)$$

将式(12)代入式(14)中经过化简可得

$$\dfrac{{\rm d}}{{\rm d}t}\pmb U _i^n \approx \sum\limits_{j_n = 0}^{N - 1} {d_n^{j_n } \pmb U _i^{j_n } }(15)$$

其中

$$d_n^{j_n } = \left\{ \begin{array}{ll} \dfrac{2\pi }{T}\cdot \dfrac{1}{2}( - 1)^{n - j_n }{\rm csc }\left(\dfrac{\pi \left( {n - j_n } \right)}{N}\right), \quad & {n \neq j_n } \\ 0, \quad & n = j_n \\ \end{array} \right. $$

若$N$为偶数, 则类似进行计算可得

$$d_n^{j_n } = \left\{ \begin{array}{ll} \dfrac{2\pi }{T}\cdot \dfrac{1}{2}( - 1)^{n - j_n }\cot \left(\dfrac{\pi \left( {n - j_n } \right)}{N}\right) , \quad & n \neq j_n \\ 0, \quad & n = j_n \\ \end{array} \right. $$

于是, 可将含有$N$个采样点在第$n$个时刻的Navier-Stokes (N-S)方程的半离散形式写为

$$\sum\limits_{j_n = 0}^{N - 1} {d_n^{j_n } V_i^{j_n } } \pmb U _i^j + \pmb R _i^n =\pmb 0,\qquad n = 0,1,{2},\cdots ,N - 1(16)$$

在式(16)中引入伪时间项, 则该方程可进一步写为

$$\dfrac{\partial }{\partial \tau _i^n }\left( {V_i^n \pmb U _i^n } \right) + \sum\limits_{j_n = 0}^{N - 1} {d_n^{j_n } V_i^j } \pmb U _i^j + \pmb R _i^n =\pmb 0,\qquad n = 0,1,{2},\cdots ,N - 1(17)$$

如果对伪时间项进行显示离散,为保证计算稳定性, 伪时间步长须限制在

$$\tau _i^n = CFL\dfrac{V_i^n }{\left\| \lambda \right\| + V_i^n k'}(18)$$

式中, $\lambda $为谱半径,$k'$为最大谐波数, 具体表达式如下

$$k' = \left\{ \begin{array}{ll} \dfrac{\pi N}{T}, \quad & {\rm mod} (N,2) = 1 \\ \dfrac{\pi \left( {N - 1} \right)}{T}, \quad & {\rm mod} (N,2) = 0 \\ \end{array} \right.$$

随着无量纲时间步长的减小, 迭代的伪时间步长会非常小,从而严重降低计算效率. 采用隐式格式可以放宽时间步长限制,因此在计算中多采用隐式格式对伪时间项进行离散.对伪时间导数进行一阶向后差分离散得到式(19)所示的迭代格式

$$\pmb A ^n\Delta \pmb U ^n = - \sum\limits_{j_n = 0}^{N - 1} {d_n^{j_n } V^{j_n }} \pmb U ^{j_n } - \pmb R ^n(19)$$

其中

$$\pmb A = \left[ \begin{array}{cccc} \dfrac{V_0 }{\Delta \tau _0 }{\pmb I} + {\pmb J}_0 \quad & {V_1 d_0^1 \pmb I } \quad & \cdots \quad & {V_{N - 1} d_0^{N - 1} \pmb I } \\ {V_0 d_1^0 \pmb I } \quad & {\dfrac{V_1 }{\Delta \tau _1 }\pmb I + {\pmb J}_1 } \quad & \cdots \quad & {V_{N - 1} d_1^{N - 1} \pmb I } \\ \vdots \quad & \vdots \quad & \vdots \quad & \vdots \\ {V_0 d_{N - 1}^0 \pmb I } \quad & {V_1 d_{N - 1}^1 \pmb I } \quad & \cdots \quad & {\dfrac{V_{N - 1} }{\Delta \tau _{N - 1} }\pmb I + {\pmb J}_{N - 1} } \\ \end{array} \right]$$

从系数矩阵$\pmb A $的表达式可以看出,时间谱方法引入了额外的非对角项, 且随着无量纲时间步长减小,非对角项的值也就越大. 这就严重削弱了系数矩阵的对角占优特性,导致了矩阵的病态.常用的伪时间推进格式比如对称高斯赛德尔迭代均要求系数矩阵对角占优,但时间谱方法引入的额外非对角项会破坏系数矩阵的对角占优特性从而导致传统的迭代方法失效.因此, 相关研究者针对这一问题展开了研究,并提出可以通过带预处理的广义极小残值(generalized minimum residual,GMRES)方法进行伪时间推进 (Su & Yin 2010),一方面通过预处理可以改善系数矩阵$\pmb A $的性质,另一方面采用GMRES算法也可以加快迭代法的收敛速度. Mundis 和Mavriplis (2013, 2014, 2015) 针对GMRES算法的预处理进行一系列改进,可以通过并行同时计算2047个采样点, 如 表1所示.

时间谱方法的优势在于将非定常周期性流动数值求解转化为$N$个时刻的全耦合定常流场的求解.对于时间响应光滑及低减缩频率的周期性流动, 其效率可以提高一个量级左右,并且时间谱方法可以实现时间与空间的同时并行计算.以亚声速翼型简谐强迫俯仰振荡为例, 从 图5可以看出只需要5个点就能与BDF吻合,图6给出了不同采样点数情况下时间谱方法的效率. 可以看出, 在这个算例中,相比于传统时域推进法, 时间谱方法在保证精度同时可以将效率提高大约一个量级.

显示原图| 下载原图ZIP| 生成PPT

图 5   亚声速算例力矩系数$C_m$ 时间谱方法TS与BDF对比

   

显示原图| 下载原图ZIP| 生成PPT

图 6   总计算时间对比

   

时间谱方法的缺点在于随着无量纲时间步长的减小, 会带来严重的矩阵病态问题,计算效率急剧下滑.这也就决定了时间谱方法一般只能适用于减缩频率低于0.4流动的数值求解,且对于由大迎角动态失速等引起的气动力随时间响应的高频分量,基于傅里叶变换的时间谱在进行数值模拟时会有明显的吉布斯效应, 如 图7所示.

显示原图| 下载原图ZIP| 生成PPT

图 7   时间谱方法中吉布斯效应$(C_l$ 为升力系数)

   

4.2 时间伪谱方法

继时间谱方法提出之后, 各种时间伪谱方法也发展起来并应用到流体数值求解当中.这些方法与时间谱方法的本质区别在于投影的谱空间不同.谱方法投影的谱空间须具有完备的周期性边界条件, 若不具备该类方法称为伪谱法.比较典型的有时间切比雪夫伪谱法 (Dong et al. 2015) 和时间有限差分法(Leffell 2016).

切比雪夫多项式前六项图形见 图8.

显示原图| 下载原图ZIP| 生成PPT

图 8   切比雪夫多项式前六项图形

   

切比雪夫伪谱法的谱矩阵如下所示

$${\pmb D}_{\rm ch} = \tilde {\pmb D}_T \pmb {GD}_T(20)$$

其中

$$\begin{array}{l} G _{ij} = \left\{ \begin{array}{ll} \dfrac{2j}{c_i } , \quad & {\rm if} \quad {i + j} \quad {\rm even} \\ 0, \quad & {\rm otherwise} \\ \end{array} \right. \\ c_i = \left\{ \begin{array}{ll} 2, \quad & {\rm if} \quad {i = 0,N} \\ 1,\quad & {\rm otherwise} \\ \end{array} \right.\\ \end{array}$$

$$\tilde {\pmb D}_T = \left[ {{\begin{array}{*{20}c} {T_0 \left( {t_0 } \right)} \quad & {T_1 \left( {t_0 } \right)} \quad & \cdots \quad & {T_N \left( {t_0 } \right)} \\ {T_0 \left( {t_1 } \right)} \quad & {T_1 \left( {t_1 } \right)} \quad & \cdots \quad & {T_N \left( {t_1 } \right)} \\ \vdots \quad & \vdots \quad & \ddots \quad & \vdots \\ {T_0 \left( {t_N } \right)} \quad & {T_1 \left( {t_N } \right)} \quad & \cdots \quad & {T_N \left( {t_N } \right)} \\ \end{array} }} \right]$$

$$\pmb D _T = \dfrac{2}{N} \left[ {{\begin{array}{*{20}c} {\dfrac{T_0 \left( {t_0 } \right)}{4}} \quad & {\dfrac{T_0 \left( {t_1 }\right)}{2}} \quad & \cdots \quad & {\dfrac{T_0 \left( {t_N }\right)}{4}} \\ {\dfrac{T_1 \left( {t_0 } \right)}{2}} \quad & {T_1 \left( {t_1 } \right)}\quad & \cdots \quad & {\dfrac{T_1 \left( {t_N } \right)}{2}} \\ \vdots \quad & \vdots \quad & \ddots \quad & \vdots \\ {\dfrac{T_N \left( {t_0 } \right)}{4}} \quad & {\dfrac{T_N \left( {t_1 } \right)}{2}} \quad & \cdots \quad & {\dfrac{T_N \left( {t_N } \right)}{4}} \\ \end{array} }} \right]$$

切比雪夫伪谱法的基函数如下:

时间切比雪夫伪谱法的好处在于其可以应用于非周期的非定常流动,但对于周期性非定常流动切比雪夫伪谱法的效果并没有时间谱方法好.

时间有限差分法中比较经典的是时间八阶中心差分法(FDMT8), 其谱矩阵如下所示

$$\dfrac{{\rm d}}{{\rm d}t}\pmb U ^n = \sum\limits_{j_n = 1}^M {d_{j_n }^p } \left( {\pmb U ^{n + j_n } - \pmb U ^{n - j_n }} \right)(21)$$

$$d^{\rm VIII} = \dfrac{\omega N}{2\pi }\left\{ {\dfrac{4}{5},-\dfrac{1}{5},\dfrac{4}{105},-\dfrac{1}{280}} \right\}$$

时间八阶中心差分法的好处在于每个时刻的时间导数仅与相邻的8个点有关,计算时间随采样点数增加近似呈线性增长.对于时间响应有高频分量的流动仍能较好地刻画出来.但缺点在于与时间谱方法相比精度较差,时间八阶中心差分法需要更多的采样点才能达到与时间谱方法相当的精度. 图9为时间有限差分法与时间谱方法的精度对比. 图10为不同采样点数下时间有限差分法与时间谱方法每一步迭代所耗用的时间.从 图10可以看出, 随着采样点数增加,时间有限差分法每一步迭代所耗用的时间基本不变, 而时间谱方法由于是全耦合,因此随着采样点数增加, 每一步迭代所用的时间会显著延长.

显示原图| 下载原图ZIP| 生成PPT

图 9   时间有限差分法FDMT与时间谱方法TS精度对比 (Leffell 2016), $\|\pmb q_{\rm ex}-\pmb q\|_2$ 为数值解与精确解的误差二范数

   

显示原图| 下载原图ZIP| 生成PPT

图 10   时间有限差分法FDMT与时间谱方法TS每步计算时间对比 (Leffell 2016)

   

4.3 时间谱元法

时间谱元法是一种比较特殊的时间离散手段, Mundis和Mavriplis (2015)首次将该方法应用到欧拉方程的时间离散当中, 并取得了很好的效果.与其他时域配点方法不同的是,时间谱元法需要先对欧拉方程沿时间积分之后进行离散. 具体步骤如下.

对守恒量进行拉格朗日插值可得

$$ \left( {V\pmb U } \right)^{\left( m \right)}\left(\varsigma \right) = \sum\limits_{k = 0}^N {\left( {V\pmb U } \right)} ^{\left( m \right)}\left( {\varsigma _k } \right)\psi _k^{\left( m \right)} \left( \varsigma \right)(22)$$

式中,$\psi _k \left( \varsigma \right)$为拉格朗日基函数,对欧拉方程沿时间进行积分可得

$$\int_{ - 1}^1 {v\left( \zeta \right)} \left[{\dfrac{{\rm d}\left( {V\pmb U } \right)^{\left( m \right)}}{{\rm d}\zeta } + \dfrac{T^{\left( m \right)}}{2}\pmb R \left( {\pmb U ^{\left( m \right)},\dot {x}^{\left( m \right)},\tilde {n}^{\left( m \right)}} \right)} \right]{\rm d}\zeta =\pmb 0(23)$$

其中权函数取为拉格朗日基函数, 通过分步积分可得

$$\left. {\left( {V\pmb U } \right)^{\left( m \right)}\psi _p^{\left( m \right)} } \right|_{ - 1}^1 - \int_{ - 1}^1 {\left( {\left( {V\pmb U } \right)^{\left( m \right)}\dfrac{{\rm d}\psi _p^{\left( m \right)} }{{\rm d}\zeta } - \dfrac{T^{\left( m \right)}}{2}\pmb R \left( {U^{\left( m \right)},\dot {x}^{\left( m \right)},\tilde {n}^{\left( m \right)}} \right)\psi _p^{\left( m \right)} } \right)}{\rm d}\zeta =\pmb 0(24)$$

对于式(24)中积分通过高斯点进行离散, 有

$$\int_{ - 1}^1 \pmb Q {\rm d}\zeta = \sum\limits_{p = 0}^N \pmb Q \left( {\zeta _p } \right)\omega _p(25)$$

整理可得

$$\pmb \varPsi ^{\left( m \right)}\left[ {{\begin{array}{*{20}c} {V\pmb U \left( {\zeta _0 } \right)} \\ \vdots \\ {V\pmb U \left( {\zeta _N } \right)} \\ \end{array} }} \right]^{\left( m \right)} = \pmb I _\omega ^{\left( m \right)} \left[ {{\begin{array}{*{20}c} {\pmb R \left( {\pmb U \left( {\zeta _0 } \right),\dot {\pmb x}\left( {\zeta _0 } \right), \tilde {\pmb n}\left( {\zeta _0 } \right)} \right)} \\ \vdots \\ {\pmb R \left( {\pmb U \left( {\zeta _N } \right), \dot {\pmb x}\left( {\zeta _N } \right), \tilde {\pmb n}\left( {\zeta _N } \right)} \right)} \\ \end{array} }} \right]^{\left( m \right)}(26)$$

$$ \pmb \varPsi ^{\left( m \right)} = \left[ {{\begin{array}{*{20}c} {\left( {\left. {\dfrac{{\rm d}\psi _0 }{{\rm d}\zeta }} \right|_{\zeta _0 } \omega _0 + 1} \right)} & {\left. {\dfrac{{\rm d}\psi _0 }{{\rm d}\zeta }} \right|_{\zeta _1 } \omega _1 } & \cdots & {\left. {\dfrac{{\rm d}\psi _0 }{{\rm d}\zeta }} \right|_{\zeta _{N - 1} } \omega _{N - 1} } & {\left. {\dfrac{{\rm d}\psi _0 }{{\rm d}\zeta }} \right|_{\zeta _N } \omega _N } \\ {\left. {\dfrac{{\rm d}\psi _1 }{{\rm d}\zeta }} \right|_{\zeta _0 } \omega _0 } \! \! \! &\! \! \! {\left. {\dfrac{{\rm d}\psi _1 }{{\rm d}\zeta }} \right|_{\zeta _1 } \omega _1 } \! \! \! &\! \! \! \cdots \! \! \! &\! \! \! {\left. {\dfrac{{\rm d}\psi _1 }{{\rm d}\zeta }} \right|_{\zeta _{N - 1} } \omega _{N - 1} } \! \! \! &\! \! \! {\left. {\dfrac{{\rm d}\psi _1 }{{\rm d}\zeta }} \right|_{\zeta _N } \omega _N } \\ \vdots \! \! \! &\! \! \! \vdots \! \! \! &\! \! \! \ddots \! \! \! &\! \! \! \vdots \! \! \! &\! \! \! \vdots \\ {\left. {\dfrac{{\rm d}\psi _{N - 1} }{{\rm d}\zeta }} \right|_{\zeta _0 } \omega _0 } \! \! \! &\! \! \! {\left. {\dfrac{{\rm d}\psi _{N - 1} }{{\rm d}\zeta }} \right|_{\zeta _1 } \omega _1 } \! \! \! &\! \! \! \cdots \! \! \! &\! \! \! {\left. {\dfrac{{\rm d}\psi _{N - 1} }{{\rm d}\zeta }} \right|_{\zeta _{N - 1} } \omega _{N - 1} } \! \! \! &\! \! \! {\left. {\dfrac{{\rm d}\psi _N }{{\rm d}\zeta }} \right|_{\zeta _N } \omega _N } \\ {\left. {\dfrac{{\rm d}\psi _N }{{\rm d}\zeta }} \right|_{\zeta _0 } \omega _0 } \! \! &\! \! {\left. {\dfrac{{\rm d}\psi _N }{{\rm d}\zeta }} \right|_{\zeta _1 } \omega _1 } \! \! &\! \! \cdots \! \! &\! \! {\left. {\dfrac{{\rm d}\psi _N }{{\rm d}\zeta }} \right|_{\zeta _{N - 1} } \omega _{N - 1} } \! \! &\! \! {\left( {\left. {\dfrac{{\rm d}\psi _N }{{\rm d}\zeta }} \right|_{\zeta _N } \omega _N - 1} \right)} \\ \end{array} }} \right] $$

$$\pmb I _\omega ^{\left( m \right)} = \dfrac{T^{\left( m \right)}}{2}\left[ {{\begin{array}{*{20}c} {\omega _0 } \quad & 0 \quad & \cdots \quad & 0 \quad & 0 \\ 0 \quad & {\omega _1 } \quad & \cdots \quad & 0 \quad & 0 \\ \vdots \quad & \vdots \quad & \ddots \quad & \vdots \quad & \vdots \\ 0 \quad & 0 \quad & \cdots \quad & {\omega _{N - 1} } \quad & 0 \\ 0 \quad & 0 \quad & \cdots \quad & 0 \quad & {\omega _N } \\ \end{array} }} \right] $$

时间谱元法巧妙地结合了谱方法与有限元方法, 对时间进行分块求解,可以应用于任意非定常流动的时间响应.当周期非定常流动时间响应中有高频分量时,时间谱元法的效果比时间谱方法的效果好. 图11 (Mundis & Mavriplis 2015) 给出了时间谱元法的时间精度验证. 从 图11}可以看出, 当$N = 1$时精度为1.85阶, 当$N = 2$时精度为2.42阶,而当$N > 2$时精度并没有明显提升. 因此,如何提升时间谱元法的精度也需要进一步研究.

显示原图| 下载原图ZIP| 生成PPT

图 11   时间谱元法精度图

   

5 时域推进与时域配点混合方法

时间谱方法只能应用于周期性非定常流动, 而很多流场并不是周期性的,而是具有准周期的特征, 如机动旋翼、气动弹性响应等.非定常流动的时间响应中含有周期性分量和时变分量但周期性分量占主导,这类流动就叫做准周期流动. 针对这类准周期流动, Yang 等 (2010)提出了混合BDF/TS方法.将准周期流场表示为周期性分量与时变分量的叠加, 如式(27)所示

$$\pmb U ( t ) = \sum\limits_{k = - \frac{N}{2}}^{\frac{N}{2} - 1} {\widehat{\pmb U }_k {\rm e}^{{\rm j}k\frac{2\pi }{T}n\Delta t}} + \overline{ \pmb U} ( t )(27)$$

式中, 等式右端首项为周期性分量, $\overline{ \pmb U} ( t )$表示时变分量. 对于时变分量的处理以一阶精度为例, 表达式如下

$$\overline{ \pmb U} ( t ) = \varphi _{12} ( t )\pmb U ^{m + 1} + \varphi _{11} ( t )\pmb U ^m(28)$$

式中, $\pmb U ^{m + 1}$和$\pmb U ^m$表示该计算周期结束时刻点与初始时刻点.$\varphi _{12} $与$\varphi _{11} $表达式如下

$$\varphi _{12} ( t ) = \dfrac{t - t^m}{T}, \quad \varphi _{11} ( t ) = \dfrac{t^{m + 1} - t}{T}$$

对式(27)中等式右端周期性分量利用3.1节中时间谱方法进行离散,等式右端时变分量用BDF方法进行离散. 因此可得

$$\dfrac{{\rm d}}{{\rm d}t}\left( {\pmb U ^n} \right) = \sum\limits_{j_n = 1}^{N - 1} {d_n^{j_n } } \widetilde{\pmb U }^{j_n } + \varphi' _{12} \left( {t_n } \right)\pmb U ^{m + 1} + \varphi' _{11} \left( {t_n } \right)\pmb U ^m(29)$$

式中,$\widetilde{\pmb U }$为周期性分量. 且有

$$t_j = t_m + \dfrac{j}{N}\left( {t_{m + 1} - t_m } \right), \quad j = 0, \cdots,N - 1(30)$$

整理可得到最终方程形式

$$ \sum\limits_{j_n = 1}^{N - 1} {d_n^{j_n } V^{j_n }\pmb U ^{j_n }} - \left( {\sum\limits_{j_n = 1}^{N - 1} {d_n^{j_n } \varphi _{12} \left( {t_{j_n } } \right)} - \varphi' _{12} \left( {t_n } \right)} \right)V^{m + 1}\pmb U ^{m + 1} $$

$$- \left( {\sum\limits_{j = 1}^{N - 1} {d_n^j \varphi _{11} \left( {t_j } \right)} - \varphi' _{11} \left( {t_n } \right)} \right)V^m\pmb U ^m + \pmb R \left( {\pmb U ^n,\overline {\pmb x} ^n,\overrightarrow {\pmb n} ^n} \right) + {\rm {\bf S}}\left( {\pmb U ^n,\overrightarrow {\pmb n} ^n} \right) =\pmb 0 $$

$$n = 0,1,2 , \cdots,N - 1$$

混合BDF/TS方法与BDF方法的区别在于BDF是以每一个时刻点进行推进,而混合BDF/TS方法是以每一个周期进行推进,在每一个周期内用时间谱方法进行求解.由于对时变分量进行BDF离散时的时间步长为一个周期,这也就决定了混合BDF/TS方法不能保证时变分量的准确度,因此只有时变分量较小且变化平缓时该方法才适用. 图12(Yang & Mavriplis 2015) 和 图13(Mundis & Mavriplis 2013)分别为混合BDF/TS方法计算直升机爬升和气动弹性问题并与BDF进行对比,从图中可看出吻合得很好.虽然混合BDF/TS方法的计算效率相比于传统BDF方法并没有太大的提升,但混合BDF/TS方法具有良好的时空并行特性,可以实现时间与空间的同时并行计算.

显示原图| 下载原图ZIP| 生成PPT

图 12   UH60A爬升时混合BDF/TS方法与BDF计算结果对比. (a)单桨叶$x$方向力响应,(b)单桨叶$y$方向力响应, (c) 单桨叶$z$方向力响应, (d) 总升力 (1 bs=4.448 2 N)

   

显示原图| 下载原图ZIP| 生成PPT

图 13   气弹响应混合BDF/TS方法与BDF计算结果对比

   

6 时间离散格式特点及选择

前几节主要介绍了各种时间离散格式及其特点, 如 表2所示.对于非定常流动的CFD计算,根据其非定常特征要合理地选择时间离散格式才能达到更好的效果. 图14 给出了不同非定常流动特征所对应的合适的时间离散格式,为在实际工程应用中非定常流动数值计算时间离散格式的选择提供了参考.

显示原图| 下载原图ZIP| 生成PPT

图 14   不同非定常流动特征时间离散格式选取参考

   

对于周期性非定常流动来说, 要充分利用其周期性特征.因此频域谐波法和时域配点法更为适合. 对于低减缩频率流场的求解,时域配点法可以通过在时域内各时刻点的耦合求解达到高效高精度的要求,与频域谐波法相比不需要时域与频域来回转换且对程序的改动相对较小,因此应选用时域配点法. 然而对于高减缩频率流场的求解, 时域配点法的矩阵病态会导致计算效率严重下滑,此时时域配点法就失去了其原有的优势,而频域谐波法计算效率随着减缩频率变化没有太大影响,仍然能保证流场求解的高效和高精度要求, 因此应选用频域谐波法.在频域谐波法中, 对于扰动大的问题,采用非线性频域法(NLFD)能够保证更高的精度和效率. 对于小扰动问题,根据脉动载荷的动力学特征的非线性大小决定采用线性或非线性谐波法.在时域配点法中, 则需要预先判断其脉动响应是否光滑,如果脉动响应是光滑的, 则可以采用时间谱方法求解.若脉动响应不光滑或含有高阶谐波分量,传统时间谱方法需要很多时刻点且容易出现吉布斯现象,此时需要采用时间谱元法或时间有限差分法.

对于准周期流动, 其流动既有周期分量也有时变分量, 但周期分量占主导,此时频域谐波法和时间谱方法由于都是基于周期性的傅里叶变换因而失效.但是混合BDF/TS方法以及切比雪夫伪谱法等仍然适用. 在减缩频率较小时,可以选用上述方法, 虽然计算量并没有明显减少,但是其并行性相比于时域推进法有较大的改善. 当减缩频率较大时,同样是由于系数矩阵病态的原因, 选用时域推进法的计算效率会更高.

对于非周期流动, 由于其时变分量占主导, 周期性特征已经很弱.频域谐波法已然失效,时域配点法中混合BDF/TS方法以及伪谱方法等虽然有效但其计算精度和效率会明显下滑甚至不如时域推进法.因此, 对于非周期无规则脉动响应, 传统的时域推进法具有优势.

以上是针对强迫运动中减缩频率已知的流动. 而对于非定常流动中的稳定性问题,其流动频率不是事先给定的, 而是待求解的参数.无论是频域谐波法还是时域配点法都需要显式表达流动频率,因此需要对流动频率进行搜索, 而搜索的效率与给定的初值和算法有关.这方面的选择还需进一步研究.

7 新型时间离散格式在工程中的应用

近年来, 以频域谐波法、时间谱方法等为代表的新型时间离散格式发展非常迅速,并以其高精度和高效率的特点在工程中得到广泛的应用.下面以叶轮机内流、动导数计算、旋翼非定常优化设计和气动弹性动力学仿真为例,阐述上述新型时间离散格式在工程中的应用和效果.

7.1 叶轮机内部流动的数值模拟

叶轮机械内部的流动在本质上是非定常的, 由于相邻叶排之间的相对周向运动,叶轮机械内流可以看成是周期性的. 对于传统的时域推进法来说,没有利用流动的周期性特征, 其计算量庞大. 计算规模越大, 计算的过渡时间越长,从而计算耗时越多, 因此很难在工程设计中应用.

目前在叶轮机的非定常计算中以频域谐波法和时域配点法为主. Chen等(2000)应用非线性谐波法求解叶轮机内流并指出在与传统时域推进法的计算结果差别小于10%的情况下,非线性谐波法可以将效率提高20倍左右. 之后Ekici 和Hall (2007, 2008)将谐波平衡法应用到叶轮机颤振和强迫响应的计算中,效率提高了一个量级以上. He (2008, 2010)改进线性谐波法并应用于叶轮机大尺度分离流动,有很高的计算效率,并总结了傅里叶快速预测在叶轮机中的应用. Frey等(2015)采用谐波平衡法求解叶轮机内流并将线化方法和时域推进法的计算结果进行比较.Gopinath等 (2007)将时间谱方法应用到了二维和三维叶轮机的数值模拟中, {\bf图15}为用时间谱计算叶轮机叶片扭转力矩的数值计算结果, 可以看出当$N> 7$时与BDF计算结果基本吻合. Dufour等 (2012)将时域谐波平衡法应用到叶轮机叶片转子、定子的数值模拟.还有Guédeney等 (2013)、Salles等 (2014) 以及Custer等(2012)均对时域配点法在叶轮机中的应用做了相关研究.

显示原图| 下载原图ZIP| 生成PPT

图 15   用时间谱TS与时域推进法BDF计算叶轮机叶片扭转力矩的结果及对比(Gopinath & Van et al. 2007 )

   

国内有Du和 Ning (2012), Su 和 Yuan (2010), Yuan 和 Su (2009) 和王丁喜等 (Rahmati & He et al. 2014, Huang & Wang 2016) 也开展了相关的研究工作. 杜鹏程(Du & Ning 2012)对比时间推进和高维谐波平衡法对涡轮机械中流场的模拟能力,苏欣荣(Su & Yuan 2010)用时间谱隐式求解叶轮机流场并对计算方法的稳定性和无反射边界条件处理方法进行研究.王丁喜对于多排叶栅频域方法的边界条件进行改进同时也对谐波平衡法在叶轮机内流计算中的快速算法进行了研究.也有许多相关研究者对于非线性谐波法在叶轮机内流中的应用进行研究(赵军和刘宝杰 2008), 在此不做赘述.

7.2 飞行器的动导数计算

动导数常规数值计算方法是采用CFD模拟飞行器的长周期微振幅强迫简谐振动.由于预测动导数需要计算周期性非定常气动响应,计算效率是CFD在动导数方面应用的瓶颈.

时域配点法刚好适用于长周期的非定常简谐振荡问题,通过求解一个周期内几个时刻的流场重建整个周期的非定常流动,具有很高的计算精度和效率,从而在动导数计算领域得到广泛的关注和应用. Murman等 (Murman et al.2004, Murman 2012) 将谐波平衡法应用于预测Finner及SDM模型的动导数;Hassan和Sicot (011) 将谐波平衡法应用于动导数的快速预测.Ronch等(Ronch et al. 2010, 2013)对比了时域推进法、线性频域方法和谐波平衡法在数值预测动导数精度、计算效率与消耗内存方面的能力.图16 为NACA0012强迫简谐振动不同时间离散方法的计算效率对比.其中加速比为1对应的是时域推进法. LFD为线性频域法,PMB-HB对应的是Liverpool大学PMB程序中的谐波平衡方法 (Badcock et al.2000); COSA-HB对应的是Glasgow大学COSA程序中的谐波平衡方法(Bonfiglioli et al. 2009).

显示原图| 下载原图ZIP| 生成PPT

图 16   NACA0012强迫简谐振动不同时间离散格式的计算效率(Ronch et al. 2010)

   

国内的研究相对较少, 谢立军等 (2015)采用时间谱方法对高超声速HBS标模和超声速Finner标模进行动导数计算,并与时域推进法的精度和效率进行对比.陈琦等采用谐波平衡法对翼型和纯锥的非定常绕流进行数值模拟 (陈琦等 2014a),并预测带翼导弹的俯仰动导数 (陈琦等 2014b).

7.3 旋翼非定常绕流数值模拟与优化设计

对于直升机旋翼非定常绕流数值模拟与设计优化, 传统的时域推进法CFD十分耗时,很难在工程设计中应用. 而直升机在前飞过程中旋翼流场具有明显的周期性特征,这就意味着频域谐波法和时域配点法有明显的优势.

Butsuntorn等 (Butsuntorn & Jameson 2008a, 2008b; Butsuntorn 2008) 将时间谱方法应用到三维直升机旋翼前飞流场的数值模拟,通过引入涡修正将多旋翼桨叶简化为一个旋翼桨叶的流场进行求解,采用时间谱方法效率相比于传统BDF方法效率可提高一个量级. Choi (Choi & Datta 2008, Yi et al. 2015)用时间谱方法求解直升机旋翼流固耦合并分析了时间谱方法的精度和效率.Yang (Yang et al. 2011, 2012)通过混合BDF/TS方法求解直升机做准周期机动飞行时的旋翼流场模拟.

旋翼非定常设计优化也非常耗时,在定常优化设计中广泛应用的伴随方法很难应用到非定常设计优化当中.随着频域谐波法和时域配点法的发展,为伴随方法在非定常设计优化中的应用搭建了桥梁. Thomas等 (2005)提出将伴随方法与二维黏性谐波平衡方程结合计算翼型的网格灵敏度.Nadaraja和Jameson (2007) 提出结合ADjoint方法与非线性频域法对跨声速振荡翼型进行优化. Choi等(2008)将时间谱方法与ADjoint方法结合进行直升机旋翼桨叶优化的梯度求解.Woodgate分别通过时间线化法和时域谐波平衡法与伴随方法结合进行旋翼设计并进行对比(Woodgate & Barakos 2012). 图17给出了谐波平衡法与时域推进法在UH60前飞时旋翼流场对比.关于伴随方法与频域谐波法或时域配点法结合进行非定常优化设计或灵敏度分析还有很多文献(Cherif et al. 2012, Ekici & Beran 2014), 这里不作赘述.

显示原图| 下载原图ZIP| 生成PPT

图 17   谐波平衡法与时域推进法计算流场对比. (a) 二阶谐波平衡法,(b) 四阶谐波平衡法,(c) 时域推进法 (Woodgate & Barakos 2012)

   

7.4 气动弹性动力学仿真

对于翼型或机翼的气动弹性问题, 无论是线性颤振还是非线性的极限环,流场的流动均有明显的准周期特征. 在气弹的数值模拟中,如果要获取某一状态下的气弹响应, 则需要采取准周期的时间离散格式;但如果要获取颤振边界, 由于在颤振边界处幅值不变,流场的非定常响应是周期的,因此可以利用周期性的时间离散格式对颤振边界进行搜索.这也是近年来新型时间离散格式在气动弹性动力学仿真中的主要应用.

针对简化的气动力模型, Liu (Liu & Dowell 2004, Liu et al. 2007) 和Lee等 (2005)使用高阶谐波平衡法研究二元机翼的Hopf分岔行为以及预测极限环的振幅与频率,Lau和Cheung (1981)还提出了增量谐波平衡法通过引入增量项来简化方程的处理难度.

对于基于CFD和CSD技术的气动弹性仿真, 如果要获取准周期的气弹响应,需要气动和结构方程均要进行相同的时间离散. Mundis等 (Mundis {\rm \&} Mavriplis 2013, Mundis et al. 2013)利用混合BDF/TS方法求解气动弹性响应如图18 (Mundis & Mavriplis 2013) 所示, 可以看出混合BDF/TS方法与BDF2吻合的很好,同时具有更好的并行性.

显示原图| 下载原图ZIP| 生成PPT

图 18   BDF/TS方法与BDF2气弹响应计算结果对比. (a)BDF/TS方法与BDF2气弹响应计算结果对比$(M_\infty = 0.875$, $V_{\rm f} = 0.4)$, (b) BDF/TS方法与BDF2气弹响应计算结果对比$(M_\infty = 0.875$, $V_{\rm f} = 0.537)$, (c)BDF/TS方法与BDF2气弹响应计算结果对比$(M_\infty = 0.875$, $V_{\rm f} = 0.65)$

   

如果要获取颤振边界或极限环响应, 其减缩频率并没有显式给出,采用频域谐波法或时域配点法需要对减缩频率进行搜寻. Thomas 等(2002, 2004)利用谐波平衡法求解流场, 在频域内通过Newton-Raphson迭代获取极限环频率和响应.这种求解方法高效并且有很好的鲁棒性. Gopinath等 (2006)将时间谱方法应用到极限环求解中并与时域推进法和谱元法对比计算精度. 此外,Besem等 (2016) 将谐波平衡法应用于圆柱涡致振动中对锁频现象进行探究.

在国内, 刘南等 (刘南等 2016, 刘南 2016)将高阶谐波平衡法与$v$-$g$法结合进行高效的颤振边界预测,通过高阶谐波平衡法求解器对各个结构模态的小扰动简谐强迫运动进行数值模拟进而获得广义力影响系数矩阵$\pmb F ( {jw} )$, 之后通过$v$-$g$法求出对应的颤振速度和颤振频率.

可见, 时域配点法在航空工程中有较为广泛的应用且效果很好, 因此,时域配点法作为一种新型的时间离散方法有很好的发展前景.

8 结论

本文总结了近年来计算流体力学中时间离散格式的研究进展,并根据离散形式和数学思想的不同本文将非定常时间离散格式分为4大类：时域推进法、频域谐波法、时域配点法和混合方法.主要包括传统时域向后差分法、线性谐波法、非线性谐波法、非线性频域法,以及近年来发展的时间谱方法、混合BDF/TS方法、时间伪谱法以及时间谱元法等.在综述过程中比较了各方法的计算精度、效率和适用范围以及对新发展的时间离散格式在航空航天工程中应用做了概述.

近年来虽然时间离散格式得到了迅速的发展, 其中时域配点法也逐渐趋于多元化.但是时间离散格式目前仍有一些难题有待解决, 也是时间离散格式下一步发展的方向：

(1)新发展的时间离散格式需要确定并拓宽其适用范围.一些时间离散格式的适用范围和效果仍不明确, 比如时间谱元法. 同时,后期发展的方法都有局限性,对于不同的非定常流动要选择对应的时间离散格式.新型时间离散格式的计算稳定性问题也一定程度限制了其应用.因此需要寻找鲁棒性更好的迭代算法并通过预处理解决时域配点法带来的矩阵病态问题.虽然Nathan在这方面做了不少工作, 但距离理想的效果还有一定距离.

(2)对于频率未知的流动稳定性问题,无论是频域谐波法还是时域配点法都需要一个快速搜索流动频率的算法,这方面还需要进一步探索和改进. 对于准周期流动问题,频域谐波法如何应用以及时域配点法如何进一步提升计算的精度和效率也需要研究.

(3)时域配点法在气动优化设计当中的应用.对于非定常问题的优化一直以来是航空领域的一大难题.而时域配点法的出现,使得采用定常优化处理手段解决非定常优化问题成为可能.比如时间谱方法与Adjoint方法结合进行非定常优化设计.因此时域配点法在这一方面的应用有很好的前景.

(4)频域谐波法与时域配点法在多学科耦合中的应用研究. 例如,在气动弹性力学、虚拟飞行等领域中,新型时间离散格式可以精确、高效地获取所需的非定常气动力,并且十分便于非定常气动力的建模.

因此, 尽管近期时间离散格式得到快速发展, 但距离实际工程应用仍有很多工作要做.需要CFD界更多的关注和共同努力,我们相信在不久的将来新型的时间离散格式在复杂非定常流动的数值仿真和优化设计上会发挥巨大的作用.

致 谢

The authors have declared that no competing interests exist.

参考文献

文献选项

• 原文顺序
• 文献年度倒序
• 文中引用次数倒序
• 被引期刊影响因子