台风强度变化最快时涡度的解析解

范景威 周伟灿 冯也骋 官元红

0 引言

台风强度变化是台风研究的一个重要课题,是当前台风观测和预报中的重难点与前沿问题．影响台风强度变化的因子可以归结为以下三类:环境气流、下垫面和台风内部结构[1-2]．随着监测手段的发展以及数值模式理论和技术的不断成熟,国内外学者针对台风强度变化的研究取得了一定成果,已有研究内容涉及诸多方面,包括环境风垂直切变、冷空气、高低空急流、西风槽、地形以及海面温度和喷沫等[3-9]．

变分法是研究泛函极值的数学方法,其主要思想是构造适当的泛函,对泛函作变分并取变分为零,研究此时满足的函数性态．变分法广泛地应用于力学、数学物理反问题、气象资料同化和大气运动的稳定性问题中．Arnold[10]和Vallis[11]运用变分原理,指出流体运动和地转平衡分别是在位涡守恒及势能守恒的约束下,使总能量达到极小．伍荣生[12]利用Finlayson[13]提出的限制性变分原理,求得了与非线性涡度方程对应的Lagrange函数的近似表达式．Barth等[14]从线性浅水方程和系统能量变分出发,推导了从高程和深度平均速度中去除惯性重力波的方法,指出变分滤波器在提前降低重力波方面的优势．黄思训等[15]通过二次变分方法研究台风流场结构,将实测风场分解成无旋场和有旋场,再将有旋流场分解成对称涡旋和非对称涡旋对,两次提取到的都是最大涡旋．Wang等[16]用变分法推导了赤道电离层Rayleigh-Taylor(RT)稳定和不稳定的充分条件,并通过系统特征值计算了RT稳定和不稳定区域．Badin等[17]将变分原理应用于波动动力学中,得到地表水波方程以及非线性问题的波频散关系．近年来,随着变分同化技术的发展与成熟,它被广泛地应用于大气和海洋等不同领域[18-21],其目的是充分利用气象四维观测资料,在动力预报模式中通过变分得到大气海洋最优状态估计,以获得更好的模式初始场及预报效果．

目前对台风强度的研究多采用统计、诊断分析、数值模拟和敏感性试验等方法,缺少理论层面的研究．此外,前人的研究多采用直接分解的方法将台风流场分解成涡旋流场和无旋流场[15],这种分解无法得到有旋场中的最大涡旋．本文通过能量泛函变分,研究台风强度变化率最大时台风能量满足的关系,在四力平衡关系的基础上提取到有旋流场中的最大涡旋,给出台风强度变化最快时涡度的解析解．该结果为预报方程提供了一定的理论指引,对研究台风发展过程中速度的变化趋势和台风的层次结构具有实用价值．

1 能量泛函变分

假设大气是均匀不可压缩的流体,由连续方程有

(1)

其中,u=u(x,y,z,t)为纬向风速,v=v(x,y,z,t)为经向风速,w=w(x,y,z,t)为垂直风速．

大气运动遵循牛顿第二运动定律．对于惯性坐标系,牛顿第二运动定律可表示为

(2)

(3)

其中,V3=(u,v,w)为三维速度矢量,Ω为地球自转角速度矢量,R为气块所在的纬圈平面内从地轴到该气块的距离矢量,其大小为R=rcosφ．作用于空气微团的气压梯度力、分子黏性力、湍流黏性力和重力的表达式分别为

(4)

(5)

(6)

(7)

其中,ρ为密度,p=p(x,y,z,t)为气压,3为三维微分算子,μ为分子黏性系数,为三维拉普拉斯算子,Az为湍流交换系数,G为引力常数,M为地球的质量,r为空气微团的位置矢量．

将式(3)—(7)代入式(2),经整理得到旋转坐标系中矢量形式的动量方程:

(8)

其中,F=F2+F3为分子黏性力与湍流黏性力之和,即摩擦力．用V3点乘式(8)两端得到旋转坐标系中的动能方程：

(9)

定义能量泛函I(u,v,w,p)∈C2(Ω),

(10)

其中,Ω是任意台风区域,用该泛函表示台风强度变化率．假设在∂Ω上u,v,w,p以及u,v,w关于x,y,z的所有一阶偏导数均为常数,此时台风在区域边界上可视为基本气流[22]．对I作变分,有

(11)

令δI=0,利用δu,δv,δw的任意性得到Euler-Lagrange方程为

(12)

(13)

(14)

矢量形式为

2F+fp-g-3K=0,

(15)

其中,fp为气压梯度力,g=(0,0,g)．式(15)表明,当台风强度变化率达到最大时,摩擦力、气压梯度力、重力和动能梯度满足四力平衡．因此,通过这四个力确定的向量可以作为台风强度的预报因子,更准确地确定台风系统强度变化率达到最大的时间点,研究台风强度的变化规律．此外,将式(15)变形成

(16)

则台风强度变化最快时的摩擦力表达式由式(16)给出,它为计算较困难且常常被忽略的摩擦力提供了一种新的计算方法．

2 台风风场分解

2.1 变分提取有旋气流

为了进一步研究台风强度变化最快时涡度的形式和流场的变化,采用黄思训等[15]的方法对台风风场进行分解．首先将风场u(x,y,z),v(x,y,z)分解成平面场u′(x,y),v′(x,y)和垂直场H(z),即

u(x,y,z)=u′(x,y)H(z),

(17)

v(x,y,z)=v′(x,y)H(z),

(18)

其中,风场u,v对应台风强度变化率最大的时刻,分别满足方程(12)和(13)．由于台风的涡旋运动受环境气流引导,因此台风流场可以分解成有旋流场和无旋流场．利用变分方法从风场u′,v′中提取无旋气流u1(x,y),v1(x,y)[15],使得泛函

(v′-v1)2]dΩ=min!

(19)

且u′,v′,u1,v1满足

(20)

(21)

(22)

这是一个条件变分问题,引入Lagrange乘子λ1(x,y,z),λ2(x,y,z),λ3(x,y,z),有

(23)

于是

(24)

假设在∂Ω上λ1,λ2,λ3以及λ2,λ3关于x,y,z的所有一阶偏导数均为0,则有

(25)

利用δu,δv,δw,δu1,δv的任意性,可得它们对应的Euler-Lagrange方程分别为

(26)

(27)

(28)

(29)

(30)

(31)

(32)

(33)

(34)

(35)

根据式(29)和(30),有旋流场(u2,v2)和无旋流场(u1,v1)的表达形式分别为

(36)

(37)

(38)

(39)

在有旋场中引入流函数ψ(x,y),满足

(40)

涡度为

(41)

将式(30)对x作微商,式(29)对y作微商,相减并结合式(33)得

(42)

将式(35)对x作微商,式(34)对y作微商,相减并结合式(42)得

(43)

将式(33)、(41)和(42)联立得

ψ=λ1.

(44)

因此,在平面风场中,二维Lagrange乘子λ1(x,y)即为有旋流场的流函数ψ(x,y)．故方程(43)可改写为

(45)

由于ζ=ζ(x,y),H=H(z),对方程(45)采用分离变量法,并设:

(46)

则:

(47)

其中,λ为常数．

由式(31)得

(48)

将式(46)、(48)代入式(26)、(27)得

(49)

(50)

给定方程(46)的边界条件为H(0)=H,H(z1)=H1,方程(46)的解为

(51)

(52)

2.2 涡度的解析解

根据黄思训等[15],台风流场一般在圆域中经变分分解能提取到最大涡旋,故将方程(47)转化成极坐标(r,θ)下的方程

0

(53)

为满足齐次边界条件,作变换

(54)

则η(0,θ)=η(r0,θ)=0．记ζ(r0,θ)=ζ0,方程(53)转化成如下形式:

(55)

令η(r,θ)=R(r)Φ(θ),ζ(r,θ)=ρ(r)Φ(θ),则有

(56)

(57)

对方程(57)作变量分离得

(58)

其中,μ为常数．结合η(r,θ)关于θ的周期性,可得

(59)

(60)

式(59)是带有周期性条件的特征值问题,其特征值和特征函数分别为

μn=n2,Φn(θ)={cosnθ,sinnθ}.

(61)

其中,n≥0．它的解为

(62)

其中,an,bn为常数．

(63)

当λ=0时,方程(63)为非齐次欧拉方程,它在边界条件下的解为

(64)

此时,结合式(56)得涡度ζ的表达式为

(65)

其中,A′n,B′n,c0为常数．

当λ>0时,采用常数变易法[24]得方程(63)的通解为

(66)

(67)

由边界条件得

(68)

(69)

因此,台风强度变化率最大时有旋流场涡度的解析解为

(70)

涡度拟能的表达式为

(71)

2.3 有旋流场的解析解

当λ=0时,流函数ψ(r,θ)满足:

(72)

根据ζ(r,θ)的表达式,ψ(r,θ)具有如下形式:

E(r)sin2nθ+F(r)cos2nθ+G(r)sinnθcosnθ.

(73)

由λ1|∂Ω=0可知ψ|∂Ω=0,故方程(72)的边界条件为

C(r0)=D(r0)=E(r0)=F(r0)=G(r0)=0.

(74)

将式(73)代入方程(72),整理得

(75)

(76)

(77)

(78)

(79)

方程(75)、(76)和(79)满足边界条件的解分别为

(80)

(81)

G(r)=0,

(82)

(83)

(84)

且根据式(74),边界条件为H(r0)=J(r0)=0．解得

(85)

J(r)=0.

(86)

由式(85)和(86)得

(87)

因此,流函数ψ(r,θ)的表达式为

(88)

根据式(40),有旋场u2(r,θ),v2(r,θ)满足:

(89)

(90)

因此,有旋流场(u2,v2)的表达式为

(91)

(92)

由变分的性质可知,分解得到的有旋气流(u2,v2)是最大的涡旋．

当λ>0时,为方便求解流函数,对涡度ζ作如下简化:

(93)

(94)

同理,ψ(r,θ)具有如下形式:

(95)

边界条件为

H(r0)=M(r0)=0.

(96)

将式(95)代入方程(94),整理得

(97)

(98)

H(r)=anGnr3-an(Gnr03-n+Fnr02n+2)rn+anFnr3n+2,

(99)

M(r)=bnGnr3-bn(Gnr03-n+Fnr02n+2)rn+bnFnr3n+2.

(100)

因此,流函数ψ(r,θ)的表达式为

(101)

根据式(89)、(90),有旋流场(u2,v2)的表达式为

(3n+2)Fnr3n+1]sinθsinnθ+

an[-3Gnr2+n(Gnr03-n+Fnr02n+2)rn-1-

(3n+2)Fnr3n+1]sinθcosnθ+

an[nGnr2-n(Gnr03-n+Fnr02n+2)rn-1+

nFnr3n+1]cosθsinnθ+

bn[-nGnr2+n(Gnr03-n+Fnr02n+2)rn-1-

nFnr3n+1]cosθcosnθ,

(102)

nFnr3n+1]sinθsinnθ+bn[-nGnr2+n(Gnr03-n+

Fnr02n+2)rn-1-nFnr3n+1]sinθcosnθ+

bn[3Gnr2-n(Gnr03-n+Fnr02n+2)rn-1+

(3n+2)Fnr3n+1]cosθsinnθ+an[3Gnr2-

n(Gnr03-n+Fnr02n+2)rn-1+

(3n+2)Fnr3n+1]cosθcosnθ.

(103)

2.4 数值试验

当λ=0时,根据有旋流场的表达式(91)、(92),取n=1,r0=1,A1=1,B1=1,c0=16,ρ(1)=2,台风强度变化最快时的流线图如图1所示．可以看出,流线图在经向和纬向上呈均匀梯度分布,符合实际流场分布[15]．因此,在有旋流场的表达式(91)、(92)中,取n=1可得有旋流场及其涡度的解析解．

图1 λ=0,n=1时台风强度变化率最大时刻的流线图Fig.1 Plot of streamline at the maximum change rate of TC intensity when λ=0,n=1

图2 λ>0,n=0时台风强度变化率最大时刻的流线图Fig.2 Plot of streamline at the maximum change rate of TC intensity when λ>0,n=0

因此,台风强度变化最快时涡度、涡度拟能和流场的一个解析解分别为

(104)

(105)

(106)

(107)

有旋流场(u2,v2)更直观清晰地反映台风每一层的流场结构,涡度和涡度拟能的变化有助于研究台风的运动机理和风场中能量的变化．台风的整体结构取决于垂直廓线H(z),而λ是决定H(z)的主要参数,因此λ决定整个台风的层次结构．

3 结论

本文从理论层面研究台风强度变化率最大时台风能量满足的关系,给出台风强度变化最快时涡度的解析解．结果表明,当台风强度变化率达到最大时,摩擦力、气压梯度力、重力和动能梯度满足四力平衡．通过这四个力确定的向量可以作为台风强度的预报因子,更准确地确定台风系统强度变化率达到最大的时间点,研究台风强度的变化规律．在四力平衡关系的基础上,进一步通过风场变分分解提取到有旋场中的最大涡旋,得到台风强度变化最快时涡度和流场的一个解析解,从理论上证明了台风发展最快时前人的研究[15]结果,对研究台风发展过程尤其是平衡过程中的运动机理和台风的空间结构具有一定的指导意义,为台风路径和强度预报提供了一定的理论指引．此外,本文的结论为计算较困难且常常被忽略的摩擦力提供了一种新的计算方法．