POD-Kriging降阶方法在串联双圆柱流场预测中的应用

王晨, 白俊强, Jan S Hesthaven, 邱亚松, 乔磊, 韩啸

(1.西北工业大学 航空学院, 陕西 西安 710072; 2.洛桑联邦理工(EPFL), 瑞士 洛桑 CH-1015)

在航空航天领域，对于复杂的带参数非定常非线性问题，基于本征正交分解方法(POD)的降阶模型(ROM)[1-2]被广泛采用以提高计算效率，其思路是通过POD对复杂非线性解空间构造一组最优降阶基以保留流场的主要特征，再将待求解空间投影到由这组基张成的线性子空间中以快速预测未知流场的低阶模型。常用的投影方法包括Galerkin投影等，但对于处理复杂的带参数非线性问题时，采用投影方法给投影系数求解带来不便，为此，目前有一批研究学者考虑采用KrigingRBF等代理模型来预测降阶模型中各POD基投影系数的降阶模型形式，如张伟伟等[3]则采用POD与径向基函数法(radial basis function,RBF)的结合预测翼型的气动响应和极限环振荡，邱亚松等[4]采用Kriging代理模型预测POD基投影系数，对由10个参数控制的翼型扰流流场进行了预测。本文考虑采用POD结合Kriging代理模型的形式构造降阶模型，预测带参数的典型流动问题——串联双圆柱构型。

在数值方法和流动机理的研究中，圆柱扰流是个经典的非定常算例。对于计算状态确定的等直径串联双圆柱，距径比L/D(圆柱柱心间距与直径的比值，L为圆柱柱心间距；D为圆柱直径)是唯一影响流动形态的参数。美国Langley实验室(LaRC)提供了L/D=3.7、L/D=1.435下basic aerodynamic research tunnel(BART)[5]和quiet flow facility(QFF)[6]2组风洞实验数据。Zdravkovich等[7-8]通过研究，总结了不同距径比下双圆柱流场的流动形态规律。由于双圆柱流动具备强烈的非定常效应，又有对流动形态起着决定性作用的几何参数，并有实验数据作支撑，是降阶模型理想的验证算例。本文试图通过降阶模型对变距径比串联双圆柱流场进行预测。

但在采用POD提取降阶基之前，降阶模型要求提供样本流场解所组成的快照数据，为此需确定双圆柱流动的数值模拟方法计算样本数据。对于串联双圆柱构型，很多学者对BART和QFF的实验工况采用不同的数值模拟方法进行了相关研究：Khorrami等[9-10]采用基于SST湍流模型的非定常雷诺平均方法(URANS)；Wessels等[11]采用格子波尔兹曼法分别对LaRC风洞实验工况进行非定常数值仿真。而相对于URANS方法，能解析非定常脉动的LES、混合RANSLES方法在模拟双圆柱流动时将更为精确，如白俊强等[12]通过采用尺度自适应模型(SAS)[13]对L/D=3.7的临界工况进行模拟，说明了SAS不仅能捕捉瞬时脉动，在时均结果的模拟精度上也优于URANS方法。本文选择SAS方法提供快照数据。

针对串联双圆柱流动，本文采用POD提取主要的降阶基模态以构造线性降阶子空间，再利用Kriging代理模型插值投影系数来封闭降阶模型以预测未知流场。通过与风洞实验及数值求解计算结果作对比，评估降阶模型在变参大分离类问题上的预测精度与效率。为保证快照的数值精度，快照数据由非定常数值求解器采集，模型选择的是基于SST[14]的尺度自适应模型(SAS)。

1 降阶基模型建立

对于带参数定常问题uθ(x)(x∈Ω),参数向量θ,系统提供N个样本uθ(i)(x)(i=1,2,…N),降阶基模型通过从样本中提取一组正交降阶基Φk(k=1,2,…,K)来预测结果:

θ(x)=u0(x)+∑Kk=1αk(θ)·Φk(x)

(1)

(Φi,Φj)=0i≠j

1i=j(i,j=1,2,…K)

(2)

αk(θ)是待求的投影系数。u0(x)通常取为快照的平均u0(x)=1N∑Ni=1uθ(i)(x),这一项的添加能够显著提高降阶基模型的预测精度。

同样对于非定常非线性问题u(x,t),降阶基模型给出类似的预测公式:

(x,t)=u0(x)+∑Kk=1αk(t)·Φk(x)

(3)

1.1 降阶基提取

(1)、(3)式中的降阶基Φk(x)可由不同的方法确定,本文选用本征正交分解(POD)方法提取Φk(x)。

POD方法提取的降阶基Φk(x)需满足快照在降阶子空间上的投影达到最大,即:

A·v=λv(5)

矩阵A是对称非负定矩阵,具有N个特征值λk和相应的特征向量vk,与之对应的基模态Φk(x)可通过k(x)=∑Ni=1vk(i)·uθ(i)(x)再对k(x)单位化后求得。为节省计算量,需要从这N个模态中选择最主要的模态来张成降阶子空间。基于j(x),j(x)〉dx=∑Nk=1λk(u′(·)=u(·)-u0(·)),Sirovich[15]提出了基于特征值λk的广义能法来选择主要模态:将特征值降序排列,定义能量百分比I(k)=∑ki=1λi∑Ni=1λi,给定能量阈值ε∈(0,1),通常取0.99,则K可通过式K=argmin(I(k)I(k)≥ε)确定,前K个特征值对应的模态就是最终获得的POD模态Φk(x)(k=1,2,…,K)。

1.2 基系数预测

Kriging代理模型[16]是通过对样本的训练建立输入、输出函数关系,对输出变量给出最佳线性无偏估计的插值方法。在本文中,Kriging模型并非直接用来预测状态变量(若将uθ)(x)(x∈Ω)直接作为输出,就需对每个网格上的状态变量建立Kriging模型,那么代理模型将会失去其原有的效率优势),而是将提取的降阶基系数αk(θ)作为输出,这样通过Kriging模型提供投影系数αk(θ),然后再通过降阶模型预测公式(1)、(3)预测待求流场。对单个状态变量来说,POD提取了K个降阶基,就需建立K个代理模型。在用POD提取基模态后,N个样本的投影系数可以通过下式得到:

αk(θ(i))=Φk(x)·(uθ(i)(x)-u0(x))

i=1,2,…,N(6)

以此作为代理模型的样本,每个Kriging代理模型的构造和预测过程如下:

设输出α,输入为p维向量θ,由1.1节得到对应的N组样本,Kriging模型对输入输出建立了如下形式的函数关系:

α(θ)=∑pi=1βifi(θ)+Z(θ)

(7)

式中，f(θ)为回归模型,Z(θ)为期望为0、方差为σ2、协方差为cov(Z(θ(i)),Z(θ(j)))=σ2R(θ(i),θ(j))的正态分布。θ(i)为第i个样本,相关函数R表征着样本两两之间在θ空间内的相关程度,其计算公式如下:

R(θ(i),θ(j))=e(-d(θ(i),θ(j)))

min:ψ(wT)=|R|1/pσ2

(θ)=f(θ)·+(R*(θ))T·R-1·(A-F)

R*(θ)=[R(θ,θ(1)),R(θ,θ(2)),…,R(θ,θ(N))]T

(10)

通过K个代理模型的建立,就可以确定参数域内任一点降阶模型的系数αk(θ),再由(1)、(3)式即可预测出对应的流场数据。

2 串联双圆柱流动简介及数值模拟方法

2.1 串联双圆柱流动简介

计算工况采用NASA双圆柱风洞实验工况(见图1):来流马赫数0.128 5(V=44 m/s),基于圆柱直径的雷诺数为1.66×105,圆柱直径D=0.057 1 m。考虑到三维非定常计算过于耗时,本文采用二维非定常计算。圆柱表面采用无滑移边界条件,远场采用无反射边界条件。物理时间步长取Δt=1.0×10-5s。

图1 串联双圆柱几何示意图

2.2 非定常数值模拟方法

本文采用双时间步,空间离散格式为二阶中心差分的有限体积求解器对双圆柱流场进行非定常模拟。

因为双圆柱的流场中存在明显的涡旋结构,而非定常雷诺平均方法(URANS)对分离流动的预测精度低,本文采用Menter在k-ωSST湍流模型的基础上发展的一种混合RANSLES模型——尺度自适应模型(SAS)方法来保证湍流预测的精度。流动方程是滤波后的Navier-Stokes方程:

∂ρ∂t+∂(ρU)∂x+∂(ρV)∂y=0

∂(ρU)∂t+∂(ρU2)∂x+∂(ρUV)∂y=-∂p∂x+

(μ+μt)∂2U∂x2+∂2U∂y2

∂(ρV)∂t+∂(ρUV)∂x+∂(ρV2)∂y=-∂p∂y+

(μ+μt)∂2V∂x2+∂2V∂y2

∂ρe+U2+V22∂t+·ρe+U2+V22U=-

(11)

式中，亚格子黏性系数μt,表征小尺度脉动对于大尺度脉动的作用,其求解仍采用了k-ωSST求解涡粘性系数的表达形式:

μt=minρkω,a1ρkΩF2

(12)

∂(ρk)∂t+ui∂(ρk)∂xi=Pk-

βkρkω+∂∂xiμ+μtσk∂k∂xi

∂(ρω)∂t+ui∂(ρω)∂xi=CωPω+QSAS-

βωρω2+∂∂xiμ+μtσω∂ω∂xi+

2ρ(1-F1)1σω21ω∂k∂xi∂ω∂xi

(13)

但在求解湍动能k和湍流比耗散率ω时,混合RANSLES通过引入比湍流模型中积分尺度更小的特征尺度,以保证模型具备识别非定常脉动的能力。其中DES类方法利用网格尺度增大湍动能(k)方程的耗散,而本文选用的SAS方法则是通过在比耗散率(ω)方程中添加由冯卡门尺度确定的源项QSAS来解析流动的脉动效应,保证模型具备捕捉非定常现象的能力。

SAS添加源项:

QSAS=ρFSASmaxζ2S2LLvk-2kσφFSST-SAS,0(14)

冯卡门尺度Lvk=κU′U″,表征流场当地速度的一阶导和二阶导的比值,它与网格尺度无关。用冯卡门尺度取代网格尺度来减小模型识别的流动特征尺度,使得SAS对网格的依赖程度要低于DES类方法。

2.3 变参数计算网格生成

对于不同距径比的双圆柱构型,在进行数值计算时,本文均采用相同的计算网格拓扑。网格采用多块结构网格,总计单元数11.2万,如图2所示:边界层内部仍采用RANS的网格要求,即物面法向第一层网格满足保证y+<1,附面层内网格物面法向增长率为1.1;附面层网格外部进行适当加密以满足非定常模拟精度。

图2 计算网格示意图

本文先以图2的网格拓扑对L/D=3.7的工况建立计算网格,当距径比变化时,通过程序生成新的双圆柱几何,再用无限代数插值(TFI)[17]动网格技术生成新构型的计算网格。

在进行计算时,CFD求解器首先将物理坐标系(oxy)转换到单位正交计算坐标系(o′x′y′)中进行计算,如图3所示。这样对于不同距径比构型采用相同网格拓扑后,尽管在物理坐标系下各网格点的空间坐标不同,但其在计算坐标系下坐标是完全相同的,本文对变距径比构型采用(1)式、(3)式建立降阶模型时采用的正是计算坐标系。采用相同的网格拓扑,和计算坐标系的使用使得降阶模型可以处理变几何参数类问题。

图3 坐标系转换示意图

3 降阶模型在串联双圆柱流动中的应用与预测

3.1 变参数时均场模拟

作为串联双圆柱构型重要的流场参数,距径比L/D的数值直接决定着双圆柱流场形态,为评估降阶模型模拟变参数流场的能力,本文试图通过一系列L/D的时均流场构造降阶模型,预测L/D=3.7.1.435这2个有实验数据的临界距比的流场,以评估模型对于预测变参类问题的精度与效率。

本文采用2.2节的非定常数值模拟方法和2.3节的网格拓扑提供20组不同L/D的时均流场作为快照数据:{1.2,1.25,1.3,1.35,1.55,1.7,1.85,2,2.25,2.5,2.75,3,3.25,3.5,3.6,3.85,4,4.25,4.5,4.75,5}。采用POD对4个状态变量(ρ,u,v,p)的快照分别提取降阶基,再通过广义能法选取主要模态:广义能百分比随基数目的变化曲线如图4所示,图中对于4个状态变量,前5个基模态占据了几乎所有的广义能;给定99.5%的能量阈值,本文最终为ρ,u,v,p分别选取5,5,4,4个降阶基。然后对这5+5+4+4=18个POD基的权重系数分别建立Kriging代理模型来预测POD基系数,这样对于任意L/D∈[1.25,5],结合降阶基和代理模型,就能预测出L/D在1.25到5范围内任一距径比的时均流场。

图4 能量百分比随基模态数目的变化曲线

1) 预测精度评估

本文重点考虑L/D=3.7和1.435的预测精度。

L/D=3.7和1.435时风洞实验、数值求解器、降阶模型三者得到的空间流线分布如图5所示。

图5 时均流线对比图

其中实验流线图来源于文献[10]。L/D=3.7时,数值模拟和降阶基预测的下游圆柱空间背风面的分离区域大于实验所探测的;L/D=1.435时,实验数据中两圆柱间是个整体的分离涡形式,而数值结果则在圆柱间预测出上下2个相互接触的漩涡结构。尽管与实验的时均流线存在差异,降阶模型与求解器的流线包括圆柱背风面分离区的大小和形态几乎完全一致,这说明时均流场的差异主要是由于数值求解方法所致,而并非降阶模型预测误差所造成。

为进一步检验降阶模型的预测精度,考察降阶模型预测的时均流场和数值求解之间的误差分布。定义降阶模型预测值与数值求解器计算值的绝对误差Δ(x)=‖θ(x)-uθ(x)‖,相对误差Δr(x)=‖θ(x)-uθ(x)‖uθ(x)。图6给出ρ,p在L/D=3.7,1.453时的相对误差分布,考虑到距径比对远离物面处流动的影响小,此时降阶模型误差小,故这里仅提供双圆柱附近流场的误差分布。在整个考察区域内降阶模型预测的相对误差处在10-3量级,说明降阶模型对ρ,p实现了精确的预测。至于速度场,在本算例中,圆柱的背风面是流动滞止区,在这个区域内的速度相对其他区域数值很小,使得降阶模型在该区域某些地方相对误差过大,故本文考虑速度场的绝对误差分布,如图7所示。图中流场各处速度的绝对误差始终处小于10-2,而在图5中降阶模型又能较精确地模拟两工况的时均流线,因此本文得出降阶模型对速度场、密度和压强均实现了全流场精确预测的结论。

物体表面压强分布是气动计算与预测的关注点之一,降阶模型预测的圆柱表面压力分布如图8、图9所示。它与求解器的计算结果除在驻点处外几乎完全重合,在L/D=3.7预测的两圆柱压力分布及L/D=1.435时上游圆柱的压强系数分布与实验也保持一致,但L/D=1.435下游圆柱与实验则相差较大,这和未能模拟出实验中此距径比下两圆柱间的分离形态相吻合。这种差距一方面是因为数值模拟方法的精度仍有待提高,另一方面是由于数值计算时没有考虑风洞实验的某些条件和外在因素所致。

图6 时均流场 图7 时均流场速度绝对误差分布图

图8 L/D=3.7圆柱表面压力分布对比 图9 L/D=1.435圆柱表面压力分布对比

通过时均流线、误差分布、表面压力分布的对比,本文的降阶模型在对20个样本提取降阶基后,能精确地模拟L/D=3.7和L/D=1.435的时均流场。

2) 计算效率评估

除预测精度之外,计算效率是评估降阶模型性能优劣的另一个重要因素,降阶模型各阶段计算耗时如表1所示:

①降阶模型对20个样本生成降阶子空间耗时22.42 s,对POD基系数建立代理模型时花费1 452 s,这两部分构成了整个训练阶段;

②在降阶模型建立后,对于单点流场,降阶模型仅需2.3 s就能较精确地预测出时均结果。而采用求解器计算时均结果至少需要两万非定常时间步,需耗时19 382 s。因此,本文为双圆柱流场构造的降阶模型利用训练阶段的基函数和代理模型能快速又精确地预测变距径比时均流场,其在计算效率上带来的收益非常可观。

表1 计算效率评估

3) 时均流场变化趋势分析

在L/D∈[1.25,5]范围内均匀选取750个测试点,用建立的降阶模型提取时均流场后,本文对変距径比串联双圆柱时均流场得出以下规律(这750个测试点的网格数据是通过对空间坐标另外建立POD-Kriging降阶模型得到的):

①在L/D∈[1.25,2.815]范围内,上游圆柱背风面的分离涡占据着圆柱间的整个区域,如图10a)所示;

②在L/D∈(2.815,5]范围内,上游圆柱背风面的分离涡仅分布在两圆柱间的部分区域,其与下游圆柱背风面并未接触,如图10b)所示。

图10 时均流场形态

3.2 定参非定常模拟

3.1节对变距径比时均流场进行了分析,在本小节中,针对L/D=3.7这个临界工况的非定常流场建立降阶基,以考察降阶模型对非定常问题的预测精度。对该工况计算5万非定常时间步后,令此时t=0 s,再在t∈[0,0.5]范围内均匀选取501个时刻的瞬时流场作为快照数据。采用POD对各状态变量(ρ,u,v,p)的快照分别提取降阶基,能量百分比随POD基数的变化曲线如图11所示,利用广义能法为ρ,u,v,p分别选取(15,15,11,13)个降阶基。

图11 能量百分比随基模态数目的变化曲线

取t=0.101 5 s、0.206 5 s、0.336 5 s、0.473 5 s这4个时刻作为测试点,显然这4个测试点均不属于样本集。图12为速度的绝对误差分布,在考察的4个时刻降阶模型预测的绝对误差始终保持在10-3量级,说明了本文建立的降阶模型在t∈[0,0.5]范围内对瞬时速度场具有较高的预测精度。图13为测试时刻的瞬时涡量图,左图中本文采用的数值模拟方法能在一定程度上反映出捕捉流场的非定常特征,而对于降阶模型来说,非定常效应往往难以精确模拟,右图中本文采用的降阶模型预测的涡量分布与左图数值计算的几乎完全一致,证明了其对非线性、非定常问题具有较高的预测精度。

图12 测试点速度绝对分布图13 瞬时涡量图

4 结 论

本文通过采用POD提取主要的降阶基模态，对流场解构造了线性降阶子空间，再利用Kriging代理模型插值投影系数来封闭降阶模型，对变距径比串联双圆柱流场进行预测，得到以下结论：

1)在与实验对比后，本文采取的非定常求解方法能较精确预测串联双圆柱非定常流动；

2)对于变距径比构型，降阶模型能较精确预测出非样本点的时均流场，相对于全阶模型计算效率提高69 470倍，其与实验数据的误差主要源于数值求解的误差；

3)通过建立的降阶模型，本文总结了双圆柱时均流场随距径比的变化趋势；

4)在特定工况下，POD-Kriging降阶模型同样能精确预测给定时刻的瞬时流场。

考虑到传统的POD-Galerkin方法在多参数问题中构造不便，本文并未给出其与传统的POD-Galerkin方法的对比而只给出了其与全阶模型精度及计算效率的对比。另外，由于三维非定常计算非常耗时，为保证取样效率以印证降阶模型对变参串联双圆柱流场的预测精度与效率，本文采用二维非定常计算，但由于忽略展向流动，这在一定程度势必会损失预测精度，后续工作将会针对三维数值模拟结果建立适应于带参数非定常问题的降阶基模型，此时相对于二维流场预测，降阶模型在计算效率上的收益预计将会更为可观。

针对变参问题的降阶模型在实际应用如气动设计、流固耦合等方向中将拥有比较广泛和诱人的应用前景，后续研究将集中在高维参数空间的强非线性非定常问题中提高降阶模型鲁棒性和自适应取样这两项工作。

参考文献：

[1] Thomas J P, Dowell E H, Hall K C. Three-Dimensional Transonic Aeroelasticity Using Proper Orthogonal Decomposition Based Reduced Order Models[J]. Journal of Aircraft, 2003, 40(3): 544-551

[2] Du J, Fang F, Pain C C, et al. POD Reduced-Order Unstructured Mesh Modeling Applied To 2D and 3D Fluid Flow[J]. Computers & Mathematics with Applications, 2013, 65(3): 362-379

[3] Zhang W, Kou J, Wang Z. Nonlinear Aerodynamic Reduced-Order Model for Limit-Cycle Oscillation and Flutter[J]. AIAA Journal, 2016, 54(10): 3304-3311

[4] Qiu Yasong, Bai Junqiang. Stationary Flow Fields Prediction of Variable Physical Domain Based on Proper Orthogonal Decomposition and Kriging Surrogate Model[J]. Chinese Journal of Aeronautics, 2015, 18(1): 44-56

[5] Jenkins L N, Khorrami M R, Choudhari M M, et al. Characterization of Unsteady Flow Structures around Tandem Cylinders for Component Interaction Studies in Airframe Noise[R]. AIAA-2005-2812

[6] Jenkins L N, Neuhart D H, McGinley C B, et al. Measurements of Unsteady Wake Interference between Tandem Cylinders[R]. AIAA-2006-3202

[7] Zdravkovich M M. Flow around Circular Cylinders Volume 2: Applications[M]. Oxford University Press, 1997

[8] Zdravkovich M M. Review of Flow Interference between Two Circular Cylinders in Various Arrangements[J]. Journal of Fluids Engineering, 1977, 99(4): 618

[9] Khorrami M R, Choudhari M M, Lockard D P, et al. Unsteady Flowfield around Tandem Cylinders as Prototype Component Interaction in Airframe Noise[J]. AIAA Journal, 2007, 45(8): 1930-1941

[10] Lockard D P, Choudhari M M, Khorrami M R, et al. Aeroacosutic Simulations of Tandem Cylinders with Subcritical Spacing[R]. AIAA-2008-2862

[11] Brès G, Wessels M, Noelting S. Tandem Cylinder Noise Predictions Using Lattice Boltzmann and Ffowcs Williams-Hawkings Methods[R]. AIAA-2010-3791

[12] 白俊强,王晨,张扬. 一种基于冯卡门尺度的湍流模式在模拟稳态和非稳态流动问题中的应用[J]. 工程力学, 2014, 31(11): 39-45

Bai Junqiang, Wang Chen, Zhang Yang. Application of a Turbulence Model Based on Von-Karmen Length Scale in Steady and Unsteady Flow Simulation[J]. Engineering Mechanics, 2014, 31(11): 39-45 (in Chinese)

[13] Menter F R. A Scale-Adaptive Simiclation Modelusing Two-Equation Medels[R]. AIAA-2005-1095

[14] Menter F R. Two-Equation Eddy-Viscosity Turbulence Models for Egineering Applications[J]. AIAA Journal, 1994, 32(8): 1598-1605

[15] Sirovich L. Turbulence and the Dynamics of Coherent Structures[J]. Quarterly of Applied Mathematics, 1987, 45(3): 561-571

[16] Simpson T W, Mauery T M, Korte J J, et al. Kriging Models for Global Approximation in Simulation-Based Multidisciplinary Design Optimization[J]. AIAA Journal, 2001, 39(12): 2233-2241

[17] Li J, Huang S, Jiang S, et al. Unsteady Viscous Flow Simulations by a Fully Implicit Method with Deforming Grid[R]. AIAA-2005-1221