Mathematical Problems in Engineering
Volume 2015 (2015), Article ID 439704, 8 pages
http://dx.doi.org/10.1155/2015/439704

Research Article

Two-Phase Nonturbulent Flow with Applications

Darko Mitrović1,2 and Andrej Novak³

¹Department of Mathematics, University of Zagreb, Bijenička Cesta 30, 10000 Zagreb, Croatia
²Department of Mathematics, University of Montenegro, Cetinjski Put bb, 81000 Podgorica, Montenegro
³Faculty of Electrical Engineering and Computing, University of Zagreb, Unska Ulica 3, 10000 Zagreb, Croatia

Received 3 December 2014; Revised 6 March 2015; Accepted 7 March 2015

Academic Editor: Yannis Dimakopoulos

Copyright © 2015 Darko Mitrović and Andrej Novak. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

We model dynamics of two almost immiscible fluids of different densities using the Stokes equations with the Dirac distribution representing the sink or source point. The equations are solved by regularizing the Dirac distribution and then using an iterative procedure based on the finite element method. Results have potential applications in water pollution problems and we present two relevant situations. In the first one, we simulate extraction of a light liquid trapped at the bottom of a pond/lake and, after being disturbed, it rises toward the surface. In the second case, we simulate heavy liquid leaking from a source and slowly dropping on an uneven bottom.

1. Introduction

Questions concerning transport of fluids are essential in environmental and industrial engineering. Among many situations of interest, we will describe here two examples which we will put in the context of pollution problems. In both situations, we have two fluids which are almost immiscible and have slightly different densities. This means that we assume almost sharp interface between the fluids and slow movements.

Such assumptions are often in different applications. Let us mention the problem of surface circulation in Northern Hemisphere lakes and marginal seas. This relates to the residual motion obtained by removing the variability related to tides and other high-frequency phenomena (primarily currents directly driven by wind, seiches, and inertial oscillations). Usually, the flow is cyclonic (i.e., counterclockwise) and it is governed by the equations that we are going to consider in the paper. Let us remark that the phenomenon is explained in [1].

Another notable application of the model to be analyzed here is in solid Earth sciences. To be more precise, recall that the mantle is a part of a terrestrial planet or other rocky bodies large enough to have differentiation by density. The interior of Earth, similar to the other terrestrial planets, is chemically divided into layers. The mantle is a layer between the crust and the outer core. The processes of gravitational instability of two types occur in the Earth’s crust and the mantle: the Rayleigh-Taylor instability of multilayered structure where one of the lower layers is less dense than the upper layers or the Rayleigh-Bernard instability of a single layer heated from below (e.g., [2, 3]). It is important to analyze the instabilities to understand the processes in the deep Earth interior. In the current contribution, the mathematical problem is described by the Stokes and density advection-diffusion equations, similar to those used in geodynamics applications [4–7].

Let us finally remark that environmental issues have already been investigated through the Stokes type equations (see, e.g., [8] and references therein).

The paper is organized as follows.

After the introduction, we provide the governing equations and explain corresponding modelling issues. In Section 3, we describe mathematical and numerical techniques that we are applying. Finally, we present simulations in two different situations.

In the first one, we have a lighter fluid which resides at the bottom of, say, a pond or lake. It is trapped there due to the pressure of surrounding water but if disturbed, it will start rising and, thus, it will cause pollution. To avoid the risk of such a situation, we would like to extract the pollutant from the bottom. If for some reason the pollutant is disturbed, it will start to move up. It is thus reasonable to rise the extraction point together with the rising fluid. To this end, we will simulate situations when the extraction point is fixed and when it moves up.

In the second situation, we will inspect a heavier pollutant leaking from some source (e.g., a ship). We will see how it spreads along through the water.

2. Governing Equations

The first equation that governs mentioned processes is (diffusive) conservation of mass (continuity equations). It has the form

where is the density of the fluid mixture, is the diffusion coefficient which is small (since we assumed that fluids are almost immiscible), and is the velocity of the fluid. The function is a source/sink term and we will model the entry/extraction point by the sink term localized around appropriate coordinates.

Since we are considering slow flow, velocities are given by the Stokes equations

where, as usual, and is the pressure. On the right-hand side we have gravity directed in the -direction and is another external force. Throughout the paper, we will assume that we do not have an external force other than gravity.

Finally, to close the system, we will assume the incompressibility type conditions (as usual, they are obtained by putting in (1)):

Next, we need to choose initial and boundary conditions for and , respectively. We will always assume that the two fluids are separated in different regions which implies that the initial condition has the following form:

where is the Heaviside function, is density of the lighter, and is density of the heavier fluid. The functions are continuous functions such that determines the (separate) interfaces between lighter fluids, for every . For instance, if we are working in the domain and we have the lighter fluid confined in the circle centered in zero and of radius 1, then the initial condition has the form .

Concerning the velocity , we need first to decide on which domain we will solve the equations. We will work on a square denoted by (on which we will distort bottoms in the case of sinking pollutant). Thus, we will have the following general conditions (mixed Dirichlet and Neumann):

where , , , and . By we denote the normal derivative on , while and are appropriate functions to be determined later. Let us briefly describe them.

We will take that the velocity (both horizontal and vertical) on the surface is zero (Dirichlet conditions implying no movement on the surface of the pond). This statement is not true in general since any small disturbance of the interface between two fluids results in a flow. However, the flow is considered to be stable and slow, and we are basically interested in the behavior of the fluids away from the surface, so we can neglect the movements there. Another reason is purely mathematical. Namely, for the other edges, we will take the Neumann conditions

that is, the horizontal velocity is not changed along the edges, and

that is, the edges do not affect the vertical direction of the velocity. In order to be able to solve corresponding elliptic equations, we must have Dirichlet conditions on a part of the boundary.

We will provide details in Section 4.

3. Solution Strategy

We will construct a sequence of approximate solutions which converge along a subsequence toward a weak solution of the problem (1), (2), (3), (4), and (5). Accordingly, fix the time interval , , and split it on even subintervals , . Remark that depends on time variable and space variables .

In order to obtain the sequences of approximate solutions and , we substitute from (4) into the Stokes equations (2), (3) and solve them with boundary conditions (5). Thus, we obtain and and assume that and for . Now, we substitute into the mass conservation equation (1) and take initial conditions (4) to compute the function for .

Accordingly, we have defined the sequences and in the intervals and , respectively. We continue in a similar manner for . Namely, we use to define and and take and for . We then use on the interval to obtain for , and so forth.

As we can see, we are solving initial-boundary value problem for system of Stokes equations and the parabolic equation. According to the results from [9, 10], we know that solutions to the Stokes equations have better regularity than the boundary function and forces given by the function . Similarly, the parabolic equation (1) representing the conservation of mass has regularizing effects with respect to initial data. Moreover, we have linear equation with respect to and , and, therefore, it is enough to prove strong convergence of the sequence .

To explain the latter more precisely, denote by the space of -functions whose Sobolev derivatives up to th order belong to as well. By we denote the space of functions which are continuous with respect to and of class with respect to .

Now, since are at least in for any fixed (see, e.g., [10]) and are continuous with respect to , since they are coefficients in (1), we conclude that is uniformly bounded in . Therefore, the sequence is strongly precompact in ; that is, it contains a subsequence which strongly converges in . Choosing a weakly converging subsequence of and indexed in the same way as the converging subsequence of , we conclude that our scheme converges toward a weak solution to (1), (2), (4), and (5).

Let us remark that the obtained solution is unique which follows from the Schauder fixed point theorem arguments. Details can be found in [11] in a similar situation. We will discuss this purely mathematical part of the paper in more detail elsewhere.

Concerning the numerical part, we will use the mixed finite element method at each step (i.e., in each time interval ) which can be found in [12]. Remark that level-set methods type techniques can be found in [13, 14] in simpler situations.

4. Simulations

In the first part of the section, we consider the problem of extracting a pollutant from the bottom of the lake represented by the square . Let us remark that the typical situation is in the frame of CO₂ sequestration where residual trapping is one of the procedures preventing the gas to go back to the surface. However, if it is disturbed by some force (e.g., an earthquake) it will start moving up. Extreme examples of such situations are, for example, Lake Nyos (Cameroon, Africa). Namely, magma which lies under the lake leaks CO₂ under the water. Part of the CO₂ is dissolved in the water changing it into carbonic acid, and part is trapped at the bottom of the lake, but in the state of unstable equilibrium. In 1986, due to landslides, Lake Nyos suddenly emitted a large amount of CO₂ which caused suffocation of over 1700 people. Other examples include tanks with different pollutants originating from human activities lying at the bottom of lakes or seas. If for some reason the tanks are opened, we will have the situation that we consider here.

Accordingly, assume that the pollutant is in an unstable equilibrium and that it is somehow disturbed. Since it is lighter than the surrounding fluid, it starts moving up.

Initial condition has the form (lighter fluid has the density and the heavier )

while the boundary conditions for the velocities are for the first three sets of simulations (recall that )

implying no movement on the surface for the first three sets of simulations;

implying no velocity change on appropriate edges.

In (1), we have chosen

In the first set of simulations (Figure 1), for the sink term we take a regularization of the -distribution

where is the distribution defined for every by

or, informally speaking,

For the numerical purposes, we replace from (12) by the regularization (see, e.g., [15] for systematic usage of -distribution in numerics)

that is,

(a)

(b)

(c)

(d)

Figure 1: Dynamics of the density distribution with fixed-position sink. After iterations (with the time step ) we have stopped the extraction since the pollutant escaped from the sink point and we see that quite large amount of pollutant escaped. More precisely, we extracted only around 25 percent of the pollutant before we stopped the sink. Here, red, yellow, green, and blue color denote densities 0.35, 0.46, 0.52, and 0.6, respectively. (a) At the initial moment , the pollutant is undisturbed at the bottom. (b) After iterations, pollutant is rising while we are extracting it at the fixed point . (c), (d) After and iterations, we observe further evolution of pollutant.

The constant multiplying the -measure in (12) describes the intensity of the sink. This means that we have a sink at point (i.e., it is a pipe extracting the liquid at point ). If the regularization is not sharp, that means that we have a wider sink (nonlocalized), while the distribution represents an infinitely thin sink. We present simulations that were obtained by implementing the proposed method.

On the next set of simulations (Figure 2), we move the sink together with the gravity. More precisely, we take for the sink

where represents marching forward in time with the time step ( is the first iteration, the second, etc.).

(a)

(b)

(c)

(d)

Figure 2: Dynamics of the density distribution of the pollutant with slowly moving sink. After iterations (with the time step ) we have stopped the extraction since the pollutant escaped from the sink point that was too slow. We extracted around 65 percent of the pollutant. Here, red, yellow, green, and blue color denote densities 0.35, 0.46, 0.52, and 0.6, respectively. (a) At the initial moment , the pollutant is undisturbed at the bottom. (b) After iterations, pollutant is rising while we are extracting it at point moving by the law . (c), (d) At and iterations we see further evolution of pollutant.

After many attempts, we reached the conclusion that the following (linear) motion of the sink will give optimal extraction results:

where, as before represents marching forward in time with the time step ( is the first iteration, the second, etc.).

The results are presented in Figure 3.

(a)

(b)

(c)

(d)

Figure 3: Dynamics of the density distribution of the pollutant moving upwards with the sink moving together with the pollutant (optimal speed). We have stopped the sink after iterations. We extracted around 85 percent of the pollutant before the sink is stopped. Here, red, yellow, green, and blue color denote densities 0.35, 0.46, 0.52, and 0.6, respectively. (a) At the initial moment, the pollutant is undisturbed at the bottom. (b), (c) After and iterations pollutant is rising while we are extracting it at point moving by the law . (d) At , we see further evolution of pollutant.

Finally, we simulate dynamics of sinking pollutant on a square with an uneven bottom (so, our domain in which we solve the equations is irregular). The pollutant originates from a source and the source is closed after a while (i.e., during an accident on the sea involving the chemical tanker). In this case, we can disregard behavior at the bottom of the pond which we denote by , and thus we put

implying no movement on the bottom  for the last set of simulations;

implying no velocity change on appropriate edges.

As for the source term, we have chosen

The densities are and . Simulations are presented in Figure 4.

(a)

(b)

(c)

(d)

Figure 4: Dynamics of the density distribution of a heavy pollutant leaking from a chemical tanker, with a hull breach. Here, red, green, and blue color denote densities 1, 0.9, and 0.7, respectively. (a) At the initial moment, pollutant starts to leak. (b) Leaking is continued after iterations and it expands along the surrounding fluid. (c) Leaking is stopped after iterations. (d) At , we observe how gravity pulls the pollutant downwards.

5. Conclusion

We proposed a mathematical model for describing extraction of one fluid from another as well as propagation of a fluid leaking into another fluid. By analyzing the problem of dynamics of two almost immiscible fluids in the context of pollution problems, we have concluded that we can control its propagation by moving the extraction point. We have also described how to predict propagation of a pollutant from, for example, a chemical tanker.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The authors would like to thank the referee for careful reading of the paper as well as for generous help in providing constructive comments and suggestions. They would also like to thank Professor Mladen Jurak for his help during writing of the paper. The research is supported by the Croatian Science Foundation’s funding of the project “Weak convergence methods and applications” no. 9780 and by the FP7 project “Micro-local defect functional and applications” (MiLDeFA) in the frame of the program Marie Curie FP7-PEOPLE-2011-COFUND.

References

M. Orlić and M. Lazar, “Cyclonic versus anticyclonic circulation in lakes and inland seas,” Journal of Physical Oceanography, vol. 39, no. 9, pp. 2247–2263, 2009.

S. Chandrasekhar, Hydrodynamic and Hydromagnetic Stability, Courier Corporation, 1961.

A. T. Ismail-Zadeh, H. E. Huppert, and J. R. Lister, “Gravitational and buckling instabilities of a rheologically layered structure: implications for salt diapirism,” Geophysical Journal International, vol. 148, no. 2, pp. 288–302, 2002.

A. Ismail-Zadeh, I. Tsepelev, C. Talbot, and A. Korotkii, “Three-dimensional forward and backward modelling of diapirism: numerical approach and its applicability to the evolution of salt structures in the Pricaspian basin,” Tectonophysics, vol. 387, no. 1–4, pp. 81–103, 2004.

G. I. Khodakovskii, V. P. Trubitsyn, M. Rabinowicz, and G. Cenleneer, “Magma migration and compaction of surrounding rocks with variable viscosity,” in Hydrodynamics and Earthquake Forecasting, vol. 30 of Computational Seismology, pp. 16–31, Nauka, Moscow, Russia, 1998.

D. Mckenzie, “The generation and compaction of partially molten rock,” Journal of Petrology, vol. 25, no. 3, pp. 713–765, 1984.

B. M. Naimark and A. T. Ismail-Zadeh, “Numerical models of a subsidence mechanism in intracratonic basins: application to North American basins,” Geophysical Journal International, vol. 123, no. 1, pp. 149–160, 1995.

P. Bauer and Z. Janour, “Mathematical modelling and numerical simulation of pollution transport in the atmospheric boundary layer,” in Proceedings of the Czech-Japanese Seminar in Applied Mathematics, Czech Technical University, Prague, Czech Republic, August 2004.

R. Temam, Navier-Stokes Equations, AMS Chelsea Publishing, 2000.

M. Mitrea, S. Monniaux, and M. Wright, “The Stokes operator with Neumann boundary conditions in Lipschitz domains,” Journal of Mathematical Sciences, vol. 176, no. 3, pp. 409–457, 2011.

S. N. Antontsev, A. M. Meirmanov, and V. A. Solonnikov, “Smooth interface in a two-component Stokes flow,” Annali dell'Università di Ferrara. Sezione XII: Scienze Matematiche XLVII, vol. 47, no. 1, pp. 269–284, 2001.

D. N. Arnold, F. Brezzi, and M. Fortin, “A stable finite element for the Stokes equations,” CALCOLO, vol. 21, no. 4, pp. 337–344, 1984.

V. G. Danilov and G. A. Omel'yanov, “Dynamics of the interface between two immiscible liquids with nearly equal densities under gravity,” European Journal of Applied Mathematics, vol. 13, no. 5, pp. 497–516, 2002.

D. Mitrovic, J. M. Nordbotten, and H. Kalisch, “Dynamics of the interface between immiscible liquids of different densities with low Froude number,” Nonlinear Analysis. Real World Applications, vol. 15, pp. 361–366, 2014.

B. Engquist, A.-K. Tornberg, and R. Tsai, “Discretization of Dirac delta functions in level set methods,” Journal of Computational Physics, vol. 207, no. 1, pp. 28–51, 2005.