含压力基Navier-Stokes方程最优动力系统建模和分析*

王金城1, 齐 进2, 吴锤结1

(1. 大连理工大学 航空航天学院, 辽宁 大连 116024;2. 北京应用物理与计算数学研究所, 北京 100088)

(我刊编委吴锤结来稿)

摘要: 研究了采用压力基函数和速度基函数的Navier-Stokes方程的最优截断低维动力系统建模理论在黏性不可压缩流体中模拟了并排三方柱绕流流场,对此流场进行了含压力基函数和速度基函数的Navier-Stokes方程的最优动力系统建模,并以此为工具分析了三方柱绕流最优动力系统的动力学特性该研究得到了如下结论:三方柱绕流的最优动力系统的动力学行为为混沌,它与双方柱绕流场的极限环动力学特性有着本质的区别,因此可以通过多柱绕流增进尾流的复杂性,从而促进流体混合

关 键 词: 压力基函数; 最优动力系统; 多方柱绕流; 混沌

引 言

最优低维动力系统理论[1](简称“最优动力系统”)作为一种降维方法,通过最优化的方法求得最优基函数最优目标泛函的选取是建立最优动力系统的核心,可以选用余项[2-3]、残差[4]和加权残数[5-6]等指标作为目标泛函王金城等[7-8]对最优动力系统理论做了进一步发展,得到了适应复杂速度边界条件的Navier-Stokes方程最优动力系统理论和Navier-Stokes方程的脉动速度方程最优动力系统理论在现有的最优动力系统理论中,速度的不可压缩条件是通过优化速度基函数来实现的对于极其复杂的流动问题,速度不可压缩条件的处理很有讲究,需要通过速度基函数和压力基函数来共同实现速度的不可压本文通过引入压力基函数来对质量守恒方程进行最优低维建模,然后将其与包含压力基函数的动量守恒方程的最优动力系统模型联立,得到了基于压力基函数和速度基函数的Navier-Stokes方程的最优动力系统

本文的安排如下:第1节提出了含压力基的Navier-Stokes方程的最优动力系统理论,采用同时基于压力基函数和速度基函数的建模方法来保证速度的不可压缩条件第2节借助三方柱绕流问题,完成了含压力基的Navier-Stokes方程的最优动力系统模型的构建第3节借助构建好的最优动力系统模型,对并排三方柱绕流最优动力系统的复杂动力学特性进行了分析最后总结并得出了结论

1 含压力基的不可压缩Navier-Stokes方程的最优动力系统建模

1.1 含压力基的最优动力系统建模理论

本文使用张量指标记号[9],指标ij为物理空间中的坐标(i, j = 1,2,3)

黏性不可压缩三维流动问题的无量纲Navier-Stokes方程、连续性方程、初始条件和边界条件为

(1)

其中,u0ip0分别为初始速度场和初始压力场,FΩ表示Dirichlet边界条件、Neumann边界条件和Robin边界条件中的任意一种

定义速度基函数ξki的函数空间为BN,压力基函数χm的函数空间为ZNp,两个函数空间的具体定义为

(2)

其中,ξkiχm二阶可微,δkl和δmq为Kronecker delta函数,HNHNp分别为N维和Np维Hilbert空间,N是速度基函数的截断阶数,Np是压力基函数的截断阶数

速度ui在函数空间BN中的低维近似表达式,以及压力p在函数空间ZNp中的低维近似表达式为

(3)

其中,uRipR为余项,akcm都是时间的函数,k=1,2,…,Nm=1,2,…,Np

将式(3)代入方程(1)中,得

(4)

其中 l,h=1,2,…,N

首先对方程组(4)第一式施加边界条件:

(5)

其中(akξki)b的具体形式由具体的边界条件决定[8]

然后再对满足边界条件的方程组(4)第一式构造Galerkin投影方程,即将式中的各项在任意速度基函数ξki上投影,来得到满足边界条件的Galerkin投影方程:

(6)

其中

(7)

再对方程组(4)中第二式构造Galerkin投影方程,即将式中的各项在任意压力基函数χm上投影,得

(8)

其中

(9)

综上所述,得到如下Galerkin投影方程组:

(10)

至此,偏微分方程定解问题(4)转化为如下的动力系统定解问题:

(11)

注意,由于方程组(11)中第一式的压力基的系数没有对时间求导,且它还是未知的,因此,此常微分方程组无法用标准的常微分方程求解方法进行求解,具体求解方法见1.3小节

需要着重说明的是,上述定解问题中,边界条件通过积分项中边界处的Ui和边界附近的导数Ui, j,Ui, jj与内部流场耦合起来动力系统方程求解完毕后,内部速度场由基函数及其系数给出,而边界处的速度由边界条件给出

1.2 含压力基最优动力系统模型

因为含压力基的动力系统中包含了连续性方程的投影方程,因此无需再将速度的不可压缩性指标作为目标泛函

下面将速度基函数的正交性作为目标泛函,压力基函数的正交性也作为目标泛函,建立关于速度基函数正交性和压力基函数正交性的广义目标泛函,得

(12)

其中

(13)

通过对方程(12)的变分运算,得到含压力基的Navier-Stokes方程的最优动力系统模型方程

1) 关于ak,t的含压力基的动力系统方程

(14)

2) 广义目标泛函梯度方程

Ω内,

(15)

至此,即可将1.3小节的常微分方程组求解方法和共轭梯度优化算法[10]结合在一起,对最优动力系统模型进行求解(亦可采用其他任意一种优化算法来代替共轭梯度优化算法)

1.3 压力基-速度基耦合动力系统方程组求解方法

采用Runge-Kutta法求解常微分方程组(14)

首先将方程组(14)第二式对时间求导,得

(16)

再将ak,t改写为al,tδlk,可得上述方程的矩阵形式:

(17)

根据Runge-Kutta法求解公式[11],有

(18)

式中由下面的线性代数方程组确定:

(19)

其中

(20)

若已知tn时刻的速度系数则按照式(19)和(20),即可求出再代入式(18),即得到tn+1时刻的速度系数压力系数可由下列方程求出:

(21)

若速度基函数的阶次和压力基函数的阶次相同,虽然速度的解较精确,但压力的解误差很大,如果压力基函数比速度基函数的阶次低,则可以得到较为准确的压力场

1.4 边界条件

基于王金城等[8]提出的任意边界条件的实现方法,并排三方柱绕流问题的边界条件FΩ如下(其具体的数值实现方法请参考文献[8],压力边界条件的处理方法与速度边界的处理方法相同,本文不再赘述):

1)左边界 速度为入流边界和无滑移边界的混合,压力零梯度;

2)右边界 速度为连续出流边界,压力给定;

3)前后上下边界 速度为滑移边界,压力零梯度

2 含压力基的最优动力系统的POD分析

2.1 并排三方柱绕流算例和POD分析

本节基于Reynolds数为200的并排三方柱绕流问题,进行含压力基的Navier-Stokes方程的最优动力系统的流场分析需要说明的是,除了特殊说明外,本计算所涉及的参数都是无量纲的三维流场的尺寸为Ω=25.0×4.0×10.0边界条件见1.4小节

需要具体明确的边界条件是,左边界由三个方柱(无滑移)和来流共同组成,给定入流速度uB=(1.0, 0.0, 0.0),方柱的特征尺寸为1.0,三个方柱是等间距的,且间距为1.0采用NASA-VOF3D程序[12]计算并排三方柱绕流算例,其粗糙网格数为50×8×20,网格尺寸为Δxyz=0.5,时间步长取为Δt=0.01图1是已经充分发展的三方柱绕流,可见流场中的周期性结构很少,且三个方柱后的尾流的摆动方向不同

(a) 速度模值云图 (b) 带状流线
(a) The contour of velocity magnitudes(b) Ribboned streamlines
图1 t=750.0时粗糙流场
Fig. 1 The coarse velocity field at t=750.0

为了解释图中的颜色,读者可以参考本文的电子网页版本,后同

表1 速度POD基能量分布

Table 1 The energy distribution of velocity POD bases

orderenergy of kth base ek/%energy of first k bases ek/%188.7288.7223.1191.8332.9494.7740.9895.7550.8996.6460.7197.3570.6297.9680.3998.3590.3798.72100.2398.95

图2 前10阶速度POD基的能量占比
Fig. 2 The energy ratios of the first 10 velocity POD bases

为了使构建出的动力系统能刻画流动的长期行为,需要使用充分发展的流动数据库来提取POD基,因此选取t=750.0~790.0的流场和压力场,并以Δt′=0.2作为提取间隔,分别获得速度POD基和压力POD基[13]

从表1和图2中的速度POD基能量分布可以看出,三方柱绕流问题中的小尺度结构极其多选取前10阶速度POD基函数来拟合速度场,以更好地捕捉三方柱绕流问题的小尺度流动结构

表2 压力POD基能量分布

Table 2 The energy distribution of pressure POD bases

orderenergy of kth base ek/%energy of first k bases ek/%199.997 15099.997 15020.000 575 699.997 72630.000 535 099.998 26140.000 491 699.998 75250.000 226 299.998 97860.000 170 199.999 14970.000 135 399.999 28480.000 132 499.999 41690.000 117 599.999 534100.000 109 399.999 643

与速度POD基不一样,前10阶压力POD基函数捕捉了压力场近乎全部的能量(表2),说明压力场没有速度场复杂我们选取前8阶压力POD基函数来拟合原始压力场

2.2 近似全局最优速度基与最优压力基

根据文献[5]中的双尺度全局最优化方法可知,初始迭代基的网格尺寸应为粗网格的一半,即Δx′=Δy′=Δz′=0.25将前10阶速度基和前8阶压力基插值到密网格上,然后代入共轭梯度优化算法进行优化,即可得到近似全局最优速度基和最优压力基图3展示了式(12)中的广义目标泛函值Jg的收敛曲线(反映了全局最优的速度基函数正交性和压力基函数正交性的变化)

图3 Jg的收敛曲线
Fig. 3 The convergence curve of Jg

第一阶最优速度基函数可以看作是时间平均速度场,而高阶最优速度基函数则用于刻画小尺度结构,三方柱绕流问题的高阶最优速度基函数在空间上几乎没有了周期性或规律性结构(图4),最优压力基函数在复杂性方面与最优速度基函数类似(图5)

(a) 第一阶 (b) 第二阶
(a) The 1st order (b) The 2nd order

(c) 第三阶 (d) 第四阶
(c) The 3rd order (d) The 4th order

(e) 第五阶 (f) 第六阶
(e) The 5th order (f) The 6th order

(g) 第七阶 (h) 第八阶
(g) The 7th order (h) The 8th order

(i) 第九阶 (j) 第十阶
(i) The 9th order (j) The 10th order
图4 流向中截面上的最优速度基模值云图
Fig. 4 The magnitude contours of optimal velocity bases on the streamwise middle section

(a) 第一阶 (b) 第二阶
(a) The 1st order (b) The 2nd order

(c) 第三阶 (d) 第四阶
(c) The 3rd order (d) The 4th order

(e) 第五阶 (f) 第六阶
(e) The 5th order (f) The 6th order

(g) 第七阶 (h) 第八阶
(g) The 7th order (h) The 8th order
图5 流向中截面上的最优压力基云图
Fig. 5 The contours of optimal pressure bases on the streamwise middle section

2.3 最优基(速度基与压力基)对流场和压力场的拟合

为了检验最优速度基在刻画流场结构方面的准确性,我们用前10阶最优速度基及其系数拟合出近似速度,并将其同数据库速度进行对比(图6),发现近似速度场非常接近数据库速度场

(a) t=750.0最优动力系统流场 (b) t=750.0数据库流场
(a) The optimal dynamical system for t=750.0 (b) The database for t=750.0
图6 含压力基最优动力系统速度模值和数据库速度模值比较
Fig. 6 The comparison of the velocity magnitudes between the optimal dynamical system and the database

(a) t=750.0最优动力系统压力场(b) t=750.0数据库压力场
(a) The optimal dynamical system for t=750.0 (b) The database for t=750.0
图7 含压力基最优动力系统压力场和数据库压力场比较
Fig. 7 The comparison of the pressure magnitudes between the optimal dynamical system and the database

同时我们用前8阶最优压力基拟合的近似压力场与数据库压力场进行对比(图7),发现近似压力场的误差也较小说明10阶最优速度基函数和8阶最优压力基函数可以较为完整地抓住流动的主要结构

3 含压力基最优动力系统的动力学特性

3.1 含压力基最优动力系统的相空间轨道

为观察含压力基最优动力系统在相空间中的轨迹,选取t=20 250.0~20 750.0无量纲时段内的最优速度基系数的轨迹进行观测首先依次观测三个相邻整数阶的最优动力系统的三维轨迹,可以发现,三维轨迹是杂乱无章的,如图8所示

(a) a1-a2-a3 (b) a4-a5-a6

(c) a6-a7-a8 (d) a8-a9-a10
图8 含压力基最优动力系统的三维相空间轨道
Fig. 8 Three-dimensional phase portaits of the optimal dynamical systems with pressure bases

(a) a1-a2(b) a3-a4

(c) a5-a6(d) a7-a8(e) a9-a10
图9 含压力基最优动力系统的二维相空间轨道
Fig. 9 Two-dimensional phase portaits of the optimal dynamical systems with pressure bases

其次,依次观测两个相邻整数阶的含压力基的最优动力系统的二维轨迹,如图9所示通过二维相轨道也可以发现,含压力基的最优动力系统的动力学行为极其复杂

(a) 相平面a1=0 (b) 相平面a2=0
(a) The phase plane of a1=0 (b) The phase plane of a2=0

(c) 相平面a3=0 (d) 相平面a4=0
(c) The phase plane of a3=0 (d) The phase plane of a4=0

(e) 相平面a5=0 (f) 相平面a6=0
(e) The phase plane of a5=0 (f) The phase plane of a6=0
图10 含压力基最优动力系统的六个典型Poincaré截面
Fig. 10 Six typical Poincaré sections of the optimal dynamical systems with pressure bases

3.2 含压力基最优动力系统的Poincaré截面

Poincaré截面可以被用于判断系统是否具有混沌特性从含压力基的最优动力系统在六个典型相平面上的Poincaré截面(图10)可见,各个Poincaré截面都是由一系列不重叠的点组成的,因此可以断定含压力基的最优动力系统在Re=200时的长期动力学特性是混沌

3.3 含压力基最优动力系统的分岔特性

下面研究最优动力系统随Reynolds数变化的分岔特性我们以1作为Reynolds数的步长,研究Reynolds数从200增加到700的过程中含压力基的最优动力系统的分岔特性,结果如图11所示其中,横轴代表Reynolds数,图11(a)纵轴代表a9,图11(b)纵轴代表a10,图中的点代表在t=6 750.0~12 750.0无量纲时段内,含压力基的最优动力系统轨迹穿过a1=15相平面时,在该Poincaré截面上的位置

(a) a9的分岔特性(b) a10的分岔特性
(a) The bifurcation characteristics of a9(b) The bifurcation characteristics of a10
图11 Re=200~700的含压力基最优动力系统分岔图
Fig. 11 The bifurcation diagram of the optimal dynamical systems for Re=200~700

值得注意的是,分岔图只显示了Re=200~500部分的图像,这是因为Re=500~700的含压力基最优动力系统是不稳定的(即发散的),不会在分岔图中留下点,因此在画图时只绘制了Re=200~500部分的图像Re=200~500区间内,随着Re的增长,含压力基的最优动力系统的动力学行为越来越不稳定,并且出现间歇性的发散

3.4 含压力基最优动力系统的功率谱分析

采用文献[8]中的功率谱分析方法来分析含压力基的最优动力系统的功率谱取无量纲时间t=750.0~20 750.0的含压力基最优动力系统时间序列(共十万个点),绘制出如图12所示的功率谱,其中的横纵坐标都是对数坐标

(a) a1的功率谱 (b) a2的功率谱
(a) The power spectrum of a1 (b) The power spectrum of a2

(c) a3的功率谱 (d) a4的功率谱
(c) The power spectrum of a3 (d) The power spectrum of a4

(e) a5的功率谱 (f) a6的功率谱
(e) The power spectrum of a5 (f) The power spectrum of a6

(g) a7的功率谱 (h) a8的功率谱
(g) The power spectrum of a7 (h) The power spectrum of a8

(i) a9的功率谱 (j) a10的功率谱
(i) The power spectrum of a9 (j) The power spectrum of a10
图12 含压力基最优动力系统的功率谱
Fig. 12 The power spectra of the optimal dynamical systems with pressure bases

图13 含压力基最优动力系统Lyapunov指数集的演化
Fig. 13 The evolution of Lyapunov exponent spectra of the optimal dynamical systems

理论上,功率谱中的尖峰代表周期运动,而宽峰和噪声背景代表混沌运动由图12可知,各阶最优动力系统的功率都具有噪声背景和宽峰,因此该最优动力系统的动力学特性为混沌(与3.2小节结论一致)

3.5 含压力基最优动力系统的Lyapunov指数

Lyapunov指数集能够更为准确地判断动力系统的长期动力学行为

采用文献[8]中的方法来分析含压力基的最优动力系统在t=750.0~20 750.0无量纲时段内的Lyapunov指数集,结果如图13所示(绘图时将时间坐标的起点平移至t=0.0)10个Lyapunov指数中6个为正,1个为零,3个为负,反映出此含压力基的最优动力系统为混沌系统,与Poincaré截面、功率谱的相关结果的混沌特性相符合

4 结 论

本文介绍了含压力基函数和速度基函数的Navier-Stokes方程的最优截断低维动力系统建模理论数值模拟了三方柱绕流问题,得到初始迭代压力基和速度基,构建了该问题的含压力基的Navier-Stokes方程最优动力系统,通过对其相空间轨道、分岔特性、功率谱、Lyapunov 指数集等动力学特性的分析,发现此最优动力系统是混沌的,它与文献[7]的双方柱绕流场的极限环动力学特性有着本质的区别,即三方柱绕流场的动力学特性更加复杂由此可见,可以通过多柱绕流增进尾流的复杂性,从而促进流体混合

参考文献:

[1] WU C J. Optimal truncated low-dimensional dynamical systems[J]. Discrete and Continuous Dynamical Systems, 1996, 2(4): 559-583.

[2] WU C J, SHI H S. An optimal theory for an expansion of flow quantities to capture the flow structures[J]. Fluid Dynamics Research, 1995, 17(2): 67-85.

[3] 吴锤结, 赵红亮. 不依赖数据库的最优动力系统建模理论及其应用[J]. 力学学报, 2001, 33(3): 289-300.(WU Chuijie, ZHAO Hongliang. A new database-free method of constructing optimal low-dimensional dynamical systems and its application[J]. Acta Mechanica Sinica, 2001, 33(3): 289-300.(in Chinese))

[4] WU C J, WANG L. A method of constructing a database-free optimal dynamical system and a global optimal dynamical system[J]. Science in China(Series G): Physics, Mechanics & Astronomy, 2008, 51(7): 905-915.

[5] PENG N F, GUAN H, WU C J. Research on the optimal dynamical systems of three-dimensional Navier-Stokes equations based on weighted residual[J]. Science China: Physics, Mechanics & Astronomy, 2016, 59(4): 644-701.

[6] PENG N F, GUAN H, WU C J. Optimal dynamical systems of Navier-Stokes equations based on generalized helical-wave bases and the fundamental elements of turbulence[J]. Science China: Physics, Mechanics & Astronomy, 2016, 59(11): 114713.

[7] 王金城, 齐进, 吴锤结. Navier-Stokes方程的脉动速度方程的最优动力系统建模和动力学分析[J]. 应用数学和力学, 2020, 41(3): 235-249.(WANG Jincheng, QI Jin, WU Chuijie. Dynamics analysis and modelling of optimal dynamical systems of the fluctuation velocity equations of incompressible Navier-Stokes equations[J]. Applied Mathematics and Mechanics, 2020, 41(3): 235-249.(in Chinese))

[8] 王金城, 齐进, 吴锤结. 不可压缩Navier-Stokes方程最优动力系统建模和分析[J]. 应用数学和力学, 2020, 41(1): 1-15.(WANG Jincheng, QI Jin, WU Chuijie. Analysis and modelling of optimal dynamical systems of incompressible Navier-Stokes equations[J]. Applied Mathematics and Mechanics, 2020, 41(1): 1-15.(in Chinese))

[9] 黄克智, 薛明德, 陆明万. 张量分析[M]. 北京: 清华大学出版社, 2003.(HUANG Kezhi, XUE Mingde, LU Mingwan. Tensor Analysis[M]. Beijing: Tsinghua University Press, 2003.(in Chinese))

[10] 叶庆凯, 王肇明. 优化与最优控制中的计算方法[M]. 北京: 科学出版社, 1986.(YE Qingkai, WANG Zhaoming. Computational Methods of Optimization and Optimum Control[M]. Beijing: Science Press, 1986.(in Chinese))

[11] 章本照. 流体力学中的有限元方法[M]. 北京: 机械工业出版社, 1986.(ZHANG Benzhao. Finite Element Method in Fluid Mechanics[M]. Beijing: China Machine Press, 1986.(in Chinese))

[12] TORREY M D, MJOLSNESS R C, STEIN L R. NASA-VOF3D: a three-dimensional computer program for incompressible flows with free surfaces[R]. NASA STI/Recon Technical Report N, 1987.

[13] SIROVICH L. Turbulence and the dynamics of coherent structures, part Ⅲ: dynamics and scaling[J]. Quarterly of Applied Mathematics, 1987, 45(3): 583-590.

(Contributed by WU Chuijie, M. AMM Editorial Board)

Modelling and Analysis of Optimal Dynamical Systems of Incompressible Navier-Stokes Equations With Pressure Base Functions

WANG Jincheng1, QI Jin2, WU Chuijie1

(1. School of Aeronautics and Astronautics, Dalian University of Technology,Dalian, Liaoning 116024, P.R.China;2. Institute of Applied Physics and Computational Mathematics,Beijing 100088, P.R.China)

Abstract: The modelling theory for optimal truncated low-dimensional dynamical systems of Navier-Stokes equations with pressure base functions and velocity base functions was studied. In the viscous incompressible fluid, the flow field around 3 square columns was simulated. According to that flow problem, the optimal dynamical systems of the Navier-Stokes equations with pressure base functions and velocity base functions were modelled and studied. The results show that, the dynamics behavior of the optimal dynamical systems around the 3 square columns is chaos, which is essentially different from the limit cycle dynamics behavior of the flow field around 2 square columns. As a result, the complexity of the wake increases in the multi-column flow, which thereby shows that the increase of number of columns can promote.

Key words: pressure base function; optimal dynamical system; flow around multiple square columns; chaos

*收稿日期: 2019-09-06;

修订日期:2020-03-04

基金项目: 国家自然科学基金(11601033;11372068);国家重点基础研究发展计划(973 计划)(2014CB744104)

作者简介:

王金城(1994—),男,硕士生(E-mail: jcwangdut@foxmail.com);

齐进(1977—),女,副研究员(E-mail: qi_jin@iapcm.ac.cn);

吴锤结(1955—),男,教授(通讯作者. E-mail: cjwudut@dlut.edu.cn).

引用格式: 王金城, 齐进, 吴锤结. 含压力基Navier-Stokes方程最优动力系统建模和分析[J]. 应用数学和力学, 2020, 41(8): 817-833.

中图分类号: O357.41

文献标志码:A

DOI: 10.21656/1000-0887.400276

Foundation item: The National Natural Science Foundation of China(11601033; 11372068); The National Basic Research Program of China (973 Program)(2014CB744104)