基于POD方法的二维方柱低雷诺数绕流流场分析研究

王掩刚, 陈俊旭, 先松川

(西北工业大学 动力与能源学院, 陕西 西安 710072)

在计算流体力学的应用中,通常求解数学模型方程组的计算规模庞大、维数较高,对计算机的容量和速度提出了极高的要求。因此如何在保证足够精度的条件下对高维或无穷维流体动力系统进行降维建模成为了计算流体力学领域研究的热点问题,国内外研究学者为此进行了广泛研究。

本征正交分解(proper orthogonal Deco-Mposition,POD)方法作为一种有效的降维手段,被广泛运用于流体动力系统的降维建模中。在国外,Favier等[1]通过POD降阶模型对圆柱绕流尾迹,以及某翼型的表面分离流动进行了分析研究,结果表明该方法可以很好地模拟所研究的流场区域;Siegel等[2]从数值模拟中提取出圆柱绕流瞬态流场在短时间内的快照(snapshots)集合,运用POD方法对该流场进行了分析研究,证实了该方法可以准确的模拟和评估瞬态流场结构;Scarano等[3]则通过PIV技术和POD方法研究了入射角对二维的方柱绕流流场的影响,结果表明该方法可以很好的展现流场的旋涡脱落过程以及入射角的影响。在国内也有学者对POD方法的进行了一定的相关研究。张震宇[4]利用POD原理设计了一种针对风力机翼型动态失速时变过程的辨识方法, 结果表明该降阶模型方法能够以明显降低的计算量精确辨识翼型的浅失速情况;倪振华、江棹荣等[5]基于POD近似的时间与空间分解来预测未知区域的风压时间序列,研究表明该方法能够有效合理地反映出风压场的时间序列。

综合国内在POD方法的相关研究来看,在非定常流场方面的应用研究相对较少,因此本文基于POD方法对二维方柱尾缘特征区域的非定常流场进行了研究,对求得的POD模态以及时间系数进行了初步分析,并用POD模态对原始流场进行了重构,证实了该方法的可行性,为后续POD方法在流场方面的研究应用提供了一定的参考。

1 计算模型及POD模型的构建

1.1 计算模型

计算采用二维方柱为研究对象,方柱的边长D

=1 m,计算域为11D×30D,上游边界距方柱中心4.5D,下游边界距方柱中心25.5D,上、下侧边界分别距方柱中心5.5D。在大量网格校验的基础上,计算网格采用H型网格,网格密度选取370×280,并在方柱周围以及下游尾迹区域加密,以准确捕捉其流动细节。

计算工质为空气,雷诺数Re=100,采用层流模型;选用压力基求解器对各控制方程进行求解,压力速度耦合算法选用SIMPLE算法,梯度项选用least squares cell based方法;动量项选用二阶迎风格式。上游边界指定为速度入口边界条件u=0.001 46 m/s,v=0 m/s;下游边界指定为压力出口边界条件;方柱表面指定为壁面无滑移边界;上下边界指定滑移边界条件u=0.001 46 m/s,v=0 m/s。

1.2 POD简介及模型的构建

本征正交分解(POD)方法由LumIey[6]于1967年首先引入到湍流相干结构研究中。该方法的基本思想是把原来在时间域及空间域上连续物理量的场,写成一个只与时间相关的函数和一个只与空间相关的函数的展开式序列,且它们在均方意义上是最优的,在展开式中只需要少量的项数就可较准确地描述该物理过程,从而可以提供具有足够高精度而自由度又较少的低维模型,以降低计算维数、节省计算时间和内存。本文将以上述方柱绕流的非定常计算模型为基础,运用POD方法来构建模型并进行研究分析。

图1中所示的点划线矩形区域即为本文进行POD分析的特征区域,区域范围为X:-9.5～-4.5;Y:-2～2。选取该区域的原因在于它是方柱绕流流动的主要特征区域，能够反映出尾迹脱落涡的整个演化过程。

图1 POD分析的特征区域

针对以上特征区域,在流场稳定周期性波动的某段时间内,本文选取了该区域中N=60个时刻的瞬态速度场快照(snapshots)集合U(x,y,ti)(i=1,2,…,N)用于建立降维的POD空间。定义这组样本的平均值和脉动量分别为:

(1)

(2)

根据POD方法的基本思想,目标是将脉动量V(x,y,ti)分解为空间模态φi(x,y)与时间系数ai(t)的函数即:

(3)

实际上,求解空间模态φi(x)等价于求解以下最大值问题:

(4)

且满足

(φ,φ)=‖φ‖2=1

(5)

式中：(·,·)和‖·‖分别表示内积空间Ω上的L2-内积和L2-范数。

利用变分法,上述最大值问题可转化为以下特征值问题:

(6)

式中:

i=1,2,…N为Vi的自相关函数,也称POD的核;可以利用原有函数空间快照脉动量的线性组合来表示空间模态,即

(7)

将上(7)式代入(6)式中得以下特征值问题:

(8)

Φ=span(φ1,φ2,…,φN)

(9)

基于上面所求的M个POD模态,原速度场的降维模式可以表示为

(10)

以上便得到了由M个POD基函数扩展成的降维空间,使其在降维空间里求解。用Galerkin逼近将模式方程投影到POD基函数扩张成的降维空间中,得到求解时间系数ai(t)的常微分方程组,便可求得时间系数ai(t),最终得到POD的低阶模型。

2 数值模拟结果及POD分析

2.1 数值模拟结果

图2和表1是通过上述计算模型得到的方柱绕流非定常数值模拟的结果。

图2 升阻力系数随无量纲化时间的变化曲线

表1 方柱的计算结果与参考文献[8-9]的对比

2.2 POD结果分析

图3 各阶POD模态的能量分布

由图3可以看出,所计算出的特征值均成对出现,而且一对特征值所占的能量基本相同。其中第一对POD模态所占总波动能量的比例分别为51.62%与44.08%,而剩下的POD模态分别所占总能量比例均不大于2%。第二对POD模态则与相对较低能量的高次谐波相关[7],而更高阶的POD模态则是包含了流体运动中的更高次谐波。

根据之前E(k)的定义可以得出前2阶POD模态占总能量的比例为95.7%,一大部分的流动动能均包含在前2阶模态当中;前4阶POD模态所占总能量的比例已高达99.4%。由此可以看出前4阶POD模态已经完全可以抓住流场流动的主要特征,因此本文仅给出前4阶POD模态的速度基函数图像,并对前4阶POD的相关结果进行了初步分析以及对非定常流场的重构。

图4分别给出了POD分析中所占能量比例较高的前2个模对相应的时间系数(a1,a2)、(a3,a4)随无量纲化时间的变化曲线。从图中可以看出:时间系数基本都是呈正余弦曲线的变化趋势,且每个模对中2个时间系数的频率及振幅基本相同,其中时间系数a1、a2的变化周期为旋涡脱落周期的1/2;由图5a)～图5c)还可以得到每个模对的2个时间系数相位差为π/2,幅值随着模对阶数提高而减小,但第一个模对的时间系数幅值占主导地位。据此,可以推测:对于方柱绕流的尾迹旋涡流动,随着旋涡的生成、发展和消亡,在微尺度上存在着多个相同频率和幅值、但旋转相位有差异的对旋涡对(counter-rotating vortex)。能量最大的第一阶涡对决定了尾迹流动的波动频率和幅值,低能量的高阶涡对则影响着流场的细微结构。

图4 前4阶时间系数随时间t′的变化

为了清楚地看到不同模对时间系数的相互关系,图5给出了POD分析中前3个模对相应的时间系数ai在60个时刻点相关的利萨如图形。不同的时间系数合成轨迹为封闭的图形,且频率比满足如下关系式:

式中：fx、fy分别为2个信号的频率,nx、ny分别为利萨如图形的外切水平线和垂直线的切点数。从图5d)～图5f)可看出:a1、a4的频率比值为1∶2;a1、a5的频率比值为1∶3;a4、a6的频率比值为2∶3。应用利萨如图形在不同频率信号叠加时其形状与相位关系可以得出:时间系数a1与a4、a5的等效相位差均为nπ/2(n=1,3,5…),而时间系数a4与a5的等效相位差nπ/4(n=1,3,5,…)。

通过以上分析,得到了前3个模对中的6个时间系数之间的频率及其相位关系。对于一个复杂的旋涡流场来说,各阶旋涡波动共存,并且相互干涉。通过获取各阶波动的关联关系,如果能够对相对较低能量的波动实施干扰,控制流动细微结构,从而实现旋涡运动的宏观改善,这对于复杂流动的主动控制技术有一定的参考价值。

图5 时间系数的利萨如图形

2.3 基于POD结果的流场重构

基于上述分析,对于本文所研究的低雷诺数条件下的方柱绕流,前4阶POD模态表征了绝大部分的流场能量,为了做进一步验证,下文通过前4阶POD模态对流场进行了重构,并比较了不同阶数的POD模态对流场的重构效果。

图6和图7给出了应用POD方法构造的速度分量前四阶POD基函数等值线云图。由图中可以看到成对的相似模式和空间变化,这与方柱绕流尾部典型旋涡脱落的对流特点有关。对于流向速度u而言,前2阶POD模态在横向出现了正负交替的反对称涡核结构,流向则是正负交替的涡核;3阶和4阶POD模态在横向呈现对称结构,流向则出现正负交替的涡核。对于横向速度v而言,前2阶POD模态显示出横向完全对称,流向交替出现的涡核结构;3、4阶POD的模态则在横向表现为反对称,流向交替出现的涡核结构。这一结果与Van Oudheusden B W等[3]、Ben Chiekh M等[7]分析结果基本一致。

图6 流向速度u的前4阶POD基函数等值线云图

图7 横向速度v的前4阶POD基函数等值线云图

图8 原始流场与不同阶数POD重构结果对比

3 结 论

本文的研究结果表明:

1) 应用本文的数值分析手段,捕捉到了方柱绕流非定常流动特征,且与公开文献结果吻合较好。

2) 对于本文的研究对象,前4阶POD模态所占波动总能量为99.4%,可以用于较准确描述其非定常流场特征。

3) 各阶POD模态的时间系数之间都有明确的频率与相位关系,其中能量最大的第一对POD模态对方柱尾迹流动的波动频率和幅值起决定作用,这对于流动分析控制有一定的参考价值。

4) 在所研究的雷诺数下,前4阶POD模态就可以很好地重构流场,这对以后低阶模型的建立有一定的指导意义。

参考文献：

[1] Favier J, Cordier L, Kourta A. Accurate POD Reduced-Order Models of Separated Flows[J]. Physics of Fluids, 2007, 8(3): 259-265

[2] Siegel S, Cohen K, Seidel J, et al. Short Time Proper Orthogonal Decomposition for State Estimation of Transient Flow Fields[C]∥43rd AIAA Airspace Sciences Meeting and Exhibit, 2005: 296

[3] Van Oudheusden B W, Scarano F, Van Hinsberg N P, et al. Phase-Resolved Characterization of Vortex Shedding in the Near Wake of a Square-Section Cylinder at Incidence[J]. Experiments in Fluids, 2005, 39(1): 86-98

[4] 张震宇. 风力机翼型动态失速的 POD 模型降阶方法[J]. 南京航空航天大学学报, 2011, 43(5): 577-580

Zhang Zhenyu. Reduced-Order POD Model for Dynamic Stall of Wind Turbine Airfoils [J]. Nanjing University of Aeronautics and Astronautics, 2011, 43(5): 577-580 (in Chinese)

[5] 倪振华, 江棹荣, 谢壮宁. 本征正交分解技术及其在预测屋盖风压场中的应用 [J]. 振动工程学报, 2007, 20(1): 1-8

Ni Zhenhua, Jiang Zhaorong, Xie Zhuangning. POD Technique and Its Application to Prediction of Wind Pressure Fields on Roof [J]. Journal of Vibration Engineering, 2007, 20(1): 1-8 (in Chinese)

[6] Lumley J L. The Structure of Inhomogeneous Turbulent Flows.∥Atmospheric Turbulence and Radio Wave Propagation[M]. Yaglom AM, Tatarski VI. Nauka, Moscow, 1967: 166-178

[7] Ben Chiekh M, Michard M, Guellouz M S, et al. POD Analysis of Momentumless Trailing Edge Wake Using Synthetic Jet Actuation[J]. Experimental Thermal and Fluid Science, 2012, 46: 89-102

[8] 余化军. 圆柱和方柱绕流及矩形柱涡激振动的二维数值分析[D]. 天津大学, 2012

Yu Huajun. Two Dimensional Numerical Analysis of Flow over a Circular and Square Cylinder and Vortex-Induced Vibration of Rectangular Cylinder[D]. Tianjin: Tianjin University, 2012 (in Chinese)

[9] Singh A P, De A K, Carpenter V K, et al. Flow Past a Transversely Oscillating Square Cylinder in Free Stream at Low Reynolds Numbers[J]. International Journal for Numerical Methods in Fluids, 2009, 61(6): 658-682