力学大会2015集Arnoldi方法在边界层稳定性问题中应用

用(1.力学系，（2.市现代工程力学，：二本流的特征值问题对传统解法而言计算量太大，于是有人使用Arnoldi方法进行研究，但传统Arnoldi方法都使用有限元方法计算扰动方程，这就很难推广到大雷诺数和可压缩动中去。为此，本文使用Arnoldi法研究边界层稳定性问题计算扰动方程，避开有限元方法的局限性。其思想是，把线性扰动方程看成一个算子，通过Arnoldi这一算子的特征值问题过向一个较小的Krylov子空间投影，将大型方阵的特征值问题转化为一个小方阵的特征值问题问题大大简化。为了验证该方法的有效性，选取了三种情况。先在一本流下计算二维扰动线性稳定性理论（LST）结果对比；然后在一本流下计算三维扰动LST结果对比；最后在二本流下计算三维扰动渐加大雷诺数并与无粘Floquet理论的结果对比。结果显示该方法是有效的，且对大雷诺数的情况完全适用。将该方法应用于展向局部自由流湍流下边界层旁路转捩的二次失稳机理的研究，发现二次失稳的贡献远大于T-1引言 ）作者简介 ,讲师 transition）和旁路转捩转捩则会发生旁路转捩，一种典型的过程是[1]：外界低频扰动边界层，形成快（低频条纹）沿展向也是变化的，所以就涉及二本流（沿法向和展向变化）自然转捩的中的T-S波的参数是通过解较成解法，如Muller法[2]。但是这些方法在用于二本流时，计算量太随着计算数学的发展，Arnoldi迭始被用于研究流动稳定性问题外形，但是在大雷诺数的情况下迭代收敛很慢甚至发散。另一方面，在可压缩流动，尤其是超音速、高超音速流动中，往往存在激波，有限元方法也难以应用。所以传统Arnoldi方法多用于不可压的低雷诺数流动。克服传统Arnoldi方法的缺陷就需要是在解扰动方程时改用有限差分。虽然有限差分不太适合处理复杂几何形状的流场，但它对雷诺数大小没有要求，且是在机理性研究中，流场的几何形状大多比较规则，有限差分完全适用。另外，有限差分还适用于可压缩、超音速及高超音速流动的问题，而且能较好的处理流场中的激波。大多数工程问题，尤其是航空航题中，雷诺数都很大。因此，很有必要对传统Arnoldi方法做一下修改，将其中的有限元方法换为有限差分，如Merle等[17所作。本文正是做了这样的修改，但是作为起点本文还只是针对不可压流体的问题，并观察大雷诺数下的表现。在本文研究成功之后，则可以推广到可压缩、超音速及高超音速问题中去。此外，在实现Arnoldi迭代的方法上，传统的ld方法[15-验证。先在一本流下（只沿法向变化）计算二维扰动，并与线性稳定性理2u

uvw u

v

w

x

z2vu

vv

vp

1(2v

w

ww

u

v

w

Re(

p表示瞬时压力。该方程已经用来流速度U和当地边界层排移厚度δ无量纲化，雷诺数为ReU*δ*/υ*，其中υ*为来流的运动学粘性系数。将瞬时量表示为二本流与扰动之和，q(x,y,z,t)U(y,z)εq'(x,y,z,t)， 其中，瞬时量为q(uvw,p)T，基本流为U(U,V,WP)T，扰动为u'v'w'

Uv'Uw'Uu'Vu'Wu

p'

(2u'2u

Re

Vv'Vw'Uv'

p'

(2v'2v

V Re V

Wv'Ww'Uw'V V

w'p

1(2w'2w'

ˆ

ˆ 1(

αu

vˆ tyvˆzˆ yWzyRe(αvˆ

wˆ

ˆ

W

ˆ

方程（5）是关于qˆ的线性方程，对于某一波数α，它定义了一个线性演化算 qˆ(y,z,T) (T,α)qˆ(y,z,0)， (T,α)仍是一个线性算子，那么它就存在特征模态，即 述 其中，ω的实部ωr表示频率，虚部ωi表示增长率。对于某一特征模态，就可以描 综合（8）式和（10）式可见，特征λ和频率ω的关系 本文后面多用复数形式的频率表示特征值下面介绍求解方程（5）时的展向和法向离散方法，以及时间格式

，

本流U(y,z)看做一维部分Ua(y)和二维部分Ub(y,z)之和，U(y,z)=Ua(y)+Ub(其中，一维部分Ua为基本流U1Lz

Ua(y)=

U(z-Lz

Lz2π/β为展向基本波长。二维部分Ub也展为Fourier级

nU(n

差分格式，须将方程变换到等间距的计算坐标η下，变换形式为ηηy)。2

fj28fj18fj1fj2 fj216fj130fj16fj1fj j靠近边界处差分格式需降阶时间使用三阶差分格刚开始时降阶启动

t

11fn18fn19fn22 qˆ从t=0到t=T的演化，是通过上述数值方法求解方程（5）实现的。而（8）式中ˆ [q,q q q以得到一有k+1个向量的Krylov序 0，归一化得到矩 ˆ

ˆ

[q,q,q,...,

k 其中为序列中每个向量的范数 N(k ˆˆTTD(k k1 D(k (k1) 其 是Kronecker函。方程（20）是Tk和Tk1Krylov序列的基本关系式，下面介绍代数上如何处理这先将Tk和Tk1做QRQR k1k1 kk的上三角方阵Qk1和Rk1的情况类似。标记为k+1的矩阵与标记为k的矩阵在同样位置的元素完全相同，如Qk1的前k列与Qk相同，只不过后面多了一个列向量qk。定义一个(k1)k的矩则式（21）可写

H(k1) k1 k Q H(k k H(k

h* 是上Hessenberg矩阵，其最后一行只有一个非零元素，记为h

。用kk的方

H H QQHh*qeˆ kk kk

ˆT

基Qk作用时的投影。所得kk的小方阵Hk是否精确捕捉到了算子(T,α的作征向量，分别对应特征值λ0λ1,λ2,...,λk1。在（24）式中将Hk对角化，得QϒQϒϒ1Hϒh*qˆTk kk k k k ˆT k k Λk为kk的对角阵，其元素为Hk的特征值ΨkNk的矩阵 特征模 的残差可以用（26）式的最后一项表示ˆTυ kk 其中标量 kj是特征向量j的最后一个元素算子迭代就是将算子(T,α投影到一个k维的Krylov子空间，因为k并不大，所以在这个子空间里求kk的方阵Hk的特征值问题并 后由Hk的特征模态算出 3结果发点，仍有必要先在一本流下进行验证。本着由简单到复杂的顺序，本文先对一本流下的二维扰动进行计算，并与线性稳定性理论结果对比；然后对一计算三维扰动，逐渐加大雷诺数并与无粘的Floquet理论的结果[10]对图3和4分别给出了Re=1000 Neutralcurve(LST)Maxgrowth(Arnoldi)Maxgrowth(LST)

0

图1二维扰动的中性曲线和最大增长曲 图2二维扰动的最大增长

yy yy0

0

（a）中性曲线上 （b）中性曲线下6y4y20

（c）最大增长曲

yy yy0

0

（a）中性曲线上 （b）中性曲线下6y4y20

（c）最大增长曲再看一本流下的三维扰动的情况图5给出了由Arnoldi方法和线性稳定性理论得到的Re=1000处的三维扰动中性曲线和最大增长曲线，图6给出了此处的最大增长率。其中只画出了β0的部分，β0的部分与之对称。可见，Arnoldi方法和线性稳定性理论的结果完全

Neutralcurve(LST)Maxgrowth(Arnoldi)Maxgrowth(LST)

Maxgrowth(LST) 0

0 图5Re=1000处三维扰动的中性曲线和最大增长曲 图6Re=1000处三维扰动的最大增长

yy yy0

0 （b）中性曲线上支处的|

yy yy0

0 0 （d）中性曲线下支处的|

yy yy0

0

（f）最大增长曲线处的|

Neutralcurve(LST)Maxgrowth(Arnoldi)Maxgrowth(LST)

MaxgrowthMaxgrowth 0

图8Re=4000处三维扰动的中性曲线和最大增长曲 图9Re=4000处三维扰动的最大增长

yy yy0

0

（b）中性曲线上支处的|

yy yy0

0 （d）中性曲线下支处的|

yy yy0

0

（f）最大增长曲线处的|表1和2分别给出了图12给出了在有限元方遇到的迭代计算量增加、甚至迭代不收敛的问题。图11文献[10]（图6b）中的基本流：流向速度U的等(无粘率(无粘率6y4y20 z 频的，时间周期很长，因此可以在各个时刻对沿法向和展向变化的二本流做分析。Zhang等[1]研究了在不同来流湍流下边界层内的条纹，并在x旁路转捩。不同时刻，最不稳定的扰动都对应流向波数相差不多，都在α0.7附近。在ϕ=37π/128时刻，最大增长率最高；ϕ=15π/128时刻，最大增长率稍低；ϕ=81π/128和121π/128时刻，最大增长率比较低。图14（b）画出了相速度Cr。只有在时间相位为ϕ=15π/128时相速度较低，应该为内模态；其它时刻相速度为内模态[13]。图16显示，ϕ=37π/128时，特征函数集中分布在中心慢条纹的“两肩”，且实部和虚部对该条纹都是称分布的，考虑到其较高的相速（个慢条纹都是称分布的，考虑到其较高的相速度（见图14（b）），模态中的正弦模态[13]。结合图14（a）可见，中心处强条纹的不稳定模态的增长率要高于两侧弱条纹。yyzyy(a)ϕ

yyzyy(c)ϕϕ

图13x*=160mm 1 (a)ω

(b)图yzyyyzyy

图yzyyyzyy

yzyyyzyy

yzyyyzyy

4结论三种情况的计算，即一本流下的二维扰动，一本流下的三维扰动，以及二本流下的扰动，与线性稳定性理论和Floquet理论结果对比，验证了该方 作者在此感谢瑞典理工学院（KTH [1]Zhang,Y.,Zaki,T.,Sherwin.S.&Wu,X.Nonlinearresponseofalaminarboundarylayertoisotropicandspanwiselocalized-streamturbulence.AIAA2011-3292.6thAIAATheoreticalFluidMechanicsConference27-30JuneHawaii,US.(2011)[2],.流动稳定性.国防工业,Klebanoff,P.Effectof-streamturbulenceonalaminarboundarylayer.Bull.Am.Phys.Soc.16.(1971)Suder,K.L.,O’Brien,J.E.&Reshotko,E.Experimentalstudyofbypasstransitioninaboundarylayer.NASATech.Mem.100913.(1998)Kendall,J.M.Boundarylayerreceptivitytostreamturbulence.AIAAPaper.85-1695.(1990)Roach,P.E.&Brierley,D.H.Theinfluenceofturbulentstreamonzeropressuregradienttransitionalboundarylayerdevelopment.PartI:TestcasesT3AandT3B.InNumericalSimulationofUnsteadyFlowsandTransitiontoTurbulence(ed.O.Pironneau,W.Rodi,I.L.Ryhming,A.M.Saville&T.V.Truong),pp.319-347.CambridgeUniversityPress.(1992)Westin,K.J.A.,BoikoA.V.,Klingmann,B.G.B.,Kozlov,V.V.&P.H.Experimentsinaboundarylayersubjected streamturbulence.1.Boundarylayerstructureandreceptivity.J.FluidMech.281,193-218.Matsubara,M.&Alfredsson,P.H.Disturbancegrowthinboundarylayerssubjectedto -streamturbulence.J.FluidMech.430,149-168.(2001)Theofilis,V.GlobalLinearInstability.Annu.Rev.FluidMech.,43,319-352.Andersson,P.,Brandt,L.,Bottaro,A.&Henningson,D.S.Onthebreakdownofboundarylayerstreaks.J.FluidMech.428,29-60.(2001)Ricco,P.,Luo,J.&Wu,X.Evolutionandinstabilityofunsteadynonlinearstreaksgeneratedby -streamvorticaldisturbances.J.FluidMech.677,1-38.Liu,Y.,Zaki,T.A.&Durbin,P.A.Boundary-layertransitionbyin ctionofdiscreteandcontinuousmodes.J.FluidMech.604,199-233.(2008)Vaughan,N.J.&Zaki,T.A.Stabilityofzero-pressure-gradientboundarylayerdistortedbyunsteadyKlebanoffstreaks.J.FluidMech.681,116-153.(2011)计算物理.24,7-12.Barkeley,D.,Blackburn,H.M.&Sherwin,S.J.Directoptimalgrowthysisfortimesteppers.Int.J.Numer.Meth.Fluids.57,1435-1458.(2008)Blackburn,H.M.,Barkeley,D.&Sher