Numerical Investigation of Two-Phase Flows in Corrugated Channel with Single and Multiples Drops

: This work aims at investigating numerically the effects of channel corrugation in two-phase ﬂows with single and multiples drops subject to buoyancy-driven motion. A state-of-the-art model is employed to accurately compute the dynamics of the drop’s interface deformation using a modern moving frame/moving mesh technique within the arbitrary Lagrangian–Eulerian framework, which allows one to simulate very large domains. The results reveal a complex and interesting dynamics when more than one drop is present in the system, leading eventually in coalescence due to the amplitude of the corrugated sinusoidal channel and distance between drops.


Introduction
Two-phase flows with variable cross section still remain a challenging subject with great interest from the scientific community and industry related to the efficient cooling of computer parts and biomedical devices. Moreover, channels with different corrugation patterns seem to considerably improve the effect of heat dissipation by modifying the drop's dynamics and thus the liquid film thickness of the two-phase system. Corrugated channel is a class of variable cross section channels that present a periodic pattern along its length and has strong impact in the fluid flow, changing its behavior by increasing recirculation zones and eventually leading to onset of boundary layer detachment. Note that the presence of a deformable interface due to the presence of more fluid phases brings an extra degree of complexity in the fluid flow system. Additionally, the periodic pattern of the corrugated channels requires very large domains to be analyzed, thus dramatically increasing the costs associated to the experiment assembling or the numerical simulations.
Single-phase flows with heat transfer were experimentally analyzed in [1,2], where the effects of several key parameters were studied such as channel wavy amplitude and flow velocity in periodic wavy passages. Numerically, several authors have investigated single-phase flows in corrugated channels. For instance, in [3] a numerical study was done for transition regimes from laminar to turbulent when convective heat transfer plays an important role in those flows. They concluded that the Nusselt number for the wavy wall channel was directly affected when compared to those for a straight channels, while in [4], a 2-dimensional laminar steady and time-dependent fluid flow solver was developed for heat transfer in periodic wavy channels with sinusoidal and arc shapes. It was concluded that the arc-shaped periodic channel enhances heat transfer due to a higher friction factor when compared to the sinusoidal channel, specially beyond the critical value of the Reynolds number for the unsteady regime. Another interesting numerical approach was suggested by the authors of [5], where the coordinate transformation method and the spline alternating-direction implicit method was used to study the rates of heat transfer for flow through a sinusoidally curved converging-diverging channel. The steady streamfunction vorticity formulation was used to model the fluid flow with the energy equation using the finite difference method. Several channel amplitude wavelengths were tested at different Reynolds number, and it was found that at a sufficiently larger value of amplitude wavelength ratio the corrugated channel increases heat transfer for large Reynolds numbers. Turbulent forced convection in a wavy channel using a two-phase model was studied in [6] for water-Al 2 O 3 nanofluids using a 2-dimensional numerical solver. It has been seen that if the channel amplitude increases, as well as the Reynolds number and volume fraction of nanoparticles, the Nusselt number is directly affected, leading to a higher value. Additionally, it was observed that the more nanoparticles added to the solution, the higher pressure drop in corrugated channels is seen.
It is already known that two-phase flows bring another level of challenges numerically and experimentally due to the presence of an interface between fluids and different phases in the system. In [7], a comparative analysis of two-phase flow in sinusoidally corrugated channels was carried out with application to polymer electrolyte membrane fuel cells. Different patterns of corrugated channels were used to numerically analyze single and multiples bubbles dynamics in air-water systems through a 3-dimensional simulation using the Volume-of-Fluid (VOF) method. The authors concluded that the two-phase flow performance in the sinusoidal channels can be substantially improved, decreasing the radius of curvature at the bends and making the confining walls extremely hydrophobic. On the other hand, in [8], experimental and numerical approaches have been used to tackle CO 2 bubble behavior flowing in a methanol solution in sinusoidally corrugated channels for fuel cell purposes. A detailed investigation of key two-phase flow parameters was undertaken, and a thermal analysis was made regarding when a single bubble travels along sinusoidal channels under pressure differences between the inlet and outlet regions of the channel. Additionally, a comparison of the same flow condition was carried out for a straight channel. The results reveals that the effect of flow disturbance structure is negligible with small values of channel amplitude and angular frequency of corrugation for the studied two-phase system. Sinusoidal channels with one and more wavy patterns were numerically simulated using the front-tracking and the finite volume methods in [9]. Two-phase 2-dimensional flow simulations were carried out to investigate single buoyant drops in channel constrictions and expansions and the interaction between two small drops in different simulation conditions with the newly implemented numerical code. Results have been successfully reported for various constricted channels.
The current study is based upon the experiments in [10], where the authors have studied experimentally the effects of buoyancy-driven motion of drops in a periodically constricted capillary. Several two-phase systems were used in the experiment such as glycerol-water, diethylene glycol, and diethylene glycol-glycerol as suspending fluid, and silicon oil, Dow Corning fluids, and UCON fluid emollient as drop fluids, and an extensive study has been reported for sinusoidal shape channels with fixed dimensionless amplitude A = 0.07 and wavelength λ = 4 for single drops. However, no further investigation has been conducted for different channel amplitude and its effect on interface break-up for single and multiple drops. Additionally, it is clear from the literature review presented above that much is still required for the complete understanding of two-phase flows in short and large periodic arrays with different corrugation levels. It is indeed of great significance that single and multiple buoyant drops in sinusoidal channels deliver another level of flow complexity due to the possible significant channel amplitude perturbation in the two-phase system, leading eventually to drop coalescence and/or interface breakup. The aim of this work is to investigate numerically two-phase flow key parameters such as the drop's rising velocity, film thickness, and interface perimeter change for one and three drops in buoyancy-driven motion for a diethylene glycol/UCON-1145 (DEG3) two-phase system in sinusoidal channels with dimensionless wavelength λ = 4 and wave amplitude A = 0.14. Additionally, a parametric study is also carried out to investigate the influence of channel wavelength and drop volume in coalescence during slug flows. A state-of-the-art model is employed to accurately compute the dynamics of the drop's interface motion using a modern moving frame/moving mesh technique within the arbitrary Lagrangian-Eulerian framework. The presented results show the drop's dynamics in large periodic channels and coalescence even for large liquid slug distances.
This article is organized with a literature review of single-and two-phase corrugated channels, followed by the mathematical description of the modeling equations in axisymmetric formulation based on the Finite Element Method (FEM). The next section provides details concerning the numerical modeling used to discretize the mathematical equations, followed by the result section where several important test cases are presented to assess and analyze the dynamics of two-phase flow in corrugated channels with single and multiples bubbles. Finally, this text ends with remarks in the conclusion section.

Two-Phase Flow Equations in the ALE Context
Momentum and continuity equations are used as the model for incompressible motion of two-phase flows in the context of a "one-fluid" description where one set of equations is used domain-wise and a jump in fluid properties is assumed at the interface. The dimensionless form of the momentum and continuity equations respectively is shown in cylindrical coordinates for axisymmetric flows with gravity g and surface tension force f accounted in as where x is the position vector with radial r and axial x coordinates. v is the velocity vector with components v r and v x ,v is the mesh velocity vector with componentsv r andv x . The velocity difference v −v is due to the arbitrary Lagrangian-Eulerian framework; otherwise, only v is found when a fixed frame (Eulerian) is established. ρ(x) is the fluid density, µ(x) is the fluid dynamic viscosity spatially distributed with a jump condition at the interface Γ, f is the surface tension with radial f st r and axial f st x components, respectively, p is pressure and t is time. Note here that the fluid properties µ(x) and ρ(x) are dimensionless and they appear in the two-phase flow equation since both are usually different among phases; however in single-phase flows both are constant and equal to the unity, thus vanishing from the momentum equation. As flow is buoyancy-drive driven, the dimensionless groups Archimedes (N) and Eötvös (Eo) appear in the momentum equation according to the standard dimensional analysis and are both defined as follows, where the subscript 0 stands for the reference values of dynamic viscosity µ, surface tension σ, gravity g = 9.81m/s 2 , and channel diameter D. ∆ρ is defined as the difference of densities among phases.

Axisymmetric Differential Operators
The dimensionless axisymmetric strain tensor in the momentum equation does not account for the mesh motion itself; however, in two-phase flows, both the velocity gradient (∇v) and the gradient transpose (∇ T v) shall be considered due to fluid properties transition at interface, resulting in nonuniform spatially variation in density ρ(x) and µ(x). Therefore, it is written according to the following expression, The divergent (∇·) and transpose gradient (∇ T ) defined as follows,

General Description of Two-Phase Flow Solver
The current numerical methodology is based on the works in [11][12][13], and a short description is given for the sake of completeness, with detailed explanation of the triangle finite element used in this work, as well as the geometrical parameters used to assemble the numerical test sections and the finite element mesh treatment. The current numerical code has been implemented in C/C++ language using object-oriented programming with aid of PETSc framework ( [14]) for the linear system solvers and preconditioners, the finite element mesh generator GMsh [15] for initial interface and boundary mesh setup, and the triangle mesh generator package Triangle [16] for adaptive remeshing support.
The classical Galerkin method has been used as approximating method to the variational formulation of the governing equations resulting in the discretization of all terms except the nonlinear advection term (v −v) · ∇v. Such a term is known to be the source of undesired spurious modes that may cause numerical instabilities if the Galerkin method is used. Thus, to overcome such an undesired issue, the semi-Lagrangian method [17,18] has been successfully employed in this work. The idea of representing the acceleration field by a Lagrangian framework instead of the commonly found Eulerian framework. Note that such a method is applied to the ALE convective velocity difference v −v within the ALE framework. For each time step the nodes move towards the flow characteristics and a re-initialization of the system coordinates is performed, thus recovering the initial mesh. The substantial derivative with the ALE velocity difference is evaluated in the strong form along the characteristic trajectory, by estimating the position of a point and solving the vectorial equation Dx/Dt backwards in time with the initial condition x(t n+1 ) = x i , where x i is the mesh node coordinate, than an integration method is used to evaluate the previous node position in the triangle unstructured mesh. A first order discretization scheme is adopted assuming the trajectory is linear.
The numerical domain is discretized using an unstructured triangle mesh with a LLBstable mini-element with details given in the next section. Once the discretization of the domain is accomplished, the system matrices are assembled and the solution of the time dependent 2-dimensional axisymmetric equations is then found by successively solving the linear system using PETSc framework in each time step for pressure and velocity. However, due to the strong coupling between pressure and velocity, the numerical procedure implemented to solve the mentioned linear system uses the Projection method based on the LU decomposition, which was first introduced by [19]. The aim of this method is to uncouple pressure and velocity and solve each quantity separately, thus reducing the large linear system size into smaller ones. The solution of the linear system for pressure is given by a direct solver based on LU decomposition, and the velocity is solved by the conjugated gradient method using the incomplete Cholesky factorization as preconditioner as the associated linear system is positive symmetric definite.

Choice of the Finite Element
The linear element with cubic bubble function uses three shape functions that are interpolated at the vertex of the element supplemented by a cubic bubble function, which is interpolated at the element centroid. As only four degrees of freedom are used for a cubic polynomial, the shape functions of this element are incomplete. This element is also know as mini element. Such a scheme gives an LBB stable element without having to introduce too many degrees of freedom. The shape functions, in terms of barycentric coordinates, are where N i are the shape functions of the corner nodes of the mini-triangle element Figure 1b, and N 4 is the shape functions of the triangle centroid node. The associated boundary element is a linear line element and its shape functions are simply the functions of the barycentric coordinates λ b 1 and λ b 2 . As can be seen in Figure 1, the combination of both triangles generates the well-known mini-element, which is LBB stable and does not require any additional term in the fluid motion equations. The linear element was used to evaluate pressure at the corner nodes and mesh update procedure, while the mini and linear boundary elements were used to evaluate velocity at all element nodes. Interpolation nodes for the triangle finite elements and boundary finite elements used in the current work. The linear element was used to evaluate pressure at the corner nodes and mesh update, while the mini and linear boundary element were used to evaluate velocity at all element nodes. The combination of these two triangle elements generates the P1 + P1 finite element.

Test Sections
Two numerical test sections were used in the present work, namely, straight and sinusoidal channels. Each test section was tested with single-drop flow and consecutive three-drop flow with different liquid slug thickness s, as can be seen in Figure 2; thus, resulting in a simplified slug flow pattern where the center drop's flow is directly affected by the front and rear drop's motion. If the amplitude A of the sinusoidal channel is set to A = 0.00, one achieves the straight channel, and the same numerical method and mesh can be used; therefore, allowing the investigation of the influence of channel wall into the drop's dynamics. The sinusoidal wall profile is set through the harmonic wave equation: where D is the channel's diameter and D/2 is its radius; A is the wave amplitude that was set to A = {0.00, 0.14} for the straight and sinusoidal channels, respectively; λ is the wavelength; and φ is its phase. The axial coordinate is x, while x re f refers to the drop's referential point. The moving boundary approach has been discussed in details in [12], where the drop or drops center of mass remain fixed in space while the boundary moves according to its position increment with time. Such a procedure has been successfully proven to be efficient and allowing large periodic domain to be simulate with low computational costs. Additionally, four dimensionless liquid slug lengths, s, have been used: s = {0.375, 0.875, 1.375, 1.875}, thus varying significantly the axial distance between drops. The drop's size factor κ d = 1.72 relative to the channel's radius (D/2) was used to set the drop's dimension for all simulations presented in this work. The chosen drop's size factor implies a drop with a radius of 1.72 times the channel's radius by considering the equivalent radius of a spherical drop with same volume. It is important to note that the initial drop's shape presents a approximated drop's shape found in two-phase flows in channels (see, e.g., in [20][21][22]).

Mesh Treatment
The numerical code solves the fluid motion equations using a moving mesh technique within the ALE framework, which requires an efficient adaptive mesh refinement when large mesh deformation occurs. It is important to note that the current code implementation moves the mesh nodes and decides whether deletion or insertion of nodes and elements is required; however, the connective matrix needed by the FEM solver is generated using the package Triangle. The spatial mesh node distribution is given through the solution of the generalized Helmholtz equation for the whole domain, and the decision whether to insert or delete nodes is performed based upon a given predefined tolerance. The solution is numerically obtained using the FE method using the classical discrete operators with the triangle linear element. Figure 3 presents a detailed view of the unstructured triangle finite element mesh using during the simulation of three drops flow in axisymmetric sinusoidal channel, where the bottom boundary stands for the axis of symmetry. It illustrates in details the resulting mesh during the corrugated channel passage of two sequential drops with large interface deformation. Furthermore, recall that an extra centroid node is used in the calculation of the velocity field that is not appearing in the current figure. As can be seen, the minimum liquid film thickness in the front drop is covered by a minimum of 5 triangle elements and 9 nodes; thus, capturing accurately the dynamics within such an important region in two-phase flows. Not less important, the unstructured triangle mesh is updated every time step due to the Lagrangian motion of the interface nodes as well as the phase nodes, both using the mesh velocityv. It is indeed noticeable that insertion and deletion of nodes as well as flipping of triangle edges are required to complete the mesh treatment and maintain consistence of the numerical algorithm thus resulting in an accurate method to track interface motion under high deformation.

Results
We present the results of single and multiples drops rising in channels with two levels of corrugation amplitude (A = {0.00, 0.14}) obtained with the ALE-FE Two-Phase flow solver for axisymmetric coordinate system. Note that such a code has been extensively validated against several single and two-phase flows problems and have been reported in different well-recognized international journals. All numerical results presented here refers to fluid the two-phase system diethylene glycol/UCON-1145 from [10] also known as DEG3 system. According to the experiments, the drop fluid (UCON-1145) presented the dynamic viscosity of µ in = 530 mPa and density of ρ in = 995 kg/m 3 . The suspending fluid (diethylene-glycol) presented dynamic viscosity and density of µ out = 28 mPa and ρ out = 1110 kg/m 3 , respectively. Gravity was set to standard condition: g = 9.81 m/s. The surface tension coefficient was set to σ = 0.332 mN/m, which is approximately 10 times lower than the original value found in the aforementioned reference (σ r = 3.2 mN/m), thus meaning the dimensionless parameters used at all simulations were the Archimedes number N = 15,417 and the Eötvös number Eo = 340.284. These dimensionless numbers were chosen to allow large drop deformation due to low surface tension force while keeping numerical stability with large time steps. The initial mesh to all test cases had approximately 16,500 nodes including centroid and 10,060 triangle elements and finished with 43,000 nodes including centroid and 28,200 triangle elements due to the adaptive remeshing described in the previous section.
A total of 10 test cases using the geometries presented in Figure 2 for the straight channel (A = 0.00) and sinusoidal channel (A = 0.14) with single and three consecutive drops were made, namely, single drop in straight channel, single drop in sinusoidal channel, slug flows in straight channel with initial slug length, and slug flows in sinusoidal channel, both with initial slug length ranging from s = {0.375, 0.875, 1.375, 1.875}; thus, the same slug lengths s were used as comparison between different channel shapes. Additionally, 14 test cases were carried out to investigate drop's coalescence time t c in slug flows, varying systematically the drop's size factor κ d = {0.86, 1.29, 1.72, 2.15} and keeping constant the wavelength λ = 4 and the initial slug length s = 1.375, and varying systematically the channel's wavelength λ = {1, 2, 4, 6, 8} and keeping constant the drop's size factor κ d = 1.72 and the initial slug length s.
The unsteady behavior of all flows in sinusoidal channel can be easily noted at all plots presented below due to their wavy response with respect to the evolution with time of drop rising velocity and minimum film thickness, which is the minimum distance between any interface node to the closest wall node. On the other hand, the straight channel flow presents negligible shape variation after a short initial transient regime and one can be noted by checking out the full straight lines in the same images colored in red and blue. Furthermore, recall that the initial bubble shape at all simulations were the equal with drop's size factor κ d = R d /(D/2) = 1.72 and a short transient period for all test cases were expected with time 0.0 ≤ t ≤ 28.0 for the straight channel with single and multiples drops and 0.0 ≤ t ≤ 21.0 for all sinusoidal channels simulations. Figure 4 shows the steady shape of a single drop in straight channel with blue color representing the suspending fluid (diethylene-glycol) and the drop fluid (UCON-1145). As can be seen, the large minimum film thickness at the tail of the drop stabilizes the flow and no further variation in its shape is seen. Additionally, due to the size of the drop relative to the channel's diameter, the channel wall constricts the drop's shape to an elongated drop with a round nose due to the hydrodynamics and compacts the tail from which the minimum liquid film thickness δ is measured. The drop's shape, drop's velocity, and the minimum liquid film thickness remain constant for t > 28.0. A more challenging test case is presented in Figure 5 where a single drop shape evolution with time in sinusoidal channel can be seen during its passage in one corrugated length where the drop undergoes severe topological change. The amplitude of the wavy channel wall is A = 0.14, and a detailed description of the geometry is presented in Figure 2. At time t = 1375.775, the drop leaves one divergent section and approaches the convergent section; however, the channel's constriction, where the minimum crosssectional area is found, stretches out the drop's shape to its minimum thickness where the interface approaches considerably the symmetry wall; therefore, suggesting a possible topological change due to the interface breakup. However, for the two-phase system numerically investigated here, no breakup has been seen for any of the presented simulations with viscosity ratio µ in /µ out = 18.92. According to the work in [10], the breakup of drops was observed in systems with lower viscosity ratios, more specifically in GW5 systems with (µ in /µ out = 0.22). Later, with time t = 1390.774, the drops shrinks to its minimum shape length and reaches the minimum film thickness to the wall before the convergent section begins. In time t = 1418.696 the drop stretches over again and a fixed periodic shape motion is noted with time t → ∞.  Slug flow is now analyzed, and the drop's shape evolution with time in straight channel with amplitude A = 0.00 and initial liquid slug length s = 1.875 is presented in Figure 6. As can be noted, the center drop (green color) approaches the front drop (red color) at time t = 908.634, where coalescence takes place. It can also be seen that before drops collision, the center drop elongates and the minimum liquid film thickness δ becomes larger if compared to time t = 584.694. Such a behavior is due to the increase of drop's velocity as result of lower pressure behind the front drop. The same physical process was also seen when the initial slug length s is smaller, i.e., s = {0.375, 0.875, 1.375}. However, the due to the smaller initial distance of the initial slug, the coalescence time is shorter; therefore, for s = {0.375, 0.875, 1.375}, the following coalescence times were computed: t c = {105.857, 228.373, 498.582}, respectively. The drop shapes of the rear and front drops remained fairly constant during the simulations.  A slug flow is analyzed in a corrugated channel with amplitude A = 0.14, and the drop's shape evolution with time is presented in Figure 7 for the initial liquid slug length s = 1.875. The center drop (green color) approaches the front drop (red color) at time t = 758.822, where coalescence takes place, which is faster than the straight channel with same initial conditions. The same physical process was also seen when the initial slug length s was smaller, i.e., s = {0.375, 0.875, 1.375}. As seen in the straight channel, due to the smaller initial distance of the initial slug, the coalescence time is shorter, therefore for s = {0.375, 0.875, 1.375} the following coalescence time were computed t c = {17.692, 113.149, 404.622} respectively. The drop shapes of the rear and front drops presented fairly the same dynamical shape of the single drop in sinusoidal channel of Figure 5.  The evolution of the drop's rising velocity with time is presented for all test cases simulated in this work in Figure 8. The subplots are divided according to the initial slug length s, and the single drop flow in straight and sinusoidal channels is used as a reference at all plots. As can be seen, the smaller is the initial slug length s, the faster the coalescence takes place, thus resulting in topological change of the two front drops. The drops flow in sinusoidal presented same wavy pattern of velocity variation with time to all test cases. The full lines were used to represent the drop velocity of single drop in straight channel (red color) and sinusoidal (blue color) and the others line types were used for the slug flows in both channels. A=0.00, single A=0.00, rear A=0.00, center A=0.00, front A=0.14, single A=0.14, rear A=0.14, center A=0.14, front A=0.00, single A=0.00, rear A=0.00, center A=0.00, front A=0.14, single A=0.14, rear A=0.14, center A=0.14, front A=0.00, single A=0.00, rear A=0.00, center A=0.00, front A=0.14, single A=0.14, rear A=0.14, center A=0.14, front A=0.00, single A=0.00, rear A=0.00, center A=0.00, front A=0.14, single A=0.14, rear A=0.14, center A=0.14, front Two more interesting results are presented in Figure 9 for initial slug length s = 1.875, where the variation of the ratio of the drop's perimeter P relative to its initial perimeter P init is investigated in two corrugated lengths in Figure 9a, and the minimum film thickness δ evolution is presented for two corrugated lengths in Figure 9b. As reference, the single drop flow in the straight channel (red color) and the sinusoidal channel (blue color) is plotted against the slug flow in sinusoidal channel. As can be noted, the drop's perimeter ratio P/P init of the slug flow differs to that of the single drop flow (blue color) when the liquid film thickness hits its maximum length; however, its evolution returns to the same behavior of the single drop flow. On the other hand, the liquid film thickness in the slug flow do not present significant variations if compared to the single drop flow in sinusoidal channel; however, it can be noted that the sinusoidal channel geometry reduces considerably the mean liquid film thickness when compared to the single drop in straight channel (red color). A=0.00, single A=0.14, single A=0.14, s=1.875D, rear A=0.14, s=1.875D, center A=0.14, s=1.875D, front A=0.00, single A=0.14, single A=0.14, s=1.875D, rear A=0.14, s=1.875D, center A=0.14, s=1.875D, front (b) slug s = 1.875 Figure 9. Single and multiple drops flow in channels with amplitudes A = 0.00 (straight channel, red color) and A = 0.14 (sinusoidal, blue color) in two corrugated lengths for (a) change in drop's perimeter ratio P relative to its initial perimeter P init and (b) film thickness evolution. The figures stand for slug flow with drop's initial slug lengths of s = 1.875. All parameters are dimensionless.
The axisymmetric streamlines ψ are plotted along with the velocity components v r = −(1/r)∂ψ/∂x and v x = (1/r)∂ψ/∂r and pressure p in Figure 10a to highlight the wake structure of center and front drop's before coalescence in slug flow with drop's initial slug lengths s = 1.875. In Figure 10b,c the velocity components, v x and v r , and pressure are interpolated axially between 9.3 ≤ x ≤ 9.7 at radius r = 0.1, respectively, showing the variation of drop's velocity and the pressure drop between the center and front drops. length [-] pressure (c) Figure 10. Wake structure of two consecutive drops before coalescence with drop's initial slug lengths of s = 1.875. (a) Axisymmetric streamlines ψ of the center and front drops, (b) velocity components v x and v r interpolated axially between 9.3 ≤ x ≤ 9.7 at radius r = 0.1, and (c) pressure distribution interpolated axially between 9.3 ≤ x ≤ 9.7 at radius r = 0.1. Figure 11 depicts the trend of the coalescence time t c relative to the initial slug length s for the straight and sinusoidal channels. As can be noted, for initial slug length s ≥ 0.875, the coalescence time increases linearly as s increases. Additionally, as noted previously, the sinusoidal channel anticipates the coalescence time for the same flow parameters of the straight channel. initial slug length s [-] straight channel sinusoidal channel Figure 11. Trend of the coalescence time relative to the initial slug length s for the straight and sinusoidal channels. For initial slug length s ≥ 0.875, the coalescence time increases linearly as s increases. All parameters are dimensionless.

Parametric study: wavelength λ
In the first parametric study, the wavelength λ has been systematically changed for values λ = {1, 2, 4, 6, 8} while keeping constant the drop's size factor κ d = 1.72 and all dimensionless numbers used at all simulations for two initial slug lengths s = {1.375, 1.875}. The standard test case for this parametric study is the straight channels with A = 0.00 in which the wavelength approaches infinity, where the coalescence time was t c = {500.268, 909.507}, respectively, for the initial slug lengths s = {1.375, 1.875}. In the sinusoidal channels, the wavelength parameters λ = {1, 2, 4, 6, 8} presented coalescence time of t c = {130.783, 771.516, 407.099, 571.755, 36.190}, respectively, indicating that the coalescence time is not linearly related to the channel's wavelength as can be seen in Figure 13.
In Figure 12, the final solution before coalescence is shown for all wavelength tested in this parametric study for initial slug length s = 1.375 and the respective coalescence time. For the wavelength λ = 2, no coalescence between the drops took place for time t < 1060; therefore, such a wavelength seems to be critical for the simulation parameters used in this work. All the others wavelengths have presented lower coalescence time relative to the straight channel with same parameters. Figure 13 reveals the coalescence time t c for the straight A = 0.00 and sinusoidal A = 0.14 channels relative to the wavelengths. Is is clear that the wavelength λ = 2 is critical for the presented simulation parameters where the coalescence time t c is much larger than the straight channels with same simulation parameters. On the other hand, all the others wavelengths anticipates the coalescence process when compared to the straight channels.  A=0.00, s=1.375 A=0.00, s=1.875 A=0.14, s=1.375 A=0.14, s=1.875 Figure 13. Trend of the coalescence time t c relative to the channel's wavelength λ for the straight A = 0.00 and the sinusoidal A = 0.14 channels. It is noted that for wavelength λ = 2, the coalescence time t c is higher than the coalescence time for the straight channel with same parameters. The others wavelengths have presented smaller coalescence time. All parameters are dimensionless.

Parametric study: drop's size factor κ d
In the last parametric study, the drop size factor is systematically changed for values κ d = {0.86, 1.29, 1.72, 2.15} while keeping constant the channel's wavelength λ = 4, initial slug length s = 1.375, and amplitude A = 0.14 for the corrugated channel as well as for the straight channel A = 0.00. In Figure 14, the center drop's perimeter ratio P/P init within two corrugated lengths is presented. As expected, the larger is the drop's size factor κ d , the larger is the perimeter ratio variation. A=0.14,κ=0.86,center drop A=0.14,κ=1.29,center drop A=0.14,κ=1.72,center drop A=0.14,κ=2.15,center drop Figure 14. Center drop's perimeter ratio P/P init within two corrugated lengths for sinusoidal with amplitude A = 0.14, initial slug length s = 1.375 and drop size factor for several drops: .86, 1.29, 1.72, 2.15}. As expected, the the larger is the drop's size factor κ d , the larger is the perimeter ratio variation. Figure 15 presents a parametric study of coalescence time t c between drops in slug flow for the sinusoidal channel with amplitude A = 0.14, initial slug length s = 1.375, and several drop's size factor κ d = 0.86, 1.29, 1.72, 2.15. For κ d = 0.86 the rear and center drops coalesce, where for the others, the coalescence takes place between the center and front drops. Before coalescence, the drop geometries are similar for κ d > 1.29 and it takes place in the same channel region. The drop's body increases and the drop's size factor increases.  Figure 15. Parametric study of coalescence time t c between drops in slug flow for the sinusoidal channel with amplitude A = 0.14, initial slug length s = 1.375 and several drop's size factor κ d = 0.86, 1.29, 1.72, 2.15. In subfigure (a), the rear and the center drops coalesce, where for the others the coalescence takes place between the center and front drops. All parameters are dimensionless. Figure 16 shows the trend of the coalescence time relative to the initial slug length s for the sinusoidal channel with amplitude A = 0.14. The drop's size factor κ d = 1.29 presented larger coalescence time t c for corrugated channels relative to the straight channels with same parameters. The others drop's size factors presented smaller coalescence time for the corrugated channel with amplitude A = 0.14 relative to the straight channel A = 0.00. Therefore, it is expected to have smaller coalescence time in corrugated channels. A=0.00 A=0.14 Figure 16. Trend of the coalescence time relative to the initial slug length s for the sinusoidal channel. The drop's size factor κ d = 1.29 presented larger coalescence time t c for corrugated channels relative to the straight channels with same parameters. The others drop's size factors presented smaller coalescence time for the corrugated channel with amplitude A = 0.14 relative to the straight channel A = 0.00. All parameters are dimensionless.

Conclusions
In this work, a numerical investigation was carried out to analyze the drop velocity, film thickness, and shape variation for one and three drops in buoyancy-driven motion for diethylene glycol/UCON-1145 (DEG3) two-phase system in straight and sinusoidal channels with dimensionless wavelength λ = 4 and wave amplitude A = 0.14. A state-ofthe-art model was employed to accurately compute the dynamics of the drop's interface motion using a modern moving frame/moving mesh technique within the arbitrary Lagrangian-Eulerian framework which allowed the simulation of very large domains. Several cases with different initial slug lengths s were simulated, and we concluded that the sinusoidal channel anticipates the coalescence of successive drops in slug flows when compared to the straight channel with same flow parameters. Additionally, two parametric studies have been carried out to evaluate drop coalescence time in slug flows varying wavelength λ and drop size factor κ d . It was observed that the sinusoidal channels with amplitude A = 0.14 is likely to decrease the coalescence time relative to the straight channels, except for wavelength λ = 2 and κ d = 1.29 where a different trend has been observed.
Funding: This research was funded by FAPERJ-Foundation for Research Support of the State of Rio de Janeiro.

Conflicts of Interest:
The author declares no conflict of interest.