Mathematical Model for Axisymmetric Taylor Flows Inside a Drop

Analytical solutions of the Stokes equations written as a differential equation for the Stokes stream function were obtained. These solutions describe three-dimensional axisymmetric flows of a viscous liquid inside a drop that has the shape of a spheroid of rotation and have a similar set of characteristics with Taylor flows inside bubbles that occur during the transfer of a two-component mixture through tubes.


Introduction
The Taylor flow regime is the most favorable of the known regimes for microchannels, since it provides, on the one hand, an intensification of mass transfer due to Taylor vortices, and on the other hand, a sufficient residence time, in view of the rather moderate velocities of phase motion (typically 0.01-0.1 m/s, and not larger than 0.8-1 m/s), which are typical for this regime [1,2]. The velocity of bubbles or droplets is characterized by capillary number Ca = µU/σ where µ is a liquid viscosity, U is adispersed phase (bubbles or droplets) velocity and σ is an interfacial tension. Usually Ca number does not exceed 0.01-0.1, and the Reynolds number for the flows in micro and mini-channels is less than 2000, corresponding thus to the laminar flow [3]. These features make it possible to achieve a high mixing intensity without significantly complicating the equipment and generate great interest in this type of flow.
In the last two to three decades, many papers were devoted to both the hydrodynamics of the Taylor flow and mass transfer in it [4][5][6][7][8][9][10]. Most of these works in fluid dynamics are based on experimental or numerical experiments. At the same time, it is important to search for independent theoretical methods that allow one to obtain a general a priori idea of the dispersed flow inside droplets or bubbles, including the shape of streamlines.
Such two-component liquid streams belong to the class of Taylor streams and are widely used in various industries [1,2].
One of the first models describing the motion of Taylor bubbles in channels was published by F. Bretherton [11]. The extended Bretherton model for Taylor bubbles at moderate capillary numbers was proposed by E. Klaseboer [12].
A number of works [13][14][15][16] were devoted to the study of Taylor flows. In [17], the method of micro Particle Image Velocimetry is used to study Taylor vortices inside an elongated bubble. In [18] the circulation patterns within the dispersed phase for regular Taylor two-phase flows were constructed using the FEATFLOW program.
It is worth noting that the actual shape of the Taylor bubbles may be far from the ideal geometric shape. Nevertheless, modeling, which considers the Taylor flow in a domain with an idealized form, is used in a number of studies [19]. In [20,21], the shape of the bubbles is approximated by a part of a cylinder with hemispherical caps or a prolate ellipsoid at ends.
In the type of flow that we are considering, a mixture of two liquids that have different viscosity coefficients is transferred through the channel and one liquid is distributed in the form of bubbles inside the other one ( Figure 1). The liquid interacting with the channel walls moves along its axis and drives the gas bubbles or the droplets of the second liquid. At a constant two-phase velocity of a mixture of liquids through the channel, this flow can be considered stationary and axisymmetric (for the case when the capillary forces dominate over gravitational forces, i.e., for liquids similar to the water at the channel diameter smaller than ≈ 3 mm).

PEER REVIEW 2 of 10
It is worth noting that the actual shape of the Taylor bubbles may be far from the ideal geometric shape. Nevertheless, modeling, which considers the Taylor flow in a domain with an idealized form, is used in a number of studies [19]. In [20,21], the shape of the bubbles is approximated by a part of a cylinder with hemispherical caps or a prolate ellipsoid at ends.
In the type of flow that we are considering, a mixture of two liquids that have different viscosity coefficients is transferred through the channel and one liquid is distributed in the form of bubbles inside the other one ( Figure 1). The liquid interacting with the channel walls moves along its axis and drives the gas bubbles or the droplets of the second liquid. At a constant two-phase velocity of a mixture of liquids through the channel, this flow can be considered stationary and axisymmetric (for the case when the capillary forces dominate over gravitational forces, i.e., for liquids similar to the water at the channel diameter smaller than  3 mm). Due to the presence of friction force at the boundary of two fluids and friction between the fluid and the channel walls, a vortex flow occurs in the inner part of the bubbles during the movement of the mixture through the channel ( Figure 2). In our previous work [22], we obtained solutions for two-dimensional flows inside an elliptical region, which demonstrated the presence of both the main toroidal vortex (similar to that shown in Figure 2) and two satellite vortices. These solutions generally reflect the nature of the droplet / bubble flows in the Taylor two-phase flow.
The system of Stokes equations is often used in mathematical models to describe the flow of a liquid with a set of characteristics similar to the Taylor flows described above. It should be noted that in most problems, the system of Stokes equations can be solved only numerically, and analytical solutions of Stokes equations have been found only for a small class of problems with a specific regular boundary shape [23,24]. In the presence of curvilinear sections of the boundary in some problems, it turns out to be expedient to use the corresponding curvilinear coordinate system. In [25], solutions were found for a circular cavity in a polar coordinate system; in [26], solutions of the Stokes equations were obtained in an elliptic coordinate system for a region bounded by two confocal ellipses. Solutions of the Stokes equations for two-dimensional flows inside an elliptic region were considered in article [22]. A similar approach is used for three-dimensional flows in cases where the boundary has a spherical or cylindrical shape [27,28].
The purpose of this work is to obtain analytical solutions of the Stokes equations for a mathematical model that describes the three-dimensional Taylor flow inside an elon-  It is worth noting that the actual shape of the Taylor bubbles may be far from the ideal geometric shape. Nevertheless, modeling, which considers the Taylor flow in a domain with an idealized form, is used in a number of studies [19]. In [20,21], the shape of the bubbles is approximated by a part of a cylinder with hemispherical caps or a prolate ellipsoid at ends.
In the type of flow that we are considering, a mixture of two liquids that have different viscosity coefficients is transferred through the channel and one liquid is distributed in the form of bubbles inside the other one ( Figure 1). The liquid interacting with the channel walls moves along its axis and drives the gas bubbles or the droplets of the second liquid. At a constant two-phase velocity of a mixture of liquids through the channel, this flow can be considered stationary and axisymmetric (for the case when the capillary forces dominate over gravitational forces, i.e., for liquids similar to the water at the channel diameter smaller than  3 mm).  In our previous work [22], we obtained solutions for two-dimensional flows inside an elliptical region, which demonstrated the presence of both the main toroidal vortex (similar to that shown in Figure 2) and two satellite vortices. These solutions generally reflect the nature of the droplet / bubble flows in the Taylor two-phase flow.
The system of Stokes equations is often used in mathematical models to describe the flow of a liquid with a set of characteristics similar to the Taylor flows described above. It should be noted that in most problems, the system of Stokes equations can be solved only numerically, and analytical solutions of Stokes equations have been found only for a small class of problems with a specific regular boundary shape [23,24]. In the presence of curvilinear sections of the boundary in some problems, it turns out to be expedient to use the corresponding curvilinear coordinate system. In [25], solutions were found for a circular cavity in a polar coordinate system; in [26], solutions of the Stokes equations were obtained in an elliptic coordinate system for a region bounded by two confocal ellipses. Solutions of the Stokes equations for two-dimensional flows inside an elliptic region were considered in article [22]. A similar approach is used for three-dimensional flows in cases where the boundary has a spherical or cylindrical shape [27,28].
The purpose of this work is to obtain analytical solutions of the Stokes equations for a mathematical model that describes the three-dimensional Taylor flow inside an elon- In our previous work [22], we obtained solutions for two-dimensional flows inside an elliptical region, which demonstrated the presence of both the main toroidal vortex (similar to that shown in Figure 2) and two satellite vortices. These solutions generally reflect the nature of the droplet/bubble flows in the Taylor two-phase flow.
The system of Stokes equations is often used in mathematical models to describe the flow of a liquid with a set of characteristics similar to the Taylor flows described above. It should be noted that in most problems, the system of Stokes equations can be solved only numerically, and analytical solutions of Stokes equations have been found only for a small class of problems with a specific regular boundary shape [23,24]. In the presence of curvilinear sections of the boundary in some problems, it turns out to be expedient to use the corresponding curvilinear coordinate system. In [25], solutions were found for a circular cavity in a polar coordinate system; in [26], solutions of the Stokes equations were obtained in an elliptic coordinate system for a region bounded by two confocal ellipses. Solutions of the Stokes equations for two-dimensional flows inside an elliptic region were considered in article [22]. A similar approach is used for three-dimensional flows in cases where the boundary has a spherical or cylindrical shape [27,28].
The purpose of this work is to obtain analytical solutions of the Stokes equations for a mathematical model that describes the three-dimensional Taylor flow inside an elongated bubble ( Figure 2). In the considered simplification of the physical model, the bubble shape is specified as the surface of a spheroid (ellipsoid of revolution), and the fluid flow inside the bubble is assumed to be axisymmetric. Thus, in comparison with work [22], in this Fluids 2021, 6, 7 3 of 10 article we simulate the flow for three-dimensional shape of the bubble boundary. Another form of boundary geometry leads to the emergence of a more complex axisymmetric pressure field inside the bubble and becomes a more accurate approximation of the real physical model.
The obtained analytical solutions of the Stokes equations make it possible to qualitatively compare the flow pattern inside the bubbles with experimental data and the results of numerical modeling for Taylor flows, and can also be used as a benchmark to assess the quality of numerical algorithms [29,30].

Problem Formulation
We are considering a three-dimensional axisymmetric Stokes flow in the region bounded by a spheroid (Figure 3). We assume that the boundary of the spheroid is impenetrable and the shape of the spheroid does not change over time. The fluid flow inside the spheroid is caused by the motion of the boundary, in which the velocity of the points on the boundary is directed tangentially to the boundary of the spheroid. In this problem, we will assume that the velocity vector at the boundary of the spheroid lies in the plane passing through the axis of symmetry of the spheroid.
R REVIEW 3 of 10 gated bubble (Figure 2). In the considered simplification of the physical model, the bubble shape is specified as the surface of a spheroid (ellipsoid of revolution), and the fluid flow inside the bubble is assumed to be axisymmetric. Thus, in comparison with work [22], in this article we simulate the flow for three-dimensional shape of the bubble boundary. Another form of boundary geometry leads to the emergence of a more complex axisymmetric pressure field inside the bubble and becomes a more accurate approximation of the real physical model. The obtained analytical solutions of the Stokes equations make it possible to qualitatively compare the flow pattern inside the bubbles with experimental data and the results of numerical modeling for Taylor flows, and can also be used as a benchmark to assess the quality of numerical algorithms [29,30].

Problem formulation
We are considering a three-dimensional axisymmetric Stokes flow in the region bounded by a spheroid (Figure 3). We assume that the boundary of the spheroid is impenetrable and the shape of the spheroid does not change over time. The fluid flow inside the spheroid is caused by the motion of the boundary, in which the velocity of the points on the boundary is directed tangentially to the boundary of the spheroid. In this problem, we will assume that the velocity vector at the boundary of the spheroid lies in the plane passing through the axis of symmetry of the spheroid. The condition specifying the tangential velocity at the boundary will be written as: The radial component of the velocity on the axis of symmetry of the spheroid equals to zero: The described type of flow inside the spheroid has a similar set of characteristics with Taylor flows inside bubbles ( Figure 2) and represents their simplified mathematical model.
Let us consider the derivation of solutions to the Stokes equations that describe fluid flows for a given system. We introduce a cylindrical coordinate system (r, z, ϕ) centered at a point that halves the segment lying between the foci of the spheroid (Figure 3). Since the flow is axisymmetric, we can write equations in two coordinates (r, z). A cross section of a spheroid is shown in Figure 4. The axis OZ divides the border of the ellipse into two parts, S 1 and S 2 . We assume that the distribution that defines the velocity at the boundary of the spheroid has the same form on any half and an ellipse (S 1 , S 2 ) in any plane passing through the axis OZ. In elliptical coordinates, the halves of the ellipse ( Figure 4) will transform into two rectangles, given by the inequalities:   Let us formulate the boundary conditions for the Stokes stream function in elliptic coordinates. Condition (4) will be written as:  The condition specifying the tangential velocity at the boundary will be written as: Fluids 2021, 6, 7 The radial component of the velocity on the axis of symmetry of the spheroid equals to zero: v r | OZ = 0 In the process of describing the fluid flow and constructing the solution of the Stokes equations in this work, we will use the Stokes stream function (it is necessary to distinguish this function from the stream function used in hydrodynamics to describe two-dimensional flows). Velocity components are related to the Stokes stream function by the following relations [23]: Let us formulate the boundary conditions for the Stokes stream function. Equation (2) will be written in the form: The normal component of the velocity at the ellipse boundary is equals to zero, therefore: The condition for the tangential velocity (1) at the boundary of the ellipse takes the form: Taking into account the specific shape of the domain boundary, the solution of Stokes equations in elliptic coordinates (ξ 1 , ξ 2 , ϕ) was derived: The Lame coefficients in elliptic coordinates are: . The Stokes equation for the Stokes stream function are written as [31]: where the operator E 2 has the form: Taking into account relation between coordinates r = dsinh(ξ 1 ) sin(ξ 2 ), we get: It should be noted that Equation (6) differs from the classical biharmonic equation for the two-dimensional stream function [26]: Formula (8) includes the last two terms containing the hyperbolic cotangent and the cotangent, which are absent in the formula for the Laplace operator in elliptic coordinates.
In elliptical coordinates, the halves of the ellipse (Figure 4) will transform into two rectangles, given by the inequalities: 0 ≤ ξ 1 ≤ ξ 0 , 0 ≤ ξ 2 ≤ π and 0 ≤ ξ 1 ≤ ξ 0 , π ≤ ξ 2 ≤ 2π ( Figure 5). The segment between the foci will correspond to the segments ξ 1 = 0, 0 ≤ ξ 2 ≤ π, π ≤ ξ 2 ≤ 2π. Parts of boundary S 1 , S 2 will transform to segments ξ 1 = ξ 0 , 0 ≤ ξ 2 ≤ π, π ≤ ξ 2 ≤ 2π. The segments between the foci and points A, B will be segments ξ 2 = π, ξ 2 = 0, 0 ≤ ξ 1 ≤ ξ 0 .  s formulate the boundary conditions for the Stokes stream function in elliptic tes. Condition (4) will be written as:  (6), we perform the following ent: replacement makes it possible to significantly simplify formula (8) for the difoperator: ill seek a solution in the form: Let us formulate the boundary conditions for the Stokes stream function in elliptic coordinates. Condition (4) will be written as: Condition (3) takes the form: We write the boundary condition for the tangential velocity (5) on the right half of the ellipse in the form: Further, we assume that the velocity at the boundary of the ellipse is symmetric about the axis OZ.

Problem Solution
In the process of constructing a solution of Equation (6), we perform the following replacement: This replacement makes it possible to significantly simplify Formula (8) for the differential operator: We will seek a solution in the form: Integration of corresponding differential equations, gives us: In elliptical coordinates, the stream function takes the form: Fluids 2021, 6, 7 6 of 10 A detailed derivation of Formula (15) is given in the Appendix A. It is worth noting that a similar technique for constructing solutions was previously used to describe the flow outside a solid body in the form of a spheroid during its motion in a viscous fluid [31].
To complete the process of constructing the Stokes stream function defining the fluid flow inside the bubble we need to find the coefficients C 1 , C 2 , C 3 in the expression (16) that satisfy the boundary conditions.
Due to the peculiarity of constructing the stream function in the form (16), boundary condition (10) will be satisfied automatically.
From the condition on the segment between the foci of the ellipse (11), we obtain linear equations for the coefficients: From condition (9) we get: where λ 0 = cosh(ξ 0 ), the value ξ 1 = ξ 0 defines the boundary of the spheroid.
Equations (17) and (18) define two linear equations for three unknown constants. Thus, one of the constants C 1 , C 2 , C 3 remains free, which makes it possible to vary the velocity of the border.
It should be noted that due to axial symmetry, the problem of constructing a solution inside an ellipse is reduced to the problem of constructing a solution in its half. Due to the peculiarity of the elliptical coordinate system, the Stokes stream functions given in different halves of the ellipse (Figure 4) will be symmetric about the axis OZ (functions for different parts will have the same modulus and the opposite sign). This will ensure the fulfillment of the boundary conditions and the continuity condition of the stream function and its derivatives on the axis OZ.

Results and Discussion
The obtained analytical solution (16) opens up an opportunity for modeling of fluid flow inside the spheroid. It is worth noting that functions of the form (16) impose specific restrictions on the velocity function at the boundary. Due to the presence of the sine of the coordinate ξ 2 in Formula (16), the flow velocity function tends to zero at the ends of the ellipse and increases at its middle part. However, this feature is natural for the Taylor flows.
Let us consider an example of such a flow, given by an analytical stream function for a certain choice of constants in Formula (16). Velocity function at the boundary S 1 is shown in Figure 6. Flow streamlines are shown in Figure 7. Streamlinesin this example correspond to the following coefficients in (16):  A flow with a similar set of parameters can be observed in the by-pass flow mode in microchannels, as well as inside a bubble moving far from the channel walls. In this type of flow, streamlines are around the bubble and are not closed, as shown in Figure 2. This feature of the flow leads to the appearance of shear stresses at the interface between liquids, acting from the nose of the bubble to its tail [32,33], which create circulation inside   A flow with a similar set of parameters can be observed in the by-pass flow mode in microchannels, as well as inside a bubble moving far from the channel walls. In this type of flow, streamlines are around the bubble and are not closed, as shown in Figure 2. This feature of the flow leads to the appearance of shear stresses at the interface between liquids, acting from the nose of the bubble to its tail [32,33], which create circulation inside the bubble (Figure 7).
The analytical solution (16) can be used to simulate a three-dimensional axisymmetric flow between two confocal spheroids. In this case, the normal velocity at the boundary of the inner spheroid should be equal to zero. Condition (17) must be replaced by a condition that ensures the equality of the stream function to zero at the boundary of the inner spheroid: where 1 cosh( )     ,the value 1    defines the boundary of the inner spheroid. It should be noted that a similar problem of the Stokes flow between two cylinders with an elliptical profile was previously considered in [26].   A flow with a similar set of parameters can be observed in the by-pass flow mode in microchannels, as well as inside a bubble moving far from the channel walls. In this type of flow, streamlines are around the bubble and are not closed, as shown in Figure 2. This feature of the flow leads to the appearance of shear stresses at the interface between liquids, acting from the nose of the bubble to its tail [32,33], which create circulation inside the bubble (Figure 7).
The analytical solution (16) can be used to simulate a three-dimensional axisymmetric flow between two confocal spheroids. In this case, the normal velocity at the boundary of the inner spheroid should be equal to zero. Condition (17) must be replaced by a condition that ensures the equality of the stream function to zero at the boundary of the inner spheroid: where λ 1 = cosh(ξ * ), the value ξ 1 = ξ * defines the boundary of the inner spheroid. Figure 8 shows the streamlines for the flow between two spheroids. The velocity function at the boundary of the spheroids is shown in Figure 9. The upper function (highlighted in blue) corresponds to the velocity on the outer spheroid, the bottom function (highlighted in red) corresponds to the velocity on the inner spheroid. Figures from this example correspond to the following coefficients in (16):  A flow with a similar set of parameters can be observed in the by-pass flow mode in microchannels, as well as inside a bubble moving far from the channel walls. In this type of flow, streamlines are around the bubble and are not closed, as shown in Figure 2. This feature of the flow leads to the appearance of shear stresses at the interface between liquids, acting from the nose of the bubble to its tail [32,33], which create circulation inside the bubble (Figure 7).
The analytical solution (16) can be used to simulate a three-dimensional axisymmetric flow between two confocal spheroids. In this case, the normal velocity at the boundary of the inner spheroid should be equal to zero. Condition (17) must be replaced by a condition that ensures the equality of the stream function to zero at the boundary of the inner spheroid:  It should be noted that a similar problem of the Stokes flow between two cylinders with an elliptical profile was previously considered in [26].  . Velocity function at S1 for outer and inner spheroids. It should be noted that a similar problem of the Stokes flow between two cylinders with an elliptical profile was previously considered in [26]. he boundary of the spheroids is shown in Figure 9. The upper function in blue) corresponds to the velocity on the outer spheroid, the bottom functed in red) corresponds to the velocity on the inner spheroid. Figures from correspond to the following coefficients in (16): be noted that a similar problem of the Stokes flow between two cylinders ical profile was previously considered in [26]. ure of the streamlines in the een outer and inner coaxial . Results of calculations. Figure 9. Velocity function at S1 for outer and inner spheroids.

Conclusions
The geometric shape of the Taylor bubble can naturally be approximated by a spheroid of revolution. Therefore, a comparison of Taylor flows with axisymmetric Stokes flows inside a spheroid is quite applicable, despite the fact that the actual shape of the bubble can be more complicated due to the presence of microwaves on its surface.
This work is devoted to the derivation of exact analytical solutions of the Stokes equations for the model of axisymmetric flows inside the Taylor bubble. These solutions can be used as a tool for assessment of the numerical algorithms quality by comparing the result of the numerical algorithm with exact analytical solution. The advantage of this evaluation method is that it excludes the possibility of error accumulation which can occur in the case of comparing numerical algorithms with each other.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Derivation of a Stokes Stream Function
Let us consider the derivation of the solution for Equation (6). Substitution of expression (14) into Equation (13), gives us: Let us introduce the function: For function E 4 Ψ we get the following expression: It is worth noting that the term in square brackets must equals to zero to execute the Equation (6). Since G is a function of a variable λ, this implies G (λ) = 0 G(λ) − λG (λ) = 0 (A4) The solution of this system is a function: Substituting this function in (A2), we obtain the inhomogeneous equation: Particular solution of this equation is a function of the following form: Another part of the solution will correspond to the solution of the following homogeneous equation: (λ 2 − 1)g (λ) − 2g(λ) = 0 (A8) Integration gives us: (λ 2 − 1)g (λ) − 2λg(λ) = C 2 (A9) By integrating this equation, we get: The general solution to Equation (A6) takes the form: For a stream function satisfying Equation (6), we obtain the following expression: