Stability analysis of an interface between immiscible liquids in Hele-Shaw flow in the presence of a magnetic field

Applied Mathematical Modelling 75 (2019) 572–588

Contents lists available at ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

Stability analysis of an interface between immiscible liquids in Hele-Shaw flow in the presence of a magnetic field Yuri Zeniti Sinzato, Francisco Ricardo Cunha∗ Department of Mechanical Engineering, Laboratory Microhydrodynamic and Rheology - VORTEX, University of Brasilia, Brasilia DF 70910 900, Brazil

a r t i c l e

i n f o

Article history: Received 22 April 2018 Revised 18 May 2019 Accepted 23 May 2019 Available online 4 June 2019 Keywords: Hydrodynamic instability Hele-Shaw cell Magnetic finger-interface Stability analysis Asymptotic solution

a b s t r a c t Hydrodynamic instabilities may occur when a viscous fluid is driven by a less viscous one through a porous medium. These penetrations are common in enhanced oil recovery, dendrite formation and aquifer flow. Recent studies have shown that the use of magnetic suspensions allow the external control of the instability. The problem is nonlinear and some further improvements of both theory and experimental observations are still needed and continue being a current source of investigation. In this paper we present a generalized Darcy law formulation in order to examine the growth of finger instabilities as a magnetic field is applied to the interface between the fluids in a Hele-Shaw cell. A new linear stability analysis is performed in the presence of magnetic effects and provides a stability criterion in terms of the non-dimensional physical parameters of the examined flow and the wavenumber of the finger disturbances. The interfacial tension inhibits small wavelength instabilities. The magnetic field contributes to the interface stability for moderate wavelength as it is applied parallel to the liquid-interface. In particular, we find an explicit expression, as a function of the susceptibility, for a critical angle between the interface and the magnetic field direction, in which its effect on the interface is neutral. We have developed a new asymptotic solution for the flow problem in a weak nonlinear regime. The first correction captures the second order nonlinear effects of the magnetic field, which tends to align the fingers with the field orientation and have a destabilizing effect. The analysis predicts that the non-linear effects at second order can counterbalance the first order stabilizing effect of a parallel magnetic field which results in a loss of effectiveness for controlling the investigated finger instabilities. The relevant physical parameters for controlling these finger instabilities are clearly identified by our non-dimensional analysis. © 2019 Elsevier Inc. All rights reserved.

1. Introduction A common problem in water pumping oil recovery methods is the hydrodynamic instabilities that occur in the interface between water and oil. Since the water is much less viscous, it penetrates through the oil, reducing drastically the sweeping efficiency of the process. Saffman and Taylor [1] and Chouke et al. [2] were the first to observe this phenomena, usually denoted as Saffman–Taylor fingers or viscous fingers. They developed a linear stability analysis for a small initial perturbation, and found a stability criteria for which the interface is unstable. Experiments were performed in a Hele-Shaw cell, showing ∗

Corresponding author. E-mail address: [email protected] (F.R. Cunha).

https://doi.org/10.1016/j.apm.2019.05.049 0307-904X/© 2019 Elsevier Inc. All rights reserved.

Y.Z. Sinzato and F.R. Cunha / Applied Mathematical Modelling 75 (2019) 572–588

Nomenclature Latin symbols a initial amplitude b gap between plates B density of magnetic flux f interface disturbance g gravity acceleration H magnetic field i imaginary unit I identity matrix k wave number Kp porous medium permeability K Boltzmann constant m intensity of particle magnetic dipole M local magnetization p pressure t time T absolute temperature u local average fluid velocity U equilibrium velocity x, y Cartesian coordinates Greek symbols  porosity ζ amplification factor η dynamic viscosity θ magnetized field angle κ curvature λ wave lenght μo vacuum magnetic permeability ρ density σ interfacial tension

stress tensor φ velocity potential ϕ volume fraction χ magnetic susceptibility ψ magnetic potential Physical non-dimensional numbers α magnetic parameter Ca capillary number ε asymptotic solution parameter NA Archimedes number ηA,B viscosity ratio Rm magnetic Reynolds number Subscripts A lower fluid B upper fluid c critical condition n normal component o cut-off condition x, y x,y-component Superscripts o equilibrium condition  disturbed condition R real part I imaginary part

573

574

Y.Z. Sinzato and F.R. Cunha / Applied Mathematical Modelling 75 (2019) 572–588

the growth of small disturbances, leading to the predominance of a single finger with half the channel width. A Hele-Shaw flow is a two dimensional analogue of porous medium single-phase flow, but not a perfect analogue in the case of viscous fingering phenomena. While a continuous interface is present in a Hele-Shaw cell, the flow in porous medium has no clear interface and the propagation of the meniscus has effects on the phenomena. Nevertheless, the Hele-Shaw analogy allows us to observe the effects of the various forces acting on the instability. A very good review on viscous fingering was carried out by Homsy [3]. Similar phenomena also occur in dendrite growth, flame propagation and aquifer flows [4]. The control and stabilization of this type of instabilities has been subject of recent research. Chen et al. [5] simulated numerically the action of centrifugal and Coriolis forces on a radial viscous fingering, detecting important changes on the morphology of the arising patterns. Fontana and Miranda [6] applied a weakly nonlinear approach to investigate interfacial pattern formation in a lifting HeleShaw cell. Al-Housseiny and Stone [7] found that the presence of a depth gradient can alter the stability of the interface offering opportunities to control and tune fingering instabilities. The optimal control of the pumping rate can also be used to prevent the phenomena, as proposed by Dias et al. [8]. Recent studies investigated the effect of chemical reactions on the interface stability [9,10]. The use of magnetic fluids to inhibit viscous fingering was first proposed by Rosensweig et al. [11,12]. They developed a linear stability analysis considering the action of an external magnetic field. The showed that a field parallel to the interface reduces the growth rate of the finger formation. Research into magnetic fluids has been extensive over the past few decades, leading to several commercial applications. This subject is largely discussed in a paper review by Rinaldi et al. [13]. Moreover, research on the use of magnetic fluid in flow control include: the process of magnetic separation under the action of an external magnetic field in the capture of oil spills [14]; drag reduction in capillaries due to the effect of a magnetic field gradient in a pressure-driven flow [15] and the behaviour of a magnetic liquid–liquid interface formed between two vertical flat plates in response to vertical magnetic fields [16]. Mostaghimi et al. [17] showed that the use of a electrically conducting fluid under the action of a external magnetic field can also inhibit the growth of viscous fingers. Some numerical simulation of magnetic suspension instability and velocity fluctuations of magnetic particles have been carried out more recently by Gontijo et al. [18]. The nonlinear dynamics of an oscillating bubble immersed in a magnetic fluid has been explored in a paper by Malvar et al. [19]. Qin and Chadam [20] obtained analytical solutions for the non-linear regime of one-phase thermal convection of ferrofluid in porous medium. The formation of the viscous fingers is essentially nonlinear and some further theoretical studies and experimental investigations are still necessary for understanding this complex phenomenon in the presence of magnetic effects. Some advances were made by Verma and Rajput [21], who obtained an analytical solution for a normal magnetic field gradient under the non-linear regime. In addition, Pacitto et al. [22] and Flament et al. [23] have explored numerically and experimentally the fully developed fingers under both perpendicular and parallel magnetic fields. During the initial finger formation, when regime is linear, a parallel field can decrease the finger growth rate. However, during the transition to the non-linear regime, it is not settled in the current literature whether the nonlinear effects imposed by the presence of a magnetic field do or do not contribute to the attenuation of the fingers’ amplitude or stability of the material interface between the viscous fluids. In the present work, we investigate theoretically the stability and interface shape transitions as a magnetic fluid drives an immiscible and more viscous fluid through a Hele-Shaw cell, under the action of an oblique magnetic field. We shall try to avoid the introduction of any parameters that do not have a clear physical meaning and are not calculable, at least in principle. A linear stability analysis is performed in order to identify the amplification factor as a function of the nondimensional physical parameters based on the appropriate characteristic scales of the flow. We wish to emphasize that the linear analysis presented in the paper by Rosensweig et al. [11] has been extended here by developing a new second order asymptotic solution in order to capture at least a weakly non-linear regime. After this, it is already possible to characterize finger transition from small to large amplitude. The theoretical approach is similar to that of Albernaz and Cunha [24] in the context of bubble dynamics. We have examined how the nonlinear effects, specially those related to the magnetic force, can affect the dynamical behavior of the finger in terms of the growth or decay rate of the disturbances amplitude and shape. An interesting finding of this work has been to show that non-linear second order effects on finger instability, in contrast with the linear prediction, have surprisingly a destabilizing effect even in the presence of a magnetic parallel field. An immediate motivation of this theoretical work has been the practical application related to the control of non-linear hydrodynamic perturbations that might occur in the oil-water interface as applied to water pumping oil recovery. These instabilities can drastically reduce the sweeping efficiency of the oil recovery process. In the present work, the relevant physical parameters for controlling these finger instabilities have also been clearly identified by our scalings and non-dimensional analysis.

2. Mathematical modeling A two dimensional porous medium can be modeled by a Hele-Shaw cell, assuming a parabolic velocity profile between plates (Fig. 1). Its permeability is given by K p = b2 /12, where b is the gap between the plates. Consider the upper fluid (fluid B) as the driven one, without magnetization. The lower fluid (fluid A) is a magnetic fluid injected in the porous medium. In a stationary reference system moving with the equilibrium position, the interface between the fluids is given by:

yi = f (x, t ),

(1)

Y.Z. Sinzato and F.R. Cunha / Applied Mathematical Modelling 75 (2019) 572–588

575

Fig. 1. Flow of the two immiscible fluids in a Hele-Shaw cell.

where f(x, t) is a small disturbance. The velocity of the equilibrium position of the interface is U. Consider an uniform applied magnetic field oblique to the interface, with its Cartesian coordinates given by:

Hx = H cos(θ )

and

Hy = H sin(θ ),

(2)

where H is the magnetic field magnitude and θ is its angle with the x axis. In the equilibrium condition, the magnetic field is uniform. As the interface is disturbed, a magnetic field gradient shall appear, leading to the action of magnetic forces on the fluid. 2.1. Magnetostatics In the absence of electrical field, polarization effects and current density flux, Maxwell’s Equations reduce simply to [25]:



∇ ·B=0 , ∇ ×H=0

(3)

where B is density of magnetic flux and H is the magnetic field, which are related by:

B = μo ( H + M ) ,

(4)

where μo is the vacuum magnetic permeability and M is the local magnetization. Combining the Eqs. (3) and (4) and writing the magnetic field as H = ∇ψ , where ψ is the magnetic potential, results that:

∇ 2 ψ = −∇ · M.

(5)

In this work we consider a superparamagnetic model [12]:

M = χ ( H ) H,

(6)

where χ (H) is the magnetic susceptibility, which is a function of the magnetic field magnitude. Substituting (6) into (5), making the variables non-dimensional by choosing appropriate typical scales of the problem and, after some algebraic manipulation, we find

⎧ ∗ ∗ ∗ ∗ ∗ ⎨∇ ∗2 ψ ∗ = − dχ ∇ ψA · ∇ |∇ ψA | A ∗ dH 1+χ ⎩ ∗2 ∗ ∇ ψB = 0

,

(7)

where ∇ ∗ = b∇ , H ∗ = H/Ms and ψ ∗ = ψ /bMs . Notice that for the fluid B (non-magnetic), χ = 0. For a dilute magnetic suspension, χ (H) is given by:

χ (H ) = L (α )

1 , with H∗

α=

mMs H ∗ . KT

(8)

Here L(α ) = coth α − 1/α denotes the Langevin function as described in [12] which depends on the non-dimensional magnetic parameter α , Ms = nm is the saturation magnetization of the fluid, n is the number of magnetic particle per unit of volume (density number), m is intensity of the particle magnetic dipole (the same for all particles), K is the Boltzmann constant and T is the absolute temperature of the fluid. The initial susceptibility χ o is the limit of χ when α → 0. In this limit L(α ) = α /3 and the susceptibility is the constant χ = mMs /(3KT ).

576

Y.Z. Sinzato and F.R. Cunha / Applied Mathematical Modelling 75 (2019) 572–588

2.2. Generalized Darcy’s Law In creeping flows through porous media, inertia effects may be neglected when compared to viscous effects. Under these conditions, the force balance per unit of volume in the fluid is given by Darcy’s Equation, with additional terms to consider gravitational and magnetic forces:

∇p = −

η Kp

u + ρ g + μoM · ∇ H,

(9)

where p is the pressure, η is the dynamical viscosity, Kp is the medium permeability, u is the local average fluid velocity, ρ is the fluid density, g is the gravity acceleration. Substituting (6) into (9) leads to:

∇p = −

η Kp

u − ρ g + μoχ H · ∇ H.

(10)

For a vertical orientation of the medium, the gravitational field can be written as g = −∇ (gy ). Now, applying the identity H · ∇ H = 1/2∇ (H · H ) − H × (∇ × H ) and Eq. (3) results that:

∇p = −

η Kp

u − ∇ (ρ gy ) + χ ∇

1

The last term can be rewritten as:

χ∇

1 2

μo H

2



 = ∇ μo

H 0

2

 μo H 2 .

(11)

χ HdH .

(12)

Now, isolating the velocity field, it can be described by a velocity potential φ :

φ=

Kp

η



−p − ρ gy + μo



H 0

χ HdH .

(13)

We define the non-dimensional variables φ ∗ = φ /Ub, p∗ = pb/ηBU , y∗ = y/b, K p∗ = K p /b2 , where ηB is the fluid B viscosity. Therefore, the non-dimensional velocity potentials are given by:

⎧ 

H∗ A ρA gb2 ∗ ⎪ ∗ ∗ ∗ ∗ ∗ ⎪ y + Rm χ H dH ⎨φA = Kp −pA − ηB U 0 

, ρB gb2 ⎪ ⎪ ⎩φB∗ = ηA,B Kp∗ −p∗B − η U y∗ B

(14)

where ηA,B (Ratio of viscosity), Ca (Capillary Number) and Rm (Magnetic Reynolds Number) represent the non-dimensional numbers ηA /ηB , ηB U/σ and μoMo2 b/ηBU . Combining the incompressibility condition ∇ · u = 0 with (13) implies that the velocity potentials must satisfy Laplace’s Equation:

 ∇ 2 φA = 0 . ∇ 2 φB = 0

(15)

2.3. Free boundary conditions Consider a reference system fixed at the interface:

Y = y − f (x, t ).

(16)

Since the interface is a material boundary, the following kinematic condition holds:

DY ∂ f Ub =− + ∇ φ · ∇ Y = 0, Dt ∂t 

(17)

with ∇ φ = ∇ φA = ∇ φB at the interface and ∇ Y = (− fx , 1 ). Making the expression (17) non-dimensional, results that:

⎧ ∗ ∗ ∂φA ∂ f ∗ ∂φA∗ ∂f ⎪ ⎪ = − ⎨ ∂t∗ ∂ y∗ ∂ x∗ ∂ x∗ y= f ∗

∗ ∂φB ∂ f ∗ ∂φB∗ ∂f ⎪ ⎪ ⎩ ∂ t ∗ = ∂ y∗ − ∂ x∗ ∂ x∗ y= f

,

(18)

where x∗ = x/b, t ∗ = tU/b e f ∗ = f /b. The magnetic field must satisfy the following boundary conditions at the interface:



ˆ · (B∗A − B∗B ) = 0 n ˆ × (H∗A − H∗B ) = 0 n

,

(19)

Y.Z. Sinzato and F.R. Cunha / Applied Mathematical Modelling 75 (2019) 572–588

577

ˆ is the unitary vector normal to the interface pointing out to the upper fluid. Another boundary where B∗ = B/Ms and n condition is that the effects of the disturbance must vanish at y → ± ∞. Therefore, the potentials must equal the equilibrium condition:



and,



limy→∞ ∇ ∗ φB∗ = 0 limy→−∞ ∇ ∗ φA∗ = 0

,

(20)

limy→∞ ∇ ∗ ψB∗ = Ho∗ B limy→−∞ ∇ ∗ ψA∗ = Ho∗ A

,

(21)

where Ho∗ is the non-disturbed non-dimensional magnetic field. 2.4. Traction jump at the interface The stress tensor in the presence of a magnetic field is given by:

1 2

= −pI − μo|H|2 I + HB.

(22)

According to Young-Laplace Equation, the traction jump at the interface is:

ˆ, ( A − B ) · nˆ = (2κσ ¯ )n

(23)

where κ¯ is the mean curvature and σ is the interfacial tension between the fluids. Substituting (4), (6) and (19) into (23) and making the variables non-dimensional gives the pressure jump on the interface:

( p∗B − p∗A ) =

∗2 Rm χ 2 HAn 2κ¯ ∗ + , Ca 2

(24)

ˆ · HA and the non-dimensional curvature is κ¯ ∗ = κ¯ b. Combining (14) and (24), results that: where HAn = n

−

φB∗ − ηA,B φA∗ K p∗2

2κ¯ ∗ = NA y + + Rm Ca ∗



HA∗

0

χ H dH + ∗

∗

∗2 χ 2 HAn



,

2

(25)

where NA is the Archimedes Number given by NA = (ρB − ρA )gb2 /ηBU . The interface average curvature is given by κ¯ = 0.5(∂ 2 f /∂ x2 )/(1 + (∂ f /∂ x )2 )3/2. The curvature between the plates is neglected, since it is constant for every x. Therefore, it does not affect the interface stability. The boundary condition must be evaluated on y∗ = y∗i = f ∗ (x, t ). Rewriting the traction jump condition (omitting the ∗ ), we obtain:



−

φB − ηA,B φA

Kp

= y= f





∂ 2 f /∂ x2



2 3/2

Ca 1 + (∂ f /∂ x )

+ NA f + Rm

HA 0

χ HdH +

2 χ 2 HAn

2

.

(26)

y= f

The interface shape is completely defined by Eqs. (7), (15), (18)–(21) and (26) combined with the initial condition. At all interface boundary conditions, there are explicit and implicit non-linearities in this problem, since the interface itself is the solution. 3. Linear stability analysis In the absence of any disturbance, the velocity potentials that satisfy the set of Eqs. (15), (18) and (20) are:



φA (x, y, t ) = GA (t ) , φB (x, y, t ) = GB (t )

(27)

where GA (t) e GB (t) are arbitrary functions that must satisfy the traction jump condition, given in Eq. (26). Also, the absence of disturbance implies that there are no magnetic field gradient. The uniform applied magnetic field Ho must satisfy ˆ = eˆ y . Then, Eq. (19). For a plane interface, n



BoAy − BoBy = 0 . o o HAx − HBx =0

(28)

We investigate now the interface stability for a small amplitude disturbance f(x, t):

f (x, t ) = aeζ t+ikx ,

(29)

where k is the wave number, ζ is the amplification factor and a is the initial amplitude. From this place, all variables will be non-dimensional and the ∗ will be omitted for simplicity. The linearized kinematic boundary condition at the interface is:

∂f ∂φ = , ∂t ∂ y y=0

(30)

578

Y.Z. Sinzato and F.R. Cunha / Applied Mathematical Modelling 75 (2019) 572–588

where φ is evaluated at the equilibrium position y = 0. The velocity potentials can be written in terms of an equilibrium potential φ o and a disturbed potential φ  , namely:

φ = φo + φ.

(31)

Solving Laplace’s equation for φ  and applying the linearized boundary condition gives:

⎧ ⎨φ  (x, y, t ) = aζ eζ t+ky+ikx A k . ⎩φ  (x, y, t ) = − aζ eζ t−ky+ikx B

(32)

k

In the same way, the interface perturbation leads to a magnetic field gradient. Therefore, the magnetic potentials can be written as:

ψ = ψ o + ψ .

(33)

Linearizing Eqs. (7), (19), (21) and (28) and solving for ψ  gives (see details in Appendix A):

 ψA (x, y, t ) = mA eζ t+kβ3 y+ikx , ψB (x, y, t ) = mB eζ t−ky+ikx

(34)

where β 3 , mA and mB are given by (56) and (58), available in the appendix of this work. Now, the linearized traction jump condition is:



−

φB − ηA,B φA



=

Kp

y=0

1 − ηA,B + NA Kp



1 f+ Ca





∂ψA ∂ψA ∂ψA ∂ f ∂2 f 2 + Rm χ Hx + Hy + χ Hy − HH (35) . ∂x ∂y ∂y ∂ x x y y=0 ∂ x2

Here we write H = HAo for simplicity. Terms that are not x-dependent were omitted, since they do not affect the interface stability. Applying Eqs. (2), (29), (32) and (34) to (35) provides a dispersion relation for ζ :



ζ=

Kp k3 kC1 + k2C2 − , (1 + ηA,B ) Ca

(36)

where:



⎧ (1 − ηA,B ) ⎪ ⎪ + NA , ⎨C1 = K p

 ⎪ ⎪ ⎩C2 = Rm χ 2 H 2 sin2 (θ ) −



1 . 1 + β3 (1 + χ )

In the linear stability analysis we consider only the real part of ζ , since the imaginary part is related to the lateral propagation of the disturbance, therefore it does not affect the interface stability. The interface will be stable if ζ R < 0 and unstable ζ R > 0 (the superscript indicates the real and imaginary parts). This result was first obtained by Rosensweig et al. [11] in a dimensional form. Fig. 2 represents the amplification factor as a function of the wave number under several conditions of magnetic field and initial magnetic susceptibility. Eq. (36) shows that the amplification factor is given by a third degree polynomial in k. For small wave numbers (i.e. large wave lengths), the term kC1 dominates the left hand side of Eq. (36). The coefficient C1 expresses the balance between viscous and gravitational forces. If the lower fluid is more dense (NA < 0), the gravitational force will tend to stabilize the interface. On the other hand, if the upper fluid is more viscous (ηA,B < 1), the viscous force will amplify the disturbance. When C1 > 0, the viscous force dominates the dynamics and the disturbance will grow exponentially (i.e. viscous fingers). For high wave numbers (i.e. s mall wave lengths), the capillary force (proportional to k3 ) dominates over viscous forces, stabilizing small scale disturbances. For intermediate wave numbers, the magnetic force becomes relevant to the disturbance growth. If the magnetic field is perpendicular to the interface (θ = π /2 or −π /2), the instability will be intensified, but if the field is parallel (θ = 0 or π ), the amplification factor will be reduced. The cut-off wave number is:



Ca ko = C2R + 2



2 C2R



C1 +4 , Ca

(37)

and the critical wave number (i. e. most unstable) is:



Ca kc = C2R + 3



2 C2R



C1 +3 . Ca

(38)

Fig. 3 presents the neutral stability curves of the interface for several wave lengths and a horizontal magnetic field. We investigate the stability as a function of the main magnetic parameters (i.e. magnetic Reynolds Number, initial magnetic susceptibility and magnitude of the field). We consider the limit case as the capillary force is negligible (i.e Ca → ∞). The

Y.Z. Sinzato and F.R. Cunha / Applied Mathematical Modelling 75 (2019) 572–588

579

Fig. 2. Amplification factor as a function of the wave number under various conditions of (a) magnetic field magnitude, (b) initial magnetic susceptibility and (c) field angle. Ca = 1, NA = 0, ηA,B = 0, K p = 1, Rm = 5. The dashed line in (c) indicates the absence of magnetic field. The susceptibility in Eq. (6) was calculated with the Langevin function L (α ) = α −1 − coth(α ).

580

Y.Z. Sinzato and F.R. Cunha / Applied Mathematical Modelling 75 (2019) 572–588

Fig. 3. Neutral stability curve relating the Magnetic Reynolds Number (Rm ), the non-dimensional magnetic field magnitude (H), the initial magnetic susceptibility (χ 0 ) for several wave numbers and a field parallel to the interface. α = 0, Ca → ∞, NA = 0, ηA,B = 0, K p = 1. The susceptibility in Eq. (6) was calculated with the Langevin function L (α ) = α −1 − coth(α ).

Y.Z. Sinzato and F.R. Cunha / Applied Mathematical Modelling 75 (2019) 572–588

581

size of the unstable region increases as the wavenumber decreases, since disturbances of smaller scale are more stable than the larger ones. Disturbances of the interface of small wavelengths are quickly attenuated by magnetic effects. For an interface of arbitrary shape, in which a wide range of wave lengths is present, those modes with k ≈ kc will grow faster and dominate. Therefore, the finger growth will be given by the maximum amplification factor ζkc (i.e amplification factor for the critical wave number). Fig. 4a and b presents ζkc as a function of the magnetic field magnitude for various field angles and initial magnetic susceptibilities. Notice that when H → ∞, the magnetic fluid reaches the saturation condition. In this limit, the amplification factor can only be changed by the Magnetic Reynolds Number, as seen in Fig. 4c. This physical parameter is associated with the fluid saturation magnetization. 4. Regular asymptotic solution As the interface behaves in a nonlinear way with disturbance of finite amplitudes, this finite effect of the finger amplitude and its consequence on the interface shape should be considered in the analysis.Therefore, we shall consider an initial disturbance with k = kc , since it dominates the instability. A regular perturbation method is applied to a weakly nonlinear regime. Firstly, the nonlinear terms are made explicit. Secondly, the variables are re-scaled in order to represent the scale of the nonlinear terms by a non-dimensional parameter. We perform the present nonlinear analysis of Eq. (18). The disturbed velocity potentials are evaluated at y = f . In order to explicit the dependence with f(x, t), we use a first order MacLaurin series as follows:

⎧ ∂φ    2  y= f = ∂φ y=0 + f ∂ φ y=0 + . . . ⎪ ⎪ ⎨ ∂x ∂x ∂ x∂ y

.

⎪  ⎪ ∂φ   ∂ 2 φ   ⎩ ∂φ  + ... y=0 + f y= f = ∂y ∂y ∂ y2 y=0

(39)

Applying the expansion (18) to (39) and neglecting higher order terms, it results that:



⎧ ∂φA ∂ f ∂φA ∂ 2 φA ∂ f ⎪ ⎪ = − +f ⎪ ⎨ ∂t ∂y ∂x ∂x ∂ y2 y=0 . 

⎪ ∂φB ∂ f ∂φB ∂ 2 φB ∂f ⎪ ⎪ ⎩ ∂ t = ∂ y − ∂ x ∂ x + f ∂ y2 y=0

(40)

Now, consider the initial condition of the interface f (x, 0 ) = a cos(kc x ) and define λc = 2π /kc as the critical wave length. The following variables have order f, φA , φB ∼ a, y, x, t ∼ λc . Therefore, we define the rescaled variables f˜ = f /a, φ˜ A = φA /a, φ˜  = φ  /a, x˜ = x/λc , y˜ = y/λc , t˜ = t/λc and substituting these variables into Eq. (40), we have: B

B



⎧ ˜ ˜ ∂ φA ∂ 2 φ˜ A ∂f ∂ f˜ ∂ φ˜ A ⎪ ˜ ⎪ = − ε − f ⎪ ⎨ ∂ t˜ ∂ y˜ ∂ x˜ ∂ x˜ ∂ y˜2 y=0 .  

 2 ˜ ⎪ ˜ ˜ ˜ ˜ ∂ φ ∂ φ ∂ φ ∂f ∂f B ˜ ⎪ B B ⎪ −ε −f ⎩ ˜= ∂ y˜ ∂ x˜ ∂ x˜ ∂ y˜2 ∂t y=0

(41)

Here the non-dimensional parameter ε = a/λc has the scale of the nonlinear terms. The same analysis has been performed for Eqs. (7), (19) and (26), as presented in Appendix B of this paper. Now, if ε 1, the nonlinear terms are weak ones so that a regular asymptotic solution O(ε ) applies:

f˜(x, t ) = f˜0 (x, t ) + ε f˜1 (x, t ).

(42)

In the same way, the disturbed magnetic and velocity potentials are given by the following asymptotic solutions:

 φ˜  (x, y, t ) = φ˜ 0 (x, y, t ) + ε φ˜ 1 (x, y, t )

ψ˜  (x, y, t ) = ψ˜ 0 (x, y, t ) + ε ψ˜ 1 (x, y, t )

,

(43)

where φ0 , φ1 , ψ0 , ψ1 must satisfy individually Eqs. (7), (15), (20) and (21). Applying the asymptotic solutions to Eqs. (41), (61), (62) and (64) and equating the coefficients of ε , gives two non-homogeneous linear systems of coupled differential equation in f˜0 (x, t ) and f˜1 (x, t ), respectively. The initial condition for the interface is:

f˜(x, 0 ) = cos(2π x˜).

(44)

We specify the initial conditions for each term of the asymptotic solutions as being:



f˜0 (x, 0 ) = cos(2π x˜) f˜1 (x, 0 ) = 0

.

(45)

582

Y.Z. Sinzato and F.R. Cunha / Applied Mathematical Modelling 75 (2019) 572–588

Fig. 4. Maximum amplification factor as a function of the (a) non-dimensional magnetic field magnitude for various initial magnetic susceptibilities, (b) magnetic field magnitude for various field angles and (c) Magnetic Reynolds Number for the fluid saturation condition. Ca = 1, NA = 0, ηA,B = 0, K p = 1. The dashed line in (a) indicates the asymptote of the curves.The susceptibility in Eq. (6) was calculated with the Langevin function L (α ) = α −1 − coth(α ).

Y.Z. Sinzato and F.R. Cunha / Applied Mathematical Modelling 75 (2019) 572–588

583

Fig. 5. Typical bifurcations of the interface shape for a weakly nonlinear case: (a) Initial condition; (b) absence of magnetic field; (c) parallel magnetic field (θ = 0, H = 0.75); (d) oblique magnetic field (θ = π /6, H = 1). - - - Absence of magnetic field. Ca = 1, NA = 0, ηA,B = 0, K p = 1/12, ε = 0.03, Rm = 10, χo = 1 .

584

Y.Z. Sinzato and F.R. Cunha / Applied Mathematical Modelling 75 (2019) 572–588

Fig. 6. Time evolution of the disturbance peak according to the asymptotic solution for various values of the small parameter  in comparison to the zero-field condition (H = 0). Ca = 1, NA = 0, ηA,B = 0, K p = 1/12, Rm = 10, χo = 1, H = 1, θ = 0.

The system in f˜0 (x, t ) correspond to the problem for ε = 0, and its solution is the same found in the linear stability analysis. The system in f˜1 (x, t ) provides the first distortion on the finger shape. Solving the differential equations, substituting into Eq. (42) and returning to the non-scaled variables gives:





f (x, t ) = Re aeζkc t eikc x +

a2 k3c K p Rm H 2 χ 2 e2 ζkc t − eζ2kc t β4 ei2kc x , 1 + ηA,B 2 ζkc − ζ2kc

(46)

where ζkc and ζ2kc are the amplification factors for k = kc and 2kc , respectively, calculated by Eq. (36). In Eq. (46) the real and imaginary parts of ζ are considered. Both terms contribute to the perturbations of the interface. The parameter β4 = β4 (H, χ , θ ), given by (65), is a complex function of the magnetic parameters (see details in Appendix B). This solution is valid for sufficiently small times, when ε| f˜1 | << | f˜0 |. The second order perturbation, associated with the effect of the magnetic field, provides the first degree of distortion in the sinusoidal shape of the disturbance. Fig. 5 presents the time evolution of an unstable interface in three conditions: absence of magnetic field, with a parallel field and in the presence of an oblique field. In the absence of magnetic field, only the first order term is present, and the interface remains sinusoidal, as shown in Fig. 5b. Non-linear terms due exclusively to hydrodynamic effects would contribute at third order only. As the field is applied, a second order term of wavelength 2kc changes the finger shape. This term can be separated into two contributions, a cosine and a sine term:



Re

e2 ζkc t − eζ2kc t β4 ei2kc x 2 ζkc − ζ2kc





= Re







e2 ζkc t − eζ2kc t e2 ζkc t − eζ2kc t β4 cos(2kc x ) − Im β4 sin(2kc x ). 2 ζkc − ζ2kc 2 ζkc − ζ2kc

(47)

The amplitude of the cosine term is positive for any field orientation. Its effect is to sharpen the tip of the finger and flatten the back, as we can see in Fig. 5c. The sine term breaks the finger symmetry with respect to the y-axis, and its amplitude is non-zero for α = 0 and α = π /2. This effect tends to align the finger with the direction of the external magnetic field. The sine term is much smaller than the cosine term, even for an oblique field. Therefore, the symmetry break is not a predominant effect of the asymptotic solution, as can be seen in Fig. 5d. In order to evaluate the practical effect of the second order term in the interface stability, the time behavior of the disturbance maximum peak (fmax ) is examined. In oil recovery process, the higher the peak is, the sooner the injected fluid will reach the extracting well. Consequently, reducing the amount of recovered oil. Fig. 6 presents the time evolution of the disturbance peak for various values of the small parameter  under a parallel field. For  = 0, corresponding to the linear stability solution, the disturbance growth under a parallel field is significantly less compared to the zero magnetic field condition. However as a small but finite value of  is considered, the second order solution surprisingly becomes large compared to the leading order. Under this condition, the degree of instabilities is accentuated while at the same time the amplitude of the disturbances in the interface between the viscous fluids amplified. It is seen that, even for values of  still small (i.e. 1), the finger growth may overcome that one corresponding to the zero magnetic field condition. 5. Conclusions In this work we have successfully examined the problem of the interface stability when a magnetic fluid drives a more viscous fluid under action of an oblique magnetic field through a porous medium. A stability criterion was obtained as a function of non-dimensional physical parameter of the flow. The regime in which the magnetic field balances viscous forces

Y.Z. Sinzato and F.R. Cunha / Applied Mathematical Modelling 75 (2019) 572–588

585

has been identified. The interfacial tension limits the range of wave numbers that lead to the instability, stabilizing small scale disturbs. A magnetic field parallel to the interface can be used to control the stability of the interface. When the disturbance is already unstable, the flow problem quickly evolves to a nonlinear regime. A regular perturbation method was successfully applied to provide an asymptotic solution of second order to the magnetic effect on the porous media flow. This solution even at second order was sufficient to characterize nonlinear magnetic effects on the finger shape caused by the application of an external magnetic field. This second order effects distorted the finger shape. In the case of an oblique field, the finger tends to align slightly with the magnetic field direction. The second order solution has a destabilizing effect and may become large compared to the leading order solution. We have concluded that, although a parallel magnetic field was always a stabilizing effect at first order (i.e. linear regime), reducing the viscous fingering growth, as the transition evolves toward a non-linear regime of the interface disturbances, the second order contribution surprisingly appeared as a destabilizing effect. In particular, the amplitude of the disturbance may be larger than the value predicted by the first order solution, even for  still small. We have argued that the global effect of a parallel magnetic field might also have a destabilizing effect even as weakly non-linear perturbations are produced on the interface by the proper field. Consequently, an increase in oil recovery would be significantly less than that predicted by the standard linear stability analysis [11]. As a final comment, we wish to state that calculations such as those described here serve as a good starting point in which a perturbation method can be applied to explore at lest weakly non-linear effects on finger instability in porous media in the presence of a parallel magnetic field without actually performing numerical simulations. Thus the method itself has the potential to be useful in the present flow problem that has been benefited from the evaluation of the first contribution of nonlinear effects, as opposed to linear studies on finger instability. Thus, the present results put in perspective the question of exploring by future works based on numerical simulations or experimental observations higher order non-linear effects on the finger instability. In this work we also have discussed on the relevant physical parameters that should be applied for controlling the dynamics of finger instabilities in the presence of magnetic fields. The control of this finger instability could be important for preserving the sweeping efficiency by using instead of pure water a dilute suspension of magnetic particles plus water pumping oil recovery processes. Acknowledgments This work was supported by the Brazilian funding agency CNPq - Ministry of Science, Technology and Innovation of Brazil (Grants: 402371/2013-5, 305764/2015-2 and 441340/2014-8). YZS also wish to acknowledge CNPq and FAPDF for supporting him with a scholarship during this work at University of Brasilia-Brazil. Appendix A. Solution for the disturbed magnetic potential We must find the magnetic potentials that satisfy Eqs. (7), (19) and (21). The solution for ψB that satisfy Laplace’s Equation and the boundary conditions is:

ψB (x, y, t ) = mA eζ t−ky+ikx .

(48)

It should be important to note that Eq. (7) for ψA is nonlinear. To simplify it, we consider the small amplitude distur-

bance:

|∇ψ  | << |∇ψ o|.

(49)

Therefore, the following approximations can be made:

(i ) ∇ 2 ψ = ∇ 2 ψ o + ∇ 2 ψ  = ∇ 2 ψ  ;



∂ 2ψ  ∂ 2 ψ  eˆ x ∂ 2ψ  ∂ 2 ψ  eˆ y (ii ) ∇ |∇ ψ| ≈ Hx 2 + Hy + Hx + Hy ; ∂ xy H ∂ xy ∂x ∂ y2 H 

2  2  1 ∂ 2ψ  2∂ ψ 2∂ ψ (iii ) ∇ |∇ ψ| · ∇ψ ≈ Hx + Hy + 2Hx Hy . H ∂ x∂ y ∂ x2 ∂ y2

(50)



(51)

(52)

Applying the approximations to (7), it is reduced to the following linear equation:



2  2  ∂ 2 ψA ∂ 2 ψA ∂ 2 ψA β1 2 ∂ ψA 2 ∂ ψA + = H + H + 2 H H . x y x y ∂ x∂ y ∂ x2 ∂ y2 H2 ∂ x2 ∂ y2

(53)

Here, the susceptibility is a function of the equilibrium magnetic field χ = χ (H o ). We define the modified first and second derivatives of χ as β 1 and β 2 , given by:

β1 = −

H dχ (1 + χ ) dH

and

β2 =

H2 d2 χ . (1 + χ ) dH 2

(54)

586

Y.Z. Sinzato and F.R. Cunha / Applied Mathematical Modelling 75 (2019) 572–588

Solving Eq. (53) for the boundary condition (21), results that:

ψA (x, y, t ) = mA eζ t+kβ3 y+ikx ,

(55)

where β 3 is given by:



β3 =

1 − β1 + iβ1 cos(θ ) sin(θ ) 1 − β1 sin (θ ) 2

.

(56)

Combining (19) and (28), linearizing and making the variables dimensionless gives the interface boundary condition for the disturbed magnetic potentials:

⎧  ∂ψA ∂ψ ⎪ = −ik f χ Hy ⎨ B− ∂x ∂x .   ⎪ ⎩ ∂ψB − (1 + χ ) ∂ψA = −ik f χ Hx ∂y ∂y

(57)

Substituting (29), (48) and (55) into (57), results that:

⎧ aχ ⎪ ⎪mA = 1 + β (1 + χ ) (iHx + Hy ) ⎨ 3 ⎪ ⎪ ⎩mB =

aχ (iHx − β3 (1 + χ )Hy ). 1 + β3 (1 + χ )

(58)

Appendix B. Details of the regular perturbation method Defining the re-scaled variables f˜ = f /a, φ˜ A = φA /a, φ˜ B = φB /a, x˜ = x/λc , y˜ = y/λc , t˜ = t/λc and applying to the Eqs. (7), (18), (19) and (26), results that

 

∂ φ˜ A ∂ 2 φ˜ A ∂ f˜ ∂ f˜ ∂ φ˜ A = −ε − f˜ + O ( ε 2 ); ∂ y˜ ∂ x˜ ∂ x˜ ∂ y˜2 ∂ t˜ y=0  

∂ φ˜ B ∂ f˜ ∂ f˜ ∂ φ˜ B ˜ ∂ 2 φ˜ B = −ε −f + O ( ε 2 ); ˜ ∂ y˜ ∂ x˜ ∂ x˜ ∂ y˜2 ∂t y=0

  2   2 

∂ ψ˜ A ∂ ψ˜ A ∂ f˜ ∂ ψ˜ A ∂ ψ˜ B ∂ ψ˜ B ∂ f˜ ∂ ψ˜ B ∂ f˜ ˜ χ Hx + − (1 + χ ) + ε f˜ − − ( 1 + χ ) f − + O ( ε 2 ) = 0; ∂ x˜ ∂ y˜ ∂ y˜ ∂ x˜ ∂ x˜ ∂ x˜ ∂ x˜ ∂ y˜2 ∂ y˜2 y=0 f˜χ Hy +

 

∂ ψ˜ B ∂ ψ˜ A ψ˜ B − ψ˜ A + ε f˜ − + O ( ε 2 ) = 0; ∂ y˜ ∂ y˜ y=0



2  2  ∂ 2 ψA ∂ 2 ψA ∂ 2 ψA β1 2 ∂ ψA 2 ∂ ψA + = H + H + 2 H H x y x y ∂ x∂ y ∂ x2 ∂ y2 H2 ∂ x2 ∂ y2  2 

   2  ˜ ∂ ψA 3 ∂ψA ∂ψA ∂ψA ∂ ψA ε + 4 β1 2 H − Hy3 + H H2+ ∂ x˜∂ y˜ x ∂ y ∂x ∂ x ∂ x2 x H  2 

 

∂ ψA ∂ 2 ψA ∂ψA ∂ψA ∂ 2 ψA ∂ψA 2 + H H + − Hy Hx H − H ∂ y ∂ y2 y ∂x y ∂y x ∂ x2 ∂ y2 



2  2  ∂ψA ∂ψA ∂ 2 ψA 2 ∂ ψA 2 ∂ ψA −β2 Hx + Hy Hx + 2Hx Hy + Hy ; ∂x ∂y ∂ y∂ x ∂ x2 ∂ y2 −

 



∂ φ˜ A ∂ φ˜ B 1 − ηA,B φ˜ B − ηA,B φ˜ A + ε f˜ − ηA,B = + NA f˜ ∂ y˜ ∂ y˜ Kp y=0  



 2 ˜ ∂ ψ˜ A ∂ f˜ ∂ ψ˜ A ∂ ψ˜ A 1 ∂ f 2 + + Rm χ Hx − H H + χ Hx + Hy Ca ∂ x˜2 ∂ y˜ ∂ x˜ x y ∂ x˜ ∂ y˜  2 2 ∂ 2 ψ˜ A 1 ∂ ψ˜ A ∂ f˜ Hx2 − Hy2 ∂ f˜ ∂ ψ˜ A ∂ f˜ ∂ ψ˜ A +ε − − Hx − Hy + f˜ H 2 ∂ y˜ ∂ x˜ 2 ∂ x˜ ∂ y˜ ∂ x˜ ∂ x˜ ∂ y˜2 y

1 Kp



(59)

(60)

(61)

(62)

(63)

Y.Z. Sinzato and F.R. Cunha / Applied Mathematical Modelling 75 (2019) 572–588





587

 2

2

∂ 2 ψ˜ A ∂ 2 ψ˜ A ∂ ψ˜ A 1 ∂ ψ˜ A ˜ +χ H + f H + + y x ∂ y˜∂ x˜ 2 ∂ x˜ ∂ y˜ ∂ y˜2 

Hy2 ∂ 2 ψ˜ A ∂ 2 ψ˜ A − f˜(1 + χ )β1 2 Hx + Hy + O ( ε 2 ), ∂ x˜∂ y˜ H ∂ y˜2 y=0 f˜

(64)

where the dimensionless parameter ε = a/λc has the scale of the nonlinear terms. The function β 4 (H, χ , θ ) described in Eq. (47) is associated with the magnetic parameter and is a complex number:

β4 (H, χ , θ ) = β4R (H, χ , θ ) + iβ4I (H, χ , θ ),

(65)

where, β4R and β4I are the real and imaginary parts of β 4 . For small magnitude fields, the magnetic susceptibility is nearly constant and the real and imaginary parts are given by:

χ 2 (cos (θ ) )2 − 2 χ (cos (θ ) )2 − 2 χ 2 − 4 (cos (θ ) )2 − 4 χ − 2 ; (χ + 2 )2

β4R = −

1 2

β4I = −

χ sin (θ ) cos (θ )(χ + 1 ) . (χ + 2 )2

(66)

(67)

For a general function of χ with H, and for small angles θ :

β4R =

 

1

2 −2 β5 + 1

2β5 (1 + (χ + 1 )β5 ) (1 − (χ + 1 )β5 ) 2

4

2



2

 4 β5 (χ + 1 )3

  χ − 5 β1 2 − 7 β1 β5     2 2 2 −β5 2 χ 2 β5 + 2 χ + 3 β1 (1 + (χ + 1 )β5 ) −χ 2 β5 + 3 β1 − 2   2   2 2 −2 8 β1 − 3 β1 + 2 β5 (χ + 1 )3 + 2 β1 + 35 β1 − 3 β2 + 2 χ 2    2 2 + 15 β1 + 25 β1 − 3 β2 + 4 χ + 13 β1 + 8 β1 ; +(3 β1 − β2 − 4 )χ 2 + 5 β1 − 2 β1 − β2 − 4

β4I =

2



sin(θ ) 2β5 (1 + (χ + 1 )β5 ) (1 − β5 (χ + 1 ) ) 3

4



2

−4 β1 + 9 β1 − β2 − 6

+4 β1 β5 (χ + 1 )3 + −12 β1 + 19 β1 − β2 − 8 4

where

β5 =



1 − β1

2

2

  χ + 10 β1 β5 2 ,

(68)

 2 χ (69)

(70)

References [1] P.G. Saffman, G.I. Taylor, The penetration of a fluid into a porous medium or Hele-Shaw cell containing a more viscous liquid, in: Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 245, The Royal Society, 1958, pp. 312–329. [2] R.L. Chuoke, P. van Meurs, C. van der Poel, The instability of slow, immiscible, viscous liquid-liquid displacements in permeable media, Pet. Trans. AIME 216 (1959) 188–194. [3] G.M. Homsy, Viscous fingering in porous media, Ann. Rev. Fluid Mech. 19 (1) (1987) 271–311. [4] J.S. Langer, Dendrites, viscous fingers, and the theory of pattern formation, in: Dynamics and Patterns in Complex Fluids, Springer, 1990, pp. 190–193. [5] C.Y. Chen, C.W. Huang, H. Gadêlha, J.A. Miranda, Radial viscous fingering in miscible Hele-Shaw flows: A numerical study, Phys. Rev. E 78 (1) (2008) 016306. [6] J.V. Fontana, J.A. Miranda, Finger competition in lifting Hele-Shaw flows with a yield stress fluid, Phys. Rev. E 88 (2) (2013) 023001. [7] T.T. Al-Housseiny, H.A. Stone, Controlling viscous fingering in tapered Hele-Shaw cells, Phys. Fluids 25 (9) (2013) 092102. [8] E.O. Dias, E. Alvarez-Lacalle, M.S. Carvalho, J.A. Miranda, Minimization of viscous fluid fingering: a variational scheme for optimal flow rates, Phys. Rev. Lett. 109 (14) (2012) 144502. [9] C. Rana, A. De Wit, Control of viscous fingering by chemical reactions, in: Proceedings of the APS Meeting Abstracts, 2016. [10] S. Arai, Y. Nagatsu, P. Shukla, A. De Wit, Stabilization of miscible viscous fingering by chemical reaction decreasing viscosity, in: Proceedings of the APS Meeting Abstracts, 2016. [11] R.E. Rosensweig, M. Zahn, T. Vogler, Stabilization of fluid penetration through a porous medium using magnetizable fluid, in: Thermomechanics of Magnetic Fluids: Theory and Applications, 1, 1978, pp. 195–211. [12] R.E. Rosensweig, Ferrohydrodynamics, Dover, 1985. [13] C. Rinaldi, A. Chaves, S. Elborai, X. He, Z. M., Magnetic fluid rheology and flows, Curr. Opin. Colloid. Interface Sci. 10 (2005) 141–157. [14] F.R. Cunha, Y.D. Sobral, Characterization of the physical parameters in a process of magnetic separation and pressure-driven flow of a magnetic fluid, Physica A Stat. Mech. Appl. 343 (2004) 36–64. [15] A.P. Rosa, R.G. Gontijo, F.R. Cunha, Laminar pipe flow with drag reduction induced by a magnetic field gradient, Appl. Math. Model. 40 (2016) 3907–30918. [16] R.G. Gontijo, S. Malvar, Y.D. Sobral, F.R. Cunha, The influence of a magnetic field on the mechanical behavior of a fluid interface, Meccanica (2016) 1–19. [17] P. Mostaghimi, M. Ashouri, B. Ebrahimi, Hydrodynamics of fingering instability in the presence of a magnetic field, Fluid Dyn. Res. 48 (5) (2016) 055504.

588

Y.Z. Sinzato and F.R. Cunha / Applied Mathematical Modelling 75 (2019) 572–588

[18] R.G. Gontijo, F.R. Cunha, Numerical simulations of magnetic suspensions with hydrodynamic and dipole-dipole magnetic interactions, Phys. Fluids 29 (2017) 062004. [19] S. Malvar, R.G. Gontijo, F.R. Cunha, Nonlinear motion of an oscillating bubble immersed in a magnetic fluid, J. Eng. Math. 108 (2018) 143–170. [20] Y. Qin, J. Chadam, A nonlinear stability problem for ferromagnetic fluids in a porous medium, Appl. Math. Lett. 8 (2) (1995) 25–29. [21] A.P. Verma, A.K. Rajput, Instabilities in displacement processes through porous media with magnetic fluid, J. Magn. Magn. Mater. 65 (2-3) (1987) 330–334. [22] G. Pacitto, C. Flament, J.C. Bacri, Viscous fingering in a magnetic fluid. ii. linear hele–shaw flow, Phys. Fluids 13 (11) (2001) 3196–3203. [23] C. Flament, G. Pacitto, J.C. Bacri, I. Drikis, A. Cebers, Viscous fingering in a magnetic fluid. i. radial Hele-Shaw flow, Phys. Fluids 10 (10) (1998) 2464–2472. [24] D.L. Albernaz, F.R. Cunha, Unsteady motion of a spherical bubble in a complex fluid: Mathematical modelling and simulation, Appl. Math. Model. 37 (2013) 8972–8984. [25] R.E. Rosensweig, Ferrohydrodynamics, Courier Corporation, 2013.