Effect of boundary roughness on nonlinear saturation of Rayleigh-Taylor instability in couple-stress fluid

1 Introduction

The growth rate of instability due to gravity associated with a dense layer overlying lighter layer, Rayleigh-Taylor instability (RTI) depends on constitutive law related to stress and strain rate. The upper fluid has a higher density compared to lower fluid, the interface is unstable to small perturbations whose amplification is well described by linear stability theory which depends on density ratio, gravity, coefficient of surface tension and viscosity (see Chandrasekhar [1]).

As the interface is destabilized by the gravity, configuration is unstable, but if there exist surface tension between two fluids, the configuration will have stabilizing effect. The abundant literature is available on linear RTI in composite two-fluid system is treated as a non-porous and fluid as Newtonian fluids. Stokes [2] has framed the theory of couple-stress for non-Newtonian fluid. The small amount of supplements in a lubricant can improve the bearing performance by increasing lubricant viscosity and producing an increase in load carrying capacity. This supplement in a lubricant decrease the coefficient of friction and also increases the temperature range in which bearing can operate. Rudraiah et al. [3] discussed the dispersion model of Gill and Sankarsubramanian [4], the effect of couple-stress is to facilitate the loss of hemoglobin, a disease known as hemolysis in blood. To comprehend the control of such situation, surface instability of RTI is studied. Sunil et al. [5] discussed the RTI of two superposed couple-stress fluids with unchanging densities in a horizontal magnetic field, they have carried out stability analysis for two highly viscous fluids of equal kinematic viscosity and equal couple-stress kinematic viscosity. Rudraiah et al. [6] analyzed the couple-stress fluid effects on the control of RTI at interface between dense fluids accelerated by lighter fluid and used the approximations to derive the growth rate of RTI. Another important application of RTI is flow past a porous layer is the design of effective porous wall insulation. Wetting of porous insulation has been noticed to persist despite many preventive devices discussed by Ognewicz and Tien [7]. This attributes to unavoidable formation of openings arising from faulty insulation and aging. The wetting phenomenon is generally believed to cause an important deterioration in the performance of wall insulation.

Although abundant literature is available on linear RTI, the work on nonlinear RTI is very scarce. The discovery of high-speed computing facility and nonlinear approach has changed significantly. If the perturbed quantities are insignificant in comparison with basic state then the product of perturbed quantities can be neglected which leads to linear stability theory. On the other hand if the perturbations are not small then their product cannot be ignored this leads to nonlinear theory. Let us give a brief review on the literature of nonlinear RTI. A mathematical model is used for both analytical and numerical approaches, of RTI [1]. McCrory et al. [8] discussed the simulation of RTI of abrasively accelerated thin-shell fusion targets and showed that nonlinear evolution exhibits spike amplitude saturation. A simple mathematical model is derived by Baker and Freeman [9] for nonlinear evolution of RTI. Babchin et al. [10] the nonlinear saturation of RTI in thin films and found that combined action of flow shear and surface tension is essence of saturation mechanism. A method of strained co-ordinates used for investigating nonlinear RTI problem given by Shivamoggi [11]. After incorporating the rectifications pointed out by Malik and Singh [12] and found a revised face for nonlinear cut-off wave number which separates the region of stability from that of instability. Mohamed and Shehawey [13] considered nonlinear electrohydrodynamic RTI in the absence of surface charges and a charge free surface separating two semi-infinite dielectric fluids inclined by a normal electric field subjected to nonlinear deformations. Allah and Yahia [14] analyzed the nonlinear RTI in the presence of magnetic field and also the heat and mass transfer using simplified formulation. Rudraiah et al. [15] discussed the various problems of nonlinear study of RTI in thin film past a porous layer.

Recently Rudraiah et al. [16] examined the non-linear study of electrohydrodynamic RTI in a composite fluid-porous layer. This problem has great impact of slip velocity at the interface between porous layer and thin film. Chavaraddi et al. [17] analyzed the hydrodynamic RTI in a composite fluid-porous layer. The main purpose of this article is to demonstrate such Rayleigh-Taylor instability in a fully nonlinear with effect of boundary roughness. We consider two-dimensional nonlinear saturation of RTI on two superposed couple-stress fluid and fluid saturated porous layer using slip condition proposed by Saffman [18] with the effect of boundary roughness.

The paper is organized as follows; the basic equations and corresponding boundary conditions are discussed in Section 2 using approximation proposed by Babchin et al. [10]. Section 3 is devoted to the study of nonlinear saturation of RTI in superposed composite couple-stress fluid and fluid saturated porous layer with the effect of boundary roughness. Finally results and discussions are drawn in Section 4.

2 Mathematical Formulation

Let us consider two superposed of fluids are compressible, immiscible and inviscid flowing down an inclined plane as shown in Figure 1. Cartesian Cartesian co-ordinate system aligned with the plane describes couple-stress fluid, x points in horizontal direction and y is vertical to it. The fluid properties of two layers are different; we use subscripts to distinguish them as ρ_f and ρ_p which denote densities in upper and lower fluids respectively.

Fig. 1

The physical configuration of the problem

The system involving of two semi-infinite inviscid, incompressible couple-stress fluid and fluid saturated porous media differentiated by an interface in the existence of surface roughness is considered. Since gravity is present; the effective acceleration is in the positive y-direction (upwards). Therefore, light fluid thrusts on heavy fluid (fluid saturated porous media). As far as the interface amid fluids remains uniform, that is completely horizontal and perpendicular to effective acceleration. The light fluid (couple stress fluid) has adequate pressure to hold heavy fluid against ceiling. However, small deviations are bound to occur at the interface. Whenever irregularities occur the portions of interface lie higher than average and therefore feel more pressure from light fluid is necessary to support water. Therefore, interface continues to increase in the spots. In the portions where the interface has dropped by small amounts below the average, more pressure is required for the support, and interface drops further. Even the smallest perturbation will cause the instability. If situation were reversed, heavy fluid layer was on the floor supporting the volume of light fluid then pressure gradients reverse. In the portions of interface have risen, the dense fluid will lower back to average, whereas the lowered portions will rise again. This qualitative analysis leads to conclusive statement that, if heavy fluid pushes the light fluid then the interface is stable and if light fluid pushes the heavy fluid then the interface is unstable and unstable case is RTI.

The upper fluid has density ρ_p, the lower one has density ρ_f and surface tension exist between two fluids is denoted by y.

To understand the physics of problem described here, it is simplified by following assumptions of Stokes and lubrication approximations (Babchin et al. [10]) as

1. The clear dense liquid is homogeneous and isotropic.

2. The thickness of film h is much lesser than thickness H of porous layer bounded above the film.

   i.e., h ≪ H

3. The Strouhal number S is assumed to be considerably small.

4. The surface elevation η is assumed to be small compared to film thickness h,

   i.e., η ≪ h

5. The fluid viscosity and thermal conductivity are assumed to be constants.

Based on the above assumptions and approximations, the basic equation leads to

where q⃗ = (u, v) is the fluid velocity, λ is the couple-stress parameter, p is the pressure, ∂p∂x is the pressure gradient, μ is the fluid viscosity and ρ is the fluid density of the fluid.

The relevant boundary conditions becomes

1. Roughness condition:

   −β1∂u∂y=uaty=0.(2.3)

2. Slip condition:

   ∂u∂y=−αk at y=h.Saffman [17](2.4)

3. Kinematic condition:

   v=∂η∂t+u∂η∂xaty=h.(2.5)

4. Dynamic condition:

   p=−∂η−γ∂2η∂x2aty=h.(2.6)

5. Couple-stress conditions:

   ∂2u∂y2=0aty=0andh.(2.7)

Let us non-dimensionalize the above equations by using

The above assumptions and approximations (i.e.Stokes and lubrication approximations), that the heavy fluid in porous layer is almost stationary because of creeping flow approximation and substituting Eq. (2.8) into Eqs. (2.1) – (2.2), we get (after neglecting the asterisks for simplicity)

where M0=λ/μh2 is the couple-stress parameter.

3 Dispersion Relation

To justify the use Stokes and lubrication approximations and derive the closed form of dispersion relation from appropriate governing equation (2.10) subjected to boundary and surface conditions.

1. Boundary roughness condition

   −β1∂u∂y=uaty=0.(3.1)

2. Saffmann slip condition

   ∂u∂y=−βuaty=1.(3.2)

   where β₁ is the roughness parameter, σp=hk the porous parameter and β = α_pσ_p.

3. Couple-stress conditions:

   ∂2u∂y2=0aty=0and1.(3.3)

4. Kinematic condition:

   v=∂η∂t+u∂η∂xaty=1.(3.4)

5. Dynamic condition:

   p=−η−1B∂2η∂x2aty=1.(3.5)

   where B = δh²/y is the Bond number and η = η(x, y, t) is the elevation of interface.

   The solution of Eq. (2.10) subject to the boundary conditions (3.1–3.3)

   u={C1cosh⁡(y/M0)+C2sinh⁡(y/M0)+12y2+M02+Ay+B}P(3.6)

   where P=∂p∂x,a1=1M0 sinh(1/M₀) + β cosh(1/M₀), a₂ = 1M0cosh⁡(1/M0)+βsinh⁡(1/M0),C1=−M02,C₂ = M02[cosh⁡(1/M0)−1]sinh⁡(1/M0)

   A={−122+β(1+2M02)+M02a1+M02[1−cosh⁡(1/M0)]sinh⁡(1/M0)a2−β1M0[1−cosh⁡(1/M0)]sinh⁡(1/M0)−M02β}/1+β(1−β1)B=M02−β1M0cosh⁡(1/M0)−1sinh⁡(1/M0)−Aβ1.

   Integrating Eq. (2.9) with respect to y between y = 0 and 1 and using Eq. (3.6), we get

   v(1)=v1=−∫01∂u∂xdy=∂2p∂x2N.(3.7)

   where N=2M03{1−cosh⁡(1/M0)}sinh⁡(1/M0)+16(1+M02)+A2+B.

   By using Eqs. (3.5) and (3.7) then Eq. (3.4) becomes

   ∂η∂t+N1∂η∂x+1B∂3η∂x3∂η∂x=−N∂2η∂x2+1B∂4η∂x4.(3.8)

   where N1=C1Cosh(1/M0)+C2sin⁡h(1/M0)+12+M02+A+B. Eq. (3.8) is a fourth order nonlinear ordinary differential equation and now let us analyzes the interface evolution by this equation. The process described here is quite different from a process wherein the film is bounded by a fluid with moving boundaries instead of porous layer discussed by Babchin et al. [10]. Eq. (3.8) is not amenable to solve by using analytical method and hence it can be solved numerically using 4th order central differences in space and time with periodic boundary conditions as explained below. For time- integration of Eq. (3.8), we use Adams-Bashforth predictor and Adams-Moulton corrector methods of 4th order [19]. For the spatial derivatives are described by the following central difference formulae of fourth order accuracy

   ∂η∂x→112Δxηi−2−8ηi−1+8ηi+1−ηi+2.(3.9)

   ∂2η∂u2→112Δx2−ηi−2+16i−1−30ηi+16ηi+1−ηi+2.(3.10)

   ∂3η∂x3→18Δx3−ηi−3−8ηi−2+13ηi−1−13i+1+8ηi+2−ηi+3.(3.11)

   ∂4η∂x4→16Δx4−ηi−3+12ηi−2−39ηi−1+56ηi−39ηi+1+12ηi+2−ηi+3.(3.12)

Here η (i − 2) stands for value of η at the position x − 2Δ x. The integer i indicates the i^th grid point. The initial condition used in numerical integration is a sine-wave with wave number ℓ and is of the form

Here the amplitude η₀ is assumed to be small. In our numerical computation η₀ = 10^{− 4}, the non-dimensional form and periodic boundary conditions have been applied in x-direction.

4 Results and Discussion

The boundary roughness effects on nonlinear saturation of Rayleigh-Taylor instability (RTI) in a couple-stress fluid have analyzed based on stability of the interface of system between two fluids. The governing fourth order nonlinear ordinary differential Eq. (3.8) is solved numerically on the interval 0<x≤π2ℓ with periodic boundary conditions using Adams-Bashforth predictor and Adams-Moulton corrector method for space and time. This method involves the approximation of spatial derivatives with central finite difference scheme for time step and consequences are depicted in Figs. 2–7. The process of study is to explain the distortion of wave profile i.e. η(x, t). It can also be represented in the following words ( in our interface interpretation of equation); every point of interface moves in horizontal x-direction with velocity which is proportional to that point and elevation η . Figures 2–4 predicts the spatial structure of growth rate of interface in terms of n_NL − nNL,Max=1η∂η∂t−1ηmax∂ηmax∂t at an early stage (initial line) before instability occurs for wave number ℓ = 0.75 and other parameters defined earlier. Also, Figures 2–4 shows six waves are contained in interval for the instance of small wave numbers and peaks reflect the position of wave modes. Figure 5 and 6 predicts that whether the full numerical solutions which deviate from a simple harmonic behaviour of η, namely η = η₀ sin(ℓ x) for ℓ = 0.75, where only a portion of the wave is shown.

Fig. 2

Shape of interface growth rate shape of interface growth rate for a roughness parameter β₁ = 0, 3.3 × 10⁻³, 3.3 × 10⁻² with fixed values of M₀ = 0.3, B = 0.02 and β = 2. (from top to bottom)

Fig. 3

Different values of M₀ = 0.2, 0.25, 0.3, and 0.35 with fixed values of B = 0.02, β = 2, ℓ = 0.75 and β₁ = 3.3 × 10^{− 6}.

Fig. 4

Different values of B = 0.01, 0.02, 0.03 and 0.04 with fixed values of β = 2, M₀ = 0.3, ℓ= 0.75 and β₁ = 3.3 × 10^–6.

Fig. 5

Comparison between full numerical wave number of ℓ̂ = 0.75 and three roughness. Parameters β₁ = 0, 3.3 × 10^{– 3}, 3.3 × 10^{– 2} with fixed values of M₀ = 0.3, B = 0.02 and β = 2. (from top to bottom)

Fig. 6

Comparison between full numerical wave number of ℓ̂ = 0.75 and the couple stress. Parameter M₀ = 0.2, 0.25, 0.3, 0.35 with fixed values of β = 2, B = 0.02 and β₁ = 3.3 × 10^{– 6}. (from top to bottom)

Fig. 7

Comparison between full numerical wave number of ℓ̂ = 0.75 and the Bond number B = 0.01, 0.02, 0.03, 0.04with fixed values of M₀ = 0.3, β₁ = 3.3 × 10^{– 6}and β = 2. (from top to bottom)

Figures 5(a)–5(c) depicts that spatial structures of interface for t>0 and Figures 2(a)–2(c) describes interface at t = 0. Initial profile is symmetric in the interval 0 ≤ x ≤ 16π . But for t > 0, we notice that symmetry is lost, since the maximum moves farther from one of zeros and closer to other one. However, this process of steepening of forward faces of interface profile does not result in breakup of interface because of the effects of surface tension, roughness parameter and couple-stress parameter. Initially these parameters may have negligible effects but as time progresses these parameters play a vital role in stabilization even though negligible effect in reducing the asymmetry of system by roughness parameter. In Figures 5(a)–5(c) only a portion of the wave is shown as it move from top to bottom and clearly notice that symmetry can be obtained with effect of roughness parameter, hence reduces the growth rate of RTI at interface considerably for t > 0.

In Figures 6(a)–6(d) it is clear that interface profile is symmetric in the interval 0 ≤ x ≤ 16π for initial time t = 0 and as it moves from Figures 5(a)–5(c) (top to bottom), it is observed that symmetry can be recovered slowly for increasing small roughness parameter β₁. Therefore, it has negligible effect on reducing asymmetry of the system because of resistance offered by surface roughness in that process a part of kinetic energy is converted into potential energy.

Also, Figures 4(a)–4(d) shows the interface is symmetric in the interval 0 ≤ x ≤ 16π for t = 0. For t > 0, the initial symmetry can be lost as shown in Fig. 6(a) and increasing Bond number the symmetry can be recovered and results are depicted in Figs. 7(b)–7(d). Hence for suitable choice of Bond number the symmetry can maintained at the interface and hence reduce the growth rate of RTI at interface. It may also be noted that full numerical solution for all the cases is not possible due to limitation of numerical scheme.