Thermocapillary Bubble Oscillations and Migration in a Vibrating Cylinder in a Zero-Gravity Environment

Bubble migration in a vibrating zero gravity environment is numerically investigated using ANSYS-FLUENT software. A 3D CFD model is developed describing the two-phase flow of a nitrogen bubble immersed in a container full of ethanol. The Volume of Fluid (VOF) method and the geometric reconstruction scheme are used to track the gas–liquid interface. The liquid in the container is vibrated horizontally in the x momentum with different frequencies from 0 to 1 Hz, and amplitudes from 0.005 to 0.1 m/s2. The vibration impact on the bubble arrival time to the top and its ensuing dynamic is analyzed. Different bubble trajectory shapes are observed, other than the conventional vertical translation induced by the temperature difference. Compared to the no vibration case, the bubble motion is slightly either accelerated or decelerated for very low vibration amplitudes. For a fixed frequency f = 1 Hz, the bubble takes more time to reach the top with the vibration amplitude increment relatively to the no vibration case. The vibration effect becomes more intense with the Marangoni number decrease when f = 0.2 Hz and Ab = 0.005 m/s2. Those results are difficult to obtain experimentally, signifying the importance of this numerical study to understand bubble motion and migration in space.


Introduction
Bubble motion inside a cavity has attracted researcher's attention for a long time, due to its importance in many industrial applications. From petroleum engineering to space, the formation of bubbles can sometimes be harmful for the enclosing equipment (cavitation issues) and beneficial when it is adopted to mass transfer operations (Thompson et al. 1980;Dijkink et al. 2006;Nagasawa et al. 2001). When gravity is taken into account, bubble motion is complicated owing to bubble deformation, coalescence and break up. Many studies were performed regarding the dynamic behavior of a bubble and multiple bubbles (Krishna and Van Baten 1999;Radulescu and Robinson 2008;Yu et al. 2011;Nie et al. 2015). By reason of buoyancy, thermocapillary force is negligible. In order to understand the influence of this force, previous studies investigated bubble motion in a microgravity environment including zero gravity conditions. Young et al. (1959) were the first to prove experimentally that a negative temperature gradient in the vertical direction in pure liquids is sufficient to hold small bubbles stationary or drive them downwards. Thus, the interfacial tension at the bubble surface will change with position because of its dependence on the temperature difference. This phenomenon is the thermocapillary migration. Subramanian (1992) discussed theoretically and experimentally thermocapillary motion of bubbles and drops. He confirmed the results obtained by Young et al. (1959) and showed that under certain conditions, a drop's pure thermocapilllary migration normal to a plane can be more rapid when the drop is near the surface than when it is far. Balasubramaniam et al. (1996) conducted experiments on isolated drop and bubble motion in a reduced gravity environment aboard the NASA Space Shuttle in orbit. They found quantitative discrepancies between the experimental results and the theoretical model of the migration velocity of air bubbles, but the qualitative tendency has 1 3 22 Page 2 of 13 been validated. The interaction between two drops was also examined. A small leading drop, with its motion unaffected, influences significantly the motion of a larger trailing drop by retarding it on account of the thermal wake behind the leading drop. Furthermore, other experiments were carried out to extend the parameter range of this last experiment, Hadland et al. (1999), Wozniak et al. (2001). In Wozniak et al. (2001), the disturbance of the imposed temperature field through bubble motion is studied and compared to a theoretical model based on the numerical model developed by Ma et al. (1999), where the steady migration of a spherical drop is simulated in a liquid under zero gravity conditions. As a consequence of the complexity of the experiments and the disturbances that can occur, the results can sometimes be difficult to interpret. In order to circumvent those issues, many numerical models have been developed as a support for the onboard researches. Herrmann et al. (2008) presented a numerical method to predict the thermocapillary motion of drops. It was found that the drop velocity decreases with Marangoni number (Ma) increase and the drops present a complex behavior at large Ma as observed by Wozniak et al. (2001). Using computational fluid dynamics (CFD) approach, Alhendal et al. (2013), analyzed a single bubble motion in stagnant fluid in 2D and 3D domains. Their results were in good agreement with the space experiments and showed that Volume of Fluid (VOF) method can be profitably employed to simulate thermocapillary flow.
Furthermore, after studying the bubble motion, other parameters were taken into consideration as the bubble size, shape and deformation (O'Shaughnessy and Robinson 2008;Colin et al. 2008;Nurse et al. 2013;Kalendar et al. 2021;Wang et al. 2004), the interaction between droplets (Alhendal et al. 2016a;Balcazar et al. 2016;Kalichetty et al. 2019) and the drop behavior under rotation, (Yamagushi, et al. 2004;Gupta and Kumar 2007;Alhendal et al. (2016b). Another interesting factor is the influence of vibration on drop motion and behavior. It is well known that in spacecraft, while floating in a microgravity environment, vibrations occur classified as g-jitter. Managing those vibrations in space is crucial as an important tool to control flows and heat/mass transfers [Kawaji et al. (2006), Garrabos et al. (2007), and Ahadi and Saghir (2013)]. It is absolutely beneficial for different applications like crystal growth, fluid management, heat exchangers and multiphase flow. Therefore, even if theoretical studies were performed such as Bleich (1956), there is a still a lack of knowledge and understanding on vibration effects in a zero gravity environment. This complicated phenomenon intrigued many researchers and became a very focused topic. Movassat et al. (2009) studied numerically the interaction of two bubbles in a square vibrating container using the VOF solver. Their results showed that bubble collision occurs faster when the vibration frequency and amplitude are important. For high accelerations, bubble circular shape deforms, while only small shape oscillations are observed for small accelerations. At larger accelerations, bubbles breakup before their collision. Shoikhedbrod (2016) developed a model to prove that controlled vibration can be used to control gas bubble motion in the vertical direction and determine conditions of gas bubbles floating, drowning and oscillations in the liquid in reduced and microgravity environments. The theoretical and numerical results obtained were in good agreement with the parabolic aircraft tests.
By simulating complex behaviors such as Marangoni flows in zero gravity conditions using the CFD approach, we can in principle have a better understanding of the ensuing physical phenomena that occur in space and can subsequently manage it in the most optimum way possible. In this work, the impact of vibration on the bubble motion and migration is investigated in a three-dimensional domain (3D) using ANSYS-FLUENT software.

Mathematical Model
The migration velocity or the YGB model, based on the linear model of Young et al. (1959), is used to define the bubble velocity: It is suitable for small Reynolds and Marangoni numbers defined respectively as: The Prandtl number is defined as the ratio of kinematic viscosity, , and thermal diffusivity, .
The velocity V T derived from the tangential stress balance at the free surface is used to scale the migration velocity where and ′ , and ′ are the dynamic viscosity and thermal conductivity of continuous phase and gas, respectively. is the density of the continuous phase fluid and r b is the radius of the bubble. The constant d ∕dT or T is the rate of change of interfacial tension and dT∕dx is the temperature gradient imposed in the continuous phase fluid. (1)

VOF Model and Computational Procedure
In this study, ANSYS-FLUENT package is employed to solve the governing continuum conversation equations. It has been shown in Alhendal et al. (2013), that the use of volume of fluid (VOF) method is very suitable to track the liquid/ gas interface and solve two-phase problems. This method deals with completely separated phases without diffusion. The geometric construction scheme based on the piece-wise linear interface calculation (PLIC) method of Youngs (1982) is chosen. Geo-reconstruction is an added module to the already existing VOF scheme that allows for a more accurate definition of the free surface, Hirt and Nichols (1981). The motion of the bubble-liquid interface is tracked based on the distribution of the bubble volume fraction, i.e. α G , in a computational cell, where the value of α G is 0 for the liquid phase and 1 for the bubble phase. Therefore, the bubble -liquid interface exists in the cell where α G lies between 0 and 1. Throughout the domain, a single momentum equation is solved and shared by all the phases, given below: where � ⃗ v is treated as the mass-averaged variable: �⃗ F vib is the vibration force expressed as: where � ⃗ A b (A b ,0,0) and f are the amplitude and the frequency of the bubble vibration, respectively. �⃗ F represents volumetric forces at the interface, resulting from the surface tension force per unit volume. The continuum surface force (CSF) model proposed by Brackbill et al. (1992) is used to compute the surface tension force for the cells containing the bubble-host liquid interface: where σ is the coefficient of surface tension, and 0 is the surface tension at a reference temperature T 0 , T is the liquid temperature, � ⃗ n is the surface normal which is estimated from the gradient of volume fraction, and k is the local surface curvature, given as: The tracking of the interface between the bubble and hostliquid is accomplished by the solution of a continuity equation for the volume fraction of bubble: The volume fraction equation is not solved for the hostliquid; rather, the liquid volume fraction is computed based on the constraint: where G and L are the volume fraction of the bubble and host-liquid phases respectively. The transport equation properties are determined by the presence of the component phases in each control volume and are calculated as volume-averaged values. The density and viscosity of each cell at the interface is computed by the application of the following equations: where , and G are the density, viscosity and volume fraction of the different phases respectively.
The energy equation is also shared among the phases: The VOF model treats energy ( E ) and temperature ( T ) as mass-averaged variables: where Eq. 16, the energy for each phase is based on the specific heat of that phase and the shared temperature. The properties and k eff (effective thermal conductivity) are also shared by the phases. j is the number of phases.

Model Setup, Boundary and Operating Conditions
The thermocapillary flow of a single bubble in a confined space was the first multiphase test case which would explain and verify the Marangoni phenomenon. After that, the liquid in the container is vibrating in x momentum with different and frequencies (f). The following is a description of the geometry (Fig. 1) and the assumptions made for 2D axisymmetric and 3D numerical simulations: 1. The width and the height of the domain are 60 mm and 120 mm, respectively.
2. Nitrogen bubble, with a diameter d = 2r b , was placed 10 mm from the bottom (cold wall) inside ethanol (Pr = 16.2) using the region adaptation setting on ANSYS-FLUENT. 3. The thermocapillary velocity is small and the flow is laminar. 4. The upper and lower surfaces are flat and non-deformable; the lateral surfaces are adiabatic and non-slip wall conditions is applied to all surfaces. 5. A steady-state temperature distribution is established as an initial condition before releasing the bubble into the unsteady motion, where the top and bottom walls are maintained at constant temperatures: T Top > T Bottom . 6. The host liquid is an incompressible Newtonian fluid, and an assumption of constant properties is applicable, except for surface tension. 7. The driving force for the flow is the variation of surface tension according to temperature, which is modelled by a linear function shown in Eq. (10). 8. For simulations, the properties of nitrogen and ethanol were taken as those given in Table 1 from Thompson et al. (1980).   9. All simulations were run at time t = 0 with an initial stationary liquid and gas, with no velocity at the inlet or the outlet, and the pressure will be taken to equal atmospheric pressure.
After creating the geometry and applying the Quad-map grid, the structured mesh is then transported into the ANSYS-FLUENT software. The calculations are carried out using the pressure-velocity formulation embodied in the PISO  algorithm (Pressure-Implicit with Splitting of Operators) that performs two corrections on the grid quality related to cell neighbor and skewness. The pressure interpolation is achieved using a pressure-staggering option (PRESTO) scheme, while the momentum and the energy equations were discretised using a second-order upwind differencing scheme. Other algorithms were also tried, such as QUICK, instead of the second-order upwind scheme. Using these alternative methods gave similar simulation results. However, using noniterative methods is generally faster than the use of iterative alternatives regarding computational time. Based on the simulation employing different operational conditions, convergence is generally obtained using a non-iterative time step of 5 × 10 −3 s. All simulation runs are executed using doubleprecision accuracy and no gravitational force is imposed.

Grid Dependency and Model Validation
The results obtained in Fig. 2a were in good agreement with those found by Thompson et al. (1980) and Alhendal et al. (2017), and stresses that the VOF Geo-Reconstruct volume fraction model with UDF of surface tension coefficient as a function of temperature is well-coded and validated correctly, suggesting that it is a suitable option for solving thermocapillary problems (Alhendal et al. (2013(Alhendal et al. ( , 2016a). As the accuracy of the simulation is mostly dependent on mesh density, the process of using different mesh sizes, time steps, convergence criteria and discretization schemes, grid tests and extending the geometry to a fully three-dimensional model was all checked (Fig. 2b) and grid independency was achieved by studying five different meshes as presented in Fig. 3a and Table 2. It has been found that for 3D simulation, a grid of 324,000 nodes corresponding to 576 cells per bubble diameter is the most suitable mesh that gives good results with less computing time and memory usage. Figure 3  performed for a better understanding the physical processes underlying many physical phenomena observed in zero gravity environments, the most important of which are vibrations, as Fig. 3c shows more clearly the effect of vibration on the bubble and fluid flow around the bubble, which we will explain in detail in this paper.

Effect of Liquid Temperature Upon Bubble Migration
Before introducing vibration in the domain, it is interesting to first show what happens to the thermocapillary bubbles when a linear temperature distribution is prescribed between the upper and lower walls in 2D axisymmetric before moving to 3D. Four different temperature differences were tested. The lower wall temperature was kept constant at 300 K, and the temperature of the upper wall was varied from 315 K, 317.5 K, 320 K, and 325 K. Figures 4a show the streamlines (top) and isothermal lines contours (bottom) for a single gas bubble N 2 (d = 10 mm) migrating in Ethanol (Pr = 16.2) at t = 9 s. These figures show the bubble's movements toward the warmer side when subjected to a temperature gradient in a zero gravity environment. Such a phenomenon is known as the Marangoni problem or the thermocapillary migration problem. Surface tension generally decreases with increasing temperature and the non-uniform surface tension at the fluid interface leads to shear stresses that act on the outer fluid by viscous forces. This causes bubbles in the fluid to move in the direction of the thermal gradient. Temperature gradients cause surface tension gradients at the liquid-gas interface, and the variation of surface or interfacial tension along the meniscus gives rise to convection. The flow from a region of low surface tension to a region of higher surface tension is referred to as Marangoni flow. This phenomenon is clearest in Fig. 4a. The motion of the bubble triggers motion of the liquid around it. This motion of the ethanol reflects off the bottom of our domain and then rebounds into the back of the  bubble. The effect of the temperature gradient is clear in these figures; for example, as the temperature gradient increases, the migration velocity increases. Figure 4b presents the relationship between migration time and migration distance for different upper wall temperatures. Increasing the temperature gradient (upper wall temperature) simultaneously increases the bubble displacement. The results show that different temperature gradients lead to different bubble migration velocities, indicating a direct relationship between temperature gradient and bubble speed. It also reveals that the bubble velocity decreases if the temperature gradient between the upper and the lower walls decreases, and vice versa.

Results and Discussion
On earth, low frequencies 0 Hz ≤ f ≤ 1 Hz are usually neglected. However, in space this range of vibration frequency can create oscillatory disturbances. At first, the vibration amplitude is fixed at A b = 0.005 m/s 2 and the frequency impact on the bubble migration is discussed. The temperatures of the bottom and the top walls are held at T cold = 300 K and T hot = 325 K respectively (Ma T = 488). Then, for f = 1 Hz, the vibration amplitude is varied ranging from 0.005 to 0.1 m/s 2 and its effect is analyzed on the bubble dynamics.

Vibration Frequency
For Ma T = 488 and A b = 0.005 m/s 2 , the effect of vibration frequency on the bubble motion is plotted in Fig. 5. When the vibration is absent, the bubble moves in a vertical trajectory from the container bottom to its top. At f = 0.01 Hz, a vertical translation of the bubble is observed at first. Then, in a position, almost equal to the half of the container height, the bubble motion deviates slightly to the left and continues deviating until it arrives at the top. A similar behavior is observed for f = 0.05 Hz with an enhancement of the trajectory deviation from a position before the half of the container height, near the bottom. The bubble arrives at the top in a position far from the center. By increasing f to 0.1 Hz, the bubble oscillates creating a trajectory in the shape of a snake and reaches the top almost in the center. For f = 0.175 Hz, The bubble's path becomes more wavy and it returns to the center as it reaches the top. When the frequency increases to f = 1 Hz, the bubble motion seems to be in a vertical translation; however, this is not the case. The bubble oscillates periodically around the container vertical axis until it arrives at the top. Figure 6 shows the bubble center position in the three directions. For no vibration case and all the frequencies, the bubble travels from the bottom to the top (Y direction) similarly in the first 2.5 s, then, beyond 2.5 s, the bubble arrives rapidly at the container top when f = 0.175 Hz and f = 1 Hz taking 10.75 s (Fig. 6a). This time is reduced by 0.25 s, 0.5 s and 0.75 s vis-à-vis the no vibration case, f = 0.01 Hz, and both f = 0.05 Hz and f = 0.1 Hz respectively. This slight difference in the arrival time can be explained by the oscillatory  (Fig. 6b and c). Those pulsations deform the bubble shape from its initial circular shape as observed by Movassat et al. (2009) as shown previously in Fig. 3c. In general, the primary force influencing the motion is thermocapillary, but choosing the right frequency can accelerate the bubble migration. For f = 0.01 Hz, f = 0.05 Hz and f = 0.1 Hz, the bubble migration is delayed compared to the no vibration case. The deviation of the bubble from the center reduced the bubble motion.

Vibration Amplitude
In this section, for Ma T = 488, the vibration frequency is fixed at f = 1 Hz and its amplitude, A b , is varied ranging from 0.005 m/s 2 to 0.1 m/s 2 . Actually, the bubble motion becomes complicated with A b increase, which consequently made it difficult to present clear plots of its ensuing behavior. For this reason, from the bubble center positions in X, Y and Z directions, the global motion can be understood (Fig. 7). Understandably, the bubble takes a longer time to arrive at the top when the amplitude increases. From 11 s in the no vibration case to 38.5 s when A b = 0.1 m/s 2 (Fig. 7a). With the amplitude increment, the oscillations become intense in both X and Z directions and change the bubble shape that slows down the bubble motion. The bubble oscillates around the Z axis in the first 10 s for all the amplitudes, then deviates (Fig. 7b). Regarding the X direction, the oscillations are noticeable from A b = 0.05 m/s 2 . For A b = 0.06 m/ s 2 , the bubble oscillations are quite symmetric on the negative side related to the X axis before reaching the top. Beyond A b = 0.06 m/s 2 , those oscillations become asymmetric and drift, displaying a transient behavior (Figs. 7c, d). The second oscillation begins with different velocity and shape than the first oscillation and the same thing behavior can be observed for the other oscillations. Additionally, similar dynamics have been observed for high frequencies and amplitudes by Movassat et al. (2009). Compared to the frequency variation effects, the amplitude influences considerably not only the bubble trajectory but also the arrival time at the top. In fact, the bubble migration decelerates or accelerates with the frequency increase as compared to a negligible difference vis-à-vis the no vibration case. However, it always decelerates with amplitude augmentation. The influence affecting the shape is greater at high amplitudes caused by the changes that also occur regarding the inertia and pressure forces on the bubble. Therefore, the bubble displays a non-linear behavior while incorporating vibration influences on the thermocapillary force.

Vibration Impact on Different Ma T
For f = 0.2 Hz and A b = 0.005 m/s 2 , the temperature of the top is changed from 310 K, 315 K, 320 K, 325, and 330 K corresponding to Ma T = 162.7, 244, 325.3, 406.7, and 488 respectively. Figure 8 shows different bubble trajectory shapes. The bubble oscillates and its shape deforms while traveling from the bottom to the top. However, in the no vibration case, regardless of Ma T , the bubble motion keeps a vertical linear trajectory. Furthermore, the bubble takes less time to reach the top whenever the Marangoni number is important (Fig. 9a). Actually, this finding is similar to the results found by Alhendal et al. (2013) for the no vibration case. Oscillations are also noticed in the Z direction (Fig. 9b). In the X direction, the bubble is in a translatory motion for the first 2.5 s when Ma T = 162.7. Then, oscillations appear drifting to the positive side when the bubble approaches the container top (Fig. 9c). For Ma T = 244 and Ma T = 325.35, the first oscillation is in x negative direction then becomes x positive direction, drifting up and finally pushing the bubble to the positive side of the container wall. Beyond those Marangoni numbers, the oscillations are asymmetrical around the X axis and pushing the bubble near the container center. In fact, the undulations, at Ma T = 162.7, are very noticeable compared to the other cases. In this studied Marangoni range, those undulations diminish with Ma T increase. Therefore, when the temperature difference rises, the thermocapillary force is less affected by vibration.

Conclusions
In a container full of ethanol, a nitrogen bubble migration under vibration, in zero gravity environment, is studied using ANSYS-FLUENT. The governing equations of the simulated model are discretized using the finite volume method and solved by the PISO algorithm. The VOF method with the geometric reconstruction scheme is employed to model the two-phase flow. A small change in the amplitude or the frequency vibration, the bubble motion changes from vertical translation to oscillatory trajectories and the bubble shape deforms as well. The bubble becomes elliptical instead of spherical. This shape variation implies the unsteady motion of the bubble in the vibrating fluid surrounding it that, at his turn, changes the bubbles orientation leading to the flow patterns observed. For a fixed amplitude and regardless of the frequency, the bubble arrival time is slightly similar to the no vibration case, despite the changes observed on the trajectory. However, when the amplitude increases, asymmetric and drift oscillations are observed making the bubble migration slower vis-à-vis the no vibration case. In general, the motion leading force is the thermocapillary force, but choosing the adequate frequency and amplitude plays an important role in either accelerating or slowing down the bubble migration. Similar to the no vibration case, increasing the Marangoni number accelerates the bubble migration. Those novel results, which are difficult to obtain by experiments, can help researchers to fathom out the thermocapillary flow and bubble motion in space and therefore, enable practices to exploit them in the best way possible in different applications as fluid management and multiphase flows.