Landweber iterative regularization method for reconstructing the unknown source of the modified Helmholtz equation

1.   Introduction

When sufficiently strong uniform shears exist between a fluid interface, the interface becomes unstable, and it wraps up into vortices (Figure 1). This phenomenon, known as the Kelvin-Helmholtz instability (KHI), appears in various areas such as geophysical fluid dynamics ^[1,2,3], internal gravity waves ^[4], and plasma physics ^[5]. In these systems, multi-interface flows appear mainly due to density stratification, and the interaction between the interfaces occurs. When the density stratification is absent or very small, the nonlinear interaction between the interfaces has been investigated numerically, giving the result that sufficiently close interfaces roll up simultaneously under strong shear flows ^[6,7,8]. This simultaneous roll-up is also found in density stratified flows when a uniform shear exists in one of the fluid regions ^[9,10]. These results suggest that interfaces in multi-layer flow with strong shears deform and roll up simultaneously.

Moore ^[11] considered the following undisturbed two-dimensional flow $(u(y), 0)$ in fixed rectangular axes $OX$, $OY$ (refer to Figure 2):

where two vortex sheets are located at $y = 0$ and $y = -d$ ($d > 0$), and $U_2 \approx -U_0$ is assumed (In his linear analysis, $|U_2|$ is assumed to be slightly larger than $|U_0|$). When $U_2 = -\theta U_0$ ($\theta \in {\mathbb R}$), putting $\epsilon = {{{\rm{e}}}}^{-kd} \ll 1$, $k$ the wavenumber, he obtained the result that the upper vortex sheet $y = {{{\rm{e}}}}^{ikx} \eta(t)$ evolvs as

under the conditions $\eta(0) = \eta_0$ and ${\dot \eta}(0) = 0$, where the constants $\alpha = O(1)$ and $\beta = O(\epsilon^2)$, and the lower sheet is assumed to be $y = -d + {{{\rm{e}}}}^{ikx} \mu(t)$, $\mu(0) = {\dot \mu}(0) = 0$. The linear growth rate $\sigma_1 \in {\mathbb C}$ and $\sigma_2 \in {\mathbb C}$ are given by $\sigma_1 = 1/2U_0 k(1 - i) + O(\epsilon^2)$ and $\sigma_2 = 1/2U_0 \theta k(1 + i) + O(\epsilon^2)$, respectively. Neglecting $O(\epsilon^2)$, $\sigma_1$ corresponds to the complex growth rate on a single vortex sheet with the same strength as the upper sheet, and $\sigma_2$ has the same significance for the lower sheet. The second term in (1.2) divided by the first is $O(R(t))$, where

When $R < 0$, i.e., when $t < t_c$, where

the evolution of the upper vortex sheet is determined only by $\sigma_1$ described by $U_0$, and it is not affected by $U_2$. This indicates that the upper vortex sheet evolves just as if the lower sheet were absent. The value $t_c$ does not depend on the wavenumber $k$. As we see from this estimate, $t_c \to \infty$ as $\theta \to 1$, i.e., as $U_2 \to -U_0$. When $t > t_c$, the linear theory fails, and the motion of the two vortex sheets is unable to calculate analytically. This analysis suggests that the two vortex sheets may not roll up simultaneously in the nonlinear region when the condition $U_2 \approx -U_0$ is satisfied. The long-time behavior of the vortex sheets with the above conditions is unknown. The purpose of this study is to examine the nonlinear motion of the two vortex sheets beyond Moore's linear analysis, taking the lower shear $U_2$ as $U_2 = -U_0$. Density stratification is not included in Moore's analysis. We also investigate the effect of density stratification on the motion of the two vortex sheets.

The long-time behavior of a vortex sheet has been investigated by various mathematical models and numerical methods such as the boundary integral method ^[12,13], the vortex (blob) method ^{[14,15,16,17,18,19,20]}, and the contour dynamics ^[21,22,23], in which the boundary integral method and the vortex method adopt a singular integral equation called the Birkhoff-Rott equation ^[24,25,26] to calculate the interfacial velocity. The vortex method is a regularization of the boundary integral method ^[27]. The boundary integral method without regularization is known to provide a weak solution to the two-dimensional Euler equation as long as the initial condition is sufficiently smooth. Then the obtained numerical solution uniformly converges to the Birkhoff-Rott equation ^[15,27,28] with spectral accuracy if we adopt an appropriate spatial integration method such as the alternate point quadrature method ^[29]. Mathematically, this method provides the most accurate numerical solution that approximates the exact solution. However, this high-accuracy calculation breaks down with the appearance of a curvature singularity (curvature divergence) on the interface ^[30] before the roll-up of the vortex sheet appears ^[13]. To avoid the occurrence of the curvature singularity and calculate the vortex sheet motion for a long time, Krasny ^[17] introduced a regularized parameter $\delta$ (often called Krasny's $\delta$) and succeeded in calculating the detailed structure of the roll-up of a vortex sheet. Extending this vortex method to the sheet motion with density stratification, Matsuoka et al. ^{[31,32,33,34,35,36,37]} succeeded to capture the complicated interfacial dynamics in the Rayleigh-Taylor instability (RTI) and the Richtmyer-Meshkov instability (RMI) ^[38,39]. This method is also applicable to multi-component density stratified systems with multiple interfaces ^[10,40], in which the merger of two vortex sheets are observed. The merging of vortex sheets is found when the compressibility is weak and the dissipation is finite. Tsoutsanis et al. capture this merging phenomenon in the double vortex pairing process using the compressible and incompressible structured-grid method ^[9]. In the boundary integral method or the vortex method, the same evolution equation (Bernoulli equation) as the analytical calculation is adopted; therefore, we can use the result of linear analysis straightforwardly as the initial condition of numerical computations. Moreover, these methods enable us to calculate various physical and mathematical quantities such as circulations, velocity fields, curvatures, and so on. In the current study, we adopt the vortex method to calculate the nonlinear motion of the interfaces.

As described above, linear analysis ^[10,40] is required to determine the initial conditions for numerical calculations. Generally, the algebra of linear ^{[41,42,43,44]} or weakly nonlinear analysis ^[45,46] for multi-component fluid systems is lengthy and analytical calculations are daunting tasks. Following the references ^[10,47], we calculate the linear solution to satisfy Moore's analysis described above by Newton's method. Using the obtained initial conditions, we perform the numerical computations by the vortex method and investigate the nonlinear development of two vortex sheets with uniform shears in opposite direction. Unlike multi-interface flows with uniform shear ^[10], the simultaneous roll-up (or merging) of two sheets is not found when $U_2 = -U_0$ in (1.1), and the initial amplitude and initial velocity of the lower vortex sheet are approximately zeros. For this case, the lower vortex sheet does not roll up, and the growth of the sheet strength, which corresponds to the magnitude of the shear velocity, is also suppressed for quite a long time. This suggests that the upper vortex sheet controls the motion of the lower sheet, as pointed out by Moore ^[11]. We also show that the suppression effect is more remarkable when a stable density stratification that the upper fluid is the lightest and the lower fluid is the heaviest exists in the system.

This paper is organized as follows. In Section 2, we perform the linear analysis and derive the dispersion relations to adopt as the initial conditions in numerical calculations. In Section 3, we overview the mathematical model and numerical methods to describe the nonlinear dynamics of two vortex sheets. In Section 4, we present some numerical results for verifying Moore's estimate in the nonlinear region. Section 5 is devoted to conclusion and discussions.

2.   Linear stability analysis

In this section, we perform a linear analysis to find the initial conditions satisfying Moore's condition. In the current study, we define Moore's condition as

under the condition that $|a_1| \gg |a_2|$, where $a_1$ and $a_2$ are the linear amplitudes of the upper and lower interfaces, respectively [refre to (2.5)]. The obtained results are used for numerical computations in Section 4. We first derive the analytical dispersion relation from the governing equations, then calculate its solution by use of the Newton's method. More detailed analyses are provided in the reference ^[10]. We consider a two-dimensional inviscid flow induced by two vortex sheets $I_1$ and $I_2$ such that a density stratification exists between them (Figure 2), where the undisturbed levels of $I_1$ and $I_2$ are set to $y = 0$ and $y = -d$ ($d > 0$), respectively. The fluid motion is assumed to be irrotational except on the interfaces $I_1$ and $I_2$. Then the velocity potential $\phi_i$ in the region $i$ ($i = 0, 1, 2$), which is related to the fluid velocity ${{\boldsymbol{ u }}}_i$ as ${{\boldsymbol{ u }}}_i = \nabla \phi_i$, satisfies the Laplace equation

and the Bernoulli equation:

where $\rho_i$ is the density of fluid $i$, $p_i$ is the pressure, $g$ is the gravitational acceleration, and $C_i(t)$ is the Bernoulli constant. Here, we assume that the Bernoulli constant does not depend on time, i.e., $C_i(t) = C_i$. By imposing the pressure continuous condition across the interfaces; $p_0 = p_1$ at $I_1$ ($y = 0$) and $p_1 = p_2$ at $I_2$ ($y = -d$), and taking into the hydrostatic pressure condition $p_i = C_i - \rho_i g y$ ($i = 0, 1, 2$, $y = 0$ at $I_1$ and $y = -d$ at $I_2$), we obtain

Then the Bernoulli equation yields

Suppose that the vortex sheets $I_1$ and $I_2$ are evaluated as the deviations $y = \eta_1(x, t)$ and $y = \eta_2(x, t)$ from the undisturbed levels $y = 0$ and $y = -d$ (Figure 2), respectively. The kinematic boundary conditions at the vortex sheets are given by

Now we assume the undisturbed velocity profile (1.1). Suppose that the upper ($I_1$) and the lower ($I_2$) vortex sheets are described by

where $k$ is the wavenumber, $\omega$ is the linear frequency, and $\Re$ denotes the real part. Then the linearized kinematic boundary condition (2.4) gives the following solution for the velocity potential $\phi_i$ ($i = 0, 1, 2$)

in which the coefficients $B_0$, $B_{11}$, $B_{12}$ and $B_2$ are given as

Linearizing the Bernoulli equation (2.3) at the undisturbed interfaces $y = 0$ and $y = -d$, and using the result of (2.7), we obtain the following dispersion relation ^[10,47]:

where

and

When the amplitude $a_1$ in (2.5) is given, the amplitude $a_2$ for $I_2$ is determined by $a_1$, $D_1$ and $\Lambda$ as

We can also describe $a_1$ using $a_2$ as

In the current study, we adopt the relation (2.9) ($a_1$ is assumed to be given).

The linear frequency $\omega$ is determined by solving the dispersion relation $D(\omega, k) = 0$ in (2.8). This calculation is also possible to perform analytically; however, the dispersion relation (2.8) is a fourth-order equation with respect to $\omega$, and the analytical representation is too complicated to understand the $k-\omega$ relations. Therefore, we adopt the Newton's method to obtain the solution in the current study. Here, we select the spacing $d$ as $d = \pi/2$. Throughout this paper, we select the initial amplitude $a_1 = 0.2$ and the gravitational acceleration $g = 1$.

We present the imaginary part of $\omega$ obtained from the dispersion relation (2.8) and the absolute value of the amplitude $a_2$ obtained from (2.9) in Figures 3-6, in which the four branches of the solution are colored by black ($\omega_1$), blue ($\omega_2$), green ($\omega_3$) and red ($\omega_4$). The fluid densities are set to the same value in Figures 3 and 4. The densities in Figure 5 are set to $\rho_0: \rho_1 : \rho_2 = 2/3 : 1: 3/2$ so that the upper fluid is the lightest ($\rho_0 < \rho_1 < \rho_2$), and the densities in Figure 6 are set to $\rho_0 : \rho_1 : \rho_2 = 3/2 : 1: 2/3$ so that the upper fluid is the heaviest ($\rho_0 > \rho_1 > \rho_2$). Here, we select the undisturbed flow in these figures as $U_2 = -U_0$ except Figure 4, in which the value of $U_2$ is given by $U_2 = -2 U_0$; i.e., Figure 4 does not satisfy Moore's condition. The value $U_2 = -U_0$ corresponds to the value that $t_c \to \infty$ in the estimate (1.3). The mode 1 (black) in Figure 5(b) becomes extremely large in the neighborhood of $k = 0.335$ due to $|\omega_r| \ll 1$ and $\omega_i = 0$ ($\omega = \omega_r + i \omega_i$) (refer to (2.9)), although it is finite.

In the current system, any of the four modes finally become unstable, but the mode having the largest positive $\omega_i$, the imaginary part of $\omega$, gives the fastest growing mode. In order to apply Moore's estimate described in the previous section, the condition $a_2 \approx 0$ is required as given by (2.1). To satisfy this condition, we adopt the mode that $\omega_i > 0$ and $|a_2| \ll 1$ at $k = 1$ (this mode is not necessarily the largest $\omega_i$, but it is required to be $\omega_i > 0$). We mention that the real part of $\omega$ is not important in investigating the nonlinear evolution of the vortex sheets. We adopt the values at $k = 1$ in the above four cases (Figures 3–6) as the initial conditions for numerical calculations in Section 4 because $a_2 \approx 0$ is realized at this wavenumber.

3.   Mathematical modeling for numerical computations

The detailed derivation of the mathematical model for numerical computations is presented in ^[40]. Here, we briefly review the numerical methods. We assume that the interfaces are $L$-periodic in the $x$ direction, where $L$ is the system size, which corresponds to the wavelength given by the wavenumber $k$ as $L = 2\pi/k$ in the current system. Suppose that the interfaces $I_i$ ($i = 1, 2$) are described by ${{\boldsymbol{ x }}} = {{\boldsymbol{ X }}}_i$, and we parameterize points on these interfaces as

using a Lagrangian parameter $e$ ($-L/2 \leq e \leq L/2$). Note that the two interfaces $I_1$ and $I_2$ are parameterized by the same Lagrangian parameter $e$. The (average) fluid velocity ${{\boldsymbol{ W }}}$ at an arbitrary point ${{\boldsymbol{ x }}} = (x, y)$ is described by the vortex induced velocities ${{\boldsymbol{ W }}}_1$ and ${{\boldsymbol{ W }}}_2$ due to the motion of interfaces $I_1$ and $I_2$ as

where ${{\boldsymbol{ W }}}_i = (W_{i, x}, W_{i, y})$ ($i = 1, 2$) is given by

in which

denotes the (true) vortex sheet strength of interface $I_i$ derived from the circulation $\Gamma_i \equiv \phi_i - \phi_{i-1}$, where $s_i$ is the arc length, and $\boldsymbol{t }_i$ is the unit tangent of the interface $I_i$, respectively. The subscript $e$ denotes the differentiation with respect to $e$ and $s_{i, e} = \sqrt{X_{i, e}^2 + Y_{i, e}^2}$. The integral (3.2) becomes a principal value integral when $(x, y) = (X_i, Y_i)$. Equations (3.1) and (3.2) correspond to the Birkhoff-Rott equation ^[24,25,26] for two interfaces, in which the velocity ${{\boldsymbol{ W }}}_i$ corresponds to the average velocity $({{\boldsymbol{ u }}}_i + {{\boldsymbol{ u }}}_{i-1})/2$ of the two fluid velocities ${{\boldsymbol{ u }}}_i$ and ${{\boldsymbol{ u }}}_{i-1}$.

The fluid velocities of the upper and the lower sides of the interfaces $I_1$ and $I_2$ are given by

where ${{\boldsymbol{ u }}}_{11}$ and ${{\boldsymbol{ u }}}_{12}$ are the velocities of the lower side of interface $I_1$ and the upper side of interface $I_2$ in the fluid region 1, respectively. There are two interfaces in the current system, and accordingly, two Lagrangian velocities exist associated with their motion. The normal component of the fluid velocity across each interface should always be continuous; however, there is an arbitrariness how to select the tangential velocity at the interfaces ^{[31,48,49,50]}. Introducing the Atwood number $A_i$ ($i = 1, 2$) describing the density ratio between fluids $i$ and $i-1$ ($i = 1, 2$) as

we define the interfacial velocity ${{\boldsymbol{ u }}}_i^+$ at each interface $I_i$ as

so that the interfacial velocities ${{\boldsymbol{ u }}}_1^+$ and ${{\boldsymbol{ u }}}_2^+$ become the weighted averages

Equating ${{\boldsymbol{ u }}}_i^+$ with the evolution of the interface $I_i$, we obtain the interfacial velocity for the Lagrangian motion as

where ${{{\rm{d}}}} /{{{\rm{d}}}} t$ denotes the Largange derivative moving with the velocity ${{\boldsymbol{ u }}}_i^+$. Rewriting the the Bernoulli equation (2.3) using the relation $\Gamma_i = \phi_i - \phi_{i-1}$ and differentiating the obtained equation with respect to $e$, we obtain the evolution equation for the sheet strength $\gamma_i$

where $T_i = \boldsymbol{t }_i \cdot {{\boldsymbol{ W }}}_i$. Equation (3.8) is the simultaneous Fredholm integral equations of the second kind. Solving (3.7) and (3.8) simultaneously taking (3.1) and (3.2) into account, we can determine the motion of two vortex sheets.

4.   Nonlinear evolution

In this section, we present the nonlinear evolution of two vortex sheets based on the linear analysis in Section 2 using the governing Eqs (3.7) and (3.8) provided in Section 3. As described in the last paragraph in Section 2, the condition $a_2 \approx 0$ is required to apply Moore's estimate (1.3). Since Moore's analysis does not depend on the wavenumber $k$ and that condition is realized in the neighborhood of $k = 1$ in our linear analysis in Figures 3–6, we adopt $k = 1$ in all quantities requited for numerical computations. Taking this into account, the initial conditions for numerical computations are selected as

where $\phi_i$ ($i = 0, 1, 2$) are given by (2.6) and (2.7). The initial conditions (4.1) are identical to those used in the reference ^[10].

In the current study, we adopt the vortex method ^{[14,15,16,17,18,19,20]} to investigate the nonlinear evolution of the vortex sheets. In numerical computations, we introduce a regularized parameter $\delta$ ^[17] in the the denominators of the singular integrals (3.2) so that the denominators do not tend to zero when ${{\boldsymbol{ x }}} = {{\boldsymbol{ X }}}_i$. We mention that this regularization is only needed for the calculation of ${{\boldsymbol{ W }}}_i$, and it is unnecessary for the calculation of ${{\boldsymbol{ W }}}_j$ when ${{\boldsymbol{ x }}} = {{\boldsymbol{ X }}}_i$, ($i \ne j$, $i, j = 1, 2$) in the velocity (3.1). The numerical results deviate from the analytical calculations performed in Section 2 when $\delta \ne 0$ ^{[10,13,32,33]}, and the deviation becomes more significant as the parameter $\delta$ becomes large. However, the curvature singularity ^[30] occurs in the limit of $\delta = 0$, and the computation breaks down before the roll-ups appear ^{[10,13,31,40]}, which does not describe the real experiments ^[51,52] or direct numerical simulations ^[35]. Krasny et al. ^[17,52] concluded that the values $\delta = 0.1 \sim 0.2$ are appropriate for the vortex method. Following these studies, we select the value of $\delta$ as $\delta = 0.1$ throughout this paper. For reference, the calculations with smaller values of $\delta$ are provided in A. In numerical computations, the interfacial coordinate $(X_i, Y_i)$ ($i = 1, 2$) is discretized using Fourier modes as ^[13,31,40]

Here, we adopt the alternate point quadrature method (alternative trapezoidal rule) ^{[13,29,31,40]} and the fourth-order Runge-Kutta scheme for the spatial integration and the temporal integration, respectively. For the spatial integration, we can also use the conventional trapezoidal rule for the current value of $\delta$. The accuracy of both methods is the same for finite $\delta$. When the discretization number $N$ is sufficiently large, the conventional trapezoidal rule becomes unstable as $\delta \to 0$. The simultaneous Fredholm equations of the second kind (3.8) are solved by iteration with tolerance level $10^{-12}$. For other techniques for numerical computations, refer to the references ^[10,40]. We select the number of grid points $N$ discretizing on the interfaces as $N = 1024$ and the time step $\triangle t$ as $\triangle t = 1.25 \times 10^{-3}$ for all calculations.

4.1.   Nonlinear evolution of vortex sheets without density stratification

In this subsection, we present the numerical results for $\rho_0 = \rho_1 = \rho_2$ ($A_1 = A_2 = 0$). Note that the gravitational term in (3.8) is negligible when $A_1 = A_2 = 0$. We consider the cases $U_2 = -U_0$ in Section 4.1.1 and $U_2 = -2 U_0$ in Section 4.1.2.

4.1.1.   Interfacial dynamics in uniform shears with the same magnitude

The initial conditions for numerical computations in this sub-subsection are provided in Figure 3, in which Moore's condition for the estimate (1.3) is almost satisfied. The linear frequency $\omega$ and the initial amplitude $a_2$ are given by

for $a_1 = 0.2$, $d = \pi/2$, and $k = 1$. These values are denoted by blue lines (mode 2) in Figure 3.

We present the interfacial structures with the colored scale of the vortex sheet strengths $\gamma_1$ and $\gamma_2$ in Figure 7, where the positive and negative strengths denote the counterclockwise and clockwise flow in the velocity field, respectively. As we see from the figure, the change of the lower interface $I_2$ is much smaller than that of $I_1$ in shape, and the roll-up is not observed even at the last computed time $t = 5$. Temporal evolution of the maximum values of $|\gamma_1|$ (blue line) and $|\gamma_2|$ (red line) are depicted in Figure 8. We see that the vortex sheet strength $|\gamma_2|$ almost does not grow when $t < 3$. The maximum sheet strength of $|\gamma_1|$ is realized in the neighborhood of vortex cores ^[53], the center of roll-ups (refer to Figure 7), around which a strong negative sheet strength (clockwise flow) is generated (refer to Figure 9). The velocity fields at $t = 0$ and $t = 5$ calculated using (3.2) are provided in Figure 9. The velocity field in the region between $I_1$ and $I_2$ at $t = 0$ is almost zero as found in the upper panel of Figure 9. A complicated velocity field is gradually induced in the flow region between the two vortex sheets due to the roll-up of the upper vortex sheet. As Moore predicted in the linear theory ^[11], the upper vortex sheet evolves as if the lower sheet were absent for quite a long time in Figure 7. This independent motion between the two vortex sheets is found for $d > \pi/2$ as well. As the distance $d$ becomes smaller, the interaction between the two vortex sheets becomes stronger, and the lower sheet also rolls up when $d \leq \pi/4$. The result for $d = \pi/4$ is presented in B.

4.1.2.   Interfacial dynamics in uniform shears with a different magnitude

The initial conditions for numerical computations in this sub-subsection are provided in Figure 4. All parameters are the same as those in Section 4.1.1 except the value of the shear, which corresponds to $U_2 = -2 U_0$. This value does not satisfy Moore's condition in the linear analysis. The linear frequency $\omega$ and the initial amplitude $a_2$ for numerical computations are given by

for $a_1 = 0.2$, $d = \pi/2$ and $k = 1$. These values are denoted by black lines (mode 1) in Figure 4. Although the value of $\omega$ in (4.4) is close to the one in (4.3), and $a_2 \approx 0$ is satisfied, the evolution of the two vortex sheets is different from that in Section 4.1.1. We see that below.

The interfacial structures with the colored scale of the vortex sheet strengths are provided in Figure 10. Unlike the evolution in Figure 7, the lower vortex sheet begins to roll up at around $t = 3$, which is also confirmed in the maximum sheet strength in Figure 11, where a sharp growth occurs in the curve of $|\gamma_2|$ (red line). The initial amplitude of $I_2$ is almost zero; however, quite a strong velocity field is induced initially in the region between $I_1$ and $I_2$ as found in the upper panel of Figure 12. This initial velocity fields cause the deformation of the lower vortex sheet, and finally, a vortex street ^[6,7,8,23] is formed as found in the panel at $t = 5$ in Figure 10. This is a typical multi-layer nonlinear KHI. A similar interfacial structure is also found in ^[10].

4.2.   Nonlinear evolution of vortex sheets with density stratification in uniform shear flows with the same magnitude

In this subsection, we present the numerical results for the motion of two vortex sheets with density stratification. All parameters except the Atwood numbers are the same in Sections 4.2.1 and 4.2.2. In Moore's linear analysis, density stratification effect is not considered. In the following sub-subsections, we investigate that effect and discuss the difference from the evolution of vortex sheets without density stratification.

4.2.1.   Interfacial dynamics when the upper fluid is the lightest

The dispersion relation in this sub-subsection is provided in Figure 5. The linear frequency $\omega$ and the initial amplitude $a_2$ for numerical computations are given by

for $a_1 = 0.2$, $d = \pi/2$ and $k = 1$. These values are denoted by black lines (mode 1) in Figure 5. Here, the density ratios are set to $\rho_0: \rho_1 : \rho_2 = 2/3 : 1: 3/2$ ($\rho_0 < \rho_1 < \rho_2$), which correspond to the Atwood numbers $A_1 = A_2 = -0.2$.

Figure 13 denotes the interfacial structures with the colored scale of the vortex sheet strengths. As we see from this figure, the lower vortex sheet does not evolve even at $t = 5$, despite the upper sheet largely deforms. The separation in motion between the two interfaces is more remarkable than that in the interfacial evolution for $A_1 = A_2 = 0$ found in Figure 7. We can also confirm that the lower interface does not evolve in Figure 14, in which the maximum value of $|\gamma_2|$ (red line) does not grow at all during the computation. The velocity fields at $t = 0$ and $t = 5$ are depicted in Figure 15. Not only the initial amplitude $a_2 \approx 0$ but also the initial velocity ${{{{\rm{d}}}} a_2}/{{{{\rm{d}}}} t} \approx 0$ holds in the current case, which almost completely satisfies the applicability conditions of Moore's estimate (1.3). Unlike Figure 9, we see that the velocity field is almost zero in the region between $I_1$ and $I_2$ even at the last computed stage $t = 5$. Although the density stratification effect is not included in Moore's analysis, his prediction that the upper vortex sheet evolves just as if the lower sheet were absent when $t < t_c$ is also true beyond the linear stage even in flows with finite Atwood numbers. Unlike the case in Section 4.1.1, when $A_1 = A_2 = -0.2$, the lower vortex sheet neither deforms nor rolls up, even for the narrower spacing $d \leq \pi/4$ (refer to B). On the other hand, the separation of motion between $I_1$ and $I_2$ as found in Figure 13 is also retained for $d > \pi/2$ as well as the case in Section 4.1.1.

4.2.2.   Interfacial dynamics when the upper fluid is the heaviest

The only one difference in Section 4.2.2 between Section 4.2.1 is in the density ratios. The dispersion relation in this sub-subsection is provided in Figure 6. The linear frequency $\omega$ and the initial amplitude $a_2$ for numerical computations are given by

for $a_1 = 0.2$, $d = \pi/2$ and $k = 1$. These values are denoted by red lines (mode 4) in Figure 6. The density ratios in the current computations are given by $\rho_0 : \rho_1 : \rho_2 = 3/2 : 1: 2/3$, which correspond to the Atwood numbers $A_1 = A_2 = 0.2$. Since $\rho_0 > \rho_1 > \rho_2$ (the upper fluid is the heaviest), KHI and RTI coexist in the current system.

The interfacial structures with the colored scale of the vortex sheet strengths is provided in Figure 16. Despite $a_2 \approx 0$, ${{{{\rm{d}}}} a_2}/{{{{\rm{d}}}} t} \approx 0$ and $U_2 = -U_0$, the lower vortex sheet grows and rolls up as found in this figure. This unstable motion is caused by RTI. As we see from Figure 17, the maximum sheet strength of $|\gamma_2|$ (red line) begins to grow at a quite early stage, and it increases up to the same value as the one of $|\gamma_1|$ (blue line) at the last computed stage $t = 5$. Figure 18 shows the velocity fields at $t = 0$ and $t = 5$. As with the case of Figure 15, the velocity field is not induced in the region between $I_1$ and $I_2$ at $t = 0$. However, a complicated velocity field is generated at $t = 5$ due to the density instability. This suggests that Moore's estimate (1.3) does not work in a system in which RTI (the density instability or the gravitational instability) occurs, even though the initial conditions $a_2 \approx 0$, ${{{{\rm{d}}}} a_2}/{{{{\rm{d}}}} t} \approx 0$ and $U_2 = -U_0$ are satisfied.

5.   Conclusions

We have investigated the nonlinear motion of two vortex sheets in a three-layer fluid with uniform shears in the opposite directions using the vortex method. When the uniform shears are equal in absolute value, and the initial amplitude and the initial velocity of the lower vortex sheet are approximately zero, the two vortex sheets evolve independently for quite a long time as predicted by Moore. Especially for the case that the fluid densities satisfying $\rho_0 < \rho_1 < \rho_2$ (the upper fluid is the lightest), the lower vortex sheet neither deformed nor rolled up while the upper sheet extensively deforms. This is a stabilization effect by density stratification. On the other hand, when the fluid densities satisfy the condition $\rho_0 > \rho_1 > \rho_2$ (the upper fluid is the heaviest), the lower vortex sheet deforms extensively at an early stage even though all conditions required for Moore's estimate are satisfied.

Generally, when the spacing $d$ is sufficiently large ($d \geq \lambda$, $\lambda$, the wavelength of the initial perturbation), the interaction between the two vortex sheets is weak, and each vortex sheet rolls up as if it is an independent sheet ^[6,23,40]. However, the results presented in Sections 4.1.1 and 4.2.1 are also true even for the spacing $d > \lambda$ ($\lambda = 2 \pi$ here). When the upper vortex sheet controls the lower sheet even if the two sheets are far apart from each other. In the current study, we calculate the case that the motion of the lower vortex sheet is suppressed. The reverse is also possible. When the initial amplitude of the upper vortex sheet $a_1$ is given by (2.10) for a given lower amplitude $a_2$, and the conditions $a_1 \approx 0$ and ${{{{\rm{d}}}} a_1}/{{{{\rm{d}}}} t} \approx 0$ are satisfied under the uniform shear $U_2 = -U_0$, the upper vortex sheet almost does not deform while the lower sheet rolls up. We mention that even though the deformation of the lower vortex sheet is extremely small, the current system is unstable, and the higher-order Fourier modes arise in the lower vortex sheet as with the upper sheet.

Acknowledgments

This work was supported by Grant-in-Aid for Scientific Research (C) (Grant No. 17K05371, Grant No. 18K03418 and Grant No. 21K03408) from the Japan Society for the Promotion of Science, the Osaka City University (OCU) Strategic Research Grant 2021 for top priority researches/basic researches/young researchers, the joint research project of ILE, Osaka University, Osaka City University Advanced Mathematical Institute (OCAMI) and the Research Institute for Mathematical Sciences, an International Joint Usage/Research Center located in Kyoto University.

Conflict of interest

The authors declare there is no conflicts of interest.

Appendix A.   Effect of the regularized parameter

In this appendix, we discuss the dependence of the regularized parameter on the interfacial structure. The validity of the vortex method and the choice of the regularized parameter $\delta$ has been discussed in various works ^{[17,31,35,52]}. The parameter $\delta$ provides a kind of numerical dissipation, and the system deviates from the perfect inviscid and incompressible flow when $\delta \ne 0$. We show the interfacial structures at $t = 3$ for $A = 0$ for various regularized parameters $\delta$ in Figure 19, where the left panels are calculated under the same conditions as those adopted in Sub-subsection 4.1.1, and the right panels are calculated under the same conditions as adopted in Sub-subsection 4.1.2; i.e., the left panels satisfy Moore's condition, while the right ones do not. As the regularized parameter $\delta$ becomes smaller, the roll-up of the interface $I_1$ becomes stronger for both conditions; however, the interface $I_2$ is hardly affected by the value of $\delta$ when Moore's condition is satisfied. This indicates that the vorticity concentration does not occur on the interface $I_2$ when Moore's condition is satisfied. Although there are some differences in the finer structure of roll-up, significant differences due to the decrease of $\delta$ are not found in asymptotic interfacial structures in Figure 19. The above tendency due to the decrease of $\delta$ is also true in finite Atwood number cases in Sub-subsections 4.2.1 and 4.2.2.

Appendix B.   Effects of density stratification and initial spacing

In this appendix, we discuss the effect of density stratification and initial spacing $d$. Figure 20 shows the interfacial structures with various Atwood numbers at $t = 5$, where $d = \pi/4$ and all figures satisfy Moore's condition. Here, the linear frequency $\omega$ and the initial amplitude $a_2$ are given by $\omega = 0.4129564 + 0.5239015 \, i$ and $a_2 = 0.0250247 + 0.0294614 \, i$ for $A_1 = A_2 = 0$ [Figure 20(a)], $\omega = 0.3152806 + 0.2517213 \, i$ and $a_2 = 0.0108192 - 0.0112804 \, i$ for $A_1 = A_2 = -0.2$ [Figure 20(b)], and $\omega = 0.6601108 + 0 \, i$ and $a_2 = 0.0125539 + 0 \, i$ for $A_1 = A_2 = -0.5$ ($\rho_0: \rho_1 : \rho_2 = 1/3 : 1: 3$) [Figure 20(c)]. As the initial spacing $d$ decreases, the lower interface without density stratification becomes unstable due to strong interaction between the two interfaces, and it also rolls up as found in Figure 20(a). However, even at this value of $d$, the shape of the lower sheet hardly changes, and the sheet strength also does not increase when density stratification exists as found in Figure 20(b), (c).