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 formwhere 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 equationswhere, 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 CO2 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 CO2 under the water. Part of the CO2 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 CO2 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 -distributionwhere 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)
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)
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)
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 putimplying 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)
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.