1. Introduction
Viscoelastic fluids exhibit a time-dependent response to deformation. They behave like elastic solids when exposed to deformations over short time scales but exhibit fluid-like properties in steady shear over long time scales. Interest around viscoelastic fluids grew as early as the 1940s and 1950s, influenced by industrial research and development in chemical engineering (Denn Reference Denn2004). Over the years, extensive efforts have been employed to combine polymers of widely dissimilar material structures to yield unique multilayered blends and properties. Currently, viscoelastic fluids find wide application in several commercial polymer processing operations like spray coating, spray drying, fibre spinning, food processing, etc. Introductory texts on this subject with the application of perturbation technique can be found in the works of Langlois & Rivlin (Reference Langlois and Rivlin1959, Reference Langlois and Rivlin1963) which investigated the effect of viscoelasticity in normal stress-driven converging flows and rectilinear flows, respectively. Liu & Durst (Reference Liu and Durst1998) and Brenn, Liu & Durst (Reference Brenn, Liu and Durst2000) reported that non-Newtonian sheets are more unstable than Newtonian sheets for both symmetric and asymmetric disturbances. Overall effects of elasticity number and viscosity were found to be weakly destabilizing and stabilizing, respectively. Brenn, Liu & Durst (Reference Brenn, Liu and Durst2001) and Liu & Liu (Reference Liu and Liu2006) considered three-dimensional disturbances in viscoelastic fluids and observed that the growth rate of two-dimensional disturbances was higher than that of three-dimensional disturbances. Yang et al. (Reference Yang, Qu, Fu and Gu2010) investigated breakup of gel propellant under the combined influence of sinuous and varicose modes of disturbance. Results were found to be consistent with previous literature which also reported that viscoelastic sheets were more unstable than Newtonian sheets. However, breakup time for the varicose mode of disturbance was found to be the same as that of Newtonian sheets. Ye et al. (Reference Ye, Yang, Fu, Ye, Yang and Fu2016) and later Alsharif (Reference Alsharif2019) performed temporal instability analysis of a viscoelastic compound jet using the Oldroyd B model. Both studies confirmed that viscoelastic sheets are more unstable than Newtonian compound jets but less unstable than inviscid compound jets. In addition to the linear studies discussed so far, the literature also contains detailed nonlinear investigations of the temporal evolution of two-dimensional disturbances in viscoelastic fluids. Nonlinear analysis of temporal instabilities in Poiseuille and Couette flow reported by Atalık & Keunings (Reference Atalık and Keunings2002) gave very interesting findings in the form of a dual effect of elasticity for the Oldroyd B model: increasing elasticity number had a destabilizing effect initially, followed by a re-stabilization effect. Wang et al. (Reference Wang, Yang, Xie and Chen2015) performed an in-depth nonlinear analysis of three rheological models, namely corotational Jeffrey's model, Oldroyd A model and Oldroyd B model, exposed to initial sinuous disturbances. The study showed that though different rheological models have different disturbance pressure, the interface displacement and disturbance velocities are the same. Nonlinear analysis of the varicose mode of disturbance by Xie et al. (Reference Xie, Yang, Wang and Qin2018) showed that the second-order displacement of the initial varicose mode was also varicose in nature, which caused sheet breakup at full wavelength interval and produced ligaments with two connected swells. However, the aforementioned nonlinear studies considered surrounding air to be quiescent. Atomizers that utilize the kinetic energy of ambient gases, such as twin-fluid atomizers, offer advantages in the form of good atomization quality at low pressure (Sovani, Sojka & Sivathanu Reference Sovani, Sojka and Sivathanu2000). Moreover, low sensitivity to fluid rheology also makes them a preferred choice in commercial operations such as spray coating, spray drying and process industries, which utilize viscoelastic fluids (Mujumdar et al. Reference Mujumdar, Huang and Dong Chen2010). The effect of non-zero unequal gas velocities on viscoelastic liquid sheets was analysed by Yang, Xu & Fu (Reference Yang, Xu and Fu2012) for para sinuous and para varicose modes. Similar to Newtonian sheets, a higher velocity difference across the interfaces increases the instability of non-Newtonian sheets. However, their study employed linear theory which can only predict the sheet behaviour at the onset of instabilities (Asare, Takahashi & Hoffman Reference Asare, Takahashi and Hoffman1981; Mitra, Li & Renksizbulut Reference Mitra, Li and Renksizbulut2001). Moreover, instabilities in Newtonian viscous flows are generally a direct consequence of the presence of nonlinear convective terms in the momentum equation. But in viscoelastic flows, the nonlinearity appears through both momentum convection and stress convection in constitutive equations, which further emphasize the need for nonlinear studies of viscoelastic fluids. The literature shows that most of the nonlinear studies that consider non-zero gas velocities either are performed within the inviscid approximation (Jazayeri & Li Reference Jazayeri and Li2000; Nath et al. Reference Nath, Mukhopadhyay, Sen and Tharakan2010; Nath, Mukhopadhyay & Datta Reference Nath, Mukhopadhyay and Datta2014; Dasgupta, Nath & Bhanja Reference Dasgupta, Nath and Bhanja2018) or involve Newtonian fluids (Yang, Chen & Wang Reference Yang, Chen and Wang2014; Dasgupta, Nath & Bhanja Reference Dasgupta, Nath and Bhanja2019). Excellent reviews on this topic can also be found in the works of Denn (Reference Denn1990), Larson (Reference Larson1992) and James (Reference James2009).
Based on the above discussion, it can be concluded that previous literature lacks detailed nonlinear analyses of breakup of viscoelastic fluids in the presence of moving ambient gas medium. The present study considers corotating convected frames in the constitutive relations and employs the second-order perturbation technique to perform a weakly nonlinear study of temporal instabilities in a viscoelastic fluid surrounded by non-zero gas velocities. The main objective is to investigate the effect of viscoelastic properties on growth rate, breakup time and sheet behaviour at different gas velocities. Ideally, the final outcome of any atomization analysis should be droplet characteristics such as droplet size and velocity distribution. Early contributions in this area can be found in the works of Rosin & Rammler (Reference Rosin and Rammler1933) and Nukiyama & Tanasawa (Reference Nukiyama and Tanasawa1939) who developed empirical relations to define droplet characteristics with sufficient accuracy. Later, Rizk & Lefebvre (Reference Rizk and Lefebvre1985) and Bhatia et al. (Reference Bhatia, Dominick, Drust and Tropea1988) employed the modified Rosin–Rammler method and the log-hyperbolic method, respectively, to achieve better data fit for larger droplet size and a wide range of experimental data. Villermaux, Marmottant & Duplat (Reference Villermaux, Marmottant and Duplat2004) observed that the atomization process for Newtonian liquids can be precisely described by a fragmentation coalescence scenario of ligament and the final droplet size distribution can be represented by gamma distributions. In the absence of any experimental data, analytical methods such as maximum entropy formulation (MEF) and discrete probability function can be useful tools to predict spray characteristics (Li, Tankin & Renksizbulut Reference Li, Tankin and Renksizbulut1990; Li & Tankin Reference Li and Tankin1992; Nath et al. Reference Nath, Datta, Mukhopadhyay, Sen and Tharakan2011; Negeed Reference Negeed2011). The present study employs the modified MEF to elucidate the effect of viscoelasticity on droplet size distribution. The MEF involves the method of Lagrange multipliers consisting of a numerical procedure for the solution of a set of nonlinear equations. The Newton–Raphson method has been believed to be a trustable numerical method. But unfortunately, due to the involvement of exponential terms and not having sufficiently close initial guess value to the root of equations, the MEF has a propensity to diverge rapidly. To eliminate this tendency to diverge, a modified Newton–Raphson method with Taylor series expansion up to second order is considered in the present study.
2. Mathematical formulation
The study considers a two-dimensional viscoelastic sheet of thickness 2h enclosed by two inviscid gas streams, flowing with equal non-zero velocities (figure 1). The sheet is subjected to an initial sinuous mode of disturbance. Liquid density, surface tension and viscosity are represented as 
${\rho _l}$, 
${\sigma _l}$ and 
${\mu _l}$, respectively. Initially, the unperturbed liquid sheet has velocity only in the x direction, represented by 
${u_l}$. Surrounding gases are considered to be inviscid in nature and gas density is denoted by 
${\rho _g}$. In the unperturbed state, gas velocities at the two interfaces are parallel to the liquid flow. The gas velocities are represented by 
${u_g}$. Both liquid and gas flows are considered to be incompressible and the effect of gravity is neglected. Non-dimensionalization of constitutive relation, governing equation and boundary conditions is performed using the following scale: [length, time, density, velocity, stress] = [
$h,\; h/{u_l},\; {\rho _l},\; {u_l},\; {\rho _l}u_l^2$]. Reynolds number is represented as 
$\; Re = {\rho _l}{u_l}h/{\mu _l}$.
Schematic of planar viscoelastic sheet of liquid enclosed within two gas streams of equal non-zero velocities under the influence of sinuous disturbances.
The tensor form for the convected Jeffrey's corotational model in its dimensionless form is expressed as
\begin{equation}\boldsymbol{\tau } + {\boldsymbol{\lambda }_1}[\tilde{\boldsymbol{\tau }} - \boldsymbol{W}\boldsymbol{\cdot }\boldsymbol{\tau } + \boldsymbol{\tau }\boldsymbol{\cdot }\boldsymbol{W}] = \frac{1}{{Re}}[\boldsymbol{\gamma } + {\lambda _2}[\tilde{\boldsymbol{\gamma }} - \boldsymbol{W}\boldsymbol{\cdot }\boldsymbol{\gamma } + \boldsymbol{W}\boldsymbol{\cdot }\boldsymbol{\gamma }]].\end{equation}
Here, 
${\lambda _1}$ and 
${\lambda _2}$ represent stress relaxation time and deformation retardation time, respectively.
Parameter τ is the extra stress tensor and its material derivative is expressed as
\begin{equation}\tilde{\boldsymbol{\tau }} = \frac{{\partial \boldsymbol{\tau }}}{{\partial t}} + (\boldsymbol{V}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{\tau }.\end{equation}Here, V represents liquid velocity vector (u, v, 0).
Parameter γ is the strain tensor such that
\begin{equation}\boldsymbol{\gamma } = \boldsymbol{\nabla }\boldsymbol{V} + {(\boldsymbol{\nabla }\boldsymbol{V})^\textrm{T}}\;\textrm{and}\;\textrm{its}\;\textrm{material}\;\textrm{derivative}\;\textrm{is}\;\tilde{\boldsymbol{\gamma }} = \frac{{\partial \boldsymbol{\gamma }}}{{\partial t}} + (\boldsymbol{V}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{\gamma }.\end{equation}Parameter W represents the vorticity tensor such that
\begin{equation}\boldsymbol{W} = {\textstyle{1 \over 2}}[\boldsymbol{\nabla }\boldsymbol{V} - {(\boldsymbol{\nabla }\boldsymbol{V})^\textrm{T}}];\quad {W_{xx}} = {W_{yy}} = 0;\quad {W_{xy}} = {W_{yx}} = {u_y} - {v_x}.\end{equation} Two characteristic numbers, elasticity number (El) and time constant ratio 
$(\lambda )$, are introduced to define the rheological property of the viscoelastic fluid. Elasticity number represents the relative magnitude of the elastic stresses as compared to inertial stresses, whereas time constant ratio signifies the magnitude of stress relaxation time as compared to deformation retardation time. Corresponding expressions are gives as
\begin{equation}El = \frac{{{\lambda _1}}}{{Re}};\quad \lambda = \frac{{{\lambda _2}}}{{{\lambda _1}}}.\end{equation}
The upper gas–liquid interface is represented by j = 1 and the lower gas–liquid interface is represented by j = 2. Gas velocities are represented by velocity potential 
$\; {\phi _{gj}}$. In the absence of any perturbation, gas velocity potential is given as 
${\phi _{gj}}{|_{t = 0}} = Ux$, where U is the non-dimensional gas velocity 
$({u_{gj}}/{u_l})$. Liquid pressure and gas pressure in their non-dimensional form are represented as 
${P_l}$ and 
${P_g}$, respectively.
Location of non-dimensional interfaces after initial perturbation is given as
\begin{equation}y(x,t) = {( - 1)^{j + 1}} + {\eta _j}(x,t),\end{equation}
where 
${\eta _j}(x,t)$ is non-dimensional surface deformation. In the following sections, partial derivatives with respect to the x direction, y direction and time t are represented by the subscripts x, y and t, respectively.
Liquid mass and momentum conservation equations are presented in (2.7) and (2.8), respectively:
\begin{gather}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{V} = 0\quad \textrm{for}\; - 1 + {\eta _2} \le y \le 1 + {\eta _1},\end{gather}
\begin{gather}{\boldsymbol{V}_t} + (\boldsymbol{V}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{V} ={-} \boldsymbol{\nabla }{P_l} + \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\tau }\quad \textrm{for}\; - 1 + {\eta _2} \le y \le 1 + {\eta _1}.\end{gather}The kinematic boundary condition for liquid is given as
\begin{equation}v - {\eta _{j,t}} - (1 + {u_l}){\eta _{j,x}} = 0\quad \textrm{for}\; - 1 + {\eta _2} \le y \le 1 + {\eta _1},\end{equation}
where ‘1’ in 
$(1 + {u_l})$ represents the non-dimensional initial liquid flow.
The inviscid nature of the gas results in zero shear force at the gas–liquid interface (Li & Tankin Reference Li and Tankin1991):
\begin{equation}(\boldsymbol{\tau }\boldsymbol{\cdot }{\boldsymbol{n}_j}) \times {\boldsymbol{n}_j} = 0\;\textrm{at}\;y = {( - 1)^{j + 1}} + {\eta _j},\;\textrm{where}\;{\boldsymbol{n}_j}\;\textrm{is}\;\textrm{the}\;\textrm{unit}\;\textrm{vector}\;\textrm{normal}\;\textrm{to}\;\textrm{the}\;\textrm{surface}\textrm{.}\end{equation}Mass conservation equation for the gas phase is given as
\begin{equation}{\varphi _{gj,xx}} + {\varphi _{gj,yy}} = 0\quad \textrm{for}\;1 + {\eta _1} \le y \le \infty ,\; - \infty \le y \le - 1 + {\eta _2}.\end{equation}Kinematic boundary condition for the gas phase is expressed as
\begin{equation}{\varphi _{gj,y}} - {\eta _{j,t}} - {\varphi _{gj,x}}{\eta _{j,x}} = 0\;\textrm{at}\;y = {( - 1)^{j + 1}} + {\eta _j}.\end{equation}
To make the gas phase bounded and finite, gas velocity is considered to be zero at 
$y ={\pm} \infty $. Thus,
\begin{equation}{\varphi _{gj,y}} = 0\;\textrm{at}\;y ={\pm} \infty .\end{equation}The unsteady Bernoulli equation is employed to obtain gas pressure and is expressed as
\begin{equation}{p_{gj}} = {\textstyle{1 \over 2}}{U^2} - \rho [{({\varphi _{gj,x}})^2} + {({\varphi _{gj,y}})^2}] - \rho {\varphi _{gj,t}}\quad \textrm{for}\;1 + {\eta _1} \le y \le \infty ,\; - \infty \le y \le - 1 + {\eta _2}.\end{equation}The dynamic boundary condition is expressed as a balance between the normal stress difference across the gas–liquid interface and the surface tension force:
\begin{equation}- {P_l} + ({\boldsymbol{n}_j}\boldsymbol{\cdot }\boldsymbol{\tau }){\boldsymbol{n}_j} + {( - 1)^{j + 1}}\frac{1}{{We}}(\boldsymbol{\nabla }\boldsymbol{\cdot }{\boldsymbol{n}_j}) + {P_{gj}} = 0.\end{equation} Here, We represents the liquid Weber number and is expressed as 
$We = {\rho _l}u_l^2h/\sigma $.
All the parameters used in governing equations and boundary conditions are presented using power series of 
${\eta _0}$ such that
\begin{equation}({\eta _j},u,v,{P_l},\boldsymbol{\tau },\boldsymbol{\gamma },{\varphi _{gj}},{P_g}) = \sum\limits_{n = 1}^\infty {\eta _0^n} [{\eta _{jn}},{u_n},{v_n},{P_{ln}},{\boldsymbol{\tau }_n},{\boldsymbol{\gamma }_n},{\varphi _{gjn}},{P_{gjn}}].\end{equation}
All the liquid and gas parameters are approximated around the unperturbed interface 
$y = {( - 1)^j}$ using Taylor's series expansion. Finally, the following set of governing equations and boundary conditions are obtained.
First-order equations 
$(\eta _o^1)$:
\begin{gather}{\boldsymbol{\tau }_1} + {\lambda _1}({\boldsymbol{\tau }_{1,t}} + {\boldsymbol{\tau }_{1,x}}) = \frac{1}{{Re}}[{\boldsymbol{\gamma }_1} + {\lambda _2}({\boldsymbol{\gamma }_{1,t}} + {\boldsymbol{\gamma }_{1,x}})],\end{gather}
\begin{gather}{u_{1,x}} + {v_{1,y}} = 0\quad \textrm{for}\; - 1 \le y \le 1,\end{gather}
\begin{gather}{u_{1,t}} + {u_{1,x}} ={-} {P_{l1,x}} + \left( {\frac{\partial }{{\partial x}}{\tau_{1,xx}} + \frac{\partial }{{\partial y}}{\tau_{1,xy}}} \right)\quad \textrm{for}\; - 1 \le y \le 1,\end{gather}
\begin{gather}{v_{1,t}} + {v_{1,x}} ={-} {P_{l1,y}} + \left( {\frac{\partial }{{\partial x}}{\tau_{1,yx}} + \frac{\partial }{{\partial y}}{\tau_{1,yy}}} \right)\quad \textrm{for}\; - 1 \le y \le 1,\end{gather}
\begin{gather}{v_1} = {\eta _{j1,t}} + {\eta _{j1,x}}\;\textrm{at}\;y = {( - 1)^{j + 1}}\;\textrm{at}\;y = {( - 1)^{j + 1}},\end{gather}
\begin{gather}{\tau _{yx1}} = 0\;\textrm{at}\;y = {( - 1)^{j + 1}},\end{gather}
\begin{gather}{\varphi _{gj1,xx}} + {\varphi _{gj1,yy}} = 0\quad \textrm{for}\;1 \le y < \infty ,\; - \infty \le y < - 1,\end{gather}
\begin{gather}{\varphi _{gj1,y}} - {\eta _{j1,t}} - U{\eta _{j1,x}} = 0\;\textrm{at}\;y = {( - 1)^{j + 1}},\end{gather}
\begin{gather}{\varphi _{gj1,y}} = 0\;\textrm{at}\;y ={\pm} \infty ,\end{gather}
\begin{gather}- {P_{l1}} + {\tau _{yy1}} + {( - 1)^j}\frac{1}{{We}}{\eta _{j1,xx}} - \rho {\varphi _{gj1,t}} - \rho U{\varphi _{gj1,x}} = 0.\end{gather}
Second-order equations 
$(\eta _o^2)$:
\begin{gather}\begin{array}{ccccc} & {\boldsymbol{\tau }_2} + {\lambda _1}\left( {\dfrac{\partial }{{\partial t}}{\boldsymbol{\tau }_2} + \dfrac{\partial }{{\partial x}}{\tau_2} + {u_1}\dfrac{\partial }{{\partial x}}{\boldsymbol{\tau }_1} + {v_1}\dfrac{\partial }{{\partial y}}{\boldsymbol{\tau }_1} - {\boldsymbol{W}_1}\boldsymbol{\cdot }{\boldsymbol{\tau }_1} + {\boldsymbol{\tau }_1}\boldsymbol{\cdot }{\boldsymbol{W}_1}} \right)\\ & \quad = \dfrac{1}{{Re}}\left[ {{\boldsymbol{\gamma }_2} + {\lambda_2}\left( {\dfrac{\partial }{{\partial t}}{\boldsymbol{\gamma }_2} + \dfrac{\partial }{{\partial x}}{\boldsymbol{\gamma }_2} + {u_1}\dfrac{\partial }{{\partial x}}{\boldsymbol{\gamma }_1} + {v_1}\dfrac{\partial }{{\partial y}}{\boldsymbol{\gamma }_1} - {\boldsymbol{W}_1}\boldsymbol{\cdot }{\boldsymbol{\gamma }_1} + {\boldsymbol{\gamma }_1}\boldsymbol{\cdot }{\boldsymbol{W}_1}} \right)} \right], \end{array}\end{gather}
\begin{gather}{u_{2,x}} + {v_{2,y}} = 0\quad \textrm{for}\; - 1 \le y \le 1,\end{gather}
\begin{gather}{u_{2,t}} + {u_{2,x}} + {u_1}{u_{1,x}} + {v_1}{u_{1,y}} ={-} {P_{l2,x}} + \left( {\frac{\partial }{{\partial x}}{\tau_{2,xx}} + \frac{\partial }{{\partial y}}{\tau_{2,xy}}} \right)\quad \textrm{for}\; - 1 \le y \le 1, \end{gather}
\begin{gather}{v_{2,t}} + {v_{2,x}} + {u_1}{v_{1,x}} + {v_1}{v_{1,y}} ={-} {P_{l2,y}} + \left( {\frac{\partial }{{\partial x}}{\tau_{2,yx}} + \frac{\partial }{{\partial y}}{\tau_{2,yy}}} \right)\quad \textrm{for}\; - 1 \le y \le 1,\end{gather}
\begin{gather}{v_2} = {\eta _{j2,t}} + {\eta _{j2,x}} + {u_1}{\eta _{j1,x}} - {\eta _{j1}}{v_{1,y}}\;\textrm{at}\;y = {( - 1)^{j + 1}},\end{gather}
\begin{gather}{\tau _{2,yx}} + {\eta _{j1}}\left( {\frac{\partial }{{\partial y}}{\tau_{1,yx}}} \right) + ({\tau _{1,yy}} - {\tau _{1,xx}}){\eta _{j1,x}} = 0\;\textrm{at}\;y = {( - 1)^{j + 1}},\end{gather}
\begin{gather}{\varphi _{gj2,xx}} + {\varphi _{gj2,yy}} = 0\quad \textrm{for}\;1 \le y < \infty , - \infty < y \le - 1,\end{gather}
\begin{gather}{\varphi _{gj2,y}} - {\eta _{j2,t}} - U{\eta _{j2,x}} = {\eta _{j1,x}}{\varphi _{gj1,x}} - {\eta _{j1}}{\varphi _{gj1,yy}}\;\textrm{at}\;y = {( - 1)^{j + 1}},\end{gather}
\begin{gather}{\varphi _{gj2,y}} = 0\;\textrm{at}\;y ={\pm} \infty ,\end{gather}
\begin{gather}\begin{array}{ccccc} & - {P_{l2}} - {\eta _{j1}}{P_{l1,y}} - [{\eta _{j1,x}}({\tau _{1,xy}} + {\tau _{1,yx}})] + \left[ {{\tau_{2,yy}} + {\eta_{j1}}\dfrac{\partial }{{\partial y}}{\tau_{1,yy}}} \right]\\ & \quad + {( - 1)^j}\dfrac{1}{{We}}({\eta _{j2,xx}}) - \rho ({\varphi _{gj2,t}} + U{\varphi _{gj2,x}}) - \dfrac{1}{2}\rho (\varphi _{gj1,x}^2 + \varphi _{gj1,y}^2)\\ & \quad - \rho {\eta _{j1}}({U_j}{\varphi _{gj1,yx}} + {\varphi _{gj1,yt}}) = 0\textrm{ at }y = {( - 1)^{j + 1}}. \end{array}\end{gather}Focusing on temporal sinuous instabilities, the initial condition for both first and second order is given as
\begin{equation}\left. {\begin{array}{*{20}{c@{}}} {{\eta_{j = 1}}{|_{t = 0}}}\\ {{\eta_{j = 2}}{|_{t = 0}}} \end{array}} \right\} = {\eta _0}\cos ({k_1}x) = \tfrac{1}{2}{\eta _0}exp (\textrm{i}{k_1}x) + \textrm{c}\textrm{.c}\textrm{.},\end{equation}
where 
${k_1}$ is the dimensionless wavenumber (Ibrahim & Jog Reference Ibrahim and Jog2008).
3. Solution procedure
3.1. First-order solutions
The first-order solutions are expressed as
\begin{equation}({u_1},{v_1},{P_{l1}},{\boldsymbol{\tau
}_1},{\boldsymbol{\gamma }_1},{\varphi _{gj1}},{\eta
_{j1}}) = [{\hat{u}_1},{\hat{v}_1},{\hat{P}_{l1}},{\hat{\boldsymbol{\tau
}}_1},{\hat{\boldsymbol{\gamma }}_1},{\hat{\varphi
}_{gj1}},{\hat{\eta }_{j1}}]exp [\textrm{i(}{k_1}x -
{\omega _1}t\textrm{)}] + \textrm{c.c}.\end{equation}
Here, 
${\omega _1}$ is the linear complex frequency whose real (α) and imaginary (β) parts signify angular frequency and growth rate, respectively. The symbol ‘
$\wedge$’ represents the components that are functions of 
$y$ only and c.c. is complex conjugate.
Substituting (3.1) in the linear constitutive relation (2.17) yields
\begin{equation}{\boldsymbol{\tau }_1} = \frac{1}{{R{e_1}}}{\boldsymbol{\gamma }_1}\quad \textrm{where}\;R{e_1} = \frac{{1 + {\lambda _1}\textrm{i(}{k_1} - {\omega _1}\textrm{)}}}{{1 + {\lambda _2}\textrm{i(}{k_1} - {\omega _1}\textrm{)}}}Re.\end{equation}
Here, 
$R{e_1}$ represents first-order effective Reynolds number.
Substitution of (3.1) and (3.2) in the first-order governing equations (2.17) to (2.20) and (2.23) yields the solution of the disturbance field with a set of integration constants. The expressions for these constants are derived using the boundary conditions described in (2.21) and (2.22), and (2.24) to (2.26). The final solutions for liquid and gas phase are given as
\begin{gather}{\hat{u}_1} = \textrm{i}{A_1}\sinh ({k_1}y) + \frac{{\textrm{i}{l_1}}}{{{k_1}}}{B_1}\sinh ({l_1}y),\end{gather}
\begin{gather}{\hat{v}_1} = {A_1}\cosh ({k_1}y) + {B_1}\cosh ({l_1}y),\end{gather}
\begin{gather}{\hat{p}_1} = \frac{{({w_1} - {k_1})}}{{{k_1}}}\textrm{i}{A_1}\sinh ({k_1}y),\end{gather} where 
${A_1} = (k_1^2 + l_1^2){\hat{\eta }_j}/R{e_1}\cosh ({k_1})$, 
${B_1} ={-} 2k_1^2{\hat{\eta }_j}/R{e_1}\cosh ({l_1})$ and
\begin{gather}l_1^2 = k_1^2 + \textrm{i}R{e_1}({k_1} - {\omega _1}), \end{gather}
\begin{gather}{\hat{\varphi }_{gj1}} = {( - 1)^{j + 1}}\left[ {\frac{{\textrm{i}{\omega_1}}}{{{k_1}}} - \textrm{i}U} \right]{\hat{\eta }_j}exp [{k_1} + {( - 1)^j}{k_1}y].\end{gather}The first-order dispersion equation is obtained by substituting (3.3) to (3.7) in the dynamic boundary condition (equation (2.26)) and is expressed as
\begin{equation}- \rho [{({\omega _1} - U{k_1})^2}] + \frac{{k_1^3}}{{We}} + \frac{{{{(l_1^2 + k_1^2)}^2}}}{{R{e_1}^2}}\tanh ({k_1}) - {\left( {\frac{2}{{R{e_1}}}} \right)^2}{l_1}k_1^3\tanh ({l_1}) = 0.\end{equation}
The disturbance growth rate requires the solution for temporal frequency 
$({\omega _1})$, which is obtained by solving (3.8).
3.2. Second-order solutions
Clark & Dombrowski (Reference Clark and Dombrowski1972), Jazayeri & Li (Reference Jazayeri and Li2000) and Yang et al. (Reference Yang, Wang, Fu, Du and Tong2013) have reported that purely time-dependent functions are not present in second-order deformation. Liquid sheet disintegration under the influence of sinuous disturbances takes place due to the combined effect of first-order growth rate, first harmonic, i.e. energy transfer from fundamental first-order to second-order disturbances, and the weak influence of inherent second-order instabilities. The first harmonic and fundamental first-order disturbances share the same frequency 
$({\omega _1})$ and correspond to wavenumber 
${k_1}$, while the second-order frequency 
$({\omega _2})$ corresponds to the homogeneous second-order equations having wavenumber 
$2{k_1}$. Hence, second-order solutions are given as
\begin{equation}\begin{array}{ccccc} & (u{}_2,{v_2},{P_{l2}},{\boldsymbol{\tau }_{21}},{\boldsymbol{\gamma }_{21}},{\varphi _{gj2}},{\eta _{j2}})\\ & \quad = [{{\hat{u}}_{21}},{{\hat{v}}_{21}},{{\hat{P}}_{l21}},{{\hat{\boldsymbol{\tau }}}_{21}},{{\hat{\boldsymbol{\gamma }}}_{21}},{{\hat{\varphi }}_{gj21}},{{\hat{\eta }}_{j21}}]exp [2\textrm{i}({k_1}x - {\omega _1}t)] + \textrm{c}\textrm{.c}\textrm{.}\\ & \quad \quad + [{{\hat{u}}_{22}},{{\hat{v}}_{22}},{{\hat{P}}_{l22}},{{\hat{\boldsymbol{\tau }}}_{22}},{{\hat{\boldsymbol{\gamma }}}_{22}},{{\hat{\varphi }}_{gj22}},{{\hat{\eta }}_{j22}}]exp [\textrm{i}(2{k_1}x - {\omega _2}t)] + \textrm{c}\textrm{.c}\textrm{.} \end{array}\end{equation}
The terms with subscript ‘21’ have the same coefficient of exponent as the non-homogeneous terms obtained as a result of products of two first-order solutions. They indicate the transfer of energy from fundamental first order to second order. The homogeneous inherent disturbances having second-order complex frequency 
$({\omega _2})$ are represented with subscript ‘22’. In this regard, it needs to be noted that the first-order solutions are present in the second-order equations in the form of products each of two terms, expressed with subscript ‘21’. The second-order solutions can be obtained based upon the non-homogeneous form of first-order solutions.
The non-homogeneous terms with subscript ‘21’in the constitutive equation (2.27) are represented as
\begin{equation}{u_1}{\gamma _{1,x}} + {v_1}{\gamma _{1,y}} - {\boldsymbol{W}_1}\boldsymbol{\cdot }{\boldsymbol{\gamma }_1} + {\boldsymbol{\gamma }_1}\boldsymbol{\cdot }{\boldsymbol{W}_1} = {\boldsymbol{D}_{21c}},\end{equation}
where 
${\boldsymbol{D}_{21c}}$ is a tensor matrix consisting of four elements. Previously obtained first-order solutions are substituted in (3.10) and the terms with coefficient 
$exp 2\,\textrm{i(}{k_1}x - {\omega _1}t\textrm{)}$ are collected to obtain the expressions for these four components as
\begin{gather}{\hat{D}_{21cxx}} ={-} 2k_1^2\hat{u}_1^2 + 2i{k_1}{\hat{v}_1}{\hat{u}_{1,y}} - \hat{u}_{1,y}^2 - k_1^2\hat{v}_1^2,\end{gather}
\begin{gather}{\hat{D}_{21cxy}} = {\hat{D}_{21cyx}} = 2i{k_1}{\hat{u}_1}{\hat{u}_{1,y}} + {\hat{v}_1}{\hat{u}_{1,yy}} + 2i{k_1}{\hat{v}_1}{\hat{v}_{1,y}} - {\hat{u}_{1,y}}{\hat{v}_{1,y}},\end{gather}
\begin{gather}{\hat{D}_{21c,yy}} = 2ik{\hat{u}_1}{\hat{v}_{1,y}} + 2{\hat{v}_1}{\hat{v}_{1,yy}} + \hat{u}_{1,y}^2 + k_1^2\hat{v}_1^2.\end{gather} Now, the first-order constitutive relation (3.2) and the non-homogeneous terms expressed in (3.10) are substituted in the second-order constitutive relation (2.27) and the terms with coefficient 
$exp 2\,\textrm{i}({k_1}x - {\omega _1}t)$ are collected to obtain
\begin{equation}\begin{array}{ccccc}
{{\hat{\boldsymbol{\tau }}}_{21}} & =
\dfrac{1}{{R{e_{21}}}}{{\hat{\boldsymbol{\gamma }}}_{21}} +
\dfrac{1}{{Re}}\dfrac{{({\lambda _2} - {\lambda
_1}){{\hat{\boldsymbol{D}}}_{21c}}}}{{[1 + {\lambda
_1}(\textrm{i}k - \textrm{i}{\omega _1})][1 + 2{\lambda
_1}(\textrm{i}k - \textrm{i}{\omega _1})]}},\\ & \quad
\textrm{where}\;R{e_{21}} = \dfrac{{1 + 2{\lambda
_1}(\textrm{i}{k_1} - \textrm{i}{\omega _1})}}{{1 +
2{\lambda _2}(\textrm{i}{k_1} - \textrm{i}{\omega _1})}}Re.
\end{array}\end{equation}
Substituting equation (3.9) in continuity equation (2.28) yields
\begin{equation}{\hat{u}_{21}} = \frac{\textrm{i}}{{2{k_1}}}{\hat{v}_{21}}.\end{equation}Similarly, the solutions for non-homogeneous terms in the expression for liquid pressure and velocities are given as
\begin{gather}{\hat{P}_{l21}} = {c_1}exp (2{k_1}y) + {c_2}exp ( - 2{k_1}y) + {\hat{S}_{21p}},\end{gather}
\begin{gather}{\hat{u}_{21}} = {c_3}exp ({l_{21}}y) \!+\! {c_4}exp ( - {l_{21}}y) \!+\! R{e_{21}}\frac{{\textrm{i}k}}{{4k_1^2 - l_{21}^2}}[2{c_1}exp (2{k_1}y) \!+\! 2{c_2}exp ( - 2{k_1}y)] + {\hat{S}_{21u}},\end{gather}
\begin{gather}{\hat{v}_{21}} = {c_5}exp ({l_{21}}y) \!+\! {c_6}exp ( - {l_{21}}y) \!+\! R{e_{21}}\frac{k}{{4k_1^2 - l_{21}^2}}[2{c_1}exp (2{k_1}y) - 2{c_2}exp ( - 2{k_1}y)] \!+\! {\hat{S}_{21v}},\end{gather} where 
$l_{21}^2 = 4k_1^2 + R{e_{21}}\textrm{i}(2{k_1} - 2{\omega _1})$ and 
${c_1},\;{c_2},\;{c_3},\;{c_4},\;{c_{5}}\;\textrm{and}\;{c_6}$ are constants of integration. The detailed solution procedure and corresponding expressions of 
${\hat{S}_{21p}}$, 
${\hat{S}_{21u}}$ and 
${\hat{S}_{21v}}$ are given in Appendix A. The expression for second-order gas velocity is obtained using a procedure similar to that for first-order analysis and is given as
\begin{equation}{\hat{\varphi }_{gj21}} = \left[ {{{( - 1)}^{j + 1}}\frac{{{{\hat{\eta }}_{j21}}}}{{{k_1}}}\textrm{i}({\omega_1} - U{k_1}) + \hat{\eta }_{j1}^2\textrm{i(}{\omega_1} - U{k_1}\textrm{)}} \right]exp [{( - 1)^j}2{k_1}y].\end{equation}
Finally, the non-homogeneous terms in the boundary conditions given in (2.31), (2.32) and (2.36) having coefficient 
$exp 2\textrm{i(}{k_1}x - {\omega _1}t\textrm{)}$ are represented as 
${\hat{E}_{21\eta }}$, 
${\hat{E}_{21\tau }}$ and 
${\hat{E}_{21d}}$. Corresponding expressions are given in Appendix A. Thus, the second-order boundary conditions given in (2.31), (2.32) and (2.36) can be represented in terms of the expression for non-homogeneous terms as
\begin{gather}{\hat{v}_{21}} - \textrm{i(}2{k_1} - {\omega _1}\textrm{)}{\hat{\eta }_{j21}} = {\hat{E}_{21\eta }},\end{gather}
\begin{gather}{\hat{\tau }_{yx21}} ={-} {\hat{E}_{21\tau }},\end{gather}
\begin{gather}- {\hat{p}_{l21}} + {\hat{\tau }_{yy21}} + 2i\rho ({\omega _1}{\hat{\varphi }_{gj21,t}} - U{k_1}{\hat{\varphi }_{gj21,xy}}) + \frac{1}{{We}}{( - 1)^{j + 1}}4k_1^2{\hat{\eta }_{j21}} ={-} {\hat{E}_{21d}}.\end{gather} With this, one arrives at the solution for liquid and gas pressure and velocity field, as well as the boundary conditions based upon the first-order solutions. The interface displacement due to energy transfer from first to second order requires the solution for the first harmonic 
$({\hat{\eta }_{j21}})$. To obtain the same, the expressions for constitutive relation, x and y components of liquid velocity, liquid pressure and gas velocity as given in (3.14), (3.15), (3.16), (3.18) and (3.19) are substituted in (3.20) to (3.22). This yields six equations that can be solved algebraically to obtain the expression for the interface displacement due to the first harmonic 
${\hat{\eta }_{j21}}$. However, the solution procedure can be further simplified following Yang et al. (Reference Yang, Wang, Fu, Du and Tong2013). Here, 
${\hat{E}_{21\eta }}$, 
${\hat{E}_{21\tau }}$, 
${\hat{D}_{21xy}}\textrm{,}\;{\hat{D}_{21yx}}$ and 
${\hat{S}_{21v}}$ are odd as they are hyperbolic functions of sine, whereas 
${\hat{D}_{c21xx}},\;{\hat{D}_{c21yy}},\;{\hat{E}_{21d}},\;{\hat{S}_{21u}}$ and 
${\hat{S}_{21p}}$ are even as they are hyperbolic functions of cosine. A parity analysis between the different terms shows that
\begin{equation}{c_1} = {c_2},\;{c_5} ={-} {c_6}\quad \textrm{and}\quad {\eta _{j21}}{|_{j = 1}}\textrm{ = } - {\eta _{j21}}{|_{j = 2}}.\end{equation}
Equation (3.23) is substituted in (3.20)–(3.22) and the set of algebraic equations are solved. Finally, the expression for 
${\hat{\eta }_{j21}}$ is obtained as
\begin{align}&{\hat{\eta }_{j21}} =-{{( -
1)}^{j + 1}}\nonumber\\ &\quad {\left. {\left\{
\begin{array}{l} \dfrac{{(l_{21}^2 + 4k_1^2)\coth
(2{k_1})}}{{2R{e_{21}}(4k_1^2 - l_{21}^2)}}\left[
\begin{array}{l} (l_{21}^2 + 4k_1^2){{\hat{E}}_{21\eta }} -
2\textrm{i}{k_1}R{e_{21}}{{\hat{E}}_{21\tau }} +
R{e_{21}}{{\hat{S}}_{21p,y}}\\ \quad + \,
{W_1}\left( { - 4\textrm{i}{k_1}{{\hat{D}}_{21cxy}} -
\dfrac{{\partial {{\hat{D}}_{21cyy}}}}{{\partial y}}}
\right) \end{array} \right]\\ \quad + \,
\dfrac{{(4\textrm{i}{k_1}{l_{21}})\coth
({l_{21}})}}{{R{e_{21}}(4k_1^2 - l_{21}^2)}}\left[
\begin{array}{l}
\dfrac{\textrm{i}}{{2{k_1}}}(R{e_{21}}{{\hat{S}}_{21p,y}})
- {W_1}\left( { - 2{{\hat{D}}_{21cxy}} +
\dfrac{\textrm{i}}{{2{k_1}}}\dfrac{{\partial
{{\hat{D}}_{21cyy}}}}{{\partial y}}}
\right)\\ \quad + \,
R{e_{21}}{{\hat{E}}_{21\tau }} +
4\textrm{i}{k_1}{{\hat{E}}_{21\eta }} -
\dfrac{{\textrm{i}(4k_1^2 -
l_{21}^2)}}{{2{k_1}}}{{\hat{S}}_{21v}} \end{array}
\right]\\ \quad +\, {{\hat{S}}_{21p}} -
\dfrac{2}{{R{e_{21}}}}{{\hat{S}}_{21v,y}} + 2{({\omega_1} -
U)^2}\rho \eta_{j1}^2 - {{\hat{E}}_{21d}} -
{W_2}{{\hat{D}}_{21cyy}} \end{array} \right\}} \right|_{y =
1}}\frac{{2k}}{{{D_{21{\mathop{\rm var}}
}}}},\end{align}
where
\begin{align}{W_1} &= \frac{{{\lambda _2} - {\lambda _1}}}{{[1 + {\lambda _1}(\textrm{i}{k_1} - \textrm{i}{\omega _1})][1 + 2{\lambda _2}(\textrm{i}{k_1} - \textrm{i}{\omega _1})]}}\quad \textrm{and}\nonumber\\ {W_2} &= \frac{{{\lambda _2} - {\lambda _1}}}{{[1 + {\lambda _1}(\textrm{i}{k_1} - \textrm{i}{\omega _1})][1 + 2{\lambda _1}(\textrm{i}{k_1} - \textrm{i}{\omega _1})]}}.\end{align} The term 
${D_{21{\mathop{\rm var}} }}$ represents the first harmonic dispersion equation and is observed to be varicose in nature. Its expression is given as
\begin{equation}{D_{21{\mathop{\rm var}} }} ={-} \rho [{(2{\omega _1} - 2U{k_1})^2}] + \frac{{8k_1^3}}{{We}} + \frac{{{{(l_{21}^2 + 4k_1^2)}^2}}}{{Re_{21}^2}}\coth ({k_1}) - {\left( {\frac{2}{{Re_{21}^2}}} \right)^2}8{l_{21}}k_1^3\coth ({l_{21}}).\end{equation} Obtaining the expression for interface deformation is enough to get the position of perturbed surfaces at different instances of time. As a result, the expressions for 
${c_1},\;{c_2},\;{c_3}$ and 
${c_4}$ are not noted. Next, the inherent second-order constitutive relation is obtained by substituting (3.9) in (2.26) and collecting terms with coefficient 
$exp [\textrm{i(}2{k_1}x - {\omega _2}t\textrm{)}]$, as shown in (3.27):
\begin{equation}{\hat{\boldsymbol{\tau }}_{22}} = \frac{1}{{R{e_{22}}}}{\hat{\boldsymbol{\gamma }}_{22}},\end{equation}
where 
$R{e_{22}}$ is the second-order effective Reynolds number expressed as
\begin{equation}R{e_{22}} = \frac{{1 + {\lambda _1}(2\,\textrm{i}{k_1} - {\omega _2})}}{{1 + {\lambda _2}({2\,\textrm{i}{k_1} - {\omega_2}} )}}Re.\end{equation} The solution for inherent second-order disturbances 
${\hat{\eta }_{j22}}$ can similarly be obtained by substituting the terms with coefficient 
$exp [\textrm{i}(2{k_1}x - {\omega _2}t)]$, as given in (3.9), in the governing equations and boundary conditions (2.27)–(2.36). However, the need for the same can be eliminated by substituting equations (2.16), (3.1) and (3.9) in equation (2.37), which yields
\begin{equation}{\hat{\eta }_{j1}} = {\textstyle{1 \over 2}};\quad {\hat{\eta }_{j21}} ={-} {\hat{\eta }_{j22}}.\end{equation}From (3.23) and (3.29), it can be deduced that
\begin{equation}{\hat{\eta }_{j22}}{|_{j = 1}} ={-} {\hat{\eta }_{j22}}{|_{j = 2}}.\end{equation}Equation (3.30) suggests that the second-order disturbance is varicose in nature. Also, the absence of non-homogeneous terms in inherent second-order solutions allows them to be expressed in a form similar to the first-order solution. As a result, the non-trivial solution of the second-order dispersion equation is given as
\begin{equation}- \rho [{({\omega _2} - 2U{k_1})^2}] + \frac{{8k_1^3}}{{We}} + \frac{{{{(l_1^2 + 4k_1^2)}^2}}}{{Re_{23}^2}}\coth (2{k_1}) - {\left( {\frac{2}{{R{e_{23}}}}} \right)^2}8{l_1}k_1^3\coth ({l_1}) = 0.\end{equation}The expression for final interface displacement is obtained from (2.16) and (3.9):
\begin{equation}{\eta _j} = \frac{1}{2}{\eta _0}exp [\textrm{i}({k_1}x - {\omega _1}t)] + {\left. {\eta_0^2\left\{ \begin{array}{l} {{\hat{\eta }}_{j21}}exp [{2\textrm{i}({k_1}x - {\omega_1}t)} ]\\ - {{\hat{\eta }}_{j21}}exp [\textrm{i}(2{k_1}x - {\omega_2}t)] \end{array} \right\}} \right|_{{{\hat{\eta }}_{j1}} = 1/2}} + \textrm{c}\textrm{.c}\textrm{.}\end{equation}4. Maximum entropy formulation for droplet size and velocity distribution using modified Newton–Raphson method
The nonlinear instability model provides breakup length and the ligament geometry following breakup. The cross-sectional area of the ligament generated after the breakup of the liquid sheet is calculated from the shape of the deformed interface and the area of an equivalent cylindrical liquid column is evaluated. The mass mean diameter of the droplets (D 30) is calculated from the generated ligament as per the Rayleigh instability of the cylindrical liquid column (Rayleigh Reference Rayleigh1878). In the second part of the comprehensive model, an MEF-based model has been developed to get the most probable droplet size distribution by maximizing the Shannon entropy, subject to the respective constraint conditions which in turn are obtained using the mean droplet diameter and the breakup length. The constraint equations are formulated based on the conservation of physical quantities, like mass, momentum and energy, and by considering the physics of the breakup process. Accordingly, the momentum and energy exchange during breakup have been accounted for. The present study closely follows the non-empirical model developed by Nath et al. (Reference Nath, Datta, Mukhopadhyay, Sen and Tharakan2011) in this regard.
Four constraint conditions, normalized, mass balance, momentum balance and energy balance conditions, have been considered, which are given as
\begin{gather}\textrm{normalization}:\quad \sum\limits_i^m {\sum\limits_j^n {{P_{i,j}}} } = {S_1},\end{gather}
\begin{gather}\textrm{mass}\;\textrm{balance:}\quad \sum\limits_i^m {\sum\limits_j^n {{P_{i,j}}D_i^3} } = {S_2},\end{gather}
\begin{gather}\textrm{momentum}\;\textrm{balance:}\quad \sum\limits_i^m {\sum\limits_j^n {{P_{i,j}}D_i^3}V_j} = {S_3},\end{gather}
\begin{gather}\textrm{energy}\;\textrm{balance:}\quad \sum\limits_i^m {\sum\limits_j^n {{P_{i,j}}D_i^3V_j^2} } + B\sum\limits_i^m {\sum\limits_j^n {{P_{i,j}}D_i^2} } = {S_4}.\end{gather} In (4.1) to (4.4), 
${S_1}$, 
${S_2}$, 
$S{}_3$ and 
${S_4}$ are source terms such that 
${S_1} = 1$; 
${S_2} = 1$; 
${S_3} = 1 + {\textstyle{1 \over 2}}\rho {(U - 1)^2}L{C_f}$; 
${S_4} = 1 + \rho {(U - 1)^3}L{C_f}$; L = non-dimensional length 
$(1/h)$; ρ = gas to liquid density ratio 
$({\rho _g}/{\rho _l})$; U = non-dimensional gas velocity 
$({u_g}/{u_l})$; Di = non-dimensional droplet diameter 
$({d_i}/{d_m})$; Vj = non-dimensional droplet velocity 
$({v_j}/{u_l})$; Pi,j = joint probability function 
$({n_{i,j}}/N)$; 
$B = 12\sigma /{\rho _l}u_l^2{d_m}$; σ = surface tension coefficient for liquid; and N represents the total number of droplets produced per unit time.
There are possibilities of an infinite probability distribution function 
$({P_{i,j}})$ which may satisfy the available constraint conditions. However, the objective is to find the distribution function that maximizes the Shannon entropy. The expression for the Shannon entropy is given as
\begin{equation}S ={-} \kappa \sum\limits_i^m {\sum\limits_j^n {{P_{i,j}}\ln {P_{i,j}}} } ,\end{equation}where κ is the Boltzmann constant. Expression for the probability distribution function after maximizing the Shannon entropy is obtained as
\begin{equation}{P_{i,j}} = exp [ - (1 + \overline{\lambda _1} + \overline{\lambda _2}D_i^3 + \overline{\lambda _3}D_i^3{V_j} + \overline{\lambda _4}(D_i^3V_j^2 + BD_i^2))].\end{equation} Substituting the expression for 
${P_{i,j}}$ in (4.1) to (4.4), the following set of nonlinear algebraic equations are obtained, which are solved numerically to obtain the values of Lagrange multipliers:
\begin{gather}f = \sum\limits_{i = 1}^m {\sum\limits_{j = 1}^k {exp [ - (1 + \overline{\lambda _1} + \overline{\lambda _2}D_i^3 + \overline{\lambda _3}D_i^3{V_j} + \overline{\lambda _4}(D_i^3V_j^2 + BD_i^2))]} } - {S_1} = 0,\end{gather}
\begin{gather}g = \sum\limits_{i = 1}^m {\sum\limits_{j = \textrm{1}}^k {exp [ - (1 + \overline{\lambda _1} + \overline{\lambda _2}D_i^3 + \overline{\lambda _3}D_i^3{V_j} + \overline{\lambda _4}(D_i^3V_j^2 + BD_i^2))]D_i^3} - {S_2}} = 0,\end{gather}
\begin{gather}h = \sum\limits_{i = 1}^m {\sum\limits_{j = \textrm{1}}^k {exp [ - (1 + \overline{\lambda _1} + \overline{\lambda _2}D_i^3 + \overline{\lambda _3}D_i^3{V_j} + \overline{\lambda _4}(D_i^3V_j^2 + BD_i^2))]D_i^3{V_j}} - {S_3} = 0} ,\end{gather}
\begin{gather}m = \sum\limits_{i = 1}^m {\sum\limits_{j = \textrm{1}}^k {exp [ - (1 \!+\! \overline{\lambda _1} \!+\! \overline{\lambda _2}D_i^3 \!+\! \overline{\lambda _3}D_i^3{V_j} \!+\! \overline{\lambda _4}(D_i^3V_j^2 \!+\! BD_i^2))]} (D_i^3{V^2} \!+\! BD_i^2) - {S_4} = 0} .\end{gather}A more detailed discussion on how the constraint conditions are derived can be found in the work of Nath et al. (Reference Nath, Datta, Mukhopadhyay, Sen and Tharakan2011). The Newton–Raphson method is employed to solve the above set of nonlinear equations. However, the solution has a propensity to diverge if the initial guess is not close to the roots. Moreover, the presence of exponential terms in the equations makes it highly sensitive to the initial guesses (Li & Tankin Reference Li and Tankin1988; Chin et al. Reference Chin, Larose, Tankin, Jackson, Stutrud and Switzer1991; Mondal, Datta & Sarkar Reference Mondal, Datta and Sarkar2003). To tackle the problem of divergence of solutions, the present paper considers a Taylor series expansion up to second order. It has been observed that inclusion of higher-order terms reduces the strict requirement of a close initial guess (Li & Li Reference Li and Li2006).
Each of the nonlinear equations is expanded in the Taylor series around the solution up to the second-order terms. For instance, the expanded normalization condition is given as
\begin{gather}f(\overline{\lambda _1^0} + \alpha ,\overline{\lambda _2^0} + \beta ,\overline{\lambda _3^0} + \gamma ,\overline{\lambda _4^0} + \varphi ) = 0,\end{gather}
\begin{gather}{f_0} + \alpha {\left. {\frac{{\partial f}}{{\partial \overline{\lambda_1}}}} \right|_0} + \alpha {\left. {\frac{{{\partial^2}f}}{{\partial \overline{\lambda_1^2}}}} \right|_0} + \beta {\left. {\frac{{\partial f}}{{\partial \overline{\lambda_2}}}} \right|_0} + \beta {\left. {\frac{{{\partial^2}f}}{{\partial \overline{\lambda_2^2}}}} \right|_0} + \gamma {\left. {\frac{{\partial f}}{{\partial \overline{\lambda_3}}}} \right|_0}\nonumber\\ + \gamma {\left. {\frac{{{\partial^2}f}}{{\partial \overline{\lambda_3^2}}}} \right|_0} + \varphi {\left. {\frac{{\partial f}}{{\partial \overline{\lambda_4}}}} \right|_0} + \varphi {\left. {\frac{{{\partial^2}f}}{{\partial \overline{\lambda_4^2}}}} \right|_0} = 0.\end{gather}The expanded normalization constraint condition (4.12) is represented in matrix form as
\begin{align}&\left[
{\begin{array}{*{20}{c}} {\dfrac{{\partial f}}{{\partial
\overline{\lambda_1}}}}&{\dfrac{{\partial f}}{{\partial
\overline{\lambda_2}}}}&{\dfrac{{\partial f}}{{\partial
\overline{\lambda_3}}}}&{\dfrac{{\partial f}}{{\partial
\overline{\lambda_4}}}} \end{array}} \right]\left[
{\begin{array}{*{20}{c}} \alpha \\ \beta \\ \gamma \\
\varphi \end{array}} \right] + \frac{1}{2}\left[
{\begin{array}{*{20}{c}} \alpha &\beta &\gamma &\varphi
\end{array}} \right]\nonumber\\\quad &\left[
{\begin{array}{*{20}{c}} {\dfrac{{{\partial^2}f}}{{\partial
\overline{\lambda_1^2}}}}&{\dfrac{{{\partial^2}f}}{{\partial
\overline{\lambda_1}\partial
\overline{\lambda_2}}}}&{\dfrac{{{\partial^2}f}}{{\partial
\overline{\lambda_1}\partial
\overline{\lambda_3}}}}&{\dfrac{{{\partial^2}f}}{{\partial
\overline{\lambda_1}\partial \overline{\lambda_4}}}}\\
{\dfrac{{{\partial^2}f}}{{\partial \overline{\lambda_2}\partial
\overline{\lambda_1}}}}&{\dfrac{{{\partial^2}f}}{{\partial
\overline{\lambda_2^2}}}}&{\dfrac{{{\partial^2}f}}{{\partial
\overline{\lambda_2}\partial
\overline{\lambda_3}}}}&{\dfrac{{{\partial^2}f}}{{\partial
\overline{\lambda_2}\partial \overline{\lambda_4}}}}\\
{\dfrac{{{\partial^2}f}}{{\partial \overline{\lambda_3}\partial
\overline{\lambda_1}}}}&{\dfrac{{{\partial^2}f}}{{\partial
\overline{\lambda_3}\partial
\overline{\lambda_2}}}}&{\dfrac{{{\partial^2}f}}{{\partial
\overline{\lambda_3^3}}}}&{\dfrac{{{\partial^2}f}}{{\partial
\overline{\lambda_3}\partial \overline{\lambda_4}}}}\\
{\dfrac{{{\partial^2}f}}{{\partial \overline{\lambda_4}\partial
\overline{\lambda_1}}}}&{\dfrac{{{\partial^2}f}}{{\partial
\overline{\lambda_4}\partial
\overline{\lambda_2}}}}&{\dfrac{{{\partial^2}f}}{{\partial
\overline{\lambda_4}\partial
\overline{\lambda_3}}}}&{\dfrac{{{\partial^2}f}}{{\partial
\overline{\lambda_4^4}}}} \end{array}} \right]\left[
{\begin{array}{*{20}{c}} \alpha \\ \beta \\ \gamma \\
\varphi \end{array}} \right] ={-}
{f_0}.\end{align}
Mass, momentum and energy constraint equations are also similarly expanded and represented in matrix form. However, they are not shown here for brevity. Since the above expression is nonlinear in nature, solution of the Lagrange multipliers requires the linearization of expressions; the simplest way is to replace the unknown values of 
$\alpha ,\beta ,\gamma $ and 
$\phi $ associated with the second-order term with the known values from previous iteration using Picard's method.
The vector consisting of first-order derivative terms and the matrix having the second-order derivative terms on the left-hand side of (4.13) are represented as F 1 and F 2, respectively. Similarly, for the mass, momentum and energy conservation constraints, first- and second-order terms are given as G 1 and G 2, H 1 and H 2, and M 1 and M 2, respectively. The solution matrix for all the nonlinear equations is obtained as
\begin{equation}\left[ {\begin{array}{*{20}{c}} {{F_1} + {F_2}}\\ {{G_1} + {G_2}}\\ {{H_1} + {H_2}}\\ {{M_1} + {M_2}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} \alpha \\ \beta \\ \gamma \\ \varphi \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} { - {f_o}}\\ { - {g_o}}\\ { - {h_o}}\\ { - {m_o}} \end{array}} \right].\end{equation}Considering the complexity of the nonlinear equations and to maintain the stability of the system, normalization constraint has been checked using expressions shown in (4.15) at the end of each iteration as suggested by Sellens & Brzustowski (Reference Sellens and Brzustowski1985):
\begin{equation}\overline{\lambda _1} = \ln \left[ {\sum\limits_{i = 1}^m {\sum\limits_{j = 1}^k { - (1 + \overline{\lambda_2}D_i^3 + \overline{\lambda_3}D_i^3{V_j} + \overline{\lambda_4}(D_i^3V_j^2 + BD_i^2))} } } \right].\end{equation} The converged solutions of 
$\overline{\lambda _1},\overline{\lambda _2},\overline{\lambda _3}$ and 
$\overline{\lambda _4}$ are used in (4.6) to obtain the probability distribution function 
${P_{i,j}}$ at different droplet diameters 
$({D_i})$ and velocities 
$({V_j})$.
Finally, the volumetric probability density for droplet size is expressed as
\begin{equation}\partial {Q_i}/\partial {d_i} = ({Q_i} - {Q_{i - 1}})/({d_i} - {d_{i - 1}}),\end{equation}
where 
${Q_i} = \sum\nolimits_{c = 1}^i {\sum\nolimits_{j = 1}^k {{P_{c,j}}D_c^3} }$ represents fraction of liquid having diameter less than 
${d_i}$.
5. Results and discussion
Stability of a viscoelastic planar liquid sheet has been studied in the past, and stabilizing and destabilizing effects of time constant ratio and elasticity number, respectively, have been established. However, a Kelvin–Helmholtz type of instability, as considered in the current study, is primarily driven by a velocity jump across the liquid–gas interfaces. So, the natural question to ask is whether the effects of viscoelasticity on the liquid sheet vary with the surrounding gas velocities. To address the same, this section features linear and nonlinear investigations of the effects of elasticity number (El) and time constant ratio 
$(\lambda )$ at different gas to liquid velocity ratios (U). Following the experimental investigation of viscoelastic jets by Carroll & Joo (Reference Carroll and Joo2006), polymer density is kept at 850 kg m−3. Gas density is considered to be 1.2 kg m−3 and, subsequently, gas to liquid density ratio (ρ) becomes 0.0014. Literature reports that when liquid Weber number (We) is varied between 200 and 500, linear analysis shows good agreement with experimental results due to the relatively weak effect of nonlinearity at medium Weber number. Yang et al. (Reference Yang, Wang, Fu, Du and Tong2013) also recommended that weakly nonlinear instability analysis should be applied to cases where 
$\rho We < 1$, such that disturbance amplitudes of orders higher than 2 are relatively very small. As a result, liquid Weber number is chosen to be 300 for the present study. Keeping consistency with Wang et al. (Reference Wang, Yang, Xie and Chen2015), liquid Reynolds number is kept at 63.54. Validation of the present model is done with the temporal instability study of viscoelastic fluid by Wang et al. (Reference Wang, Yang, Xie and Chen2015) and is presented in table 1.
Comparison of breakup time obtained in the present study with that of Wang et al. (Reference Wang, Yang, Xie and Chen2015) for different elasticity number (El) keeping We = 600, Re = 63.54, ρ = 0.0012, λ = 0 and k = 0.46.
Table 1 shows the effect of elasticity number on breakup time predicted by the present study and it shows excellent agreement with the results of Wang et al. (Reference Wang, Yang, Xie and Chen2015). For the purpose of validation, the wavenumber is kept fixed for all values of El. However, in the subsequent discussion of results, the dominant wavenumber for calculating breakup time is obtained from maximum linear growth rate, and hence may change with fluid properties.
5.1. First-order results
Solution of the first-order dispersion relation, as given in (3.8), provides the first-order complex root 
$({\omega _1})$. The imaginary part (β) of the complex root represents disturbance growth rate. Variation of β with wavenumber 
$({k_1})$ provides the maximum growth rate 
$({\beta _{max}})$. The wavenumber corresponding to 
${\beta _{max}}$ is considered as the most dominating wavenumber 
$({k_{dominant}})$ since it registers the highest growth rate (Rayleigh Reference Rayleigh1878).
Figures 2(a) and 2(b) show variation of maximum growth rate 
$({\beta _{max}})$ and corresponding 
${k_{dominant}}$ with elasticity number (El), respectively, for different gas to liquid velocity ratios (U) keeping We = 300 and ρ = 0.0014.
Variation of (a) maximum growth rate and (b) dominant wavenumber with El at different U keeping 
$\lambda = 0$, We = 300 and ρ = 0.0014.
An increase in elasticity number increases the stress relaxation time 
${\lambda _1}$ and consequently produces a larger effective Reynolds number. An effective Reynolds number is nothing but a similarity parameter to link the behaviour exhibited by a non-Newtonian fluid to that of a Newtonian fluid. An increase in effective Reynolds number enhances the linear growth rate. As a result, 
${\beta _{max}}$ increases with El at all gas to liquid velocity ratios (U). However, the overall effect of elasticity on growth rate is small. In general, the viscous component of viscoelastic fluid tends to weaken the instability by decreasing the disturbance growth rate, whereas the liquid elasticity counterpart promotes enhancement of instability. Theoretical investigation of viscoelastic fluid has shown that under the combined action of liquid viscosity and elasticity effects, the growth rate curve of non-Newtonian liquid sheets lies between those of inviscid and Newtonian sheets (Liu & Durst Reference Liu and Durst1998; Brenn et al. Reference Brenn, Liu and Durst2000). The area between the growth rate curves of the inviscid and non-Newtonian liquid sheets induced by the interaction of the liquid viscosity and elasticity effects is known as the viscoelasticity-induced region, whereas the area between the growth rate curves of Newtonian and non-Newtonian liquid sheets is generally termed the elasticity-enhanced region. The potential increase in growth rate with increase in El is restricted by the difference in area between the elasticity-enhanced region and viscoelasticity-induced region, which itself is small. This is the reason why even a significant increase in elasticity number does not lead to a correspondingly strong effect on growth rate. However, the effect of El on 
${\beta _{max}}$ is very minor at low gas to liquid velocity ratios (U = 2, 2.15) and increases slightly with increase in U (U = 2.50, 2.75). It should also be noted that the effect of elasticity is stronger at low to intermediate El, as 
${\beta _{max}}$ almost becomes insensitive to El at high values. Figure 2(b) shows that similar to 
${\beta _{max}}$, variation of the most unstable wavenumber 
$({k_{dominant}})$ with El at low U (U = 2, 2.15) is almost negligible. But as U increases (U = 2.50, 2.75), an increase in El shifts 
${k_{dominant}}$ towards higher values. It is well established that an increase in gas velocity broadens the range of unstable wavenumber. This allows an increase in elasticity number to displace the dominating wavenumber towards higher frequency, thereby producing more distinct amplified growth rate, as observed in figure 2(a).
The most dominant wavenumber from linear analysis is chosen as the characteristic wavenumber for the second-order analysis (Jazayeri & Li Reference Jazayeri and Li2000; Nath et al. Reference Nath, Mukhopadhyay, Sen and Tharakan2010). The interface position at different instants of time is obtained by substituting equation (3.32) in equation (2.6). The total time from the moment the sheet leaves the nozzle until the distance between the two interfaces becomes negligible is considered as the breakup time. Figure 3 shows variation of breakup time with El at different U keeping We = 300 and ρ = 0.0014. At larger velocity ratios (U = 2.50, 2.75), an initial increase in El reduces breakup time, thereby indicating a destabilizing effect of elasticity at low to intermediate El. This effect lasts for El < 20, beyond which breakup time almost becomes saturated and shows negligible variation with change in El. However, at low velocity ratios (U = 2 and 2.15), an initial increase in El stabilizes the sheet and delays breakup time until breakup time reaches a maximum value.
Variation of breakup time with El at different U keeping 
$\lambda = 0$, We = 300 and ρ = 0.0014.
An increase in El above the critical point reduces breakup time, which indicates the presence of an elasticity-induced destabilizing regime. It must be noted that the value of El corresponding to the critical point that marks the dual effect of elasticity is not constant and depends upon the gas to liquid velocity ratio(El = 4 for U = 2; El = 1 for U = 2.15).
To provide an explanation for this non-monotonic effect of elasticity number on sheet breakup, previous theoretical studies of viscoelastic fluids are considered. It is revealed that the complex dual effect of elasticity as observed in the present work has also been reported in the past. Atalık & Keunings (Reference Atalık and Keunings2002) studied temporal instabilities in viscoelastic plane channel flows. The study employed a spectral method to investigate two-dimensional temporally evolving disturbances in Poiseuille and Couette flows of viscoelastic fluids. In the case of Poiseuille flow of Oldroyd B fluid, wave amplitude initially increased with increase in El until it reached a critical value, thereby exhibiting a destabilizing effect. A further increase in El reduced the disturbance amplitude, which was referred to as a re-stabilizing zone. Nonlinear investigation of a viscoelastic liquid sheet in quiescent gas medium by Wang et al. (Reference Wang, Yang, Xie and Chen2015) reported only a destabilizing effect of elasticity. However, they considered a high Weber number (We = 600) for their study. Considering surface tension to be constant, a high Weber number represents a more significant effect of inertia forces, which is similar to the effect created by high gas to liquid velocity ratio. Hence, the observation made at U = 2.50 and 2.75 agrees with the study of Wang et al. (Reference Wang, Yang, Xie and Chen2015). Among past literature, Clark & Dombrowski (Reference Clark and Dombrowski1972) and Jazayeri & Li (Reference Jazayeri and Li2000) revealed that 
${\eta _{j1}}$ alone cannot cause breakup of a liquid sheet subjected to initial sinuous disturbances as the distance between the two interfaces in the y direction remains unchanged. As a result, the interfaces are displaced in the same direction by exactly the same amount. For an initially sinuous surface disturbance, the thinning and subsequent breakup of the liquid sheet take place only due to nonlinear effects with the generation of higher harmonics. The second-order surface disturbances in the works of Clark & Dombrowski (Reference Clark and Dombrowski1972) and Jazayeri & Li (Reference Jazayeri and Li2000) are presented as
\begin{equation}\begin{array}{ccccc} {\eta _{j2}} & = \textrm{exp(}2\textrm{i}{k_1}x\textrm{)[}{{\hat{\eta }}_{j21}}exp ( - 2\textrm{i}{\omega _1}t) + {{\hat{\eta }}_{j22}}exp ( - 2\textrm{i}{{\bar{\omega }}_1}t) + {{\hat{\eta }}_{j23}}exp (\textrm{i}{\omega _1}t - \textrm{i}{{\bar{\omega }}_1}t)\textrm{]}\\ & \quad + {{\hat{\eta }}_{j24}}exp [ - \textrm{i}({\omega _2}t)] + {{\hat{\eta }}_{j25}}exp [ - \textrm{i}({{\bar{\omega }}_2}t)], \end{array}\end{equation}
where 
${\hat{\eta }_{j21}}$, 
${\hat{\eta }_{j22}}$ and 
${\hat{\eta }_{j23}}$ represent the energy transfer from first to second order, and 
${\hat{\eta }_{j24}}$ and 
${\hat{\eta }_{j25}}$ are the inherent second-order disturbances. A comparison with the present study shows that 
${\hat{\eta }_{j21}}$ and 
${\hat{\eta }_{j24}}$ correspond to 
${\hat{\eta }_{j21}}$ and 
${\hat{\eta }_{j22}}$ of the current study. In the present paper, only the positive roots of the dispersion equations are considered. Hence, 
${\eta _{j2}}$ in the present study does not contain terms such as 
${\hat{\eta }_{j22}}$, 
${\hat{\eta }_{j23}}$ and 
${\hat{\eta }_{j25}}$. Also, the growth rates 
${\bar{\omega }_1}$, 
$({\omega _1} - {\bar{\omega }_1})$ and 
${\bar{\omega }_2}$ corresponding to the disturbance components 
${\hat{\eta }_{j22}}$, 
${\hat{\eta }_{j23}}$ and 
${\hat{\eta }_{j25}}$, respectively, are either negative or too small. Hence, their contribution to sheet instability is negligible. Lastly, 
${\omega _2}$ is obtained by solving the second-order dispersion equation which is varicose in nature, and hence it is inherently small. Also, 
${\omega _2}$ corresponds to 
$2{k_1}$, which is close to or more than the cut-off wavenumber. As a result, the disturbance component 
${\hat{\eta }_{j24}}$ that contains 
${\omega _2}$ also does not contribute to sheet instability.
In view of the above discussion, it can be concluded that 
${\hat{\eta }_{j21}}$ is the only component of second-order disturbance that contributes to sheet thinning and subsequent breakup. The effect of any fluid property on 
${\hat{\eta }_{j21}}$ can provide a possible explanation as to how it affects the overall sheet instability. As a result, the present study attempts to explain the dual effect of rheological properties of the fluid by studying their corresponding effects on 
${\hat{\eta }_{j21}}$. Since 
${\hat{\eta }_{j21}}$ is in proportion to the square of 
${\eta _{j1}}$, 
${\hat{\eta }_{j2}}$ expressed as 
${\hat{\eta }_2} = \textrm{real}({\hat{\eta }_{j21}}/\hat{\eta }_{j1}^2)$ is used to represent the influence of the first harmonic. The expression for 
${\hat{\eta }_{j21}}$ is given by (3.24). Similar approach was also adopted by Yang et al. (Reference Yang, Wang, Fu, Du and Tong2013, Reference Yang, Chen and Wang2014) and Wang et al. (Reference Wang, Yang, Xie and Chen2015) for temporal investigation of sinuous disturbances in viscous and viscoelastic fluids, respectively.
Figure 4 shows the effect of El on second-order amplitude 
$({\hat{\eta }_2})$ at different gas to liquid velocity ratios (U) keeping 
$\lambda = 0$, We = 300 and ρ = 0.0014. A comparison between figures 2 and 4 shows that the effect of elasticity on 
${\hat{\eta }_2}$ is much stronger than its effect on linear growth rate. Hence the effect of second-order amplitude associated with nonlinearity dominates the breakup process.
Variation of second-order amplitude with El at different U keeping 
$\lambda = 0$, We = 300 and ρ = 0.0014.
It can be observed that 
${\hat{\eta }_2}$ shows an upward trend with increasing El for all values of U. However, for low gas to liquid velocity ratios (U = 2 and 2.15), 
${\hat{\eta }_2}$ is initially negative for 0.1 < El < 4 at U = 2 and 0.1 < El < 1 at U = 2.15. In this range of El, an increase in elasticity increases 
${\hat{\eta }_2}$ but causes a reduction in its absolute value 
$(|{\hat{\eta }_2}|)$. Since 
${\hat{\eta }_{j21}}$ is directly proportional to 
${\hat{\eta }_{j2}}$, a decrease in 
${\hat{\eta }_{j2}}$ implies a reduction in 
${\hat{\eta }_{j21}}$ as well. Hence at this initial range of El, an increase in El dampens the very component of disturbance that primarily drives sheet breakup, thereby leading to an increase in breakup time. The stabilizing effect of El persists until 
${\hat{\eta }_2}$ finally achieves a positive value. The singularities in figure 3 represent those particular values of El at which 
${\hat{\eta }_2}$ becomes positive (El = 4 for U = 2 and El = 1 for U = 2.15). In the range of positive values of 
${\hat{\eta }_2}$, an increase in El only causes an increase in 
${\hat{\eta }_2}$, which represents an increase in 
${\hat{\eta }_{j21}}$ as well. As a result, beyond El = 4 and El = 1 for U = 2 and 2.15, respectively, an increase in elasticity produces a destabilizing effect on the liquid sheet and reduces breakup time. In the case of higher gas to liquid velocity ratio (U = 2.50 and 2.75), 
${\hat{\eta }_2}$ is positive in the entire range of El. An increase in El only enhances 
${\hat{\eta }_2}$, before it becomes saturated at its high values. Thus, elasticity promotes nonlinearity and destabilizes the sheet, ultimately resulting in a reduction in breakup time for U = 2.50 and 2.75. A comparison of different gas velocity ratios also shows that 
${\hat{\eta }_2}$ increases with an increase in U and the effect of elasticity on 
${\hat{\eta }_2}$ is more amplified and distinct in the case of higher velocity ratios (U = 2.50 and 2.75).
In view of the above discussion involving linear and nonlinear findings, it can be concluded that an increase in elasticity number monotonically increases the first-order effective Reynolds number and hence enhances the linear growth rate at all gas to liquid velocity ratios. However, sheet breakup under the influence of initial sinuous waves essentially takes place only when the first harmonic is added to the fundamental. The effect of first harmonic is expressed using second-order amplitude. Elasticity produces a much stronger effect on second-order amplitude as compared with linear growth rate. It has been observed that the relation between elasticity and second-order amplitude is non-monotonic and depends upon the range of elasticity number and gas to liquid velocity ratio. This highly coupled non-monotonic interaction between elasticity number and gas to liquid velocity ratio produces a complex effect on the second-order amplitude, and is responsible for the dual effect of elasticity on sheet breakup.
The effect of elasticity on sheet behaviour at the instant of breakup is analysed next. Figure 5 shows the effect of El on final sheet profile for different gas to liquid velocity ratios keeping 
$\lambda = 0$, We = 300 and ρ = 0.0014. In the case of low velocity ratios (U = 2 and U = 2.15), 
${\hat{\eta }_2}$ is negative at El = 0.1 (see figure 4). Hence, sheet breakup for El = 0.1 occurs at 1/4 and 3/4 of the fundamental wave, as shown in figures 5(a) and 5(b). However, as El increases (El = 10 and 100), the breakup point tends to shift towards half and full wavelength, owing to the corresponding positive value of 
${\hat{\eta }_2}$. It has already been observed that the effect of second-order amplitude is less for low velocity ratios. As a result, the sheet experiences sufficient deformation and high wave amplitude at the time of breakup for U = 2 and 2.15. In the case of high velocity ratios (U = 2.5 and 2.75), 
${\hat{\eta }_2}$ is positive for the entire range of El (see figure 4), resulting in breakup near half and full wavelength. Moreover, wave amplitude is less due to the stronger effect of nonlinearity. The final sheet appears to be more distorted and takes a saw tooth shape.
Final sheet profile for different El for (a) U = 2, (b) U = 2.15, (c) U = 2.50 and (d) U = 2.75 keeping 
$\lambda = 0$, We = 300 and ρ = 0.0014.
Subsequent discussion examines the effect of time constant ratio 
$\lambda = {\lambda _2}/{\lambda _1}$. It represents the effect of deformation retardation time for constant value of El. Figures 6(a) and 6(b) show the variation of maximum growth rate 
$({\beta _{max}})$ and dominant wavenumber with 
$\lambda $, respectively, for different gas to liquid velocity ratios keeping El = 10, We = 300 and ρ = 0.0014. An increase in 
$\lambda $ increases the deformation retardation time and causes a reduction in first-order effective Reynolds number 
$(R{e_1})$. A smaller effective Reynolds number leads to dampening of linear growth rate, and hence 
${\beta _{max}}$ reduces with increasing 
$\lambda $. However, the effect of 
$\lambda $ on 
${\beta _{max}}$ is very weak at U = 2 and 2.15 and increases slightly with an increase in U (U = 2.50 and 2.75). Figure 6(b) shows that the dominant wavenumber almost remains invariant to change in 
$\lambda $ at U = 2 and 2.15. However, at high U (U = 2.50 and 2.75), an increase in 
$\lambda $ reduces the dominant wavenumber, suggesting formation of longer waves. The overall effect of time constant ratio on linear growth rate and dominant wavenumber appears to be opposite to that of elasticity number.
Variation of (a) maximum growth rate and (b) dominant wavenumber with 
$\lambda $ at different U keeping El = 10, We = 300 and ρ = 0.0014.
Figure 7 shows the effect of time constant ratio 
$(\lambda )$ on breakup time at different gas to liquid velocity ratios (U) keeping El = 10, We = 300 and ρ = 0.0014. It is observed that the influence of 
$\lambda $ on breakup time changes significantly with change in velocity ratio. Breakup time decreases monotonically with an increase in 
$\lambda $ at U = 2, whereas it shows a monotonic increase with an increase in 
$\lambda $ at U = 2.50 and 2.75. This signifies a destabilizing and stabilizing effect of 
$\lambda $ at low (U = 2) and relatively high (U = 2.50 and 2.75) velocity ratio, respectively. However, at U = 2.15, 
$\lambda $ exhibits both stabilizing and destabilizing effects depending upon the range of 
$\lambda $. A critical value of 
$\lambda $ 
$(\lambda = 0.4)$ is identified, below and above which 
$\lambda $ increases and decreases the breakup time, respectively. Thus, an increase in time constant ratio results in two regimes at U = 2.15: the first regime indicates a stabilizing effect of 
$\lambda $ for 0.1 < λ < 0.40 and the second regime indicates a destabilizing effect of 
$\lambda $ for 0.40 < λ < 1. Similar to elasticity number, the dual behaviour of time constant ratio can be explained by studying its effect on second-order amplitude 
$({\hat{\eta }_2})$.
Variation of breakup time with 
$\lambda $ at different U keeping El = 10, We = 300 and ρ = 0.0014.
Variation of second-order amplitude 
$({\hat{\eta }_2})$ with 
$\lambda $ at different gas to liquid velocity ratios is presented in figure 8. A comparison between figures 6 and 8 shows that similar to elasticity, the effect of time constant ratio on 
${\hat{\eta }_2}$ is much greater than that on linear growth rate for all velocity ratios.
Variation of second-order amplitude with 
$\lambda $ at different U keeping El = 10, We = 300 and ρ = 0.0014.
At U = 2, 
${\hat{\eta }_2}$ is negative in the entire range of 
$\lambda $. An increase in 
$\lambda $ reduces 
${\hat{\eta }_2}$ but causes a weak enhancement of its absolute value 
$(|{\hat{\eta }_2}|)$. Hence, time constant ratio enhances nonlinearity and exhibits a destabilizing effect on the liquid sheet at U = 2. The effect of 
$\lambda $ on 
${\hat{\eta }_2}$ at higher velocity ratios U = 2.50 and 2.75 becomes opposite. It can be seen that 
${\hat{\eta }_2}$ is positive for the entire range of 
$\lambda $. An increase in 
$\lambda $ decreases 
${\hat{\eta }_2}$ and hence dampens nonlinearity. Thus, breakup time increases with time constant ratio for U = 2.50 and 2.75. At U = 2.15, 
${\hat{\eta }_2}$ is initially positive. An increase in 
$\lambda $ decreases 
${\hat{\eta }_2}$ until it becomes negative at 
${\lambda } = 0.4$. Hence, in this regime, 
$\lambda$ suppresses nonlinearity and makes the sheet more stable, resulting in longer breakup time. However, beyond 
${\lambda} = 0.4$, 
${\hat{\eta }_2}$ is negative and an increase in 
$\lambda $ increases 
$|{\hat{\eta }_2}|$. This leads to a second regime, where 
$\lambda $ enhances nonlinearity, and thereby produces faster breakup. This provides an explanation for the dual effect of 
$\lambda $ at U = 2.15.
Figure 9 shows the effect of 
$\lambda $ on final sheet profile at different gas to liquid velocity ratios keeping El = 10, We = 300 and ρ = 0.0014. At U = 2, 
${\hat{\eta }_2}$ is negative in the entire range of 
$\lambda $. Hence, breakup always takes place at 1/4 and 3/4 of the fundamental wave. Also, the effect of 
$\lambda $ on sheet profile is less and the sheet remains sinuous in shape for all values of 
$\lambda $. At U = 2.15, disturbance amplitude initially increases as 
${\lambda }$ increases from 0.1 to 0.4 but reduces with a further increase in 
$\lambda $ from 0.4 to 0.9.
Final sheet profile for 
${\lambda}$ for (a) U = 2, (b) U = 2.15, (c) U = 2.50 and (d) U = 2.75 keeping El = 10, We = 300 and ρ = 0.0014.
The reason can be attributed to initial increase and subsequent decrease in nonlinearity between 0.1 < λ < 0.4 and 0.4 < λ < 0.9, respectively, as observed earlier. Also, the breakup point is near half and full wavelength for λ = 0.1 due to corresponding positive values of 
${\eta _2}$. But at 
${\lambda } = 0.4$ and 0.9, the breakup point shifts to 1/4 and 3/4 of the fundamental wave, due to corresponding negative values of 
${\hat{\eta }_2}$. For U = 2.50 and 2.75, 
${\hat{\eta }_2}$ is positive for the entire range of 
${\lambda }$. Hence, sheet breakup always takes place at half and full wavelength. Also, the distortions appear to reduce and disturbance amplitude to increase with an increase in 
$\lambda $.
Figure 10 shows the effect of El on droplet size distribution at different velocity ratios (U) keeping 
$\lambda = 0$, We = 300 and ρ = 0.0014. It can be observed that an increase in El shifts the droplet distribution curve towards finer droplets. This effect is more prominent for 0.1 < El < 10 at higher velocity ratios U = 2.50, 2.75. The reason can be attributed to an increase in effective Reynolds number with an increase in El.
Effect of El on droplet size distribution for (a) U = 2, (b) U = 2.15, (c) U = 2.50 and (d) U = 2.75 keeping 
${\lambda } = 0$, We = 300 and ρ = 0.0014.
At a constant Weber number, an increase in Reynolds number implies a reduction in fluid viscosity (Yang et al. Reference Yang, Wang, Fu, Du and Tong2013). In most practical conditions, a decrease in fluid viscosity promotes faster deformation of ligaments between two breakup points into droplets, so that atomization takes place well upstream where the velocity of liquid is comparatively high. This leads to formation of finer droplets.
However, experimental studies of dilute polymeric solutions reveal that an increase in viscoelasticity inhibits fragmentation and results in larger mean droplet diameter (Thompson & Rothstein Reference Thompson and Rothstein2007). In this regard, it must be emphasized that the actual dynamics of viscoelastic fluid near the pinch-off region is very different from that of Newtonian fluid. Theoretical investigation of the process of fragmentation reported by Villermaux et al. (Reference Villermaux, Marmottant and Duplat2004) indicated that the atomization of Newtonian spray is controlled by the size and geometry of initial ligaments. The study provides a detailed analysis as to why gamma distributions are better in fitting drop size distribution data, as compared with log-normal or Poisson distribution. This approach was later followed by Kooij et al. (Reference Kooij, Sijs, Denn, Villermaux and Bonn2018) for fitting drop size distribution data for a Newtonian fluid, under the influence of different nozzle operating conditions and liquid properties. A good agreement between predicted and experimentally measured distribution was obtained, confirming that final droplet size distribution is governed by ligament sizes and distribution. Keshavarz et al. (Reference Keshavarz, Houze, Moore, Koerner and McKinley2016) performed an experimental investigation of ligament-driven fragmentation dynamics of viscoelastic fluid. The study revealed that, as opposed to the case of Newtonian fluid where rapid droplet pinch-off happens due to capillary thinning, atomization of viscoelastic fluids is governed by markedly different series of events. Through a montage of images, it was observed that, as the flow squeezed near the neck of ligaments, there was a large variation in axial diameter of the localized blobs, resulting in very corrugated ligaments. The resulting ligament profile gradually deviated from purely cylindrical shape to wavelike modulation of arbitrary shape. The final distribution consisted of a greater number of large and small droplets ultimately leading a larger average droplet size and broader distribution range. The study identified enhanced extensional viscosity of viscoelastic fluid as the primary reason for elongated ligaments near the neck region. Several other experimental studies have also shown that for viscoelastic ligaments, there is the occurrence of high strain rate in the ligament neck which resists pinch-off and leads to formation of thin cylindrical threads of initially uniform and later non-uniform thickness, known as beads-on-string structure (Mun, Byars & Boger Reference Mun, Byars and Boger1998; Christanti & Walker Reference Christanti and Walker2001; Cooper et al. Reference Cooper, Fagan, Tirtaatmadja, Lester and Boger2002; Wagner et al. Reference Wagner, Amarouchene, Bonn and Eggers2005). At later stages of thinning of the ligament, when the threads come to their full extension and are of the order of several micrometres in radius, tiny beads may begin to appear resulting in a structure known as blistering (Sattler, Wagner & Eggers Reference Sattler, Wagner and Eggers2008). Experimental investigation of pearling instability in viscoelastic thread by Deblais, Velikov & Bonn (Reference Deblais, Velikov and Bonn2018) showed that the bead-on-string structure is controlled by an interplay between capillary and elastic forces, whereas the blistering instability is due to a dynamical phase separation that takes place in the elongational flow. However, neither of these phenomena are accessible in the context of the present study, which captures the sheet dynamics from the moment waves start growing on the surface of the interfaces until the distance between the two interfaces becomes a minimum. It considers finite time singularity at breakup, like that of Newtonian fluid, where rapid thinning of the ligament neck takes place due to formation of capillary waves, which ultimately leads to droplet pinch-off. Tracking the geometry of ligaments after inception of pinch-off points is a visualization challenge and is beyond the scope of this analytical study. However, the finding of this study is important from the standpoint of understanding the effect of viscoelastic properties on the growth of instabilities as well as probable droplet size distribution based on initial ligament area.
Figure 11 shows the effect of 
$\lambda $ on droplet size distribution at different velocity ratios (U) keeping El = 10, We = 300 and ρ = 0.0014. An increase in time constant ratio reduces the effective Reynolds number. As discussed earlier, a decrease in Reynolds number delays breakup and produces larger droplets. Hence, the size distribution curve tends to shift towards larger diameter with an increase in 
$\lambda $. However, similar to El, the effect of 
$\lambda $ on droplet size is very minor at low velocity ratios (U = 2, 2.15) and increases with increasing U (U = 2.50 and 2.75).
Effect of 
${\lambda }$ on droplet size distribution for (a) U = 2, (b) U = 2.15, (c) U = 2.50 and (d) U = 2.75 keeping El = 10, We = 300 and ρ = 0.0014.
6. Conclusion
The study involves a weakly nonlinear temporal analysis of a planar viscoelastic liquid sheet sandwiched between two inviscid gas streams of equal velocities. The rheological model of the liquid sheet is considered as corotational Jeffrey's model. The sheet is subjected to initial sinuous mode of disturbance. Perturbation analysis is employed to obtain the second-order interface displacement. Effects of viscoelastic properties such as elasticity number (El) and time constant ratio (λ) are discussed for different gas to liquid velocity ratios. Linear analysis predicts that elasticity number and time constant ratio monotonically enhance and dampen the linear growth rate, respectively, at all velocity ratios. However, the effect of these properties on sheet breakup is non-monotonic and depends upon the range of the viscoelastic properties as well as the surrounding gas velocities. At relatively low velocity ratios, elasticity produces a stabilizing effect for 0.1 < El < 4 at U = 2 and 0.1 < El < 1 at U = 2.15. A further increase in elasticity beyond this initial range of El produces a destabilizing effect. However, at U = 2.50 and 2.75, elasticity only produces a destabilizing effect on the liquid sheet. Time constant ratio exhibits a destabilizing effect at U = 2, whereas it stabilizes the sheet at U = 2.50 and 2.75. At intermediate gas liquid velocity ratio U = 2.15, two regimes are identified representing stabilizing and destabilizing effects of time constant ratio for 
$\lambda < 0.40$ and 
$\lambda > 0.40$, respectively. The reason for such dual behaviour of elasticity and time constant ratio is obtained by studying their corresponding effect on second-order amplitude. Liquid sheet breakup under the influence of initial sinuous waves essentially takes place only when the first harmonic is added to the fundamental. The effect of first harmonic is represented using second-order amplitude 
${\hat{\eta }_2}$. It is observed that the effect of elasticity number and time constant ratio on 
${\hat{\eta }_2}$ is dominant over their corresponding effect on linear growth rate. A highly coupled nonlinear interaction of elasticity number with gas velocities or time constant ratio with gas velocity changes the property of 
${\hat{\eta }_2}$ and is responsible for their dual effect on sheet breakup. Practical applications involving non-Newtonian fluids usually employ an air assist atomizer which uses the kinetic energy of surrounding gases to disintegrate the liquid sheet. Hence, the findings of this work are significant as they throw light on the effect of viscoelastic properties on sheet behaviour in the presence of different gas velocities.
Size distribution of droplets generated after primary breakup is also investigated using an MEF-based non-empirical model, which employs a modified Newton–Raphson method with Taylor series expansion up to second order to eliminate the chances of divergence of solutions. It has been observed that an increase in elasticity number and time constant ratio produces finer and larger droplets, respectively. Their corresponding influence on effective Reynolds number has been identified as the primary reason for such an effect on droplet size distribution. However, the present analytical study does not address the effect of viscoelastic properties after the inception of pinch points. Moreover, breakup of interconnected ligaments into secondary droplets is also not addressed by the present MEF model. Hence, the present MEF model may not be suitable for studying the isolated effects of viscoelastic properties near the pinch region. A more complete modelling of secondary breakup in future or adoption of direct numerical simulations may provide better prediction of droplet size.
Declaration of interests
The authors report no conflict of interest.
Appendix A
Terms having 
$exp [2\textrm{i(}{k_{1x}} - {\omega _1}t\textrm{)}]$ as coefficient in (2.29) to (2.32) and (2.36) are given as 
${\hat{E}_{21u}} = {u_1}{u_{1,x}} + {v_1}{u_{1,y}}$; 
${\hat{E}_{21v}} = {u_1}{v_{1,x}} + {v_1}{v_{1,y}}$; 
${E_{21\eta }} = {u_1}{\eta _{j1,x}} - {\eta _{j1}}{v_{1,y}}$; 
${\hat{E}_{21\tau }} = {\eta _{j1}}(\partial /\partial y){\tau _{1yx}} + ({\tau _{yy1}} - {\tau _{xx1}}){\eta _{j1,x}}$; and 
${\hat{E}_{21d}} ={-} {\eta _{j1}}{P_{l1,y}} - [{\eta _{j1,x}}({\tau _{1,xy}} + {\tau _{1,yx}})] + {\eta _{j1}}(\partial /\partial y){\tau _{1,yy}}$
$- \rho {\eta _{j1}}({U_j}{\varphi _{gj1,yx}} + {\varphi _{gj1,yt}}) - {\textstyle{1 \over 2}}\rho (\varphi _{gj1,x}^2 + \varphi _{gj1,y}^2)$.
Corresponding expressions are obtained as
\begin{gather}{\hat{E}_{21u}} =
\frac{\textrm{i}}{{2{k_1}}}[2(A_1^2k_1^2 + B_1^2l_1^2) +
{A_1}{B_1}{({l_1} - {k_1})^2}\cosh ({l_1}y + {k_1}y)\nonumber\\ \hspace{-2.2pc} +\,
{A_1}{B_1}{({l_1} + {k_1})^2}\cosh ({l_1}y -
{k_1}y)],\end{gather}
\begin{gather}{\hat{E}_{21v}} = 0,\end{gather}
\begin{gather}{\hat{E}_{21\eta }} ={-} 2{\hat{\eta }_{j1}}[{A_1}{k_1}\sinh ({k_1}y) + {B_1}{l_1}\sinh ({l_1}y)],\end{gather}
\begin{gather}{\hat{E}_{21\tau }} = {\hat{\eta }_{j1}}\left[ {6i{A_1}k_1^2\sinh ({k_1}y) + \frac{{{l_1}}}{{{k_1}}}(l_1^2 + 5k_1^2)i{B_1}\sinh ({l_1}y)} \right],\end{gather}
\begin{gather}{\hat{E}_{21d}} = \left[ {\frac{{(4l_1^2 + 2k_1^2)}}{{Re}}{B_1}\cosh ({l_1}y) + \frac{{(l_1^2 + 5k_1^2)}}{{Re}}{A_1}\cosh ({k_1}y)} \right]{\hat{\eta }_{j1}}\nonumber\\ \hspace{-1.95pc}+\, \eta _{j1}^2\rho {({\omega _1} - U{k_1})^2}exp [{k_1} + {( - 1)^j}{k_1}y].\end{gather}The details of the particular solutions for liquid pressure and liquid velocity as mentioned in (3.16)–(3.18) are given below.
Term 
${\hat{S}_{21p}}$ is the particular solution of
\begin{align}\begin{aligned}
&\left( {\dfrac{{{\partial^2}}}{{\partial
{y^2}}}{{\hat{P}}_{l21}} - 4k_1^2{{\hat{P}}_{l21}}} \right)
={-} 2\textrm{i}{k_1}{{\hat{E}}_{21u}} +
\dfrac{1}{{Re}}\dfrac{{{\lambda _2} - {\lambda _1}}}{{[1 +
{\lambda _1}(\textrm{i}{k_1} - \textrm{i}{\omega _1})][1 +
2{\lambda _1}(\textrm{i}{k_1} - \textrm{i}{\omega _1})]}}\\
&\quad \left[
{\dfrac{{{\partial^2}{{\hat{D}}_{21c,yy}}}}{{\partial
{y^2}}} + 4\textrm{i}{k_1}\dfrac{{\partial
{{\hat{D}}_{21c,yx}}}}{{\partial y}} -
4k_1^2{{\hat{D}}_{21c,xx}}} \right].
\end{aligned}\end{align}
The expression for 
${\hat{S}_{21p}}$is obtained as
\begin{align}
{{\hat{S}}_{21p}} & ={-} \dfrac{1}{{2k_1^2}}(A_1^2{k_1} +
B_1^2{l_1}) + {A_1}{B_1}\left[ \begin{array}{l}
\dfrac{{{l_1} + {k_1}}}{{({l_1} - 3{k_1})}}\cosh ({l_1}y -
{k_1}y)\nonumber\\ \quad + \dfrac{{{l_1} - {k_1}}}{{({l_1} +
3{k_1})}}\cosh ({l_1}y + {k_1}y) \end{array} \right]\nonumber\\ &
\quad + \frac{1}{{Re}}\frac{{({\lambda _2} - {\lambda
_1})}}{{[1 + {\lambda _1}(\textrm{i}{k_1} -
\textrm{i}{\omega _1})][1 + 2{\lambda _1}(\textrm{i}{k_1} -
\textrm{i}{\omega _1})]}}\nonumber\\ & \quad \left[ \begin{array}{l}
- 2l_1^2B_1^2 - 2k_1^2A_1^2 - \dfrac{1}{2}k_1^2B_1^2 -
\dfrac{1}{{2k_1^2}}l_1^4B_1^2 - \dfrac{{{{(l_1^2 -
k_1^2)}^2}}}{{2k_1^2}}B_1^2\cosh (2{l_1}y)\nonumber\\ \quad -
\dfrac{{{{({l_1} - {k_1})}^2}(2{l_1} + 4{k_1})}}{{({l_1} +
3{k_1})}}{A_1}{B_1}\cosh ({l_1}y - {k_1}y)\\ \quad -
\dfrac{{{{({l_1} + {k_1})}^2}(2{l_1} - 4{k_1})}}{{({l_1} +
3{k_1})}}{A_1}{B_1}\cosh ({l_1}y - {k_1}y) \end{array}
\right]. \end{align}
Terms 
${\hat{S}_{21u}}$ and 
${\hat{S}_{21v}}$ are the particular solutions of
\begin{gather}\begin{array}{ccccc} & {{\hat{u}}_{21,yy}} - l_{21}^2{{\hat{u}}_{21}} = R{e_{21}}[[2\textrm{i}{k_1}{c_1}exp (2{k_1}y) - 2{k_1}{c_2}exp ( - 2{k_1}y)] + {{\hat{S}}_{21p}}]\\ & \quad - \dfrac{{R{e_{21}}}}{{Re}}\dfrac{{{\lambda _2} - {\lambda _1}}}{{[1 + {\lambda _1}(\textrm{i}{k_1} - \textrm{i}{\omega _1})][1 + 2{\lambda _1}(\textrm{i}{k_1} - \textrm{i}{\omega _1})]}}\left( {2i{k_1}{{\hat{D}}_{21c,xx}} + \dfrac{\partial }{{\partial y}}{{\hat{D}}_{21c,yy}}} \right), \end{array}\end{gather}
\begin{gather}\begin{array}{ccccc} & {{\hat{v}}_{21,yy}} - l_{21}^2{{\hat{v}}_{21}} = R{e_{21}}[[2{k_1}{c_1}exp (2{k_1}y) - 2{k_1}{c_2}exp ( - 2{k_1}y)] + {{\hat{S}}_{21p,y}}]\\ & \quad - \dfrac{{R{e_{21}}}}{{Re}}\dfrac{{{\lambda _2} - {\lambda _1}}}{{[1 + {\lambda _1}(\textrm{i}{k_1} - \textrm{i}{\omega _1})][1 + 2{\lambda _1}(\textrm{i}{k_1} - \textrm{i}{\omega _1})]}}\left( {2\textrm{i}{k_1}{{\hat{D}}_{21c,xy}} + \dfrac{\partial }{{\partial y}}{{\hat{D}}_{21c,yy}}} \right). \end{array}\end{gather}
Corresponding expressions for 
${\hat{S}_{21u}}$ and 
${\hat{S}_{21v}}$are given as
\begin{gather}{\hat{S}_{21u}} = \frac{\textrm{i}}{{2{k_1}}}{S_{21v,y}}, \end{gather}
\begin{gather}
{{\hat{S}}_{21v}} = R{e_{21}}{A_1}{B_1}(l_1^2 -
k_1^2)\left\{ \begin{array}{l} \dfrac{{\sinh ({l_1}y +
{k_1}y)}}{{({l_1} + 3{k_1})[{{({l_1} + {k_1})}^2} -
l_{21}^2]}}\\ + \dfrac{{\sinh ({l_1}y -
{k_1}y)}}{{({l_1} - 3{k_1})[{{({l_1} + {k_1})}^2} -
l_{21}^2]}} \end{array} \right\}\nonumber \\ \hspace{0.4pc} +
\frac{{{\lambda _2} - {\lambda _1}}}{{[1 + {\lambda
_1}(\textrm{i}{k_1} - \textrm{i}{\omega _1})][1 + 2{\lambda
_2}(\textrm{i}{k_1} - \textrm{i}{\omega
_1})]}}{A_1}{B_1}\nonumber\\ \hspace{4.7pc}\times\left[ \begin{array}{l} - \dfrac{{(l_1^2 -
k_1^2)}}{{({l_1} + 3{k_1})\{ {{({l_1} + {k_1})}^2} -
l_{21}^2\} }}\sinh ({l_1}y + {k_1}y)\\ - \dfrac{{(l_1^2 -
k_1^2)}}{{({l_1} - 3{k_1})\{ {{({l_1} - {k_1})}^2} -
l_{21}^2\} }}\sinh ({l_1}y - {k_1}y) \end{array} \right].
\end{gather}
