英语原文共 25 页

数值和实验数据的动态模态分解

Peter Schmid

描述流体的耦合特性对我们理解流体动力学和输运关系至关重要。在本篇文章中，我们将介绍一种能够从流场中提取动态信息的方法，这种方法的基础数据要么通过直接的数值模拟获得，要么通过在物理实验中的观测或者测量来获得。我们提取的动态模态能被理解为是全局稳定模型的概括，能够用来描述被捕获在数据序列中的最根本的物理机制，或者将大规模的问题投射到一个自由度小得多的动力系统上。在相关动力学描述的流场子域上的集中，这种方法允许将复杂的流动分解到局部不稳定现象的区域，并进一步说明了该方法的灵活性，正如在空间框架内描述动力学一样。这个方法应用的例子，包括平面槽流、二维空腔流动、柔性膜后尾流和两柱间射流。

1、介绍和动机

在数值模拟和物理实验中，精确描述复杂几何中的干扰行为，以及提取和量化这种行为对计算算法提出了很大的挑战。与此同时，在许多的工业应用中，如燃烧室内的流动，我们可以从对根本的转变和不稳定机制的更深入了解中得到很大的好处。

对于复杂几何中的流动进行全局稳定性分析变得越来越普遍，但是由此产生的巨大稳定矩阵对于计算资源有着巨大的压力。对于简单问题，我们可以选择直接法，但是选择这个方法运算的代价是十分昂贵。因此，我们必须使用迭代方法来提取全局稳定模型。Arnoldi方法（参考Edwards et al. 1994）在这方面特别成功。它通过投影到低维的Krylov子空间来得到高维系统矩阵的近似。这样，整个系统的主要特征值（以及相应的特征向量）能够被相当有效地计算出来。我们可以通过位移、映射、重新开始、锁定和清除技术来提高这种方法收敛性，但是在大多数情况下，由于需要额外的矩阵求逆，这种计算的完成通常是以牺牲计算效率来实现的。

通过全局稳定分析来提取动力学特征已经是一个几乎专门应用于数值模拟的工具。这是因为数值模拟各自的算法需要基础流场的系统矩阵来构件一个（人工）流场序列，而数值方法的收敛依赖于这些流场序列。在物理实验中，这些系统矩阵是得不到的，然而在数值模拟中，右边的子程序调用是可以直接自适应的，即满足系统矩阵的，而在实验中却满足不了这样的要求。相反，在实验设置中，我们能得到的唯一的输入就是流场本身，不论是粒子图像速度测量（PIV）的形式，还是被动示踪剂可视化的形式。因此，如果想要在实验数据中辨认出拟序结构，算法的设计就只需要依赖这些测量值。

识别拟序结构的一种常用方法是正交分解法（POD）（Lumley 1970;Sirovich 1987;Berkooz, Holmes amp;Lumley 1993）。这种方法能够从流场快照（snapshots）中提取信息，因此同样能够适用于实验数据（参考 Herzog 1986; Bonnet et al. 1994;Noack et al. 2003）。该方法通过将从流场快照中计算得出的空间相关矩阵对角化，来确定流场中具有最大能量的结构。但是使用这个方法有两个默认的缺陷：（1）在所有情况下，能量也许不是最正确的测量方式来对流体结构进行排列；（2）由于此方法选择了二阶统计量作为分解的基础，有价值的相位信息丢失了。第一个缺点已经得到广泛承认，Noack等人（2008）对动态相关零能量模态的存在进行了解释。当我们选择更侧重于流场特定分量的权函数或者更活跃的流动区域，可以改善我们对总动能扰动能的关注。然而，第二个缺点更难去克服。双正交分解（BOD）（参考Aubry 1991;Hemon amp;Santi 2007）通过对流场快照矩阵进行奇异值分解，同时产生时间和空间结构，我们称之为“chronos”和 “topos”；当计算相关张量的时候，通过选择时间或者空间的平均值，我们可以得到相同的结果；或者分别地使用经典的POD方法（Lumley 1970）与snapshot POD方法（Sirovich 1987）也可以得到。时间结构（chronos）代表了空间上平均的时间相关矩阵的特征向量，而空间结构（topos）表示了时间上平均的空间相关矩阵的特征向量。在任何一种情况下，产生二阶统计量的平均过程都会造成信息的丢失，这些信息在对流场快照中包含的动态过程进行分类时可能非常重要。

在这篇文章中，我们将提出一种仅仅基于流场快照的方法（参考Schmid amp; Sesterhenn 2008），但是它仍然会产生精确描述流体运动的流动结构。对于线性化流（即稳定基流的小扰动流），我们所提取的模态与整体稳定性分析的结果相当；对于非线性流动，结果产生了与底层流动近似的线性切线结构，并描述了表达数据序列中主要动态行为的流体单元。从经典的线性全局稳定性分析出发，我们将提取的流动结构称为动态模态，将分解技术称为动态模态分解。这一项技术是基于Koopman对于非线性动态系统的分析（参考Lasota amp;Mackey 1994; Meziacute;c 2005），该模型最近被Rowley等人（2009）应用于横流中射流的大规模模拟。我们将指出新方法与常用分解方法之间的联系和不同，这将会有助于用熟悉的描述拟序结构的技术正确地看待这一新方法。我们将验证此方法并将其应用于三种不同的流动情况：其中一个是基于用Navier-Stokes代码生成的数值数据，后面两个是基于从实验流动测量中提取的数据。我们强调在任一种情况下，只利用流场来提取相关的动态信息。

2、利用流场数据进行动态模态分解

动态模态分解（DMD）的推导将会与迭代方法中用于计算线性特征值或者其他线性代数问题解的参数密切相关。此外，熟悉正交模态分解（Lumley 1970; Sirovich 1987）和双正交模态分解（Aubry 1991; Hemon amp; Santi 2007）的读者将会注意到共性；然而，重大的差异也确实存在，在下文我们将会指出。

如前所述，我们的目标是寻找一种同样适用于实验和数值流场数据的方法；因此，我们尝试使用一种仅仅依赖收集到的输入数据而忽略底层系统矩阵任何信息的无矩阵公式。实际上，我们避免使用基于模型的方法来提取动态信息，相反的，我们更专注于基于数据的过程。

2.1大体的描述

我们将从通过提取直接数值模拟或者实验数据中任意一个采样收集的流场的一般描述开始。为了消除固有测量噪声，对实验数据进行预处理是必要的。数据应该以流场快照序列的形式表示，由矩阵给出，

（2.1）

式中代表第i个流场。在上述定义中，下标1表示序列中的第一个元素，上标N表示序列中的最后一项，即分别表示矩阵的第一列和最后一列。我们进一步假设一组由恒定采样时间分隔的有序数据序列；采样时间在两个连续流场快照之间的选择将在做进一步讨论。

第一步，我们假设一个线性的映射A将流场和随后的流场联系起来，即

（2.2）

我们假设这个映射在整个采样区间是近似相同的。如果流场起源于一个非线性过程，那么这个假设就等于线性切线近似。对于缓慢变化的系统，多尺度参数能够为上述假设提供基础。对于纯线性过程的特殊情况，假设一个常数映射不会引起近似。无论如何，流场快照之间的常数映射的假设将允许我们将流场序列表示为Krylov序列。（例如参考 Greenbaum 1997;Trefethen amp; Bau 1997）。

（2.3）

然后，我们的目标是提取由A描述的基于序列的动态过程的动态特性。（本征值，本征向量，伪本征值，能量放大，共振行为等等）

随着流场快照数量的增加以及由给出的数据序列捕获底层物理过程的主要特征，可以合理地进行假设，在流场快照的临界数量之外，（2.2）式给出的向量是线性相关的。换句话说，在数据序列中增加更多的流场并不会改善所形成的向量空间。当这个限制达到时，我们能够将向量表示为先前流场向量的线性组合，这些流场向量是相互线性无关的，例如

（2.4）

或者用矩阵的形式

（2.5）

其中，r作为残余的向量。我们继续跟随Ruhe（1984）的想法，并且这样改写：

（2.6）

或者用矩阵的形式

（2.7）

其中作为第（N-1）个单位向量。

简单的计算表明，矩阵S是这个样子的:

（2.8）

它的子对角线元素反映了这样一个事实，通过我们的设计，当时，向量的第i列和向量 的第（i 1）列相同。S中唯一的未知数是系数{}，它们是先前样本线性表示最后一个样本的线性系数。

S的特征值近似于A的一些特征值。著名的Arnoldi方法（参考例如Greenbaum 1997;Trefethen amp; Bau 1997）和上述的分解方法密切相关，但是互相正交的向量，我们可以得到分解形式的Hessenberg矩阵。同样的，矩阵H的特征值也近似于矩阵A的一些特征值。实际上，用Arnoldi方法将矩阵A还原成Hessenberg形式，并不是简单地通过的QR分解来实现的，而是通过一系列对连续Krylov子空间的投影来完成的。这将产生一种稳定的算法，但是对于这些投影，矩阵A必须是可得到的，这使得经典的Arnoldi方法对我们没有吸引力。相反，为了得到一种完全基于流场的数值技术，并且同样适用于实验数据和大规模的数值模拟，我们需要解决算法稳定性（收敛性）较差的问题。

S矩阵的计算过程如下：给定数据序列的最后一个元素可以表示为序列前几个元素的线性组合，如（2.5）所述，对于一个满秩矩阵他的最小二乘解为：

（2.9）

以作为数据序列的经济规模的QR分解。向量a的（N-1）个元素形成了矩阵S的最后一列。

即使上述的分解方法基于矩阵S在数学上是正确的，而且它常被用来证明完全Arnoldi 方法的收敛性质（Ruhe 1984），但是实际上它会产生一个病态的算法，这个算法通常不能够提取超过第一阶，或者前两阶的主要动态模式。当数据来源于实验并受到噪声和其他不确定因素污染时，这一点尤其正确。基于这个原因，我们选择一个更加正确的实现方法，它通过相似变换得到一个与S矩阵相关的完整矩阵。我们通过使用奇异值分解的预处理步骤使数据序列来实现鲁棒性，即。将上述奇异值分解结果代入式（2.7），经过整理，我们可以得到。通过辨认矩阵U所包含的数据序列的正交模型，上述操作相当于一个线性算子A在POD基上的投影。（更多关于这些联系的细节会在2.2节中给出）。除了对A的低纬度表示进行更稳定的计算之外，这个方法的另一个优点是，通过sum;的非零奇异值给出的有限投影基U的限制来解释在数据序列中的非满秩问题（或者以高于规定阈值的奇异值表示）。在这篇文章中我们已经使用了基于矩阵的方法得到了满意的结果。

在某种程度上，从矩阵中提取模态结构的方法类似于从标准Arnoldi方法的海森堡矩阵H的特征向量中恢复全局模态。在我们的例子中，我们有如下的动态模式表示：

（2.10）

其中以是的第i个特征向量，例如，同时U是流场快照序列的右奇异向量。上述分解方法，无论是基于伴随矩阵的数学形式，还是基于满秩矩阵的方法，都只能从一系列流场数据中提取拟序结构。

2.2与其他分解方法的关系

表达式（2.7）清楚的显示了与其他常用分解技术的关系。正如前面所说的，Arnoldi方法通过连续投影到正交Krylov序列上，将矩阵A降维成一个更低维度的海森堡阵。在我们的例子中没有采取正交化的步骤；相反地，我们直接从流体快照中生成了一个更低维度的系统矩阵。这个特性的代价是降低了算法的稳定性和收敛性，但反过来，它又允许了一种无模型方法的应用。

大气或海洋数据分析中常用的一种方法，是将数据分解成主相互模式（PIPS）或者是主振荡模式（POPs），同时利用了一种统计方法，将相干模式在一个时间步长内，从时间互相关的比率恢复成初始流场的自相关状态（(Hasselmann 1988; von Storch et al. 1995）。首先假设，我们通过一个线性过程来获得数据矩阵，即

（2.11）

通过两边同时右乘数据向量的转置，并在时域上取平均值，

资料编号：[4903]