Global stability analysis of two-fluid flows

Author: Luis Miguel Carrión Matamoros Supervisor: Miguel Ángel Herrada Gutiérrez Doctoral Thesis Departmento de Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla 2019 Global stability analysis of two-fluid flows Author: Luis Miguel Carrión Matamoros Supervisor: Miguel Ángel Herrada Gutiérrez i Doctoral Thesis Global stability analysis of two-fluid flows Author: Luis Miguel Carrión Matamoros Supervisor: Miguel Ángel Herrada Gutiérrez A thesis submitted in fulfillment of the requirements for the degree of Doctor of Philosophy Departamento de Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2019 ii Tesis doctoral: Global stability analysis of two-fluid flows Autor: Luis Miguel Carrión Matamoros Tutor: Miguel Ángel Herrada Gutiérrez El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2019 El Secretario del Tribunal iii A Dios, por su bendición. A mi esposa Daniela. A mis padres y hermanos. iv Acknowledgements Me gustaría agradecer a mi tutor y director Dr. Miguel Ángel Herrada Gutiérrez por brindarme la oportunidad de ser su estudiante de doctorado en esta prestigiosa institución. Gracias por impartir sus valiosos conocimientos en el mundo de la mecánica de fluidos y su simulación computacional, así como también por la paciencia para atender a cualquier duda que surgía en el transcurso de las investigaciones. Igualmente, a su familia, por su cálida atención y momentos agradables compartidos. Al Dr. Vladimir Shtern, por sus conocimientos en el estudio de la dinámica no-lineal y de estabilidad de los flujos de dos fluidos en rotación. Ha sido muy enriquecedor aprender de su experiencia científica. De la misma manera, quiero expresar mi agradecimiento al Dr. José María Montanero Fernández, por su apoyo en la investigación de la hidrodinámica y estabilidad de los puentes líquidos. Al Dr. José López Herrera, por las tertulias científicas en modelación matemática y simulación de flujos multifásicos. A mi esposa Daniela, por la paciencia y las palabras de apoyo que me ha brindado a cada momento. Gracias por estar a mi lado. v Resumen La presente tesis doctoral es un compendio de siete estudios dedicados al análisis numérico de la dinámica no lineal y la estabilidad lineal global de flujos de dos fluidos separados por una interfaz continua. Para esto se implementa la técnica numérica eficiente computacional desarrollada por Herrada y Montanero [1] para el estudio de la dinámica de sistemas capilares de un solo fluido. En vez de usar técnicas de aproximación para simplificar las ecuaciones de Navier-Stokes de flujos de dos fluidos obteniendo flujos básicos que dependen de una sola dirección, la aplicación de este método tiene la capacidad de modelar el conjunto completo de ecuaciones de gobierno de ambos fluidos para obtener el comportamiento hidrodinámico y el análisis de estabilidad lineal global. El primer capítulo es una introducción para explicar los diferentes tipos de flujos multifásicos y su importancia en la industria y la tecnología. Aquí, se presta especial atención a trabajos anteriores para el análisis de estabilidad lineal global de flujos de dos fluidos, en particular aquellos que contienen una interfaz libre. También se revisan los métodos numéricos para simularlos, así como los tipos de discretización espacial. El Capítulo 2 se centra en el estudio de un flujo de aire y agua como una aplicación en biorreactores aéreos. Este flujo estable de dos fluidos y axisimétrico es impulsado por un disco superior giratorio en un cilindro vertical sellado. La investigación demuestra la aparición de una burbuja de rotura de vórtice que aparece en el centro inferior y que se expande hacia la interfaz a medida que aumenta el número de Reynolds. Además, regiones de circulación meridional en el sentido de las agujas del reloj aparecen en el agua y en el aire, las cuales están separadas por una delgada capa de circulación en sentido contrario a las agujas del reloj adyacente a la interfaz en la región del agua. La novedad de este estudio muestra que el flujo se vuelve inestable para números de Reynolds mayores que aquellos en los que emerge la burbuja de rotura de vórtice y la capa delgada de circulación. En el capítulo 3 es analizado la topología del flujo y el comportamiento de estabilidad que se ocultan en el trabajo experimental desarrollado por Tsai et al. [2]. El disco giratorio superior de un cilindro vertical sellado, lleno de agua y aceite de soja, genera la circulación meridional y el flujo de rotación. El análisis numérico demuestra la forma de la interfaz de superficie plana y la aparición de la ruptura del vórtice en ambos fluidos al aumentar la fuerza de rotación. También se encontró que la inestabilidad del flujo estable axisimétrico que emerge en el agua es de tipo cizalla. Con respecto al capítulo 4, al igual que en el anterior, se analiza el flujo de rotación de dos fluidos que fue estudiado experimentalmente por Fujimoto y Takeda [3], con el objetivo de descubrir las topologías de flujo y las inestabilidades hidrodinámicas. En este caso, el contenedor cilíndrico vertical se llena con agua y aceite de silicona y sus movimientos se establecen mediante la tapa giratoria superior. El estudio numérico reproduce con precisión las complejas formas de la interfaz, llamadas en la investigación experimental como joroba, cúspide, monte Fuji y campana. Además de que la inestabilidad del flujo axisimétrico es de tipo cizalla, se encuentra que la energía de perturbación se concentra en la profundidad del agua, lo que explica la no dependencia de la forma de la interfaz en la inestabilidad. Los patrones de múltiples vórtices y la estabilidad de un flujo de rotación se investigan en el capítulo 5. A diferencia de los casos anteriores, se utiliza un disco giratorio inferior y aire y agua como fluidos de trabajo. En base a la fracción de volumen de agua, se encuentran diferentes escenarios topológicos con características interesantes, como cadenas de vórtices de aire y agua y vi vórtices de aire con forma de anillo de burbujas. Además, se demuestra que los movimientos de múltiples vórtices son estables y las inestabilidades aparecen cuando la interfaz se coloca cerca de la parte inferior o superior del contenedor. En relación al capítulo 6, este se dedica al análisis de un sistema de aire y agua axisimétrico y estable, impulsado por un disco superior giratorio en un contenedor semiesférico. En este estudio se demuestra que el aumento de la fuerza de rotación produce un patrón de tres vórtices, donde la circulación de aire y la co-circulación de agua están separadas por una región delgada de agua que gira en sentido contrario. Además, se encuentra que la perturbación de energía crítica se localiza en la región del aire el cual es de tipo cizalla. Finalmente, en los dos últimos capítulos, se estudia un puente líquido bajo los efectos de vibraciones de pequeña amplitud, surfactantes insolubles superficiales, gradientes térmicos y fuerza gravitacional. Específicamente, el capítulo 7 analiza la interacción entre las no homogeneidades de concentración de surfactante interfacial debido al flujo de transmisión y la convección de Marangoni. Como consecuencia, se encuentra que los patrones de flujo tienen un comportamiento no monótono al aumentar el número de Marangoni. Además, se ve que la intensidad del flujo de transmisión en la interfaz produce dos regímenes definidos cuando la amplitud de la fuerza se intensifica. En el capítulo 8, se estudia la estabilidad de la convección térmica para altos números de Prandtl con deformación de la superficie libre. Se muestra que las deformaciones dinámicas de la superficie libre del flujo base son muy parecidos a los experimentos previos y colapsan en una curva al normalizarse con su respectivo número capilar. Además, un efecto mínimo en el análisis de estabilidad es causado por la deformación interfacial. vii Abstract The present thesis is a compendium of seven studies devoted to the numerical analysis of nonlinear dynamics and global linear stability of two-fluid flows separated by a continuous interface. The computationally efficient numerical technique developed by Herrada and Montanero [1] is implemented to study the dynamics of capillary one-fluid systems. Instead of using approximation techniques to simplify the two-fluid flows Navier-Stokes equations, the application of this method has the capability of modelling the complete set of governing equations for both fluids in order to obtain the hydrodynamics behavior and the global linear stability analysis. The first chapter is an introduction to explain the different types of multiphase flows and their importance in the industry and technology. Here, special attention is paid to previous works for global linear stability analysis of two-fluid flows, in particular those with interfaces. It is also reviewed the numerical methods for simulating them as well as the spatial discretization techniques. Chapter 2 is focused on the study of an air-water swirling flow as an application in aerial bioreactors. This steady axisymmetric two-fluid flow is driven by a rotating top disk in a vertical sealed cylinder. It is found that a vortex breakdown bubble appears at the bottom center and spread out to the interface as Reynolds number increases. Then, clockwise meridional circulation regions appear in water and air, which are separated by a thin anticlockwise circulation layer adjacent to the interface in water. The novelty of this study shows that the flow becomes unstable for larger Reynolds numbers than those at which vortex breakdown bubble and thin circulation layer emerge. In chapter 3, the topological flow structure and stability behaviour hidden in the experimental work developed by Tsai et al. [2] is analysed. The upper rotating disk of a sealed vertical cylinder, filled with water and soybean oil, generates the meridional circulation and swirl flow. The numerical analysis proves the flattop interface shape and the vortex breakdown emergence in both fluids by increasing the rotation strength. It is also shown that the instability of the steady axisymmetric flow emerging in the water fluid is of the shear-layer type. With regard to chapter 4, as with to the previous one, the two-fluid swirling flow experimentally studied by Fujimoto and Takeda [3] is analyzed in order to unveil the flow topologies and hydrodynamic instabilities. In this case, the vertical cylindrical container is filled with water and silicone oil and their motions are set by the upper rotating lid. The numerical study accurately reproduces the complex interface shapes, called in the experimental investigation as hump, cusp, Mt. Fuji and bell. Besides the instability of the axisymmetric flow is of shear-layer type, it is found that the disturbance energy is concentrated in the water depth, explaining the no-dependency of interface shape on the instability. Multi-eddy patterns and stability of a whirlpool flow is investigated in chapter 5. Unlike the preceding cases, use is made of a bottom rotating disk and air and water as working fluids. Relying on the water volume fraction, it is found different topological scenarios with interesting characteristics such as zipper-like chains of air and water eddies and bubble-ring air eddies. In addition, it is shown that multi-eddy motions are stable and the instabilities appear when the interface is placed near either bottom or top of the container. Concerning chapter 6, it is dedidated to the analysis of a steady axisymmetric air-water driven by a rotating top disk in a semispherical container. It is found that the increasing of rotation strength produces a three-eddy pattern, where the air circulation and the water co-circulation are separated by a thin region of water counter-circulation. Moreover, the critical energy disturbance is xiv Figure 3.9 Distribution of swirl velocity at the interface at Reo = 700 (dashes), 850 (solid curve) and 1082 (dots). ......................................................................................................................................... 33 Figure 3.10 Distribution of oil velocity at the axis at Reo values shown near the curves. ............ 33 Figure 3.11 (color online) Patterns of meridional motion at Reo = 966 (a) and 1033 (b) show vortex breakdown cell in oil. .................................................................................................................... 34 Figure 3.12 (color online) (a) Meridional motion and (b) contours of constant energy of critical disturbances at Reo = 810 and m = 2. ..................................................................................................... 35 Figure 3.13 (a) z-distributions of bases-flow radial velocity vr and critical disturbance energy Ed .................................................................................................................................................................... 36 Figure 4.1 Schematic of water-spout model. .................................................................................... 39 Figure 4.2 Patterns of creeping flow at water height Hw = 0.9 (a), 0.93 (b) and 0.95 (c). ............ 42 Figure 4.3 Flow patterns at water height Hw = 1 and Reo = 1 (a), 50 (b) and 350 (c). ................. 43 Figure 4.4 Radial distribution at the interface of normalized (a) tangential-meridional vtn= vt/vtm and (b) swirl vn=v/vm velocities at Reo values shown on the plots. ................................ 44 Figure 4.5(a) Distribution of velocity at the axis for Reo = 350, (b) is a close up of (a) near the interface, z=zi = 1.0431. ........................................................................................................................ 45 Figure 4.6 Dependence on z of swirl (dashes) and radial (solid curve) velocities in water at r = 0.4 and Reo = 350. .................................................................................................................................... 45 Figure 4.7 Streamlines of basic flow (a)-(d) and contours of disturbance energy (e)-(f) for parameters listed in Table 3: Reo = 324 (a, e), 368 (b, f), 400 (c, g) and 614 (d, h). ........................... 47 Figure 4.8 Radial distribution of the basic-flow axial velocity (solid curve) and energy (dashed curve) of critical disturbance at z = 0.405, Reo = 324 and m = 4. ........................................................ 47 Figure 4.9 Surfaces of constant (a) azimuthal vorticity and (b) pressure of the critical disturbance at Reo = 324 and m = 4. ........................................................................................................................... 48 Figure 5.1 Geometry of the whirlpool model. .................................................................................. 52 Figure 5.2 Streamline patterns at Re = 350 (a) and 360 (b) show appearance of VB region CR3. .................................................................................................................................................................... 55 Figure 5.3 Distributions of velocity on the axis at Re = 350-400 illustrate (a) development of VB regions in water (CR3) and air (CR4) and (b) development of CR5 and its merging with CR4......... 56 Figure 5.4 Distributions of meridional velocity at the interface (vt) for Re from 360 up to 1300 illustrate (a) the first and (b) the second vt reversals near the axis. ....................................................... 56 Figure 5.5 Streamline patterns at Re = 380 (a) and 390 (b) show appearance of VB regions CR4 (a) and CR4* (b) in the air flow. .............................................................................................................. 57 Figure 5.6 Streamline patterns at Re = 400 (a) and 800 (b) show radial expansion of regions CR3 and CR4. .................................................................................................................................................... 57 Figure 5.7 Streamline patterns at Re = 1000 (a) and 1800 (b) correspond to the transformation of region CR5 from a bubble (a) into a bubble-ring (b). ............................................................................. 58 Figure 5.8 Schematics showing the flow topology in (a) Figure 5.7(a) and (b) in Figure 5.7(b). . 59 Figure 5.9 Distribution of velocity on the axis at Re = 2200 (solid curve) 2300 (dashes) and 2400 (dots). ......................................................................................................................................................... 60 Figure 5.10 Streamline patterns at Re = 2200 (a), 2400 (b), 2500 (c) and 2600 (b) show reduction of VB regions. ........................................................................................................................................... 61 xv Figure 5.11 Streamline pattern (a) and disturbance energy contours (b) at critical Re = 2629 for Hw= 0.8. .................................................................................................................................................... 62 Figure 5.12 Profiles of the disturbance energy (solid curve) and base-flow velocity (dashed curve) at z = zm. Both energy and velocity are normalized by their maximum values. ..................... 62 Figure 5.13 Energy contours for marginal disturbance with m = 2 (a) and 4 (b) at Hw = 0.8. ..... 63 Figure 5.14 Streamline patterns at Hw = 0.5 for Re = 250 (a) and 300 (b) show the outmost Moffatt eddy ME and the emergence of circulation regions CR3 and CR4. ....................................... 64 Figure 5.15 Distributions of velocity w at the axis, r = 0, for Re = 250 (dashed curve) and 300 (solid curve) show the emergence of circulation regions CR3 and CR4, see Fig. 14(b), in the water (air) region 0 < z < zi (zi < z < 1); zi = 0.497 (0.495) at Re = 250 (300). ............................................. 64 Figure 5.16 Profiles of meridional velocity on the interface vt at Re = 250, 300, and 400 illustrate the emergence and radial expansion of regions CR3 and CR4 (where vt > 0) as Re increases. ......... 65 Figure 5.17 Streamline patterns at (a) Re = 350 and (b) 400 shows (a) the emergence of circulation region CR4* and (b) merger of CR4 and CR4*. ................................................................. 65 Figure 5.18 Distributions of velocity w at the axis at Re = 300, 350, and 400 illustrate the emergence of circulation region CR4* (curves 300 and 350) and merging of CR4 and CR4* (curves 350 and 400, see also Fig. 17). ................................................................................................................. 66 Figure 5.19 Streamline patterns at (a) Re = 450 and (b) 1000 show the separation of region CR3 from axis and the emergence of region CR5. .......................................................................................... 67 Figure 5.20 Distribution of meridional velocity on the interface vt at Re = 400 (dotted curve), 450 (solid curve), and 1000 (dashed curve) shows the emergence (solid curve) and expansion (dashed curve) of region CR5 shown in Fig. 19(b). ............................................................................... 67 Figure 5.21Schematic of circulation regions CR1-CR6 at Re > Re5m = 1183 (see Figure 5.22). The bold line denotes the interface, the thin curve is separatrix, and the arrows show flow directions. .................................................................................................................................................................... 69 Figure 5.22 Streamline patterns at Re = 1500 corresponds to schematic in Figure 5.21. .............. 69 Figure 5.23Distribution of meridional velocity at the interface vt at Re = 1500 characterizes local circulation regions, CR3 (vt > 0), CR5 (vt < 0), and CR6 (vt> 0), adjacent to interface (Figure 5.21 and Figure 5.22). ....................................................................................................................................... 69 Figure 5.24 Streamline patterns at Re = 1700 (a), 1800 (b), 1900 (c) and 1950 (d). ..................... 70 Figure 5.25 Energy contours indicate where the critical disturbance mostly grows in the flow depicted in Fig. 24(d). ............................................................................................................................... 71 Figure 5.26 Streamline patterns at Re = 150 (a) and 200 illustrate the emergence of CR3. .......... 72 Figure 5.27 Streamline patterns at Re = 700 (a) and 900 (b) illustrate how CR2 transforms from bubble (Figure 5.26) in bubble-ring (a) and then in ring (b). Emergence of off-axis VB region CR4 is also shown. ................................................................................................................................................ 73 Figure 5.28 Profiles of velocity w on axis show how w > 0 range shrinks as Re increases from 200 (solid curve) to 250 (dashed curve). ................................................................................................. 74 Figure 5.29 (a) Velocity on axis at Re = 700, (b) close up of interface (zi) vicinity. .................... 75 Figure 5.30 (a) Distributions of swirl v and meridional vt velocities at the interface at Re = 800 (dotted curves), 900 (solid curves), and 1100 (dashed curves); (b) distributions of vt1/3 reveal the appearance of vt > 0 range, CR4 and CR5 as Re increases. .................................................................. 75 Figure 5.31 Schematics show the flow topology at Re = 829 (a) and Re > 829 (b). The bold lines denote the interface and the arrows indicate the flow direction............................................................. 76 xvi Figure 5.32 Distribution of meridional velocity on the interface at Re = 1400 and 1450 shows reversal of vt near axis and appearance of region CR6. ......................................................................... 77 Figure 5.33 Distribution of meridional velocity on the interface at Re = 1500 shows the appearance of region CR7. ....................................................................................................................... 78 Figure 5.34 Streamline patterns at Re = 1300 (a) and 1500 (b) show appearance of circulation region CR6; (c) is close-up of (b) visualizing region CR5. .................................................................... 78 Figure 5.35 Schematic shows the zipper-like topology of shallow whirlpool for Re  1500 (see Figure 5.33). The bold line is the interface and arrows indicate flow direction. .................................. 79 Figure 5.36 (a) Base-flow streamlines and (b) energy contours of critical disturbance at Hw = 0.3 and Re = 1547. ........................................................................................................................................... 80 Figure 5.37 (a) Global CR2 and (b) local CR5 kinds of bubble-ring cells. The bold curves depict the interface and the bold lines depict the sidewall. ................................................................................ 80 Figure 5.38 Dependence on Re of the maximal magnitudes of meridional vm and swirl vs velocities at the interface. Curve zi shows the interface height at the axis. Hw = 0.5. ....................... 82 Figure 6.1 Schematic of the problem. The lid only rotates. ............................................................ 86 Figure 6.2 Mesh generated using  w=25,  a=25,ns=20 with G(s) = 0.5-cos( s) and F(s)= f0s. ................................................................................................................................................... 89 Figure 6.3 (a) The base-flow pattern and (b) critical-disturbance energy distribution at Hw = 0.8. .................................................................................................................................................................... 90 Figure 6.4 Finding the critical Re by interpolation. The crosses present the numerical results. The curve and formula show the fitting polynomial. ..................................................................................... 90 Figure 6.5 (a) The z-profiles of the base-flow radial velocity u (solid curve) and Ed (dashed curve) near the lid at r = 0.618 Hw = 0.8. The dotted curve depicts the Kármán solution for u. (b) Contours of the disturbance axial velocity at the interface (color online). ............................................ 91 Figure 6.6 (a) The base-flow streamlines and (b) critical disturbance energy contours at Hw = 0.6. .............................................................................................................................................................. 92 Figure 6.7 (a) The z-profiles of the base-flow radial velocity u (solid curve) and Ed (dashed curve) in the air domain at r = 0.76 and Hw = 0.6. The dotted curve depicts the Kármán solution for u. (b) Contours of the disturbance axial velocity at the interface (color online). ................................ 92 Figure 6.8 (a) The base-flow streamlines and (b) critical disturbance energy contours at Hw = 0.4. .............................................................................................................................................................. 93 Figure 6.9 The z-profiles of the base-flow radial velocity u (solid curve) and Ed (dashed curve) in the air domain at r = 0.469 and Hw = 0.4. (b) Contours of the disturbance axial velocity at the interface (color online). ............................................................................................................................. 93 Figure 6.10 The r-profiles of the base-flow azimuthal v and meridional vt velocities at the interface for Hw = 0.4. .............................................................................................................................. 94 Figure 6.11 (a) The base-flow streamlines and (b) critical disturbance energy contours at Hw = 0.2. .............................................................................................................................................................. 94 Figure 6.12 (a) The r-profiles of the base-flow normalized axial velocity w (solid curve) and Ed (dashed curve) in the air domain at z = 0.289 and Hw = 0.2. (b) Contours of the disturbance axial velocity at the interface (color online). .................................................................................................... 95 Figure 7.1 Liquid bridge configuration ............................................................................................. 99 Figure 7.2 Time evolution of the axial component w of the velocity numerically calculated (solid xvii line) and of the axial component wL of the Lagrangian mean velocity (dots). These quantities were evaluated at the liquid bridge position (r = 1,z = 3Λ/2). The component wL was obtained from the difference between the fluid element positions at the instants t and t - T [Eq. (7.24)]. The simulation was conducted for Λ = 1.25, V = 1, C = 0.01, Ma = 0, ω = 2.1, and B = 0.005. .............. 101 Figure 7.3 (a) Streamlines of the mean Lagrangian velocity field.The cyan lines (first, third, and fifth vortex starting from the left) and gray lines (second, fourth, and sixth vortex starting from the left) indicate clockwise and anticlockwise rotations, respectively. The arrows are the mean Lagrangian velocities calculated from Eq. (7.24). (b) Magnification of plot (a) near the point (r = 0.85,z = 0.9).The black lines indicate the time-dependent trajectories of the fluid particles A and B calculated from the integration of Eqs. (7.23). (c and d) Magnification of plot (b) to show the trajectories of the fluid particles A and B. The subindexes i and f indicate the initial and final positions of the fluid particles in a cycle, respectively. The maximum value of the velocity field is vmax = 1.77 × 10−4. The velocity magnitude at the points P, A, and B are 0.92vmax, 0.98vmax, and 0.43vmax, respectively. The values of the governing parameters are the same as those in Figure 7.2, namely, Λ = 1.25, V = 1, C = 0.01, γ0 = 0.5, Ma = 0, ω = 2.1, and B = 0.005. ................................. 105 Figure 7.4 Mean free surface deflection fE(z) (left) and mean surfactant concentration γE(z) (right). The values of the governing parameters are the same as those in Fig. 2, namely, Λ = 1.25, V = 1, C = 0.01, γ0 = 0.5, Ma = 0, ω = 2.1, and B = 0.005. ....................................................................... 106 Figure 7.5 Counterpart of Fig. 3 for different Marangoni numbers as indicated by the labels in the figure. The values of the governing parameters are Λ = 1.25, V = 1, C = 0.01, γ0= 0.5, and ω = 2.1 .................................................................................................................................................................. 108 Figure 7.6 Mean surfactant concentration γE(z). The colors (and types of lines) correspond to the different Marangoni numbers as indicated in the left-hand graph. The values of the rest of governing parameters are Λ = 1.25, V = 1, C = 0.01, γ0 = 0.5, and ω = 2.1 ..................................... 108 Figure 7.7 Intensities Ioscbulk (a), Imfbulk (b), and Imfinterf (c) of the oscillatory flow, the mean flow in the bulk and the mean flow in the interface, respectively, for the Marangoni numbers considered in Figure 7.5. The values of the rest of governing parameters are Λ = 1.25, V = 1, C = 0.01, γ0 = 0.5, and ω = 2.1. .................................................................................................................. 109 Figure 7.8 Imfbulk (left) and Imfinterf (right) vs. Ioscbulkfor the indicated Marangoni numbers. The values of the rest of governing parameters are Λ = 1.25, V = 1, C = 0.01, γ0 = 0.5, and ω = 2.1. .................................................................................................................................................... 110 Figure 7.9 Counterpart of Figure 7.8 for the indicated oscillation frequencies. The values of the rest of governing parameters are Λ = 1.25, V = 1, C = 0.01, γ0 = 0.5, and Ma = 0.2. ...................... 110 Figure 7.10 Counterpart of Figure 7.8 for the indicated values of C. The values of the rest of governing parameters are Λ = 1.25, V = 1, γ0 = 0.5, Ma = 0.2, and ω = 2.1. .................................... 111 Figure 7.11 Ioscbulk and Imfbulk as a function of the dynamical Bond number B starting from equilibrium, as well as starting from the previous point while increasing (up) and decreasing (down) progressively the dynamical Bond number. The symbols overlap. The values of the governing parameters are Λ = 1.25, V = 1, C = 0.01, γ0 = 0.5, Ma = 0.2, and ω = 2.1. ..................................... 112 Figure 8.1 Liquid bridge configuration. The dynamical free surface deformation has been exaggerated. ............................................................................................................................................. 117 Figure 8.2 Details of the grid used in the simulations .................................................................... 121 Figure 8.3 Experimental (solid symbols) and numerical (open symbols) values of the critical Marangoni number Ma* (left-hand graph) and the oscillation frequency ωr* (right-hand graph) for liquid bridges with different aspect ratios AR under microgravity conditions (Fr=Gr=0) (see Fig. 6 of Ref. (Nishino et al., [204])). The circles and triangles correspond to the numerical results calculated with Bi = 0.15 and 0.3, respectively. ................................................................................... 123 xviii Figure 8.4 Streamlines (left) and temperature isolines (right) for AR=0.615,Pr=67,αν= 5.892,Fr=1.973×104,Θamb=0,Bi=0.15 and V=0.8265,Gr=326.9,Ma= 1.082x104,Ca=0.0508,T0*=1.645(left) and V=0.8298,Gr=468,Ma=1.549x104 Ca=0.0727,T0*=1.149 (right). The blue and cyan lines in the left-hand figure indicate clockwise and anticlockwise rotations, respectively. ........................................................................... 124 Figure 8.5 Experimental (open symbols) and numerical (solid symbols) free Surface deformation Fb-F0 of the base flow AR=0.615,Pr=67,αν=5.892,Fr=1.973×104,Θamb=0,Bi= 0.15 and V=0.8200,Gr=117.4,Ma=3.885x103,Ca=0.0182,T0*=4.579 (diamonds), V=0.8241,Gr=237.7,Ma=7.863x103,Ca=0.0369,T0*=2.262 (down-triangles), V=0.8265,Gr=326.9,Ma=1.082x104,Ca=0.0508,T0*=1.645 (up-triangles), V= 0.8298,Gr=468,Ma=1.549x104,Ca=0.0727,T0*=1.149 (circles) and V=0.8405,Gr= 544.8,Ma=1.803x104,Ca=0.0846,T0*=0.987 (squares). .................................................... 125 Figure 8.6 (Left) Free Surface deformation normalized with the Capillary number, Fb-F0/Ca, for the configurations considered in Figure 8.5: AR=0.615,Pr=67,αν=5.892,Fr=1.973× 104,Θamb=0,Bi=0.15 and V=0.8200,Gr=117.4,Ma=3.885x103,Ca=0.0182,T0*= 4.579 (diamonds), V=0.8241,Gr=237.7,Ma=7.863x103,Ca=0.0369,T0*=2.262 (down-triangles), V=0.8265,Gr=326.9,Ma=1.082x104,Ca=0.0508,T0*=1.645 (uptriangles), V=0.8298,Gr=468,Ma=1.549x104,Ca=0.0727,T0*=1.149 (circles) and V=0.8405,Gr=544.8,Ma=1.803x104,Ca=0.0846,T0*=0.987 (squares). ................... 126 Figure 8.7. (Left-hand image) Isolines of kinetic energy 1/2U2+V2+W2 (left) and magnitude of p of reduced pressure (right). (Right-hand image) Isolines of magnitude of the temperature field, Θ(r,z) (left), and amplitude F(z) of the free surface oscillation. The latter has been magnified to show clearly its axial dependence. The results correspond to the eigenmode m=1 responsible for instability at AR=0.615,Pr=67,αν=5.892,Fr=1.973×104,Θamb=0,Bi=0.15 and V=0.828,Gr=350.2,Ma=1.159x104,Ca=0.0544,T0*=1.535 ....................................... 126 Figure 8.8 Eigenvalues characterizing the m=1 (circles) and 2 (triangles) linear modes for the configurations considered in Figure 8.5: AR=0.615,Pr=67,αν=5.892,Fr=1.973× 104,Θamb=0,Bi=0.15 and V=0.8200,Gr=117.4,Ma=3.885x103,Ca=0.0182,T0*= 4.579 , V=0.8241,Gr=237.7,Ma=7.863x103,Ca=0.0369,T0*=2.262 , V= 0.8265,Gr=326.9,Ma=1.082x104,Ca=0.0508,T0*=1.645, V=0.8298,Gr= 468,Ma=1.549x104,Ca=0.0727,T0*=1.149 and V=0.8405,Gr=544.8,Ma= 1.803x104,Ca=0.0846,T0*=0.987. The open and solid symbols were calculated with dynamical free surface deformation and the equilibrium contour, respectively. ................................ 127 Figure 8.9 Spatial dependence of the temperature field perturbation and free surface oscillation amplitude for marginally stable liquid bridges with different aspect ratios. The calculations were conducted for 5-cSt silicone oil liquid bridges with L = 3.69 mm and V = 1. The free surface oscillation amplitude has been magnified to show clearly its axial dependence. ............................... 128 Figure 8.10 Critical Marangoni number Ma (circles) and oscillation frequency ωr(triangles) as a function of the aspect ratio AR for 5-cSt silicone oil liquid bridges with L = 3.69 mm and V = 1. The results were calculated with the dynamical free surface deformation (open symbols) and the rigid free surface approximation (solid symbols). ................................................................................................. 128 xix List of tables Table 2.1 Dependence of critical parameters on heavy fluid fraction Hh for fluids used by Brady at al. ............................................................................................................................................................ 15 Table 2.2 Dependence on heavy-fluid fraction Hh of radial rm and axial zm coordinates of absolute maximum of disturbance energy Ed. ........................................................................................ 17 Table 2.3 Dependence of marginal parameters for m = 1 on water fraction Hh in the air-water flow. ........................................................................................................................................................... 18 Table 2.4 Dependence on water fraction Hh of radial rm and axial zm coordinates of absolute maximum of disturbance energy Ed......................................................................................................... 18 Table 2.5. Dependence of marginal and critical (bold) parameters on water fraction Hh and azimuthal wave number m. ...................................................................................................................... 21 Table 2.6 Dependence on water fraction Hh of radial rm and axial zm coordinates of absolute maximum of disturbance energy Ed for critical Re (see Table 5). ......................................................... 23 Table 3.1 Marginal (i = 0) and critical (bold) stability characteristics. ......................................... 35 Table 4.1 Dependence on Reo of interface height zi, width of CRW1 at the axis zizs, maximal swirl vmi and minimal tangential-meridional vtm velocities at the interface and the Reynolds number of water flow Rew. ..................................................................................................................... 44 Table 4.2 Dependence of ω=ωr+iωi on Reo for the most dangerous disturbance at m = 4.. 46 Table 4.3 Critical (m = 4) and marginal (m = 1, 2 and 3) values of the Reynolds number (Reo) and frequency (ωr) at ωi = 0. .................................................................................................................. 46 Table 5.1 Dependence on Re of zi, zs, and CR5 thickness on axis (zs-zi) reveals that CR5 is a bubble-ring for Red6 ≤ Re ≤ Red7. ........................................................................................................ 59 Table 5.2 Critical (bold) and marginal parameters of the flow instability at Hw = 0.8. ................ 61 Table 5.3 Values of z for the interface (zi) and separatrix (zs) on the axis, and axial range (zs-zi) of CR5 for Re4m < Re < Re5m and of CR6 for Re > Re5m. ............................................................. 68 Table 5.4 Critical (bold) and marginal parameters of the flow instability at Hw = 0.5. ................ 71 Table 5.5 Axial values for the interface-axis (zi) and separatrix-axis (zs) intersections. Value of ( zs-zi) is the z-extent of CR2 at the axis for Re < 829. ......................................................................... 73 Table 5.6 Critical (bold) and marginal parameters of the flow instability at Hw = 0.3. ................ 79 Table 5.7 Changes in flow topology as Re increases at Hw = 0.8. ................................................. 81 Table 5.8 Dependence on water fraction of critical parameters of the whirlpool instability. ........ 83 1 1. INTRODUCTION ultiphase flows are systems in which different phases in solid, liquid or gas state coexist and interact simultaneously. These phenomena can be observed in a variety of industrial and technological applications, such as fluidized beds, cavitating flows, pollution controls, heat exchangers, aerosols, plasma spray coating, synthesis of nanophase materials, abrasive water-jet cutting, electrophotographic processes, to mention a few. Environmental phenomena are also typical multiphase flows, for instance, rain, snow, fog, avalanches, mudslides, sediment transport debris flows and countless more. Most biological and medical processes like blood flow, cell-fluid interface interaction, microbe locomotion, gas exchange, etc, have multiphase flows inherent characteristics. As mentioned above, multiphase flows are presented in a variety of processes, but their complex behavior makes the analysis and study a challenge task, due to mainly the interface interaction, type of fluid, change of phases, chemical reactions and size of particles. In order to explore these phenomena, it is necessary to determine the topologies [4] involved in, which can be disperse flows and separated flows. In disperse flows [5], small particles, bubbles or droplets are spread in a continuous fluid, whereas separated flows are composed by different continuous phases limited each other by interfaces [6], [7]. Figure 1.1 depicts the different types of multiphase flows and their industrial and technological applications. The definition of the topology is an important previous step to model the physical behavior of multiphase flows, which can be studied experimentally, theoretically and computationally. In the first one, it is used a laboratory equipment that is instrumented adequately so that the main physical characteristics can be obtained [8]; generally, empirical correlations are acquired from this experimental process. In the second approach, mathematical mechanistic equations [9], which are generally based on the main characteristics of the phenomenon, are utilized in order to model the physical behavior. Finally, computational approach studies the multiphase phenomena by applying numerical methods so as to solve the differential equations that govern them [10]. 1.1 Hydrodynamic linear stability of multiphase flows One of the most important branches of the fluid dynamics is the hydrodynamic stability, in which perturbations around a basic state flow are studied in order to analyze the transition from a stable to unstable condition[11], [12], . Basically, there are five means [13] to establish a hydrodynamic stability study: through the natural phenomena and laboratory experiments, numerical experiments, application of linear and weakly nonlinear theory, qualitative theory analysis of bifurcation and chaos and strongly nonlinear theory approach. The study in the present thesis is devoted to the analysis of modal linear stability of two-fluid flows with a continuos interface. In this approach, the variables of the flow system such as velocity, pressure or temperature are decomposed into a steady (base flow) and an unsteady part M 2 Global stability analysis of two-fluid flows (perturbation), the latter being approximated by small amplitude and phase functions. These anzats, represented by 𝐯(x,t)=𝐯(𝐱)+ε𝐯(x,t) for the velocity with ε≪1 and x,t as the spatial coordinates and time respectively, are replaced into the governing equations. The resulting equations are subtracted from the base flow and the terms in ε2 are dropped in order to yield the linearized perturbation equations. It is worth mentioning that the perturbation is considered as a superposition of many waves that do not interfere with each other. Figure 1.1 Types of multiphase flows [14] Based on the flow’s response to the perturbation or impulse, it can be identified both a stable and Introduction 3 an unstable flow. In the former one, the perturbation decays toward the base state in every direction from the source point of impulse. Conversely, if there is a growing effect of the perturbation along at least one ray emanating from the point of impulse, the flow is considered unstable. The modal linear stability theories can be classified as: local, nonlocal and global analysis [15]. Local stability analysis focuses on flows that vary slowly in one spatial direction or parallel flows [16], [17], where the basic flow is considered homogeneous along two spatial directions [18][13]. The first studies of hydrodynamics stability were focused in local linear analysis theory and some assuptioms were imposed in order to solve the differential equations that model the physical phenomenon. For instance, the pioneering work of Kelvin-Helmholtz (1871), which used the potential flow approach, analysed the stability of two fluids of different density in two horizontal parallel infinite streams; one of the natural phenomena that can be explained by using this theory is the generation of waves. From this study, by applying an additional upward vertical acceleration to the system, the Rayleigh-Taylor instability emerged and served to explain certain processes involved in salt domes, weather inversions and in phenomena related with astrophysics and electrohydrodynamics. Rayleigh [19], in other interesting work presented in 1878, used the Euler’s equations for inviscid fluids to study the linear stability leading to the break-up process of a round liquid jet under the surface tension effects. The phase function of the perturbation is composed of the phase velocity and the wave numbers. If this numbers are real or complex, then it can be distinguish three ways [20] of leading the local, nonlocal and global stability analysis, namely, (a) through temporal approach, from which the first classical hydrodynamics theories were developed [21][22], that uses the phase velocity as a complex number and real wave numbers; (b) by using spatial approach, where use is made of complex wave numbers and real phase velocity; and (c) by means of the spatio-temporal approach, where both the wave numbers and phase velocity are complex numbers. With respect to spatio-temporal stability analysis, special attention is given to unstable open flows, where a perturbation never passes the same point more than once. In this case, if the perturbation response grows along some rays from the point of impulse, the flow is called convectively unstable. On the other hand, if the perturbation grows at the point where it was applied, the flow is absolutely unstable. The reader interested in this subject is referred to Huerre and Monkewitz [23] for more details. The Figure 1.2 depicts a graphical description differentiation between stable and convectively/absolutely unstable flows. Figure 1.2 (a) Stable flow, (b) convectively unstable flow and (c) absolutely unstable flow [18]. In the nonlocal linear stability analysis, the amplitude of the perturbation, which depends only on one spatial direction for the local case, varies slowly in other one. This analysis emerged due to the fact that the local stability does not take into account both the growth of boundary layer and the 11 2. INSTABILITY ANALYSIS OF WATER-SPOT FLOW his chapter studies the hydrodynamic linear stability of a air-water swirling flow as a potential application in aerial bioreactors. This steady axisymmetric two-fluid flow is driven by the rotating top disk in a vertical sealed cylinder. A vortex breakdown bubble (VBB) appears at the bottom center and spread out to the interface as Reynolds number (Re) increases. Clockwise meridional circulation regions appear in water and air, which are separated by a thin anticlockwise circulation layer (TCL) adjacent to the interface in water. This study shows that the flow becomes instable for larger Re than those at which VBB and TCL emerge. Additionally, the instability is of shear-layer type in the air region but becomes centrifugal when the air volume fraction is small. 2.1. Introduction The study of the behavior and metamorphosis that occur in swirling flow is important because of industrial and natural applications [105], [106], [107], specially for the emergence of a local circulation cells called as vortex breakdown (VB). This phenomenon has been a focus of research owing to its relevance to delta-wing aircraft, furnaces and combustion chambers, water turbines, compressors, tornadoes, and other applications [108], [109], [110]. The first studies begun in 1957 with a pionner work of Peckham and Atkison [111] in which analysed the flow over a gothic wing in a low speed wind tunnel. Since then, numerous interesting hypotheses [110] has been developed in order to uncover the mechanism behind VB phenomenon. Recently, Shtern and Borissov [112] worked in swirl-decay mechanism so that to explain the VB emergence and the ways to be controlled it [113], [114]. Unlike VB studies in single fluid flows, the reseach in two-fluid flows have attracted attention in recent years because of applications in aerial bioreactors [51], [115], where air-water flow is driving by a propeller and or a rotating lid. The air motion expands towards the reactor axis, near the interface, and produces the meridional circulation and hence the rotation of water. The air flow delivers oxygen to the interface and helps tissue culture growth and water circulation improves the mixing process of the dissolved oxygen with other ingredients. The process of tissue growth riquires a considerable interval of time under small shear stresses and slow motion of ingredients. Because of the tissue fraction is small compared to that of water, it is neglected in the flow pattern studies [51], [115]. VB phenomena is found to be important in this type of biological process [51]. The interfacial interactions occurring in two-fluid swirling flows produce interesting characteristics which are not present in single-fluid flows. One of them is a thin anticlockwise circulation layer [116], [117] (TCL) adjacent to the interface. This remarkable phenomenon, having a large surface-area-to-volume ratio, enhances significantly the heat and mass transfer and can be important in bioreactors applications. Herrada et al. [118] studied a two-fluid flow driven by the rotating lid in a vertical sealed cylindrical container with same volumetric fractions of air and water, where the meridional circulation of air is clockwise type, which diverges from the axis near the lid and converges to the axis near the interface. Conversely, when the rotation lid is slow, the meridional circulation in the T 12 Global stability analysis of two-fluid flows water region is anticlockwise, converging to the axis near the interface and diverging from the axis near the bottom. The rising of speed rotation produces a vortex breakdown bubble (VBB) at the bottom center of the water region and spread out to the interface up to the thin layer anticlockwise circulation (TCL)[118]. There are vast studies related with the stability analysis of single-fluid VBB flows. For instance, Gelfgat et al. [119], [120] described the emergence of VBB in the steady axisymmetric flow as a consequence of local flow reversal without generating instabilities; therefore, there is not a direct relation between instability and the appareance of vortex breakdown. Escudier[121] and Sorensen et al. [122]–[124] found that, the emergence of VB in a steady axisymmetric flow, increasing Reynolds number, occurs for the container height-to-radius ratio 𝐻 less than 3.2. For larger values of H, the flow becomes unstable with respect to 3D time oscillatory disturbances, with 𝑚=3 for 3.2 < 𝐻 < 4.3, 𝑚=2 for 4.3 < 𝐻 < 5.2 and 𝑚=4 for 5.2 < 𝐻 < 5.5, being m the azimuthal wave number. Recently, Herrada et al. [113] argued that the instability nature of these flows is of the shear layer type. However, although there are detailed and rigorous instability studies of VBB single flows, no attention has been paid to the instability of two-fluid VB flow. Our study reveals that the disturbances only grow for Reynolds number being larger than those at which VBB and TCL emerge. To this end, a numerical technique has been developed for the stability investigation which can linearize efficiently the conditions at the interface. 2.2. Problem Formulation 2.2.1. Flow geometry Figure 2.1 Schematic of water-spout model. Figure 2.1 is a schematic of the problem. The lower part, 0 < 𝑧 < ℎ𝑤, of the cylindrical container of radius R and height h is filled with a heavy fluid (e.g., water), the upper part, ℎ𝑤 < 𝑧 < ℎ, is filled with a light fluid (e.g., air). With no motion, the interface is flat, 𝑧 = ℎ𝑤 (thin horizontal line in Fig. 1). When the top disk (at 𝑧 = ℎ) rotates with angular velocity  while the other walls are stationary, the interface becomes deformed as the curve in Figure 2.1 schematically shows. This flow mimics a water spout where swirling air raises water near the rotation axis. One control parameter is aspect ratio 𝐻 = ℎ/𝑅 which is fixed here: 𝐻 = 1. The dimensionless height of the interface at rest is 𝐻ℎ = ℎ𝑤/𝑅 being also a volumetric fraction of the heavy fluid which varies here. Other control parameters are the Reynolds number, 𝑅𝑒=𝜔𝑅2/𝜈𝑤, characterizing the swirl strength, the Froude number, 𝐹𝑟=𝜔2𝑅/𝑔=𝑎𝑅𝑒2, which is a centrifugal- to-gravity acceleration ratio, and the Weber number, 𝑊𝑒=𝜌𝑤𝜔2𝑅3/𝜎=𝑏𝑅𝑒2, characterizing the Instability analysis of water-spot flow 13 effect of surface tension 𝜎 at the interface; 𝑔 = 9.81 m2/s is the gravity acceleration, 𝜌𝑤 is the density of heavy fluid, 𝑎=𝜈𝑤2/(𝑔𝑅3) and 𝑏=𝜌𝑤𝜈𝑤2/(𝜎𝑅). As the light and heavy fluids are specified, coefficients a and b are fixed numbers while Re varies. We first consider the fluids, studied by Brady et al [116],[117], and imagine a physical experiment where 𝑅 = 5 mm, 𝑎 = 5.77105 and 𝑏 = 1.42105. Then we focus on the air-water flow[118] where 𝜈𝑎 = 15106 m2/s and 𝜈𝑤 = 106 m2/s are the kinematic viscosities, 𝜌𝑎 = 1.22 kg/m3 and 𝜌𝑤 = 1000 kg/m3 are densities; subscripts “𝑎” and “𝑤” means “air” and “water” as well as the light and heavy fluids. Pressure on the interface at rest is atmospheric, temperature 𝑇 = 300K and 𝜎 = 0.0715 kg/s2. Here we imagine a physical experiment where 𝑅 = 1 mm, 𝑎 = 1.02104 and 𝑏 = 1.39105. 2.2.2. Governing equations Using 𝑅, 𝜔𝑅, and 𝜌𝑤𝜔2𝑅2 as scales for length, velocity, and pressure, respectively, renders all variables dimensionless. We consider a flow of two viscous incompressible immiscible fluids governed by the Navier-Stokes equations [125], 1 𝑟𝜕(𝑟𝑢) 𝜕𝑟 +1 𝑟𝜕𝑣 𝜕𝜙+𝜕𝑤 𝜕𝑧=0, 2.1 𝜕𝑢 𝜕𝑡+𝑢𝜕𝑢 𝜕𝑟+𝑣 𝑟𝜕𝑢 𝜕𝜙+𝑤𝜕𝑢 𝜕𝑧−𝑣2 𝑟=−𝜌𝑛𝜕𝑝 𝜕𝑟+𝜈𝑛 𝑅𝑒(𝛻2𝑢−𝑢 𝑟2−2 𝑟2𝜕𝑣 𝜕𝜙) 2.2 𝜕𝑣 𝜕𝑡+𝑢𝜕𝑣 𝜕𝑟+𝑣 𝑟𝜕𝑣 𝜕𝜙+𝑤𝜕𝑣 𝜕𝑧+𝑢𝑣 𝑟=−𝜌𝑛 𝑟𝜕𝑝 𝜕𝜙+𝜈𝑛 𝑅𝑒(𝛻2𝑣−𝑣 𝑟2+2 𝑟2𝜕𝑢 𝜕𝜙) 2.3 𝜕𝑤 𝜕𝑡+𝑢𝜕𝑤 𝜕𝑟+𝑣 𝑟𝜕𝑤 𝜕𝜙+𝑤𝜕𝑤 𝜕𝑧=−𝜌𝑛𝜕𝑝 𝜕𝑧+𝜈𝑛 𝑅𝑒𝛻2𝑤 2.4 where 𝛻2≡1 𝑟𝜕 𝜕𝑟(𝑟𝜕 𝜕𝑟)+1 𝑟2𝜕2 𝜕𝜙2+𝜕2 𝜕𝑧2 is the Laplace operator for a scalar field, (𝑢,𝑣,𝑤) are the velocity components in cylindrical coordinates (𝑟,𝜙,𝑧), 𝑡 is time, and 𝑝 is pressure. The coefficients, 𝜌𝑛 and 𝜈𝑛, are both equal 1 at 𝑛 = 1 (in the heavy fluid) while 𝜌𝑛=𝜌𝑤/𝜌𝑎 and 𝜈𝑛=𝜈𝑎/𝜈𝑤 at 𝑛 = 2 (in the light fluid). We denote the list (𝑢,𝑣,𝑤,𝑝) as 𝑽, and look for a solution of the system 2.1-2.4 in the form 𝑽=𝑽𝑏(𝑟,𝑧)+𝜀𝑽𝑑(𝑟,𝑧)𝑒(𝑖𝑚𝜙−𝑖𝜔𝑡)+𝑐.𝑐., 2.5 where subscripts “b” and “d” denote the base flow and a disturbance, respectively; c.c. denotes the complex conjugate of the preceding term;  << 1 is an amplitude; integer 𝑚 is an azimuthal wave number; and 𝜔=𝜔𝑟+𝑖𝜔𝑖 is a complex number to be found, with frequency 𝜔𝑟 and growth rate of disturbance i. For a decaying (growing) disturbance, 𝜔𝑖 is negative (positive). The equations governing the base flow result from substituting 2.5 in system 2.1-2.4 and setting  = 0. The terms of order 𝑂(  ) constitute equations governing infinitesimal disturbances. 2.2.3. Boundary conditions Equations 2.1-2.4 are solved under the following boundary conditions: (i) Regularity at the axis, 0<𝑧<𝐻 , 𝑟=0: 14 Global stability analysis of two-fluid flows a) 𝑢=𝑣=0,𝜕𝑤/𝜕𝑟=0 (basic flow and 𝑚 = 0 disturbances), b) 𝑤𝑑=0,𝑢𝑑+𝑚𝑣𝑑=0,𝜕𝑢𝑑/𝜕𝑟=0 (𝑚 = 1 disturbances) c) 𝑤𝑑=𝑢𝑑=𝑣𝑑=0 (𝑚 >1 disturbances) (ii) No-slip at the walls: 𝑢=𝑣=𝑤=0 at the still disk, 0<𝑟<1, 𝑧=0, and at the sidewall, 0<𝑧<𝐻, 𝑟=1; 𝑢=𝑤=0, 𝑣=𝑟 at the rotating disk, 0<𝑟<1, 𝑧=𝐻. (iii) Continuity of all the velocity and stress components at the interface, 𝑧=𝐹(𝑟,𝜙,𝑡). In particular, the balance for the normal stresses yields that 𝑝𝑤−𝑝𝑎=1 𝑊𝑒𝛻.𝒏− 1 𝑅𝑒𝒏.(𝝉𝑤−𝜇𝑟𝝉𝑎).𝒏− 1 𝐹𝑟(1−𝜌𝑟)𝑧, 2.6 where n is the unit vector normal to the interface, 𝝉𝑤 and 𝝉𝑎 are tensors of the viscous stresses in the heavy and light fluids, respectively, 𝜇𝑟 and 𝜌𝑟 are the light-to-heavy fluid ratios of the dynamic viscosities and densities, respectively. (iv) The kinematic equation for the interface shape, 𝑧=𝐹(𝑟,𝜙,𝑡), yields that 𝜕𝐹 𝜕𝑡−𝑤+𝑢𝜕𝐹 𝜕𝑟+𝑣 𝑟𝜕𝐹 𝜕𝜙=0, 2.7 and we look for a solution in the form 𝑧=𝐹𝑏(𝑟,𝑧)+𝜀𝐹𝑑(𝑟,𝑧)𝑒(𝑖𝑚𝜙−𝑖𝜔𝑡)+𝑐.𝑐.. 2.3. Numeric technique To simulate the nonlinear problem for the axisymmetric basic flow and the generalized eigenvalue problem for disturbances, we use a numerical technique which is a variation of that described in detail in Ref. [1]. First, the heavy-fluid (e.g. water) and light-fluid (e.g. air) regions are mapped onto the standard square domain (0 ≤ (w,a) ≤ 1, 0 ≤  ≤ 1) by means of the coordinate transformations (a) w = z/F and  = r and (b)ηa = (z−F)/(H−F) and  = r for the water and air, respectively. Then, each variable (velocities, pressure field and the interface shape) and all its spatial and temporal derivatives, which appear in the transformed equations, are composed as a single symbolic vector. For example, for the axial velocity in the water flow we create a vector having 11 components: xw = [ww, ww/w, ww/, 2ww/2, 2ww/2, 2ww/, ww/, 2ww/2, 2ww/, 2ww/, ww/t]. The next step is to use a symbolic toolbox to calculate the analytical Jacobians of all the equations with respect to all the symbolic vectors. Using these analytical Jacobians we generate functions which then are evaluated point by point in the square domains. In this procedure, we used the MATLAB procedure matlabFunction to convert the symbolic Jacobians in MATLAB functions. Then, we carry out the spatial and temporal discretization of the problem. The water and air domains are discretized using a set of nw and na Chebychev spectral collocation points in the axial direction (along the w and a axes, respectively)[126]. The water and air domains are discretized using a set of n Chebychev spectral collocation points in the radial direction, . The second-order backward finite differences are used to compute the time derivatives for the basic flow. Since the basic flow is axisymmetric, all the azimuthal derivatives are set to zero. For disturbances, we obtain the temporal and azimuthal derivatives using the representation 2.5. The final step is to set up the numerical matrices allowing us to solve the problem by using a Newton procedure for the basic steady flow and by solving a generalized eigenvalue problem for disturbances. Details of this procedure are described in Ref. [1]. Instability analysis of water-spot flow 15 To summarize, the numeric procedure includes the mapping of water and air regions, the proper spatial and temporal discretization creating the discrete Jacobian matrix for the Newton procedure for the basic flow and two more matrices for the generalized eigenvalue problem for disturbances. For the basic flow, we get the final steady solution though an unsteady process. Starting from the rest and selecting a time step, dt, the solution is advanced in time until a steady state is reached. Since the nonlinear procedure used to compute the basic flow is fully implicit, dt can be taken sufficiently large to quickly reach the steady solution. Once the base flow is computed, and given an azimuthal wavenumber m, we use MATLAB subroutine eigs to calculate the eigenvalues () of the system of discrete linear equations. Most of the simulations presented here are done using nw = 25, na = 25, and n = 30 (standard grid), but some runs for flows, having small circulation regions, are performed with more fine grids. Since the Chebyshev grid points concentrate near the interface from both sides, the approach is adequate to resolve thin circulation layers, located near the interface, even at moderate values of nw, na, and n. The numerical code was validated by the agreement with the results of the paper [118] where the same base flow was simulated using a different code. The stability results where validated in the paper [1] by comparison with experimental studies. 2.4. Instability of liquid-liquid flow We first address the media studied by Brady et al. [116], where the light-to-heavy fluid density ratio is 𝜌𝑟 = 0.5284 and dynamic viscosity ratio is 𝜇𝑟 = 0.2. These values are larger than those for air and water by orders of magnitude that makes the interface remarkably deformed even at moderate Re values. We analyze stability features as the fraction of the heavy fluid 𝐻ℎ varies. We checked the disturbances with 𝑚 = 0, 1, 2, 3, and 4 and examined 𝜔𝑖. For small 𝑅𝑒, 𝜔𝑖 is negative. As 𝑅𝑒 increases, 𝜔𝑖 first become positive for 𝑚 = 1 disturbances for all 𝐻ℎ listed in Table 2.1 which presents the results for the critical (𝜔𝑖 = 0) values of the Reynolds number 𝑅𝑒, frequency 𝜔𝑟 and azimuthal wave number 𝑚 of the leading disturbances. To better understand the instability nature, we consider the distribution of disturbance kinetic energy 𝐸𝑑(𝑟,𝑧)=〈|𝑢𝑑|2+|𝑣𝑑|2+ |𝑤𝑑|2〉, normalized by its maximal value; < > denote averaging with respect 𝑡 and 𝜙. 𝑯𝒉 0.1 0.2 0.3 0. 4 0.5 0. 6 0. 7 0 .8 𝑹𝒆 97 0 134 6 105 0 9 16 85 0 7 79 6 90 6 30 𝝎𝒓 0.0 26 0.3 0.19 1 0. 221 0.2 47 0. 268 0. 298 0 .34 𝒎 1 1 1 1 1 1 1 1 Table 2.1 Dependence of critical parameters on heavy fluid fraction 𝐻ℎ for fluids used by Brady at al. Table 2.1 shows that as 𝐻ℎ decreases from 0.8 down to 0.2, 𝑅𝑒𝑐𝑟 grows while 𝜔𝑟 decreases. The physical reasoning for this trend likely is the following. The instability develops in the upper fluid Figure 2.2. As 𝐻ℎ decreases at a fixed 𝑅𝑒, the upper-fluid meridional motion weakens because its volume increases while its driving does not. This stabilizes the flow resulting in the 𝑅𝑒𝑐𝑟 growth. The limited container volume terminates this trend as the interface approaches the bottom. 16 Global stability analysis of two-fluid flows Figure 2.2 (color online) Streamline patterns at critical 𝑅𝑒 (see Table1) for heavy fluid fraction 𝐻ℎ = 0.1 (a), 0.5 (b) and 0.8 (c). The bold curves show the interface, light (dark) contours show clockwise (anticlockwise) meridional circulation. Contours of critical disturbance energy 𝐸𝑑 for 𝐻ℎ = 0.1 (d), 0.5 (e) and 0.8 (f); z-dependence of 𝐸𝑑, base-flow axial 𝑤𝑏 and radial 𝑢𝑏 velocities at 𝑟 = 𝑟𝑚 (Table 2.2) for 𝐻ℎ = 0.1 (g), 0.5 (h) and 0.8 (i). At small 𝐻ℎ, the jet-like boundary layer develops near the top disk, sidewall, interface and axis (Figure 2.2(a)). The shear-layer instability typical of jets results in the 𝑅𝑒𝑐𝑟 decrease at 𝐻ℎ = 0.1. The 𝜔𝑟 trend likely relates to the weakening of the upper-fluid rotation near the interface and axis as 𝐻ℎ decreases. This trend reverses at small 𝐻ℎ for the same reason as that for for Recr. The trend makes r negative at 𝐻ℎ = 0.2, i.e., the critical single-helix disturbances changes the direction of its rotation. For all 𝐻ℎ, the most dangerous disturbance is single-helix, 𝑚 = 1. The simulations at 𝐻ℎ = 0.9 reveal that the flow remains stable until the interface touches the rotating lid at 𝑟 = 0 as 𝑅𝑒 increases. Table 2.1 does not include the 𝐻ℎ = 0.9 case because our numeric technique does not allow 0 0.2 0.4 0.6 0.8 1 0 0.5 1 0 0.2 0.4 0.6 0.8 1 0 0.5 1 -0.75 -0.5 -0.25 0 0.25 0.5 0.75 1 0 0.5 1 (g) (h) (i) z z z Ed Ed Ed wb wb wb ub Instability analysis of water-spot flow 17 describing the interface disconnection from the container axis. The radial 𝑟𝑚 and axial 𝑧𝑚 locations of 𝐸𝑑 absolute maximum help understand the instability nature by indicating where disturbance energy 𝐸𝑑 focuses. Table 2.2 lists 𝑟𝑚 and 𝑧𝑚 values for the small (𝐻ℎ = 0.1), middle (𝐻ℎ = 0.5) and large (𝐻ℎ = 0.8) fractions of the heavy fluid. To help understand the instability nature, Figure 2.2 depicts the base-flow streamlines (a)-(c), contours of disturbance energy 𝐸𝑑 (d)-(f), and 𝐸𝑑 and base-flow velocity 𝑧-profiles at 𝑟 = 𝑟𝑚 (g)-(i) for 𝐻ℎ = 0.1 (the first column), 0.5 (the second column) and 0.8 (the third column) at critical 𝑅𝑒. The disturbance energy is normalized by its maximal value. 𝑯𝒉 0.1 0. 5 0. 8 𝒓𝒎 0.0 029 0 0. 898 𝒛𝒎 0.9 2 0. 92 0. 965 Table 2.2 Dependence on heavy-fluid fraction 𝐻ℎ of radial 𝑟𝑚 and axial 𝑧𝑚 coordinates of absolute maximum of disturbance energy 𝐸𝑑. As the locations of the absolute maximum of 𝐸𝑑 indicate in Figure 2.2(d)- Figure 2.2(f), the instability develops in the upper fluid. This is physically reasonable since (i) the light-fluid viscosity is smaller than the heavy-fluid viscosity and (ii) the light-fluid flow, being adjacent to the rotating lid, is faster than the heavy-fluid flow since the swirl velocity rapidly decays downward[113]. The velocity difference is especially large for the shallow heavy fluid at 𝐻ℎ = 0.1. Figure 2.2(g) shows that the axial velocity of base flow 𝑤𝑏 (normalized by its maximal value, 0.011) is significantly smaller in the lower fluid, 0 < 𝑧 < 𝑧𝑖 = 0.242, than in the upper fluid 𝑧𝑖 < 𝑧 < 1. In the upper fluid, 𝑤𝑏 is positive according to the clockwise circulation (Figure 2.2(a)) while in the lower fluid there is the small range, 0.225 𝑧 < 𝑧𝑖, where 𝑤𝑏 is negative in the anticlockwise circulation region (where contours are dark). The peak locations for disturbance energy 𝐸𝑑 and 𝑤𝑏 are very close in Figure 2.2(g) that indicates that the instability is likely related to the near-axis flow whose inflection point in the 𝑤𝑏(r) profile causes the shear-layer instability. This feature is similar to that in the single-fluid Vogel- Escudier flow[113]. The shear-layer instability also occurs at 𝐻ℎ = 0.5 as Figure 2.2(b), Figure 2.2(e) and Figure 2.2(h) illustrate. There are two peaks of 𝐸𝑑 located at the container axis in Figure 2.2(e) and two local peaks of 𝑤𝑏 which are close to the 𝐸𝑑 peaks in Figure 2.2(h). The character of instability changes as 𝐻ℎ increases up to 0.8 as Figure 2.2(c), Figure 2.2(f) and Figure 2.2(i) illustrate. The location of the 𝐸𝑑 absolute maximum shifts away from the axis and becomes close to the sidewall in Figure 2.2(f) with 𝑟𝑚 = 0.898 and 𝑧𝑚 = 0.965 (Table 2.2). This peak is located within the near-lid radial jet (Kármán’s boundary layer[127]) as Figure 2.2(i) illustrates where the radial velocity 𝑢𝑏(z) is depicted by the dashed curve (𝑢𝑏 and 𝑤𝑏 are normalized by maximal 𝑢𝑏 = 0.136). Figure 2.2(f) also reveals a local peak of 𝐸𝑑 which is close to the sidewall in the middle part of the upper-fluid region. This peak is likely associated with the centrifugal instability[22] which also occurs in the single-fluid Vogel-Escudier flow[113]. 2.5. Instability of air-water flow 2.5.1. Instability to m=1 disturbances Now we study the stability of the air-water flow where r = 0.0012 and r = 0.018 and therefore 18 Global stability analysis of two-fluid flows the interface deformation is smaller than that for the fluids used by Brady et al.[116] at close Re values. H h 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0. 8 0. 9 R e 62 718 67 260 99 045 84 165 57 492 43 140 18 960 10 950 2 520 R ea 41 82 44 84 66 03 56 11 38 33 28 76 12 64 73 0 1 68 R ew 40 6 52 8 88 2 70 1 51 0 40 7 20 4 12 3 6 0  r 0. 3282 0. 373 0.0 338 0. 488 0.0 146 0.0 194 0.0 372 0. 064 0. 263 m 1 1 1 1 1 1 1 1 1 Table 2.3 Dependence of marginal parameters for m = 1 on water fraction Hh in the air-water flow. Hh 0.1 0.5 0.9 rm 0.5 809 0.0 029 0,9 883 zm 0.2 060 0,6 754 0,9 894 Table 2.4 Dependence on water fraction Hh of radial rm and axial zm coordinates of absolute maximum of disturbance energy Ed. Figure 2.3, Table 2.3 and Table 2.4 show the characteristics of the air-water flow which are similar to those shown in Figure 2.2, Table 2.1 and Table 2.2 respectively. Following Sorensen at al.[123] we use terms “marginal” and “critical”. Neutral disturbances corresponding to the minimal Re at prescribed m are marginal. The marginal disturbance corresponding to the minimal Re among all m values is critical. Corresponding parameters Re, r and m are also named as marginal and critical. In Table 2.3, Re and Rea are based on the water and air viscosities respectively and Rew = ʋbimRe where ʋbim is the maximal value of the base-flow swirl velocity at the interface; Rea (Rew) characterizes the strength of the air (water) base flow. The air-water base-flow and its stability characteristics at m = 1 are similar to those for the fluids used by Brady et al.[116] at Hh = 0.1 and 0.5. The instability is likely of the shear-layer kind and develops in the upper fluid (air). Comparison of plots (g) and (h) in Fig. 3 as well as Rea and Rew values in Table 2.3 show that the water flow is drastically weaker than the air flow. In contrast, Figure 2.2 shows that the meridional flow strengths are of the same order of magnitude at marginal Re values for the heavy and light fluids studied by Brady et al.[116]. Figure 2.3 and Table 2.3 show that the Hh = 0.9 case qualitatively differs from the Hh = 0.1 and 0.5 cases with respect to both the basic flow pattern and its stability. At Hh = 0.1 (Figure 2.3(a)) and 0.5 (Figure 2.3(b)), the bulk air meridional flow is one-cellular and clockwise. It is separated from the water clockwise circulation by a thin layer of anticlockwise circulation (TCL) adjacent to the interface from below. In contrast, the air flow is two-cellular and the TCL is adjacent to the interface from above at Hh = 0.9 in Figure 2.3(c). Instability analysis of water-spot flow 19 Figure 2.3 Streamline patterns of air-water flow at critical Re values (see Table 3) for water fraction Hh = 0.1 (a), 0.5 (b) and 0.9 (c). The bold curve shows the interface, light (dark) contours show clockwise (anticlockwise) meridional circulation. Contours of marginal disturbance energy Ed for Hh = 0.1 (d), 0.5 (e) and 0.9 (f); z-dependence of Ed, base-flow axial wb and swirl ʋb velocities at r = rm (see Table 4) for Hh = 0.1 (g), 0.5 (h) and 0.9 (i). The Hh = 0.8 curve in Figure 2.4 shows that the radial velocity of base flow is negative, ub < 0, that corresponds to a one-cell air flow, as Figure 2.5(a) schematically depicts. As Hh increases, ub becomes positive near the axis, r = 0, up to point S, as curve Hh = 0.825 illustrates. Interpolation indicates that the ub > 0 region emerges at Hh = 0.8115. This corresponds to a change in the flow topology. As Figure 2.5(b) schematically depicts, the air flow becomes two-cellular with the anticlockwise circulation in region ACa and clockwise circulation in region CCa while the water anticlockwise circulation region ACw becomes separated from the axis. All flow cells meet at the saddle point S of the meridional air-water motion. 0 0.2 0.4 0.6 0.8 1 0 0.5 1 z ʋb 0 0.2 0.4 0.6 0.8 1 00.5 1 z ʋb 0 0.2 0.4 0.6 0.8 1 0.8 0.9 1 z wb ʋb Ed Ed Ed wb d wb d (g) (i) (h) 26 Global stability analysis of two-fluid flows 3.2.2. Governing equations Using 𝑅, 𝜔𝑅, and 𝜌𝑤𝜔2𝑅2 as scales for length, velocity, and pressure, respectively, renders all variables dimensionless. We consider a flow of two viscous incompressible immiscible fluids governed by the Navier-Stokes equations [17], 1 𝑟𝜕(𝑟𝑢) 𝜕𝑟 +1 𝑟𝜕𝑣 𝜕𝜙+𝜕𝑤 𝜕𝑧=0, 3.1 𝜕𝑢 𝜕𝑡+𝑢𝜕𝑢 𝜕𝑟+𝑣 𝑟𝜕𝑢 𝜕𝜙+𝑤𝜕𝑢 𝜕𝑧−𝑣2 𝑟=−𝜌𝑛𝜕𝑝 𝜕𝑟+𝜈𝑛 𝑅𝑒(𝛻2𝑢−𝑢 𝑟2−2 𝑟2𝜕𝑣 𝜕𝜙), 3.2 𝜕𝑣 𝜕𝑡+𝑢𝜕𝑣 𝜕𝑟+𝑣 𝑟𝜕𝑣 𝜕𝜙+𝑤𝜕𝑣 𝜕𝑧+𝑢𝑣 𝑟=−𝜌𝑛 𝑟𝜕𝑝 𝜕𝜙+𝜈𝑛 𝑅𝑒(𝛻2𝑣−𝑣 𝑟2+2 𝑟2𝜕𝑢 𝜕𝜙), 3.3 𝜕𝑤 𝜕𝑡+𝑢𝜕𝑤 𝜕𝑟+𝑣 𝑟𝜕𝑤 𝜕𝜙+𝑤𝜕𝑤 𝜕𝑧=−𝜌𝑛𝜕𝑝 𝜕𝑧+𝜈𝑛 𝑅𝑒𝛻2𝑤, 3.4 where 𝛻2≡1 𝑟𝜕 𝜕𝑟(𝑟𝜕 𝜕𝑟)+1 𝑟2𝜕2 𝜕𝜙2+𝜕2 𝜕𝑧2 is the Laplace operator for a scalar field, (𝑢,𝑣,𝑤) are the velocity components in cylindrical coordinates (𝑟,𝜙,𝑧), 𝑡 is time, and 𝑝 is pressure. The coefficients, 𝜌𝑛 and 𝜈𝑛, are both equal to 1 at 𝑛 = 1 (in water) while 𝜌𝑛=𝜌𝑤/𝜌𝑜 and 𝜈𝑛=𝜈𝑜/𝜈𝑤 at 𝑛 = 2 (in oil). We denote the list (𝑢,𝑣,𝑤,𝑝) as 𝑽, and look for a solution of the system 3.1 - 3.4 in the form 𝑽=𝑽𝑏(𝑟,𝑧)+𝜀𝑽𝑑(𝑟,𝑧)𝑒(𝑖𝑚𝜙−𝑖𝜔𝑡)+𝑐.𝑐., 3.5 where subscripts “𝑏” and “𝑑” denote the base flow and a disturbance, respectively; 𝑐.𝑐. denotes the complex conjugate of the preceding term;  << 1 is an amplitude; integer 𝑚 is an azimuthal wave number; and 𝜔=𝜔𝑟+𝑖𝜔𝑖 is a complex number to be found, with frequency 𝜔𝑟 and growth rate of disturbance 𝜔𝑖. For a decaying (growing) disturbance, 𝜔𝑖 is negative (positive). The equations governing the base flow follow from substituting 3.5 in system 3.1 - 3.4 and setting  = 0. The terms of order 𝑂(  ) constitute equations governing infinitesimal disturbances. 3.2.3. Boundary conditions Equations 3.1 - 3.4 are solved under the following boundary conditions: (i) Regularity at the axis, 0<𝑧<𝐻 , 𝑟=0: a) 𝑢=𝑣=0,𝜕𝑤/𝜕𝑟=0 (base flow and 𝑚 = 0 disturbances), b) 𝑤𝑑=0,𝑢𝑑+𝑣𝑑=0,𝜕𝑢𝑑/𝜕𝑟=0 (𝑚 = 1 disturbances) c) 𝑤𝑑=𝑢𝑑=𝑣𝑑=0 (𝑚 >1 disturbances) (ii) No-slip at the walls: 𝑢=𝑣=𝑤=0 at the still disk, 0<𝑟<1, 𝑧=0, and at the sidewall, 0<𝑧<𝐻, 𝑟=1; 𝑢=𝑤=0, 𝑣=𝑟 at the rotating disk, 0<𝑟<1, 𝑧=𝐻. (iii) Continuity of all the velocity and tangential stress components at the interface, 𝑧= 𝐹(𝑟,𝜙,𝑡). The balance for the normal stresses yields that 𝑝𝑤−𝑝𝑜=1 𝑊𝑒𝛻.𝒏− 1 𝑅𝑒𝒏.(𝝉𝑤−𝜇𝑟𝝉𝑜).𝒏− 1 𝐹𝑟(1−𝜌𝑟)𝑧, 3.6 Topology and stability of water-soybean-oil swirling flow 27 where 𝒏 is the unit vector perpendicular to the interface, 𝝉𝑤 and 𝝉𝑜 are tensors of the viscous stresses in the water and oil, respectively, and 𝜇𝑟 and 𝜌𝑟 are the light-to-heavy fluid ratios of the dynamic viscosities and densities, respectively. (iv) The kinematic equation for the interface shape, 𝑧=𝐹(𝑟,𝜙,𝑡), yields that 𝜕𝐹 𝜕𝑡−𝑤+𝑢𝜕𝐹 𝜕𝑟+𝑣 𝑟𝜕𝐹 𝜕𝜙=0, 3.7 and we look for a solution in the form 𝑧=𝐹𝑏(𝑟,𝑧)+𝜀𝐹𝑑(𝑟,𝑧)𝑒(𝑖𝑚𝜙−𝑖𝜔𝑡)+𝑐.𝑐.. 3.3. Numeric Technique To simulate the nonlinear problem for the axisymmetric base flow and the generalized eigenvalue problem for disturbances, we use a numerical technique which is a variation of that described in detail in Ref. [1]. First, the water and oil regions are mapped onto the standard square domain (0 ≤ (w,o) ≤ 1, 0 ≤  ≤ 1) by means of the coordinate transformations (a) w = z/F and  = r and (b) o = (z −F)/(H−F) and  = r for the water and oil, respectively. Then, each variable (velocities, pressure field, and the interface shape) and all its spatial and temporal derivatives, which appear in the transformed equations, are composed as a single symbolic vector. For example, for the axial velocity in the water flow, we create a vector having 11 components: xw = [ww, ww/w, ww/, 2ww/2, 2ww/2, 2ww/, ww/, 2ww/2, 2ww/, 2ww/, ww/t]. The next step is to use a symbolic toolbox to calculate the analytical Jacobians of all the equations with respect to all the symbolic vectors. Using these analytical Jacobians, we generate functions which then are evaluated point by point in the square domains. In this procedure, we used the MATLAB procedure matlabFunction to convert the symbolic Jacobians in MATLAB functions. Then, we carry out the spatial and temporal discretization of the problem. The water and oil domains are discretized using a set of nw and no Chebychev spectral collocation points in the axial direction (along the w and o axes, respectively) [126]. The water and oil domains are discretized using a set of n Chebychev spectral collocation points in the radial direction, . The second-order backward finite differences are used to compute the time derivatives for the base flow. Since the base flow is axisymmetric, all the azimuthal derivatives are set to zero. For disturbances, we obtain the temporal and azimuthal derivatives using the representation (3.5). This helps overcome the technical difficulty of linearizing relation (3.6) for the linear stability study of two-fluid flows. The final step is to set up the numerical matrices, allowing us to solve the problem by using a Newton procedure for the base steady flow and by solving a generalized eigenvalue problem for disturbances. Details of this procedure are described in Ref. [1] To summarize, the numeric procedure includes the mapping of water and oil regions, the proper spatial and temporal discretization creating the discrete Jacobian matrix of the Newton procedure for the base flow, and two more matrices for the generalized eigenvalue problem for disturbances. For the base flow, we get the final steady solution through an unsteady process. Starting from rest and selecting a time step, dt, the solution is advanced through time until a steady state is reached. Since the nonlinear procedure used to compute the base flow is fully implicit, a sufficiently large dt can be taken to quickly reach the steady solution. Once the base flow is computed and given an azimuthal wavenumber m, we use MATLAB subroutine eigs to calculate the eigenvalues () for the system of discrete linear equations. Most of the simulations presented here are done using nw = 60 and no = 25 in the axial direction and n = 40 in the radial direction (standard grid), but due to their small circulation regions, some simulations for flows at large Re, are performed with finer grids. The larger nw compared with no is used because the oil viscosity is significantly larger than the water viscosity. 28 Global stability analysis of two-fluid flows Since the Chebyshev grid points concentrate near the interface from both sides, the approach is adequate to resolve thin circulation layers located near the interface, even at moderate values of nw, no, and n. The Chebyshev grid points also concentrate near the axis and all walls, helping to resolve the small flow cells located near them. To analyze the flow topology, we use the Stokes stream function , w = r1/r and u = r1/z, and plot streamlines of the meridional motion, i.e., contours  = constant. 3.4. Topologycal metamorphoses of water flow 3.4.1. Moffatt and near-interface eddies Figure 3.2. (color online) Patterns of meridional motion at Reo = 1 (a), 40 (b) and 60 (c). It is instructive to start with a creeping flow and to gradually increase 𝑅𝑒𝑜. Figure 3.2 depicts streamlines at 𝑅𝑒𝑜 = 1 (a), 40 (b), and 60 (c). The bold (red online) curve denotes the interface. The clockwise (anticlockwise) circulation is shown hereafter by light (dark) contours which are blue online. Figure 2(a) depicts the pattern of the creeping flow. The rotating lid pushes the upper fluid to periphery near the lid and thus drives the clockwise circulation in CRO; CR is an abbreviation for “circulation region” and O denotes the oil. Near the interface, the oil moves toward the container axis. The water flow is driven by two factors: (i) swirl and (ii) meridional velocities at the interface, both provided by the oil motion. The factors (i) and (ii) push the water in opposite directions; the swirl (meridional) velocity tends to drive the clockwise (anticlockwise) meridional circulation. For 𝑅𝑒𝑜 = 1, the swirl velocity at the interface is very small, as discussed in more detail below. Therefore, factor (ii) dominates over (i) and drives the bulk anticlockwise circulation marked as CRW1 in Figure 3.2(a), where W denotes the water. Topology and stability of water-soybean-oil swirling flow 29 This circulation induces the Moffatt eddies near the sidewall-bottom intersection, as Figure 3.2(a) illustrates. Moffatt [132] revealed that a flow in a corner (between two inclined planes with a sufficiently small angle) has an infinite set of eddies whose strength and scale diminish to zero as the edge is approached. Since the outmost eddy size is small compared with the cylinder radius, see ME in Fig. 2(a), the sidewall can be approximated by a tangential plane thus reducing the cylindrical geometry to the corner between orthogonal planes. The corner angle, 90, is sufficiently small for the Moffatt eddies to occur [132]. Our numerical grid only resolves the outmost Moffatt eddy, denoted as ME in Figure 3.2(a). The Moffatt eddies also develop near the interface-sidewall intersection [133]. Our numerical grid only resolves the outmost near-interface eddy in water, marked as CRW2 in Figure 3.2(a). CRW2 drives the small counter-circulation cell right above CRW2 in the oil flow. In turn, this cell induces one, in Fig. 2(a), or two, in Figure 3.2(b) and Figure 3.2(c), near-sidewall small eddies. As 𝑅𝑒𝑜 increases, ME and CRW2 expand in z-direction, touch each other, as Figure 3.2(b) illustrates, and merge. The merged cell is denoted again as CRW2 in Figure 3.2(c). 3.4.2. Vortex breakdown in water flow The water moves upward in the bulk flow near the sidewall and thus blocks the downward transfer of the angular momentum from the interface. The angular momentum is transported toward the axis along the interface and then downward near the axis in the water domain. Therefore, the swirl effect is the strongest near the axis-bottom intersection where the VB cell emerges as Reo increases and factor (i) strengthens. The VB cell is denoted as CRW3 in Figure 3.3(c). Figure 3.3 (Color online) Patterns of meridional motion at Reo = 300 (a), 310 (b) and 370 (c) show formation of VB cell CRW3 in water. Figure 3.4 shows the distribution of velocity w at the axis in water at 𝑅𝑒𝑜 = 250 (dashes), 300 30 Global stability analysis of two-fluid flows (solid curve), and 305 (dots). To better observe small-magnitude velocity near the bottom, z = 0, Figure 3.4 plots 𝑤1/3. As 𝑅𝑒𝑜 increases, the 𝑤 > 0 range that corresponds to the VB water cell appears. Interpolating 𝜕𝑤/𝜕𝑧 at 𝑟 = 𝑧 = 0 as a function of 𝑅𝑒𝑜 yields that CRW3 emerges at 𝑅𝑒𝑜  300. Figure 3.5 depicts the dependence of the maximal value of swirl velocity at the interface, 𝑣𝑠𝑖, on 𝑅𝑒𝑜. This dependence helps understand why the VB in water occurs at 𝑅𝑒𝑜 = 300. At small 𝑅𝑒𝑜, the water rotation is very weak, 𝑣𝑠𝑖 = 0.0015. The 𝑣𝑠𝑖 value starts to significantly grow for 𝑅𝑒𝑜 > 100 and becomes nearly 0.07 at 𝑅𝑒𝑜 = 300, which corresponds to the inflection point of the curve in Figure 3.5 Figure 3.4 Distribution of velocity 𝑤 at the axis in the water domain at 𝑅𝑒𝑜 = 250 (dashes), 300 (solid curve), and 305 (dots). Figure 3.5 Dependence of maximal swirl velocity at the interface vsi on the Reynolds number. The Reynolds number, based on 𝑣𝑠𝑖 and the water viscosity 𝑣𝑤, is 𝑅𝑒𝑤 = 1370 at 𝑅𝑒𝑜 = 300, i.e., the water rotates sufficiently fast for the nonlinear terms to become important in equations (3.2)- (3.4). As the rotating fluid meets a normal wall, the secondary flow develops [134]. This flow radially converges to the rotation axis near the wall and goes away from the wall near the axis. The driving mechanism for the secondary flow is the radial gradient of pressure. It develops away from the wall to balance the centrifugal force: 𝜕𝑝/𝜕𝑟=𝜌𝑣2/𝑟. Near the wall, 𝜕𝑝/𝜕𝑟 is not balanced by the centrifugal force (because it has the second-order zero since 𝑣 = 0 at the wall) and pushes the fluid radially inward. This tornado-like secondary motion becomes sufficiently strong to overcome -0.1 -0.08 -0.06 -0.04 -0.02 0 0.02 0 0.05 0.1 0.15 0.2 0.25 𝑤1/3 𝑧 0.001 0.01 0.1 1 10 100 1000 𝑣𝑠𝑖 𝑅𝑒𝑜 Topology and stability of water-soybean-oil swirling flow 31 the anticlockwise circulation of water near the axis-bottom intersection at 𝑅𝑒𝑜 = 300. 3.4.3. Formation of thin anticlockwise circulation layer in water As 𝑅𝑒𝑜 further increases, the vortex breakdown cell CRW3 and the outmost Moffatt eddy ME, shown in Figure 3.3(c), merge and form the clockwise circulation, denoted again as CRW3 in Figure 3.6(a). The anticlockwise circulation CRW1 separates from the bottom, as Figure 3.6(a) illustrates at 𝑅𝑒𝑜 = 400, and becomes a thin layer, as Figure 3.6(b) illustrates at 𝑅𝑒𝑜 = 600. The zextent of CRW1 at the axis, 𝑟 = 0, becomes very small: 0.007 at 𝑅𝑒𝑜 = 600. Figure 3.6 (color online) Patterns of meridional motion at 𝑅𝑒𝑜 = 400 (a) and 600 (b) show formation of thin anticlockwise circulation layer in water. 3.4.4. Growth of the interface height at the axis This reduction results in the significant rise of the interface near the axis, driven by the water upflow in region CRW3, as Figure 3.6(b) illustrates. Figure 3.7 depicts this trend in more detail by showing the dependence of the interface height at the axis, 𝑧𝑖, on 𝑅𝑒𝑜. 32 Global stability analysis of two-fluid flows Figure 3.7 Dependence of interface height at the axis 𝑧𝑖 on the rotation speed 𝑅𝑒𝑜. 3.4.5. Formation of a bulge in the anticlockwise layer Figure 3.8 (color online) Patterns of meridional motion at 𝑅𝑒𝑜 = 700 (a) and 933 (b) show formation of a bulge in anticlockwise circulation layer in water. The interface rises in the axis vicinity but remains nearly flat for 𝑟 > 0.3, as Figure 3.8 illustrates. This causes an intriguing effect: the formation of a bulge in the anticlockwise water circulation. At 𝑅𝑒𝑜 = 700, the water motion is pretty fast, since 𝑅𝑒𝑤 = 3742. A strong jet forms near the interface. The interface peak at the axis makes this jet inclined downward. The inertia force separates the jet from the interface, where the interface becomes nearly flat, resulting in the bulge observed in Figure 3.8(b). The bulge divides the clockwise water circulation in the near-axis and near-sidewall parts in Figure 3.8(b). Figure 3.8(b) also shows that oil streamlines have a swelling near the axis for 𝑧 close to 1.2. This swelling is a precursor of vortex breakdown in the oil flow discussed below. 0.2 0.25 0.3 0.35 0.4 0 500 1000 𝑧𝑖 𝑅𝑒 𝑜 Topology and stability of water-soybean-oil swirling flow 33 3.5. Vortex breakdown in oil flow 3.5.1. Accumulation of swirl near the axis-interface intersection Figure 3.8 indicates the formation of a jet radially converging toward the axis near the interface in the oil flow. This intensifies the swirl near the axis-interface intersection. The dashed (Reo = 700) and solid (Reo = 850) curves in Figure 3.9 show that the swirl velocity v at the interface increases with Reo and its peak location shifts closer to the axis. This significantly increases the centrifugal acceleration, v2/r, and consequently the radial gradient of pressure, resulting in a local pressure minimum near the axis-interface intersection. Figure 3.9 Distribution of swirl velocity at the interface at 𝑅𝑒𝑜 = 700 (dashes), 850 (solid curve) and 1082 (dots). Figure 3.10 Distribution of oil velocity at the axis at 𝑅𝑒𝑜 values shown near the curves. 0 0.02 0.04 0.06 0.08 0.1 0 0.2 0.4 0.6 0.8 1 𝑣 𝑟 -0.02 -0.01 0 0.01 0.02 0.03 0.3 0.9 1.5 2.1 𝑤 𝑧 866 900 932 966 1032 34 Global stability analysis of two-fluid flows 3.5.2. Distribution of velocity at the axis in the oil flow The reduced pressure sucks the downstream oil, decelerating the upward oil motion near the axis. Figure 3.10 illustrates this trend by depicting the distribution of oil velocity at the axis at the 𝑅𝑒𝑜 values shown near the curves. A local minimum of 𝑤 forms, its location approaches the interface, and its value decreases as 𝑅𝑒𝑜 increases. 3.5.3. Emergence and development of vortex breakdown in the oil flow Interpolation of the local-minimal w value as a function of 𝑅𝑒𝑜 yields that vortex breakdown in oil occurs at 𝑅𝑒𝑜 = 941. The VB cell in oil emerges at the axis near 𝑧 = 1 and then extends in both radial and axial directions. Interpolating of 𝜕𝑤/𝜕𝑧 at the interface as a function of 𝑅𝑒𝑜 yields that the vortex breakdown cell reaches the interface at 𝑅𝑒𝑜 = 1000. Figure 3.11 depicts the vortex breakdown patterns at 𝑅𝑒𝑜 = 966 (a) and 1033 (b). As the VB oil cell reaches the interface, it meets the VB water cell. Next, the stagnation saddle point S separates from the axis, as Figure 3.11(b) illustrates. The anticlockwise circulation region in water detaches from the axis as well. All bulk flow regions meet at S for 𝑅𝑒𝑜 > 1000. The vortex breakdown patterns in the oil flow, depicted in Figure 3.11, agree with the experimental vortex breakdown patterns shown in Figure 3.2(b) and (b*c) of Ref. [2]. In contrast, the striking pattern of the water flow is not resolved in Ref. [2]. An advantage of numerical simulations is that they help investigate a flow in fine details. In particular, our study reveals the complicated structure of the water flow, which plays an important role in the development of the flat-top shape (named as “Mt. Fuji” shape in Ref. [3]) of the oil-water interface. Figure 3.11 (color online) Patterns of meridional motion at 𝑅𝑒𝑜 = 966 (a) and 1033 (b) show vortex breakdown cell in oil. Topology and stability of water-soybean-oil swirling flow 35 3.6. Stability analysis m 𝜔𝑟 Re o Re w vsi vtm rm zm 1 0.073 85 4 48 29 0.0827 -0.0342 0 0.0261 2 0.0097 81 0 43 71 0.0867 -0.0315 0.714 0.153 3 0.0241 84 2 46 89 0.0854 -0.0329 0.874 0.106 4 0.0388 84 9 47 58 0.0859 -0.0332 0.874 0.111 Table 3.1 Marginal (i = 0) and critical (bold) stability characteristics. We refer a reader to Refs. [1] and [135] for technical details of the stability study and discuss the results summarized in Table 3.1. It lists the characteristics corresponding to the marginal disturbances at the azimuthal wave number 𝑚 = 1, 2, 3, and 4. The marginal disturbances are neutral, 𝜔𝑖 = 0, corresponding to 𝑅𝑒𝑜 at which all other eigenvalues have 𝜔𝑖 < 0 for a prescribed 𝑚. Figure 3.12 (color online) (a) Meridional motion and (b) contours of constant energy of critical disturbances at 𝑅𝑒𝑜 = 810 and 𝑚 = 2. The critical disturbance is the marginal one, corresponding to the minimal 𝑅𝑒𝑜 for all 𝑚. The critical characteristics are shown bold in Table 3.1. The Reynolds number of the water flow, 𝑅𝑒𝑤, is based on the water viscosity, 𝜈𝑤, and the maximal azimuthal velocity at the interface, whose dimensionless value, 𝑣𝑠𝑖, is also presented in Table 3.1. The minimal meridional tangential velocity at the interface, 𝑣𝑡𝑚, is negative in Table 3.1 because it is directed toward the axis as Figure 3.12(a) illustrates. The location of peak marginal-disturbance squared velocity modulus, 𝐸𝑑, corresponds to 𝑟 = 𝑟𝑚 and 𝑧 = 𝑧𝑚 shown in the last two columns of Table 3.1. Here, 𝐸𝑑=〈|𝑢𝑑|2+|𝑣𝑑|2+ |𝑤𝑑|2〉 where the brackets denote averaging with respect to time and azimuthal angle. 𝐸𝑑 is normalized by its maximal value in the entire domain. A distribution of 𝐸𝑑 indicates where the disturbance magnitude is maximal indicating where it most production occurs. Other 𝑚 values not presented in Table 3.1 correspond to larger 𝑅𝑒𝑜. Thus, the flow becomes unstable with respect to time-oscillating three-dimensional disturbances for 𝑅𝑒𝑜 > 810. 42 Global stability analysis of two-fluid flows wavenumber m, we use MATLAB subroutine eigs to calculate the eigenvalues () of the system of discrete linear equations. Most of the simulations presented here are done using nw = 50, no = 25, and n = 40 (standard grid), but some runs for flows, having small circulation regions, are performed with more fine grids. Since the Chebyshev grid points concentrate near the interface from both sides, the approach is adequate to resolve thin circulation layers, located near the interface, even at moderate values of nw, no, and n. To analyze the flow topology, we use the Stokes stream function , w = r1/r and u = r1/z, and plot streamlines of the meridional motion, i.e., contours  = constant. 4.4. Topological metamorphoses of water flow 4.4.1. Creeping flow It is instructive to start with a creeping flow, where the flow topology only depends on the fluid fractions. Figure 4.2 depicts streamlines at 𝑅𝑒 = 1, 𝐻 = 2.3 and 𝐻𝑤 = 0.9 (a), 0.93 (b), 0.95 (c). Figure 4.2 Patterns of creeping flow at water height 𝐻𝑤 = 0.9 (a), 0.93 (b) and 0.95 (c). Figure 4.2(a) shows the simplest flow pattern. The rotating lid pushes the silicon oil to periphery and thus drives the clockwise circulation in CRO; CR is an abbreviation for “circulation region” and O denotes the oil. Near the interface, the oil moves toward the container axis. The water flow is driven by two factors: (a) swirl and (b) meridional velocities at the interface, both provided by the oil motion. The factors (a) and (b) pushes the water in opposite directions: the swirl (meridional) velocity tends to drive the clockwise (anticlockwise) circulation. For small 𝐻𝑤, factor (b) dominates [114], and the bulk water circulation is anticlockwise as Figure 4.2(a) shows in CRW1; here W denotes the water flow. The water circulation induces the Moffatt eddies (ME) near the sidewall-bottom intersection. Moffatt [132] revealed that a flow in a corner (between two inclined planes with a sufficiently small Metamorphosis in a water-oil swirling flow 43 angle) has an infinite set of eddies whose strength and scale diminish to zero as the edge is approached. Since the outmost eddy size is small compared with the cylinder radius, e.g., see ME in Figure 4.2(a), the sidewall can be approximated by a tangential plane thus reducing the cylindrical geometry to the corner between orthogonal planes. The corner angle, 90, is sufficiently small for the Moffatt eddies to occur. Our numerical grid only resolves the outmost Moffatt eddy: ME in Figure 4.2(a). The water moves upward to the interface near the sidewall and thus blocks the downward transfer of the angular momentum there. The angular momentum is transported toward the axis near the interface and then downward near the axis. Therefore, the swirl effect is the strongest near the axisbottom intersection where, as 𝐻𝑤 increases, factor (a) generates CRW2 shown in Figure 4.2(b). Interpolating the velocity distribution at the axis yields that CRW2 emerges at 𝐻𝑤 being very close to 0.93 that explains why CRW2 is so small in Figure 4.2(b). As 𝐻𝑤 further increases, CRW2 and ME expand. They merge into a single region denoted again as CRW2 in Figure 4.2(c). This feature is similar to that described in Ref. [130]. The merger occurs at 𝐻𝑤 = 0.938 as interpolation of w(z) yields at r = 0.4. 4.4.2. Topological changes as rotation speeds up Now we fix 𝐻𝑤 = 1 (as in Ref. [3]) and increase 𝑅𝑒. For convenient comparison with the results of Ref. [3], we introduce 𝑅𝑒𝑜=𝜔𝑅2/𝜈𝑜 = 𝑅𝑒/137.5 based on the oil viscosity. Figure 4.3 depicts the flow patterns at 𝑅𝑒𝑜 = 1 (a), 𝑅𝑒𝑜 = 50 (b) and 𝑅𝑒𝑜 = 350 (c). As 𝑅𝑒 increases, CRW2 expands upward and CRW1 shrinks toward the interface as Figure 4.3 illustrates. Table 4.1 lists the dependence on 𝑅𝑒𝑜 of the following flow characteristics: the interface height at the axis 𝑧𝑖, the CRW1 width at the axis 𝑧𝑖𝑧𝑠 , see Figure 4.3(b), the maximal swirl velocity at the interface 𝑣𝑚𝑖, the minimal tangential-meridional velocity at the interface 𝑣𝑡𝑚, and the Reynolds number of water flow 𝑅𝑒𝑤=𝑣𝑚𝑖 𝑅𝑒; 𝑅𝑒𝑤 is based on the maximal swirl velocity at the interface and water viscosity. Figure 4.3 Flow patterns at water height 𝐻𝑤 = 1 and 𝑅𝑒𝑜 = 1 (a), 50 (b) and 350 (c). 44 Global stability analysis of two-fluid flows 𝑅𝑒𝑜 𝑧𝑖 𝑧𝑖𝑧𝑠 𝑣𝑚𝑖 𝑣𝑡𝑚 𝑅𝑒𝑤 0 1 0.47719 0.01022 -0.00025 0 50 1.0004 0.23842 0.01040 -0.00793 73 100 1.0011 0.12863 0.01235 -0.00933 173 150 1.0021 0.07607 0.01647 -0.01115 346 200 1.0044 0.04343 0.02347 -0.01429 657 250 1.0097 0.02296 0.03423 -0.01847 1198 300 1.0212 0.00495 0.04817 -0.02333 2023 350 1.0431 0 0.06296 -0.02692 3085 400 1.0777 0 0.07571 -0.03451 4240 Table 4.1 Dependence on 𝑅𝑒𝑜 of interface height 𝑧𝑖, width of CRW1 at the axis 𝑧𝑖𝑧𝑠, maximal swirl 𝑣𝑚𝑖 and minimal tangential-meridional 𝑣𝑡𝑚 velocities at the interface and the Reynolds number of water flow 𝑅𝑒𝑤. Figure 4.4 depicts the radial profiles of (a) the tangential-meridional and (b) swirl velocities at the interface, normalized by their respectively minimal and maximal values, for the 𝑅𝑒𝑜 values indicated on the plots. Figure 4.4(a) reveals that the profile of the meridional velocity is nearly Reindependent. Figure 4.4(b) shows that the location of the maximal swirl velocity shifts closer to the axis for large 𝑅𝑒𝑜. Figure 4.4 Radial distribution at the interface of normalized (a) tangential-meridional 𝑣𝑡𝑛= 𝑣𝑡/𝑣𝑡𝑚 and (b) swirl 𝑣𝑛=𝑣/𝑣𝑚 velocities at 𝑅𝑒𝑜 values shown on the plots. 4.4.3. Formation of a robust bubble-ring A striking feature is the shape of CRW1 for large Re. Figure 4.3(c) shows that CRW1 becomes a thin circulation layer (TCL) adjacent to the interface. The third column in Table 4.1 lists values of 𝑧𝑖𝑧𝑠 which is the 𝑧-extent at the axis of CRW1. As 𝑅𝑒 increases, 𝑧𝑖𝑧𝑠 reduces, becomes zero at 𝑅𝑒𝑜 = 314 and remains zero for larger 𝑅𝑒𝑜. Therefore, TCL topologically is a bubble for 𝑅𝑒𝑜 < 314 and a bubble-ring for 𝑅𝑒𝑜  314. Typically, a bubble-ring is a transient shape. A small variation of a control parameter transforms a bubble-ring either in a bubble or in a ring. However here, CRW1 is a 0 0.2 0.4 0.6 0.8 1 00.2 0.4 0.6 0.8 1 𝑣𝑡𝑛 r 50 100 200 400 (a) 0 0.2 0.4 0.6 0.8 1 0.0 0.2 0.4 0.6 0.8 1.0 𝑣𝑛 r 50 100 200 400 (b) Metamorphosis in a water-oil swirling flow 45 bubble-ring in the range of 𝑅𝑒𝑜 > 314 shown in Table 4.1, i.e., TCL is a robust bubble-ring (RBR). Figure 4.5 confirms that TCL is a bubble-ring at 𝑅𝑒𝑜 = 350 by depicting the distribution of velocity at the axis. The dots show w(z) in the water region and the solid curve shows 𝑤(𝑧) in the oil region. Figure 4.5(b) is a close-up of Fig. 5(a) near the interface. Figure 4.5 shows that 𝑤(𝑧) is zero but not tangent to the 𝑤 = 0 line at 𝑧=𝑧𝑖. The derivative, 𝜕𝑤/𝜕𝑧, has a jump at 𝑧=𝑧𝑖 being negative (positive) in the water (oil) region. The jump is due to different densities and viscosities of the oil and water. Figure 4.5(a) Distribution of velocity at the axis for 𝑅𝑒𝑜 = 350, (b) is a close up of (a) near the interface, 𝑧=𝑧𝑖 = 1.0431. Figure 4.6 Dependence on z of swirl (dashes) and radial (solid curve) velocities in water at 𝑟 = 0.4 and 𝑅𝑒𝑜 = 350. The TCL is a region of large shear as Figure 4.6 illustrates depicting profiles of swirl 𝑣(𝑧) (dashes) and radial 𝑢(𝑧) (sold curve) velocities in water at 𝑟 = 0.4 and 𝑅𝑒𝑜 = 350. TCL, shown in Fig. 3(c), is located between the peak and valley of 𝑢(𝑧) near 𝑧=𝑧𝑖 = 1.0431 in Figure 4.6. The swirl velocity also significantly changes within TCL. The magnitudes of 𝜕𝑢/𝜕𝑧 and 𝜕𝑣/𝜕𝑧 are maximal at the interface, 𝑧=𝑧𝑖. These large shear stresses likely contribute to the striking shapes of the interface shown in Fig. 2 of Ref. [3]. However, before these shapes develop, the steady flow can suffer from an instability, because the Reynolds number of the water flow 𝑅𝑒𝑤 is rather large at 𝑅𝑒𝑜  350 according to the values of 𝑅𝑒𝑤 in the last column of Table 4.1. In addition, the water flow includes a jet forming near the interface as Figure 4.3(c) illustrates at 𝑅𝑒𝑜 = 350. This jet forms below TCL, 0 0.02 0.04 0.06 0.08 0.1 0123 𝑤 𝑧 water oil (a) 0 0.002 0.004 0.006 0.008 0.01 1 1.025 1.05 1.075 1.1 𝑤 𝑧 water oil (b) -0.03 0.00 0.03 0.06 00.35 0.7 1.05 𝑢, 𝑣 𝑧 u v 46 Global stability analysis of two-fluid flows goes to the sidewall and then downward near the sidewall and toward the axis near the bottom. These features are similar to those observed in Ref. [135]. The solid curve in Figure 4.6 has a number of inflection points that typically causes the shear-layer instability. This motivates us to perform a stability study, whose results are reported below. 4.5. Stability analysis 4.5.1. Marginal and critical characteristics Our simulations show that the most dangerous disturbance corresponds to the azimuthal wave number 𝑚 = 4. Table 4.2 lists eigenvalues, 𝜔=𝜔𝑟+𝑖𝜔𝑖, for a few values of 𝑅𝑒𝑜 at 𝑚 = 4. Interpolation yields the critical values: 𝑅𝑒𝑜 = 324, 𝜔𝑟 = 0.0391 and 𝜔𝑖 = 0. Disturbances with 𝑚 = 1, 2, and 3 becomes growing for larger 𝑅𝑒𝑜 as Table 4.3 shows listing the marginal (corresponding to 𝜔𝑖 = 0) values. For 𝑚 > 4, marginal 𝑅𝑒𝑜 values also are larger than the critical Reo at 𝑚 = 4. 𝑅𝑒𝑜 300 336 368 400 𝜔𝑟 0.0358 0.0407 0.0441 0.0465 𝜔𝑖 0.000464 0.000245 0.000927 0.001536 Table 4.2 Dependence of 𝜔=𝜔𝑟+𝑖𝜔𝑖 on 𝑅𝑒𝑜 for the most dangerous disturbance at 𝑚 = 4. 𝑅𝑒𝑜 324 368 400 614 𝜔𝑟 0.0391 0.0318 0.0381 0.0244 𝑚 4 3 2 1 Table 4.3 Critical (𝑚 = 4) and marginal (𝑚 = 1, 2 and 3) values of the Reynolds number (𝑅𝑒𝑜) and frequency (𝜔𝑟) at 𝜔𝑖 = 0. 4.5.2. Basic-flow streamlines and disturbance energy contours Metamorphosis in a water-oil swirling flow 47 Figure 4.7 Streamlines of basic flow (a)-(d) and contours of disturbance energy (e)-(f) for parameters listed in Table 3: 𝑅𝑒𝑜 = 324 (a, e), 368 (b, f), 400 (c, g) and 614 (d, h). Figure 4.7 depicts streamlines of the basic meridional motion, (a)-(d), and contours of disturbance energy 𝐸𝑑 = constant for critical, (a, e), and marginal parameter values listed in Table 4.3. Here 𝐸𝑑=〈|𝑢𝑑|2+|𝑣𝑑|2+|𝑤𝑑|2〉 where the brackets denote averaging with respect to time and azimuthal angle and 𝐸𝑑 is normalized by its maximal value in the entire domain. At critical 𝑅𝑒𝑜 = 324, the energy peak is located near the middle of the water domain as Figure 4.7(e) shows. The outmost energy contour in the water domain corresponds to 𝐸𝑑 = 0.1, i.e., the disturbance energy is focused between the water down-flow, going along the sidewall, and the water up-flow, going near 𝑟  0.5 as the pair of arrows indicate in Figure 4.7(a). 4.5.3. Shear-layer character of instability Figure 4.8 Radial distribution of the basic-flow axial velocity (solid curve) and energy (dashed curve) of critical disturbance at 𝑧 = 0.405, 𝑅𝑒𝑜 = 324 and 𝑚 = 4. -1.5 -1 -0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 𝑤𝑏, 𝐸𝑑 𝑟 𝑤𝑏 𝐸𝑑 48 Global stability analysis of two-fluid flows Figure 4.8 shows this in more detail by depicting the radial profiles of the basic-flow axial velocity 𝑤𝑏(𝑟) (solid curve) and the energy of critical disturbances 𝐸𝑑(𝑟) (dashed curve) at 𝑧 = 𝑧𝑚 = 0.405 (𝑧𝑚 is the axial coordinate of the 𝐸𝑑 peak in Figure 4.7(e)), 𝑅𝑒𝑜 = 324 and 𝑚 = 4. The velocity is normalized by its maximal value for the convenient comparison with 𝐸𝑑(𝑟) in Figure 4.8. The radial location of 𝐸𝑑 peak, 𝑟 = 𝑟𝑒 = 0.714, is close to the inflection point of 𝑤𝑏(𝑟) located at 𝑟 = 𝑟𝑖 = 0.76. This indicates that the instability likely is of shear-layer type, corresponding to the inflection point in the counterflow marked in Figure 4.7(a) by the arrows in the water domain. A smaller energy peak exists in the oil domain near the interface, which likely can be related to the radial counterflow in the thin circulation layer shown by the dark curves in Figure 4.7(a). As 𝑅𝑒𝑜 increases, the 𝐸𝑑 peaks shift toward the axis both in the water and oil domains. An interesting feature is the development of vortex breakdown bubble (VBB) in the oil domain shown by the dark contours in Fig. 7(d). This VBB development looks similar to that experimentally observed by Tsai et al. (Fig. 2 in Ref. [2]). 4.5.4. Iso-surfaces of critical disturbance vorticity and pressure Figure 4.9 Surfaces of constant (a) azimuthal vorticity and (b) pressure of the critical disturbance at 𝑅𝑒𝑜 = 324 and 𝑚 = 4. Figure 4.9 depicts surfaces of constant (a) azimuthal vorticity and (b) pressure for the critical disturbance (the first column in Table 4.3). The light (dark) color corresponds to the water (oil) domain and 0.2 (0.02) values of maximal (a) vorticity and (b) pressure values. The surfaces visualize the four-branch patterns of the 𝑚 = 4 disturbance fields. 4.6. Concluding remarks Our numerical study indicates that the striking shapes of the interface, observed by Fujimoto & Takeda [3], develop in unsteady three-dimensional flow states. The instability develops at 𝑅𝑒𝑜 = 324 at which the interface deformation is small and has the “hump” shape typical for a slow motion. This instability likely is of the shear-layer type, corresponds to time-oscillating disturbances with the azimuthal wave number 𝑚 = 4, and emerges in the lower-fluid (water) domain. It is also revealed that even at this rather small 𝑅𝑒𝑜, the water flow has a non-trivial topology, consisting of the bulk meridional circulation of the same direction, as the upper-fluid (oil) flow has, and a thin counter-circulation region adjacent to the interface. It is shown how this topology develops as 𝑅𝑒𝑜 increases starting from the creeping (Reo << 1) Metamorphosis in a water-oil swirling flow 49 flow. The creeping flow is also multicellular and its topology depends of fluid fractions. The physical reason for the two-cell bulk water flow is in the competing driving factors: (i) the oil meridional circulation which tends to induce the counter-circulation in water and (ii) the oil swirl which tends to induce the co-circulation in water. For small water fraction (𝐻𝑤) and slow rotation (𝑅𝑒), the meridional driving (i) dominates. As 𝐻𝑤 and 𝑅𝑒 increases, the swirl driving (ii) intensifies, causing the water co-circulation. It emerges near the bottom center, then expands and occupies nearly the entire water domain except a thin layer of counter-circulation adjacent to the interface. The layer 𝑧-thickness near the axis diminishes to zero at the axis. The water up-flow, resulting from the co-circulation, contributes to the significant rise of the interface near the container axis observed in the experiment [3]. 50 5. PATTERNS AND STABILITY OF A WHIRLPOOL FLOW his numerical study reveals stable multi-eddy patterns of a steady axisymmetric air-water flow driven by the rotating bottom disk in a vertical sealed cylindrical container. As rotation strength Re increases, eddies emerge, coalesce, separate, and disappear in both air and water. The topological scenario varies with water volume fraction 𝐻𝑤 according to the results obtained for 𝐻𝑤 = 0.3, 0.5, and 0.8. Interesting features are (a) zipper-like chains of air and water eddies forming near the bent interface and (b) bubble-ring air eddies existing in the Re ranges specified in the paper. The stability analysis, performed with the help of a novel efficient technique for two-fluid flows, shows that these multi-eddy motions are stable. The shear-layer instability develops as the interface approaches either the top or bottom of container and some eddies vanish. The physical reasoning behind the eddy formation and the flow instability is provided. The results are of fundamental interest and can have applications in bioreactors. 5.1. Introduction An intriguing and important fluid-mechanics phenomenon is the emergence of a local circulation cell in a swirling flow, often referred to as vortex breakdown (VB). VB applications include deltawing aircraft, where VB is dangerous causing an abrupt change in lift and drag, combustion chambers, where VB is beneficial stabilizing flame, and natural swirling jets like tornadoes, where VB decreases the twister strength. Escudier[108] performed a comprehensive review of early VB studies. More recent works, including VB control strategies, are discussed in Shtern [109]. Vogel [136] and Escudier[121] initiated fundamental VB studies in a sealed cylindrical container with one end disk rotating. An advantage is the closed domain with well-defined and controlled boundary conditions allowing for meaningful comparisons of experimental and numerical results. They well agree as was first shown by Lopez [137]. The analysis of the Vogel-Escudier flow helps understand the VB nature. A recent view is that VB develops via the swirl-decay mechanism (Shtern et al. [109]). While one-fluid VB flows have been studied rather in detail, two-fluid VB flows have not attracted much attention until recent time. The situation changed with the development of aerial bioreactors where air-water flows are used for the growth of tissue culture [138]. The air flow transports the oxygen, required for tissue growth, to the interface and the water circulation enhances mixing of the dissolved oxygen with other ingredients. The tissue fraction is small compared with that of water and is neglected in the studies of flow patterns. The bioreactor applications stimulated the experimental investigations by Lo Jacono et al. [115]and the numerical simulations by Liow et al. [51], [128]. Early numerical studies modeled the gas-liquid interface as a symmetry plane [139] and as a deformable stress-free surface [140]. The first work, which is free from these idealizations of the interface, was performed by Brady et al. [116], [117]. Two-fluid VB flows have a number of unusual features absent in single-fluid flows. One striking feature is the existence of a thin circulation layer (TCL) adjacent to the interface. A TCL attached to the entire interface develops in a water-spout flow [118]. Another striking feature is the emergence of an off-axis VB ring in the depth of a lower fluid [141]. Since eddies arise in both fluids, their T Patterns and stability of a whirlpool flow 51 variety is rich and their transformations are numerous. For example, eighteen topological metamorphoses follow each other as the water volume fraction increases in a truncated conical container where a creeping air-water flow is induced by the slowly rotating top disk [129]. Multiple eddies in both upper and lower fluids in a whirlpool flow driven by the rotation of bottom disk in a sealed vertical cylindrical container have not been studied. Our study partially fills this gap describing and explaining the development of VB cells, TCL, and other eddies. We reveal two more interesting features: (a) zipper-like chains of near-interface eddies and (b) robust bubblerings (RBRs). A bubble-ring typically is a transitional state between bubble and ring shapes. A variation of a control parameter transforms a bubble-ring into either a bubble or a ring. A RBR first was found in a creeping water-spout flow [130]. There RBR occurs in some range of water fraction. The current study shows that RBRs also develop as the rotation strength Re increases in the whirlpool flow. These RBRs exist in wide ranges of Re, listed in the study. The resulting rather complicated topology raises a question about the flow stability. The stability of one-fluid VB flows in a cylindrical container has been investigated rather in detail. The numerical study of Gelfgat et al. [119], [120] showed that the Vogel-Escudier flow can become unstable at either smaller or larger Re than that, at which VB emerges, depending on the length-to-radius ratio, 𝐻. The experimental and numerical studies by Escudier [121]and Sorensen et al. [122]–[124] documented that, as Re increases, the steady axisymmetric VB bubble first develops for 𝐻 < 3.2. For larger 𝐻, the flow first becomes unstable with respect to 3D time-oscillatory disturbances with 𝑚 = 3 for 3.2 < 𝐻 < 4.3, 𝑚 = 2 for 4.3 < 𝐻 < 5.2, and 𝑚 = 4 for 5.2 < 𝐻 < 5.5; 𝑚 is the azimuthal wave number. Herrada et al. [113] found that this instability is of the shear-layer type developing for 𝐻 > 5.5 as well. Here we consider a two-fluid flow and 𝐻 = 1 because this 𝐻 value is close to that used in bioreactors [138]. A technical difficulty of studying the stability of a two-fluid flow is the linearization of a rather complicated relation describing the balance of normal stresses at the bent interface. Herrada & Montanero [1]elaborated an efficient routine which resolves this problem and facilitates numerical simulations. The routine includes (i) mappings converting the time-dependent upper and lower fluid regions onto fixed squared domains, (ii) a symbolic toolbox to calculate the analytical Jacobians, and (iii) the Chebyshev grid in both radial and axial directions. The Chebyshev grid, which concentrates near the interface from both sides, helps observe small bubbles, bubble-rings, and TCLs. The grid also concentrates near the axis and all walls. This helps resolve small circulation cells, emerging near the axis and walls, by using a moderate number of grid points. Herrada & Montanero [1]proved the method efficiency in their study of liquid-bridge dynamics. Here this numerical technique is modified and applied for the whirlpool problem that helps investigate and understand the base-flow topology and the instability nature. Our study provides physical explanations of different topological scenarios as Re increases at 𝐻𝑤 = 0.3, 0.5, and 0.8; 𝐻𝑤 is the water volume fraction. It is found that the instability emerges in a jetlike flow forming in the water. This jet originates in the Kármán boundary layer near the rotating disk, turns upward along the sidewall and then goes toward the axis near the interface. As the jet meets a VB region, it decelerates, diverges and becomes unstable. 58 Global stability analysis of two-fluid flows and CR4. 5.4.4. Merging of near-axis cells in air As 𝑅𝑒 increases, regions CR4 and CR4* merge at 𝑅𝑒 = 𝑅𝑒𝑑4= 394. This is the fourth change in the flow topology. Figure 5.3(b) illustrates this merger: the 𝑤 > 0 range, located near 𝑧 = 1 and observed in the 𝑅𝑒 = 390 curve, is absent in the 𝑅𝑒 = 400 curve. We denote the merged region again as CR4. Figure 5.6 depicts CR4 at 𝑅𝑒 = 400 (a) and 800 (b). The Bödewadt pumping and the 𝑣𝑡 > 0 water flow at the interface work together, driving the anticlockwise circulation of air in CR4, resulting in its expansion in the radial direction as 𝑅𝑒 increases (Figure 5.6). 5.4.5. Separation of water VB bubble from axis The fifth topological metamorphosis is the separation of CR3 from the axis. The boundary between CR1 and CR3 moves up near the axis as 𝑅𝑒 increases (Fig. 6), reaches the interface at 𝑅𝑒 = 𝑅𝑒𝑑5= 873 and moves away from the axis for larger 𝑅𝑒. This transformation occurs together with the second reversal of the meridional velocity at the interface: 𝑣𝑡 becomes negative for small 𝑟 as curve 𝑅𝑒 = 900 in Fig. 4(b) illustrates. Figure 7(a) depicts the resulting streamline pattern at 𝑅𝑒 = 1000. Figure 5.7 Streamline patterns at 𝑅𝑒 = 1000 (a) and 1800 (b) correspond to the transformation of region CR5 from a bubble (a) into a bubble-ring (b). The scenario of this VB development in water is similar to that occurring in the one-fluid flow with the stress-free fixed upper surface, described and explained by BrØns et al. [143]. The difference is that here the VB development in the air flow occurs as well. This development, in particular, includes the appearance of new circulation region CR5 in air depicted in Figure 5.7(a) and induced by the second reversal of meridional velocity at the interface shown in Figure 5.4(b). Region CR5 emerges as a bubble and then transforms into a bubble-ring, as discussed below in more detail. Patterns and stability of a whirlpool flow 59 5.4.6. Emergence of the robust bubble-ring Figure 5.8 Schematics showing the flow topology in (a) Figure 5.7(a) and (b) in Figure 5.7(b). 𝑅𝑒 𝑧𝑖 𝑧𝑠 𝑧𝑠−𝑧𝑖 𝑅𝑒𝑑5 = 873 0.7716 0.7716 0 950 0.7654 0.7874 0.0220 1000 0.7609 0.7927 0.0318 1100 0.7511 0.7939 0.0428 1200 0.7400 0.7853 0.0454 1300 0.7274 0.7699 0.0424 1500 0.6980 0.7246 0.0266 1700 0.6620 0.6693 0.0073 𝑅𝑒𝑑6 = 1779 0.6459 0.6459 0 1900 0.6180 0.6180 0 2000 0.5923 0.5923 0 2200 0.5344 0.5344 0 𝑅𝑒𝑑7 = 2348 0.482 0.482 0 2400 0.4574 0.4631 0.0057 2500 0.4088 0.4256 0.0168 Table 5.1 Dependence on 𝑅𝑒 of 𝑧𝑖, 𝑧𝑠, and CR5 thickness on axis (𝑧𝑠−𝑧𝑖) reveals that CR5 is a bubble-ring for 𝑅𝑒𝑑6 ≤ 𝑅𝑒 ≤ 𝑅𝑒𝑑7. Though Figure 5.7(a) and Figure 5.7(b) look similar, the important difference is that the flow pattern, depicted in Figure 5.7(a), corresponds to the schematic shown in Figure 5.8(a) while the flow, depicted in Figure 5.7(b), corresponds to the schematic shown in Figure 5.8(b), as it follows from Table 5.1. Accordingly, CR5 is a bubble in Figure 5.7(a) and a bubble-ring in Figure 5.7(b). Table 5.1 lists 𝑧-values on the axis for (a) the interface (𝑧𝑖), (b) the boundary, separating regions CR4 and CR5 (𝑧𝑠), see Figure 5.7(a), and (c) thickness of CR5 at the axis (𝑧𝑠−𝑧𝑖) for 𝑅𝑒 shown in the first column. It is clear from Table 1 that region CR5 is a bubble in the range, 873 = 𝑅𝑒𝑑5 < 𝑅𝑒 < 𝑅𝑒𝑑6 = 1779, and is a bubble-ring in the range, 𝑅𝑒𝑑6 ≤ 𝑅𝑒 ≤ 𝑅𝑒𝑑7 = 2348, because its thickness (𝑧𝑠−𝑧𝑖) is zero at the axis. The transformation of CR5 from a bubble into a bubble-ring is the sixth (b) z zs zi r r z zi ,zs (a) 60 Global stability analysis of two-fluid flows topological change. Figure 5.9 Distribution of velocity on the axis at 𝑅𝑒 = 2200 (solid curve) 2300 (dashes) and 2400 (dots). Figure 5.9 confirms the existence of robust bubble-ring by depicting the profile of velocity 𝑤 at the axis for 𝑅𝑒 = 2200 (solid curve) and 𝑅𝑒 = 2300 (dashed curve). The velocity 𝑤 is negative, i.e., directed downward both in water, 0 < 𝑧 < 𝑧𝑖, and in air, 𝑧𝑖 < 𝑧 < 1; being zero only at the bottom, 𝑧 = 0, interface, 𝑧 = 𝑧𝑖, and top, 𝑧 = 1. This corresponds to the topological pattern shown in Figure 5.8(b). As increasing 𝑅𝑒 passes 𝑅𝑒𝑑7= 2348, the thin 𝑤 > 0 region emerges near the interface in the air domain as the dotted curve depicts at 𝑅𝑒 = 2400 in Figure 5.9. We use 𝑤1/3 to conveniently observe the 𝑤 > 0 region in Figure 5.9 where 𝑤 is very small. Therefore, CR5 again becomes a bubble corresponding to the schematic shown in Figure 5.8(a). 5.4.7. Reversal of topological transformations The disappearance of the robust bubble-ring at 𝑅𝑒 = 𝑅𝑒𝑑7 = 2348 is the seventh event and the first back step among the changes in the flow topology. The eighth change is the coalescence of regions CR2 and CR5, shown in Figure 5.10(a), into new region CR2, which extends up to the axis as Figure 5.10(b) depicts. The transition from the pattern in Figure 5.10(a) to that in Figure 5.10(b) has a number of intermediate changes in the flow topology which were systematically described by Brøns [144]. We refer a reader to this instructive paper for details. The ninth change is the splitting of CR4 into CR4a and CR4b as Figure 5.10(c) illustrates at Re = 2500. The tenth change is the disappearance of CR4b; see Figure 5.10(d). -0.25 -0.2 -0.15 -0.1 -0.05 0 0.05 0 0.2 0.4 0.6 0.8 1 w1/3 z Patterns and stability of a whirlpool flow 61 Figure 5.10 Streamline patterns at 𝑅𝑒 = 2200 (a), 2400 (b), 2500 (c) and 2600 (b) show reduction of VB regions. 5.4.8. Instability nature Table 5.2 summarizes the instability results for the flow at 𝐻𝑤 = 0.8. It shows parameter values corresponding to the marginal 𝑅𝑒 (minimal at a fixed 𝑚 and 𝜔𝑖 = 0) and the critical 𝑅𝑒 (shown bold and minimal among marginal). The critical 𝑅𝑒 corresponds to single-helix disturbances, i.e., 𝑚 = 1. No instability was found for 𝑚 = 0. Disturbances with other 𝑚 are less dangerous than those presented in Table 2. Figure 5.11 depicts the streamline pattern at critical 𝑅𝑒 = 2629 and contours of critical-disturbance squared velocity modulus, 𝐸𝑑 = constant; 𝐸𝑑 is averaged with respect to time and azimuthal angle and normalized by its maximal value, 𝐸𝑑𝑚. 𝑅𝑒 2629 2664 2715 2667 𝑚 1 2 3 4 𝜔𝑟 0.204 0.615 1.054 1.376 Table 5.2 Critical (bold) and marginal parameters of the flow instability at 𝐻𝑤 = 0.8. 62 Global stability analysis of two-fluid flows Figure 5.11 Streamline pattern (a) and disturbance energy contours (b) at critical 𝑅𝑒 = 2629 for 𝐻𝑤= 0.8. The 𝐸𝑑 peak is located in the water domain at 𝑟 = 𝑟𝑚 = 0.298 and 𝑧 = 𝑧𝑚 = 0.437. 𝐸𝑑 is localized near the peak: the outmost 𝐸𝑑 contour corresponds to 𝐸𝑑𝑚/10 in Figure 5.11(b). These features help understand the instability nature. The flow pattern in Figure 5.11(a) visualizes the jet development near the rotating bottom disk where streamlines are packed. This jet forms in the Kármán boundary layer, turns upward and goes near the sidewall, reaches the interface and turns toward the axis along the interface. The convergence to the axis accelerates the jet that compensates the entrainment-caused widening thus making the jet width nearly constant as Figure 5.11(a) shows near the interface middle. The situation changes as the jet approaches the VB bubble. The jet turns down and away from the axis and diverges where the disturbance energy peaks. We guess that these jet deceleration and divergence near the VB region cause the shear-layer instability typical of jet flows [109]. Figure 5.12 Profiles of the disturbance energy (solid curve) and base-flow velocity (dashed curve) at 𝑧 = 𝑧𝑚. Both energy and velocity are normalized by their maximum values. Figure 5.12 supports this conjecture by depicting the 𝑟-distribution of the disturbance energy -1 -0.5 0 0.5 1 0 0.2 0.4 0.6 0.8 1 Ed,  w r Patterns and stability of a whirlpool flow 63 𝐸𝑑 (solid curve) and the base-flow 𝑤-velocity (dashed curve) at 𝑧 = 𝑧𝑚. Since the jet is directed downward near the interface, we depict 𝑤, normalized by its maximal magnitude, for convenient discussion of Figure 5.12. The 𝑤-profile indeed has a jet-like shape near 𝑟 = 0.4 with inflection points located near 𝑟 = 0.3 and 𝑟 = 0.5 in Figure 5.12. The disturbance energy peaks near the 𝑟 = 0.3 inflection point and has a local maximum near the 𝑟 = 0.5 inflection point. These features indicate that the instability likely is of the shear-layer type. The multi-helix instability with 𝑚 = 2 and 4 develop for higher 𝑅𝑒 than that with 𝑚 = 1 according to Table 5.2, but the marginal 𝑅𝑒 values are close to the critical 𝑅𝑒. This can lead to polygonal patterns observed in experiments ([145], [146], see also references there). For this reason, Figure 5.13 depicts the Ed distribution for the marginal disturbances with 𝑚 = 2 and 4 presented in Table 5.2. The energy of these disturbances is less localized (especially at 𝑚 = 2) than that for the 𝑚 = 1 disturbance as comparison of Figure 5.11(b) and Figure 5.13 illustrates. Nevertheless, the marginal disturbance also located in the water domain, where the jet “reflects” from the VB region in Figure 5.11(a). To summarize, the stability results reveal that the topological metamorphoses occur in the stable steady axisymmetric flow. The instability develops as the interface becomes strongly deformed and approaches the top disk of container (Figure 5.11(a)). The analysis of the base-flow pattern and the energy distribution at the critical 𝑅𝑒 indicates that the instability likely is of the shearlayer type. Figure 5.13 Energy contours for marginal disturbance with 𝑚 = 2 (a) and 4 (b) at 𝐻𝑤 = 0.8. 5.5. Moderately deep whirlpool Now we describe changes in the flow topology occurring as 𝑅𝑒 increases at the equal volumes of water and air, 𝐻𝑤 = 0.5, and examine the flow stability. 5.5.1. Flow reversal near the interface Figure 5.14 depicts the streamlines patterns at 𝑅𝑒 = 250 (a) and 300 (b). The pattern shown in Figure 5.14(a) is similar to that shown in Figure 5.2(a). The difference is that the outmost Moffatt eddy ME is observed in Figure 5.14while it is not seen in Figure 5.2 being unresolved by our standard grid. Moffatt [132] revealed that a flow in a corner (between two inclined planes with a sufficiently small angle) has an infinite set of eddies whose strength and scale diminish to zero as the edge is 64 Global stability analysis of two-fluid flows approached. Since the eddy size is small compared with the cylinder radius (see ME in Figure 5.14), the sidewall can be approximated by a tangential plane thus reducing the cylindrical geometry to the corner between orthogonal planes. The corner angle, 90, is sufficiently small for the Moffatt eddies to occur. The air circulation in CR2 generates the Moffatt corner eddies. Among them, only the outmost eddy is resolved by our standard grid and denoted as ME in the upper-right corner of the air flow in Figure 5.14. Figure 5.14 Streamline patterns at 𝐻𝑤 = 0.5 for 𝑅𝑒 = 250 (a) and 300 (b) show the outmost Moffatt eddy ME and the emergence of circulation regions CR3 and CR4. Figure 5.15 Distributions of velocity 𝑤 at the axis, 𝑟 = 0, for 𝑅𝑒 = 250 (dashed curve) and 300 (solid curve) show the emergence of circulation regions CR3 and CR4, see Fig. 14(b), in the water (air) region 0 < 𝑧 < 𝑧𝑖 (𝑧𝑖 < 𝑧 < 1); 𝑧𝑖 = 0.497 (0.495) at 𝑅𝑒 = 250 (300). The Moffatt eddies can also develop in a flow inside a cone if cone half-angle c < 80.9 [147]. A cone becomes a plane at c = 90. This angle is too large for the Moffatt eddies to occur. Accordingly, no eddy exists near the top disk and axis intersection. -0.01 -0.005 0 0.005 0 0.25 0.5 0.75 1 w z CR3 CR4 300 250 Patterns and stability of a whirlpool flow 65 Figure 5.16 Profiles of meridional velocity on the interface 𝑣𝑡 at 𝑅𝑒 = 250, 300, and 400 illustrate the emergence and radial expansion of regions CR3 and CR4 (where 𝑣𝑡 > 0) as 𝑅𝑒 increases. Figure 5.14(a) remains topologically invariant as 𝑅𝑒 decreases down to zero. In contrast as 𝑅𝑒 increases, a local flow reversal occurs near the axis-interface intersection point, 𝑟 = 0 and 𝑧 = 𝑧𝑖 = 0.497, and a two-fluid VB bubble, CR3 and CR4, emerges as Figure 5.14(b) shows. The topological scenarios differ: the VB occurs in the water depth for 𝐻𝑤 = 0.8, but at the interface for 𝐻𝑤 = 0.5. Figure 5.15 depicts the distribution on the axis of velocity 𝑤 for 𝑅𝑒 values shown near the curves. At 𝑅𝑒 = 250, the velocity is negative in the range, 0 < 𝑧 < 𝑧𝑖, that corresponds to water region CR1, and positive in the range, 𝑧𝑖 < 𝑧 < 1 that corresponds to air region CR2 in Figure 5.14(a). At 𝑅𝑒 = 300, 𝑤 changes its signs near the interface in both water, 0.33 < 𝑧 < 𝑧𝑖 = 0.495, and air, 𝑧𝑖 < 𝑧 < 0.551, flows. These 𝑧-ranges respectively correspond to CR3 and CR4 shown in Figure 5.14(b). Interpolating yields that CR3 and CR4 simultaneously emerge at 𝑅𝑒 = 𝑅𝑒1𝑚= 257. The leading event is the appearance of region CR3 in water which is a corner bubble. The emergence of region CR4 in air is a sequence of 𝑣𝑡 reversal near the axis which accompanies the CR3 development. Figure 5.16 depicts 𝑣𝑡1/3 to better observe small positive values of 𝑣𝑡 near the axis at 𝑅𝑒 = 300 (curve 2) and 400 (curve 3). Velocity 𝑣𝑡 becomes positive for 𝑅𝑒 > 257 and the 𝑣𝑡 > 0 range, corresponding to CR3 and CR4 in Figure 5.14(b), expands as 𝑅𝑒 increases (Figure 5.16). The emergence of the two-fluid VB bubble is the first change in the flow topology at 𝐻𝑤 = 0.5. 5.5.2. Flow reversal near the upper wall Figure 5.17 Streamline patterns at (a) 𝑅𝑒 = 350 and (b) 400 shows (a) the emergence of circulation region CR4* and (b) merger of CR4 and CR4*. -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0 0.2 0.4 0.6 0.8 1 vt1/3 r 250 300 400 66 Global stability analysis of two-fluid flows Figure 5.18 Distributions of velocity 𝑤 at the axis at 𝑅𝑒 = 300, 350, and 400 illustrate the emergence of circulation region CR4* (curves 300 and 350) and merging of CR4 and CR4* (curves 350 and 400, see also Fig. 17). The second change in the flow topology as 𝑅𝑒 increases is the development of local countercirculation of air near the center of the top disk at 𝑅𝑒 = 𝑅𝑒2𝑚= 347. Figure 5.17(a) depicts this region (CR4*) at 𝑅𝑒 = 350. The physical reasoning for CR4* emergence is the Bödewadt pumping similar to that in the 𝐻𝑤= 0.8 case. Figure 5.18 depicts the distribution of velocity 𝑤 at the axis at 𝑅𝑒 = 300, 350 and 400. We plot 𝑤1/3 to better observe circulation region CR4 (in the range 0.495 < 𝑧 < 0.632 for curve 350) and the emergence of circulation region CR4* (in the range 0.947 < 𝑧 < 1 for curve 350) where |𝑤| is small. Comparison of curves 300 and 350 in Figure 5.18 helps understand how region CR4* emerges. 5.5.3. Merging of anticlockwise circultation regions in air flow The third change in the flow topology as Re increases is the merging of regions CR4 and CR4* at 𝑅𝑒 = 𝑅𝑒3𝑚 = 393. The two saddle stagnation points in Figure 5.17(a), located on the axis between regions (i) CR2 and CR4 and (ii) CR2 and CR4*, coalesce at 𝑅𝑒 = 393 in one saddle point S, which separates from the axis for 𝑅𝑒 > 393 as Figure 5.17(b) shows at 𝑅𝑒 = 400. The merged region, denoted again as CR4, extends from the interface up to the top disk for 𝑅𝑒 > 393. Comparison of Figure 5.17(a) and Figure 5.17(b) shows how the streamline pattern changes. Figure 5.18 depicts the corresponding change in the velocity distribution at the axis: ranges CR4 and CR4*, where 𝑤 < 0 in the air flow (curve 350), merge into the range, 𝑧𝑖 < 𝑧 < 1 (curve 400), as 𝑅𝑒 increases. The anticlockwise circulation in the merged region CR4 in Figure 5.19(a) is driven by two factors: (i) positive meridional velocity at the interface, 𝑣𝑡 > 0, and (ii) Bödewadt pumping. In contrast to the opposite effects of 𝑣𝑡 and swirl in region CR2 for 𝑅𝑒 < 325, factors (i) and (ii) cooperate in region CR4 for 𝑅𝑒 > 303 and result in the radial expansion of CR4 as Figure 5.19 illustrates. Saddle S and center C stagnation points of the meridional flow in Figure 5.19(a) merge and disappear in Figure 5.19(b) via the saddle-center bifurcation [144]. -0.3 -0.15 0 0.15 0 0.25 0.5 0.75 1 w1/3 z 300 400 350 CR4 CR4* Patterns and stability of a whirlpool flow 67 Figure 5.19 Streamline patterns at (a) 𝑅𝑒 = 450 and (b) 1000 show the separation of region CR3 from axis and the emergence of region CR5. 5.5.4. Separation of water VB region from the axis The forth change in the flow topology is the separation of water VB region CR3 from the axis at 𝑅𝑒 = 𝑅𝑒4𝑚 = 434. Comparison of Figure 5.17(a) and Figure 5.17(b) shows that the boundary between regions CR1 and CR3 moves upward near the axis as 𝑅𝑒 increases. The boundary reaches the interface at 𝑅𝑒 = 434 and the separation occurs. Figure 5.19(a) depicts the streamline pattern and separated region CR3 at 𝑅𝑒 = 450. Figure 5.20 Distribution of meridional velocity on the interface 𝑣𝑡 at 𝑅𝑒 = 400 (dotted curve), 450 (solid curve), and 1000 (dashed curve) shows the emergence (solid curve) and expansion (dashed curve) of region CR5 shown in Fig. 19(b). Figure 5.20 depicts the distribution of meridional velocity at the interface 𝑣𝑡 at 𝑅𝑒 = 400 (dotted curve), 450 (solid curve), and 1000 (dashed curve), and illustrates one more reversal of 𝑣𝑡 which becomes negative near the axis for 𝑅𝑒 > 434. The 𝑣𝑡 > 0 range of curve 2 in Figure 5.20 corresponds to region CR3 in Figure 5.19(a) and the left 𝑣𝑡 < 0 range verifies that region CR3 is indeed separated from the axis at 𝑅𝑒 = 450. This near-axis 𝑣𝑡 < 0 range corresponds to new region CR5 which radially expands as 𝑅𝑒 -0.6 -0.4 -0.2 0 0.2 0 0.2 0.4 0.6 0.8 1 vt1/3 r CR3 CR5 74 Global stability analysis of two-fluid flows 𝑧𝑠−𝑧𝑖 becomes (does not become) negative as 𝑅𝑒 increases according to Table 5.3 (Table 5.5). Physical reasoning behind this difference is that an air (water) flow cannot (can) easily reverse a water (air) flow due to the small density ratio, 𝜌𝑟 = 0.00122. For 𝑅𝑒 slightly exceeding 690, the air flow in region CR3, shown in Figure 5.27(a), is not capable to reverse the interface meridional motion, which is driven by the water flow. For this reason, region CR3 touches the interface only at a single point, namely at the axis, where 𝑣𝑡 = 0. The Bödewadt pumping extends region CR3 from the axis along the interface in Figure 5.27(a), but cannot overcome the converging motion at the interface, leaving a thin layer of CR2 region, where the air meridional velocity reverses from being directed toward the axis at the interface into being directed away from the axis at the CR3 lower boundary. To mathematically investigate, why region CR2 is thin between CR1 and CR3 in Figure 5.27(a), consider the Stokes stream function, , near the axis-interface intersection. With no loss of generality, we can apply the representation,  = r2[zF(r)]1(r,z), 5.9 where 1 is a bounded function. Presentation (5.9) takes into account the boundary conditions at the axis, 𝑟 = 0, and the interface, 𝑧 = 𝐹(𝑟). Since the axis is a streamline, we can put  = 0 at 𝑟 = 0. The regularity requirement yields that /r = 0 at 𝑟 = 0 and, therefore,  is proportional to 𝑟2. The interface, 𝑧 = 𝐹(𝑟), is a stream surface, where  = 0, that explains the presence of expression in the rectangular brackets in (8). The boundary between regions CR2 and CR3, e. g., in Figure 5.27(a), z = Fs(r), is also a stream surface where  = 0. Therefore for 𝑅𝑒 > 𝑅𝑒𝑠1 = 174, we can use the presentation, 1(r,z) = [zFs(r)]2(r,z), where 2 is a bounded function. Due to the axial symmetry, both F(r) and Fs(r) are even functions of r and their power-law expansions are F(r) = zi+air2+O(r4) for the interface and Fs(r) = zs+asr2+O(r4) for separatrix near the axis, where r << 1. Here subscripts “i” and “s” denote the interface and separatrix, respectively. Both the interface and separatrix are normal to the axis at the intersection points, 𝑟 = 0, 𝑧 = 𝑧𝑖 and 𝑧 = 𝑧𝑠, where 𝑤 = 0. This feature remains valid after merging of 𝑧𝑖 and 𝑧𝑠, i.e., the interface and separatrix are tangential at r = 0. This tangency is a mathematical reason why region CR2 is thin near the axis in Figure 5.27(a). Figure 5.28 Profiles of velocity 𝑤 on axis show how 𝑤 > 0 range shrinks as 𝑅𝑒 increases from 200 (solid curve) to 250 (dashed curve). -0.003 -0.0015 0 0.0015 00.2 0.4 0.6 0.8 1 w z 200 250 zi zs Patterns and stability of a whirlpool flow 75 Figure 5.28 and Figure 5.29 depict how the distribution of velocity on the axis, 𝑤(𝑧) at 𝑟 = 0, changes as the roots 𝑧𝑖 and 𝑧𝑠 merge. At 𝑅𝑒 = 200 (the solid curve in Figure 5.28), 𝑤 is positive in the air flow, corresponding to the CR2 circulation, except in the 𝑤 < 0 range, 0.9134 = 𝑧𝑠 < 𝑧 < 1, corresponding to the CR3 circulation in Figure 5.26(b). Figure 5.29 (a) Velocity on axis at 𝑅𝑒 = 700, (b) close up of interface (𝑧𝑖) vicinity. Figure 5.30 (a) Distributions of swirl 𝑣 and meridional 𝑣𝑡 velocities at the interface at 𝑅𝑒 = 800 (dotted curves), 900 (solid curves), and 1100 (dashed curves); (b) distributions of 𝑣𝑡1/3 reveal the appearance of 𝑣𝑡 > 0 range, CR4 and CR5 as Re increases. The boundary, separating CR2 and CR3 (referred to as “separatrix” here), intersects the axis at 𝑧 = 𝑧𝑠. At Re = 250 (dashed curve in Figure 5.28), the 𝑤 > 0 range, 𝑧𝑖 < 𝑧 < 𝑧𝑠 (corresponding to CR2), is small compared with the w < 0 range, 𝑧𝑠 < 𝑧 < 1 (corresponding to CR3). The 𝑤 > 0 range collapses at 𝑅𝑒 = 𝑅𝑒𝑠2 = 690. Figure 5.29 depicts the velocity distribution on the axis at 𝑅𝑒 = 700 where 𝑧𝑖 and 𝑧𝑠 are merged (Table 5.5). The dashed (solid) curves correspond to the water (air) flow in Figure 5.29. Curve w(z) looks tangential to line 𝑤 = 0 at 𝑧 = 𝑧𝑖 in Figure 5.29(a). However, the close up of the interface vicinity, depicted in Figure 5.29(b), reveals that the curve and the line are indeed not tangential. Despite the curve does touch the line at 𝑧 = 𝑧𝑖, the derivative, w/z, is not zero and has a jump at 𝑧 = 𝑧𝑖. -0.025 -0.02 -0.015 -0.01 -0.005 0 0 0.2 0.4 0.6 0.8 1 w z (a) zi -0.0003 -0.0002 -0.0001 0 0.0001 0.2 0.22 0.24 0.26 w z (b) zi -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0 0.2 0.4 0.6 0.8 1 v,vt vt v r (a) -0.6 -0.4 -0.2 0 0.2 0 0.2 0.4 0.6 0.8 1 (b) vt1/3 CR4 CR2 CR5 r 76 Global stability analysis of two-fluid flows The jump follows from the normal stress balance (equation 5.6) and the different properties of air and water. A jump in w/z at the interface is a specific feature occurring in a two-fluid flow where the fluid density and viscosity are not continuous at the interface. The presence of region CR2 (or other flow cells described below) is necessary in order to separate CR1 and CR3 which have the same circulation direction and therefore can only touch each other at a saddle point on the interface. One such point is located at the axis in Figure 5.27(a). Other saddle points, located away from the axis for 𝑅𝑒  𝑅𝑒𝑠3 = 829, are discussed below. 5.6.3. Off-axis vortex breakdown in water flow The third topological transformation is the VB emergence in water away from the axis at 𝑅𝑒 = 𝑅𝑒𝑠3 = 829. This topological change occurs due to the air flow which locally reverses the meridional velocity at the interface, 𝑣𝑡. Figure 5.30 depicts the distributions of swirl 𝑣 and meridional 𝑣𝑡 velocities at the interface at 𝑅𝑒 = 800 (dotted curves), 900 (solid curves), and 1100 (dashed curves). For 𝑅𝑒 < 829, the meridional velocity is directed toward the axis: 𝑣𝑡 < 0 in 0 < 𝑟 < 1, as the dotted curves in Figure 5.30 show. For 𝑅𝑒 > 829, a range of 𝑣𝑡 > 0 exists. We reveal this by plotting 𝑣𝑡1/3 in Figure 5.30(b) in order to better observe regions CR4 and CR5, located in the range 0 < 𝑟 < 0.4, where 𝑣𝑡 is very small in Fig. 30(a). Figure 5.31 Schematics show the flow topology at 𝑅𝑒 = 829 (a) and 𝑅𝑒 > 829 (b). The bold lines denote the interface and the arrows indicate the flow direction. Figure 5.30(a) shows that the distribution of swirl velocity 𝑣 is nearly linear (solid-body-rotation) and meridional velocity 𝑣𝑡 is negligibly small compared with 𝑣 in the range, 0 < 𝑟 < 0.4. Therefore, this part of the interface acts like a rotating disk for adjacent air. The Kármán pumping pushes air away from the axis near the interface. The air meridional velocity is maximal near 𝑟 = 0.35, where air-flow streamlines converge and are maximally close to the interface in region CR3 (Figure 5.27). This explains why the air flow reverses 𝑣𝑡 and a small portion of the water flow near 𝑟 = 0.35 as 𝑅𝑒 increases. Figure 5.30(b) shows that there is a local maximum of 𝑣𝑡 near 𝑟 = 0.35, whose value is negative at 𝑅𝑒 = 800 and positive at 𝑅𝑒 = 900. Interpolation yields that the maximum value is zero at 𝑅𝑒 = 𝑅𝑒𝑠3 = 829. At this 𝑅𝑒, the boundary between regions CR2 and CR3 touches the interface at point SP as Figure 5.31(a) schematically shows. SP is a saddle point where the meridional motion stagnates, 𝑣𝑡 = 0, but 𝑣  0. CR2 becomes divided into two regions: (i) near-sidewall part (denoted again as CR2) which has the ring shape, as Figure 5.27(b) shows, and (ii) near-axis part denoted as CR2a in Figure 5.31(a) which inherits the bubble-ring shape. At 𝑅𝑒 > 829, SP splits in two saddle points SP1 and SP2, located on the interface as Figure 5.31(b) shows. The interval between SP1 and SP2, where the flow goes away from the axis in Figure 5.31(b), corresponds to the 𝑣𝑡 > 0 range in Figure 5.30. The 𝑣𝑡 reversal results in the appearance of two new cells: CR4 in water and CR5 in air existing for Re > 829 as Figure 5.31(b) S P Patterns and stability of a whirlpool flow 77 illustrates. Thus, an off-axis VB occurs resulting in the emergence of VB ring CR4 which expands as 𝑅𝑒 increases being well observed at 𝑅𝑒 = 900 in Figure 5.27(b). Figure 5.32 Distribution of meridional velocity on the interface at 𝑅𝑒 = 1400 and 1450 shows reversal of 𝑣𝑡 near axis and appearance of region CR6. In contrast, region CR5 (former region CR2a) remains thin; it is a near–axis part separated from former region CR2 shown in Figure 5.27(a). Region CR5 is depicted large in Figure 5.31(b) just for convenient observation. Indeed, region CR5 is so thin that it is not visible at 𝑅𝑒 = 900 in Figure 5.27(b). The CR5 existence follows from the above analysis of flow topology for 𝑅𝑒 > 829 and from the 𝑣𝑡 distribution depicted in Figure 5.30(b), where region CR5 corresponds to the near-axis range where 𝑣𝑡 < 0. Thus, the third change in the flow topology transforms CR2 from a bubble-ring into a ring and results in the emergence of two new circulation regions: (i) VB ring CR4 in water and (ii) bubble-ring CR5 in air. 5.6.4. Chain of VB regions The fourth and fifth topology transformations occur as the interface approaches the bottom disk. Figure 5.32 depicts the second reversal of the meridional velocity on the interface, which occurs near the axis at 𝑅𝑒 = 𝑅𝑒𝑠4  1440. This reversal results in the appearance of the 𝑣𝑡 > 0 near-axis range, CR6, in Figure 5.32 and the second VB bubble, CR6, in water. As Re exceeds 𝑅𝑒𝑠5, which is slightly less than 1500, the third reversal of the meridional velocity separates region CR6 from the axis, where the new 𝑣𝑡 < 0 range, CR7, emerges, as Figure 5.33 depicts. Figure 5.34 shows the streamline patterns at (a) 𝑅𝑒 = 1300 (as a reference) and (b) 𝑅𝑒 = 1500 where the new-born VB ring CR6 is observed in the water flow. -0.6 -0.4 -0.2 0 0.2 0 0.2 0.4 0.6 0.8 1 1450 1400 vt1/3 r CR6 CR5 CR4 78 Global stability analysis of two-fluid flows Figure 5.33 Distribution of meridional velocity on the interface at 𝑅𝑒 = 1500 shows the appearance of region CR7. Thin circulation layers CR5 and CR7, located in air, are not visible in Figure 5.34(b). They are much thinner than water VB cells CR4 and CR6 due to the small air-to-water density ratio, 𝜌𝑟 = 0.00122. Being a close-up of Fig. 34(b), Figure 5.34(c) reveals layer CR5, but its resolution is poor despite we use the fine grid here. Region CR7 is too thin to be observed even in Figure 5.34(c); its existence follows from the 𝑣𝑡 distribution depicted in Figure 5.33. Figure 5.34 Streamline patterns at 𝑅𝑒 = 1300 (a) and 1500 (b) show appearance of circulation region CR6; (c) is close-up of (b) visualizing region CR5. Based on the above analysis, we interpret the flow pattern in the left-lower corner of Figure 5.34(b) as the zipper-like near-interface chain of VB regions schematically shown in Figure 5.35. The chain develops as the interface approaches the bottom disk similar to the chain development at 𝐻𝑤 = 0.5 (Figure 5.21 and Figure 5.22). -0.6 -0.4 -0.2 0 0.2 0 0.2 0.4 0.6 0.8 1 vt1/3 r CR2 CR4 CR6 CR7 CR5 Patterns and stability of a whirlpool flow 79 Figure 5.35 Schematic shows the zipper-like topology of shallow whirlpool for 𝑅𝑒  1500 (see Figure 5.33). The bold line is the interface and arrows indicate flow direction. 5.6.5. Reduction of VB cells by interface deformation As the interface deformation becomes remarkable, it reduces vortex breakdown cells in both fluids. Figure 5.34 illustrates this effect: region CR4 shrinks as 𝑅𝑒 increases from 1300 in Figure 5.34(a) to 1500 in Figure 5.34(b). The physical reasoning behind this trend is the swirl decay mechanism (Shtern et al. [114]). The interface deformation enlarges the axial extent of both air and water domains. The enlarged axial extent increases the trajectory length and travel time of fluid particles transporting the angular momentum from the rotating disk to the interface in the water flow and from the interface to the axis in the air flow. This enhances the momentum losses due to viscous diffusion and friction, thus weakening the mechanism of vortex breakdown development. 5.6.6. Instability of Hw = 0.3 flow Similarly to the 𝐻𝑤 = 0.5 case, the 𝐻𝑤 = 0.3 flow has the critical instability with respect to singlehelix (𝑚 = 1) disturbances and the marginal instability with respect to two-helix (𝑚 = 2) disturbances, according to Table 5.6. For larger Re, the interface reaches the bottom near the axis that makes our numerical technique inapplicable. Figure 5.36 shows (a) base-flow streamlines and (b) energy contours for the critical disturbance at critical 𝑅𝑒 = 1547. The energy peak is located at 𝑟 = 0.383 and 𝑧 = 0.1295 that is very close to the interface where the jet-like water flow, converging toward the axis along the interface, meets the VB bubble (Figure 5.36). This feature is similar to that at 𝐻𝑤 = 0.8 and 0.5. 𝑅𝑒 1547 1583 𝑚 1 2 𝜔𝑟 0.673 1.297 Table 5.6 Critical (bold) and marginal parameters of the flow instability at 𝐻𝑤 = 0.3. 80 Global stability analysis of two-fluid flows Figure 5.36 (a) Base-flow streamlines and (b) energy contours of critical disturbance at 𝐻𝑤 = 0.3 and Re = 1547. 5.7. Conclusions 5.7.1. Robust bubble-rings An interesting finding of this paper is that robust bubble-ring circulation cells (RBRs) emerge in a whirlpool flow as the rotation speeds up. In general, a bubble-ring is a transitional shape. A variation of a control parameter (e.g., 𝑅𝑒) typically transforms a bubble-ring into a bubble or into a ring. In contrast in the whirlpool flow, we have revealed a few RBRs, existing in the ranges of 𝑅𝑒 listed below. The bubble-rings are robust here because of (a) the physical limitation that the interface cannot separate from the axis away from the end disks and (b) an air flow cannot easily reverse a water flow due to small air-to-water density ratio. The RBRs are found in the deep (𝐻𝑤 = 0.8) and shallow (𝐻𝑤 = 0.3) whirlpools. In the 𝐻𝑤 = 0.8 whirlpool, bubble-ring CR5, schematically shown in Figure 5.8(a), exists for 1779 < 𝑅𝑒 < 2348 (Table 5.1). In the 𝐻𝑤 = 0.3 whirlpool, the first RBR, CR2 in Figure 5.27(a), is observed in the range, 690 < 𝑅𝑒 < 829 (Table 5.5). The second RBR (CR5) emerges at 𝑅𝑒 = 829 being separated from CR2 by a local reversal of meridional velocity on the axis (Figure 5.31 and Figure 5.32). The third RBR (CR7) emerges near the interface-axis intersection at 𝑅𝑒 = 𝑅𝑒𝑠5  1500 (Figure 5.33 and Figure 5.35). All these RBRs (CR2, CR5, and CR7) occur in the air flow and are adjacent to the interface. Figure 5.37 (a) Global CR2 and (b) local CR5 kinds of bubble-ring cells. The bold curves depict the interface and the bold lines depict the sidewall. Patterns and stability of a whirlpool flow 81 Two kinds of RBRs are observed: (a) global, extending from the axis up to the sidewall, CR2 in Figure 5.27(a) and Figure 5.37(a), and (b) local, bounded by the interface and separatrix, CR5 in Figure 5.37(b). The bold curves (vertical lines) denote the interface (sidewall) and the thin curves denote the boundaries between CR2 and CR3 in Figure 5.37(a), and CR4 and CR5 in Figure 5.37(b). The local RBR is shown large for convenient observation in Figure 5.37(b), but the radial extent of a local RBR can be small, e.g., 0 < 𝑟 < 0.011 for region CR7 at 𝑅𝑒 = 1500 (Figure 5.33 and Figure 5.35). 5.7.2. Multiple changes in flow topology N Re Event Figure 1 359 Emergence of VB bubble CR3 deep in water Figure 5.2 2 365 Emergence of two-fluid VB bubble CR3-CR4 Figure 5.5(a) 3 377 Emergence of VB bubble CR4* in air Figure 5.5(b) 4 394 Merger of CR4 and CR4* in air Figure 5.6(a) 5 873 Separation of CR3 from the axis Figure 5.7(a) 6 1779 Emergence of RBR CR5 in air Figure 5.7(b) and Figure 5.8 7 2348 Transformation of RBR CR5 into a bubble Figure 5.9 8  2400 Merger of CR2 and CR5 Figure 5.10(b) 9  2450 Splitting of CR4 into CR4a and CR4b Figure 5.10(c) 10  2600 Disappearance of CR4b Figure 5.10(d) 2629 Instability development Figure 5.11 Table 5.7 Changes in flow topology as 𝑅𝑒 increases at 𝐻𝑤 = 0.8. The second interesting finding is a rich variety of topological metamorphoses in the whirlpool flow. Topological scenarios depend on the water volume fraction 𝐻𝑤. Vortex breakdown (VB) first emerges in (i) the air flow near the top disk center at 𝐻𝑤 = 0.3, (ii) at the interface-axis intersection at 𝐻𝑤 = 0.5, and (iii) deep in water near the axis at 𝐻𝑤 = 0.8. Table 5.7 Changes in flow topology as 𝑅𝑒 increases at 𝐻𝑤 = 0.8. lists the ten flow topological transformations found as Re increases in the deep whirlpool, 𝐻𝑤 = 0.8. We also found the nine (five) transformations at 𝐻𝑤 = 0.5 (𝐻𝑤 = 0.3); the number of topological changes increases with 𝐻𝑤. The plurality of metamorphoses is due to flow cells emerge, change, and disappear in both fluids. 5.7.3. Chain of VB regions The third interesting finding is the development of VB cell chain adjacent to the interface as the water depth near the axis becomes small. At 𝐻𝑤 = 0.3, there are four VB cells arranged in a zipperlike pattern attached to the interface near the axis (Figure 5.35). At 𝐻𝑤 = 0.5, the chain consists of three VB cells (Figure 5.22). At 𝐻𝑤 = 0.8, the chain consists of two VB cells observed in Figure 5.10(a). 82 Global stability analysis of two-fluid flows 5.7.4. Effect of large deformation of the interface Figure 5.38 Dependence on 𝑅𝑒 of the maximal magnitudes of meridional 𝑣𝑚 and swirl 𝑣𝑠 velocities at the interface. Curve 𝑧𝑖 shows the interface height at the axis. 𝐻𝑤 = 0.5. The forth interesting finding is that the interface deformation (ID), when it becomes remarkable, suppresses VB in both fluids. Figure 5.10 illustrates this effect: VB region CR3 shrinks as ID increases. The physical reasoning behind this feature is the swirl decay mechanism (SDM) (Shtern et al. [114]). SDM plays an important role in the development of vortex breakdown, counterflows, and double counterflows in vortex devices ([109]). According to SDM, a value of Re, at which VB emerges, increases together with the aspect ratio. The aspect-ratio and ID effects are analogous because ID enlarges the axial extent of both lower and upper fluid domains. The enlarged axial extent increases the trajectory length and travel time of fluid particles transporting the angular momentum from the rotating disk to the interface near the sidewall. This enhances the momentum losses due to viscous friction and diffusion, thus reducing the maximal velocity magnitude at the interface and weakening the mechanism of VB development. The largest increase in the axial extent of air and water domains occurs at 𝐻𝑤 = 0.5. Figure 5.38 illustrates the ID effect at 𝐻𝑤 = 0.5 by depicting the dependence on the Reynolds number of the maximal magnitudes of meridional 𝑣𝑚 and swirl 𝑣𝑠 velocities at the interface as well as the interface height at the axis 𝑧𝑖. As 𝑅𝑒 increases from zero, 𝑣𝑚 and 𝑣𝑠 initially grow because the convective transport intensifies while viscous losses remain nearly the same. As 𝑅𝑒 further increases and 𝑧𝑖 becomes remarkably reduced compared with its no-flow value 0.5, the amplified ID makes larger the trajectory lengths and, therefore, viscous losses. For this reason, first 𝑣𝑚 and then 𝑣𝑠 start to decrease. 5.7.5. Instability The instability develops in the water flow because the kinematic viscosity of air is fifteen times that of water and the maximal velocity of air motion (achieved at the interface) is significantly smaller than that of water (achieved at the rotating disk) as Figure 5.38 illustrates. It is striking that the multi-cellular flows, studied here are stable. The instability develops when the number of flow cells decreases as the interface becomes significantly deformed and approaches either the top (at 𝐻𝑤 = 0.8) or bottom (at 𝐻𝑤 = 0.3 and 0.5) disk. The instability occurs with respect to helical disturbances with the azimuthal wave number 𝑚 = 1, 2, 3 and 4. The critical Re values (Table 5.8) correspond to 𝑚 = 1 disturbances. 0 0.1 0.2 0.3 0.4 0.5 0 800 1600 2400 zi vs vm Re vm , vs , zi Patterns and stability of a whirlpool flow 83 Table 5.8 Dependence on water fraction of critical parameters of the whirlpool instability. As Table 5.8 shows, the critical Reynolds number grows and the disturbance frequency drops as the water volume fraction 𝐻𝑤 increases. The instability develops in a jet-like boundary layer. It originates near the rotating disk, propagates along the sidewall and goes to the axis near the interface, e.g., see Figure 5.11(a). The critical disturbance starts to grow where the jet meets VB region, decelerates and diverges. The disturbance energy peaks are located near the jet inflection points (Figure 5.12). This indicates that the instability likely is of shear-layer type. The obtained results seem of fundamental interest and might be useful for aerial bioreactors. 𝐻𝑤 0.3 0.5 0.8 𝑅𝑒 1547 1950 2628 𝜔𝑟 0.673 0.259 0.204 𝑚 1 1 1 90 Global stability analysis of two-fluid flows eddy pattern of the creeping flow transforms into the three-eddy pattern, as the disk rotation speeds up. Figures 9(c) and 9(d) of Ref. [148] show this pattern at 𝑅𝑒 = 10000 and 18000 respectively. Figure 6.3(a) of our study depicts the three-eddy pattern at 𝑅𝑒 = 55700, at which the interface is remarkably deformed upward near the axis and downward near the sidewall, and the flow becomes unstable. Figure 6.3 (a) The base-flow pattern and (b) critical-disturbance energy distribution at 𝐻𝑤 = 0.8. The critical characteristics are found using the interpolation as Figure 6.4 illustrates for the critical 𝑅𝑒. The cross symbols correspond to the numerical results, the curve and formula show the fitting polynomial. The critical value is determined as the polynomial root where 𝜔𝑖 = 0. A similar procedure is used to find the critical value of 𝜔𝑟. The obtained critical characteristics are: 𝑚 = 2, 𝑅𝑒 = 55700, 𝜔𝑟 = 0.0176, 𝜔𝑖 = 0, 𝑅𝑒𝑎 = 3713, and 𝑅𝑒𝑤= 813. In addition to 𝑅𝑒, which is convenient to use in the numerical procedure, we calculate the Reynolds number, 𝑅𝑒𝑎 = 𝑅𝑒/𝜈𝑛, based on the air viscosity and therefore characterizing the strength of the air flow; 𝜈𝑛 = 15 is the air-to-water kinematic viscosity ratio. To characterize the strength of the water flow, we also calculate 𝑅𝑒𝑤 = 𝑣𝑠𝑖𝑅𝑒; 𝑣𝑠𝑖 is the maximal (dimensionless) swirl velocity at the interface. Here 𝑣𝑠𝑖 = 0.0146 is small compared with 1, which is the maximum velocity of the rotating disk. Figure 6.4 Finding the critical 𝑅𝑒 by interpolation. The crosses present the numerical results. The curve and formula show the fitting polynomial. y = 9.753E-18x3 - 2.489E-12x2 + 2.141E-07x - 5.894E-03 -0.0008 -0.0006 -0.0004 -0.0002 0 0.0002 40000 50000 60000 70000 Re 𝜔𝑖 Stability of an air-water flow in a semispherical container 91 Figure 6.5 (a) The 𝑧-profiles of the base-flow radial velocity 𝑢 (solid curve) and 𝐸𝑑 (dashed curve) near the lid at 𝑟 = 0.618 𝐻𝑤 = 0.8. The dotted curve depicts the Kármán solution for 𝑢. (b) Contours of the disturbance axial velocity at the interface (color online). This small value of 𝑣𝑠𝑖 results in that 𝑅𝑒𝑎 > 𝑅𝑒𝑤 that helps explain why the critical disturbance is focused in the air domain, as Figure 6.3(b) illustrates. It depicts contours of squared velocity magnitude, 𝐸𝑑(𝑟,𝑧) = constant, calculated for the critical disturbance; 𝐸𝑑=<|𝑢𝑑|2+|𝑣𝑑|2+ |𝑤𝑑|2> is normalized by its maximal value in the air and water domains respectively; brackets, < >, denote averaging with respect to 𝑡 and 𝜙. We depict a few contours of Ed with the outermost contour corresponding to 𝐸𝑑=0.1𝐸𝑑𝑚𝑎𝑥. This contour encompasses a flow region where the disturbance energy mostly is focused. In Figure 6.3(b), 𝐸𝑑 reaches its peak value at 𝑟 = 𝑟𝑚 = 0.618 and 𝑧 = 𝑧𝑚 = 0.957, i.e., near the rotating lid. At 𝑟 = 𝑟𝑚, Figure 6.5(a) depicts the axial distribution of the base-flow radial velocity 𝑢 (solid curve) and of 𝐸𝑑 (dashed curve) near the lid. For convenient presentation, 𝑢 is normalized by its maximal value, 𝑢𝑚 = 0.1065, in Figure 6.5(a), which reveals that both 𝑢 and 𝐸𝑑 are very small in the water flow (for 𝑧 < 0.8) compared with their maximal values in the air flow. The profile of 𝑢 in Figure 6.5(a) agrees with the streamline pattern in Figure 6.3(a). The rotating disk pushes the air to the periphery (the Kármán pump) forming a jet-like motion near the lid. The dotted curve in Figure 6.5(a) corresponds to the Kármán boundary layer solution as 𝑅𝑒𝑎   [127]. The near-lid parts (up to the 𝑢 peak) of the boundary layer and numerical solutions are close. For smaller 𝑧, the effect of the air circulation (the counterflow in Figure 6.5(a)) makes the difference. The fact that the 𝐸𝑑 peak is located near the inflection point of 𝑢(𝑧) in Figure 6.5(a) indicates that the instability likely is of the shear-layer type. The 𝐸𝑑 distribution in Figure 6.5(a) suggests that the air-flow instability does not necessarily disturb the water flow. A reason could be the large water-to-air density ratio. The 𝐸𝑑 contours shown in Figure 6.3(b) reveal that the 𝐸𝑑 peak in the water domain is located far away from the 𝐸𝑑 peak in the air domain. The 𝐸𝑑 peak in the water is likely induced by the counterflow in the thin water cell observed in Figure 6.3(a). Figure 6.5(b) depicts contours of critical-disturbance axial velocity at the interface. The bright (orange online) and dark (blue online) regions correspond to the interface deformation upward and downward, respectively. The locations of maximal |wd| in Figure 6.5(b) agree with the location of 𝐸𝑑 peak in the water in Figure 6.3(b). According to 𝑚 = 2, Figure 6.5(b) corresponds to two bumps and two dimples at the interface. As time increases, the picture in Figure 6.5(b) makes precession around the axis with the angular velocity r/m in the direction opposite to the base-flow rotation. 92 Global stability analysis of two-fluid flows 6.4.2. Flow at Hw = 0.6 Figure 6.6 (a) The base-flow streamlines and (b) critical disturbance energy contours at 𝐻𝑤 = 0.6. Figure 6.7 (a) The 𝑧-profiles of the base-flow radial velocity 𝑢 (solid curve) and 𝐸𝑑 (dashed curve) in the air domain at 𝑟 = 0.76 and 𝐻𝑤 = 0.6. The dotted curve depicts the Kármán solution for 𝑢. (b) Contours of the disturbance axial velocity at the interface (color online). At 𝐻𝑤 = 0.6, the critical parameters are 𝑅𝑒 = 76430, 𝜔𝑟 = 0.01295, 𝜔𝑖 = 0, 𝑚 = 2, 𝑅𝑒𝑎 = 5095, and 𝑅𝑒𝑤 = 876. The critical Reynolds numbers increase at 𝐻𝑤 = 0.6 compared to those at 𝐻𝑤 = 0.8. Figure 6.6(a) depicts the base-flow pattern at the critical 𝑅𝑒 and reveals that the Kármán boundary layer near the lid, where streamlines are packed, is now well separated, by a large region of slowly ascending air, from the near-interface boundary layer, corresponding to the radially converging jetlike air flow, where streamlines also are parked. The 𝑧-profile of 𝑢 in Figure 6.7(a) has the jet-like peak, corresponding to the Kármán boundary layer (dotted curve) near the lid and a jet-like valley, corresponding to the near-interface boundary layer. Figure 6.6(b) reveals that 𝐸𝑑 has peaks in both boundary layers near their inflection points of the u profile that indicates the shear-layer nature of instability. The near-interface jet diverges near the 𝐸𝑑 peak that stimulates the shear-layer instability [109]. Closer to the axis, this divergence is compensated by the radially converging motion. This explains why the 𝐸𝑑 peak in the near-interface jet is located close to the sidewall. Figure 6.7(b) depicts contours of critical-disturbance axial velocity at the interface whose features are similar to those in Figure 6.5(b). Stability of an air-water flow in a semispherical container 93 6.4.3. Flow at Hw = 0.4 Figure 6.8 (a) The base-flow streamlines and (b) critical disturbance energy contours at 𝐻𝑤 = 0.4. At 𝐻𝑤 = 0.4, the critical parameters are 𝑅𝑒 = 100620, 𝜔𝑟 = 0.0113, 𝜔𝑖 = 0, 𝑚 = 2, 𝑅𝑒𝑎 = 6708, and 𝑅𝑒𝑤 = 1078. The critical 𝑅𝑒 numbers increase at 𝐻𝑤 = 0.4 compared to those at 𝐻𝑤 = 0.6. Figure 6.8(a) depicts the base-flow pattern at the critical Re. It reveals that the Kármán boundary layer near the lid, where streamlines are packed, is now even more, than in the 𝐻𝑤 = 0.6 case, separated from the near-interface boundary layer. There is the axially elongated region in Figure 6.8(a) where the air streamlines are nearly parallel to the 𝑧 = axis. This motion is induced by the Kármán boundary layer [127], which sucks the ambient air. Figure 6.9 The 𝑧-profiles of the base-flow radial velocity 𝑢 (solid curve) and 𝐸𝑑 (dashed curve) in the air domain at 𝑟 = 0.469 and 𝐻𝑤 = 0.4. (b) Contours of the disturbance axial velocity at the interface (color online). Figure 6.9(a) depicts the 𝑧-profile of radial velocity 𝑢, normalized by its maximal value, 𝑢𝑚 = 0.0676, at 𝑟 = 0.469 and 𝐻𝑤 = 0.4. Near the interface, at z > 0.41, the multi-layer counterflow of air is formed where 𝑢 oscillates. The nature of these oscillations is similar to that for the secondary flow near a wall normal to a rotating fluid, first described by Bödewadt [134]. 94 Global stability analysis of two-fluid flows Figure 6.10 The 𝑟-profiles of the base-flow azimuthal 𝑣 and meridional 𝑣𝑡 velocities at the interface for 𝐻𝑤 = 0.4. Since both azimuthal and meridional velocities are very small at the interface compared with their maximal magnitudes, as Figure 6.10 illustrates, the interface acts as a wall for the air flow. Away from the interface, the radial pressure gradient is balanced by the centrifugal force, p/r = ʋ2/r. Near the interface, the centrifugal force is drastically reduced while the pressure gradient is not and, being unbalanced by the centrifugal force, pushes the fluid to the axis. The radial convergence of air accumulates the swirl vorticity near the axis that develops a local minimum of pressure. The minimal pressure sucks the ambient fluid, decelerates the downstream air flow, and reverses the radial velocity (see the peak of 𝑢 near 𝑧 = 0.5 in Figure 6.9(a)), i.e., the swirl-decay mechanism [113] works, generating the oscillations of 𝑢. The jet-like air flow near 𝑧 = 0.5 in Figure 6.9(a) has (i) inflection points and (ii) radially diverges. These two features, (i) and (ii), cause the strong shear-layer instability [109] that helps explains the location of the 𝐸𝑑 peak in Figure 6.8(b), where 𝑟𝑚 = 0.469 and 𝑧𝑚 = = 0.572, and in Figure 6.9(a), where the 𝐸𝑑 peak is located near the inflection point where 𝑢 is positive. Figure 6.9(b) depicts contours of critical-disturbance axial velocity at the interface whose features are similar to those in Figure 6.5(b). 6.4.4. Flow at Hw = 0.2 Figure 6.11 (a) The base-flow streamlines and (b) critical disturbance energy contours at 𝐻𝑤 = 0.2. In the shallow-water case, 𝐻𝑤 = 0.2, the disturbance with 𝑚 = 1 becomes the most dangerous with the critical values: 𝑅𝑒 = 95350, 𝜔𝑟 = 0.842, 𝜔𝑖 = 0, 𝑅𝑒𝑎 = 6357, and 𝑅𝑒𝑤 = 713. -0.0015 -0.001 -0.0005 0 0.0005 0 0.2 0.4 0.6 0.8 𝑟 𝑣,𝑣t 𝑣t 𝑣 Stability of an air-water flow in a semispherical container 95 Disturbances with 𝑚 = 2 and 𝑚 = 3 also become growing, but for larger values of 𝑅𝑒. The critical 𝑅𝑒 numbers slightly decrease at 𝐻𝑤 = 0.2 compared to those at 𝐻𝑤 = 0.4. Figure 6.11(a) depicts the base-flow pattern at the critical 𝑅𝑒. Figure 6.12 (a) The r-profiles of the base-flow normalized axial velocity 𝑤 (solid curve) and 𝐸𝑑 (dashed curve) in the air domain at 𝑧 = 0.289 and 𝐻𝑤 = 0.2. (b) Contours of the disturbance axial velocity at the interface (color online). Figure 6.11(b) depicts the 𝐸𝑑 distribution and reveals that the growing disturbance now is focused in the region where the near-wall jet-like air flow approaches the interface and decelerates. The deceleration enhances the shear-layer instability [109]. The 𝐸𝑑 peak is located at 𝑟 = 𝑟𝑚 = 0.638 and 𝑧 = 𝑧𝑚 = 0.289. Figure 6.12(a) shows the radial distribution of 𝐸𝑑 (dashed curve) and the axial velocity 𝑤, normalized by its minimal value, 𝑤𝑚𝑖𝑛 = 0.0532, (solid curve) at 𝑧 = 𝑧𝑚. The 𝐸𝑑 peak is located near the inflection point of 𝑤(𝑟) that indicates the shear-layer nature of the instability. Figure 6.12(b) depicts contours of critical-disturbance axial velocity at the interface. This pattern is spiral with the most high and low points located near the axis. 6.5. Concluding remarks The goal of this study is to explore the instability nature of the air-water flow in the semispherical container driven by the rotating lid. The study reveals that the instability develops in the air domain. The disturbance energy 𝐸𝑑 in the water is very small compared with the 𝐸𝑑 maximal magnitude in the air flow. At critical values of the Reynolds number, the flow has the well-developed three-eddy pattern: (i) the global meridional circulation of air is separated from (ii) the water co-circulation by (iii) the thin layer of water counter-circulation, which serves as a liquid bearing. The air circulation is concentrated in the jet-like boundary layer, (a) radially diverging near the lid (the Kármán flow), (b) descending near the sidewall, and (c) radially converging near the interface. The instability is caused by the inflection points in the boundary-layer velocity profile enhanced by the jet divergence and deceleration. The location of 𝐸𝑑 peak depends on the water volume fraction characterized by the depth-to-radius ratio 𝐻𝑤. For the deep water, the instability focuses near the rotating lid (Section 6.4.1). As 𝐻𝑤 decreases, the 𝐸𝑑 peak shifts downward (Sections 6.4.2 and 6.4.3) and concentrates near the interface-sidewall intersection for the shallow water (Section 6.4.4). The instability is timeoscillatory and helical with the azimuthal wave number 𝑚 = 2 except in the shallow water case, where 𝑚 = 1. Since 𝐸𝑑 drastically drops below the interface, it is not expected that the air-flow instability can significantly affect the bulk water motion. The results could have possible applications for aerial bio-reactors. 96 7. MEAN FLOW PRODUCED BY SMALLAMPLITUDE VIBRATIONS OF A LIQUID BRIDGE WITH ITS FREE SURFACE COVERED WITH AN INSOBLUBLE SURFACTANT s is well known, confined fluid systems subject to forced vibrations produce mean flows, called in this context streaming flows. These mean flows promote an overall mass transport in the fluid that has consequences in the transport of passive scalars and surfactants, when these are present in a fluid interface. Such transport causes surfactant concentration inhomogeneities that are to be counterbalanced by Marangoni elasticity. Therefore, the interaction of streaming flows and Marangoni convection is expected to produce new flow structures that are different from those resulting when only one of these effects is present. The present chapter focuses on this interaction using the liquid bridge geometry as a paradigmatic system for the analysis. Such analysis is based on an appropriate post-processing of the results obtained via direct numerical simulation of the system for moderately small viscosity, a condition consistent with typical experiments of vibrated millimetric liquid bridges. It is seen that the flow patterns show a nonmonotone behavior as the Marangoni number is increased. In addition, the strength of the mean flow at the free surface exhibits two well-defined regimes as the forcing amplitude increases. These regimes show fairly universal power-law behaviors. 7.1. Introduction Oscillatory flows in confined and low-viscosity fluid systems develop oscillatory boundary layers near solid walls and interfaces. Reynolds stresses exhibit nonzero temporal mean values in these layers, where they drive a mean flow. This interesting effect was first shown in a pioneering work by Lord Rayleigh [152] when explaining the anomalous sand accumulation at the bottom of a vertically vibrated container, observed 50 years earlier by Faraday [153]. Rayleigh also noted that the same mechanism was seemingly responsible for nonuniform dust accumulation (known as Kundt figures) at the walls of sound tubes [154] and for the formation of steady vortical structures in vibrated soap films, already observed by Sedley Taylor [155]. This latter problem has been more recently studied by Vega et al. [156]. They concluded that, under reasonableassumptions and according to Rayleigh’s insight, the observed vortical structures are most seemingly due to the streaming flow in the air surrounding the soap film. Rayleigh made a careful analysis of the (nowadays called) boundary layer attached to a (no-slip) solid boundary, obtaining the steady mean flow responsible for sand or dust nonuniform accumulation near the boundary. Rayleigh’s work was completed 50 years later by Schlichting [157] (see also Ref. [158]), who noticed that the streaming flow velocity was nonzero at the internal edge of the boundary layer, which is a driving mechanism for the generation of mean flows in the fluid bulk. Longuet-Higgins [159] made a similar analysis for the streaming flow within the boundary layer attached to a free surface, which again was nonzero at the internal edge of the layer and thus also responsible for the generation of a global mean flow in the bulk. An important well-known property is that the time-averaged velocity or stress converges to a generally nonzero A Mean flow produced by small-amplitude vibrations of a liquid bridge with its free surface covered with an insobluble surfactant 97 value as viscosity goes to zero, while it would vanish if viscosity were identically zero (i.e., when the boundary layer is absent) [125]. This seemingly surprising effect is explained by the singular perturbation character of the inviscid limit. More recently, the mean flow (also called steady or acoustic streaming [160], [161]) has been studied in various contexts. In particular, the streaming flow induced by no-slip boundaries has been studied in connection with flows in blood vessels [162], generation of mean motions in the ear [163], interaction of sound waves with bodies [164], flows around vibrating bodies [165], and streamingflow-jets produced by acoustic waves [166]. On the other hand, the streaming flow produced in the boundary layer attached to a vibrating free surface has an effect on the dynamics of oscillating bubbles [167], is of interest in water wave theory [168]–[172], and has been shown to play a role in the instability of the ocean Langmuir circulations [173]. In different contexts, the streaming flow has been studied in microfluidic systems [174] and in conjunction with thermal effects for various purposes, including the control of thermocapillary convection at low viscosity [175]–[178] and cooling in narrow channels [179]. For nearly resonant systems at small viscosity, the streaming flow is not just a byproduct of the primary waves, but generally affects the dynamics of those waves, as repeatedly shown for various flow configurations and vibrating devices [180]–[184]. This is because the streaming flow velocity scales with the square of the wave amplitude and produces an effect in the amplitude equations governing the dynamics of the primary waves that is proportional to the product of the primary waves amplitude and the streaming flow velocity. Thus, this latter effect scales as the cube of the primary waves amplitude, which is the order of magnitude of the cubic nonlinearity that saturates the dynamics in these systems [185]. The streaming flow is responsible for the effective mass transport in the fluid system and, in particular, for the overall motion of passive scalars and surfactants. Therefore, the mean flow is expected to exhibit a nontrivial interaction with Marangoni convection. The understanding of this interaction has been pursued in a number of papers. For instance, Martin and Vega [186] included the effect of Marangoni elasticity in a previous theory for clean free surfaces [182] to study the drift stability of standing Faraday waves in annular containers, previously found experimentally by Douady et al. [187]. More recently, Strickland et al. [188] have pointed out that materials (such as crude oil, biogenic slicks, or industrial and medical surfactants) absorbed at the fluid free surface are expected to move in response to surface waves. They have experimentally studied such effect for Faraday waves in a shallow cylindrical container with an insoluble surfactant monolayer. They distinguished the effects of the meniscus waves near the boundaries and the travelling waves in the far field, which resulted in very complex spatiotemporal patterns. The case of Faraday waves in a rectangular container was considered in Ref. [189], where no specific surfactant was added. In this study, surface contamination was seemingly present due to the use of tap water, which is easily contaminated. In fact, water contamination was also most probably responsible for some striking mean flow effects due to wave-maker oscillations that have been recently reported [190]. The streaming flow theory and systems mentioned above all deal with the low-viscosity limit, in which the thickness of the oscillatory boundary layers is very small compared to the characteristic lengths of the primary surface waves. This sets an upper limit to the surface waves oscillation frequency. If that condition does not hold, the streaming flow can still be analyzed by appropriate time-averaging of the trajectories of fluid elements calculated via direct numerical simulation. This task has not been addressed in the literature, to our knowledge. As mentioned above, Faraday waves in vibrating containers may produce meniscus oscillations that further complicate the dynamics, and may mask the interaction between the streaming flow and Marangoni convection. To overcome this difficulty, we consider a vibrating liquid bridge with the contact lines attached to the edges of the supporting disks, as we have done in other related works dealing with surfactant monolayer dynamics [191], [192]. This geometry has the additional advantage that direct numerical simulation is feasible at a reasonable computational price. The liquid bridge geometry has already been considered in connection with the streaming flow produced by vibrations, focusing on weakly nonlinear oscillations [193] and the effect of free surface waves [194]. However, these studies dealt with the low viscosity limit. On the contrary, viscosity is not 98 Global stability analysis of two-fluid flows necessarily small in the present study. Our analysis is based on direct numerical simulation for arbitrary values of that parameter. Concerning the interaction of the streaming flow with Marangoni convection, we may anticipate here that there is a tradeoff between both. Without Marangoni elasticity (the surfactant is treated as a passive scalar), the surfactant is transported by the streaming flow, which gives rise to nonuniform passive scalar concentrations at the free surface. However, Marangoni elasticity reacts against surfactant concentration gradients, decreasing the tangential velocity at the free surface and counterbalancing the effect of the streaming flow alone. Thus, the resulting streaming flow patterns are expected to be highly affected by Marangoni stresses, at least near the free surface. The main object of this investigation is to analyze the interaction between the streaming flow and Marangoni elasticity in a vibrated liquid bridge using direct numerical simulation. For simplicity, we consider an insoluble surfactant monolayer on top of the interface of a vibrating liquid bridge. We ignore gravity, consider the symmetric case in which the disks supporting the bridge have equal radii, and neglect surface viscosities. The streaming flow velocity and the time average of both the liquid bridge shape and the surfactant distribution will be calculated in two ways, one asymptotic assuming that the oscillatory flow intensity goes to zero, and another “exact” that does not take that assumption into account. In the former case, the streaming flow velocity is calculated as the sum of the mean Eulerian velocity and the Stokes drift, while in the latter case that velocity is obtained through the time-averaging of the trajectories of the fluid elements. These two approaches rely on the appropriate post-processing of the results obtained by direct numerical simulation, and their outcomes will be used for cross-checking. 7.2. Theoretical background As anticipated, the analysis of the interaction between the streaming flow and Marangoni stress will be conducted by post-processing “exact results” obtained upon direct numerical simulation. The theoretical background in Sec. 7.2.1 includes the problem formulation, the description of the numerical solver, and the nondimensional parameter range considered. The post-processing formulas to compute the streaming flow is derived and discussed in Sec. 7.2.2. 7.2.1. Governing equations and numerical solver We consider a liquid bridge (sketched in Figure 7.1) of length 𝐿 and volume 𝒱, held between two circular disks of equal radius 𝑅, and surrounded by another fluid medium of negligible density and viscosity. Nondimensionalization of the governing equations is made by using R and 𝑡0 = (𝜌𝑅3/𝜎0)1/2 as units for length and time, respectively, where 𝜌 is the liquid density and 𝜎0 is the equilibrium (initial) surface tension. For axisymmetric flows, the incompressible Navier-Stokes equations are (𝑟𝑢)𝑟+𝑟𝑤𝑟=0 (7.1) 𝑢𝑡+𝑢𝑢𝑟+𝑤𝑢𝑧=−𝑝𝑟+𝐶[𝑢𝑟𝑟+(𝑢 𝑟)𝑟+𝑢𝑧𝑧] (7.2) 𝑤𝑡+𝑢𝑤𝑟+𝑤𝑤𝑧=−𝑝𝑧+𝐶[𝑤𝑟𝑟+𝑤𝑟 𝑟+𝑤𝑧𝑧] (7.3) where r and z are the radial and axial cylindrical coordinates, respectively, t is the time variable, the Mean flow produced by small-amplitude vibrations of a liquid bridge with its free surface covered with an insobluble surfactant 99 Figure 7.1 Liquid bridge configuration subscripts 𝑟, 𝑧, and 𝑡 denote hereinafter partial derivatives with respect to those variables, 𝑢 and 𝑤 are the radial and axial velocity components, respectively, 𝑝 is the (hydrostatic) reduced pressure, and 𝐶 = 𝜇(𝜌𝜎0𝑅)−1/2 is the volumetric Ohnesorge number defined in terms of the bulk viscosity 𝜇. As boundary conditions, we consider the regularity conditions 𝑢=𝑤𝑟= 𝑝𝑟=0 at the symmetry axis 𝑟=0, as well as the kinematic compatibility condition 𝑓𝑡+𝑓𝑧𝑤−𝑢=0 (7.4) and the equilibrium of normal and tangential stresses, 𝑝−𝐵(𝑡)𝑧−𝜎𝜅−2𝐶[𝑢𝑟−𝑓𝑧(𝑤𝑟+𝑢𝑧)+𝑓𝑧2𝑤𝑧] 1+𝑓𝑧2=0 (7.5) 𝐶[1−𝑓𝑧2](𝑤𝑟+𝑢𝑧)+2𝑓𝑧(𝑢𝑟−𝑤𝑧) 1+𝑓𝑧2= 𝜎𝑧 √1+𝑓𝑧2 (7.6) at the free surface, 𝑟 = 𝑓 (𝑧,𝑡). It must be noted that these boundary conditions do not apply if overturning takes place, which may occur just before the pinch-off for vibration amplitudes larger than those considered in this paper. In the above equations, 𝜅=1+𝑓𝑧2−𝑓𝑓𝑧𝑧 𝑓(1+𝑓𝑧2)3/2 (7.7) is (twice) the mean curvature of the free surface, 𝐵(𝑡) is the unsteady Bond number defined in Eq. (7.12), and 𝜎= 𝜎/𝜎0 is the ratio of the local surface tension value 𝜎 to its equilibrium value 𝜎0. The dependence of the surface tension on the surfactant concentration obeys differentrelationships, depending on the surfactant. Here, we shall use the Szyskowski-Frumkin equation of state [195], 𝜎=1+𝑀𝑎 𝑙𝑛(1−𝛾 1−𝛾0) (7.8) where 𝛾 = 𝛾/𝛾∞ and 𝛾0 = 𝛾0/𝛾∞ are the local surface coverage and the equilibrium surface coverage, respectively. These quantities are defined in terms of the local surfactant surface concentration 𝛾, its equilibrium value 𝛾0, and its maximum value 𝛾∞ corresponding to the maximum 106 Global stability analysis of two-fluid flows Figure 7.4 Mean free surface deflection 𝑓𝐸(𝑧) (left) and mean surfactant concentration 𝛾𝐸(𝑧) (right). The values of the governing parameters are the same as those in Fig. 2, namely, Λ = 1.25, 𝑉 = 1, 𝐶 = 0.01, 𝛾0 = 0.5, 𝑀𝑎 = 0, 𝜔 = 2.1, and 𝐵 = 0.005. The simulation illustrated in Figs. Figure 7.3 and Figure 7.4 was conducted for Ma = 0. Therefore, the surfactant does not alter the surface tension and behaves as a passive scalar in that case. If the scalar were a surface active substance (𝑀𝑎 > 0), then this concentration distribution would cause Marangoni convection in the direction opposite to that of the streaming flow, which would alter significantly both the mean surfactant concentration distribution and the mean flow pattern. This effect will be analyzed in the next section. 7.3. Results on the interaction of the streaming flow and Marangoni convection All the results shown in this section were obtained for the values Eq. (7.14) of the governing parameters. The strength of the oscillatory and mean flows are measured using RMS values instead of values at particular points, which would be less robust and less appropriate to obtain scaling laws when different configurations are considered. Specifically, the intensity of the oscillatory flow is defined as 𝐼𝑜𝑠𝑐 𝑏𝑢𝑙𝑘=√1 Λ𝑉∫ ∫ [|𝑈(𝑟,𝑧)|2+|𝑊(𝑟,𝑧)|2]𝑟𝑑𝑟𝑑𝑧 𝑓𝐸 0 2Λ 0 (7.33) where 𝑈 and 𝑊 are given by Eq. (7.18). The counterpart of Eq. (7.33) for the streaming flow in the bulk is defined as 𝐼𝑚𝑓 𝑏𝑢𝑙𝑘=√1 Λ𝑉∫ ∫ [|𝑢𝐿(𝑟,𝑧)|2+|𝑤𝐿(𝑟,𝑧)|2]𝑟𝑑𝑟𝑑𝑧 𝑓𝐸 0 2Λ 0 (7.34) In addition, the intensity of the mean flow in the free surface is defined as 𝐼𝑚𝑓 𝑖𝑛𝑡𝑒𝑟𝑓=√1 𝐿𝑠∫[|𝑢𝐿|2+|𝑤𝐿|2]𝑟=𝑓𝐸𝑑𝑧 2Λ 0 (7.35) where Mean flow produced by small-amplitude vibrations of a liquid bridge with its free surface covered with an insobluble surfactant 107 𝐿𝑠=∫ √1+|𝑓𝑧𝐸|2𝑑𝑧 2Λ 0 (7.36) is the length of the time-averaged free surface. Finally, the strength of the mean flow along the axial coordinate is given by 𝐼𝑚𝑓 𝑎𝑥𝑖𝑎𝑙(𝑧)=√2∫ [|𝑢𝐿(𝑟,𝑧)|2+|𝑤𝐿(𝑟,𝑧)|2]𝑟𝑑𝑟 𝑓𝐸 0 (7.37) The components uL and wL of the mean Lagragian velocities in Eqs. (7.34), (7.35), and (7.37) are given by Eq. (7.29). In the following two subsections, we consider the influence on the streaming flow of the two main parameters: the Marangoni number and the strength of the oscillatory flow. As will be seen, the Marangoni number affects the shape of the flow patterns, while the strength of the oscillatory flow influences quantitatively the intensity of the mean flow. We aim at uncovering the main physical effects of Marangoni convection on the streaming flow, rather than a full account of the results obtained as the values of the governing parameters are varied. 7.3.1. Influence of the Marangoni number on the flow patterns In this section, we consider the dependence of the streaming flow on the Marangoni number for the parameter values listed in Eq. (7.14). We are mainly interested in illustrating how Marangoni stresses completely change the streaming flow pattern, without intending an exhaustive description of those patterns. For this reason, we will restrict ourselves to the Ohnesorge number 𝐶 = 0.01 and the forcing frequency 𝜔 = 2.1. Smaller values of 𝐶 would give flow patterns that are more concentrated near the disks and the free surface. The forcing frequency 𝜔 = 2.1 is close (but not equal) to the first resonance frequency of the liquid bridge for 𝑀𝑎 = 0, which means that the oscillatory flow intensity is very sensitive to the value of the dynamical Bond number 𝐵 . After some calibration, 𝐵 has been chosen such that the maximum over a cycle of 𝜀𝑊𝑒𝑖𝜔𝑡+𝑐.𝑐. at the center of the liquid bridge is 0.05. This choice makes the intensity of the streaming flow for 𝑀𝑎 = 0 comparable to the strength of the Marangoni convection in the considered range of Ma. It must be noted that the values of 𝐵 resulting from this criterion are physically realistic since they correspond to accelerations that can be easily produced in experiments [197]. Figure 7.5 shows the counterpart of Figure 7.3 for the indicated values of the Marangoni number. As can be seen, all the flow patterns are reflection-symmetric about the mid plane of the liquid bridge. However, and surprisingly at first sight, the streaming flow patterns do not show a monotone behavior as 𝑀𝑎 increases. Instead, the six counter-rotating eddies that are present for 𝑀𝑎 = 0 first merge for 𝑀𝑎 = 0.001, 0.005, and 0.05, and then give rise to just two counter-rotating eddies for 𝑀𝑎 = 0.1 and 0.2. As Ma increases further, new eddies appear leading to a flow pattern for 𝑀𝑎 = 1 (very strong Marangoni stress) that is qualitatively similar to that obtained for amuch smaller value 𝑀𝑎 = 0.001. This unexpected behavior is a consequence of the trade off between two mechanisms: nonuniformities in the surfactant concentration at the free surface, and the intensity of the oscillatory flow. 108 Global stability analysis of two-fluid flows Figure 7.5 Counterpart of Fig. 3 for different Marangoni numbers as indicated by the labels in the figure. The values of the governing parameters are Λ = 1.25, 𝑉 = 1, 𝐶 = 0.01, 𝛾0= 0.5, and 𝜔 = 2.1 Figure 7.6 Mean surfactant concentration 𝛾𝐸(𝑧). The colors (and types of lines) correspond to the different Marangoni numbers as indicated in the left-hand graph. The values of the rest of governing parameters are Λ = 1.25, 𝑉 = 1, 𝐶 = 0.01, 𝛾0 = 0.5, and 𝜔 = 2.1 The surfactant concentration along the free surface is shown in Figure 7.6 for the same values of Mean flow produced by small-amplitude vibrations of a liquid bridge with its free surface covered with an insobluble surfactant 109 the Marangoni number considered in Figure 7.5. As can be seen, the surfactant concentration becomes almost spatially uniform for very small Marangoni stress. This is due to the strong stabilizing effect of this stress, and is consistent with our guess at the end of Sec. 7.2.2. This behavior explains the streaming flow patterns shown in Figure 7.5 for 𝑀𝑎 = 0.001, 0.005, and 0.05. However, as 𝑀𝑎 is increased further, the strength of the oscillatory flow increases, peaks at 𝑀𝑎 ≅ 0.2, and decreases again [Figure 7.7(a)]. This is due to the fact that the forcing frequency matches the resonance frequency of the contaminated liquid bridge for 𝑀𝑎≅ 0.2. Since the intensity of the streaming flow depends on that of the oscillatory flow, they follow similar trends, as the comparison between Fig. Figure 7.7(a) and Figs. Figure 7.7(b) and Figure 7.7(c) shows. Figure 7.7 Intensities 𝐼𝑜𝑠𝑐 𝑏𝑢𝑙𝑘 (a), 𝐼𝑚𝑓 𝑏𝑢𝑙𝑘 (b), and 𝐼𝑚𝑓 𝑖𝑛𝑡𝑒𝑟𝑓 (c) of the oscillatory flow, the mean flow in the bulk and the mean flow in the interface, respectively, for the Marangoni numbers considered in Figure 7.5. The values of the rest of governing parameters are Λ = 1.25, 𝑉 = 1, 𝐶 = 0.01, 𝛾0 = 0.5, and 𝜔 = 2.1. Summarizing the above, the Marangoni stress has a great quantitative and qualitative influence on the streaming flow patterns. As 𝑀𝑎 increases, these patterns do not behave monotonously due to the fairly complex behavior of both the surfactant concentration at the free surface and the strength of the nearly resonant oscillation. 7.3.2. Influence of the oscillatory flow intensity In this section, we analyze the influence of the oscillatory flow intensity on the streaming flow intensity for the parameter values listed in Eq. (7.14). This analysis is more involved than that presented in the previous section due to the larger number of parameters considered, and the fact that hysteresis could be present because viscous effects are small and the forcing frequency is close to the 110 Global stability analysis of two-fluid flows resonant one. The results are obtained by varying the dynamical Bond number. However, and to uncover the underlying scaling laws for the intensity of streaming flow in the bulk and the interface, these quantities are plotted in terms of the oscillatory flow intensity. Figure 7.8 𝐼𝑚𝑓 𝑏𝑢𝑙𝑘 (left) and 𝐼𝑚𝑓 𝑖𝑛𝑡𝑒𝑟𝑓 (right) vs. 𝐼𝑜𝑠𝑐 𝑏𝑢𝑙𝑘for the indicated Marangoni numbers. The values of the rest of governing parameters are Λ = 1.25, 𝑉 = 1, 𝐶 = 0.01, 𝛾0 = 0.5, and 𝜔 = 2.1. Figure 7.8 shows the influence of 𝐼𝑜𝑠𝑐 𝑏𝑢𝑙𝑘 on 𝐼𝑚𝑓 𝑏𝑢𝑙𝑘 and 𝐼𝑚𝑓 𝑖𝑛𝑡𝑒𝑟𝑓for the indicated Marangoni numbers, 𝐶 = 0.01 and 𝜔 = 2.1. As can be seen, 𝐼𝑚𝑓 𝑏𝑢𝑙𝑘 scales as (𝐼𝑜𝑠𝑐 𝑏𝑢𝑙𝑘)2 for all values of 𝑀𝑎. This quadratic scaling law can be expected in view of Eq. (7.29). The curves 𝐼𝑚𝑓 𝑏𝑢𝑙𝑘 versus 𝐼𝑜𝑠𝑐 𝑏𝑢𝑙𝑘 are parallel to each other. The maximum and minimum values of the scaling factor differ by one order of magnitude. This factor does not behave monotonously as 𝑀𝑎 increases, but first decreases and then increases. This is consistent with the non-monotone behavior of the streaming flow patterns as 𝑀𝑎 increases (Figure 7.5). The curves 𝐼𝑚𝑓 𝑖𝑛𝑡𝑒𝑟𝑓versus 𝐼𝑜𝑠𝑐 𝑏𝑢𝑙𝑘 clearly show two power-law regimes: one for small 𝐼𝑜𝑠𝑐 𝑏𝑢𝑙𝑘 and another for large 𝐼𝑜𝑠𝑐 𝑏𝑢𝑙𝑘, with a fairly well defined threshold at 𝐼𝑜𝑠𝑐 𝑏𝑢𝑙𝑘≅0.02. Morover, the curves for 𝑀𝑎 > 0 are very close to each other, and clearly deviate from that for 𝑀𝑎 = 0. Figure 7.9 Counterpart of Figure 7.8 for the indicated oscillation frequencies. The values of the rest of governing parameters are Λ = 1.25, 𝑉 = 1, 𝐶 = 0.01, 𝛾0 = 0.5, and 𝑀𝑎 = 0.2. Mean flow produced by small-amplitude vibrations of a liquid bridge with its free surface covered with an insobluble surfactant 111 Figure 7.10 Counterpart of Figure 7.8 for the indicated values of 𝐶. The values of the rest of governing parameters are Λ = 1.25, 𝑉 = 1, 𝛾0 = 0.5, 𝑀𝑎 = 0.2, and 𝜔 = 2.1. The same behavior is found for other forcing frequencies, as shown in Figure 7.9 for 𝑀𝑎 = 0.2. As in the previous cases, this figure shows a unique scaling law for 𝐼𝑚𝑓 𝑏𝑢𝑙𝑘 and two scaling laws for 𝐼𝑚𝑓 𝑏𝑢𝑙𝑘 with a well defined threshold in 𝐼𝑜𝑠𝑐 𝑏𝑢𝑙𝑘. This threshold is very close to that in Figure 7.8. In addition, the scaling laws almost coincide for the three values of 𝜔, which suggests that these laws show a fairly weak dependence on the forcing frequency. The scaling behavior found in Figs. Figure 7.8 and Figure 7.9 is also present for smaller values of 𝐶, as shown in Figure 7.10. As can be seen, the scaling laws for 𝐼𝑚𝑓 𝑏𝑢𝑙𝑘 and 𝐼𝑚𝑓 𝑖𝑛𝑡𝑒𝑟𝑓 coincide with their counterparts in the previous cases. Moreover, the curve 𝐼𝑚𝑓 𝑏𝑢𝑙𝑘 versus 𝐼𝑜𝑠𝑐 𝑏𝑢𝑙𝑘 almost coincide for the two values of 𝐶, in spite of the fact that these are disparate. On the other hand, the curves 𝐼𝑚𝑓 𝑖𝑛𝑡𝑒𝑟𝑓 versus 𝐼𝑜𝑠𝑐 𝑏𝑢𝑙𝑘 slightly deviate from each other. The transition between the two scaling behaviors is less sharp for the smallest value of 𝐶. Summarizing the above, the intensity of the streaming flow in the interface shows two well defined power-law behaviors, one for small oscillatory flow intensity and another for larger values of this quantity. These behaviors are fairly universal, since they show the same exponents and the same transition threshold for the considered values of the Marangoni number, forcing frequency, and Ohnesorge number. Such regimes are reminiscent of those found in Ref. [189] (see Fig. 4 in that paper) for the streaming flow produced by Faraday waves in a rectangular container with tap water. However, we cannot ensure that there is a plain analogy between our problem and those experiments for two reasons. First, the data in Ref. [[189], Fig. 4] were not plotted on log-log scale, meaning that the power law for small oscillatory flow intensity could not be uncovered in that paper. In addition, tap water was most seemingly contaminated with a soluble surfactant in those experiments. Finally, we aim at elucidating whether hystheresis is present as the dynamic Bond number is varied. As already mentioned, this could occur since viscous effects are weak and the forcing frequency is close to the first resonant frequency. Figure 7.11 shows the intensities of the oscillatory and mean flows in the bulk obtained when the periodic regime is reached starting from equilibrium, as well as starting from the previous point while increasing or decreasing progressively the dynamical Bond number 𝐵 . As can be seen, the liquid bridge oscillations do not not exhibit hystheresis, i.e., they do not depend on the liquid bridge history. 112 Global stability analysis of two-fluid flows Figure 7.11 𝐼𝑜𝑠𝑐 𝑏𝑢𝑙𝑘 and 𝐼𝑚𝑓 𝑏𝑢𝑙𝑘 as a function of the dynamical Bond number 𝐵 starting from equilibrium, as well as starting from the previous point while increasing (up) and decreasing (down) progressively the dynamical Bond number. The symbols overlap. The values of the governing parameters are Λ = 1.25, 𝑉 = 1, 𝐶 = 0.01, 𝛾0 = 0.5, 𝑀𝑎 = 0.2, and 𝜔 = 2.1. 7.4. Conclusions We have developed a consistent theory of the interaction between the streaming flow produced by axial vibrations in a liquid bridge and Marangoni convection due to the presence of an insoluble surfactant monolayer. In contrast to previous papers in the field, this theory does not rely on any nearly inviscid formulation, but only on the assumption that the forcing amplitude is sufficiently low for the steepness of the free surface deflection to be small. The streaming flow intensities have been calculated from fully nonlinear direct numerical simulations. Two sets of formulas have been derived and used for cross-checking: one “exact” and based on the average along the trajectories of the fluid elements and another asymptotic and based on the addition of the Eulerian mean velocity and Stokes drift. After validation, the (computationally cheaper) asymptotic formulas have been used in the rest of computations. The main results of this investigation were obtained for a realistic set of parameter values. In particular, the forcing frequency was very close to the first resonant frequency of the liquid bridge. The following general conclusions can be established: (i) The streaming flow pattern is highly affected by Marangoni stress even for very small values of Ma. Moreover, the flow patterns show a nonmonotonous behavior as 𝑀𝑎 is increased, which is seemingly due to two main effects: the tendency of Marangoni stresses to oppose any surfactant inhomogeneities in the free surface and the nonmonotonous oscillatory response of the liquid bridge due to resonance peaks whose frequencies depend on 𝑀𝑎. (ii) The strength of the streaming flow shows a “universal” behavior with a unique power law for the intensity in the bulk and a transition between two power laws for the intensity in the free surface. This could be seen as somewhat similar to some recent findings in the literature [189]. The above conclusions are expected to apply to other vibrating systems despite the peculiarities of the primary oscillating flow producing the mean flow. In fact, we have selected the liquid bridge geometry only because it shows fairly simple oscillatory flow patterns when the first resonant mode is excited. This allowed us to focus on the main involved physical mechanisms. The scaling laws obtained in the present analysis require that the forcing amplitude be small. For larger amplitudes, the dynamics are fully nonlinear and much more complex. In this case, the Mean flow produced by small-amplitude vibrations of a liquid bridge with its free surface covered with an insobluble surfactant 113 oscillatory and mean flow components of the state variables are not expected to be related by simple laws. On the other hand, the somewhat small viscous effects considered in this paper might also play a role. Although the quadratic dependence of the mean flow on the vibrating amplitude also holds for more viscous fluids, the interaction with Marangoni convection may change that dependence for not too small vibrating amplitudes. The analysis of this possibility is beyond the scope of this study, where only realistic values of the parameters have been considered. Most of surface-active substances in nature are water-soluble. Surfactant solubility makes our analysis less general and does affect the results of the study quantitatively (in particular, the scaling laws). However, there is a large set of surfactants characterized by adsorption and desorption time scales much larger than the liquid bridge capillary time. The results presented here can be extended to such substances. In any case, the fact that Marangoni elasticity opposes surfactant inhomogeneities produced by the streaming flow remains valid for soluble surfactants as well. 7.5. Appendix: Computation of the lagrangian mean velocities in the numerical domain As explained in Sec. 7.2.1, the computational spatial coordinates are 𝜂 ≡ 𝑟/𝑓(𝑧,𝑡) and 𝑧. The counterparts of Eqs. (7.15), (7.18), and (7.20) for the velocity components are (𝑢,𝑤)=𝜀[(𝑈 (𝜂,𝑧),𝑊 (𝜂,𝑧))𝑒𝑖𝜔𝑡+𝑐.𝑐]+𝜀2[(𝑈 2(𝜂,𝑧),𝑊 2(𝜂,𝑧))𝑒2𝑖𝜔𝑡+𝑐.𝑐+ (𝑢𝑚(𝜂,𝑧),𝑤𝑚(𝜂,𝑧))]+(𝑂𝜀3), (7.38) 𝜀(𝑈 (𝜂,𝑧),𝑊 (𝜂,𝑧))=1 𝑇∫(𝑢(𝜂,𝑧,𝑡),𝑤(𝜂,𝑧,𝑡))𝑒−𝑖𝜔𝑡𝑑𝑡 𝑇 0 (7.39) 𝜀2(𝑈 2(𝜂,𝑧),𝑊 2(𝜂,𝑧))=1 𝑇∫(𝑢(𝜂,𝑧,𝑡),𝑤(𝜂,𝑧,𝑡))𝑑𝑡 𝑇 0 (7.40) where the tilde indicates that the quantity is evaluated in the fixed numerical domain. In addition, Eq. (7.16) implies that 𝜂=𝑟 𝑓=𝑟 𝑓0−𝜀[𝑟𝐹 𝑓02𝑒𝑖𝜔𝑡+𝑐.𝑐.]+(𝑂𝜀2) (7.41) The Eulearian mean velocity at a fixed point (𝑟,𝑧) can be calculated from the velocity field in the numerical domain as (𝑢𝐸(𝑟,𝑧),𝑤𝐸(𝑟,𝑧))=1 𝑇∫(𝑢(𝑟/𝑓(𝑡),𝑧,𝑡),𝑤(𝑟/𝑓(𝑡),𝑧,𝑡))𝑑𝑡 𝑇 0 (7.42) This calculation involves the evaluation of the velocity field at a moving point (r/f (t),z) of the numerical domain in the course of the oscillation. To avoid this evaluation, the above integral is expanded in the following way 114 Global stability analysis of two-fluid flows 1 𝑇∫(𝑢(𝑟/𝑓(𝑡),𝑧,𝑡),𝑤(𝑟/𝑓(𝑡),𝑧,𝑡))𝑑𝑡 𝑇 0=1 𝑇∫ (𝑢(𝑟 𝑓0,𝑧,𝑡),𝑤(𝑟 𝑓0,𝑧,𝑡))𝑑𝑡 𝑇 0 −𝜀2[𝑟𝐹(𝑧) 𝑓02(𝑈𝜂(𝑟 𝑓0,𝑧),𝑊𝜂(𝑟 𝑓0,𝑧))+𝑐.𝑐]+(𝑂𝜀3) (7.43) Therefore, the Eulearian mean velocity at a fixed point (𝑟,𝑧) can be obtained as (𝑢𝐸,𝑤𝐸)=𝜀2[(𝑢𝑚(𝜂0,𝑧),𝑤𝑚(𝜂0,𝑧))−𝜂0𝐹(𝑧) 𝑓0(𝑧) (𝑈 𝜂(𝜂0,𝑧),𝑊 𝜂(𝜂0,𝑧))+𝑐.𝑐] +(𝑂𝜀3) (7.44) where 𝜂0≡𝑟/𝑓0(𝑧) is a fixed coordinate in the numerical domain. To leading order, 𝑈 (𝑟/𝑓,𝑧) = 𝑈 (𝑟/𝑓0,𝑧) and 𝑊 (𝑟/𝑓,𝑧)=𝑊 (𝑟/𝑓0,𝑧). Therefore, the Lagrangian mean velocity at a fixed point (𝑟,𝑧) can be calculated as (𝑢𝐿,𝑤𝐿)=𝜀2(𝑢𝑚,𝑤𝑚) +𝜀2[−𝜂0𝐹 𝑓0(𝑈 𝜂,𝑊 𝜂)+1 𝑖𝜔(𝑈 𝜂𝑈 +𝑈  𝑧𝑊 ,𝑊  𝜂𝑈 +𝑊  𝑧𝑊 )+𝑐.𝑐.] +𝑂(𝜀3) (7.45) where all the functions are evaluated at (𝜂0,𝑧). Finally, the kinematic Eqs. (7.23) are rewritten as 𝑑𝜂 𝑑𝑡=𝑢(𝜂,𝑧,𝑡) 𝑓(𝑧,𝑡) −𝜂𝑓𝑧(𝑧,𝑡)𝑤(𝜂,𝑧,𝑡) 𝑓(𝑧,𝑡)−𝜂𝑓𝑡(𝑧,𝑡) 𝑓(𝑧,𝑡), 𝑑𝑧 𝑑𝑡=𝑤(𝜂,𝑧,𝑡) (7.46) 115 8. INFLUENCE OF THE DYNAMICAL FREE SURFACE DEFORMATION ON THE STABILITY OF THERMAL CONVECTION IN HIGH-PRANDTL- NUMBER LIQUID BRIDGES e analyze theoretically the stability of the thermal convection in high-Prandtl-number liquid bridges. The steady axisymmetric base flow, as well as its corresponding linear non-axisymmetric eigenmodes, are calculated taking into account the free surface deformation caused by both that flow and the perturbations. The stability limits and the oscillation frequencies obtained from the linear stability analysis satisfactorily agree with previous experimental data. The dynamical free surface deformation produced by the base flow approximately coincides with that measured in the experiments. When the deformations are normalized with their corresponding values of the Capillary number, they collapse onto a single curve. The dependence of the free surface oscillation amplitude with respect to the axial coordinate approximately also coincides with that measured in previous experiments under similar conditions. Our results show that the dynamical free surface deformation has very little effect on the eigenvalues, and, therefore, on the stability limits. 8.1. Introduction Convective flows in systems with free surfaces have attracted attention both at the fundamental and practical levels over the last decades. A good example of this is the thermocapillary convection in liquid bridges under microgravity conditions, where flow sets in due to variations of the surface tension along the free surface caused by the temperature difference ∆𝑇 applied between the two supporting rods. For small ∆𝑇, an axisymmetric steady flow evolves driven by thermocapillarity. For ∆𝑇 > ∆𝑇𝑐 and high-Prandtl-number liquids, this flow bifurcates into a three-dimensional oscillatory one (Velten et al., [198]). There have been numerous both theoretical and experimental studies of thermocapillary convection in cylindrical liquid bridges. Attention has mainly been paid to the examination of the onset of instability and, more recently, the formation of Particle Accumulation Structures (PAS) (Schwabe et al., [199]). Direct numerical simulations of the three-dimensional, time-dependent problem (Shevtsova et al., [200]) have enabled the description of the supercritical regimes arising as ∆T increases above its critical value. These studies have been extended to normal gravity conditions and non-cylindrical volumes (𝒱≠1,𝒱 being the liquid bridge volume divided by that of the cylinder enclosed by the supporting rods). In almost all the cases, the numerical simulations have been conducted with a statically deformed free surface calculated for ∆𝑇 = 0 from the equilibrium Young-Laplace equation. For high-Prandtl-number liquids, one typically obtains a stability diagram ∆𝑇𝑐(𝒱 ) which consists of two branches that approach each other at very high values of the critical temperature difference ∆𝑇𝑐 when the liquid bridge volume 𝒱 tends to a value around unity. W 122 Global stability analysis of two-fluid flows grid points, and ℱ(𝑝) represents the inhomogeneous part of the equations. Both the elements of ℐ(𝑝,𝑞) and ℱ(𝑝) are computed analytically before running the simulation using a symbolic software. The resulting functions are evaluated numerically over the Newton-Raphson iterations, which greatly reduces the CPU time. In each of those iterations, one first evaluates the Jacobian ℐ𝑏(𝑝,𝑞)= ℐ(𝑝,𝑞)(𝒰𝑏 (𝑞)) and the inhomogeneous part ℱ𝑏(𝑝)(𝒰𝑏 (𝑞)) for the updated value of 𝒰𝑏 (𝑞). Then, one calculates the inverse matrix ℐ𝑏−1(𝑝,𝑞) with an algorithm which takes advantage of the high sparsity of the matrix. Finally, the correction vector 𝛿𝒰𝑏 (𝑞) is obtained as 𝛿𝒰𝑏 (𝑞)=−ℐ𝑏−1(𝑝,𝑞)ℱ𝑏(𝑝). The above numerical procedure can be easily adapted to solve the eigenvalue problem which determines the linear modes of the system. The spatial dependence of the linear perturbation 𝒰(𝑞) is the solution to the generalized eigenvalue problem ℐ𝑏(𝑝,𝑞)𝒰(𝑞)=𝑖𝜔𝒬(𝑝,𝑞)𝒰(𝑞), where ℐ𝑏(𝑝,𝑞) is the Jacobian of the system evaluated with the base solution 𝒰𝑏 (𝑞), and 𝒬(𝑝,𝑞) accounts for the temporal dependence of the problem. This matrix is calculated with essentially the same procedure as that for ℐ𝑏(𝑝,𝑞). Therefore, the numerical algorithm developed for the base flow problem can also be applied to the linear stability analysis. More details of the numerical method can be found elsewhere (Herrada and Montanero, [1]). We have validated the numerical method from comparison with previous simulation results for a rigid free surface. The average difference between the critical Reynolds numbers and oscillation frequencies shown in Tables 1 and 3 of Ref. Shevtsova et al. [200] and our corresponding values are 0.2% and 0.6%, respectively. The difference between the critical Rayleigh number shown in Table 1 of Ref. Kuhlmann and Rath [217] for the natural convection in a liquid bridge with 𝐴𝑅=0.5 and that predicted by our simulation is less than 0.5%. Finally, the difference between the critical Reynolds number shown in Table 1 of Ref. (Prange et al., [224]) in the absence of a magnetic field and our result is less than 0.5%. 8.4. Results It is believed that the onset of instability under microgravity conditions is determined by the Marangoni 𝑀𝑎=𝑅𝑒𝑃𝑟=𝜎𝑇∆𝑇𝐿/(𝜌0𝜈0𝛼). and Capillary 𝐶𝑎=𝑅𝑒𝑊𝑒=𝜎𝑇∆𝑇/𝜎0 numbers rather than by the Reynolds and Weber numbers, which naturally arise from the characteristic quantities introduced in Sec. 8.2. For this reason, we will adopt {𝒱,𝐴𝑅,𝑃𝑟,𝐹𝑟,𝐵𝑖,𝛼𝜈,Θ𝑎𝑚𝑏,𝐺𝑟,𝑀𝑎,𝐶𝑎,𝑇0∗} (8.35) as the set of governing parameters. The last five numbers involve the temperature difference ∆𝑇. Now, the Capillary number is the single parameter that includes the equilibrium surface tension 𝜎0, and, therefore, it can be regarded as a measure of the free surface rigidity. In fact, the limit 𝐶𝑎→0 (𝜎0→∞) corresponds to the rigid free surface approximation commonly made in the theoretical analyses. Given the large dimension of the parameter space, we will restrict ourselves to liquid bridges of 5-cSt silicone oil, characterized by a high Prandtl number (𝑃𝑟=67) and very frequently used in experiments. The physical properties of this liquid measured at 𝑇0 = 25 oC by Nishino et al. ([206]): 𝜌0 = 915 kg/m3, 𝛽 = 1.09 × 10−3 1/K, 𝜈0= 5.0 × 10−6 m2/s, 𝛼𝜈 = 5.892, 𝜎0 = 19.7 × 10−3 N/m, 𝜎𝑇 = 6.58 × 10−5 N/(m·K), 𝜅 = 0.12 W/(m·K) and 𝛼 = 7.46 × 10−8 m2/s. Figure 8.3 shows a comparison between the experimental values and our numerical predictions for the critical Marangoni number 𝑀𝑎∗=𝜎𝑇∆𝑇(𝐷 2)/(𝜌𝜈𝛼) and the dimensionless oscillation frequency 𝜔𝑟∗=𝜔𝑟(𝐷 2)2/(𝛼√𝑀𝑎∗) of liquid bridges with different aspect ratios under microgravity conditions (𝐹𝑟=𝐺𝑟=0) (Fig. 6 of Ref. Nishino et al. ([206])). Here, 𝜈 and 𝜔𝑟 are the Influence of the dynamical free surface deformation on the stability of thermal convection in high-prandtl-number liquid bridges 123 kinematic viscosity for 𝑇=𝑇0+∆𝑇/2 and the dimensional oscillation frequency, respectively. In our simulations, the instability for 𝐴𝑅= 0.17 and 0.325 was caused by the growth of the 𝑚 = 2 mode, while the dominant mode for the rest of aspect ratios was that with 𝑚 = 1. There is good agreement between the theoretical and experimental results for the Biot numbers 𝐵𝑖 = 0.15 and 0.3 proposed by Nishino et al. [206]. This agreement is similar to that exhibited by the linear stability analysis presented by in that work. In particular, the simulation reproduces the sharp decrease of the oscillation frequency for 𝐴𝑅= 1.25. Figure 8.3 Experimental (solid symbols) and numerical (open symbols) values of the critical Marangoni number 𝑀𝑎∗ (left-hand graph) and the oscillation frequency 𝜔𝑟∗ (right-hand graph) for liquid bridges with different aspect ratios AR under microgravity conditions (𝐹𝑟=𝐺𝑟=0) (see Fig. 6 of Ref. (Nishino et al., [206])). The circles and triangles correspond to the numerical results calculated with 𝐵𝑖 = 0.15 and 0.3, respectively. Montanero et al. [205] studied experimentally the dynamical free surface deformation arising in the thermal convection of 5-cSt silicone oil liquid bridges. In those experimental runs, the liquid bridge was formed at the ambient temperature 𝑇𝑎𝑚𝑏= 25oC, and then the upper rod was progressively heated up while keeping the lower rod temperature 𝑇0 = 25oC constant, so that a certain range of ∆T values was covered. The dimensionless number {𝐴𝑅=0.615,𝑃𝑟=67,𝛼𝜈=5.892,𝐹𝑟= 1.973×104,Θ𝑎𝑚𝑏=0} were kept constant, while {𝒱,𝐺𝑟,𝑀𝑎,𝐶𝑎,𝑇0∗} varied with ∆𝑇. The free surface location was accurately measured in the course of the experiment with a sub-pixel resolution technique. The liquid thermal expansion produced variations of the liquid bridge volume ∆𝒱 ∼ 10−2, which caused free surface deformations comparable to those produced by the flow. We will devote the rest of this section to examine numerically the first experimental run conducted in that work (the other produced very similar results). In that case, the instability was caused by the 𝑚 = 1 mode for ∆𝑇 = 16.8oC, which corresponds to a critical Marangoni number 𝑀𝑎= 1.196 × 104. In order to take into account thermal expansion, our simulations were conducted with the liquid bridge volumes measured in the experiment for the corresponding value of ∆𝑇 (see Fig. 10b of Ref. (Montanero et al., [205])). Figure 8.4 shows the base flows corresponding to the experimental run I in Ref. Montanero et al. [205] with ∆𝑇 = 15.20 and 21.76oC, which correspond to subcritical and supercritical conditions, respectively. In both cases, the streamlines show the formation of four vortices. The two anticlockwise vortices near the free surface are driven by the Marangoni stress, while the other two clockwise rotations are caused by the buoyancy force. A thin recirculation cell arises in contact with the upper rod, where the temperature gradient reaches higher values. As can be observed, a large number of grid points are required in that region and next to the free surface to reproduce the flow accurately. The only qualitative difference between the flow patterns in the subcritical and supercritical cases is the growth of the vortex next to the center of the lower rod as the temperature difference increases. 124 Global stability analysis of two-fluid flows We have verified that neither the base flow nor its linear stability are significantly affected by the Biot number. Specifically, the critical Marangoni number 𝑀𝑎 = 1.159 × 104 calculated for 𝐵𝑖=0.15 differed in less than 0.5% from that calculated for 𝐵𝑖=0.3. Those values deviated in less than 3.6% from that measured in the experiment. We will set 𝐵𝑖=0.15 in the rest of the simulations. As occurred in the experiment, instability was caused by the 𝑚 = 1 mode. Figure 8.4 Streamlines (left) and temperature isolines (right) for {𝐴𝑅=0.615,𝑃𝑟=67,𝛼𝜈= 5.892,𝐹𝑟=1.973×104,Θ𝑎𝑚𝑏=0,𝐵𝑖=0.15} and {𝒱=0.8265,𝐺𝑟=326.9,𝑀𝑎= 1.082𝑥104,𝐶𝑎=0.0508,𝑇0∗=1.645}(left) and {𝒱=0.8298, 𝐺𝑟=468,𝑀𝑎=1.549𝑥104 𝐶𝑎=0.0727,𝑇0∗=1.149} (right). The blue and cyan lines in the left-hand figure indicate clockwise and anticlockwise rotations, respectively. As mentioned above, Montanero et al. [205] described the free surface deformation of a 5-cSt silicone oil liquid bridge when a temperature difference ∆T was applied between two solid supports of the same diameter. For ∆𝑇 = 0, the equilibrium liquid bridge contour is characterized by the function 𝐹0(𝑧) which obeys the Young-Laplace equation (8.34). The deformation 𝐹0(𝑧)− (2 𝐴𝑅)−1 with respect to the cylinder delimited by the two supports results from the competition between gravity and surface tension, and from the difference between the liquid bridge volume and that of the cylinder. An axisymmetric steady flow arises in the liquid bridge when a temperature difference ∆𝑇 is applied between the supports. The corresponding free surface contour 𝐹𝑏(𝑧) differs from the equilibrium shape 𝐹0(𝑧) by the perturbation 𝐹𝑏(𝑧)−𝐹0(𝑧), which is at least one order of magnitude smaller than the deformation 𝐹0(𝑧)−(2 𝐴𝑅)−1 at equilibrium under normal gravity conditions (Montanero et al., [205]). This perturbation is caused by the liquid thermal expansion, the surface tension variation over the free surface, and the hydrodynamic effects on that surface. The three effects are commensurate with each other. When ∆𝑇 exceeds a certain critical value ∆𝑇𝑐, the axisymmetric steady flow evolves towards a non-axisymmetric oscillatory regime. In this regime, the free surface contour 𝐹(𝜃,𝑧;𝑡) can be seen as the sum of those corresponding to the underlying base flow, 𝐹𝑏(𝑧), and to the oscillations. In the experiments, the oscillation contribution 𝐹(𝜃,𝑧;𝑡)− 𝐹𝑏(𝑧) was a periodic or quasi-periodic function of time, and at least one order of magnitude smaller than 𝐹𝑏(𝑧)−𝐹0(𝑧) for the range of ∆𝑇 explored (Ferrera et al., [206]). The underlying base flow contour 𝐹𝑏(𝑧) for ∆𝑇>∆𝑇𝑐 can be approximately measured by averaging 𝐹(𝜃,𝑧;𝑡) over an oscillation cycle (Montanero et al., [205]). To summarize, one must bear in mind the following hierarchy of deformations for the thermal convection in liquid bridges: |𝐹(𝜃,𝑧;𝑡)−𝐹𝑏(𝑧)|≪|𝐹𝑏(𝑧)−𝐹0(𝑧)|≪|𝐹0(𝑧)−(2 𝐴𝑅)−1|≪(2 𝐴𝑅)−1 (8.36) Influence of the dynamical free surface deformation on the stability of thermal convection in high-prandtl-number liquid bridges 125 Figure 8.5 Experimental (open symbols) and numerical (solid symbols) free Surface deformation 𝐹𝑏−𝐹0 of the base flow {𝐴𝑅=0.615,𝑃𝑟=67,𝛼𝜈=5.892,𝐹𝑟=1.973×104,Θ𝑎𝑚𝑏= 0,𝐵𝑖=0.15} and {𝒱=0.8200,𝐺𝑟=117.4,𝑀𝑎=3.885𝑥103,𝐶𝑎=0.0182,𝑇0∗=4.579} (diamonds), {𝒱=0.8241,𝐺𝑟=237.7,𝑀𝑎=7.863𝑥103,𝐶𝑎=0.0369,𝑇0∗=2.262} (downtriangles), {𝒱=0.8265,𝐺𝑟=326.9,𝑀𝑎=1.082𝑥104,𝐶𝑎=0.0508,𝑇0∗=1.645} (uptriangles), {𝒱=0.8298,𝐺𝑟=468,𝑀𝑎=1.549𝑥104,𝐶𝑎=0.0727,𝑇0∗=1.149} (circles) and {𝒱=0.8405,𝐺𝑟=544.8,𝑀𝑎=1.803𝑥104,𝐶𝑎=0.0846,𝑇0∗=0.987} (squares). Figure 8.5 shows 𝐹𝑏(𝑧)−𝐹0(𝑧) both measured experimentally (Montanero et al., [205]) and obtained from our numerical simulations. As mentioned above, the simulation was conducted with the liquid bridge volume determined in the experiment for the corresponding value of ∆𝑇 to account for thermal expansion. In both the experiments and simulations, the function 𝐹0(𝑧) was calculated from the integration of the Young-Laplace equation (8.34) with that volume and the surface tension value 𝜎0. Therefore, 𝐹𝑏(𝑧)−𝐹0(𝑧) indicates the free surface deformation caused by the surface tension variation over the free surface and the hydrodynamic effects on that surface. There is good agreement between the theoretical predictions and experiments. In all the cases, the deviation of the numerical contour from its experimental counterpart is smaller than the (dimensionless) size 2.3 x 10- 3 of the pixels in the images acquired in the experiments (Montanero et al., [205]). Part of that deviation can be attributed to the small differences between the nominal geometry and the real one (eccentricity, difference between the diameters of the rods, . . .), which may commensurate with the dynamical perturbation 𝐹𝑏(𝑧)−𝐹0(𝑧) [Eq. (8.36)]. The free surface deformation induced by the base flow exhibits similar shapes for the temperature differences considered in Figure 8.5. In fact, the maximum and minimum deformations are reached at practically the same heights. This can be clearly appreciated in Figure 8.6, where the free surface deformation has been divided by the Capillary number. As can be seen, the maximum free surface deformation is one order of magnitude smaller than the corresponding value of the Capillary number. The normalized deformations do not perfectly match probably due to the influence of the liquid bridge equilibrium shape 𝐹0(𝑧). This shape is different in each case due to the thermal expansion, as quantified by the volume . To verify this hypothesis, we have run all the simulations keeping constant that parameter, and, therefore, the equilibrium contour 𝐹0(𝑧). The results are shown in the right-hand graph of Figure 8.6. The normalized free surface deformations fall on a single curve ℱ(𝑧), which indicates that 𝐹𝑏(𝑧)−𝐹0(𝑧)=𝐶𝑎 ℱ for 𝐶𝑎≲0.1 (provided that the volume remains constant). This means that the dynamical free surface deformation can be seen as a linear steady perturbation with respect to the liquid bridge equilibrium shape. We have explored the validity of the above linear relationship by conducting simulations with larger values of the Capillary number 126 Global stability analysis of two-fluid flows (applied temperature difference). The above linear law underestimates the dynamical free surface deformation for 𝐶𝑎≳0.15, which corresponds to the supercritical regime. tr Figure 8.6 (Left) Free Surface deformation normalized with the Capillary number, (𝐹𝑏−𝐹0)/ 𝐶𝑎, for the configurations considered in Figure 8.5: {𝐴𝑅=0.615,𝑃𝑟=67,𝛼𝜈=5.892,𝐹𝑟= 1.973×104,Θ𝑎𝑚𝑏=0,𝐵𝑖=0.15} and {𝒱=0.8200,𝐺𝑟=117.4,𝑀𝑎=3.885𝑥103,𝐶𝑎= 0.0182,𝑇0∗=4.579} (diamonds), {𝒱=0.8241,𝐺𝑟=237.7,𝑀𝑎=7.863𝑥103,𝐶𝑎= 0.0369,𝑇0∗=2.262} (down-triangles), {𝒱=0.8265,𝐺𝑟=326.9,𝑀𝑎=1.082𝑥104,𝐶𝑎= 0.0508,𝑇0∗=1.645} (up-triangles), {𝒱=0.8298,𝐺𝑟=468,𝑀𝑎=1.549𝑥104,𝐶𝑎= 0.0727,𝑇0∗=1.149} (circles) and {𝒱=0.8405,𝐺𝑟=544.8,𝑀𝑎=1.803𝑥104,𝐶𝑎= 0.0846,𝑇0∗=0.987} (squares). Figure 8.7. (Left-hand image) Isolines of kinetic energy 1/2(𝑈 2+𝑉2+𝑊 2) (left) and magnitude of |𝑝| of reduced pressure (right). (Right-hand image) Isolines of magnitude of the temperature field, |Θ  (𝑟,𝑧)| (left), and amplitude |𝐹(𝑧)| of the free surface oscillation. The latter has been magnified to show clearly its axial dependence. The results correspond to the eigenmode 𝑚=1 responsible for instability at {𝐴𝑅=0.615,𝑃𝑟=67,𝛼𝜈=5.892,𝐹𝑟=1.973× 104,Θ𝑎𝑚𝑏=0,𝐵𝑖=0.15} and {𝒱=0.828,𝐺𝑟=350.2,𝑀𝑎=1.159𝑥104,𝐶𝑎=0.0544,𝑇0∗= 1.535} Figure 7 shows the spatial structure of the dominant eigenmode for the marginally stable base flow of the experimental run analyzed in this section. The perturbation kinetic energy sharply increases in the region next to the free surface where the liquid flows between the two vortices driven by the Marangoni stress (Figure 8.4). The magnitude of the reduced pressure perturbation increases close to the triple contact line anchored to the upper rod. The upper part of the free surface oscillates more than the lower part. The shape of the free surface oscillation amplitude is very akin to that measured experimentally by Ferrera et al. [206] under similar conditions (see Fig. 9 of that Influence of the dynamical free surface deformation on the stability of thermal convection in high-prandtl-number liquid bridges 127 reference). It also resembles the experimental results obtained more recently by Yano et al. [207] for the oscillatory Marangoni convection under microgravity conditions (see Fig. 10 of that reference). It must be pointed out that the linear stability analysis cannot predict the magnitude of the free surface oscillations, but only the relative dependence of the amplitude with respect to the axial coordinate. Figure 8.8 Eigenvalues characterizing the m=1 (circles) and 2 (triangles) linear modes for the configurations considered in Figure 8.5: {𝐴𝑅=0.615,𝑃𝑟=67,𝛼𝜈=5.892,𝐹𝑟=1.973× 104,Θ𝑎𝑚𝑏=0,𝐵𝑖=0.15} and {𝒱=0.8200,𝐺𝑟=117.4,𝑀𝑎=3.885𝑥103,𝐶𝑎= 0.0182,𝑇0∗=4.579} , {𝒱=0.8241,𝐺𝑟=237.7,𝑀𝑎=7.863𝑥103,𝐶𝑎=0.0369,𝑇0∗=2.262} , {𝒱=0.8265,𝐺𝑟=326.9,𝑀𝑎=1.082𝑥104,𝐶𝑎=0.0508,𝑇0∗=1.645}, {𝒱=0.8298,𝐺𝑟= 468,𝑀𝑎=1.549𝑥104,𝐶𝑎=0.0727,𝑇0∗=1.149} and {𝒱=0.8405,𝐺𝑟=544.8,𝑀𝑎= 1.803𝑥104,𝐶𝑎=0.0846,𝑇0∗=0.987}. The open and solid symbols were calculated with dynamical free surface deformation and the equilibrium contour, respectively. Figure 8.8 shows the eigenvalues characterizing the 𝑚 = 1 and 2 linear modes for the liquid bridge in Figure 8.5 at the marginal stability. The solid symbols are the results calculated when the free surface contour is assumed to be rigid and the same as that at equilibrium (∆𝑇 = 0). As can be observed, there is no signicant influence of the dynamical free surface deformation on the liquid bridge stability: the growth factors of both the 𝑚 = 1 and 𝑚 = 2 modes are hardly affected by that deformation. In order to extend our analysis to a variety of liquid bridge shapes, we simulated the same experimental run as that analyzed in this section but with different diameters of the supporting rods (aspect ratios) and the same dimensionless volume 𝒱 = 1 (neglecting the increase of volume due to thermal expansion). Figure 8.9 shows the spatial dependence of both the temperature field perturbation and the amplitude of the free surface oscillation for marginally stable liquid bridges. The location of the temperature perturbation maximum moves away from the upper rod as the aspect ratio decreases for 𝐴𝑅≥0.75 (lower graphs). However, the height of that point exhibits a non-monotonous behavior for 𝐴𝑅≤0.615. The amplitude of the free surface oscillation in the liquid bridge upper part increases with 𝐴𝑅, i.e., with the deformation of the base flow free surface with respect to the cylinder. In fact, the upper part of the free surface oscillates more than the lower part for welldeformed equilibrium shapes (large values of 𝐴𝑅). Finally, we verified whether the dynamical free surface deformation affects the thermal convection stability for different aspect ratios. Figure 8.10 shows the critical Marangoni number and oscillation frequency obtained with both the dynamical free surface deformation and the rigid free surface approximation. As can be observed, there is little influence of the dynamical free surface deformation on the liquid bridge stability for the all the aspect ratios considered. 128 Global stability analysis of two-fluid flows Figure 8.9 Spatial dependence of the temperature field perturbation and free surface oscillation amplitude for marginally stable liquid bridges with different aspect ratios. The calculations were conducted for 5-cSt silicone oil liquid bridges with 𝐿 = 3.69 mm and 𝒱 = 1. The free surface oscillation amplitude has been magnified to show clearly its axial dependence. Figure 8.10 Critical Marangoni number 𝑀𝑎 (circles) and oscillation frequency 𝜔𝑟(triangles) as a function of the aspect ratio 𝐴𝑅 for 5-cSt silicone oil liquid bridges with 𝐿 = 3.69 mm and 𝒱 = 1. The results were calculated with the dynamical free surface deformation (open symbols) and the rigid free surface approximation (solid symbols). Influence of the dynamical free surface deformation on the stability of thermal convection in high-prandtl-number liquid bridges 129 8.5. Conclusions In this investigation, we have studied numerically the thermal convection in liquid bridges with high Prandtl numbers. Both the base flow and its eigenmodes have been calculated accounting for the effect of the dynamical free surface deformation. This includes not only the steady axisymmetric deformation caused by the base flow, but also the small-amplitude non-axisymmetric oscillations produced by the growth/decay of the linear modes. We have restricted ourselves to the parameter space corresponding to 5-cSt silicone oil liquid bridges due to its popularity among experimentalists. The comparison between the linear stability analysis predictions and the experimental results measured on board of the ISS (Nishino et al., [222]) shows reasonable agreement, similar to that exhibited by the numerical calculations of Nishino et al. [222]. Our efforts have mainly focused on reproducing a representative experimental run among those conducted in Ref. (Montanero et al., [205]), where the dynamical free surface deformation was accurately measured. We find good agreement between the numerical and experimental results for both the dynamical free surface deformation and the critical Marangoni number. If the deformations caused by the base flows are normalized with the corresponding values of the Capillary number, they approximately collapse onto a single curve. This means that the dynamical free surface deformation can be seen as a steady perturbation with respect to the equilibrium shape, which increases linearly with the Capillary number. We have described the complex spatial structure of the dominant eigenmode responsible for instability. The amplitude of the free surface oscillation is very similar to that measured by Ferrera et al. [206]. We have described the dependence of both the temperature field perturbation and free surface oscillation amplitude upon the liquid bridge aspect ratio. One of the major conclusions of the present work is the fact that the dynamical free surface deformation has little effect on the eigenvalues characterizing the linear modes, and, therefore, on the thermal convection stability. Considering the liquid bridge free surface as a rigid boundary constitutes an accurate approximation. This conclusion must not be necessarily expected given the complex structure of the base flow next to the free surface, especially in the vicinity of the upper hot rod. 131 BIBLIOGRAPHY [1] M. A. Herrada and J. M. Montanero, “A numerical method to study the dynamics of capillary fluid systems,” J. Comput. Phys., vol. 306, pp. 137–147, Feb. 2016. [2] J.-C. Tsai et al., “Vortex-induced morphology on a two-fluid interface and the transitions,” Phys. Rev. E, vol. 92, no. 3, p. 031002, Sep. 2015. [3] S. Fujimoto and Y. Takeda, “Topology changes of the interface between two immiscible liquid layers by a rotating lid,” Phys. Rev. E, vol. 80, no. 1, p. 015304, Jul. 2009. [4] C. E. Brennen, Fundamentals of multiphase flow. Cambridge: Cambridge University Press, 2005. [5] C. T. Crowe, J. D. Schwarzkopf, M. Sommerfeld, and Y. Tsuji, Multiphase Flows with droplets and particles, Second Ed. Boca Raton: CRC Press, 2012. [6] Y. Taitel and D. Barnea, “Two-Phase Slug Flow,” Adv. Heat Transf., vol. 20, no. C, pp. 83– 132, 1990. [7] G. F. Hewitt and N. S. Hall-Taylor, Annular Two-phase Flow, First Ed. Edinburgh: Pergamon Press, 1970. [8] W. D. Bachalo, “Experimental methods in multiphase flows,” Int. J. Multiph. Flow, vol. 20, no. SUPPL. 1, pp. 261–295, 1994. [9] O. Shoham and Society of Petroleum Engineers of AIME., Mechanistic modeling of gasliquid two-phase flow in pipes. Society of Petroleum Engineers, 2006. [10] A. Prosperetti and G. Tryggvason, Computational Methods of Multiphase Flow. New York: Cambridge University Press, 2007. [11] S. Chandrasekhar, Hydrodynamic and hydromagnetic stability. New York: Dover Publications, 2013. [12] W. Criminale, T. Jackson, R. Joslin, and R. Joslin, Theory and computation of hydrodynamic stability. Cambridge: Cambridge University Press, 2003. [13] P. Drazin and W. Reid, Hydrodynamic stability, 2nd edn. Cambridge: Cambridge University Press, 2004. [14] M. Ishii and T. Hibiki, Thermo-fluid dynamics of two-phase flow. Springer, 2011. [15] M. P. Juniper, A. Hanifi, and V. Theofilis, “Modal Stability Theory,” Appl. Mech. Rev., vol. 66, no. 2, p. 024804, Mar. 2014. [16] H. F. FASEL, U. RIST, and U. KONZELMANN, “Numerical investigation of the threedimensional development in boundary-layer transition,” AIAA J., vol. 28, no. 1, pp. 29–37, Jan. 1990. [17] L. Kleiser and T. A. Zang, “Numerical Simulation of Transition in Wall-Bounded Shear Flows,” Annu. Rev. Fluid Mech., vol. 23, no. 1, pp. 495–537, Jan. 1991. [18] F. Charru and P. de Forcrand-Millard, Hydrodynamic Instabilities. Cambridge: Cambridge University Press, 2011. [19] Lord Rayleigh, “On The Instability Of Jets,” Proc. London Math. Soc., vol. s1-10, no. 1, pp. 4–13, Nov. 1878. [20] P. Huerre, “Spatio-Temporal Instabilities in Closed and Open Flows,” in Instabilities and