1. Introduction
1.1. Literature review
The instability of planar liquid sheets is a classical problem of fluid mechanics, with broad practical applications, such as atomization, combustion and paper-making. A good understanding of it is of both scientific interest and pragmatic usage. Therefore, a number of researchers have studied it in various aspects.
Squire (Reference Squire1953) and Hagerty & Shea (Reference Hagerty and Shea1955) investigated the linear instability of inviscid sheets. They obtained the dispersion relation and found that the instability of the sheet is caused mainly by its interaction with the surrounding gas; increasing the Weber number and the gas-to-liquid density ratio extensively enhances instability. There are two unstable modes in linear instability, i.e. the sinuous mode (antisymmetric waves) and the varicose or dilatational mode (symmetric waves), schematically shown in figure 1. Both their theories and experiments indicated that, in common conditions, the sinuous waves grow much faster than the varicose correlate.
Figure 1. Linear sinuous mode (a) and varicose mode (b).
Linear analysis only predicts the onset of instability, and depicting the final breakup of the sheet requires nonlinear analysis. The asymptotic analytical perturbation solution for weakly nonlinear instability of the linear sinuous mode was investigated in the study of Clark & Dombrowski (Reference Clark and Dombrowski1972), Jazayeri & Li (Reference Jazayeri and Li2000) and Yang et al. (Reference Yang, Wang, Fu, Du and Tong2013). It was found that the first harmonic of the linear sinuous mode is varicose, which explains the disintegration of sinuous waves. The harmonics of the linear varicose mode were not addressed in these studies, because the growth rate of the varicose wave is much smaller and the sinuous mode was believed to be the major reason for the disintegration. The nonlinear instability of the linear varicose mode was studied analytically or semi-analytically by Matsuuchi (Reference Matsuuchi1974, Reference Matsuuchi1976) and Mehring & Sirignano (Reference Mehring and Sirignano1999). They found that, although a sheet in a vacuum is stable according to linear analysis, it can be modulationally unstable when the nonlinearity is taken into account. However, partly because they neglected the ambient gas, which plays a vital role in the growth of a disturbance, their studies could not show the sheet disintegration caused by a small initial disturbance. The demands for very small disturbance amplitude (Matsuuchi Reference Matsuuchi1974) and very long waves (Matsuuchi Reference Matsuuchi1976; Mehring & Sirignano Reference Mehring and Sirignano1999) are also their limitations. The numerical approach to this issue was undertaken by Rangel & Sirignano (Reference Rangel and Sirignano1991) and Kan & Yoshinaga (Reference Kan and Yoshinaga2007). By means of boundary integral methods, they simulated the nonlinear distortion and disintegration for both the sinuous and varicose modes.
Because the linear sinuous mode grows much faster than the varicose counterpart in common conditions, many researchers believe that the sinuous mode is the dominant mode that determines the instability behaviour, including the disturbance growth and breakup. In the theoretical analysis of Squire (Reference Squire1953), Clark & Dombrowski (Reference Clark and Dombrowski1972), Asare, Takahashi & Hoffman (Reference Asare, Takahashi and Hoffman1981), Jazayeri & Li (Reference Jazayeri and Li2000) and Yang et al. (Reference Yang, Wang, Fu, Du and Tong2013), only the linear sinuous mode was investigated, and in the experiments of Squire (Reference Squire1953), Hagerty & Shea (Reference Hagerty and Shea1955), Crapper, Dombrowski & Pyott (Reference Crapper, Dombrowski and Pyott1975), Asare et al. (Reference Asare, Takahashi and Hoffman1981) and Tammisola et al. (Reference Tammisola, Sasaki, Lundell, Matsubara and Söderberg2011), the waves were evidently characterized by sinuous characters. However, Mitra, Li & Renksizbulut (Reference Mitra, Li and Renksizbulut2001) offered a different opinion. In their dual-mode linear instability analysis, the linear varicose mode was the only reason for disintegration because the linear sinuous mode does not contribute to the thinning. Since using the linear analysis to study the breakup is certainly dubious and inappropriate, their argument has rarely been adopted by other scholars. But their opinion is right on at least two points: first, at the onset of instability where nonlinear amplitudes are negligible, the varicose mode has a dominant effect in causing the thinning; and second, the linear sinuous mode itself cannot lead to the breakup anyway, however large the amplitude to which it grows. Regarding the breakup, the core of the controversy lies mainly in the question of which is more important in causing the attenuation: the varicose harmonics of the linear sinuous mode, which grow very fast but have very small initial amplitudes, or the linear varicose mode, which has a much larger initial amplitude but grows slowly. Such a problem cannot be answered by previous studies, and therefore it will be addressed in the present study of nonlinear dual-mode instability.
1.2. Present work
The main objective of the present work is to study the nonlinear instability of liquid sheets leading to breakup whose linear instability contains both sinuous and varicose modes. This issue has rarely been studied so far: only Mitra et al. (Reference Mitra, Li and Renksizbulut2001) have considered the sheet’s dual-mode linear instability, but since the linear instability is not justified to study the breakup, their results are not genuinely convincing. The topic is important in several respects. For practical application, where any forms of disturbance can exist, that the initial disturbance is purely sinuous or purely varicose is only an ideal assumption, and if the initial phase difference between the sheet’s two surfaces is neither 0 nor
${\rm\pi}$
, both modes will exist. So the present study, which considers both modes, is an important attempt to approach real situations. More merits of it lie in its scientific interest. For the sheet breakup, as noted in the literature review, there have been controversial opinions on the roles of the sinuous mode and varicose mode. Which is dominant? While many previous studies believe that considering only the fast-growing sinuous mode is enough, our current analysis indicates that this is not appropriate, at least in the conditions considered in this paper. By examining the case where the linear sinuous and varicose modes have identical initial amplitudes, we demonstrate that the linear varicose mode plays a vital role in the breakup – more important than the sinuous mode for most conditions – through both itself and its first harmonic. Another important point of this issue is the nonlinear coupling between linear unstable modes, characterized by the first harmonic generated by the coupling between the linear sinuous and varicose modes. This unstable mode exists only when both the linear sinuous and varicose modes exist, indicating that the two linear modes do not propagate independently when nonlinearity is taken into account. While it is sinuous in nature and does not contribute to breakup, it has an important function in modulating the shapes of wave profiles, leading to many interesting phenomena.
Another important task of this study is the comparison between weakly nonlinear analysis and numerical simulation. This issue also has not been addressed by previous studies for the aerodynamic instability of sheets. Clark & Dombrowski (Reference Clark and Dombrowski1972), Jazayeri & Li (Reference Jazayeri and Li2000) and Yang et al. (Reference Yang, Wang, Fu, Du and Tong2013) only undertook the weakly nonlinear analysis using a perturbation technique, so the question exists whether their results are appropriate to study the breakup when nonlinearity may become significant. The numerical simulations of Rangel & Sirignano (Reference Rangel and Sirignano1991) and Kan & Yoshinaga (Reference Kan and Yoshinaga2007) were compared only to linear instability analysis, so the breakup mechanism could not be explained by the analytical solution. Only Mehring & Sirignano (Reference Mehring and Sirignano1999) compared their analytical and numerical results, but their study was for the liquid sheet in a vacuum, in which the disturbance amplitude does not grow and the sheet breaks up only at large initial amplitudes. This is different from the aerodynamic instability considered in our study, where small initial amplitudes grow to large ones and cause the breakup. While the numerical simulation provides a relatively precise solution, the weakly nonlinear analysis gives results quickly and reliably, and, more importantly, it explains the mechanism of instability. By comparing their results, we show essential agreements between the two methods, and thus the validity of weakly nonlinear analysis is justified and the mechanism of numerical solution is explained.
2. Mathematical solution
2.1. Problem description
A two-dimensional planar liquid sheet moving through ambient gas is considered. Both the liquid and the gas are assumed to be inviscid and incompressible. The surface tension of the liquid sheet is
${\it\sigma}$
; the densities of the gas and liquid are
${\it\rho}_{g}$
and
${\it\rho}_{l}$
, respectively. For the basic flow, the sheet has a constant thickness
$2a$
and a uniform velocity
$U$
, and the ambient gas is stationary. Gravity is neglected, as we consider the sheet to be moving at a high velocity. To generalize the results, all the variables and parameters are non-dimensionalized. The non-dimensional Weber number and gas-to-liquid density ratio are defined as
$\mathit{We}={\it\rho}_{l}U^{2}a/{\it\sigma}$
and
${\it\rho}={\it\rho}_{g}/{\it\rho}_{l}$
. The time, length, velocity and velocity potential are scaled with
$a/U$
,
$a$
,
$U$
and
$Ua$
, respectively. It is noted that the viscosity of liquid was considered in our previous weakly nonlinear analysis for the pure sinuous mode (Yang et al.
Reference Yang, Wang, Fu, Du and Tong2013), but in the current paper, principally because of the limitation of numerical simulation (see § 2.3), the liquid is assumed inviscid as well.
Figure 2. Sketch diagram of a moving liquid sheet.
A sketch diagram of the sheet is schematically shown in figure 2. In Cartesian coordinate, the basic flow direction is parallel to the
$x$
-direction and the
$y$
-direction is normal to the undisturbed gas-to-liquid interface. The
$x$
-axis is the centreline of the undisturbed sheet, and thus the surfaces of the sheet for the basic flow are located at
$y=\pm 1$
. When the flow is disturbed, the non-dimensional interface displacements in the
$y$
-direction are defined as
${\it\eta}_{j}(x,t)$
, where
$j=1$
represents the upper interface and
$j=2$
represents the lower, and therefore the locations of the disturbed interfaces are
$y=(-1)^{j+1}+{\it\eta}_{j}(x,t)$
.
In the present study, only temporal instability is considered, so the sheet is assumed to be infinite, and at
$t=0$
an initial interface displacement is given as
$$\begin{eqnarray}{{\it\eta}_{j}|}_{t=0}={\it\eta}_{0}\hat{{\it\eta}}_{j}\exp (\text{i}kx)+\text{c.c.},\end{eqnarray}$$
where
${\it\eta}_{0}$
is a small number representing the initial disturbance amplitude,
$k$
is a real and positive number referred to as the wavenumber and ‘c.c.’ represents the complex conjugate. Here
$\hat{{\it\eta}}_{j}~(j=1,2)$
are complex numbers indicating the difference between the upper and lower interfaces. Both
$\hat{{\it\eta}}_{1}$
and
$\hat{{\it\eta}}_{2}$
can be independently given with any arbitrary values, and if they are not exactly the same (sinuous mode), or opposite (varicose mode), both linear sinuous and varicose modes exist in the sheet’s linear instability. Note however that, because the initial velocity is not given, the solution of the problem is not fixed. The focus of our study is one particular solution whose linear part can be expressed by the sum of two normal modes (2.21); the instability behaviour of other solutions may be studied in future research.
Two different approaches are used to solve this problem: the weakly nonlinear analysis using a perturbation expansion, and numerical simulation using a boundary integral method.
2.2. Weakly nonlinear analysis
A weakly nonlinear analysis using a second-order perturbation expansion is undertaken to obtain an approximate analytical solution. The weakly nonlinear solution can be precise for most of the time of disturbance evolution. It gives results quickly and reliably and, more importantly, it provides insight into the mechanism of instability and breakup.
2.2.1. Governing equations
Assuming that both the liquid and gas are initially irrotational, the non-dimensional velocity potentials of the liquid and gas phases are defined as
${\it\phi}_{l}$
and
${\it\phi}_{gj}$
, where
$j=1$
represents the gas above the sheet and
$j=2$
represents the gas below.
The mass conservation equations of the liquid and gas are
$$\begin{eqnarray}\displaystyle & {\it\phi}_{l,xx}+{\it\phi}_{l,yy}=0\quad \text{for }-1+{\it\eta}_{2}\leqslant y\leqslant 1+{\it\eta}_{1}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & {\it\phi}_{gj,xx}+{\it\phi}_{gj,yy}=0\quad j=1,2\text{ for }y\geqslant 1+{\it\eta}_{1},y\leqslant -1+{\it\eta}_{2}. & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & {\it\phi}_{l,y}-{\it\eta}_{j,t}-{\it\phi}_{l,x}{\it\eta}_{j,x}=0\quad \text{at }y=(-1)^{j+1}+{\it\eta}_{j}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & {\it\phi}_{gj,y}-{\it\eta}_{j,t}-{\it\phi}_{gj,x}{\it\eta}_{j,x}=0\quad \text{at }y=(-1)^{j+1}+{\it\eta}_{j}. & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle \frac{1}{2}+{\it\rho}{\it\phi}_{gj,t}-{\it\phi}_{l,t}+\frac{1}{2}{\it\rho}({\it\phi}_{gj,x}^{2}+{\it\phi}_{gj,y}^{2})-\frac{1}{2}({\it\phi}_{l,x}^{2}+{\it\phi}_{l,y}^{2})-\frac{(-1)^{j}{\it\eta}_{j,xx}}{\mathit{We}(1+{\it\eta}_{j,x}^{2})^{1.5}}=0\nonumber\\ \displaystyle & & \displaystyle \quad \text{at }y=(-1)^{j+1}+{\it\eta}_{j}.\end{eqnarray}$$
A second-order weakly nonlinear analysis is used to solve the nonlinear equations. Assuming that the nonlinearity is weak, all the variables are expressed by a second-order expansion in power series of
${\it\eta}_{0}$
:
$$\begin{eqnarray}({\it\phi}_{l},{\it\phi}_{gj},{\it\eta}_{j})=(x,0,0)+(\text{}^{1}{\it\phi}_{l},\text{}^{1}{\it\phi}_{gj},\text{}^{1}{\it\eta}_{j}){\it\eta}_{0}+(\text{}^{2}{\it\phi}_{l},\text{}^{2}{\it\phi}_{gj},\text{}^{2}{\it\eta}_{j}){\it\eta}_{0}^{2}+O({\it\eta}_{0}^{3}).\end{eqnarray}$$
Note that the terms in the first bracket on the right-hand side of (2.7) represent the basic flow field. For the variables on the disturbed boundaries, their values are obtained by a Taylor series using their values on the undisturbed boundaries, for example,
$$\begin{eqnarray}{{\it\phi}_{l,t}|}_{y=(-1)^{j+1}+{\it\eta}_{j}}={{\it\phi}_{l,t}|}_{y=(-1)^{j+1}}+{\it\eta}_{j}{{\it\phi}_{l,yt}|}_{y=(-1)^{j+1}}+\cdots \,.\end{eqnarray}$$
Substituting (2.7) and (2.8) into (2.1)–(2.6), collecting the terms with
${\it\eta}_{0}$
and
${\it\eta}_{0}^{2}$
, the first-order and the second-order governing equations are achieved.
The first-order equations are
$$\begin{eqnarray}\displaystyle & {\text{}^{1}{\it\eta}_{j}|}_{t=0}=\hat{{\it\eta}}_{j}\exp (\text{i}kx)+\text{c.c.}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \text{}^{1}{\it\phi}_{l,xx}+\text{}^{1}{\it\phi}_{l,yy}=0\quad \text{for }-1\leqslant y\leqslant 1, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \text{}^{1}{\it\phi}_{gj,xx}+\text{}^{1}{\it\phi}_{gj,yy}=0\quad j=1,2\text{ for }y\geqslant 1,y\leqslant -1, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \text{}^{1}{\it\phi}_{l,y}-\text{}^{1}{\it\eta}_{j,t}-\text{}^{1}{\it\eta}_{j,x}=0\quad \text{at }y=(-1)^{j+1}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \text{}^{1}{\it\phi}_{gj,y}-\text{}^{1}{\it\eta}_{j,t}=0\quad \text{at }y=(-1)^{j+1}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\it\rho}\text{}^{1}{\it\phi}_{gj,t}-\text{}^{1}{\it\phi}_{l,t}-\text{}^{1}{\it\phi}_{l,x}-(-1)^{j}\frac{\text{}^{1}{\it\eta}_{j,xx}}{\mathit{We}}=0\quad \text{at }y=(-1)^{j+1}. & \displaystyle\end{eqnarray}$$
The second-order equations are
$$\begin{eqnarray}\displaystyle & {\text{}^{2}{\it\eta}_{j}|}_{t=0}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \text{}^{2}{\it\phi}_{l,xx}+\text{}^{2}{\it\phi}_{l,yy}=0\quad \text{for }-1\leqslant y\leqslant 1, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \text{}^{2}{\it\phi}_{gj,xx}+\text{}^{2}{\it\phi}_{gj,yy}=0\quad j=1,2\text{ for }y\geqslant 1,y\leqslant -1, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \text{}^{2}{\it\phi}_{l,y}-\text{}^{2}{\it\eta}_{j,t}-\text{}^{2}{\it\eta}_{j,x}=\text{}^{1}{\it\phi}_{l,x}\text{}^{1}{\it\eta}_{j,x}-\text{}^{1}{\it\eta}_{j}\text{}^{1}{\it\phi}_{l,yy}\quad \text{at }y=(-1)^{j+1}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \text{}^{2}{\it\phi}_{gj,y}-\text{}^{2}{\it\eta}_{j,t}=\text{}^{1}{\it\phi}_{gj,x}\text{}^{1}{\it\eta}_{j,x}-\text{}^{1}{\it\eta}_{j}\text{}^{1}{\it\phi}_{gj,yy}\quad \text{at }y=(-1)^{j+1}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle {\it\rho}\text{}^{2}{\it\phi}_{gj,t}-\text{}^{2}{\it\phi}_{l,t}-\text{}^{2}{\it\phi}_{l,x}-(-1)^{j}\frac{\text{}^{2}{\it\eta}_{j,xx}}{\mathit{We}}\nonumber\\ \displaystyle & & \displaystyle \quad =\text{}^{1}{\it\eta}_{j}(-{\it\rho}\text{}^{1}{\it\phi}_{gj,yt}+\text{}^{1}{\it\phi}_{l,yt}+\text{}^{1}{\it\phi}_{l,yx})-\frac{1}{2}{\it\rho}(\text{}^{1}{\it\phi}_{gj,x}^{2}+\text{}^{1}{\it\phi}_{gj,y}^{2})+\frac{1}{2}(\text{}^{1}{\it\phi}_{l,x}^{2}+\text{}^{1}{\it\phi}_{l,y}^{2})\nonumber\\ \displaystyle & & \displaystyle \qquad \text{at }y=(-1)^{j+1}.\end{eqnarray}$$
2.2.2. First-order solutions
The first-order equations are linear in
$x$
and
$t$
; therefore, the sinuous mode and the varicose mode can be added together directly to give
$$\begin{eqnarray}\displaystyle (\text{}^{1}{\it\phi}_{l},\text{}^{1}{\it\phi}_{gj},\text{}^{1}{\it\eta}_{j}) & = & \displaystyle (\text{}^{1,s}\hat{{\it\phi}}_{l},\text{}^{1,s}\hat{{\it\phi}}_{gj},\text{}^{1,s}\hat{{\it\eta}}_{j})\exp (\text{i}kx-\text{i}{\it\omega}_{1s}t)\nonumber\\ \displaystyle & & \displaystyle +\,(\text{}^{1,v}\hat{{\it\phi}}_{l},\text{}^{1,v}\hat{{\it\phi}}_{gj},\text{}^{1,v}\hat{{\it\eta}}_{j})\exp (\text{i}kx-\text{i}{\it\omega}_{1v}t)+\text{c.c.},\end{eqnarray}$$
where
${\it\omega}_{1s}$
and
${\it\omega}_{1v}$
are the first-order sinuous and varicose complex frequencies; their real parts represent the oscillation frequencies and their imaginary parts represent the temporal growth rates. The superscript ‘
$1,s$
’ represents the linear sinuous mode and ‘
$1,v$
’ represents the linear varicose mode. The symbol ‘
$\hat{~}$
’ denotes the disturbance amplitudes, which are functions of
$y$
only and, in particular,
$\text{}^{1,s}\hat{{\it\eta}}_{j}$
and
$\text{}^{1,v}\hat{{\it\eta}}_{j}$
are constants.
The two unstable modes are solved using standard linear analysis combining (2.21) and (2.10)–(2.14). Here only the results are presented. The upper and lower interface displacements are identical for the sinuous mode and opposite for the varicose mode:
$$\begin{eqnarray}\text{}^{1,s}\hat{{\it\eta}}_{1}=\text{}^{1,s}\hat{{\it\eta}}_{2},\quad \text{}^{1,v}\hat{{\it\eta}}_{1}=-\text{}^{1,v}\hat{{\it\eta}}_{2}.\end{eqnarray}$$
The dispersion relations must be satisfied in linear instability to make the solution non-trivial:
$$\begin{eqnarray}D_{s}({\it\omega}_{1s},k)=0,\quad D_{v}({\it\omega}_{1v},k)=0.\end{eqnarray}$$
Here
$D_{s}$
and
$D_{v}$
are defined according to the character equations for the linear sinuous and varicose modes:
$$\begin{eqnarray}D_{s}({\it\omega},k)=-{\it\rho}{\it\omega}^{2}+\frac{k^{3}}{\mathit{We}}-({\it\omega}-k)^{2}\tanh k,\quad D_{v}({\it\omega},k)=-{\it\rho}{\it\omega}^{2}+\frac{k^{3}}{\mathit{We}}-({\it\omega}-k)^{2}\coth k.\end{eqnarray}$$
The velocity potentials are
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \text{}^{1,s}\hat{{\it\phi}}_{l}=\text{i}\left(1-\frac{{\it\omega}_{1s}}{k}\right)\frac{\sinh ky}{\cosh k}\text{}^{1,s}\hat{{\it\eta}}_{1},\quad \text{}^{1,s}\hat{{\it\phi}}_{gj}=(-1)^{j+1}\frac{\text{i}{\it\omega}_{1s}}{k}\exp [k+(-1)^{j}ky]\text{}^{1,s}\hat{{\it\eta}}_{1},\\ \displaystyle \text{}^{1,v}\hat{{\it\phi}}_{l}=\text{i}\left(1-\frac{{\it\omega}_{1v}}{k}\right)\frac{\cosh ky}{\sinh k}\text{}^{1,v}\hat{{\it\eta}}_{1},\quad \text{}^{1,v}\hat{{\it\phi}}_{gj}=\frac{\text{i}{\it\omega}_{1v}}{k}\exp [k+(-1)^{j}ky]\text{}^{1,v}\hat{{\it\eta}}_{1}.\end{array}\!\!\right\}\qquad & & \displaystyle\end{eqnarray}$$
Equations (2.22)–(2.25) are the well-known results of standard linear analysis and they are identical to those presented by Squire (Reference Squire1953) and Hagerty & Shea (Reference Hagerty and Shea1955). Note that (2.24) and (2.25) are valid only at
$k>0$
;
$k<0$
renders different expressions for the solutions but the same instability behaviour governed by the real parts of the complex solutions; at
$k=0$
, all the temporal frequencies and velocity potentials are zero. The unique solution of the linear homogeneous equations is determined by the initial condition. Substituting (2.21) and (2.22) into the initial condition (2.9), the disturbance amplitudes of the first-order sinuous and varicose modes are derived as
$$\begin{eqnarray}\text{}^{1,s}\hat{{\it\eta}}_{j}={\textstyle \frac{1}{2}}(\hat{{\it\eta}}_{1}+\hat{{\it\eta}}_{2}),\quad \text{}^{1,v}\hat{{\it\eta}}_{j}=(-1)^{j+1}{\textstyle \frac{1}{2}}(\hat{{\it\eta}}_{1}-\hat{{\it\eta}}_{2}).\end{eqnarray}$$
To measure the proportion of initial sinuous amplitude in the linear part, we define
$$\begin{eqnarray}r_{s}=\frac{|\text{}^{1,s}\hat{{\it\eta}}_{j}|}{|\text{}^{1,s}\hat{{\it\eta}}_{j}|+|\text{}^{1,v}\hat{{\it\eta}}_{j}|}.\end{eqnarray}$$
2.2.3. Second-order solutions
The second-order solutions are assumed as
$$\begin{eqnarray}\displaystyle (\text{}^{2}{\it\phi}_{l},\text{}^{2}{\it\phi}_{gj},\text{}^{2}{\it\eta}_{j}) & = & \displaystyle (\text{}^{21,ss}\hat{{\it\phi}}_{l},\text{}^{21,ss}\hat{{\it\phi}}_{gj},\text{}^{21,ss}\hat{{\it\eta}}_{j})\exp (2\text{i}kx-2\text{i}{\it\omega}_{1s}t)\nonumber\\ \displaystyle & & \displaystyle +\,(\text{}^{21,s\overline{s}}\hat{{\it\phi}}_{l},\text{}^{21,s\overline{s}}\hat{{\it\phi}}_{gj},\text{}^{21,s\overline{s}}\hat{{\it\eta}}_{j})\exp (-\text{i}{\it\omega}_{1s}t+\text{i}\overline{{\it\omega}}_{1s}t)\nonumber\\ \displaystyle & & \displaystyle +\,(\text{}^{21,vv}\hat{{\it\phi}}_{l},\text{}^{21,vv}\hat{{\it\phi}}_{gj},\text{}^{21,vv}\hat{{\it\eta}}_{j})\exp (2\text{i}kx-2\text{i}{\it\omega}_{1v}t)\nonumber\\ \displaystyle & & \displaystyle +\,(\text{}^{21,v\overline{v}}\hat{{\it\phi}}_{l},\text{}^{21,v\overline{v}}\hat{{\it\phi}}_{gj},\text{}^{21,v\overline{v}}\hat{{\it\eta}}_{j})\exp (-\text{i}{\it\omega}_{1v}t+\text{i}\overline{{\it\omega}}_{1v}t)\nonumber\\ \displaystyle & & \displaystyle +\,(\text{}^{21,sv}\hat{{\it\phi}}_{l},\text{}^{21,sv}\hat{{\it\phi}}_{gj},\text{}^{21,sv}\hat{{\it\eta}}_{j})\exp (2\text{i}kx-\text{i}{\it\omega}_{1s}t-\text{i}{\it\omega}_{1v}t)\nonumber\\ \displaystyle & & \displaystyle +\,(\text{}^{21,s\overline{v}}\hat{{\it\phi}}_{l},\text{}^{21,s\overline{v}}\hat{{\it\phi}}_{gj},\text{}^{21,s\overline{v}}\hat{{\it\eta}}_{j})\exp (-\text{i}{\it\omega}_{1s}t+\text{i}\overline{{\it\omega}}_{1v}t)\nonumber\\ \displaystyle & & \displaystyle +\,(\text{}^{22,s}\hat{{\it\phi}}_{l},\text{}^{22,s}\hat{{\it\phi}}_{gj},\text{}^{22,s}\hat{{\it\eta}}_{j})\exp (2\text{i}kx-\text{i}{\it\omega}_{2s}t)\nonumber\\ \displaystyle & & \displaystyle +\,(\text{}^{22,\overline{s}}\hat{{\it\phi}}_{l},\text{}^{22,\overline{s}}\hat{{\it\phi}}_{gj},\text{}^{22,\overline{s}}\hat{{\it\eta}}_{j})\exp (-\text{i}{\it\omega}_{2\overline{s}}t)\nonumber\\ \displaystyle & & \displaystyle +\,(\text{}^{22,v}\hat{{\it\phi}}_{l},\text{}^{22,v}\hat{{\it\phi}}_{gj},\text{}^{22,v}\hat{{\it\eta}}_{j})\exp (2\text{i}kx-\text{i}{\it\omega}_{2v}t)\nonumber\\ \displaystyle & & \displaystyle +\,(\text{}^{22,\overline{v}}\hat{{\it\phi}}_{l},\text{}^{22,\overline{v}}\hat{{\it\phi}}_{gj},\text{}^{22,\overline{v}}\hat{{\it\eta}}_{j})\exp (-\text{i}{\it\omega}_{2\overline{v}}t)+\text{c.c.}\end{eqnarray}$$
The modes with superscript ‘21’ have the same exponent coefficients as the inhomogeneous terms formed by the first-order products in (2.18)–(2.20), and they represent the transfer of perturbation from the first order to the second order. Note that
$\overline{{\it\omega}}_{1s}$
and
$\overline{{\it\omega}}_{1v}$
represent the complex conjugates of
${\it\omega}_{1s}$
and
${\it\omega}_{1v}$
, but
$\overline{s}$
and
$\overline{v}$
are merely symbols indicating that the disturbances are partly generated by the complex conjugates of first-order disturbances. The modes with superscript 22 have second-order frequencies:
${\it\omega}_{2s}$
and
${\it\omega}_{2\overline{s}}$
are the second-order sinuous frequencies, and
${\it\omega}_{2v}$
and
${\it\omega}_{2\overline{v}}$
are the second-order varicose frequencies. Thus they represent the second-order inherent disturbances. Their corresponding wavenumbers are chosen to be
$2k$
and 0 (identical to the ‘21’ modes) in order to satisfy the initial condition (2.15).
In (2.28), modes ‘
$21,ss$
’ and ‘
$21,s\overline{s}$
’ are generated from first-order sinuous waves, so they represent the first harmonics of the linear sinuous mode. Similarly modes ‘
$21,vv$
’ and ‘
$21,v\overline{v}$
’ represent the first harmonics of the linear varicose mode. Modes ‘
$21,sv$
’ and ‘
$21,s\overline{v}$
’ are the first harmonics generated by the first-order sinuous and varicose waves together; therefore, they represent the coupling between the linear sinuous and varicose modes in the second order. They are produced only when the linear sinuous and varicose waves are both present, so their existence signifies that the sinuous and varicose waves do not propagate independently when the nonlinearity is taken into account.
The solutions for ‘
$21$
’ modes are derived by solving inhomogeneous equations. Take the mode ‘
$21,ss$
’ as an example. Substituting (2.28) into the mass conservation equations (2.16) and (2.17), and collecting terms with coefficient
$\exp (2\text{i}kx-2\text{i}{\it\omega}_{1s}t)$
, one arrives at
$$\begin{eqnarray}\displaystyle & \text{}^{21,ss}\hat{{\it\phi}_{l}^{\prime \prime }}-4k^{2}(\text{}^{21,ss}\hat{{\it\phi}}_{l})=0\quad \text{for }-1\leqslant y\leqslant 1, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \text{}^{21,ss}\hat{{\it\phi}_{gj}^{\prime \prime }}-4k_{}^{2}(\text{}^{21,ss}\hat{{\it\phi}}_{gj})=0\quad j=1,2\text{ for }y\geqslant 1,y\leqslant -1. & \displaystyle\end{eqnarray}$$
$y\rightarrow \pm \infty$
, yields 
$$\begin{eqnarray}\displaystyle & \text{}^{21,ss}\hat{{\it\phi}}_{l}=A\exp (2ky)+B\exp (-2ky)\quad \text{for }-1\leqslant y\leqslant 1, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \text{}^{21,ss}\hat{{\it\phi}}_{gj}=C_{j}\exp [(-1)^{j}2ky]\quad j=1,2\text{ for }y\geqslant 1,y\leqslant -1, & \displaystyle\end{eqnarray}$$
$A$
,
$B$
and
$C_{j}~(j=1,2)$
are integral constants. Substituting (2.28) into boundary conditions (2.18)–(2.20), the coefficients of
$\exp (2\text{i}kx-2\text{i}{\it\omega}_{1s}t)$
produce 
$$\begin{eqnarray}\displaystyle & \text{}^{21,ss}\hat{{\it\phi}}_{l}^{\prime }+2\text{i}({\it\omega}_{1s}-k)\text{}^{21,ss}\hat{{\it\eta}}_{j}=-k^{2}(\text{}^{1,s}\hat{{\it\phi}}_{l})\text{}^{1,s}\hat{{\it\eta}}_{j}-\text{}^{1,s}\hat{{\it\eta}}_{j}\text{}^{1,s}\hat{{\it\phi}}_{l}^{\prime \prime }\quad \text{at }y=(-1)^{j+1}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \text{}^{21,ss}\hat{{\it\phi}}_{gj}^{\prime }+2\text{i}{\it\omega}_{1s}\text{}^{21,ss}\hat{{\it\eta}}_{j}=-k^{2}(\text{}^{1,s}\hat{{\it\phi}}_{gj})\text{}^{1,s}\hat{{\it\eta}}_{j}-\text{}^{1,s}\hat{{\it\eta}}_{j}\text{}^{1,s}{\hat{{\it\phi}}_{gj}}^{\prime \prime }\quad \text{at }y=(-1)^{j+1}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle 2\text{i}[-{\it\rho}{\it\omega}_{1s}\text{}^{21,ss}\hat{{\it\phi}}_{gj}+{\it\omega}_{1s}\text{}^{21,ss}\hat{{\it\phi}}_{l}-k(\text{}^{21,ss}\hat{{\it\phi}}_{l})]+(-1)^{j}\frac{4k^{2}(\text{}^{21,ss}\hat{{\it\eta}}_{j})}{\mathit{We}}\nonumber\\ \displaystyle & & \displaystyle \quad =\text{i}\text{}^{1,s}\hat{{\it\eta}}_{j}[{\it\rho}{\it\omega}_{1s}\text{}^{1,s}\hat{{\it\phi}}_{gj}^{\prime }-{\it\omega}_{1s}\text{}^{1,s}{\hat{{\it\phi}}^{\prime }}_{l}+k(\text{}^{1,s}{\hat{{\it\phi}}^{\prime }}_{l})]-\frac{1}{2}{\it\rho}[\text{}^{1,s}\hat{{\it\phi}}_{gj}^{\prime 2}-k^{2}(\text{}^{1,s}\hat{{\it\phi}}_{gj}^{2})]\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\frac{1}{2}[\text{}^{1,s}\hat{{\it\phi}}_{l}^{\prime 2}-k^{2}(\text{}^{1,s}\hat{{\it\phi}}_{l}^{2})]\quad \text{at }y=(-1)^{j+1}.\end{eqnarray}$$
Substituting (2.31) and (2.32) into (2.33)–(2.35), six equations in six unknown constants
$A$
,
$B$
,
$C_{j}~(j=1,2)$
and
$\text{}^{21,ss}\hat{{\it\eta}}_{j}~(j=1,2)$
are obtained.
The solvability of the inhomogeneous equations is discussed next according to the Fredholm alternative. In the second-order equations (2.16)–(2.20), the forms of the second-order disturbances, i.e. the linear parts on the left-hand side of the equations, are exactly the same as the forms of the first-order disturbances in the first-order equations (2.10)–(2.14), so the coefficient matrix of
$A$
,
$B$
,
$C_{j}~(j=1,2)$
and
$\text{}^{21,ss}\hat{{\it\eta}}_{j}~(j=1,2)$
, referred to as
$\unicode[STIX]{x1D648}(2{\it\omega}_{1s},2k)$
, is the same as the coefficient matrix of the corresponding constants in the first-order linear analysis,
$\unicode[STIX]{x1D648}({\it\omega}_{1s},k)$
, as long as
$2{\it\omega}_{1s}$
and
$2k$
are replaced by
${\it\omega}_{1s}$
and
$k$
, respectively. The linear analysis indicates that the singularity of
$\unicode[STIX]{x1D648}({\it\omega},k)$
, which guarantees the non-trivial solutions of the linear homogeneous equations, is satisfied when and only when the dispersion relation
$D_{s}({\it\omega},k)D_{v}({\it\omega},k)=0$
is satisfied. The property of
$D_{s}$
and
$D_{v}$
governed by (2.24) determines that
$D_{s}({\it\omega}_{1s},k)=0$
(the linear dispersion relation required by (2.23)) automatically results in
$D_{s}(2{\it\omega}_{1s},2k)\neq 0$
and
$D_{v}(2{\it\omega}_{1s},2k)\neq 0$
, and therefore
$\unicode[STIX]{x1D648}(2{\it\omega}_{1s},2k)$
is non-singular and thus invertible. So according to the first Fredholm alternative,
$A$
,
$B$
,
$C_{j}~(j=1,2)$
and
$\text{}^{21,ss}\hat{{\it\eta}}_{j}~(j=1,2)$
have one and only one solution.
The solvability of the other inhomogeneous modes can be examined by the same procedure. Similar to mode ‘
$21,ss$
’, the modes ‘
$21,vv$
’, ‘
$21,sv$
’ and ‘
$21,s\overline{v}$
’ are solvable in the whole wavenumber range
$k>0$
, but the modes ‘
$21,s\overline{s}$
’ and ‘
$21,v\overline{v}$
’ are solvable only at wavenumbers smaller than the cutoff wavenumber
$k_{c}$
(the wavenumber corresponding to zero temporal growth rate). This is because, when
$k>k_{c}$
,
${\it\omega}_{1s}$
and
${\it\omega}_{1v}$
become real and this makes their corresponding coefficient matrices become
$\unicode[STIX]{x1D648}^{\prime }(0,0)$
. (The form of the coefficient matrix for zero wavenumber modes
$\unicode[STIX]{x1D648}^{\prime }({\it\omega},0)$
is different from the non-zero wavenumber counterpart
$\unicode[STIX]{x1D648}({\it\omega},k)$
because the velocity potential is in the form of
$Ay+B$
rather than the exponential terms in (2.31) and (2.32).) The matrix
$\unicode[STIX]{x1D648}^{\prime }(0,0)$
is singular and not invertible, and the vector of the inhomogeneous terms generated by the first order, referred to as
$\boldsymbol{b}$
, does not guarantee the second Fredholm alternative, that is,
$\boldsymbol{b}$
is orthogonal to all eigenvectors of
$\unicode[STIX]{x1D648}^{\prime }(0,0)^{\text{T}}$
with zero eigenvalue, so the inhomogeneous equations are not solvable. It can be confirmed from (A 4) and (A 8) in appendix A that the solution of
$\text{}^{21,s\overline{s}}\hat{{\it\phi}}_{gj}$
and
$\text{}^{21,v\overline{v}}\hat{{\it\phi}}_{gj}$
does not exist when
${\it\omega}_{1s}$
and
${\it\omega}_{1v}$
are real, which happens at
$k>k_{c}$
. Therefore, for the present study, the solvability condition requires that the wavenumber
$k$
should be smaller than the cutoff wavenumber
$k_{c}$
. This is also desirable in the current analysis because, in this wavenumber range, the temporal growth rates
$\text{Im}({\it\omega}_{1s})$
and
$\text{Im}({\it\omega}_{1v})$
are positive, representing the instability of the flow field. In order to extend the wavenumber range to
$k>k_{c}$
, other measures, such as expanding the cutoff wavenumber
$k_{c}$
in power series of
${\it\eta}_{0}$
, may be required (Yuen Reference Yuen1968), but that is not the focus of the present work.
Back to the solution of mode ‘
$21,ss$
’. Knowing that the solution is unique, the following parity analysis is helpful. If the relations
$$\begin{eqnarray}\text{}^{21,ss}\hat{{\it\eta}}_{1}=-\text{}^{21,ss}\hat{{\it\eta}}_{2},\quad A=B,\quad C_{1}=C_{2},\end{eqnarray}$$
exist, then the boundary conditions (2.33)–(2.35) will be identical at
$y=1$
and
$y=-1$
, and solving only those at
$y=1$
is enough. Thus
$\text{}^{21,ss}\hat{{\it\eta}}_{1}$
is solved as
$$\begin{eqnarray}\text{}^{21,ss}\hat{{\it\eta}}_{1}=\frac{k[({\it\omega}_{1s}-k)^{2}(3-\tanh ^{2}k-4\tanh k\coth 2k)+2{\it\rho}{\it\omega}_{1s}^{2}]}{D_{v}(2{\it\omega}_{1s},2k)}\text{}^{1,s}\hat{{\it\eta}}_{1}^{2}.\end{eqnarray}$$
The solutions for
$\text{}^{21,ss}\hat{{\it\phi}}_{l}$
and
$\text{}^{21,ss}\hat{{\it\phi}}_{gj}$
are given in appendix A.
Using the same solution procedure for other ‘21’ modes, it is easy to deduce the relations between the upper and lower interface displacements via parity analysis:
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\text{}^{21,s\overline{s}}\hat{{\it\eta}}_{1}=-\text{}^{21,s\overline{s}}\hat{{\it\eta}}_{2},\quad \text{}^{21,vv}\hat{{\it\eta}}_{1}=-\text{}^{21,vv}\hat{{\it\eta}}_{2},\quad \text{}^{21,v\overline{v}}\hat{{\it\eta}}_{1}=-\text{}^{21,v\overline{v}}\hat{{\it\eta}}_{2},\\ \text{}^{21,sv}\hat{{\it\eta}}_{1}=\text{}^{21,sv}\hat{{\it\eta}}_{2},\quad \text{}^{21,s\overline{v}}\hat{{\it\eta}}_{1}=\text{}^{21,s\overline{v}}\hat{{\it\eta}}_{2}.\end{array}\!\!\right\} & & \displaystyle\end{eqnarray}$$
The solutions for the upper interface are
$$\begin{eqnarray}\displaystyle & \text{}^{21,s\overline{s}}\hat{{\it\eta}}_{1}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \text{}^{21,vv}\hat{{\it\eta}}_{1}=\frac{k[({\it\omega}_{1v}-k)^{2}(3-\coth ^{2}k-4\coth k\coth 2k)+2{\it\rho}{\it\omega}_{1v}^{2}]}{D_{v}(2{\it\omega}_{1v},2k)}\text{}^{1,v}\hat{{\it\eta}}_{1}^{2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \text{}^{21,v\overline{v}}\hat{{\it\eta}}_{1}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \text{}^{21,sv}\hat{{\it\eta}}_{1} & = & \displaystyle \frac{2k}{D_{s}({\it\omega}_{1s}+{\it\omega}_{1v},2k)}\{\!({\it\omega}_{1s}-k)^{2}+({\it\omega}_{1v}-k)^{2}-({\it\omega}_{1s}+{\it\omega}_{1v}-2k)\nonumber\\ \displaystyle & & \displaystyle \times \,[({\it\omega}_{1s}-k)\tanh k+({\it\omega}_{1v}-k)\coth k]\tanh 2k+2{\it\rho}{\it\omega}_{1s}{\it\omega}_{1v}\!\}\text{}^{1,s}\hat{{\it\eta}}_{1}\text{}^{1,v}\hat{{\it\eta}}_{1},\qquad\end{eqnarray}$$
$$\begin{eqnarray}\text{}^{21,s\overline{v}}\hat{{\it\eta}}_{1}=0.\end{eqnarray}$$
Their corresponding velocity potentials are given in appendix A. It can be seen that
$\text{}^{21,ss}\hat{{\it\eta}}_{j}$
and
$\text{}^{21,s\overline{s}}\hat{{\it\eta}}_{j}$
are identical to
$\hat{{\it\eta}}_{j21}$
and
$\hat{{\it\eta}}_{j22}$
respectively (the corresponding harmonics of linear sinuous mode for temporal instability) in the study of Yang et al. (Reference Yang, Wang, Fu, Du and Tong2013) when viscosity vanishes.
One important characteristic of the solutions is that the zero wavenumber interface displacements, i.e.
$\text{}^{21,s\overline{s}}\hat{{\it\eta}}_{j}$
,
$\text{}^{21,v\overline{v}}\hat{{\it\eta}}_{j}$
and
$\text{}^{21,s\overline{v}}\hat{{\it\eta}}_{j}$
, are zero, indicating that up to the second order there is no unstable mode that makes the sheet dilate or attenuate uniformly everywhere. This is a special property of the temporal instability of planar liquid sheets, which can also be proved by the study of Clark & Dombrowski (Reference Clark and Dombrowski1972), and it is very different from the nonlinear instability of cylindrical jets (Yuen Reference Yuen1968). Another important characteristic is that the first harmonic of the linear sinuous mode
$\text{}^{21,ss}\hat{{\it\eta}}_{j}$
and the first harmonic of the linear varicose mode
$\text{}^{21,vv}\hat{{\it\eta}}_{j}$
are varicose, but the first harmonic generated by the coupling between the linear sinuous and varicose modes,
$\text{}^{21,sv}\hat{{\it\eta}}_{j}$
, is sinuous.
Since
$\text{}^{21,ss}\hat{{\it\eta}}_{1}$
,
$\text{}^{21,vv}\hat{{\it\eta}}_{1}$
and
$\text{}^{21,sv}\hat{{\it\eta}}_{1}$
are respectively in proportion to
$\text{}^{1,s}\hat{{\it\eta}}_{1}^{2}$
,
$\text{}^{1,v}\hat{{\it\eta}}_{1}^{2}$
and
$\text{}^{1,s}\hat{{\it\eta}}_{1}\text{}^{1,v}\hat{{\it\eta}}_{1}$
, from which they are generated, it is useful define the ‘relative second-order disturbance amplitudes’ as
$$\begin{eqnarray}\text{}^{2r,ss}\hat{{\it\eta}}=\frac{\text{}^{21,ss}\hat{{\it\eta}}_{1}}{\text{}^{1,s}\hat{{\it\eta}}_{1}^{2}},\quad \text{}^{2r,vv}\hat{{\it\eta}}=\frac{\text{}^{21,vv}\hat{{\it\eta}}_{1}}{\text{}^{1,v}\hat{{\it\eta}}_{1}^{2}},\quad \text{}^{2r,sv}\hat{{\it\eta}}=\frac{\text{}^{21,sv}\hat{{\it\eta}}_{1}}{\text{}^{1,s}\hat{{\it\eta}}_{1}\text{}^{1,v}\hat{{\it\eta}}_{1}}.\end{eqnarray}$$
Here
$\text{}^{2r,ss}\hat{{\it\eta}}$
,
$\text{}^{2r,vv}\hat{{\it\eta}}$
and
$\text{}^{2r,sv}\hat{{\it\eta}}$
represent the ability of each linear mode to generate the corresponding nonlinear harmonic; their norms are a good indicator of the initial (at
$t=0,$
so the effect of temporal growth rate is excluded) extent of nonlinearity, especially when
$|\text{}^{1,s}\hat{{\it\eta}}_{j}|=|\text{}^{1,v}\hat{{\it\eta}}_{j}|$
.
The equations of modes ‘
$22$
’ do not contain the first-order inhomogeneous products, so they are solved using the same method as standard linear analysis. The second-order temporal frequencies are determined by the dispersion relations
$$\begin{eqnarray}D_{s}({\it\omega}_{2s},2k)=0,\quad D_{s}({\it\omega}_{2\overline{s}},0)=0,\quad D_{v}({\it\omega}_{2v},2k)=0,\quad D_{v}({\it\omega}_{2\overline{v}},0)=0.\end{eqnarray}$$
Similar to linear analysis, the sinuous and varicose disturbance amplitudes yield
$$\begin{eqnarray}\text{}^{22,s}\hat{{\it\eta}}_{1}=\text{}^{22,s}\hat{{\it\eta}}_{2},\quad \text{}^{22,\overline{s}}\hat{{\it\eta}}_{1}=\text{}^{22,\overline{s}}\hat{{\it\eta}}_{2},\quad \text{}^{22,v}\hat{{\it\eta}}_{1}=-\text{}^{22,v}\hat{{\it\eta}}_{2},\quad \text{}^{22,\overline{v}}\hat{{\it\eta}}_{1}=-\text{}^{22,\overline{v}}\hat{{\it\eta}}_{2}.\end{eqnarray}$$
Substituting (2.28), (2.36a
), (2.38) and (2.46) into the second-order initial condition (2.15),
$\text{}^{22,s}\hat{{\it\eta}}_{j}$
,
$\text{}^{22,\overline{s}}\hat{{\it\eta}}_{j}$
,
$\text{}^{22,v}\hat{{\it\eta}}_{j}$
and
$\text{}^{22,\overline{v}}\hat{{\it\eta}}_{j}$
are solved as
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\text{}^{22,s}\hat{{\it\eta}}_{j}=-\text{}^{21,sv}\hat{{\it\eta}}_{j},\quad \text{}^{22,\overline{s}}\hat{{\it\eta}}_{j}=-\text{}^{21,s\overline{v}}\hat{{\it\eta}}_{j}=0,\\ \text{}^{22,v}\hat{{\it\eta}}_{j}=-\text{}^{21,ss}\hat{{\it\eta}}_{j}-\text{}^{21,vv}\hat{{\it\eta}}_{j},\quad \text{}^{22,\overline{v}}\hat{{\it\eta}}_{j}=-\text{}^{21,s\overline{s}}\hat{{\it\eta}}_{j}-\text{}^{21,v\overline{v}}\hat{{\it\eta}}_{j}=0.\end{array}\!\!\right\} & & \displaystyle\end{eqnarray}$$
The solutions of the velocity potentials of the ‘22’ modes are given in appendix A.
2.2.4. Summary
Removing the modes whose interface displacements are zero and neglecting the remainder term
$O({\it\eta}_{0}^{3})$
, the second-order approximate interface displacement is summarized as
$$\begin{eqnarray}\displaystyle {\it\eta}_{j}=\text{}^{1,s}{\it\eta}_{j}+\text{}^{1,v}{\it\eta}_{j}+\text{}^{21,ss}{\it\eta}_{j}+\text{}^{21,vv}{\it\eta}_{j}+\text{}^{21,sv}{\it\eta}_{j}+\text{}^{22,s}{\it\eta}_{j}+\text{}^{22,v}{\it\eta}_{j}, & & \displaystyle\end{eqnarray}$$
with
$$\begin{eqnarray}\displaystyle \text{}^{1,s}{\it\eta}_{j} & = & \displaystyle \text{}^{1,s}\hat{{\it\eta}}_{j}{\it\eta}_{0}\exp (\text{i}kx-\text{i}{\it\omega}_{1s}t)+\text{c.c.},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \text{}^{1,v}{\it\eta}_{j} & = & \displaystyle \text{}^{1,v}\hat{{\it\eta}}_{j}{\it\eta}_{0}\exp (\text{i}kx-\text{i}{\it\omega}_{1v}t)+\text{c.c.},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \text{}^{21,ss}{\it\eta}_{j} & = & \displaystyle \text{}^{21,ss}\hat{{\it\eta}}_{j}{\it\eta}_{0}^{2}\exp (2\text{i}kx-2\text{i}{\it\omega}_{1s}t)+\text{c.c.},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \text{}^{21,vv}{\it\eta}_{j} & = & \displaystyle \text{}^{21,vv}\hat{{\it\eta}}_{j}{\it\eta}_{0}^{2}\exp (2\text{i}kx-2\text{i}{\it\omega}_{1v}t)+\text{c.c.},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \text{}^{21,sv}{\it\eta}_{j} & = & \displaystyle \text{}^{21,sv}\hat{{\it\eta}}_{j}{\it\eta}_{0}^{2}\exp (2\text{i}kx-\text{i}{\it\omega}_{1s}t-\text{i}{\it\omega}_{1v}t)+\text{c.c.,}\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \text{}^{22,s}{\it\eta}_{j} & = & \displaystyle \text{}^{22,s}\hat{{\it\eta}}_{j}{\it\eta}_{0}^{2}\exp (2\text{i}kx-\text{i}{\it\omega}_{2s}t)+\text{c.c.},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \text{}^{22,v}{\it\eta}_{j} & = & \displaystyle \text{}^{22,v}\hat{{\it\eta}}_{j}{\it\eta}_{0}^{2}\exp (2\text{i}kx-\text{i}{\it\omega}_{2v}t)+\text{c.c.}\end{eqnarray}$$
2.3. Numerical simulation
The numerical simulation is performed using a boundary integral method; it gives the solution of the full problem and verifies the validity of the weakly nonlinear analysis.
In the inviscid potential flow considered in this paper, the boundary integral method represents the flow field by two vortex sheets located at the upper and lower interfaces. (Note that the ‘sheets’ in the ‘vortex sheets’ are different from those in the ‘liquid sheets’; the ‘vortex sheet’ here represents the velocity discontinuity between inviscid gas and inviscid liquid on the gas-to-liquid interface.) The motions of these two vortex sheets are induced by the vorticity on them. The vortex sheet strength changes with time, owing to the baroclinic generation of vorticity.
Because this boundary integral method requires that vorticity is confined to the two interfaces, it is only capable of potential flows, and this is the very reason why both the liquid and the gas must be assumed inviscid. However, the sacrifice of viscosity brings about the benefit of accuracy: since any vortex on the interface induces a continuous and irrotational velocity field, the governing equations of internal and external flow, (2.2) and (2.3), are always precisely satisfied, though the discretization of the vortex sheet can cause numerical error in the boundary conditions. As a result, this boundary integral method can give very precise numerical solutions of interface evolution. It is noted that, apart from the inviscid potential flow in this paper, other types of boundary integral methods are applied to other problems: for example, in the context of high viscosity and small disturbance amplitude, respectively, the nonlinear Navier–Stokes equation is linearized into the Stokes equation and the unsteady Stokes equation, which make it possible to derive boundary integral solutions via the Green’s function. These two types of boundary integral methods are very different from the one in this paper; a comprehensive introduction to them can be found in the textbook of Pozrikidis (Reference Pozrikidis1992).
The following is the formulation of the boundary integral simulation. We represent the upper and lower interfacial curves by complex variables parametrically in
${\it\xi}$
,
$$\begin{eqnarray}z_{1}({\it\xi},t)=x_{1}({\it\xi},t)+\text{i}y_{1}({\it\xi},t),\quad z_{2}({\it\xi},t)=x_{2}({\it\xi},t)+\text{i}y_{2}({\it\xi},t),\end{eqnarray}$$
where the subscript ‘1’ represents the upper interface and ‘2’ represents the lower interface. Equations (2.50) are similar to the interface denoted by Baker & Beale (Reference Baker and Beale2004), although they only considered the motion of one interface.
The arclengths along the interfaces are
$s_{1}({\it\xi},t)$
and
$s_{2}({\it\xi},t)$
. Similar to the study of Hou, Lowengrub & Shelley (Reference Hou, Lowengrub and Shelley1994), we use
${\it\gamma}_{1}({\it\xi},t)$
and
${\it\gamma}_{2}({\it\xi},t)$
to denote the ‘unnormalized vortex sheet strengths’, and thus
${\it\gamma}_{1}({\it\xi},t)/s_{1,{\it\xi}}({\it\xi},t)$
and
${\it\gamma}_{2}({\it\xi},t)/s_{2,{\it\xi}}({\it\xi},t)$
are the corresponding ‘true vortex sheet strengths’, i.e. the tangential jumps in velocity across the interfaces. As temporal instability is considered, both interfaces are
$2{\rm\pi}$
-periodic in
${\it\xi}$
, that is,
$x_{j}({\it\xi}+2{\rm\pi})=2{\rm\pi}/k+x_{j}({\it\xi})$
,
$y_{j}({\it\xi}+2{\rm\pi})=y_{j}({\it\xi})$
and
${\it\gamma}_{j}({\it\xi}+2{\rm\pi})={\it\gamma}_{j}({\it\xi})$
(
$j=1,2$
).
The motion of the interfaces is governed by the Birkhoff–Rott equation, which indicates that the velocity of the fluid is induced by all the vorticity on the interfaces:
$$\begin{eqnarray}\displaystyle \frac{\partial \bar{z}_{1}}{\partial t}({\it\xi},t) & = & \displaystyle \frac{k}{4{\rm\pi}\text{i}}\text{P.V.}\int _{0}^{2{\rm\pi}}{\it\gamma}_{1}({\it\xi}^{\prime },t)\cot \left(k\frac{z_{1}({\it\xi},t)-z_{1}({\it\xi}^{\prime },t)}{2}\right)\text{d}{\it\xi}^{\prime }\nonumber\\ \displaystyle & & \displaystyle +\,\frac{k}{4{\rm\pi}\text{i}}\int _{0}^{2{\rm\pi}}{\it\gamma}_{2}({\it\xi}^{\prime },t)\cot \left(k\frac{z_{1}({\it\xi},t)-z_{2}({\it\xi}^{\prime },t)}{2}\right)\text{d}{\it\xi}^{\prime },\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \frac{\partial \bar{z}_{2}}{\partial t}({\it\xi},t) & = & \displaystyle \frac{k}{4{\rm\pi}\text{i}}\int _{0}^{2{\rm\pi}}{\it\gamma}_{1}({\it\xi}^{\prime },t)\cot \left(k\frac{z_{2}({\it\xi},t)-z_{1}({\it\xi}^{\prime },t)}{2}\right)\text{d}{\it\xi}^{\prime }\nonumber\\ \displaystyle & & \displaystyle +\,\frac{k}{4{\rm\pi}\text{i}}\text{P.V.}\int _{0}^{2{\rm\pi}}{\it\gamma}_{2}({\it\xi}^{\prime },t)\cot \left(k\frac{z_{2}({\it\xi},t)-z_{2}({\it\xi}^{\prime },t)}{2}\right)\text{d}{\it\xi}^{\prime },\end{eqnarray}$$
$k$
is included in these equations, because we do not force the wavelength to be
$2{\rm\pi}$
, while Baker & Beale (Reference Baker and Beale2004) did. It should be noted that the equations in Baker & Beale (Reference Baker and Beale2004) have additional terms involving a weighting parameter
${\it\alpha}$
. Here we choose that parameter to be zero, i.e. the Lagrangian formulation (Hou, Lowengrub & Shelley Reference Hou, Lowengrub and Shelley1997), in order to simplify the equations.
The change of vortex sheet strengths is caused by the baroclinicity on the gas-to-liquid interface, governed by the following equation:
$$\begin{eqnarray}\frac{\partial {\it\gamma}_{j}}{\partial t}=-2\mathit{At}_{j}\left[\text{Re}\left\{z_{j,{\it\xi}}\frac{\partial ^{2}\bar{z}_{j}}{\partial t^{2}}\right\}+\frac{1}{8}\left(\frac{{\it\gamma}_{j}^{2}}{|z_{j,{\it\xi}}|^{2}}\right)_{\!{\it\xi}}\right]+\frac{2}{\mathit{We}({\it\rho}+1)}{\it\kappa}_{j,{\it\xi}},\quad j=1,2,\end{eqnarray}$$
where
${\it\kappa}_{j}$
is the curvature of the interfacial curve, and
$\mathit{At}_{j}$
is the Atwood number,
$$\begin{eqnarray}\mathit{At}_{j}=(-1)^{j}\frac{{\it\rho}-1}{{\it\rho}+1}.\end{eqnarray}$$
Equation (2.53) is a modification of equation (2.5a) of Baker & Beale (Reference Baker and Beale2004) to include the effect of surface tension, reflected in the term that contains
$\mathit{We}$
. Equation (10) of Hou et al. (Reference Hou, Lowengrub and Shelley1997) contains a similar but slightly different term, as the definition of the Weber number in that paper is different from that used here.
The governing equations (2.51), (2.52) and (2.53) are integrated explicitly in time by the fourth-order Runge–Kutta method. As pointed out by Hou et al. (Reference Hou, Lowengrub and Shelley1994), boundary integral methods for surface tension flows suffer from a severe time-step stability constraint when the grid spacing is small, and they developed a small-scale decomposition technique to break this constraint. However, our numerical experiments indicate that, for the problems here, where the gas-to-liquid density ratio is of the order of 0.1, the two surfaces of the liquid sheets always contact before any small-scale features develop on the fluid interfaces. Therefore, it is not necessary to apply the technique of Hou et al. (Reference Hou, Lowengrub and Shelley1994).
For the spatial discretization, discrete markers on the interfaces are used: in (2.50),
${\it\xi}$
is given with discrete values of
$$\begin{eqnarray}{\it\xi}=2{\rm\pi}n/N_{j},\quad n=0,1,2,\dots ,N_{j}-1,j=1,2,\end{eqnarray}$$
where
$N_{j}$
is the number of markers within one wavelength on the upper interface (for
$j=1$
) and lower interface (for
$j=2$
).
As shown in Beale, Hou & Lowengrub (Reference Beale, Hou and Lowengrub1996), with regard to boundary integral methods for interfacial flows, special care must be taken to avoid numerical instability, and one way is to use the pseudospectral method along with a dealiasing filter. Therefore, in our numerical method, the following measures are taken. First, the principal value integrals in (2.51) and (2.52) are evaluated by the alternating trapezoidal rule, while the integrals other than the principal value integral are evaluated by the normal trapezoidal rule. This ensures that all integrals in these two equations are evaluated in spectral accuracy. Then, all derivatives with respect to
${\it\xi}$
in (2.53) are evaluated by the pseudospectral derivative operator. Finally, we adopt the 25th-order Fourier filter in Hou et al. (Reference Hou, Lowengrub and Shelley1994) to control aliasing errors. Details of these procedures can be found in the study of Hou et al. (Reference Hou, Lowengrub and Shelley1994) and Beale et al. (Reference Beale, Hou and Lowengrub1996).
A problem associated with the Lagrangian formulation is the accumulation and separation of interfacial markers during the course of simulation (Rangel & Sirignano Reference Rangel and Sirignano1991; Hou et al.
Reference Hou, Lowengrub and Shelley1997). We circumvent this problem by a redistribution technique similar to the method of Rangel & Sirignano (Reference Rangel and Sirignano1991). For each interface, when the maximum separation of markers exceeds 120 % of the initial separation or when the minimum separation drops below 80 % of the initial separation, that interface is rediscretized to make the computational points equally spaced in arclength. In contrast to the approach of Rangel & Sirignano (Reference Rangel and Sirignano1991), which used linear interpolation to rediscretize the interface, we use trigonometric interpolation, which gives higher accuracy. During the rediscretization operation, the number of markers
$N_{j}$
may increase, since both interfaces are usually elongated over time.
The initial conditions for the numerical simulation are obtained by evaluating the analytical solution at
$t=0$
as follows:
$$\begin{eqnarray}\displaystyle & x_{j}({\it\xi},0)={\it\xi}/k,\quad j=1,2, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & y_{j}({\it\xi},0)=(-1)^{j+1}+{{\it\eta}_{j}|}_{t=0,x=x_{j}({\it\xi},0)},\quad j=1,2, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle {\it\gamma}_{j}({\it\xi},0) & = & \displaystyle (-1)^{j}x_{j,{\it\xi}}({\it\xi},0)({\it\phi}_{gj,x}-{\it\phi}_{l,x})|_{t=0,x=x_{j}({\it\xi},0),y=y_{j}({\it\xi},0)}\nonumber\\ \displaystyle & & \displaystyle +\,(-1)^{j}y_{j,{\it\xi}}({\it\xi},0)({\it\phi}_{gj,y}-{\it\phi}_{l,y})|_{t=0,x=x_{j}({\it\xi},0),y=y_{j}({\it\xi},0)},\quad j=1,2.\end{eqnarray}$$
Here the initial interface displacement
${{\it\eta}_{j}|}_{t=0}$
is determined by the initial condition (2.1), and the velocity potentials
${\it\phi}_{l}$
and
${\it\phi}_{gj}$
are given according to the second-order approximate solution of weakly nonlinear analysis:
$$\begin{eqnarray}({\it\phi}_{l},{\it\phi}_{gj})=(x,0)+(\text{}^{1}{\it\phi}_{l},\text{}^{1}{\it\phi}_{gj}){\it\eta}_{0}+(\text{}^{2}{\it\phi}_{l},\text{}^{2}{\it\phi}_{gj}){\it\eta}_{0}^{2}.\end{eqnarray}$$
Discrete values of
${\it\xi}$
are given according to (2.55) with
$N_{1}=N_{2}=N$
, where
$N$
is the initial resolution. Hence, (2.56) makes the initial
$x$
-coordinates of the markers evenly distributed within one wavelength. Note that (2.58) is based on the fact that the vortex sheet strength is equal to the jump of tangential velocity across the interface.
Although the vortex sheet strengths are derived according to the solution of weakly nonlinear analysis via (2.58) and (2.59), the initial velocity field induced by the vortex sheets is not necessarily completely identical to that of the weakly nonlinear analysis. The main reason for this discrepancy is that the second-order weakly nonlinear solution is only an approximate solution; the kinematic boundary condition that the velocity normal to the interface should be continuous is not strictly satisfied. Since this boundary condition is certainly precisely guaranteed by the vortex sheets, disparities of the initial condition emerge. In order to mitigate this discrepancy, the initial velocity field of weakly nonlinear analysis must be very precise. The measure we undertake here is to choose small values of initial disturbance amplitude:
${\it\eta}_{0}\leqslant 0.2$
at
$|\hat{{\it\eta}}_{1}|=|\hat{{\it\eta}}_{2}|=0.5$
. But sometimes it may not be enough: in order that the initial second-order disturbances are very small, the chosen condition should also guarantee the demand that the relative second-order amplitudes given in (2.44) should not be too big; an example regarding this issue can be found in § 3.2.1.
To verify the accuracy of our numerical method, we have conducted a spatial convergence test, with three different initial resolutions,
$N=128$
,
$256$
and
$512$
, while the time steps are all
${\rm\Delta}t=1\times 10^{-3}$
. The condition of the test is
$\mathit{We}=5$
,
${\it\rho}=0.1$
,
$k=0.25$
and
${\it\eta}_{0}=0.2$
, and the initial disturbance is given according to (3.1) and (3.2), which appear subsequently with initial phase difference
${\it\theta}={\rm\pi}/2$
; the temporal evolution of the sheet is presented in figure 4(e–h). We measure the error by comparing the two lower-resolution calculations to the
$N=512$
calculation:
$$\begin{eqnarray}\displaystyle {\it\epsilon}_{N}(t) & = & \displaystyle \max \left[\vphantom{\frac{a}{b}}\!\right.\max _{0\leqslant {\it\xi}<2{\rm\pi}}|x_{1}^{N}({\it\xi},t)-x_{1}^{512}({\it\xi},t)|,~\max _{0\leqslant {\it\xi}<2{\rm\pi}}|x_{2}^{N}({\it\xi},t)-x_{2}^{512}({\it\xi},t)|,\nonumber\\ \displaystyle & & \displaystyle \qquad \max _{0\leqslant {\it\xi}<2{\rm\pi}}|y_{1}^{N}({\it\xi},t)-y_{1}^{512}({\it\xi},t)|,~\max _{0\leqslant {\it\xi}<2{\rm\pi}}|y_{2}^{N}({\it\xi},t)-y_{2}^{512}({\it\xi},t)|\!\left.\vphantom{\frac{a}{b}}\right]\!.\end{eqnarray}$$
The error is plotted in a negative logarithm scale in figure 3. For most of the time, the error is less than
$10^{-8}$
for both the
$N=128$
and
$N=256$
calculations, so the accuracy of our numerical simulation is very high. The
$N=128$
calculation exhibits a rapid loss of accuracy after approximately
$t=45$
, while the
$N=256$
case exhibits a similar loss after approximately
$t=55$
. By examining the sequence of interface positions shown in figure 4(e–h), it can be concluded that loss of accuracy occurs when the sheet is close to breakup. This is consistent with the conclusion of Hou et al. (Reference Hou, Lowengrub and Shelley1994), who found that the loss of accuracy occurs when the interfaces nearly pinch. But even at
$t=57.5$
, when the sheet tends to breakup, the error is still very small, smaller than
$10^{-1}$
and
$10^{-3}$
for
$N=128$
and
$N=256$
, respectively, and such accuracies are enough to analyse the breakup property of the sheet.
For all the numerical simulations in § 3, initially we use a resolution of
$N=128$
, which gives a higher calculation speed; when the sheet is close to breakup, we switch to
$N=256$
and
$N=512$
successively in order to obtain a higher accuracy.
Figure 3. Error plot:
$\mathit{We}=5$
,
${\it\rho}=0.1$
,
$k=0.25$
,
${\it\eta}_{0}=0.2$
,
${\rm\Delta}t=1\times 10^{-3}$
. The initial disturbance is (3.2) with
${\it\theta}={\rm\pi}/2$
.
Figure 4. Temporal evolution of liquid sheets at
${\it\theta}={\rm\pi}/2$
: (a–d) the results of the weakly nonlinear analysis and (e–h) the results of numerical simulation; (a,e)
$t=20$
, (b, f)
$t=45$
, (c,g)
$t=53$
and (d,h)
$t=57.5$
.
3. Results and discussion
3.1. Instability property and breakup mechanism
The temporal instability of sheets under a typical initial disturbance is studied: in (2.1),
$$\begin{eqnarray}\hat{{\it\eta}}_{1}={\textstyle \frac{1}{2}},\quad \hat{{\it\eta}}_{2}={\textstyle \frac{1}{2}}\exp (\text{i}{\it\theta})\end{eqnarray}$$
are set, and thus the initial disturbances for the upper and lower interfaces are
$$\begin{eqnarray}{{\it\eta}_{1}|}_{t=0}={\it\eta}_{0}\cos (kx),\quad {{\it\eta}_{2}|}_{t=0}={\it\eta}_{0}\cos (kx+{\it\theta}).\end{eqnarray}$$
Here
${\it\theta}$
represents the initial phase difference between the two interfaces. As
${\it\theta}$
changes from 0 to
${\rm\pi}$
, the initial disturbance gradually varies from sinuous to varicose, and the initial sinuous proportion
$r_{s}$
alters from 1 to 0. The instability behaviour of sheets for various values of
${\it\theta}$
is studied in this section, with an emphasis on the case of
${\it\theta}={\rm\pi}/2$
,
$r_{s}=0.5$
, in which the linear sinuous and varicose modes have equal initial amplitudes. The parameters in this section are given as follows: the condition is
$\mathit{We}=5$
,
${\it\rho}=0.1$
,
${\it\eta}_{0}=0.2$
,
$k=0.25$
; the complex temporal frequencies are
${\it\omega}_{1s}=0.178+0.0617\text{i}$
,
${\it\omega}_{1v}=0.244+0.0267\text{i}$
,
${\it\omega}_{2s}=0.5$
,
${\it\omega}_{2v}=0.5$
; and the relative second-order disturbance amplitudes (defined in (2.44)) are
$\text{}^{2r,ss}\hat{{\it\eta}}=-0.0103-0.0247\text{i}$
,
$\text{}^{2r,vv}\hat{{\it\eta}}=1.19+1.02\text{i}$
,
$\text{}^{2r,sv}\hat{{\it\eta}}=0.704+0.136\text{i}$
.
The comparison between the weakly nonlinear analysis and numerical simulation is discussed first. Figure 4 gives an example of sheet’s surface evolution (
${\it\theta}={\rm\pi}/2$
), which finally leads to breakup at
$t=57.5$
. Up till the time
$t=45$
, the results of the two methods are nearly the same, and this indicates that the weakly nonlinear analysis can precisely predict the shape of the sheet for most of the time of disturbance evolution. This also holds for other initial phase differences. Small discrepancies between the two methods emerge at
$t=53$
, when distinct thinning on the sheet takes place just before breakup. The shapes of the sheets when they nearly break up at various initial phase differences are given in figure 5. There are some subtle differences between the two methods: the wave amplitudes predicted by the weakly nonlinear analysis are slightly larger for most conditions; the surface from the numerical simulation is smoother in most regions; and the numerical simulation predicts a protuberance near the point of breakup. The main reason for these differences is that, because the
$n$
th-order disturbance has growth rates of
$n\,\text{Im}({\it\omega}_{1s})$
and
$n\,\text{Im}({\it\omega}_{1v})$
, the remainder term
$O({\it\eta}_{0}^{3})$
, which contains higher-order harmonics, becomes non-negligible at large values of
$t$
. Particularly, the small-scale protuberance near the location of breakup is an indicator of strong local nonlinearity, which is beyond the competence of the weakly nonlinear analysis. Despite these differences, however, the weakly nonlinear analysis still provides representative characteristics of the wave profile. The general shapes of the two methods are very similar, and, more importantly, except for the case of
${\it\theta}={\rm\pi}$
, the weakly nonlinear analysis gives very good prediction of the breakup time and breakup location, which are most important when studying the breakup property. Therefore, the qualitative agreement with numerical simulation justifies the weakly nonlinear analysis in explaining the mechanism of instability.
Figure 5. Shapes of sheets when they nearly break up at various initial phase differences: (a,b)
${\it\theta}=0$
,
$r_{s}=1$
, (c,d)
${\it\theta}={\rm\pi}/4$
,
$r_{s}=0.707$
, (e, f)
${\it\theta}={\rm\pi}/2$
,
$r_{s}=0.5$
, (g,h)
${\it\theta}=3{\rm\pi}/4$
,
$r_{s}=0.293$
, (i, j)
${\it\theta}={\rm\pi}$
,
$r_{s}=0$
; (a)
$t=61$
, (b)
$t=65$
, (c,d)
$t=59$
, (e, f)
$t=57.5$
, (g,h)
$t=57$
, (i)
$t=55$
, (j)
$t=63$
. (a,c,e,g,i) are the results of weakly nonlinear analysis, and (b,d,f,h,j) are the results of numerical simulation.
Figure 6. Temporal evolution of amplitudes of various modes: (a)
${\it\theta}=0$
,
$r_{s}=1$
; (b)
${\it\theta}={\rm\pi}/4$
,
$r_{s}=0.707$
; (c)
${\it\theta}={\rm\pi}/2$
,
$r_{s}=0.5$
; (d)
${\it\theta}=3{\rm\pi}/4$
,
$r_{s}=0.293$
; (e)
${\it\theta}={\rm\pi}$
,
$r_{s}=0$
.
Figure 6 shows the temporal evolution of the amplitudes of various unstable modes in weakly nonlinear analysis, which explains the instability behaviour in figure 5. The modes ‘
$22,s$
’ and ‘
$22,v$
’ are not plotted because their zero temporal growth rates make their amplitudes negligible when sheets break up. The linear sinuous mode ‘
$1,s$
’ has a much larger temporal growth rate than the linear varicose mode ‘
$1,v$
’, and has much larger initial amplitudes (amplitudes at
$t=0$
) than the second-order nonlinear modes, so its amplitude is predominant among all the other modes except for the case of
${\it\theta}={\rm\pi}$
,
$r_{s}=0$
, where it does not exist. As a result, in figure 5(a–h), the wave profiles are dominated by sinuous characteristics, though the decrease of initial sinuous proportion
$r_{s}$
from 1 to 0.293 leads to a significant reduction of wave amplitude. This also explains the phenomenon that, in most experiments, the sheets exhibit sinuous properties. The temporal growth rate of the linear varicose mode ‘
$1,v$
’ is much smaller, so, although its amplitude is prominent among various modes initially at
$t=0$
(except for
${\it\theta}=0$
,
$r_{s}=1$
), when the sheet breaks up at large values of
$t$
, its amplitude does not predominate as the linear sinuous mode does. However, the linear varicose mode’s role in the breakup is vital. Among the modes ‘
$1,v$
’, ‘
$21,ss$
’ and ‘
$21,vv$
’, which are varicose in nature and thus can contribute to breakup, the amplitude of the linear varicose mode ‘
$1,v$
’ is significant except at
${\it\theta}=0$
,
$r_{s}=1$
. Moreover, its ‘
$k$
’ wavenumber (in contrast to the ‘
$2k$
’ wavenumber of harmonics ‘
$21,ss$
’ and ‘
$21,vv$
’) makes the ligaments in figure 5(c–j) interspaced by one wavelength, while in figure 5(a,b), where the breakup is caused by the harmonic ‘
$21,ss$
’, the ligament is interspaced by one half-wavelength. The harmonics ‘
$21,ss$
’, ‘
$21,vv$
’ and ‘
$21,sv$
’ have small initial amplitudes, but, because they also have large temporal growth rates, their amplitudes become pronounced near breakup, indicating that the nonlinearity becomes important. In particular, for figure 6(b–d), where both linear modes exist, it is the mode ‘
$21,sv$
’ that finally has the largest amplitude among the nonlinear harmonics, and this indicates that the nonlinear coupling between the linear sinuous and varicose modes is significant. The first harmonics of the linear sinuous mode and linear varicose mode separately, ‘
$21,ss$
’ and ‘
$21,vv$
’, are varicose, and thus they contribute to the breakup process. On the contrary, the first harmonic generated by the coupling between the linear sinuous and varicose modes, ‘
$21,sv$
’, is sinuous, and it modulates the wave profile. The detailed functions of these nonlinear harmonics will be discussed subsequently.
In order to further probe the mechanism of instability and breakup, a Fourier analysis is performed. For the periodic wave with wavelength
$2{\rm\pi}/k$
, the interface displacement can be expressed by a Fourier series
$$\begin{eqnarray}{\it\eta}_{j}=\mathop{\sum }_{n=0}^{+\infty }\,\text{}^{nk}{\it\eta}_{j}=\mathop{\sum }_{n=0}^{+\infty }\,\text{}^{nk}\hat{{\it\eta}}_{j}\exp (\text{i}nkx)+\text{c.c.},\end{eqnarray}$$
with
$$\begin{eqnarray}\text{}^{nk}\hat{{\it\eta}}_{j}=\left\{\begin{array}{@{}l@{}}\displaystyle \frac{k}{4{\rm\pi}}\int _{0}^{2{\rm\pi}/k}{\it\eta}_{j}\,\text{d}x,\quad n=0,\\ \displaystyle \frac{k}{2{\rm\pi}}\int _{0}^{2{\rm\pi}/k}{\it\eta}_{j}\exp (-\text{i}nkx)\,\text{d}x,\quad n\neq 0,\end{array}\right.\end{eqnarray}$$
where
$\text{}^{nk}{\it\eta}_{j}$
is the displacement component with wavenumber
$nk$
and
$\text{}^{nk}\hat{{\it\eta}}_{j}$
is its corresponding complex amplitude. In the context of liquid sheets,
$\text{}^{nk}{\it\eta}_{j}$
is further analysed into a sinuous displacement
$\text{}^{nk,s}{\it\eta}_{j}$
and a varicose displacement
$\text{}^{nk,v}{\it\eta}_{j}$
:
$$\begin{eqnarray}\text{}^{nk}{\it\eta}_{j}=\text{}^{nk,s}{\it\eta}_{j}+\text{}^{nk,v}{\it\eta}_{j},\end{eqnarray}$$
with
$$\begin{eqnarray}\text{}^{nk,s}{\it\eta}_{j}={\textstyle \frac{1}{2}}(\text{}^{nk}{\it\eta}_{1}+\text{}^{nk}{\it\eta}_{2}),\quad \text{}^{nk,v}{\it\eta}_{j}=(-1)^{j+1}{\textstyle \frac{1}{2}}(\text{}^{nk}{\it\eta}_{1}-\text{}^{nk}{\it\eta}_{2}).\end{eqnarray}$$
The components of the weakly nonlinear analysis have more explicit expressions. According to (2.48) and (2.49), one arrives at
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\text{}^{0k,s}{\it\eta}_{j}=0,\quad \text{}^{0k,v}{\it\eta}_{j}=0,\quad \text{}^{1k,s}{\it\eta}_{j}=\text{}^{1,s}{\it\eta}_{j},\quad \text{}^{1k,v}{\it\eta}_{j}=\text{}^{1,v}{\it\eta}_{j},\\ \text{}^{2k,s}{\it\eta}_{j}=\text{}^{21,sv}{\it\eta}_{j}+\text{}^{22,s}{\it\eta}_{j},\quad \text{}^{2k,v}{\it\eta}_{j}=\text{}^{21,ss}{\it\eta}_{j}+\text{}^{21,vv}{\it\eta}_{j}+\text{}^{22,v}{\it\eta}_{j}.\end{array}\!\!\right\} & & \displaystyle\end{eqnarray}$$
According to (3.7), there are no zero wavenumber components in the second-order weakly nonlinear analysis. The numerical simulation agrees well with the weakly nonlinear analysis on this issue. Our results of numerical simulation indicate that the amplitudes of the ‘
$0k,s$
’ and ‘
$0k,v$
’ components are of the order of
$10^{-6}$
, smaller than the possible numerical error limit. Whether the precise solution should be exactly zero may require a further analysis, but that is beyond the scope of this paper. Because the remainder term
$O({\it\eta}_{0}^{3})$
also contains ‘
$1k$
’ and ‘
$2k$
’ components (albeit small), the displacements of ‘
$1k$
’ and ‘
$2k$
’ obtained from the two methods are not exactly the same, but their qualitative agreement presented subsequently is adequate and satisfactory to study the mechanism of instability.
Figure 7. Interface displacements of various Fourier components for
${\it\theta}={\rm\pi}/2$
,
$r_{s}=0.5$
at
$t=57.5$
: (a,b) shape of sheet, (c,d) ‘
$1k,s$
’ component, (e, f) ‘
$2k,s$
’ component, (g,h) ‘
$1k,v$
’ component, (i, j) ‘
$2k,v$
’ component; (a,c,e,g,i) weakly nonlinear analysis, and (b,d, f,h, j) numerical simulation.
Figure 7 shows the displacements of ‘
$1k$
’ and ‘
$2k$
’ components, which explain the instability mechanism more clearly. Only the case of
${\it\theta}={\rm\pi}/2$
,
$r_{s}=0.5$
is discussed, because in this condition the linear sinuous and varicose modes have equal initial amplitudes, which is more likely to happen in practical situations and serves as a good foundation for the comparison between the effects of various modes. The sinuous components have major influences on the wave profile. The ‘
$1k,s$
’ component, which represents the linear sinuous mode, has a much larger wave amplitude than all the other components, and this determines the general sinuous shape of the sheet. The ‘
$2k,s$
’ component modulates the wave profile: its peak near
$x=35$
shifts the peak of ‘
$1k,s$
’ leftwards; its upgrade between
$x=29$
and
$x=35$
strengthens the upgrade of ‘
$1k,s$
’ while its upgrade between
$x=41$
and
$x=47$
weakens the downgrade of ‘
$1k,s$
’, and this is the reason why in figure 7(a,b), the slope on the left side of the peak is steeper than the slope on the right. Since the ‘
$22,s$
’ mode is negligible near breakup, according to (3.7), the ‘
$2k,s$
’ component is mainly determined by mode ‘
$21,sv$
’, the mode generated by linear sinuous and varicose modes together, so this demonstrates the significance of the nonlinear coupling between the linear sinuous and varicose modes.
The varicose components contribute to the breakup of the sheet. Both ‘
$1k,v$
’ and ‘
$2k,v$
’ components contract near
$x=50$
, making the two interfaces contact. The function of the ‘
$1k,v$
’ component, which represents the linear varicose mode ‘
$1,v$
’, is more important not only because its amplitude is larger, but also because it primarily determines the location of breakup. Its period is twice as long as that of ‘
$2k,v$
’, so the thinning caused by ‘
$2k,v$
’ near
$x=37$
is counteracted by the dilation of ‘
$1k,v$
’, and only near the location where ‘
$1k,v$
’ contracts does the sheet break up, making the ligaments interspaced by one wavelength. The ‘
$2k,v$
’ component plays the role of accelerating the thinning and disintegration. According to (3.7), noting that the mode ‘
$22,v$
’ is also negligible, the ‘
$2k,v$
’ component is determined by modes ‘
$21,ss$
’ and ‘
$21,vv$
’, the first harmonics of the linear sinuous and varicose modes. According to figure 6(c), ‘
$21,ss$
’ and ‘
$21,vv$
’ have parallel amplitudes at
$t=57.5$
, so they have similar contributions to breakup. It should be pointed out that the first harmonic of the linear sinuous mode ‘
$21,ss$
’ does not have a predominant impact on the breakup process, as most previous researchers have conjectured, because, although it has the largest temporal growth rate, its initial amplitude represented by
$|\text{}^{2r,ss}\hat{{\it\eta}}|$
is too small, making its effect less important than the influence caused by varicose modes ‘
$1,v$
’ and ‘
$21,vv$
’.
Figure 8. Shape of sheet for varicose mode at
$t=40$
: (a) weakly nonlinear analysis; (b) numerical simulation.
Finally, since the nonlinear instability of the pure varicose mode is a very fundamental problem and its property differs significantly from other situations, the instability mechanism for
${\it\theta}={\rm\pi}$
is briefly discussed. Figure 5(i, j) indicates that the discrepancies between the weakly nonlinear analysis and numerical simulation become significant when the sheet breaks up. This may be because, as figure 6(e) shows, the second-order nonlinear amplitude becomes proportional to that of the linear first order at the final stage; this indicates strong nonlinearity and contradicts the basis of the weakly nonlinear analysis. Therefore, we only analyse the sheet at
$t=40$
when the weakly nonlinear analysis still agrees well with numerical simulation, which is shown in figure 8. Since in this condition
$r_{s}=0$
, no sinuous mode exists and the sheet has a symmetric profile. As figure 8(a) shows, the dilatational region of the first harmonic of the linear varicose mode ‘
$21,vv$
’ overlaps the contracting region of the linear varicose mode ‘
$1,v$
’, and thus their interference produces a long thread. The weakly nonlinear analysis can explain only the formation of this thread, while the numerical simulation in figure 5(i, j) shows that it is maintained till the sheet breaks up.
3.2. Effect of flow parameters
The effects of wavenumber
$k$
, initial disturbance amplitude
${\it\eta}_{0}$
, Weber number
$\mathit{We}$
and gas-to-liquid density ratio
${\it\rho}$
on the nonlinear instability of sheets are studied in this section. We only consider the case of initial phase difference
${\it\theta}={\rm\pi}/2$
, in which
$r_{s}=0.5$
, so that the linear sinuous and varicose modes have identical disturbance amplitudes initially.
3.2.1. Effect of wavenumber
The complex temporal frequency as a function of wavenumber
$k$
, i.e. the dispersion relation, is given in figure 9. For most of the wavenumber range, the linear sinuous mode has a much larger temporal growth rate
$\text{Im}({\it\omega})$
, and the phase velocity represented by
$\text{Re}({\it\omega})/k$
is different for the linear sinuous and varicose modes. Figure 9 also reveals that the linear varicose mode has greater temporal growth rates at very large wavenumbers, and that it has a larger cutoff wavenumber (the wavenumber corresponding to zero temporal growth rate) for this mode.
Figure 9. Dispersion relation:
$\mathit{We}=5$
,
${\it\rho}=0.1$
.
Figure 10. Relative second-order disturbance amplitudes at different wavenumbers:
$\mathit{We}=5$
,
${\it\rho}=0.1$
.
The relative second-order disturbance amplitudes, which indicate the second-order initial amplitudes, are given in figure 10. For most wavenumbers,
$|\text{}^{2r,vv}\hat{{\it\eta}}|$
is the largest and
$|\text{}^{2r,ss}\hat{{\it\eta}}|$
is the smallest, so the first harmonic of the linear varicose mode ‘
$21,vv$
’ has the largest initial amplitude and the first harmonic of the linear sinuous mode ‘
$21,ss$
’ has the smallest initial amplitude. Figure 10 also shows that, while
$|\text{}^{2r,ss}\hat{{\it\eta}}|$
increases monotonically with increase of wavenumber
$k$
, indicating that short waves have stronger nonlinear effects, both
$|\text{}^{2r,vv}\hat{{\it\eta}}|$
and
$|\text{}^{2r,sv}\hat{{\it\eta}}|$
have a peak at a certain wavenumber.
Our examination shows that, unfortunately, in the wavenumber range
$0<k<0.39$
where both modes grow (required by the solvability condition discussed in § 2.2.3), only at approximately
$k>0.24$
does the weakly nonlinear analysis agree well with numerical simulation when the sheet breaks up. For
$0<k<0.1$
, the small values of
$\text{Im}({\it\omega}_{1v})$
and
$|\text{}^{2r,ss}\hat{{\it\eta}}|$
make the breakup time very long, resulting in extremely large wave amplitudes (dominated by the linear sinuous mode ‘
$1,s$
’). This makes the nonlinearity significantly strong and the weakly nonlinear analysis is no longer appropriate to study the breakup. For
$0.1<k<0.24$
, the large values of
$|\text{}^{2r,sv}\hat{{\it\eta}}|$
and
$|\text{}^{2r,vv}\hat{{\it\eta}}|$
bring about large initial second-order velocities of the weakly nonlinear analysis, making it less precise at
$t=0$
. Consequently, the initial vortex sheet strengths of the numerical simulation given by (2.58) and (2.59) do not induce an initial velocity field that agrees well with that of the weakly nonlinear analysis, and this certainly brings about significant disparities between the analytical and numerical results in the subsequent disturbance evolution. Only at
$k>0.24$
are all the first-order and second-order amplitudes not very large when the sheet breaks up, and the weakly nonlinear analysis agrees qualitatively well with numerical simulation. The present study only focuses on this range.
Figure 11. Shapes of sheets at different wavenumbers: (a,b) the results of numerical simulation and (c,d) the results of weakly nonlinear analysis;
$\mathit{We}=5$
,
${\it\rho}=0.1$
,
${\it\eta}_{0}=0.2$
; (a,c)
$k=0.25$
,
$t=57.5$
and (b,d)
$k=0.35$
,
$t=60$
.
Table 1. Disturbance amplitudes of various modes at different wavenumbers
$k$
when sheets break up (predicted by weakly nonlinear analysis);
$\mathit{We}=5$
,
${\it\rho}=0.1$
,
${\it\eta}_{0}=0.2$
.
The shapes of sheets at two wavenumbers
$k=0.25$
and
$k=0.35$
are given in figure 11. The most prominent characteristic is that an increase of wavenumber results in a reduction of wave amplitude. The reason is explained in the weakly nonlinear analysis: the amplitude of the linear sinuous mode ‘
$1,s$
’ becomes significantly smaller at
$k=0.35$
(see table 1). This can be easily understood from the dispersion relation shown in figure 9: the increase of
$k$
from 0.25 to 0.35 leads to a significant reduction of sinuous growth rate
$\text{Im}({\it\omega}_{1s})$
but a slight increase of varicose growth rate
$\text{Im}({\it\omega}_{1v})$
. For the amplitudes of various modes given in table 1, while the linear varicose mode ‘
$1,v$
’ is enhanced significantly due to the increase of
$\text{Im}({\it\omega}_{1v})$
, modes ‘
$21,ss$
’ and ‘
$21,vv$
’ do not change very much, as the rise of
$|\text{}^{2r,ss}\hat{{\it\eta}}|$
and the reduction of
$|\text{}^{2r,vv}\hat{{\it\eta}}|$
with the increase of
$k$
(see figure 10) offset the corresponding change in temporal growth rates. Therefore, the breakup mechanism does not vary a lot with the increase of wavenumber.
The cutoff wavenumbers are generally in proportion to
${\it\rho}\mathit{We}$
. In the following analysis, we only consider the case of
$k=0.5{\it\rho}\mathit{We}$
. At this wavenumber, the instability behaviour is conspicuous and the weakly nonlinear analysis agrees well with numerical simulation when the sheet breaks up.
3.2.2. Effect of initial disturbance amplitude
The shapes of the sheets at two different initial disturbance amplitudes,
${\it\eta}_{0}=0.2$
and
${\it\eta}_{0}=0.01$
, are given in figure 12. It is easy to understand that the sheet with a smaller initial amplitude requires a longer time for the disturbance to grow to one that is large enough to cause breakup, so the breakup time for
${\it\eta}_{0}=0.01$
is significantly longer than that for
${\it\eta}_{0}=0.2$
. Apart from this obvious difference, the instability behaviour of
${\it\eta}_{0}=0.01$
exhibits many other different properties from that of
${\it\eta}_{0}=0.2$
which was studied in detail in § 3.1: it has a larger wave amplitude, its wave’s upgrade and downgrade have nearly the same slopes, and it has two locations of attenuation in one wavelength.
Figure 12. Shapes of sheets at different initial disturbance amplitudes: (a,b) the results of numerical simulation and (c,d) the results of weakly nonlinear analysis;
$\mathit{We}=5$
,
${\it\rho}=0.1$
,
$k=0.25$
; (a,c)
${\it\eta}_{0}=0.2$
,
$t=57.5$
, (b)
${\it\eta}_{0}=0.01$
,
$t=119$
, and (d)
${\it\eta}_{0}=0.01$
,
$t=114$
.
Figure 13. Effect of initial amplitude on the temporal evolution of the amplitudes of various modes:
$\mathit{We}=5$
,
${\it\rho}=0.1$
,
$k=0.25$
; (a)
${\it\eta}_{0}=0.2$
and (b)
${\it\eta}_{0}=0.01$
.
The reason for these differences can be explained by the weakly nonlinear analysis. As figure 13 shows, with
${\it\eta}_{0}$
reducing from 0.2 to 0.01, at the later stage of disturbance evolution, while the amplitudes of modes ‘
$1,v$
’, ‘
$21,sv$
’ and ‘
$21,vv$
’ decrease significantly, the amplitude of ‘
$21,ss$
’ increases extensively. The large amplitude of mode ‘
$1,s$
’ is not covered by figure 13, but at the figure’s upper limit,
${\it\eta}_{a}=2$
, the slope of
${\it\eta}_{0}=0.01$
is much steeper, indicating that it grows much faster, and so, when the sheet breaks up, the amplitude of mode ‘
$1,s$
’ is 8.019 for
${\it\eta}_{0}=0.01$
, much larger than 4.192 for
${\it\eta}_{0}=0.2$
. The increase of the ‘
$1,s$
’ mode’s amplitude results in a larger wave amplitude; the decrease of the ‘
$21,sv$
’ mode’s amplitude reduces its function in modulating the waves, making the slopes of upgrade and downgrade nearly identical. As for the attenuation, the domination of the ‘
$21,ss$
’ mode makes the thinning interspaced by one half-wavelength (see figure 12
d), the one-wavelength mode ‘
$1,v$
’ significantly reducing its amplitude and weakening its effect of eliminating one of the attenuations.
Figure 14. Amplification factor of modes ‘
$21,ss$
’ and ‘
$21,sv$
’:
$\mathit{We}=5$
,
${\it\rho}=0.1$
,
$k=0.25$
.
We give a further explanation of the dominance of modes ‘
$21,ss$
’ and ‘
$1,s$
’ at small initial amplitudes. The main reason for this phenomenon is that different modes grow at different speeds. Figure 14 gives the amplification factors of modes ‘
$21,ss$
’ and ‘
$21,sv$
’, denoted as
$\text{}^{21,ss}f$
and
$\text{}^{21,sv}f$
:
$$\begin{eqnarray}\text{}^{21,ss}f=\exp [2\,\text{Im}({\it\omega}_{1s})t],\quad \text{}^{21,sv}f=\exp [\text{Im}({\it\omega}_{1s})t+\text{Im}({\it\omega}_{1v})t].\end{eqnarray}$$
At
$t=57.5$
, the breakup time for
${\it\eta}_{0}=0.2$
,
$\text{}^{21,ss}f$
and
$\text{}^{21,sv}f$
have parallel values, but at a much longer breakup time,
$t=114$
, rendered by a much smaller initial amplitude
${\it\eta}_{0}=0.01$
,
$\text{}^{21,ss}f$
becomes far larger than
$\text{}^{21,sv}f$
, the reason being that its dominance in exponential growth rate becomes extensive at longer times. This explains why the reduction of
${\it\eta}_{0}$
leads to the domination of mode ‘
$21,ss$
’ over modes ‘
$21,sv$
’ and ‘
$21,vv$
’. For the first-order modes ‘
$1,s$
’ and ‘
$1,v$
’, this analysis is not applicable, as their dependences on
${\it\eta}_{0}$
are different: the first-order modes are proportional to
${\it\eta}_{0}$
, while the second-order modes are proportional to
${\it\eta}_{0}^{2}$
. The amplitudes of modes ‘
$1,s$
’ and ‘
$21,ss$
’, denoted as
$\text{}^{1,s}A$
and
$\text{}^{21,ss}A$
, have a simple relation: according to their expressions
$$\begin{eqnarray}\text{}^{1,s}A=2{\it\eta}_{0}|\text{}^{1,s}\hat{{\it\eta}}_{j}|\exp [\text{Im}({\it\omega}_{1s})t],\quad \text{}^{21,ss}A=2{\it\eta}_{0}^{2}|\text{}^{21,ss}\hat{{\it\eta}}_{j}|\exp [2\,\text{Im}({\it\omega}_{1s})t],\end{eqnarray}$$
using the definition of
$\text{}^{2r,ss}\hat{{\it\eta}}$
given in (2.44), we have
$$\begin{eqnarray}\text{}^{1,s}A=\sqrt{\frac{2(\text{}^{21,ss}A)}{|\text{}^{2r,ss}\hat{{\it\eta}}|}}.\end{eqnarray}$$
Therefore, the increase of the amplitude of mode ‘
$21,ss$
’ results in a corresponding amplification of mode ‘
$1,s$
’, explaining the large amplitude of the linear sinuous mode at a small initial disturbance amplitude. For the linear varicose mode ‘
$1,v$
’, such a simple relation does not exist. We can only explain a little vaguely that, because its growth rate
$\text{Im}({\it\omega}_{1v})$
is small, the increment of its amplification factor
$\text{}^{1,v}f=\exp [\text{Im}({\it\omega}_{1v})t]$
as a result of a longer breakup time is not enough to compensate the effect of the reduction in
${\it\eta}_{0}$
, leading to the decrease of its amplitude.
3.2.3. Effect of Weber number and gas-to-liquid density ratio
The effects of Weber number and gas-to-liquid density ratio have been extensively studied by previous researchers. It is now well known that an increase in Weber number or gas-to-liquid density ratio significantly increases temporal growth rates, thereby prominently enhancing instability, and that the instability is mainly caused by the interaction between liquid and gas. In the present study, while the typical discussion, i.e. the effect of
$\mathit{We}$
and
${\it\rho}$
on the temporal growth rates, second-order amplitudes and breakup time, is briefly noted, we also reveal some interesting properties of unstable waves, and give explicit explanations for them through the weakly nonlinear analysis.
Figure 15. Shapes of sheets at different Weber numbers: (a,b) the results of numerical simulation and (c,d) the results of weakly nonlinear analysis;
${\it\rho}=0.1$
,
$k=0.5{\it\rho}\mathit{We}$
,
${\it\eta}_{0}=0.2$
; (a,c)
$\mathit{We}=3$
,
$t=124.8$
, (b)
$\mathit{We}=6$
,
$t=46$
, and (d)
$\mathit{We}=6$
,
$t=45$
.
Table 2. Effect of Weber number on temporal growth rates and relative second-order disturbance amplitudes;
${\it\rho}=0.1$
,
$k=0.5{\it\rho}\mathit{We}$
.
Table 3. Effect of gas-to-liquid density ratio on temporal growth rates and relative second-order disturbance amplitudes;
$\mathit{We}=3$
,
$k=0.5{\it\rho}\mathit{We}$
.
Figure 16. Shapes of sheets at different gas-to-liquid density ratios: (a,b) the results of numerical simulation and (c,d) the results of weakly nonlinear analysis;
$\mathit{We}=3$
,
$k=0.5{\it\rho}\mathit{We}$
,
${\it\eta}_{0}=0.2$
; (a,c)
${\it\rho}=0.1$
,
$t=124.8$
, (b)
${\it\rho}=0.13$
,
$t=87$
, and (d)
${\it\rho}=0.13$
,
$t=83.5$
.
The effects of Weber number and gas-to-liquid density ratio on the major parameters, i.e. temporal growth rates and relative second-order disturbance amplitudes, are given in tables 2 and 3. An increase in Weber number and gas-to-liquid density ratio results in an extensive increase in the temporal growth rates, and this leads to a much shorter breakup time (see figures 15 and 16). For the second-order amplitudes, the amplitude of mode ‘
$21,ss$
’ increases significantly, while the other two modes only change slightly. But this does not bring about the predominance of ‘
$21,ss$
’ when the sheet breaks up, because the effect of this increase in second-order initial amplitude is somewhat marginal compared to the prominent effect of growth rates. A further examination shows that such a change in Weber number and gas-to-liquid density ratio does not result in a qualitative change in the breakup mechanism.
The numerical simulation shows that, as the Weber number and gas-to-liquid density ratio increase, the protuberance near the point of breakup becomes larger. This indicates that a more extensive interaction between gas and liquid can significantly enhance the local nonlinearity and distort the sheet more severely. The protuberances on the sheets have already been observed in the experiments of Asare et al. (Reference Asare, Takahashi and Hoffman1981) and Tammisola et al. (Reference Tammisola, Sasaki, Lundell, Matsubara and Söderberg2011); their experiments also showed that these protuberances are very likely to break up into smaller drops. But the mechanism of their formation and evolution is not very clear so far, and we look forward to future studies on this issue.
Next, we use the weakly nonlinear analysis to explain two interesting phenomena in figures 15 and 16. The first phenomenon is that, in figure 15, at
$\mathit{We}=3$
, the sheet breaks up on its wave crest, but at
$\mathit{We}=6$
, the sheet breaks up near its wave trough. The main reason is the shift of the linear varicose mode ‘
$1,v$
’, as the discussion of figure 7 has indicated that the sheet breaks up near the contraction of mode ‘
$1,v$
’ (see figure 15
c,d). At
$\mathit{We}=3$
, the contraction of mode ‘
$1,v$
’ is near the wave crest, but at
$\mathit{We}=6$
, its contraction moves rightwards, resulting in the corresponding change in breakup location.
The second phenomenon is that, in figure 16, using
$y=0$
as the origin, at
${\it\rho}=0.1$
, the depth of the wave trough is larger than the height of the wave crest, but at
${\it\rho}=0.13$
, the situation is the opposite, the height of the wave crest being larger than the depth of the wave trough. This is mainly caused by the first harmonic generated by the coupling between the linear sinuous mode and the varicose mode, i.e. mode ‘
$21,sv$
’. At
${\it\rho}=0.1$
, the troughs of mode ‘
$21,sv$
’ overlap both the peak and trough of the basic sinuous wave, shifting both of them downwards; while at
${\it\rho}=0.13$
, it is the crests of mode ‘
$21,sv$
’ that are in the positions of the peak and trough of the basic sinuous wave, and thus shift both of them upwards. This phenomenon has again highlighted the importance of the nonlinear coupling between the linear sinuous and varicose modes.
3.3. Remarks and future work
The parameter range in the present study is small Weber numbers (
$\mathit{We}\leqslant 6$
) and large gas-to-liquid density ratios (
${\it\rho}\sim 0.1$
). The present study neglects gravity; given a fixed small Weber number, this is applicable for very thin liquid sheets, so that the Froude number
$\mathit{Fr}=U^{2}/ga=\mathit{We}{\it\sigma}/{\it\rho}_{l}ga^{2}$
, which describes the ratio of inertia to gravity, is very large. For example, using the properties of water, that is,
${\it\rho}_{l}=1000~\text{kg}~\text{m}^{-3}$
and
${\it\sigma}=0.073~\text{N}~\text{m}^{-1}$
, if
$\mathit{We}=5$
,
$a=4~{\rm\mu}\text{m}$
, then
$U=9.55~\text{m}~\text{s}^{-1}$
and
$\mathit{Fr}=2.33\times 10^{6}$
. Such thin films have appeared in the study of Squire (Reference Squire1953): one of his experiments was in the condition
$a=4~{\rm\mu}\text{m}$
,
$U=20~\text{m}~\text{s}^{-1}$
, and thus
$\mathit{Fr}=1.02\times 10^{7}$
; the Froude number in his experiment was even larger. High gas densities exist in situations such as the combustion chamber of a rocket engine (the gas pressure can be as high as 10 MPa).
A more common situation is that the sheet is not very thin: the thickness is of the order of 1 mm, and the gas has the density of air. In this condition, the Weber number is relatively large and the gas-to-liquid density ratio is very small. For example, in the experiment of Tammisola et al. (Reference Tammisola, Sasaki, Lundell, Matsubara and Söderberg2011), which also used water, the conditions were
$a=0.5~\text{mm}$
,
$U=7.15~\text{m}~\text{s}^{-1}$
,
${\it\rho}_{g}=1.3~\text{kg}~\text{m}^{-3}$
, and thus
$\mathit{We}=350$
,
${\it\rho}=0.0013$
,
$\mathit{Fr}=10433$
; note that gravity also had little effect in their study. Similar conditions are broadly used in industrial applications, such as paper-making, atomization and sprays, and thus they have been imposed with more focus by previous studies. Unfortunately, the present study cannot address such situations. The weakly nonlinear analysis can easily give solution for these conditions, but the numerical simulation encounters some difficulties. Our numerical experiments indicate that, in these conditions, small-scale features develop on the fluid interfaces before the sheet’s two interfaces contact, similar to the ‘Kelvin–Helmholtz spirals’ obtained in Hou et al. (Reference Hou, Lowengrub and Shelley1997). Our current numerical code is not capable of handling these small-scale features, and in order to overcome this difficulty, local grid refinement is necessary. This issue is left for our future research.
4. Conclusions
The nonlinear temporal instability of an inviscid planar liquid sheet moving in an inviscid gas is studied, with an emphasis on the case in which both the linear sinuous mode and the linear varicose mode exist simultaneously. Two different approaches, the weakly nonlinear analysis using a second-order perturbation expansion and numerical simulation using a boundary integral method, have been employed. While the numerical solution is relatively precise, the approximate analytical solution gives important explanations to it. There are good agreements between them. For most of the time of disturbance evolution, the shapes predicted by the two methods are nearly the same. When the sheet nearly breaks up, some disparities between them emerge, but the weakly nonlinear analysis’s predictions of the sheet’s general shape, breakup time and breakup location are still very close to those of the numerical simulation in most conditions. The Fourier analysis also shows that the two methods have good agreements in the fundamental wave and its first harmonic when the sheet breaks up. These comparisons justify the validity of the weakly nonlinear analysis in studying the mechanism of instability and breakup.
The weakly nonlinear analysis shows that both the first harmonic of the linear sinuous mode and the first harmonic of the linear varicose mode are varicose; they contribute to the sheet’s breakup. But the first harmonic generated by the coupling between the linear sinuous and varicose modes is sinuous; it plays an important role in modulating the shapes of the waves. In some conditions, it makes the upgrade and downgrade of the wave have different slopes, while in some other conditions, it brings about a difference between the height of the wave crest and the depth of the wave trough. The numerical simulation shows that a protuberance appears near the location of breakup; its size increases with the increase of Weber number and gas-to-liquid density ratio. Such a strong local nonlinear effect cannot be explained by the present weakly nonlinear analysis.
We examined the temporal evolution of disturbances leading to the final breakup with various initial phase differences between the upper and lower interfaces, representing various proportions of initial sinuous and varicose amplitudes. The shapes of sheets are dominated by the linear sinuous mode’s characteristics when their amplitudes grow large, except for the case of initial varicose disturbance, in which the linear sinuous mode does not exist and the sheet maintains a symmetric profile. When the initial phase difference is
${\rm\pi}/2$
and the linear sinuous and varicose modes have identical initial amplitudes, the breakup is mainly caused by three modes according to the second-order weakly nonlinear analysis: the linear varicose mode, the first harmonic of the linear sinuous mode and the first harmonic of the linear varicose mode. The role of the linear varicose mode is especially important because, in most conditions, its amplitude is relatively large and it primarily determines the location of breakup.
The effects of flow parameters, including the wavenumber, initial disturbance amplitude, Weber number and gas-to-liquid density ratio, have been investigated. Our examination shows that, when sheets break up, the weakly nonlinear analysis agrees well with numerical simulation only at relatively large wavenumbers, and in this range increasing the wavenumber results in a reduction of wave amplitude, as this will reduce the sinuous growth rate. Reducing the initial disturbance amplitude significantly increases the breakup time, and consequently makes the first harmonic of the linear sinuous mode the dominant mode in breakup, because the effect of its fast growth rate predominates at longer times. Increasing the Weber number or gas-to-liquid density ratio significantly increases the linear sinuous and varicose growth rates, and therefore extensively reduces the breakup time, significantly enhancing instability.
Acknowledgements
First of all the authors would like to thank the three anonymous referees, who gave valuable comments and suggestions on the manuscript, which significantly improved its quality. A number of researchers have helped us with our boundary integral simulation; they are (in time sequence of our communication) Professor X. Jiang, Professor W. A. Sirignano, Professor R. H. Rangel, Dr C. Mehring, Professor L. Ying, Professor T. Beale and Professor M. Shelley. Their help was essential for our work and we thank all of them very much. We also thank Professor W. A. Sirignano for his advice on our discussion, the most important one being the suggestion to use Fourier analysis. Finally, we thank Professor H. A. Stone for giving us important information on the application of the boundary integral method. This work was supported by China’s National Nature Science Funds (support no. 11272036).
Appendix A. Solution for second-order velocity potential
The solution is as follows:
$$\begin{eqnarray}\displaystyle & \displaystyle \text{}^{21,ss}\hat{{\it\phi}}_{l}=\frac{\text{i}({\it\omega}_{1s}-k)[k\tanh k\text{}^{1,s}\hat{{\it\eta}}_{1}^{2}-\text{}^{21,ss}\hat{{\it\eta}}_{1}]}{k\sinh 2k}\cosh 2ky, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \text{}^{21,ss}\hat{{\it\phi}}_{gj}=\text{i}{\it\omega}_{1s}\left(\frac{\text{}^{21,ss}\hat{{\it\eta}}_{1}}{k}+\text{}^{1,s}\hat{{\it\eta}}_{1}^{2}\right)\exp [2k+(-1)^{j}2ky], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \text{}^{21,s\overline{s}}\hat{{\it\phi}}_{l}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \text{}^{21,s\overline{s}}\hat{{\it\phi}}_{gj}=\frac{\text{i}[-{\it\rho}(|{\it\omega}_{1s}|^{2}-{\it\omega}_{1s}^{2})+\frac{1}{2}|{\it\omega}_{1s}-k|^{2}(1+\tanh ^{2}k)-({\it\omega}_{1s}-k)^{2}]|\text{}^{1,s}\hat{{\it\eta}}_{1}|^{2}}{{\it\rho}({\it\omega}_{1s}-\overline{{\it\omega}}_{1s})},\quad & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \text{}^{21,vv}\hat{{\it\phi}}_{l}=\frac{\text{i}({\it\omega}_{1v}-k)[k\coth k\text{}^{1,v}\hat{{\it\eta}}_{1}^{2}-\text{}^{21,vv}\hat{{\it\eta}}_{1}]}{k\sinh 2k}\cosh 2ky, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \text{}^{21,vv}\hat{{\it\phi}}_{gj}=\text{i}{\it\omega}_{1v}\left(\frac{\text{}^{21,vv}\hat{{\it\eta}}_{1}}{k}+\text{}^{1,v}\hat{{\it\eta}}_{1}^{2}\right)\exp [2k+(-1)^{j}2ky], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \text{}^{21,v\overline{v}}\hat{{\it\phi}}_{l}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \text{}^{21,v\overline{v}}\hat{{\it\phi}}_{gj}=\frac{\text{i}\{-{\it\rho}[|{\it\omega}_{1v}|^{2}-{\it\omega}_{1v}^{2}]+\frac{1}{2}|{\it\omega}_{1v}-k|^{2}(1+\coth ^{2}k)-({\it\omega}_{1v}-k)^{2}\}|\text{}^{1,v}\hat{{\it\eta}}_{1}|^{2}}{{\it\rho}({\it\omega}_{1v}-\overline{{\it\omega}}_{1v})},\qquad & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \text{}^{21,sv}\hat{{\it\phi}}_{l} & = & \displaystyle \{\!\text{2i}k[({\it\omega}_{1s}-k)\tanh k+({\it\omega}_{1v}-k)\coth k]\text{}^{1,s}\hat{{\it\eta}}_{1}\text{}^{1,v}\hat{{\it\eta}}_{1}\nonumber\\ \displaystyle & & \displaystyle -\,\text{i}({\it\omega}_{1s}+{\it\omega}_{1v}-2k)\text{}^{21,sv}\hat{{\it\eta}}_{1}\!\}\frac{\sinh 2ky}{2k\cosh 2k},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \text{}^{21,sv}\hat{{\it\phi}}_{gj}=(-1)^{j+1}\text{i}({\it\omega}_{1s}+{\it\omega}_{1v})\left(\frac{\text{}^{21,sv}\hat{{\it\eta}}_{1}}{2k}+\text{}^{1,s}\hat{{\it\eta}}_{1}\text{}^{1,v}\hat{{\it\eta}}_{1}\right)\exp \{2k[1+(-1)^{j}y]\}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \text{}^{21,s\overline{v}}\hat{{\it\phi}}_{l}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \text{}^{21,s\overline{v}}\hat{{\it\phi}}_{gj}=-(-1)^{j+1}\frac{(1-{\it\rho})}{{\it\rho}}\text{i}({\it\omega}_{1s}-\overline{{\it\omega}}_{1v})\text{}^{1,s}\hat{{\it\eta}}_{1}\text{}^{1,v}\overline{\hat{{\it\eta}}}_{1}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \text{}^{22,s}\hat{{\it\phi}}_{l}=\text{i}\left(1-\frac{{\it\omega}_{2s}}{2k}\right)\frac{\sinh (2ky)}{\cosh (2k)}\text{}^{22,s}\hat{{\it\eta}}_{1}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \text{}^{22,s}\hat{{\it\phi}}_{gj}=(-1)^{j+1}\frac{\text{i}{\it\omega}_{2s}}{2k}\exp [2k+(-1)^{j}2ky]\text{}^{22,s}\hat{{\it\eta}}_{1}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \text{}^{22,\overline{s}}\hat{{\it\phi}}_{l}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \text{}^{22,\overline{s}}\hat{{\it\phi}}_{gj}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \text{}^{22,v}\hat{{\it\phi}}_{l}=\text{i}\left(1-\frac{{\it\omega}_{2v}}{2k}\right)\frac{\cosh (2ky)}{\sinh (2k)}\text{}^{22,v}\hat{{\it\eta}}_{1}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \text{}^{22,v}\hat{{\it\phi}}_{gj}=\frac{\text{i}{\it\omega}_{2v}}{2k}\exp [2k+(-1)^{j}2ky]\text{}^{22,v}\hat{{\it\eta}}_{1}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \text{}^{22,\overline{v}}\hat{{\it\phi}}_{l}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \text{}^{22,\overline{v}}\hat{{\it\phi}}_{gj}=0. & \displaystyle\end{eqnarray}$$
