0 前言
目前,对三种典型的射流——平面液膜、圆射流和环状液膜气液相界面的数理模化和碎裂机理的研究已经有所进展,但尚未完善[1]。液体表面波线性不稳定性理论（或者称为线性稳定性理论）是以气、液体的质量、动量守恒为基础,以连续性方程和纳维-斯托克斯（NAVIER-STOKES）方程组作为控制方程组,代入运动学边界条件、附加边界条件和动力学边界条件,考虑到气液体速度、密度、液体的表面张力和黏性及气体可压缩性等影响,推导得到色散关系式（dispersion relation）。它是一个复数指数方程。表面波增长率随表面波的波数或波长的变化关系是隐含给出的。波数k与波长λ的关系为$k=2\text{ }\pi\text{ }/\lambda $（m-1）。由于色散关系式很复杂,无法得到其解析解。故应用穆勒（MULLER）方法[2]编制FORTRAN语言程序,可以求得色散关系式的数值解,得到表面波增长率随表面波数的变化曲线或者曲面。线性不稳定性理论模型将不考虑雷诺应力的影响,这样方程组就是封闭的,不需要补充模型。在推导连续性方程和纳维-斯托克斯方程组时,并没有限制流动状态是层流还是湍流,因而其对层流和湍流同样成立[3]。对于环境气流马赫数$Ma\le 1$的小扰动,可采用线性不稳定性理论进行研究,实际上大多数喷雾应用都属于此范畴;但对于$Ma>1$的超声速强湍流,就要基于雷诺方程,采用非线性不稳定性理论,并考虑激波和气体的可压缩性进行分析,研究接受性问题,其数值解还可多支分叉,牵涉混沌问题。虽然目前已有基于雷诺方程的解析解研究,但讨论的是定常流进入静止气体环境中的简单模型[3]。对于液体射流,要探讨多级表面波、气液交界面相互垂直的瑞利波（Rayleigh wave）和泰勒波（Taylor wave）、表面波的时间（temporal）、空间（spatial or convective）和时空（absolute）模式、扰动的维数、线性和非线性不稳定性,逐步建立起液体喷射不稳定性的理论体系。
1 液体喷射的多级表面波
在对柴油机缸内超临界喷雾油束空间形态的观察中发现,在油束的边缘有一些梳子状结构[4],亚临界射流与超临界射流在空间形态上存在区别,超临界油束边缘的梳子状结构是由于亚临界油束边缘的卫星油滴气化造成的。要从理论上解释油束边缘卫星油滴的形成过程,则需要对射流的碎裂机理进行更为深入的研究。射流的多级表面波和瑞利-泰勒联合表面波理论的共同应用以期用来解释喷雾液束边缘析出卫星液滴的物理现象[5]。射流表面波是多级的,由此使得亚临界喷雾液束的边缘呈现出粗糙的表面。多级表面波是叠加而成的,第二级波附着在第一级波之上,而第三级波又附着在第二级波之上。该方式是由NAYFEH[6]、JAZAYERI等[7]建立的。多级叠加之后的气液交界面变得不再光滑,而是错落有致,粗糙不平。而粗糙不平表面的进一步碎裂最终将会导致卫星液滴的析出。目前,对于多级表面波的研究是采用非线性不稳定性理论进行的,根据表面波振幅解数值计算得到的叠加波形图还能够预测射流的碎裂长度和碎裂时间。由于非线性不稳定分析的推导过程异常繁复,仅第三级波的动力学边界条件已经有50余项,微分方程的求解就更加困难,目前最高研究到三级波。
2 气液交界面的瑞利波和泰勒波
射流流体团块在周围气体的扰动下会在气液体交界面处形成表面波,沿射流喷射方向的表面波称为瑞利波或者R波,沿横向或者旋转方向的表面波称为泰勒波或者T波。对于平面液膜沿液体喷射的$x$方向的是瑞利波,沿$y$方向的是泰勒波。通常,泰勒波的振幅要比瑞利波的小得多,因此目前大多数学者都忽略了射流气液交界面的泰勒波,而仅研究瑞利波。在这种情况下,对于平面液膜,沿$x$方向液膜是一个波动的曲面,而沿$y$方向液膜则是一个未经扰动的面。假设$a$为平面液膜在喷嘴出口处的半厚度（m）,$\xi $为表面波的振幅（m）。则$a+\xi $为瑞利波的波峰面,$a-\xi $为瑞利波的波谷面。当考虑有泰勒波时,沿$y$方向液膜是一个经过扰动的波形曲面。对于仅考虑瑞利波的情况,表面波为一光滑曲面,没有卫星液滴的形成。液膜仅在顶端裂开,形成带状断裂带,断裂带随后在液体表面张力的作用下聚集成棒状或线状,再碎裂成大量的离散液滴。对于圆射流,沿液体喷射的$z$方向的是瑞利波,沿旋转$\theta $方向的是泰勒波。与平面液膜同样,圆射流的泰勒波的振幅也要比瑞利波的小得多,沿θ方向是一个未经扰动的圆面,而沿$z$方向则是波动的曲面。假设a为圆射流在喷嘴出口处的半径,$\xi $为表面波的振幅。则圆面$a+\xi $为瑞利波的波峰面,圆面$a-\xi $为瑞利波的波谷面。当考虑有泰勒波时,沿θ向是一个经过扰动的波形圆面。对于环状液膜,沿θ方向是两个内环半径为${{r}_{\text{i}}}$、外环半径为${{r}_{\text{o}}}$的未经扰动的圆面,而沿$z$方向则是两个波动的内外环曲面。其中：圆面${{r}_{\text{i}}}+\xi $、${{r}_{\text{o}}}+\xi $为瑞利波的波峰面,圆面${{r}_{\text{i}}}-\xi $、${{r}_{\text{o}}}-\xi $为瑞利波的波谷面。当考虑有泰勒波时,沿θ方向是内外环两个经过扰动的波形圆面。环状液膜受环境气体的扰动作用,在喷嘴出口处就产生了波动,其碎裂长度比平面液膜的短。当不考虑泰勒波时,环状液膜的内外环均为光滑的曲面,没有卫星液滴形成。液膜在顶端碎裂形成环形断裂带,随后再碎裂成大量的细小液滴。RAYLEIGH[8]认为环形断裂带的厚度就等于液膜碎裂时顶端的厚度,宽度等于一个波长。在我们对瑞利-泰勒联合波的理论研究中发现,如果没有了瑞利波,泰勒波将不复存在;但如果没有了泰勒波,瑞利波仍然存在。说明瑞利波是主波,在液体射流的碎裂过程中起主导作用。
3 表面波的时间、空间和时空模式
对射流表面波的研究分为时间模式、空间模式和时空模式。时间模式表面波扰动表达式的形式与空间和时空模式的不同,而空间模式和时空模式扰动表达式的形式相同,只是各参数的含义不同。对于线性不稳定性理论,目前大多数学者研究的只是时间模式,对于空间模式和时空模式的研究较少。对于非线性不稳定性理论,则采用时空模式以及共轭复数方式进行研究。
对于一个复数指数函数$\exp \left( a+ib \right)$,如果虚部$b=0$,称${{\text{e}}^{a}}$为该指数函数的增长项。如果$a=0$,则称${{\text{e}}^{ib}}$为该指数函数的波动项。如果用$\exp \left( a+ib \right)$来描述液体射流气液交界面处的表面波行为,那么当$a$为正时,表面波振幅持续波动增长,射流变得越来越不稳定;当$a$为负时,表面波振幅持续波动减小,射流变得越来越稳定。波动项${{\text{e}}^{ib}}$的值介于 -1 ~ 1之间,即$-1\le {{\text{e}}^{ib}}\le 1$。当$b=0$时,${{\text{e}}^{ib}}=1$;而当$b=\text{ }\pi\text{ }$时,${{\text{e}}^{ib}}=-1$。如果增长项${{\text{e}}^{a}}$仅与时间$t$相关,即$a=f\left( t \right)$,则称之为时间模式;如果增长项${{\text{e}}^{a}}$仅与位移$x$相关,即$a=f\left( x \right)$,则称之为空间模式;如果增长项${{\text{e}}^{a}}$与时间$t$和位移$x$均相关,即$a=f\left( x,t \right)$,则称之为时空模式。在时间模式与空间模式中,$a$与$b$中的一个为复数,则另一个必须为实数。对于线性不稳定性分析,时间模式瑞利波的扰动表达式的一般形式为$\exp \left( \omega t+ikx \right)$。式中：对于直角坐标系的平面液膜,$\omega $称为复数特征频率（eigenfrequency）（r/s或s-1）,$\omega ={{\omega }_{r}}+i{{\omega }_{i}}$;$t$为时间（s）;$k$为实数波数（wave number）,简称波数（m-1）;$x$为位移（m）。将$\omega ={{\omega }_{r}}+i{{\omega }_{i}}$代入$\exp \left( \omega t+ikx \right)$,得扰动表达式$\exp \left( {{\omega }_{r}}t+i{{\omega }_{i}}t+ikx \right)$。式中：${{\omega }_{r}}$为时间轴表面波增长率（s-1）;$k$为空间轴波数,简称波数（m-1）;${{\omega }_{i}}$为时间轴波数,即实数特征频率（s-1）;推导得到的色散关系式为$f\left( {{\omega }_{r}},k \right)=0$和$f\left( {{\omega }_{r}},{{\omega }_{i}} \right)=0$,研究多采用前者。空间模式瑞利波的扰动表达式的一般形式为$\exp \left( i\omega t+ikx \right)$,式中：$\omega $为实数特征频率,$k$为复数波数。将实数$\omega $和复数$k={{k}_{r}}+i{{k}_{i}}$代入$\exp \left( i\omega t+ikx \right)$,得扰动表达式$\exp \left( -{{k}_{i}}x+i\omega t+i{{k}_{r}}x \right)$。则 $-{{k}_{i}}$为空间轴表面波增长率,由于${{k}_{i}}$本身为负,因此$-{{k}_{i}}$为正;${{k}_{r}}$为波数;$\omega $为实数特征频率,色散关系式为$f\left( -{{k}_{i}},{{k}_{r}} \right)=0$和$f\left( -{{k}_{i}},\omega \right)=0$,研究多采用前者。时空模式瑞利波扰动表达式的一般形式与空间模式相同,为$\exp \left( i\omega t+ikx \right)$,式中：$\omega $为复数特征频率,$k$为复数波数。将复数特征频率$\omega ={{\omega }_{r}}+i{{\omega }_{i}}$和复数波数$k={{k}_{r}}+i{{k}_{i}}$代入,得扰动表达式$\exp \left( -{{\omega }_{i}}t-{{k}_{i}}x+i{{\omega }_{r}}t+i{{k}_{r}}x \right)$。式中：$-{{k}_{i}}$为空间轴表面波增长率,$-{{\omega }_{i}}$为时间轴表面波增长率,同样因${{\omega }_{i}}$本身为负,故$-{{\omega }_{i}}$为正;${{\omega }_{r}}$为特征频率（s-1）,${{k}_{r}}$为波数。色散关系式分别取$f\left( -{{k}_{i}},{{\omega }_{r}},{{k}_{r}} \right)=0$和$f\left( -{{\omega }_{i}},{{\omega }_{r}},{{k}_{r}} \right)=0$,研究多采用前者。非线性不稳定性分析采用$\exp \left( i\omega t+ikx \right)$形式,与线性不稳定性分析对比时采用$\exp \left( -{{\omega }_{i}}t-{{k}_{i}}x+i{{\omega }_{r}}t+i{{k}_{r}}x \right)$形式,色散关系式通常取$f\left( -{{\omega }_{i}},{{k}_{r}} \right)=0$进行研究。时间模式仅有时间轴表面波增长率${{\omega }_{r}}$,而没有空间轴表面波增长率,表达了射流气液交界面表面波的时域特征;空间模式仅有空间轴表面波增长率$-{{k}_{i}}$,而没有时间轴表面波增长率,表达了表面波的空域特征。时空模式既有时间轴表面波增长率$-{{\omega }_{i}}$,又有空间轴表面波增长率$-{{k}_{i}}$;既有特征频率${{\omega }_{r}}$,又有波数${{k}_{r}}$,表达了表面波的时空域特征。
GASTER、LI、史绍熙的时间模式和空间模式扰动表达式不同,但时空模式的扰动表达式均相同,为$\exp \left( i\omega t+ikx \right)$。对于时间模式、空间模式和时空模式扰动表达式有三种划分方法：①GASTER[9]的时间模式和空间模式扰动表达式均为$\exp \left( i\omega t+ikx \right)$,则时间模式、空间模式和时空模式扰动表达式都是完全相同的;②LI等[10,11,12,13]的时间模式扰动表达式为$\exp \left( \omega t+ikx \right)$,空间模式和时空模式的是$\exp \left( i\omega t+ikx \right)$;③CHEN等[14]、史绍熙等[15]的时间模式和空间模式扰动表达式为$\exp \left( \omega t+ikx \right)$,时空模式的是$\exp \left( i\omega t+ikx \right)$。史绍熙还选用时间模式$f\left( {{\omega }_{r}},{{\omega }_{i}} \right)=0$与空间模式$f\left( -{{k}_{i}},{{\omega }_{r}} \right)=0$的数值计算结果进行了比较,显示正反对称波形的时间模式与空间模式没有明显的差别。其实,对于模式的划分是可以相互转化的,只是参数符号代表的含义不同而已,大多数学者采用的是第2种划分方法。GASTER和史绍熙分别从不同的表面波扰动表达式出发,分析了时空模式以及时间模式和空间模式之间的参数变换关系。
3.1 GASTER变换
假设扰动如实际射流表现的那样,既在时间上扰动与发展,又在空间上扰动与发展。以平面液膜表面波为例,将纵向位移,即振幅记为
$\xi ={{\xi }_{0}}\exp \left( i\omega t+ikx \right)$ (1)
式中：${{\xi }_{0}}$为喷嘴出口处的初始扰动振幅（m）。将时间模式记为$\left( \text{T} \right)$,空间模式$\left( \text{S} \right)$,时空模式$\left( \text{T,S} \right)$。则时空模式扰动表达式可以满足柯西-黎曼关系式
${{\omega }_{r,{{k}_{r}}}}={{\omega }_{i,{{k}_{i}}}} 和 {{\omega }_{r,{{k}_{i}}}}=-{{\omega }_{i,{{k}_{r}}}}$ (2)
式中：偏导数记为${{\omega }_{r,{{k}_{r}}}}=\frac{\partial {{\omega }_{r}}}{\partial {{k}_{r}}}$,后文同理。时间模式与空间模式的相互变换为
$k\left( \text{T} \right)={{k}_{r}}\left( \text{S} \right)$ (3)
${{\omega }_{r}}\left( \text{T} \right)=\omega \left( \text{S} \right)$ (4)
和
${{\omega }_{r,{{k}_{r}}}}=-\frac{{{\omega }_{i}}\left( \text{T} \right)}{{{k}_{i}}\left( \text{S} \right)}$ (5)
3.2 史绍熙变换
史绍熙时间模式和空间模式扰动振幅表达式采用
$\xi ={{\xi }_{0}}\exp \left( \omega t+ikx \right)$ (6)
则时空模式扰动表达式可以满足柯西-黎曼关系式（2）。时间模式与空间模式的相互变换为
$\omega \left( \text{S} \right)=i{{\omega }_{i}}\left( \text{S} \right)$ (7)
${{\omega }_{i}}\left( \text{T} \right)={{\omega }_{i}}\left( \text{S} \right)$ (8)
${{\omega }_{i,{{k}_{r}}}}=\frac{{{\omega }_{r}}\left( \text{T} \right)}{{{k}_{i}}\left( \text{S} \right)}$ (9)
应该注意的是,史绍熙空间模式对于参数的定义与GASTER的不同。其${{\omega }_{r}}\left( \text{S} \right)$对应于GASTER的$-{{\omega }_{i}}\left( \text{S} \right)$,而${{\omega }_{i}}\left( \text{S} \right)$对应于GASTER的${{\omega }_{r}}\left( \text{S} \right)=\omega \left( \text{S} \right)$。此外,对于三种模式中大多数学者所采用的第2类划分方法,时空模式的时间轴表面波增长率$-{{\omega }_{i}}$对应于时间模式的时间轴表面波增长率${{\omega }_{r}}$,即$-{{\omega }_{i}}\left( \text{T,S} \right)\tilde{\ }{{\omega }_{r}}\left( \text{T} \right)$;时空模式的波数${{k}_{r}}$对应于时间模式的波数$k$,即${{k}_{r}}\left( \text{T,S} \right)\tilde{\ }k\left( \text{T} \right)$;时空模式的特征频率${{\omega }_{r}}$对应于时间模式的特征频率${{\omega }_{i}}$,即${{\omega }_{r}}\left( \text{T,S} \right)\tilde{\ }{{\omega }_{i}}\left( \text{T} \right)$。时空模式的空间轴表面波增长率$-{{k}_{i}}$对应于空间模式的空间轴表面波增长率$-{{k}_{i}}$,即$-{{k}_{i}}\left( \text{T,S} \right)\tilde{\ }-{{k}_{i}}\left( \text{S} \right)$;时空模式的波数${{k}_{r}}$对应于空间模式的波数${{k}_{r}}$,即${{k}_{r}}\left( \text{T,S} \right)\tilde{\ }{{k}_{r}}\left( \text{S} \right)$;时空模式的特征频率${{\omega }_{r}}$对应于空间模式的特征频率$\omega $,即${{\omega }_{r}}\left( \text{T,S} \right)\tilde{\ }\omega \left( \text{S} \right)$。
实际上,在进行时间模式或空间模式色散准则关系式和稳定极限推导时,都不会将$\omega $或$k$写成$\left( {{\omega }_{r}}+i{{\omega }_{i}} \right)$或$\left( {{k}_{r}}+i{{k}_{i}} \right)$的复数形式。只是在编制数值计算程序时,时间模式要将$\omega $定义为双精度型复数,$k$定义为双精度型实数。空间模式要将$\omega $定义为双精度型实数,$k$定义为双精度型复数。时空模式要将$\omega $和$k$写成$\left( {{\omega }_{r}}+i{{\omega }_{i}} \right)$和$\left( {{k}_{r}}+i{{k}_{i}} \right)$的复数形式,并且要将${{\omega }_{r}}$、${{\omega }_{i}}$、${{k}_{r}}$、${{k}_{i}}$分别定义为双精度型实数。对于数值计算的复数结果,要根据模式选择复数的实部和/或虚部进行输出。
4 扰动的维数
平面液膜采用$x$、$y$、$z$方向的直角坐标系;圆射流和环状液膜采用$r$、$\theta $、$z$方向的圆柱坐标系。三维扰动的质量守恒方程、动量守恒方程、运动学边界条件、附加边界条件、动力学边界条件,推导得到扰动振幅、扰动压力和扰动速度表达式都是三维的。对于平面液膜,由于液膜很薄,即沿$z$方向上的厚度远远小于沿$x$方向上的长度以及沿$y$方向上的宽度。因此,液膜碎裂主要受到振幅在$z$方向上表面波的影响。对于圆射流和环状液膜,液束沿$r$方向的变形比沿$\theta $方向和$z$方向变形大得多[16,17,18]。一维扰动有两个动量守恒方程,一个质量守恒方程、一个运动学边界条件、一个附加边界条件、一个动力学边界条件,推导得到一个一维振幅、一个一维压力、一个一维速度表达式。三维扰动可简化为一个三维质量守恒方程、三个三维动量守恒方程、一个一维运动学边界条件、两个二维附加边界条件、一个一维动力学边界条件,推导得到一个一维振幅、一个一维压力、三个三维速度表达式。
5 液体喷射的线性和非线性不稳定性理论
5.1 线性和非线性不稳定性的区分
通常,流体运动基本方程还要加上初始条件和边界条件才能构成流体力学的定解问题[19]。流体运动所遵循的控制方程组是普适的,因此流动的个性就体现在初始条件和边界条件的差异上。初始条件是对不恒定流动指定初始时刻流场的某些流动参数。也就是说,能够满足某流场初始条件的流体流动形态可能有多种,流体流动方程可能有很多个,并不是唯一的。这些流体流动方程在某一设定的初始时刻均受限于流场的初始条件,但这些流体流动方程却是彼此不相同的,这就构成了满足初始条件的不恒定流动问题。边界条件是指运动方程的解在流场的边界上必须满足的运动学、附加和动力学条件。线性不稳定性理论对于控制方程组的定解只需要引入边界条件,即运动学边界条件、附加边界条件和动力学边界条件,而不必设置初始条件就能够解决流体不恒定流动的定解问题。由LI和JAZAYERI[7]建立的共轭复数模式非线性不稳定性理论对于连续性和动能守恒控制方程组的定解则除了要引入边界条件,即运动学和动力学边界条件以外,还要设置初始条件才能够解决流体不恒定流动的定解问题,否则定解的条件就不够。从他们的研究角度说明,非线性不稳定性理论不如线性不稳定性理论那样普适。而由曹建明和研究生王德超、舒力、张凯妹等建立的时空模式非线性不稳定性理论对于连续性和动能守恒控制方程组的定解仅引入边界条件,即运动学和动力学边界条件,而不必设置初始条件也能够求得气液相微分方程的定解,尽管时空模式的扰动振幅也同样能够满足共轭复数模式的初始条件。而且时空模式非线性不稳定性理论能够很好地与线性不稳定性理论相衔接。
线性不稳定性理论的研究对象是小扰动,对于小扰动来说,纳维-斯托克斯方程组中的非线性项可以忽略不计。动量守恒方程组中扰动速度对于位移或者旋转角度坐标一阶偏导数前面的系数如果是变量,则认为是非线性项;如果是常数,则认为是线性项。线性不稳定性理论即将动量守恒方程组中的非线性项直接删除,而保留线性项。该推导过程称为对动量守恒方程组的线性化。此外,在动力学边界条件中,表面波扰动振幅$\xi $对于顺流方向位移坐标$x$的一阶偏导数${{\xi }_{,x}}$表示$\xi $在$x$方向变化曲线的斜率。在射流喷射的不稳定性理论中,当${{\xi }_{,x}}=1$时,表面波的切线为一条45°直线。如果${{\xi }_{,x}}\le 1$,表面波波动曲线比较平缓,${{\xi }_{,x}}$可以被忽略,则认为扰动为小扰动,其表达式可以简化为线性的;如果${{\xi }_{,x}}>1$,$\xi $在$x$方向的曲线的斜率较大,则认为扰动不再是小扰动,而是大扰动,其表达式为非线性的。非线性表面波将比线性表面波更加不稳定,也更接近于液体喷射与雾化的实际情况。因此,${{\xi }_{,x}}=1$是区分小扰动与大扰动的界限,也是采用线性不稳定性分析与非线性不稳定性理论进行分析的界限。在推导动力学边界条件时,将${{\xi }_{,x}}$忽略掉的推导过程称为线性化。
5.2 时间和空间模式线性不稳定分析
通过对复数形式色散关系式的数值计算,可以对计算结果进行不稳定性分析。对于时间模式,色散关系式将分别得到$f\left( {{\omega }_{r}},k \right)=0$和$f\left( {{\omega }_{r}},{{\omega }_{i}} \right)=0$曲线,由于两条曲线几乎重合,通常取前者进行研究。曲线最高点所对应的表面波增长率为时间轴支配表面波增长率${{\omega }_{r-\text{dom}}}$,所对应的波数为支配波数${{k}_{\text{dom}}}$。${{\omega }_{r-\text{dom}}}-{{k}_{\text{dom}}}$就是射流的最不稳定工况点,该点处射流最易碎裂,其所具备的流动条件就是射流碎裂的必要条件。曲线与横坐标的远端交点为截断波数${{k}_{0}}$,或者称为稳定极限,它表示波数的范围。
对于空间模式,色散关系式为$f\left( -{{k}_{i}},{{k}_{r}} \right)=0$和$f\left( -{{k}_{i}},\omega \right)=0$,两条曲线也几乎重合,通常取前者进行研究。空间模式的支配工况点为$\left( -{{k}_{i-\text{dom}}} \right)-{{k}_{r-\text{dom}}}$。
5.3 时空模式线性不稳定分析
5.3.1 三维曲面图
对于时空模式,色散关系式为$f\left( -{{\omega }_{i}},{{\omega }_{r}},{{k}_{r}} \right)=0$和$f\left( -{{k}_{i}},{{\omega }_{r}},{{k}_{r}} \right)=0$,可以绘制出两个三维曲面图。每个曲面均有一个峰值点,其支配表面波增长率：时间轴为$-{{\omega }_{i-\text{dom}}}$,空间轴为$-{{k}_{i-\text{dom}}}$。支配波数分别是${{\left( {{\omega }_{r-\text{dom}}},{{k}_{r-\text{dom}}} \right)}_{\text{T}}}$或者${{\left( {{\omega }_{r-\text{dom}}},{{k}_{r-\text{dom}}} \right)}_{\text{S}}}$。如果绘制的是$f\left( {{\omega }_{i}},{{\omega }_{r}},{{k}_{r}} \right)=0$和$f\left( {{k}_{i}},{{\omega }_{r}},{{k}_{r}} \right)=0$三维曲面图,那么求取的就不是一个峰值点,而是鞍点,即曲面的最低点。以前可以绘制三维曲面图的等高线图（contour map）,再寻求等高线图的夹点（pinch point）,LI等[13]在进行时空模式不稳定分析中采用的就是这种方法,计算效率很低。现在计算机的运行速度已经飞速发展,无需采用绘制等高线图的方法来求取鞍点的数值解,直接绘制三维曲面图就可以了。然而,尽管计算机技术已经十分先进,但是数值计算仍然会面临不小的困难。在数值计算中,要预设一个公差数Nt。对于时间模式或空间模式的一维计算,预设Nt = 10-4。当采用穆勒方法前后两次相邻计算的相对误差小于等于Nt时,计算即被中止。也就是说,当两次计算的相对误差在万分之一以内时,就可以认为计算已经相当精确地得到了数值解。对于时空模式的三维计算,计算网格划分得越细,则越接近于数值解。要想达到与一维数值计算同样的计算精度,三维网格数就要划分为${{10}^{13}}\tilde{\ }{{10}^{14}}$,少了会使得网格过粗,极有可能漏掉数值解;如果采用网格搜索方法缩小网格范围,网格分区的计算量也很庞大。即使对于鞍点有所预估而减少分区的搜索数目,也有可能搜索到谬根。如此巨大的计算量,目前的计算机很难完成。因此,三维曲面图很难绘制。我们将采用穆勒方法,得到数值解的三维空间曲线图,以代替三维空间曲面图。BERS[20]认为三维曲线数值解与三维曲面解等效,可以用来分析时空模式不稳定性。
5.3.2 波动的相速度与群速度
指数函数可以与三角函数在复数范围相互转换,并且以正弦函数和余弦函数表示。对于正/余弦波,相速度（phase velocity）的定义为：单一频率波的位相面在介质中的传播速度。对于扰动振幅为$\xi ={{\xi }_{0}}\exp \left[ i\left( \omega t+kx \right) \right]$的液体射流,相速度定义式为
${{v}_{\text{p}}}={{x}_{,t}}=\frac{\omega }{k}$ (10)
群速度（group velocity）定义为波包的包络在介质中的增长速度。群速度定义式为
${{v}_{\text{g}}}={{\omega }_{,k}}$ (11)
将方程（11）代入方程（10）,可得群速度与相速度之间的关系
${{v}_{\text{g}}}=\frac{{{v}_{\text{p}}}}{1-\frac{\omega }{{{v}_{\text{p}}}}{{v}_{\text{p},\omega }}}$ (12)
按照色散类型可以分为：①${{v}_{\text{p,}\omega }}<0$,则${{v}_{\text{g}}}<{{v}_{\text{p}}}$,称为正常色散;②${{v}_{\text{p},\omega }}>0$,则${{v}_{\text{g}}}>{{v}_{\text{p}}}$,称为反常色散;③${{v}_{\text{p},\omega }}=0$,则${{v}_{\text{g}}}={{v}_{\text{p}}}$,称为无色散。当群速度和相速度均为变量时,为正常色散或者反常色散;当群速度与相速度均为固定的常数时,那么${{v}_{\text{p},\omega }}=0$。根据方程（12）,则${{v}_{\text{g}}}={{v}_{\text{p}}}=\text{c}$（式中$\text{c}$为常数）,为无色散。例如,对于群速度的均值${{\bar{v}}_{\text{g}}}$和相速度的均值${{\bar{v}}_{\text{p}}}$,必然有${{\bar{v}}_{\text{g}}}={{\bar{v}}_{\text{p}}}=\text{c}$。
5.3.3 三维曲线图
时空模式的三维曲线图的纵坐标为时间轴表面波增长率或者空间轴表面波增长率$-{{\omega }_{i}}/-{{k}_{i}}$,横坐标为特征频率${{\omega }_{r}}$和波数${{k}_{r}}$。在BERS[20]的论文中,纵坐标为波的增长率${{{\omega }'}_{0i}}$,横坐标为速度${{V}_{x}}$和${{V}_{y}}$。此处角标“0”表示支配参数“dom”。${\omega }'=\omega -kV$,在“dom”点,${{v}_{\text{g}0}}={{v}_{\text{p}0}}=\frac{{{\omega }_{0i}}}{{{k}_{0}}}=\text{c}$,则$kV={{k}_{0}}{{V}_{0}}={{\omega }_{0i}}$。那么,${{{\omega }'}_{\text{0}i}}={{\omega }_{i}}-{{\omega }_{\text{0}i}}$。式中：${{\omega }_{\text{0}i}}$为常数,${{{\omega }'}_{\text{0}i}}$和${{\omega }_{i}}$均为变量。可以进行纵坐标变换为$-{{\omega }_{i}}=-\left( {{{{\omega }'}}_{\text{0}i}}+{{\omega }_{\text{0}i}} \right)$,同理,也可以对${{{k}'}_{0i}}$进行纵坐标变换为$-{{k}_{i}}$。横坐标${{V}_{x}}$和${{V}_{y}}$可以为时间坐标和空间坐标。对于射流,$V={{v}_{\text{g}}}$。由于相速度与群速度具有函数关系（12）,即${{v}_{\text{p}}}=f\left( {{v}_{\text{g}}} \right)$,因此横坐标${{v}_{\text{g}t}}$和${{v}_{\text{g}x}}$可以变换为${{v}_{\text{p}t}}$和${{v}_{\text{p}x}}$。由于特征频率${{\omega }_{r}}$与相速度${{v}_{\text{p}t}}$成正比,波数${{k}_{r}}$与相速度${{v}_{\text{p}x}}$正相关。因此,横坐标又可以从${{v}_{\text{p}t}}$和${{v}_{\text{p}x}}$再变换为${{\omega }_{r}}$和${{k}_{r}}$。曲线$f\left( -{{\omega }_{i}},{{\omega }_{r}},{{k}_{r}} \right)=0$和曲线$f\left( -{{k}_{i}},{{\omega }_{r}},{{k}_{r}} \right)=0$的初终点截距分别为均速${{\bar{v}}_{\text{g}1}}={{\bar{v}}_{\text{p}1}}=\frac{{{\omega }_{r1}}}{{{k}_{r1}}}$和${{\bar{v}}_{\text{g}2}}={{\bar{v}}_{\text{p}2}}=\frac{{{\omega }_{r2}}}{{{k}_{r2}}}$。射流的数值计算结果显示,通常初点的${{\omega }_{r1}}={{k}_{r1}}\approx 0$,则${{\bar{v}}_{\text{g}1}}={{\bar{v}}_{\text{p}1}}\approx 0$。终点坐标与时间轴截断特征频率${{\omega }_{r0}}$和空间轴截断波数${{k}_{r0}}$有关,即${{\bar{v}}_{\text{g}2}}=\sqrt{\omega _{r0}^{2}+k_{r0}^{2}}$。
对于$f\left( -{{\omega }_{i}},{{\omega }_{r}},{{k}_{r}} \right)=0$图,每固定一个$-{{k}_{i}}$,就可获得一个$\left( -{{\omega }_{i-\text{dom}}} \right)-{{\left( {{\omega }_{r-\text{dom}}},{{k}_{r-\text{dom}}} \right)}_{\text{T}}}$点。改变$-{{k}_{i}}$的值,就可以得到一组$\left( -{{\omega }_{i-\text{dom}}} \right)-{{\left( {{\omega }_{r-\text{dom}}},{{k}_{r-\text{dom}}} \right)}_{\text{T}}}$点,绘制出$\left( -{{\omega }_{i-\text{dom}}} \right)-\left( -{{k}_{i}} \right)$曲线图。图中曲线有一个特殊点,即$-{{k}_{i}}=0$的点,它是时间模式的点。也就是说,时间模式是时空模式的一个特例。同理,对于$f\left( -{{k}_{i}},{{\omega }_{r}},{{k}_{r}} \right)=0$图,每固定一个$-{{\omega }_{i}}$,就可以得到一个$\left( -{{k}_{i-\text{dom}}} \right)-{{\left( {{\omega }_{r-\text{dom}}},{{k}_{r-\text{dom}}} \right)}_{\text{S}}}$点。改变$-{{\omega }_{i}}$的值,就可以得到一组$\left( -{{k}_{i-\text{dom}}} \right)-{{\left( {{\omega }_{r-\text{dom}}},{{k}_{r-\text{dom}}} \right)}_{\text{S}}}$点,绘制出$\left( -{{k}_{i-\text{dom}}} \right)-\left( -{{\omega }_{i}} \right)$曲线图。图中曲线也有一个特殊点,即$-{{\omega }_{i}}=0$的点,它是空间模式的点。也就是说,空间模式也是时空模式的一个特例。 $\left( -{{\omega }_{i-\text{dom}}} \right)-\left( -{{k}_{i}} \right)$和$\left( -{{k}_{i-\text{dom}}} \right)-\left( -{{\omega }_{i}} \right)$曲线图的峰值点就是时空模式的最不稳定工况点。通过改变计算过程中的流动参数,如液流雷诺数$R{{e}_{\text{l}}}$、韦伯数$W{{e}_{_{\text{l}}}}$、欧拉数$E{{u}_{_{\text{l}}}}$、欧尼索数$O{{h}_{_{\text{l}}}}$、马赫数$M{{a}_{\text{l}}}$、气流马赫数$M{{a}_{\text{g}}}$、气液流速比$U$、气液密度比$\rho $、气液压力比$P$等,可以分析这些参数对时间模式、空间模式和时空模式碎裂过程的影响。
5.4 非线性不稳定分析
非线性不稳定性分析是要对第一级波、第二级波和第三级波的支配参数${{\bar{\omega }}_{i1/2/3-\text{dom}}}$、${{\bar{\omega }}_{r1/2/3-\text{dom}}}$、${{\bar{k}}_{r1/2/3-\text{dom}}}$进行计算,得到第一级波、第二级波和第三级波的扰动振幅初始函数表达式,再分别乘以${{\bar{\xi }}_{0}}$、$\bar{\xi }_{0}^{2}$和$\bar{\xi }_{0}^{3}$,即可得到第一级波、第二级波和第三级波的扰动振幅值。根据各时间点三级波的扰动振幅值,绘制第一级波、第二级波和第三级波的波形图,以及三级波的波形叠加图,得到射流的碎裂长度和碎裂时间。
6 结语
国际上对于喷雾理论的研究已经有100多年了,发表的论文也不少,但系统全面的研究尚未见报道。理论上,数值计算可以适用于任何工况,除非是计算设备不能满足计算需求,例如存储空间和/或运行速度不够,以及计算数据溢出等情况。其数值计算范围要比试验研究广泛得多。由于自由射流碎裂过程的影响因素很多,且对一些影响因素十分敏感。因此,理论研究应将重点放在揭示射流规律性的普适原理之上,以起到指导实践的深层次作用。本文在阅读了大量文献的基础上,集二十余年的不断积累,做了一点归纳总结,希望能够供国内外从事喷雾理论的研究者参考。