## Investigation of Surfactant-Laden Bubble Migration Dynamics in Self-Rewetting Fluids using Lattice Boltzmann Method

Bashir Elbousefi, William Schupbach, Kannan N. Premnath, Samuel W.J. Welch ∗ Department of Mechanical Engineering University of Colorado Denver 1200 Larimer Street, Denver, CO 80204, U.S.A.

October 10, 2024

## Abstract

Self-rewetting fluids (SRFs), such as aqueous solutions of long-chain alcohols, show anomalous nonlinear (quadratic) variations of surface tension with temperature involving a positive gradient in certain ranges, leading to different thermocapillary convection compared to normal fluids (NFs). They have recently been used for enhancing thermal transport, especially in microfluidics and microgravity applications. Moreover, surfaceactive materials or surfactants can significantly alter interfacial dynamics by their adsorption on fluid interfaces. The coupled effects of temperature- and surfactant-induced Marangoni stresses, which arise due to surface tension gradients, on migration bubbles in SRFs remain unexplored. We use a robust lattice Boltzmann (LB) method based on central moments to simulate the two-fluid motions, capture interfaces, and compute the transport of energy and surfactant concentration fields, and systematically study the surfactant-laden bubble dynamics in SRFs. When compared to motion of bubbles in normal fluids, in which they continuously migrate without a stationary behavior, our results show that they exhibit dramatically different characteristics in SRFs in many different ways. Not only is the bubble motion directed towards the minimum temperature location in SRFs, but, more importantly, the bubble attains an equilibrium position. In the absence of surfactants, such an equilibrium position arises at the minimum reference temperature occurring at the center of the domain. The addition of surfactants moves the equilibrium location further upstream, which is controlled by the magnitude of the Gibbs elasticity parameter that determines the magnitude of the surface tension variation with surfactant concentration. The parabolic dependence of surface tension in SRF is parameterized by a quadratic sensitivity coefficient, which modulates this behavior. The lower this quantity, the greater is the role of surfactants modifying the equilibrium position of the bubble in SRF. Furthermore, the streamwise gradient in the surfactant concentration field influences the transient characteristics in approaching the terminal state of the bubble. These findings provide new means to potentially manipulate the bubble dynamics, and especially to tune its equilibrium states, in microchannels and other applications by exploiting the interplay between surfactants and SRFs.

## 1 Introduction

Dispersed two-phase flows, such as those involving bubbles or drops, are common in nature, everyday life, and industrial applications, and represent a topic of major research interests. Because surface tension forces can alter how fluids behave at the interface, they are

*Corresponding author (Email: Bashir.Elbousefi@ucdenver.edu)

important in multiphase and heat transport processes [1] and, in particular, variable surface tension flows have found various applications related to microfluidics and microgravity transport phenomena. Such technological applications rely on Marangoni stresses [2], which are caused by variations in the surface tension (i.e., the surface tension gradients) at the interface between two immiscible fluids and arise due to variations in the local interfacial temperature or due to the use of surfactants or due to the coupling between the two. These stresses cause convective motions to occur close to the interfaces due to the viscous effects of the fluids [3]. Surface tension driven flows due to temperature, or the thermocapillary convection, and the associated migration of bubbles or drops in fluids have been the subject of many studies over the years (see e.g., [4-8]). In the case of micro-electro-mechanicalsystems, the thermocapillary convection is utilized to manipulate fluid streams and thermal transport phenomena, as well as the motion of bubbles and drops in microchannels as interfacial forces become dominant at the small device scales (see, e.g., [9,10]).

Surface tension is a common property of fluids that decreases with increasing temperatures in most fluids. However, some fluids exhibit anomalous behavior in which the surface tension deviates from this linear relationship and instead exhibits a non-linear parabolic dependence on temperature with a range involving a positive gradient. Aqueous solutions of high carbon content-based alcohols (e.g., 1-butanol, 1-pentanol, and 1-heptanol) exhibit an increase in surface tension as the temperature increases in certain operating ranges. Over the years, various experiments on such special classes of fluids conducted by different researchers have measured and demonstrated this behavior [11-14]. It may be noted that in addition to mixtures of certain alcohols with water, some liquid metallic alloys and nematic liquid crystals also exhibit anomalous behavior in which the surface tension decreases with temperature up to a certain point, after which it begins to increase. Abe et al. [15] named these fluids 'self-rewetting' fluids (SRFs). They experimentally revealed that SRFs have much potential for transferring heat because of how easily they can change phases, which, compared to common or normal fluids (NFs), exhibit a significantly altered thermocapillary convection that promotes a desired wetting effect especially over hotspots.

Marangoni stresses, in particular, cause fluids near the interfaces to move towards higher temperatures (where the fluids have higher surface tensile strength) in SRFs (as opposed to NFs). For example, the schematic diagram in Fig. 1 shows how the two kinds of fluids behave differently with temperature and the resulting Marangoni stress in SRF around a bubble in motion. In self-rewetting fluids (SRFs), the surface tension increases with temperature.

Fig. 1. ( a ) Surface tension variation with respect to temperature T and surfactant concentration ψ for a SRF ( M 2 = 5) and ( ) for a NF ( b M 1 = -2 5). . ( c ) Schematic diagram shows the directions of both surfactant-induced and thermal-induced Marangoni stresses for SRF acting on the bubble surface. See Eq. (10) for definitions of M 1 and M 2 .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00146v2-figure-caf02607-4616-4e43-bddc-dc66f7a6902e.jpg)

As a result, when a temperature gradient is present at the interface of two fluids, the SRF flows from the cold to the hot region, in contrast to what occurs in normal fluids (NFs), where the fluid flows from hot to cold regions due to the decrease in surface tension with increasing temperature [16-18]. Effectively, this causes fluid currents towards higher temperature zones in SRF, which can be exploited to enhance transport. As such, SRFs have been used as working fluids in a variety of thermal management applications in both terrestrial and microgravity environments [15,19], such as in heat pipes [20-30], flow boiling [31] and evaporation [32] in microchannels, pool boiling processes [33-37], and two-phase heat transfer devices [38]. Moreover, some experiments [39, 40] have recently studied how bubble migration in SRFs has certain unique characteristics unlike what has been observed in common fluids.

Apart from temperature variations, surfactants can also significantly modify the local surface tension in fluid interfaces. Surfactants, also known as surface-active substances, are molecules that preferentially absorb at the interface between two phases. They can be soluble or insoluble in the bulk fluids and their net adsorption/desorption characteristics

on interfaces effectively cause a local reduction in surface tension as well generating its tangential gradients or the Marangoni stresses. The latter arises due to a nonuniform concentration of surfactants on the interfaces. For example, a bubble in motion involves an adsorption at its trailing edge, while a part of its front region remains free of surfactants and the resulting Marangoni stress causes the bubble to move towards areas with lower surface tension. See the schematic Fig. 1c illustrating this effect. Surfactants are used widely in many areas related to technical, chemical, and biological applications [41], including in microfluidic systems [42-46].

Multiphase flows arise in a wide range of complex situations (see e.g., [47,48]). To provide insights into the underlying physics associated with the transport processes in such cases, computational studies that are based on appropriate physical models are essential [49,50]. It would be interesting to investigate the coupled effects of both thermal- and surfactantinduced Marangoni stresses on bubble migration in self-rewetting fluids (SRFs). However, much of the prior research on SRFs in this regard involving either analytical [51-55] or numerical [56-60] approaches do not include considerations of surfactants, including our recent work on analytical and numerical study on thermocapillary convection of superimposed SRF layers [61]. One notable exception is that of [62], which studied the SRF droplet spreading over a heated film in the presence of surfactants. Nevertheless, most of existing studies that combine thermocapillary effects with surfactant behavior are all directed towards normal fluids (NFs), where the surface tension varies linearly with temperature (see e.g., [45, 63, 64]), and have not explored them for SRFs, especially involving bubble dynamics.

Thus, the aim of this work is to systematically investigate the coupled effects of both thermal- and surfactant-induced Marangoni stresses on bubble migration in self-rewetting fluids (SRFs) in order to provide a more complete physical understanding of the behavior of the dispersed phases in such configurations, which can, in turn, inform novel ways to potentially manipulate them. More specifically, previous studies have examined the effects of thermal-induced Marangoni stresses on bubble migration in SRFs, but have neglected the additional influence of surfactant-induced Marangoni stresses which, as this work will show, leads to certain interesting outcomes for the bubble dynamics. One primary objective is to present and apply a numerical simulation approach based on the lattice Boltzmann method (LBM) in conjunction with a phase field model for capturing of interfaces along with the attendant transport models for energy and surfactant dynamics. The LBM is a numerical approach inspired by the kinetic theory of gases [65-67], and its natural parallelization features, efficient and simpler implementation algorithms, and ability to handle complex geometries and boundary conditions make it a powerful and versatile technique for simulating fluid motions under a variety of situations. As such, the LBM has been found to be a versatile approach to tackle a variety of complex fluid flows and beyond. The recent advances of this method have been discussed in a number of reviews for various applications, such as thermal and aerodynamic flows [68], turbulent and multiphase systems of complex fluids [69], phase change problems including melting and solidification [70], non-ideal fluids [71] and reacting flows [72]. Over the years, the LBM attracted attention especially for simulation of multiphase flows [73-78]. As a result, the LBM has become an important numerical simulation tool for multiphase flows alongside more conventional techniques such as volume of fluid, level set method, and front tracking method in multiphase flows. LB

methods have also been used to model thermocapillary flows (see e.g., [58, 59, 61, 79, 80]). In recent years, new collision models based on central moments have been developed that enhance the robustness of the LBM [81] and allow for the application of the LBM method to simulate high-density ratio multiphase flows that are subject to Marangoni stress [77]. Such collision models have been used in the LBM for simulation of thermocapillary flows in SRFs in our recent work [61], which will be further extended by incorporating surfactant dynamics, including its adsorption/desorption process and bulk solubility effects, and its contribution to the surface tension equation of state.

Thus, in this work, we will use a robust central moment LB approach that involves the computing of four different distribution functions. One is used to compute the two-fluid motion and includes the effect of Marangoni stresses based on surface tension gradients, another one to capture the interfacial motions represented by the conservative Allen-Cahn equation (ACE), and, the third one to compute the energy transport and, finally, the fourth one for the surfactant concentration field. Local surface tension is then modulated by using a parabolic dependence on temperature for self-rewetting fluids coupled with the effect of the surfactant concentration via the Langmuir isotherm. As another key objective, we present the simulation results of surfactant-laden bubble dynamics and its terminal states in SRFs by studying the effect of variations of different characteristic parameters associated with this flow configuration. We will compare and contrast the behavior of bubbles in SRFs with those in NFs with the addition of surfactants.

The remaining sections of the paper are structured as follows. In the upcoming section (Section 2), we will explain the problem setup of the migration of bubbles laden with surfactants through thermocapillary effects in a SRF. We will also discuss the governing equations that describe the incompressible two-fluid motion, the transport of energy and surfactant concentration, and the equation of state for interfacial surface tension. The diffuse-interface computational model equations for the LBM of surfactant-laden self-rewetting fluids are provided in Section 3. We summarize the discretized central moment LB algorithms for simulating multiphase flows of SRFs in Section 4, while Appendices B,C,D, and E discuss the LB schemes for interfacial dynamics, two-fluid motion, energy transport, and surfactant concentration transport, respectively. Section 5 presents a numerical validation of our computational approach by comparing it to established benchmark problems. The results and discussion regarding the influence of various characteristic parameters on the physics of surfactant-laden thermocapillary bubble migration in SRF using LB schemes are presented in Section 6. Finally, we provide an overview of the main conclusions and contributions of our research in Section 7, and additional supporting details can be found in the appendix sections.

## 2 Problem setup, governing equations, and modeling of surface tension effects

The objective of the study is to simulate the surfactant-laden thermocapillary bubble migration in self-rewetting fluid (SRF). The simulations are conducted within a rectangular computational domain with dimensions ( L × H ), where L = 600 and H = 0 5 . L are the lengths of the domain in the x and y directions, respectively. The problem is schematically

represented in Fig. 2, which shows the computational domain and the bubble with a diameter of D = 0 1 . H is initially located at (0 25 . L, 0 5 . H ) where the origin of the coordinate system is chosen to be at the bottom left corner of the domain. The lattice units are used for all the values as is typical for the implementation of the LB method (see e.g. Ref. [82]). No slip boundary conditions are imposed on the left and right walls, while periodic boundary conditions are used on the top and bottom sides. As noted below, the side walls are imposed with different thermal conditions to drive the bubble motion via thermocapillary effects. Since they are generally solid surfaces, the no-slip boundary conditions are appropriate. Moreover, periodic conditions are imposed on the sides to ensure accurate representation of bubble dynamics without boundary interference, which is a common approach used in such dispersed flow simulations. We initialize this problem with a constant temperature gradient applied in the x direction, with a value of | ∇ T | = ( T H -T C ) /L such that a fixed hot temperature on the left side ( T H = 1) and a constant cold temperature, ( T C = 0) on the right side, and then we maintain the reference temperature ( T ref = 0 5) to occur at the . middle of the domain. Moreover, we introduce a background surfactant gradient such that the gradient in the dimensionless surfactant concentration ψ decreases from the left side ( ψ = ψ L ) to the right side ( ψ = ψ R ) (see the paragraph above Eq. (2) for the definition of ψ ). The initial bubble placement and surfactant gradients were chosen to adequately study the interesting features of dynamics of the bubble in the SRF. Given the temperature boundary conditions noted above with the temperature having a minimum in the middle of the domain and increasing on both sides, we placed the bubble halfway between the left side and the middle so that it accommodates the peculiar features of its motion in the SRF under thermocapillary effects. This ensures it sufficiently captures the transient effects before the bubble is expected to equilibrate in the vicinity of the center of the domain. A surfactant gradient is introduced in the background to generate additional Marangoni stresses to manipulate the dynamics and the final equilibrium position of the bubble. For convenience, the various symbols used in this paper are listed in Appendix A.

Fig. 2. Schematic of the initial configuration of the surfactant-laden thermocapillary bubble migration in self-rewetting fluid (SRF).

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00146v2-figure-87130f27-52bf-420d-bae7-6db464cce4cd.jpg)

Here, and in what follows, we use a subscript notation with 'A' for the ambient fluid and 'B' for bubble. The thermocapillary convection in the SRFs including surfactant obeys the equations of mass and momentum (i.e., the Navier-Stokes equations (NSE)), energy transport, and the equation of surfactant concentration transport. They can be respectively

written as follows:

$$\nabla \cdot u = 0,$$

$$\rho \left ( \frac { \partial u } { \partial t } + u \cdot \nabla u \right ) = - \nabla p + \nabla \cdot \left [ \mu ( \nabla u + \nabla u ^ { \dagger } ) \right ] + F _ { e x t },$$

$$\frac { \partial T } { \partial t } + \mathbf u \cdot \mathbf \nabla T = \mathbf \nabla \cdot \left ( \alpha \mathbf \nabla T \right ),$$

$$\frac { \partial \psi } { \partial t } + u \cdot \nabla \psi = \nabla \cdot \left ( M _ { \psi } \nabla \mu _ { \psi } \right ),$$

where ρ , µ , and α are the fluid density, dynamic viscosity, and thermal diffusivity of the fluid, respectively, with α = k/ ρc ( p ) based on the thermal conductivity k and specific heat c p . In the above, u , p , and T denote the velocity, pressure, and temperature fields of the fluids, respectively, and the superscript symbol † represents taking the transpose of the dyadic velocity gradient ∇ u . Also, F ext is any external body force. While the mass, momentum, and energy equations are standard, the modeling of surfactants as represented by Eq. (1d) calls for further discussion. Few different approaches to model surfactants within the phase field framework for its use with the LBM have been proposed (see e.g., [83-87]).

In Eq. (1d), the surfactant concentration ψ is non-dimensionalized by the maximum surfactant concentration on the interface ψ max , which is a material property that determines its interfacial adsorption capacity, such that ψ = ψ/ψ max and M ψ = m ψ ψ (1 -ψ ) is the local surfactant mobility, with m ψ being the scale for the mobility parameter. Hence, the nondimensional local surfactant concentration will be restricted between 0 and 1 (0 ≤ ψ ≤ 1); µ ψ is the chemical potential whose gradients drive the diffusion-adsorption dynamics of the advecting surfactant concentration field and can be written as [83,84]

$$\mu _ { \psi } = \lambda \ln \left ( \frac { \psi } { 1 - \psi } \right ) - \frac { s } { 2 } | \nabla \phi | ^ { 2 } + \frac { w } { 2 } ( \phi - \phi _ { m } ) ^ { 2 }.$$

Here, ϕ is the order parameter of the phase field variable (see Eq. (12) for its governing equation based on a diffuse interface model) that is used to capture the dynamics of interfaces nominally located at ϕ m = ( ϕ A + ϕ B ) / 2, where ϕ A and ϕ B are the values of ϕ in the bulk fluids A and B. The preferred tendency of the surfactant adsorption on the interface is represented by the second term on the right hand side (RHS) of the above chemical potential Eq. (2) represented in terms of a delta function through a square gradient of the order parameter. However, some recent studies [88, 89] suggested replacing it with a gradient-free regularized delta function formulation obtained using the hyperbolic tangent profile across the interface in the normal direction for the phase field variable ϕ in order to achieve well-posedness with better numerical properties for a wider range of parameter choices. By using the latter approach and taking ϕ o = ( ϕ A -ϕ B ) / 2, it is possible to rewrite the chemical potential utilized in this work as

$$\mu _ { \psi } = \lambda \ln \left ( \frac { \psi } { 1 - \psi } \right ) - \frac { s } { 2 } \frac { 4 } { \phi _ { o } ^ { 2 } W ^ { 2 } } \left [ \phi _ { o } ^ { 2 } - ( \phi - \phi _ { m } ) ^ { 2 } \right ] ^ { 2 } + \frac { w } { 2 } ( \phi - \phi _ { m } ) ^ { 2 }.$$

The model parameters λ , s , and w here represent the relative strengths of the several competing processes as explained below.

Each term in Eq. (3) on the RHS of the chemical potential has the following possible interpretation [88, 89]: The first term, known as the entropy term, bounds the surfactant concentration between 0 and 1 and indicates an increase in the system's entropy when the substance is evenly distributed everywhere. Higher concentrations λ drive stronger diffusion, which tends to redistribute the substance more evenly throughout the domain. The parameter s may be used to adjust the strength of the second term, which is also referred to as the adsorption term. It indicates the surfactant's energetic propensity to get adsorbed on interfaces. The final term, also known as the bulk term, is a penalty term that represents the solubility effect and penalizes the presence of surfactants in bulk fluids. Its magnitude is controlled by the parameter w . The surfactant concentration profile surrounding interfaces is sharpened by these latter two terms, which oppose the diffusive process (first term). The surfactant concentration field ψ evolution equation may be further reformulated by substituting Eq. (3) into Eq. (1d) and rearranging to get

$$\frac { \partial \psi } { \partial t } + u \cdot \nabla \psi = \nabla \cdot \left ( \lambda m _ { \psi } \nabla \psi \right ) + \nabla \cdot \left ( m _ { \psi } \psi ( 1 - \psi ) R \right ),$$

where the second term in the RHS of the above equation can be considered as a flux term related to the surfactant's adsorption and solubility, and R is given by the following equation [90]

$$R = \nabla \left ( - \frac { s } { 2 } \frac { 4 } { \phi _ { o } ^ { 2 } W ^ { 2 } } \left [ \phi _ { o } ^ { 2 } - ( \phi - \phi _ { m } ) ^ { 2 } \right ] ^ { 2 } + \frac { w } { 2 } ( \phi - \phi _ { m } ) ^ { 2 } \right ).$$

## 2.1 Surface tension equation of state including surfactant and thermocapillary effects for self-rewetting fluids

In the case of normal fluids (NFs), the effect of inclusion of surfactants was modeled by Luo et al. [45] who reported a form for the surface tension equation of state based on a linear term related to its variation with temperature which is combined with an additional term associated with surfactant effects. Our work extends this for self-rewetting fluids (SRFs) with surfactants. To impose an equation for the surface tension at the interface, we need to relate it to the variations in temperature T and surfactant concentration ψ . To combine both effects, we use a nonlinear (quadratic) dependence of surface tension on temperature for SRF and Langmuir isotherm for surfactant effects which are given in the following equation:

$$\sigma ( T, \psi ) = \sigma _ { 0 } + \sigma _ { 0 } \beta \left ( \frac { T } { T _ { r e f } } \right ) \ln ( 1 - \psi ) + \sigma _ { T } ( T - T _ { r e f } ) + \sigma _ { T T } ( T - T _ { r e f } ) ^ { 2 },$$

where σ 0 denotes the value of the surface tension at a reference temperature T ref which is in the middle of the domain according to the quadratic functional distribution given above; σ T = dσ dT ∣ ∣ T ref and σ TT = 1 2 d σ 2 dT 2 ∣ ∣ T ref are the sensitivity coefficients of the surface tension with linear and quadratic dependencies on temperature, respectively. In Eq. (6), β is the Gibbs

elasticity parameter, which is a dimensionless quantity that determines the influence of the surfactant on the surface tension variations and is given by

$$\beta = \frac { R _ { g } T _ { r e f } \psi _ { m a x } } { \sigma _ { 0 } },$$

where R g is the ideal gas constant and ψ max is the maximum surfactant concentration on the interface. Thus, β models the magnitude of the variations in the surface tension due to the addition of surfactants which is a material property of the latter.

When the above governing equations are nondimensionalized using a reference velocity scale U and a length scale L corresponding to the radius of the bubble R , the additional dimensionless groups that arise in the study of the thermocapillary motion of a bubble in SRF are as follows: The Reynolds number Re, Marangoni number Ma, and capillary number Ca are defined by

$$\text{Re} = \frac { U R } { \nu _ { A } }, \quad \text{Ma} = \frac { U R } { \alpha _ { A } }, \quad \text{and} \quad C \text{a} = \frac { U \mu _ { A } } { \sigma _ { 0 } },$$

respectively, where ν = µ/ρ is the kinematic viscosity of the fluid, and the rest of the symbols are already defined below Eq. (1d). Also, the following ratios of the material properties govern the thermocapillary-driven motion of a bubble in a SRF:

$$\tilde { \rho } = \frac { \rho _ { A } } { \rho _ { B } }, \quad \tilde { \mu } = \frac { \mu _ { A } } { \mu _ { B } }, \quad \tilde { k } = \frac { k _ { A } } { k _ { B } }, \text{ \ and \ } \tilde { c } _ { p } = \frac { c _ { p _ { A } } } { c _ { p _ { B } } }.$$

Moreover, the dimensionless forms of the linear and quadratic sensitivity coefficients of the surface tension σ T and σ TT (see Eq. (6)) can be, respectively, represented as

$$M _ { 1 } = \left ( \frac { \Delta T } { \sigma _ { 0 } } \right ) \sigma _ { T }, \ \ M _ { 2 } = \left ( \frac { \Delta T ^ { 2 } } { \sigma _ { 0 } } \right ) \sigma _ { T T },$$

where ∆ T is the characteristic temperature difference in the domain. Based on the above, we can also rewrite the surface tension equation of state given in Eq. (6) in the following dimensionless form:

$$\frac { \sigma ( T, \psi ) } { \sigma _ { 0 } } = 1 + \beta \left ( \frac { T } { T _ { r e f } } \right ) \ln ( 1 - \psi ) + M _ { 1 } \frac { ( T - T _ { r e f } ) } { \Delta T } + M _ { 2 } \frac { ( T - T _ { r e f } ) ^ { 2 } } { \Delta T ^ { 2 } }.$$

In addition, the non-dimensional numbers involving the surfactant dynamics based on Eqs. (4) and (5) include Pi = λW / 2 (8 σ 0 ), which indicates the proportionate contribution of surfactant diffusion, and Ex = 4 s/ wW ( 2 ), which represents the relative intensity of adsorption and solubility effects where the parameters λ , s , and w , represent the relative strengths of the various competing processes. Finally, We can estimate the scale for the reference velocity U of thermocapillary bubble migration in a fluid used in the above via U o ∼ ( σ /µ 0 a )( R/L )( M 1 + M 2 ) by balancing the scale for the viscous shear stress or µ U/R a with that of the primary Marangoni stress due to the surface tension gradient arising from temperature variations or | dσ/dT | (∆ T/L ) + ( d σ/dT 2 2 )(∆ T /L 2 ).

## 3 Diffuse-interface computational modeling of surfactant-laden self-rewetting fluids

In the upcoming section, we will explore a modeling approach that utilizes the LBM to simulate thermocapillary convection in SRFs. To accurately capture interfacial dynamics and maintain the segregation of two immiscible fluids, we will employ a phase-field lattice Boltzmann approach based on the conservative Allen-Cahn equation (ACE) [91]. This approach is an improvement over an earlier model [92] that utilized a counter-term approach [93]. In this model, the binary fluids are distinguished by an order parameter or the phase field variable ϕ , with fluid A being identified by ϕ = ϕ A and fluid B by ϕ = ϕ B . Ultimately, we will utilize the interface-tracking equation based on the conservative ACE to solve for the phase field variable which is given as

$$\frac { \partial \phi } { \partial t } + \nabla \cdot ( \phi u ) = \nabla \cdot [ M _ { \phi } ( \nabla \phi - \theta n ) ].$$

The fluid velocity is denoted by u , the mobility by M ϕ , and the unit normal vector by n . The unit normal vector may be computed using the order parameter ϕ as follows: n = ∇ ϕ/ | ∇ ϕ | . In this case, θ = -4 ( ϕ -ϕ A ) ( ϕ -ϕ B ) / W [ ( ϕ A -ϕ B )] may be used to define the parameter θ , where W is the interface width. In Eq. (12), the term M θ ϕ n functions as the interface sharpening term which counteracts the diffusive flux -M ϕ ∇ ϕ that occurs after ϕ is advected by the fluid velocity. A hyperbolic tangent profile across the diffuse interface is the order parameter that the conservative ACE reduces to at equilibrium. It can be found as follows: ϕ ζ ( ) = 1 2 ( ϕ A + ϕ B ) + 1 2 ( ϕ A -ϕ B ) tanh (2 ζ/W ), where ζ is a spatial coordinate along the normal with the origin at the interface.

In a single-field formulation depicting the motion of binary fluids, the interfacial surface tension effects may now be easily implemented inside a diffuse interface by the use of a distributed or smoothed volumetric force term. Thus, the equivalent Navier-Stokes equations for binary fluids as a single-field formulation for incompressible flow may be expressed as

$$\nabla \cdot u = 0,$$

$$\rho \left ( \frac { \partial u } { \partial t } + u \cdot \nabla u \right ) = - \nabla p + \nabla \cdot \left [ \mu ( \nabla u + \nabla u ^ { \dagger } ) \right ] + F _ { s } + F _ { e x t },$$

where any exterior body force is denoted by F ext and the surface tension force by F s . Here, as surface tension varies with temperature and surfactant concentration, and the surface tension force effectively acts in both the normal and tangential directions to the interface. The continuous surface force approach [94] is a geometric technique that can be used to model these features. It is represented by the following equation, which uses the Dirac delta function δ s :

$$F _ { s } = \left ( \sigma \kappa n + \nabla _ { s } \sigma \right ) \delta _ { s }.$$

Here, n = ∇ ϕ/ | ∇ ϕ | and κ = ∇ · n are the unit normal vector and interface curvature, respectively. The normal or capillary force operating on the interface is represented by the

first term on the right side of Eq. (15), while the tangential or Marangoni force caused by surface tension gradients is represented by the second term, which involves the surface gradient operator ∇ s . The delta function δ s is necessary to satisfy the property that ∫ + ∞ -∞ δ dy s = 1 since surface tension only affects the interface. Here, δ s = 1 5 . W | ∇ ϕ | 2 is a formulation of δ s that localizes the smoothed surface tension force suitable for use within the phase-field modeling framework. In addition, ∇ s = ∇ -n n ( · ∇ ) is the surface gradient ∇ s in Eq. (15). Consequently, it is possible to define the Cartesian components of the surface tension force in Eq. (15) as

$$F _ { s x } \ & = \ \ - \sigma ( T, \psi ) | \nabla \phi | ^ { 2 } ( \nabla \cdot n ) n _ { x } + | \nabla \phi | ^ { 2 } \left [ ( 1 - n _ { x } ^ { 2 } ) \partial _ { x } \sigma ( T, \psi ) - n _ { x } n _ { y } \partial _ { y } \sigma ( T, \psi ) \right ], \\ F _ { s y } \ & = \ \ - \sigma ( T, \psi ) | \nabla \phi | ^ { 2 } ( \nabla \cdot n ) n _ { y } + | \nabla \phi | ^ { 2 } \left [ ( 1 - n _ { y } ^ { 2 } ) \partial _ { y } \sigma ( T, \psi ) - n _ { x } n _ { y } \partial _ { x } \sigma ( T, \psi ) \right ]. \ ( 1 6 )$$

Here, the combined effect of the nonlinear (parabolic) component and the Langmuir isotherm term, respectively, yield the functional dependence of the surface tension on temperature for the SRF and on surfactant concentration (see Eq. (6) for σ = σ T, ψ ( )). This study uses an isotropic finite differencing approach [95] to determine the needed spatial gradients ∂ σ T, ψ x ( ) and ∂ σ T, ψ y ( ) in Eq. (16) for numerical implementations. The temperature field T and the surfactant concentration field ψ are determined by solving the equations for the energy transport and surfactant concentration, which are previously provided in Eq. (1c) and Eq. (1d), respectively. Finally, it is possible to represent the variations in fluid properties across the interface, including density and viscosity, as a continuous function of the phase field variable. This expression is then utilized in Eq. (14). In this work (see e.g., [96]), we employ the following linear interpolation to account for the interfacial variations of fluid properties:

$$\rho = \rho _ { B } + \frac { \phi - \phi _ { A } } { \phi _ { A } - \phi _ { B } } \left ( \rho _ { A } - \rho _ { B } \right ), \quad \mu = \mu _ { B } + \frac { \phi - \phi _ { A } } { \phi _ { A } - \phi _ { B } } \left ( \mu _ { A } - \mu _ { B } \right ),$$

where ρ A , ρ B and µ A , µ B are the densities and the dynamic viscosities in the fluid phases, respectively and denoted by ϕ A and ϕ B . An equation similar to Eq. (17) will also be utilized for distributing the interfacial jump in the thermal conductivity in solving the energy equation. In this study, we use ϕ A = 0 and ϕ B = 1.

## 4 Central moment lattice Boltzmann schemes for simulation of surfactant-laden self-rewetting fluid motions

This section will introduce a numerical LB approach based on more robust collision models involving central moments [61,77,81,97,98] for solving the equations of the phase-field model for tracking the interface (Eq. (12)) and the binary fluid motions (Eqs. (13)-(16)) given in the previous section, along with the transport of energy and the surfactant concentration presented in Eqs. (1c) and Eqs. (4) and (5), respectively, earlier. In order to solve these four equations, four different distribution functions need be evolved on the standard twodimensional square lattice (D2Q9). This is done by performing a collision step that is based on the relaxation of the different central moments of such distribution functions to their respective equilibria. The distribution functions then undergo lock-step advection to their adjacent nodes along the characteristic directions in the streaming step . The moments of

the corresponding distribution functions are then used to derive the macroscopic variables, which include the order parameter, the fluid pressure and velocity, as well as the temperature and the surfactant concentration fields. Note that this requires the use of appropriate mappings that transform between these quantities pre- and post-collision steps since the streaming step is carried out by means of the distribution functions, while the collision step is carried out using central moments. When compared to the other collision models in the LB framework, the central moment LB approaches are demonstrated to be more robust (e.g., enhanced numerical stability) (see [61,77,99,100] for recent examples). Although an orthogonal moment basis was used to develop the original central moment LB scheme for two-fluid interfacial flows [77], we will adopt an improved structure employing the nonorthogonal moment basis in the following.

Brief descriptions of the central moment LB schemes and attendant details fof the capturing of interfaces based on the conservative Allen-Cahn equation, two-fluid motion, energy transport, and surfactant concentration transport, respectively, are presented Appendices B, C, D, and E, respectively. While these methods are applicable for general class variable surface-tension driven flows, in this work, they will be mainly applied to study the effect of various characteristic parameters on the flow patterns and the intensity of thermocapillary convention in migrating surfactant-laden bubble in self-rewetting fluids (SRFs). Before proceeding with this, we will first validate our numerical approach given above for some standard thermocapillary benchmark problems for normal fluids for which analytical solutions are available in the literature in the next section.

## 5 Model validation

## 5.1 Laplace-Young law of a static drop

The static drop test case is first used to validate the capability of this model in the absence of body force effects or surfactants. In this test, the analytical solution from the LaplaceYoung's relation for a 2D drop at rest relates the pressure difference inside and outside of the drop ( △ P ) of radius R to the surface tension σ and its radius of curvature 1 /R via △ P = σ/R , which will be used as a comparison. To do the simulation, we consider a drop of density ρ A is placed at the center of a square domain that is divided into 200 × 200 lattice nodes and filled with an ambient fluid of density ρ B . Periodic boundary conditions are specified for all walls. For four different surface tensions of σ = 1 × 10 -3 , 2 × 10 -3 , 3 × 10 -3 , and 4 × 10 -3 , we ran the simulations at a density ratio of 1000. The estimated pressure differences between the drop and ambient fluid are displayed in Fig. 3 in comparison to the Laplace-Young relation equation. It is evident that there is a predicted linear relationship between △ P and 1 /R ; a good quantitative agreement between the analytical solution and the computed results is seen.

Fig. 3. Comparison of the computed pressure differences (symbols) against the analytical predictions using the Laplace-Young relation for various values of the drop curvature 1/R with four different values of surface tension.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00146v2-figure-eeb6a325-a701-4a1c-8b58-f903d05ba828.jpg)

## 5.2 Equilibrium distribution of surfactant concentration on a flat interface

To further validate the model quantitatively in the presence of surfactants, we will next compare the numerical results with an analytical solution representing the equilibrium surfactant concentration profile at a planar interface. A schematic figure of this problem is shown in Fig. 4.

Fig. 4. Schematic of the initial and final distributions of the surfactant concentration field around a flat interface illustrating the tendency of interfacial surfactant adsorption.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00146v2-figure-f3b1326b-df2c-4979-ab4f-acfb84ef4025.jpg)

This analytical solution can be derived when the chemical potential µ ψ (see Eq. (3)) is uniformly distributed throughout the domain. By setting µ ψ = µ ψ,b , where µ ψ,b corresponds to the value in the bulk region that can be determined by µ ψ,b = λ ln ψ b + w ϕ ( A + ϕ B ) 2 / 8, we obtain the analytical solution for the equilibrium surfactant concentration profile as [90]

$$\psi ^ { e q } ( \zeta ) = \frac { \psi _ { b } } { \psi _ { b } + \psi _ { c } ( \zeta ) }.$$

Here, ψ b is the fixed bulk surfactant concentration loading, and ψ c ( ζ ) in an auxiliary function determined by

$$\psi _ { c } ( \zeta ) = \exp \left \{ - \frac { s } { 2 \lambda } \left [ \frac { 2 } { W } \phi _ { m } \, \text{sech} ^ { 2 } \left ( \frac { 2 \zeta } { W } \right ) \, \right ] ^ { 2 } + \frac { w } { 2 \lambda } \left [ \phi ( \zeta ) - \phi _ { m } \right ] ^ { 2 } - \frac { w } { 8 } ( \phi _ { A } + \phi _ { B } ) ^ { 2 } \right \}.$$

Here, ϕ ζ ( ) = ϕ m + ϕ o tanh(2 ζ/W ) is the corresponding equilibrium hyperbolic tangent profile of the phase variable, ϕ m = ( ϕ A + ϕ B ) / 2, ϕ o = ( ϕ A -ϕ B ) / 2 and ζ is a coordinate along the normal direction that originates at the interface in the above equation. For this test, we select a computational domain with a grid resolution of 5 × 101 and an interface width of W = 5. The values of the parameters s , w , and λ may be computed using the non-dimensional numbers Ex = 0 025 and Pi = 50. . Other parameters are set as follows: σ 0 = 0 01 . , M ϕ = M ψ = 0 02. . The profiles of equilibrium surfactant concentration are obtained through numerical simulations for three different bulk surfactant concentration values: ψ b = 0 05 . , ψ b = 0 15 and . , ψ b = 0 225. . These profiles are displayed alongside the analytical solution for comparison in Fig. 5.

Fig. 5. Profiles of the surfactant concentration ψ for a planar interface for various bulk surfactant concentrations ψ b = 0 05 0 15 and 0 225. . , . , . The numerical solution and the analytical prediction are represented by symbols and solid lines, respectively.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00146v2-figure-d0fa28df-9895-4086-bd8a-6607185bc192.jpg)

It is evident that the interface is where the surfactant concentration is highest. This represents the preferred surfactant adsorption at the interface, which is resisted by diffusion. Furthermore, the estimated equilibrium profiles and the analytical solution match up nicely, further confirming our choice to solve the surfactant concentration problem using the central moment LB method accurately.

## 5.3 Buoyancy-driven rising bubble with large density ratio

Finally, we will simulate a rising air bubble in water with a high density ratio and a complex interface change that includes the surface tension force, buoyancy, and viscous effect in order to further test our computational approach; this case study serves as a precursor to studying the physics of the surfactant-laden bubble migration in SRFs that will be discussed in the next section. We study, under various parametric conditions, a bubble with a diameter of D and a density of ρ B rising due to buoyancy forces in an ambient fluid with a density of ρ A , where the density ratio ρ A /ρ B equals 1000 and the dynamic viscosity ratio is given as µ A /µ B = 100. Our objective is to evaluate the capability of the central moment LB approach to reproduce the temporal history of the bubble route quantitatively and to represent the various shape changes the bubble experiences due to the balance between competing forces.

The non-dimensional parameters specified in this problem are the Reynolds number Re and the Eotvos number Eo, which can be defined as follows:

$$\text{Re} = \frac { \rho _ { A } U _ { 0 } D } { \mu _ { A } }, \quad \text{Eo} = \frac { \rho _ { A } U _ { 0 } ^ { 2 } D } { \sigma },$$

where g is the gravitational acceleration and σ denotes the surface tension coefficient. The

length scale L = D , the velocity scale U 0 = √ gD , which denotes a characteristic gravitational velocity, and the time scale T = L/U 0 are taken into consideration in this problem. In this benchmark problem, a body force of F text = -( ρ -ρ A ) g j is applied to set the bubble in motion. In a rectangular domain of 161 × 481 lattice nodes, a gas bubble with a diameter of D = 80 is centered at a location (80 120). , Periodic boundary conditions are used on the left and right sides and no-slip boundary conditions are imposed on the upper and lower walls. For reporting time histories, we use the dimensionless time which is denoted by t ∗ = t/T . Depending on the magnitudes of these dimensionless groups, the bubble experiences intricate interfacial shape changes, leading to a variety of forms, including spherical-cap, dimpled ellipsoidal-cap, and skirted (see [101] for more details). In this work, we perform buoyancy-driven bubble rise simulations at four Eotvos numbers ( i.e., Eo = 0 1 10 50 and 125) to illustrate the interface shapes of rising bubbles in different . , , , flow regimes. We set ρ A /ρ B = 1000 and µ A /µ B = 100 at a fixed Reynolds number Re = 35.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00146v2-figure-2296ef2a-60af-4e8e-af7b-f3ecf4e6bd78.jpg)

(d) Eo = 125

Fig. 6. Evolution of the interface for the bubble rising at Re = 35 and ( a ) Eo = 0 1, ( ) . b Eo = 10, ( ) Eo = 50, and ( c d ) Eo = 125.

Figure 6 shows the computed evolution of the interface of the rising bubble at the above

four values of Eo. Due to the strong surface tension force, the form of the gas bubble remains unchanged for the rising bubble with Eo = 0 1. . When the surface tension force plays a relatively large role compared to the other forces, i.e., when the Eotvos number is low (Eo = 10), the bubble goes through smaller deformation that is beginning at its rear end, which then leads in a flattening of that side as the bubble rises. When Eo = 50, the driving buoyancy force exceeds the surface tension under the dominant viscous force, causing a substantially larger deformation via stretching that results in the formation of tails that extend with time. The gas bubble shape significantly changes at larger Eotvos numbers (Eo = 125) due to decreased surface tension forces, and two extended tails can be seen behind the main bubble. The findings obtained through other approaches [102-104] are quite comparable to these estimated shape changes with varying Eo at different times.

In order to evaluate the current approach with more confidence, a quantitative comparison is required. The mass center of the rising bubble as its shape changes with time is shown in Fig. 7 for the case of Re = 35 and Eo = 125. The computational results are compared with numerical results available in the literature [103] used as a reference. The fact that the current results are in good quantitative agreement with the range of the reference data indicates that the current approach is accurate in solving this problem with a high density ratio.

Fig. 7. Comparison of the non-dimensional center of mass for the rising bubble with the reference numerical results given in Wang et al. [103] for Re = 35 and Eo = 125.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00146v2-figure-7448a492-8032-49e4-a590-226a2d7c911c.jpg)

Finally, we note here that the present LB numerical approach was further tested and validated against a new analytical solution for thermocapillary convection due to Marangoni stresses in two superimposed SRF liquid layers subject to a nonuniform heating in our recent study [61].

## 6 Results and Discussion

̸

̸

The goal of this work is to model and study the physics of the migration of surfactantladen thermocapillary bubbles in a self-rewetting fluid (SRF) by investigating the effects of different characteristic parameters associated with this problem using the LB schemes, which were validated in the previous section. In this regard, we will study the cases where the dimensionless surface tension quadratic sensitivity coefficient is non-zero, i.e., σ TT = 0 or M 2 = 0 to illustrate the role of self-rewetting fluids (SRFs), and also to compare and contrast the results with the normal fluids (NFs), where only the linear coefficient of the surface tension (i.e., σ T or M 1 ) exists and M 2 is zero. The surfactant-laden thermocapillary flow strengths are controlled by the choice of the dimensionless linear and quadratic sensitivity coefficients of the surface tension variation with temperature, i.e., M 1 and M 2 , respectively, and also by the choice of the Gibbs elasticity parameter β . In this study, unlike otherwise stated, the following parameters are kept constant: T H = 1, T C = 0, T ref = 0 5, . D = 0 1 . H , ˜ = 1000, ρ ˜ = 100, µ ˜ = 1, Re = 5, and Ma = 0 25. α . Thus, in particular, in all the cases considered in what follows, the left end of the domain is maintained hotter at T H than the right end which is at T C (see Fig. 2 for a schematic of the problem setup). The reference surface tension for the clean interface case used in this study is assumed to be σ 0 = 5 × 10 -3 . We used W = 5 and M ϕ = 0 1 for the model parameters associated with the conservative . ACE for interface tracking. In addition, in the surfactant model, the non-dimensional numbers Ex and Pi have been selected to be as Ex = 0 16 and Pi = 8 33, similar to the . . values used by previous studies in Refs. [88,89] with M ψ = 0 01 and were kept fixed in all . the simulations. During the simulations, we determine the migration velocity of the bubble ( U ) by calculating the average velocity inside the bubble. The reference time is specified by the scale T o = R/U o where R is the radius of the bubble ( R = D/ 2), and therefore, the dimensionless time used is t ∗ = t/T o . Here, the characteristic thermocapillary velocity U o can be obtained from balancing the surface tension force for the clean interface with the viscous force, which reads

$$U _ { o } = \left ( \frac { \sigma _ { 0 } } { \mu _ { o } } \right ) \left ( \frac { R } { L } \right ) ( M _ { 1 } + M _ { 2 } ) \,.$$

## 6.1 Effect of Gibbs elasticity parameter β for normal fluids at fixed M 1

To begin with, in order to establish a reference point for comparison, we will initially present results for the velocity and position of a bubble in a specific scenario involving NFs, as depicted in Fig. 8. This particular case considers M 1 = -2 5 and . M 2 = 0 0. . Furthermore, the selection of other fluid characteristics is such that the Gibbs elasticity parameter is varied in the range β = 0 0 2 0 3 and 0 4, , . , . , . while the concentration of surfactant on the left and right sides are ψ L = 0 1 and . ψ R = 0 4, . respectively. As we can see in Fig. 8 for the clean interface case, i.e., without surfactants or β = 0, the bubble moves, after an initial transient phase, at a constant velocity completely to the end of the domain in one direction and adding surfactants via varying β such that β = 0 2 0 3 and 0 4 causes the . , . , . rate of motion of the bubble to decrease. Figure 10 shows the contours of the bubble at selected times for the case when β = 0 and β = 0 3, where, for brevity in presentation, the . snapshots are taken by only focusing on the relevant section of the bubble path by slicing

̸

̸

horizontally from the full domain presented in Fig. 9. The motion of the bubble from right to left in the surfactant-laden case ( β = 0) leads to a nonuniform distribution of the surfactant concentration field (see the contours shown in Figure 10), which in turn generates a tangential surface tension gradient or the Marangoni stress that opposes the motion of the bubble. Thus, the greater the sensitivity of surface tension on the surfactant field or β , the longer it takes for the moving bubble to cover a given distance (see Fig. 8b), since it generates progressively greater tangential Marangoni stresses due to surfactant adsorption on its interface that retard the rate of bubble motion which manifests as a decrease in | ¯ U | as β is increased (see Fig. 8a). In effect, the bubble migrating in a normal fluid does not settle to an equilibrium position and monotonically and continuously moves from right to left and it attains a finite and non-zero terminal velocity ¯ = 0 and having a negative sign ( ¯ U U &lt; 0) indicating the motion is directed opposite to the positive direction of the coordinate system.

̸

In contrast to Fig. 8 where the surfactant concentration at the right and the left ends are ψ R = 0 4 and . ψ L = 0 1, we now increase the concentration at the right end to . ψ R = 0 6 . while keeping ψ L = 0 1, as before, . in Fig. 11. Effectively, this increases the streamwise gradient magnitude of the surfactant concentration | dψ/dx | . As a result, comparing Figs. 8 and 11 for the same β , it can be seen that the magnitude of the bubble velocity decreases relative to the clean interface case ( β = 0) further when ψ R = 0 6 when compared to . ψ R = 0 4. . By contrast, when the right end of the domain is loaded to ψ R = 0 2, such that . | dψ/dx | is lower compared to the previous two cases, the bubble travels faster for β = 0 (see Fig. 12). We emphasize that these observations are applicable for surfactant-laden thermocapillary bubble migration in normal fluid only.

Fig. 8. Bubble velocity( a ) and its location ( b ) at different values of the Gibbs elasticity parameter β for surfactant-laden bubble in a normal fluid (NF). In this figure, Re = 5, Ma = 0 25, . M 1 = -2 5, . M 2 = 0, D/H = 0 1, . ψ L = 0 1, and . ψ R = 0 4. .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00146v2-figure-1ff6397d-7d1d-48fb-8240-2442b2c09a8b.jpg)

Fig. 9. This figure is used to indicate the horizontal slice of domain, using the black dashed lines, as seen in the time sequence snapshots in Fig. 10 below. The color represents the magnitude of the normalized surfactant concentration field for the baseline case of parameters presented in Fig. 8. Similar cases are shown without color that do not have surfactant, but have the same slice of data shown.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00146v2-figure-4c7c0fec-6849-43e8-8912-b674e097b8ca.jpg)

Fig. 10. Time sequence snapshots of bubble migration in a normal fluid with and without surfactant. The time stamps from top to bottom are t ∗ = [0 50 100 150 200 250] for , , , , , β = 0 (left) and β = 0 3 (right), . and the slice of the domain for each snapshot is indicated in Fig. 9. Some relevant parameters for this problem are as follows: Re = 5, Ma = 0 25, . M 1 = -2 5, . M 2 = 0, D/H = 0 1, . ψ L = 0 1, and . ψ R = 0 4. .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00146v2-figure-8871734b-ed43-4b8e-9ac0-6a021ea267e3.jpg)

Fig. 11. Bubble velocity( a ) and its location ( b ) at different values of the Gibbs elasticity parameter β for surfactant-laden bubble in a normal fluid (NF). In this figure, Re = 5, Ma = 0 25, . M 1 = -2 5, . M 2 = 0, D/H = 0 1, . ψ L = 0 1, and . ψ R = 0 6. . Effective magnitude of the streamwise gradient of the surfactant concentration field | dψ/dx | is higher here when compared to Fig. 8.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00146v2-figure-4e96b411-5d16-4031-9a15-ee063f8dd21a.jpg)

Fig. 12. Bubble velocity( a ) and its location ( b ) at different values of the Gibbs elasticity parameter β for surfactant-laden bubble in a normal fluid (NF). In this figure, Re = 5, Ma = 0 25, . M 1 = -2 5, . M 2 = 0, D/H = 0 1, . ψ L = 0 1, and . ψ R = 0 2. . Effective magnitude of the streamwise gradient of the surfactant concentration field | dψ/dx | is lower here when compared to Fig. 8.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00146v2-figure-12406fe1-7abb-45e0-87ae-bcbd1fba57f8.jpg)

## 6.2 Effect of Gibbs elasticity parameter β for self-rewetting fluids at larger M 2

̸

On the other hand, when we turn off the linear coefficient of surface tension (i.e., σ T = 0 or M 1 = 0) and only retain a non-zero quadratic coefficient, i.e., σ TT = 0, we conduct

simulations to explore the migration of surfactant-laden thermocapillary bubbles in selfrewetting fluids (SRFs), while keeping the same values for the property ratios and other parameters as before. Specifically, in dimensionless form, we set M 1 = 0 and M 2 = 5 as the selected values for the SRFs' baseline case to show the effect of the Gibbs elasticity parameter, β . We impose a surfactant concentration gradient by setting ψ L = 0 4 and . ψ R = 0 1. .

̸

The computed results, in terms of the bubble's velocity and position for different β , are depicted in Fig. 13. First, for β = 0 (i.e., with no surfactant sensitivity), the bubble initially speeds up and then slows down and continues to slow down until it reaches zero velocity ( ¯ U = 0) attaining an equilibrium position at the reference temperature location. This behavior is dramatically different from the normal fluid case, where, as seen earlier, there was no equilibrium position of the bubble which executes a continuous motion from the right end to the left end of the domain, i.e., from the cold to the hot sides; by contrast, the bubble in the self-rewetting fluid moves from left to right or from the hot to the cold sides eventually attaining a stationary position with ¯ U = 0. Moreover, observe that the presence of surfactant (i.e., β = 0) causes the bubble velocity to decrease for this case, since the initial bubble location ( x/L ) initial = 0 25 is . in the higher concentration part of the surfactant gradient. During the transient phases, the greater the Gibbs elasticity parameter β , the larger is the reduction in the bubble velocity, but in all cases at long times they are come to rest, i.e., attain zero velocity (see Fig. 13a). Interestingly, the combined effects of the Marangoni stresses arising from the tangential gradients in the surface tension due to non uniform temperature field and the surfactant concentration field (see Eq. (6)), causes the bubble to attain an equilibrium position at long times that is sensitive to the presence of surfactant field and the degree of its propensity to get adsorbed in the interface. In particular, as seen in Fig. 13b, only when there is no surfactant effect (i.e., the clean interface case) with β = 0, the equilibrium bubble position x eq is the same as the reference temperature location at x = 0 5 . L . In all other cases, as the surfactant's sensitivity or β increases the equilibrium position of the bubble x eq is located upstream, i.e., x eq /L &lt; 0 5; in other words, the greater the Gibbs elasticity parameter . β , the further the bubble's equilibrium location is pushed upstream or towards regions with higher surfactant concentration. From a physical point of view, when the thermally-induced Marangoni stresses for SRFs (which depends on M 2 ) is equal to the counteracting surfactant-induced Marangoni stresses (which depends on β ), the bubble stops moving and thereby attaining an equilibrium position that is dependent on β for a fixed M 2 .

Figure 13 shows these effects with the baseline case with ψ L = 0 4 and . ψ R = 0 1 for . setting up the surfactant concentration gradient along the streamwise direction. We now explore the role of changing this streamwise gradient | dψ/dx | by decreasing ψ L = 0 2 while . keeping ψ R = 0 1 (see Fig. 14) and the opposite case with higher . ψ L = 0 6 for the same . ψ R (see Fig. 15). Clearly, when | dψ/dx | is increased as in Fig. 15, there is greater reduction in the bubble velocity for a specified β when compared to the baseline case during the initial transients while the bubble halting completely at long times in each case; more importantly, comparing Figs. 13b and 15b, it can be seen that for a specified β , the bubble equilibrates at a location further upstream for the later case with greater | dψ/dx | . Thus, ( x eq /L ) G s,higher &lt; ( x eq /L ) G s,lower , where G s = | dψ /dx | is the magnitude of the imposed streamwise gradient of the surfactant concentration field. For example, when β = 0 4, the .

bubble's equilibrium position ( x eq /L ) is shifted from about 0 45 to about 0 4 when . . | dψ/dx | is increased via changing ψ L from 0 4 to 0 6. . . The contour plots for the bubble at selected times for the clean interface case ( β = 0) and the surfactant-laden case ( β = 0 3) presented . in Fig. 17 corroborates these findings.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00146v2-figure-ae1d5e6c-b68b-4dbe-82d8-3f6404437489.jpg)

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00146v2-figure-58162f3b-e575-4421-81ee-5a53a0da4695.jpg)

(b)

Bubble Centroid Location

Fig. 13. Bubble velocity( a ) and its location ( b ) at different values of the Gibbs elasticity parameter β for surfactant-laden bubble in a self-rewetting fluid (SRF). In this figure, Re = 5, Ma = 0 25, . M 1 = 0, M 2 = 5, D/H = 0 1, . ψ L = 0 4, and . ψ R = 0 1. .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00146v2-figure-3e08db38-a0a0-4aa3-aae7-94cde1b5b85a.jpg)

Fig. 14. Bubble velocity( a ) and its location ( b ) at different values of the Gibbs elasticity parameter β for surfactant-laden bubble in a self-rewetting fluid (SRF). In this figure, Re = 5, Ma = 0 25, . M 1 = 0, M 2 = 5, D/H = 0 1, . ψ L = 0 2, and . ψ R = 0 1. . Effective magnitude of the streamwise gradient of the surfactant concentration field | dψ/dx | is lower here when compared to Fig. 13.

Fig. 15. Bubble velocity( a ) and its location ( b ) at different values of the Gibbs elasticity parameter β for surfactant-laden bubble in a self-rewetting fluid (SRF). In this figure, Re = 5, Ma = 0 25, . M 1 = 0, M 2 = 5, D/H = 0 1, . ψ L = 0 6, and . ψ R = 0 1. . Effective magnitude of the streamwise gradient of the surfactant concentration field | dψ/dx | is higher here when compared to Fig. 13.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00146v2-figure-e158a9a1-5f80-46a2-968f-c6344bcb000a.jpg)

Fig. 16. This figure is used to indicate the horizontal slice of domain, using the black dashed lines, as seen in the time sequence snapshots in Figs. 17, 20, and 24 below. The color represents the magnitude of the normalized surfactant concentration field for the baseline case of parameters presented in Fig. 13. Similar cases are shown without color that do not have surfactant, but have the same slice of data shown.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00146v2-figure-acc17a39-658d-42fc-b8b9-3587c20b63aa.jpg)

Fig. 17. Time sequence snapshots of bubble migration in a self-rewetting fluid with and without surfactant. The time stamps from top to bottom are t ∗ = [0 50 100 200 500 1200] , , , , , for β = 0 (left) and β = 0 3 (right), . and the slice of the domain for each snapshot is indicated in Fig. 16. Some relevant parameters for this problem are as follows: Re = 5, Ma = 0 25, . M 1 = 0, M 2 = 5, D/H = 0 1, . ψ L = 0 4, and . ψ R = 0 1. .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00146v2-figure-0ae0531a-0f29-4a3e-986c-8e9c3584d21e.jpg)

## 6.3 Effect of Gibbs elasticity parameter β for self-rewetting fluids at smaller M 2

Moving forward, we now investigate the implications of reducing the dimensionless surface tension quadratic sensitivity coefficient from M 2 = 5 0 for the baseline case discussed above . to M 2 = 1 25 and its impact on the average velocity of the bubble and its centroid location . for different Gibbs elasticity parameter β . The results of this simulation case are presented in Fig. 18. As in the higher M 2 case, here the bubble in the SRF moves from the left to the right sides while experiencing a deceleration during the initial transients as β is increased since it is launched from an initial position ( x/L ) initial = 0 25 which is in the . higher concentration region of the surfactant field. However, the results reveal that as M 2 decreases, the magnitude of the thermocapillary effects is reduced which causes the bubble to reduce its speed and, as a result, the impact of surfactants becomes even more prominent compared to the baseline case considered in the previous section. Specifically, when the bubble moves at a slower rate due to lower M 2 , the thermally-induced contribution to the Marangoni stress is smaller and so the relative contribution of this tangential stress from surfactant's interfacial adsorption for different and higher Gibbs elasticity parameter leads to progressively larger effects.

It would be interesting to see the bubble dynamics in SRF when it is launched from a different initial position compared to the above, i.e., from ( x/L ) initial = 0 75 just like in the . NF case discussed earlier. Figure 19 shows the results of the bubble velocity and its centroid location as a function of time. Since the minimum temperature occurs at the center of the domain corresponding to the low surface tension region, the bubble seeks to move towards

it. Thus, in this case, it moves from right to left. Moreover, since here the initial location of the bubble is in a lower surfactant concentration region, during the initial transients it accelerates unlike that in the case discussed above (see Fig. 18). However, in both cases, unlike in the NF case, the bubble equilibrates, and its final position is determined by the balance of the Marangoni stresses induced by thermocapillary and surfactant effects. The previous work on bubble migration in SRF in the absence of surfactants [56-59] showed it always attains a final equilibrium position of x eq /L = 0 5 . , where the temperature is minimum . On the other hand, this study reveals a potentially new approach to modulate the bubble equilibrium position in SRF by adding surfactants in such a way that the Marangoni effects due to thermocapillary effects are counteracted by the solutal effects. In particular, we obtain the general result that regardless of how the bubble is initially launched, after initial transients, the bubble in SRF attains equilibrium at a location different from x eq /L = 0 5 . (corresponding to the clean interface case) that can be controlled by the Gibbs elasticity parameter β . This can be seen by comparing Fig. 18(b) and 19(b).

These observations for the bubble dynamics in SRF in the presence of surfactants at M 2 = 1 25 for . different β at ψ L = 0 4 and . ψ R = 0 1 can be clearly . seen in Fig. 18 by contrasting it with Fig. 13. In particular, the equilibrium position of the bubble with M 2 = 1 25 is shifted further towards the hotter, upstream side under the weakened . thermally-induced Marangoni stress to a position where it balances that from the surfactantgenerated Marangoni stress. For example, the bubble equilibrium position is further shifted to x eq /L = 0 24 for . β = 0 4 at . M 2 = 1 25 when compared to . x eq /L = 0 45 for the same . β at M 2 = 5 0 in Fig. 13, which corresponds to a significant relative change in the equi-. librium position by about 87 5%. . Thus, the inclusion of surfactants provides an ability to manipulate the equilibrium position of the bubble, contingent upon the chosen value of β related to sensitivity of surface tension on surfactant concentration, which is a material property, and the dimensionless quadratic thermocapillary coefficient M 2 . Figure 20 shows the contours of a bubble during its migration in SRF from left to its equilibrium location for different dimensionless times for β = 0 (left) and β = 0 3 (right) and for . M 2 = 1 25. We . observe a behavior similar to Fig. 17 but with a considerably stronger effect on modulating the bubble equilibrium location where the bubble always migrates from the hot side to the cold side, which is opposite in direction with using NF, towards the reference temperature location occurring at the center of the domain based on the setup of the problem.

Fig. 18. Bubble velocity ( a ) and its location ( b ) at different values of the Gibbs elasticity parameter β for surfactant-laden bubble in a self-rewetting fluid (SRF). In this figure, Re = 5, Ma = 0 25, . M 1 = 0, M 2 = 1 25, . D/H = 0 1, . ψ L = 0 4, and . ψ R = 0 1. .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00146v2-figure-cf9c2c08-5ddd-4d05-b970-03667fa02bb4.jpg)

Fig. 19. Bubble velocity ( a ) and its location ( b ) at different values of the Gibbs elasticity parameter β for surfactant-laden bubble in a self-rewetting fluid (SRF). In this figure, Re = 5, Ma = 0 25, . M 1 = 0, M 2 = 1 25, . D/H = 0 1, . ψ L = 0 4, and . ψ R = 0 1. . Here, the initial location of the bubble is at ( x/L ) initial = 0 75 compared that in Fig. 18 where the . bubble is launched from the initial location ( x/L ) initial = 0 25. .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00146v2-figure-ee2ce230-d9cf-48a5-a371-fbb2064df11a.jpg)

Fig. 20. Time sequence snapshots of bubble migration in a self-rewetting fluid with and without surfactant. The time stamps from top to bottom are t ∗ = [0 50 100 200 500 1200] , , , , , for β = 0 (left) and β = 0 3 (right), . and the slice of the domain for each snapshot is indicated in Fig. 16. Some relevant parameters for this problem are as follows: Re = 5, Ma = 0 25, . M 1 = 0, M 2 = 5, D/H = 0 2, . ψ L = 0 4, and . ψ R = 0 1. .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00146v2-figure-ecf27c1d-affc-428b-902a-36ca82e03482.jpg)

In order to encompass a wider range of situations that accommodates greater role of the surfactant effects when compared to thermocapillarity in SRFs, we now reduce the quadratic sensitivity coefficient of surface tension M 2 even further by a factor of 2.5 when compared to the baseline case reported in Fig. 18, which uses M 2 = 1 25. . Thus, Fig. 21 shows the results of the bubble velocity and its centroid location in SRF at different values of the Gibbs elasticity parameter with M 2 = 0 5. . It can be observed that the bubble again attains equilibrium at long times that is sensitive to the Gibbs elasticity parameter similar to the baseline case shown in Fig. 18 but with an enhanced effect due to the presence of surfactants that results in greater deviations from the equilibrium location at x eq /L = 0 5 . with increase in β when compared to the case with smaller M 2 . As such, the general trends and findings are similar to the other cases considered earlier.

Fig. 21. Bubble velocity ( a ) and its location ( b ) at different values of the Gibbs elasticity parameter β for surfactant-laden bubble in a self-rewetting fluid (SRF). In this figure, Re = 5, Ma = 0 25, . M 1 = 0, M 2 = 0 5, . D/H = 0 1, . ψ L = 0 4, and . ψ R = 0 1. . Here, the choice of M 2 = 0 5 represents a reduction by a factor of 2 5 when compared to the baseline . . case shown in Fig. 18 with M 2 = 1 25. .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00146v2-figure-9f79a280-4cf7-4229-834d-42c35555b6de.jpg)

In the above cases, the gradient in surfactant concentration to achieve counteracting Marangoni stresses to thermocapillary effects is introduced by specifying the surfactant concentrations on the left and right sides of the domain to be at ψ L and ψ R , respectively, following our previous work [105], as representative of the presence of their sources at the respective ends of the domain. Nevertheless, we now investigate an alternative approach by specifying ψ L and ψ R as the initial condition at those locations and then allowing them to vary via the Neumann condition. The results with this approach for the baseline case with M 2 = 1 25 for the bubble velocity and its centroid location are shown in Fig. 22. . Compared to the baseline case shown in Fig. 18, the results are generally qualitatively similar with somewhat reduced effect of surfactants at higher β in modulating the equilibrium location of the bubble. Overall, the main conclusions made earlier regarding the influence of surfactants on the bubble dynamics in SRFs are still valid.

Fig. 22. Bubble velocity ( a ) and its location ( b ) at different values of the Gibbs elasticity parameter β for surfactant-laden bubble in a self-rewetting fluid (SRF). In this figure, Re = 5, Ma = 0 25, . M 1 = 0, M 2 = 1 25, . D/H = 0 1, . ψ L = 0 4, and . ψ R = 0 1. . Here, ψ L and ψ R are set only initially and then allowed to vary through the Neumann condition compared that in Fig. 18 where ψ L and ψ R are set as boundary conditions at the respective boundaries.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00146v2-figure-799ef7d9-3ec1-4ba9-a1fe-7c9de9547d4f.jpg)

## 6.4 Effect of bubble size or diameter D

Next, we explore the effect of bubble size on its dynamics when moving in a SRF for different surfactant sensitivities. In all the previous cases, we considered the bubble size D to be such that D/H = 0 1. In this section, we double the bubble size by setting . D/H = 0 2 and study . its influence on its velocity and position. The fluid properties for this study are similar to those mentioned previously and are also shown in the figures captions. Smaller bubbles have a higher curvature, which can affect the distribution of surfactants on their surface. This can lead to more pronounced surface tension gradients, causing stronger Marangoni flows. On the other hand, larger bubbles have a smaller curvature and could result in a more uniform distribution of surfactants around the interface. The surface tension gradients are less pronounced compared to smaller bubbles, thereby resulting in weaker Marangoni-driven flows. The size of the bubble thus influences its migration velocity under Marangoni forces.

Due to stronger Marangoni flows, smaller bubbles (see Fig. 13) migrate faster and slow down more quickly compared to larger bubbles (see Fig. 23), especially if the effect is primarily driven by surface tension gradients due to solely thermocapillary effects in the absence of surfactants in SRFs (i.e., β = 0). In other words, increasing the bubble size causes a significant reduction in its velocity during the initial transients while the final equilibrium position is not significantly influenced by this change. Moreover, the bubble contours during its migration in a SRF from left to its equilibrium location for different dimensionless times for the larger bubble size presented in Fig. 24, shows that the bubble reaches its equilibrium location relatively faster than the case with the smaller bubble (see Fig. 17).

Fig. 23. Bubble velocity( a ) and its location ( b ) at different values of the Gibbs elasticity parameter β for surfactant-laden bubble in a self-rewetting fluid (SRF). In this figure, Re = 5, Ma = 0 25, . M 1 = 0, M 2 = 5, D/H = 0 2, . ψ L = 0 4, and . ψ R = 0 1. .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00146v2-figure-fe3e02d7-3311-4cbf-8443-199ca0769666.jpg)

Fig. 24. Time sequence snapshots of bubble migration in a self-rewetting fluid with and without surfactant. The time stamps from top to bottom are t ∗ = [0 50 100 200 500 1000] , , , , , for β = 0 (left) and β = 0 3 (right), . and the slice of the domain for each snapshot is indicated in Fig. 16. Some relevant parameters for this problem are as follows: mboxRe = 5, Ma = 0 25, . M 1 = 0, M 2 = 5, D/H = 0 2, . ψ L = 0 4, and . ψ R = 0 1. .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00146v2-figure-13eea859-ef1f-44b6-9e0d-b63eafb55165.jpg)

## 6.5 Effect of Reynolds number Re

We now study the effect of the Reynolds number Re on the bubble velocity and its position during its migration by varying Re in the range 0 1 1 5 10. . , , , The rest fluid properties for

this study are similar to those specified in the previous case and are also shown in the captions of the associated figures. The baseline case shown in Fig. 13 already presented the case Re = 5 and Figs. 25, 26, and 27 show the results for Re = 0 1 . , 1 0 . , and 10, respectively. As Fig. 25 illustrates, bubbles with lower Re reach a higher peak velocity during the initial transients than those with higher Re, as Fig. 27 shows. The reason why bubbles in low Reynolds number environments generally move more quickly than those at higher Reynolds numbers in such thermocapillary-driven motions of the fluids and interfaces is that the greater viscous effects in low Reynolds number situations tend to transmit effects of the surface tension gradients or Marangoni stresses from the interfaces to the bulk fluid more effectively. The presence of surfactants reduces the concentration on the trailer edge of the bubble and, in the case of a lower Re, alter its velocity during the transients the bubble velocity and shift the equilibrium location slightly upstream toward a higher surfactant gradient location. Generally, for the cases studied, it is found that the overall effect of Re is not very significant compared to the influence of other characteristic parameters such as M 2 and β .

Fig. 25. Bubble velocity( a ) and its location ( b ) at different values of the Gibbs elasticity parameter β for surfactant-laden bubble in a self-rewetting fluid (SRF). In this figure, Re = 0 1, Ma = 0 25, . . M 1 = 0, M 2 = 5, D/H = 0 1, . ψ L = 0 4, and . ψ R = 0 1. .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00146v2-figure-9efd7c60-66c0-4e4a-958f-d9e49ce73cf1.jpg)

Fig. 26. Bubble velocity( a ) and its location ( b ) at different values of the Gibbs elasticity parameter β β for surfactant-laden bubble in a self-rewetting fluid (SRF). In this figure, Re = 1, Ma = 0 25, . M 1 = 0, M 2 = 5, D/H = 0 1, . ψ L = 0 4, and . ψ R = 0 1. .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00146v2-figure-034fc073-8b0c-40c9-af7d-4c79aa8088ca.jpg)

Fig. 27. Bubble velocity( a ) and its location ( b ) at different values of the Gibbs elasticity parameter β for surfactant-laden bubble in a self-rewetting fluid (SRF). In this figure, Re = 10, Ma = 0 25, . M 1 = 0, M 2 = 5, D/H = 0 1, . ψ L = 0 4, and . ψ R = 0 1. .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00146v2-figure-ce0206ff-f228-4f9b-86d2-84de42d118ec.jpg)

## 6.6 Effect of Marangoni number Ma

Finally, we investigate the effect of changing the Marangoni number Ma on the bubble migration. The baseline case (see Fig. 13) used Ma = 0 25 and we now consider two . additional cases, one with a smaller value (Ma = 0 1) and two larger values (Ma = 1 0 . . , 10 0) . than the above reference value. All the other parameters are kept fixed as before. According to Eq. (8), the Marangoni number Ma in the context of bubble migration is a dimensionless number that represents the effect of Marangoni convection on bubble motion relative to thermal diffusion. A larger Ma thus means greater effect of the convective motion compared to the thermal diffusion. As a result, when comparing Fig. 13 (Ma = 0 25) and Fig. 28 .

̸

(Ma = 0 1) with Fig. 29 (Ma = 1 0) and Fig. 30 (Ma = 10 0), it is seen that the bubble . . . undergoes rapid acceleration and deceleration during the initial transients especially at higher Ma, while by contrast, the presence of surfactants or β = 0 has a greater role at lower Ma. However, as in the case with varying Re, changing Ma does not significantly influence the bubble equilibrium location unlike the self-rewetting surface tension sensitivity coefficient on temperature M 2 or the surfactant parameters.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00146v2-figure-f0816c28-b3f7-47cf-85c5-7540e31cb8dc.jpg)

(b)

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00146v2-figure-d5c5f798-64ef-4180-bf83-8978a1f09975.jpg)

Bubble Centroid Location

Fig. 28. Bubble velocity( a ) and its location ( b ) at different values of the Gibbs elasticity parameter β for surfactant-laden bubble in a self-rewetting fluid (SRF). In this figure, Re = 5, Ma = 0 1, . M 1 = 0, M 2 = 5, D/H = 0 1, . ψ L = 0 4, and . ψ R = 0 1. .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00146v2-figure-fc533736-daa2-4f38-85e5-b69bb7f69772.jpg)

Fig. 29. Bubble velocity( a ) and its location ( b ) at different values of the Gibbs elasticity parameter β for surfactant-laden bubble in a self-rewetting fluid (SRF). In this figure, Re = 5, Ma = 1, M 1 = 0, M 2 = 5, D/H = 0 1, . ψ L = 0 4, and . ψ R = 0 1. .

Fig. 30. Bubble velocity( a ) and its location ( b ) at different values of the Gibbs elasticity parameter β for surfactant-laden bubble in a self-rewetting fluid (SRF). In this figure, Re = 5, Ma = 10, M 1 = 0, M 2 = 5, D/H = 0 1, . ψ L = 0 4, and . ψ R = 0 1. .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00146v2-figure-47a3bc68-8902-4a30-9910-dc781d24df79.jpg)

## 7 Summary and Conclusions

Surface tension, a property generally influenced by temperature, plays a vital role in interfacial transport phenomena within fluids. While normal fluids (NFs) follow a linear relationship between surface tension and temperature, self-rewetting fluids (SRFs) display atypical non-linear (quadratic) dependency on temperature. This relationship is characterized by a quadratic dependence with a minimum point and a positive gradient. Consequently, SRFs possess distinctive features, such as interfacial fluid motion towards higher temperature regions. These unique characteristics make SRFs highly promising for various applications, including microfluidics, both in microgravity conditions and terrestrial settings. On the other hand, however, the inclusion of surface-active agents or surfactants result in their adsorption on fluid interfaces and can significantly alter the behavior of the latter. The combined effects of SRFs and surfactants on bubble migration were not explored previously. Thus, the primary objective of this work is to simulate the thermocapillary migration of bubbles laden with surfactants in SRFs. The computational approach based on LBM using central moment was initially validated against some benchmark problems to ensure its validity.

̸

Simulations resulted in a number of interesting conclusions. It is found that the bubble migrates in SRF towards the minimum temperature location, and it equilibrates at the reference temperature location located at the middle of the domain for the clean interface case where the Gibbs elasticity parameter β = 0. Additionally, including surfactants ( β = 0) is shown to slow the bubble down and move the equilibrium location upstream toward a higher surfactant gradient location. Also, the terminal bubble location in a SRF at long times can be shifted by increasing β . As a result, β can be used to control the strength of the surfactant. The equilibrium position of the bubble is determined by the competing tangential Marangoni stresses due to thermocapillarity in SRFs and from the addition of surfactants. On the other hand, in NF, there is no equilibrium position for the bubble,

which migrates completely from the cold side to the hot side.

Upon analyzing the influence of the other parameters associated with the surface tension, it was found that variations in the value of M 2 or the quadratic surface tension sensitivity coefficient on temperature, significantly impact the dominant thermocapillary force, leading to observable alterations in the velocity and equilibrium position of the bubble. When higher values of M 2 are present, the bubble experiences a more pronounced surface tension gradient along its path. As a result, the bubble undergoes greater acceleration, affecting both its velocity and final equilibrium position. Moreover, the surfactants exert greater effect in shifting the equilibrium position at lower M 2 . Finally, the streamwise gradient in the imposed surfactant concentration field modulates the transient velocity before attaining the final resting state of the bubble. These findings provide potentially new pathways to manipulate the dynamics of bubbles in microchannels and other applications by exploiting the combined effects of surfactants and anomalous thermocapillary behaviour in SRFs.

## Acknowledgements

The authors would like to acknowledge the support of the US National Science Foundation (NSF) for research under Grant CBET-1705630. The third author (KNP) would also like to thank the NSF for support of the development of a computer cluster 'Alderaan' hosted at the Center for Computational Mathematics at the University of Colorado Denver under Grant OAC-2019089 (Project 'CC* Compute: Accelerating Science and Education by Campus and Grid Computing'), which was used in performing the simulations.

## A List of Symbols

## Nomenclature

L domain length

H domain height

D bubble diameter

T H hot temperature

T C cold temperature

T ref reference temperature

∆ T = T H - T C temperature difference

ψ dimensionless surfactant concentration (0 ≤ ψ ≤ 1)

k thermal conductivity

u velocity field

p pressure field

T temperature field

Re = UR/ν

a

the Reynolds number

Ma = UR/α

a the Marangoni number

Ca = Uµ /σ a 0

the capillary number

Eo = ρ U D/σ 2

a 0 the Eotvos number

U migration velocity of the bubble

U o characteristic thermocapillary velocity

## Greek Symbols:

- ρ fluid density
- µ dynamic viscosity
- α thermal diffusivity of the fluid
- ν = µ/ρ

kinematic viscosity of the fluid

- ˜ = ρ ρ /ρ a b density ratio
- ˜ = µ µ /µ a b dynamic viscosity ratio
- ˜ = k k /k a b thermal conductivity ratio
- β Gibbs elasticity parameter
- σ 0 the value of the surface tension at a reference temperature T ref
- σ T surface tension linear sensitivity coefficient
- M 1 = (∆ T/σ 0 ) σ T dimensionless surface tension linear sensitivity coefficient
- M 2 = ∆ ( T 2 /σ 0 ) σ TT dimensionless surface tension quadratic sensitivity coefficient
- M ψ local surfactant mobility
- m ψ the scale for the mobility parameter µ ψ chemical potential
- ϕ
- order parameter of the phase field variable
- λ , s , and w model parameters
- Pi indicates the proportionate contribution of surfactant diffusion
- Ex represents the relative intensity of adsorption and solubility effects
- M ϕ local mobility for interface tracking (in LBM)
- e α particle velocity directions (in LBM)
- f α distribution function for interface tracking (in LBM)
- g α distribution function for two-fluid flow (in LBM)
- h α distribution function for temperature field (in LBM)
- q α distribution function for surfactant concentration field (in LBM)

## Subscripts:

- A ambient fluid
- B bubble
- L left
- σ TT surface tension quadratic sensitivity coefficient
- R right

## B LBM for phase-field based interface capturing

In order to solve the conservative ACE provided in Eq. (12), we will now discuss a central moment LB approach that involves evolving a distribution function f α on the D2Q9 lattice, where α = 0 1 2 , , , . . . , 8 denotes the discrete particle directions. During the collision, the distribution functions f = ( f , f 0 1 , f 2 , . . . , f 8 ) † generally rela-x to the corresponding equilibrium distribution functions given by f eq = ( f eq 0 , f eq 1 , f eq 2 , . . . , f eq 8 ) † . These distribution functions must then be implemented via their central moments in the following.

First, in this regard, the following vectors, expressed in standard Dirac's bra-ket nota-

tion, may be used to describe the components of the particle velocities of this lattice as

$$| e _ { x } \rangle = ( 0, 1, 0, - 1, 0, 1, - 1, - 1, 0 ) ^ { \dagger },$$

$$| e _ { y } \rangle = ( 0, 0, 1, 0, - 1, 1, 1, - 1, - 1 ) ^ { \dagger }.$$

We also need the following 9-dimensional vector to define the zeroth moment of f α :

$$| \mathbf 1 \rangle = ( 1, 1, 1, 1, 1, 1, 1, 1, 1 ) ^ { \dagger }.$$

In other words, the order parameter ϕ of the phase-field model should result from its inner product with a set of distribution functions ⟨ f | 1 ⟩ . Next, the central moment LB method will be built using the nine non-orthogonal basis vectors shown below (which diffrents from methodology described in [77]):

$$\left | P _ { 0 } \right \rangle = \left | \mathbf 1 \right \rangle, \ \left | P _ { 1 } \right \rangle = \left | e _ { x } \right \rangle, \ \left | P _ { 2 } \right \rangle = \left | e _ { y } \right \rangle,$$

$$| P _ { 3 } \rangle = | e _ { x } ^ { 2 } + e _ { y } ^ { 2 } \rangle \,, \ \ | P _ { 4 } \rangle = | e _ { x } ^ { 2 } - e _ { y } ^ { 2 } \rangle \,, \ \ | P _ { 5 } \rangle = | e _ { x } e _ { y } \rangle \,,$$

$$| P _ { 6 } \rangle = \left | e _ { x } ^ { 2 } e _ { y } \right \rangle, \ \ | P _ { 7 } \rangle = \left | e _ { x } e _ { y } ^ { 2 } \right \rangle, \ \ | P _ { 8 } \rangle = \left | e _ { x } ^ { 2 } e _ { y } ^ { 2 } \right \rangle.$$

Symbols such as ∣ ∣ e 2 x e y 〉 = | e x e x e y ⟩ denote a vector that results from the element-wise vector multiplication of vectors | e x ⟩ | , e x ⟩ and | e y ⟩ . Using the moment basis vectors mentioned above, they may be arranged into the following matrix, which maps the distribution functions to the raw moments:

$$\mathbf P = \left [ \langle P _ { 0 } | \,, \langle P _ { 1 } | \,, \langle P _ { 2 } | \,, \langle P _ { 3 } | \,, \langle P _ { 4 } | \,, \langle P _ { 5 } | \,, \langle P _ { 6 } | \,, \langle P _ { 7 } | \,, \langle P _ { 8 } | \, \right ].$$

It should be emphasized that by shifting the particle velocity e α by the fluid velocity u , the central moments are obtained from the distribution moments. Using this, we can then define the raw moments of the distribution function and equilibrium f eq α explicitly as

$$\left ( \begin{array} { c } \kappa _ { m n } ^ { \prime } \\ \kappa _ { m n } ^ { \prime \, e q } \end{array} \right ) = \sum _ { \alpha = 0 } ^ { 8 } \left ( \begin{array} { c } f _ { \alpha } \\ f _ { \alpha } ^ { e q } \end{array} \right ) e _ { \alpha x } ^ { m } e _ { \alpha y } ^ { n },$$

and the corresponding central moments as

$$\left ( \begin{array} { c } \kappa _ { m n } \\ \kappa _ { m n } ^ { e q } \end{array} \right ) = \sum _ { \alpha = 0 } ^ { 8 } \left ( \begin{array} { c } f _ { \alpha } \\ f _ { \alpha } ^ { e q } \end{array} \right ) ( e _ { \alpha x } - u _ { x } ) ^ { m } ( e _ { \alpha y } - u _ { y } ) ^ { n }.$$

The raw moment of order ( m + ) is thus represented by n κ ′ mn , and the corresponding central moment is κ mn . The two vectors that follow can be used to easily group together all the possible raw moments and central moments for the D2Q9 lattice:

$$\kappa ^ { ^ { \prime } } = ( \kappa ^ { ^ { \prime } } _ { 0 0 }, \kappa ^ { ^ { \prime } } _ { 1 0 }, \kappa ^ { ^ { \prime } } _ { 0 1 }, \kappa ^ { ^ { \prime } } _ { 2 0 }, \kappa ^ { ^ { \prime } } _ { 0 2 }, \kappa ^ { ^ { \prime } } _ { 1 1 }, \kappa ^ { ^ { \prime } } _ { 2 1 }, \kappa ^ { ^ { \prime } } _ { 1 2 }, \kappa ^ { ^ { \prime } } _ { 2 2 } ),$$

$$\kappa \, = \, ( \kappa _ { 0 0 }, \kappa _ { 1 0 }, \kappa _ { 0 1 }, \kappa _ { 2 0 }, \kappa _ { 0 2 }, \kappa _ { 1 1 }, \kappa _ { 2 1 }, \kappa _ { 1 2 }, \kappa _ { 2 2 } ).$$

It should be noted that one can readily map from the distribution functions to the raw moments via κ ′ = P f , which can then be transformed into the central moments through κ = F κ ′ , where the frame transformation matrix F follows readily from binomial expansions of ( e αx -u x ) m ( e αy -u y ) n to relate to e m x α e n αy etc. Similarly, the inverse mappings from central moments to raw moments, from which the distribution functions can be recovered via the matrices F -1 and P -1 , respectively. All these mapping relations are explicitly listed in Appendix F.

As previously indicated, an important aspect of our approach is carrying out the collision step in a way that causes those different central moments shown above to relax to their their corresponding central moment equilibria. The discrete central moment equilibria κ eq mn defined above can be obtained by matching them to the corresponding central moments of the continuous Maxwell distribution function after replacing the density ρ with the order parameter ϕ ; Moreover, M θn ϕ x and M θn ϕ y [77,90] must be added to the first order central moment equilibrium components in order to account for the interface sharpening flux terms in the conservative ACE (Eq. (12)). Therefore, we have

$$\kappa _ { 0 0 } ^ { e q } & = \phi, \quad \kappa _ { 1 0 } ^ { e q } = M _ { \phi } \theta n _ { x }, \quad \kappa _ { 0 1 } ^ { e q } = M _ { \phi } \theta n _ { y }, \\ \kappa _ { 2 0 } ^ { e q } & = c _ { s \phi } ^ { 2 } \phi, \quad \kappa _ { 0 2 } ^ { e q } = c _ { s \phi } ^ { 2 } \phi, \quad \kappa _ { 1 1 } ^ { e q } = 0, \\ \kappa _ { 2 1 } ^ { e q } & = 0, \quad \kappa _ { 1 2 } ^ { e q } = 0, \quad \kappa _ { 2 2 } ^ { e q } = c _ { s \phi } ^ { 4 } \phi, \\. & \quad \partial _ { x x }.$$

where c 2 sϕ = 1 3. /

We can now summarize the central moment LB algorithm for solving the conservative ACE for a time step ∆ , starting from t f α = f α ( x , t ), as follows, based on the above considerations, which were inspired by the algorithmic implementation presented in [106] (see also [99,100] for its further extensions):

- · Compute pre-collision raw moments from distribution functions via κ ′ = P f (see Eq. (45) in Appendix F for P )
- · Compute pre-collision central moments from raw moments via κ = F κ ′ (see Eq. (46) in Appendix F for F )
- · Perform collision step via relaxation of central moments κ mn to their equilibria κ eq mn :

$$\tilde { \kappa } _ { m n } = \kappa _ { m n } + \omega _ { m n } ^ { \phi } ( \kappa _ { m n } ^ { e q } - \kappa _ { m n } ),$$

where ( mn ) = (00) (10) (01) (20) (02) (11) (21) (12), and (22), and , , , , , , , ω ϕ mn is the relaxation parameter for moment of order ( m + n ). Here, the implicit summation convention of repeated indices is not assumed. The relaxation parameters of the first order moments ω ϕ 10 = ω ϕ 01 = ω ϕ are related to the mobility coefficient M ϕ in Eq. (12) via M ϕ = c 2 sϕ ( 1 ω ϕ -1 2 ) ∆ , and the rest of the relaxation parameters are typically set t to unity, i.e., ω ϕ mn = 1 0, where ( . m + ) n ≥ 2. The results of Eq. (27) are then grouped in ˜ κ .

- · Compute post-collision raw moments from post-collision central moments via ˜ κ ′ = F -1 ˜ κ (see Eq. (47) in Appendix F for F -1 )

- · Compute post-collision distribution functions from post-collision raw moments via ˜ f = P -1 ˜ κ ′ (see Eq. (48) in Appendix F for P -1 )
- · Perform streaming step via f α ( x , t +∆ ) = ˜ ( t f α x -e α ∆ ), where t α = 0 1 2 , , , ..., 8.
- · Update the order parameter ϕ of the phase-field model for interface capturing through

$$\phi = \sum _ { \alpha = 0 } ^ { 8 } f _ { \alpha }.$$

## C LBMfor two-fluid flow with capillary and Marangoni forces

The motion of binary fluids with interfacial forces represented in Eqs. (13)-(14) will then be solved using a central moment LBM by evolving another distribution function g α , where α = 0 1 2 , , , . . . , 8. Our approach is based on discretizing the modified continuous Boltzmann equation and using a matching principle with their continuous counterparts to obtain the discrete central moment equilibria and central moments of the source terms for the body forces, as explained in Ref. [77]. Nevertheless, in the following, we consider the simpler, non-orthogonal moment basis vectors as provided previously in Eq. (23), in contrast to Ref. [77,90], where an orthogonal moment basis is applied resulting in the so-called cascaded LB method.

The distribution function g α , its equilibrium g eq α , and the source term S α are defined first, as in the previous section. The latter accounts for the surface tension, body forces, and those resulting from the application of a transformation to simulate flows at high-density ratios in the incompressible limit (see [73,77]).

$$\left ( \begin{array} { c } \eta _ { m n } ^ { \prime } \\ \eta _ { m n } ^ { \prime \ e q } \\ \sigma _ { m n } ^ { \prime } \end{array} \right ) = \sum _ { \alpha = 0 } ^ { 8 } \left ( \begin{array} { c } g _ { \alpha } \\ g _ { \alpha } ^ { \ e q } \\ S _ { \alpha } \end{array} \right ) e _ { \alpha x } ^ { m } e _ { \alpha y } ^ { n },$$

$$\left ( \begin{array} { c } \eta _ { m n } \\ \eta _ { m n } ^ { e q } \\ \sigma _ { m n } \end{array} \right ) = \sum _ { \alpha = 0 } ^ { 8 } \left ( \begin{array} { c } g _ { \alpha } \\ g _ { \alpha } ^ { e q } \\ S _ { \alpha } \end{array} \right ) ( e _ { \alpha x } - u _ { x } ) ^ { m } ( e _ { \alpha y } - u _ { y } ) ^ { n }.$$

For convenience, we can group the source term, equilibrium, and components of the distribution function for the D2Q9 lattice as the following vectors: g eq = ( g eq 0 , g eq 1 , g eq 2 , . . . , g eq 8 ) † , S = ( S , S 0 1 , S 2 , . . . , g 8 ) † . Additionally, we use the following to group the possible raw moments and the central moments already defined above for the D2Q9 lattice:

$$\eta ^ { \prime } & = ( \eta ^ { ^ { \prime } } _ { 0 0 }, \eta ^ { ^ { \prime } } _ { 1 0 }, \eta ^ { ^ { \prime } } _ { 0 1 }, \eta ^ { ^ { \prime } } _ { 2 0 }, \eta ^ { ^ { \prime } } _ { 0 2 }, \eta ^ { ^ { \prime } } _ { 1 1 }, \eta ^ { ^ { \prime } } _ { 2 1 }, \eta ^ { ^ { \prime } } _ { 1 2 }, \eta ^ { ^ { \prime } } _ { 2 2 } ),$$

$$\eta = ( \eta _ { 0 0 }, \eta _ { 1 0 }, \eta _ { 0 1 }, \eta _ { 2 0 }, \eta _ { 0 2 }, \eta _ { 1 1 }, \eta _ { 2 1 }, \eta _ { 1 2 }, \eta _ { 2 2 } ),$$

and similarly for raw moments and the central moments the equilibrium and the source term.

The collision step will be performed such that different central moments shown above relax to their corresponding central moment equilibria, which are augmented by changes in the central moments due to the net forces; the latter is given by sum the surface tension force F s = ( F sx , F sy ), which can have contributions from both the capillary and Marangoni forces as represented in Eq. (16), and any external force F ext = ( F ext,x , F ext,y ), i.e., F t = F s + F ext or ( F tx , F ty ) = ( F sx + F ext,x , F sy + F ext,y ). Moreover, the use of an incompressible transformation as mentioned above leads to a pressure-based formulation, involving the incorporation of a net pressure force F p arising from φ ρ ( ) = p -ρc 2 s , i.e., F p = -∇ φ , or ( F px , F py ) = ( -∂ φ, x -∂ φ y ) (see [77] for details). Then, the discrete central moment equilibria η mn defined above can be obtained by matching them to the corresponding continuous central moments of the equilibrium that arise from the incompressible transformation, and similarly for the central moments of the source term σ mn , which then results in the following expressions for the D2Q9 lattice [77]:

$$\eta _ { 0 0 } ^ { e q } & = p, \quad \eta _ { 1 0 } ^ { e q } = - \varphi ( \rho ) u _ { x }, \quad \eta _ { 0 1 } ^ { e q } = - \varphi ( \rho ) u _ { y }, \quad \eta _ { 2 0 } ^ { e q } = p c _ { s } ^ { 2 } + \varphi ( \rho ) u _ { x } ^ { 2 }, \\ \eta _ { 0 2 } ^ { e q } & = p c _ { s } ^ { 2 } + \varphi ( \rho ) u _ { y } ^ { 2 }, \quad \eta _ { 1 1 } ^ { e q } = \varphi ( \rho ) u _ { x } u _ { y }, \quad \eta _ { 2 1 } ^ { e q } = - \varphi ( \rho ) ( u _ { x } ^ { 2 } + c _ { s } ^ { 2 } ) u _ { y }, \\ \eta _ { 1 2 } ^ { e q } & = - \varphi ( \rho ) ( u _ { y } ^ { 2 } + c _ { s } ^ { 2 } ) u _ { x }, \quad \eta _ { 2 2 } ^ { e q } = c _ { s } ^ { 6 } \rho + \varphi ( \rho ) ( u _ { x } ^ { 2 } + c _ { s } ^ { 2 } ) ( u _ { y } ^ { 2 } + c _ { s } ^ { 2 } ).$$

and

$$\sigma _ { 0 0 } & = \Gamma _ { 0 0 } ^ { p }, \quad \sigma _ { 1 0 } = c _ { s } ^ { 2 } F _ { t x } - u _ { x } \Gamma _ { 0 0 } ^ { p }, \quad \sigma _ { 0 1 } = c _ { s } ^ { 2 } F _ { t y } - u _ { y } \Gamma _ { 0 0 } ^ { p }, \\ \sigma _ { 2 0 } & = 2 c _ { s } ^ { 2 } F _ { p x } u _ { x } + ( u _ { x } ^ { 2 } + c _ { s } ^ { 2 } ) \Gamma _ { 0 0 } ^ { p }, \quad \sigma _ { 0 2 } = 2 c _ { s } ^ { 2 } F _ { p y } u _ { y } + ( u _ { y } ^ { 2 } + c _ { s } ^ { 2 } ) \Gamma _ { 0 0 } ^ { p }, \\ \sigma _ { 1 1 } & = c _ { s } ^ { 2 } ( F _ { p x } u _ { y } + F _ { p y } u _ { x } ) + u _ { x } u _ { y } \Gamma _ { 0 0 } ^ { p }, \quad \sigma _ { 2 1 } = 0, \quad \sigma _ { 1 2 } = 0, \quad \sigma _ { 2 2 } = 0,$$

where Γ p 00 = ( F px u x + F py u y ).

We can now summarize the central moment LB algorithm for computing the two-fluid motion with interfacial forces for a time step ∆ t starting from g α = g α ( x , t ) using the above developments as follows:

- · Compute pre-collision raw moments from distribution functions via η ′ = P g (see Eq. (45) in Appendix F for P )
- · Compute pre-collision central moments from raw moments via η = F η ′ (see Eq. (46) in Appendix F for F )
- · Perform collision step via relaxation of central moments η mn to their equilibria η eq mn and augmented with the source terms σ mn :

In order to allow for an independent specification of the shear viscosity ν from the bulk viscosity ζ , the trace of the second order moments η 20 + η 02 should be evolved independently from the other second-order moments. To accomplish this, prior to collision, we combine the diagonal parts of the second-order moments as follows (see e.g., [99, 100, 106]):

$$\eta _ { 2 s } & = \eta _ { 2 0 } + \eta _ { 0 2 }, \quad \eta _ { 2 s } ^ { e g } = \eta _ { 2 0 s } ^ { e g } + \eta _ { 0 2 } ^ { e g }, \quad \sigma _ { 2 s } = \sigma _ { 2 0 s } + \sigma _ { 0 2 }, \\ \eta _ { 2 d } & = \eta _ { 2 0 } - \eta _ { 0 2 }, \quad \eta _ { 2 d } ^ { e g } = \eta _ { 2 0 s } ^ { e g } - \eta _ { 0 2 } ^ { e g }, \quad \sigma _ { 2 d } = \sigma _ { 2 0 s } - \sigma _ { 0 2 },$$

and thus η 2 s and η 2 d will be evolved independently under collision. Then, the postcollision central moments under relaxation and augmentation due to the forces can be computed via

$$\tilde { \eta } _ { m n } = \eta _ { m n } + \omega _ { m n } \left ( \eta _ { m n } ^ { e q } - \eta _ { m n } \right ) + \left ( 1 - \omega _ { m n } / 2 \right ) \sigma _ { m n } \Delta t,$$

where ω mn is the relaxation time corresponding to the central moment η mn , and ( mn ) = (00) (10) (01) (2 , , , s , ) (2 d , ) (11) (21) , , (12) and (22). , , Here, the relaxation parameter ω 2 s is related to the bulk viscosity via ζ = c 2 s (1 /ω 2 s -1 2) ∆ , while the relax-/ t ation parameters ω 2 d and ω 11 are related to shear viscosity via ν = c 2 s (1 /ω ij -1 2) ∆ / t where ( ij ) = (2 ) (11). d , Typically, c 2 s = 1 3. / In view of Eq. (17) it should be noted that if the bulk fluid properties are different, the relaxation parameters ω 2 d and ω 11 will then vary locally across the interface. The rest of the relaxation parameters of central moments are generally set to unity, i.e., ω ij = 1 0, where ( . ij ) = (00) (10) , , (01) (2 , s , ) (21) (12) , , (22).

Also, the combined forms of the post-collision central moments ˜ η 2 s and ˜ η 2 d are then segregated in their individual components ˜ η 20 and ˜ η 02 via

$$\tilde { \eta } _ { 2 0 } = \frac { 1 } { 2 } \left ( \tilde { \eta } _ { 2 s } + \tilde { \eta } _ { 2 d } \right ), \quad \tilde { \eta } _ { 0 2 } = \frac { 1 } { 2 } \left ( \tilde { \eta } _ { 2 s } - \tilde { \eta } _ { 2 d } \right ).$$

Finally, the results of Eq. (33) by accounting for the above segregation are then grouped in ˜ η .

- · Compute post-collision raw moments from post-collision central moments via ˜ η ′ = F -1 ˜ η (see Eq. (47) in Appendix F for F -1 )
- · Compute post-collision distribution functions from post-collision raw moments via ˜ g = P -1 ˜ η ′ (see Eq. (48) in Appendix F for P -1 )
- · Perform streaming step via g α ( x , t +∆ ) = ˜ ( t g α x -e α ∆ ), where t α = 0 1 2 , , , ..., 8.
- · Update the pressure field p and the components of the fluid velocity u = ( u , u x y ) through

$$p = \sum _ { \alpha } g _ { \alpha } + \frac { 1 } { 2 } F _ { p } \cdot u \Delta t, \quad \rho c _ { s } ^ { 2 } u = \sum _ { \alpha } g _ { \alpha } e _ { \alpha } + \frac { 1 } { 2 } c _ { s } ^ { 2 } F _ { t } \Delta t.$$

## D LBM for energy transport equation

The energy transport equation (Eq. (1c)) can be solved using a central moment LB technique by evolving a third distribution function h α , where α = 0 1 2 , , , . . . , 8, on the D2Q9 lattice. As an advection-diffusion equation, Eq. (1c) is constructed in a manner similar to the LB scheme for the conservative ACE previously discussed, but without including a term like the interface sharpening flux term that is included in the latter case. As previously, we first

establish the following raw moments and central moments, respectively, of the distribution function h α , as well as its equilibrium h eq α :

$$\left ( \begin{array} { c } \chi ^ { \prime } _ { m n } \\ \chi ^ { \prime \, e q } _ { m n } \end{array} \right ) = \sum _ { \alpha = 0 } ^ { 8 } \left ( \begin{array} { c } h _ { \alpha } \\ h _ { \alpha } ^ { e q } \end{array} \right ) e _ { \alpha x } ^ { m } e _ { \alpha y } ^ { n },$$

$$\left ( \begin{array} { c } \chi _ { m n } \\ \chi _ { m n } ^ { e q } \end{array} \right ) = \sum _ { \alpha = 0 } ^ { 8 } \left ( \begin{array} { c } h _ { \alpha } \\ h _ { \alpha } ^ { e q } \end{array} \right ) ( e _ { \alpha x } - u _ { x } ) ^ { m } ( e _ { \alpha y } - u _ { y } ) ^ { n }.$$

For convenience, we list the components of the distribution function and its equilibrium, respectively, using h = ( h , h 0 1 , h 2 , . . . , h 8 ) † and h eq = ( h eq 0 , h eq 1 , h eq 2 , . . . , h eq 8 ) † , and analogously for the raw moments and central moments via

$$\chi ^ { ^ { \prime } } = ( \chi ^ { ^ { \prime } } _ { 0 0 }, \chi ^ { ^ { \prime } } _ { 1 0 }, \chi ^ { ^ { \prime } } _ { 0 1 }, \chi ^ { ^ { \prime } } _ { 2 0 }, \chi ^ { ^ { \prime } } _ { 0 2 }, \chi ^ { ^ { \prime } } _ { 1 1 }, \chi ^ { ^ { \prime } } _ { 2 1 }, \chi ^ { ^ { \prime } } _ { 1 2 }, \chi ^ { ^ { \prime } } _ { 2 2 } ),$$

$$\chi = ( \chi _ { 0 0 }, \chi _ { 1 0 }, \chi _ { 0 1 }, \chi _ { 2 0 }, \chi _ { 0 2 }, \chi _ { 1 1 }, \chi _ { 2 1 }, \chi _ { 1 2 }, \chi _ { 2 2 } ).$$

Similar to Sec. B, to construct a central moment-based collision model for solving the energy equation, we obtain the discrete equilibrium central moments from the corresponding continuous counterpart of the Maxwellian by replacing the density ρ with the temperature T , and the results read as

$$\chi _ { 0 0 } ^ { e q } & = T, \quad \chi _ { 1 0 } ^ { e q } = 0, \quad \chi _ { 0 1 } ^ { e q } = 0, \\ \chi _ { 2 0 } ^ { e q } & = c _ { s T } ^ { 2 } T, \quad \chi _ { 0 2 } ^ { e q } = c _ { s T } ^ { 2 } T, \quad \chi _ { 1 1 } ^ { e q } = 0, \\ \chi _ { 2 1 } & = 0, \quad \chi _ { 1 2 } ^ { e q } = 0, \quad \chi _ { 2 2 } ^ { e q } = c _ { s T } ^ { 4 } T, \\ & \wedge \quad.$$

where, typically, c 2 sT = 1 3. / Then, the computational procedure for solving the energy equation for a time step ∆ t starting from h α = h α ( x , t ) can be summarized as follows:

- · Compute pre-collision raw moments from distribution functions via χ ′ = P h (see Eq. (45) in Appendix F for P )
- · Compute pre-collision central moments from raw moments via χ = F χ ′ (see Eq. (46) in Appendix F for F )
- · Perform collision step via relaxation of central moments χ mn to their equilibria χ eq mn :

$$\tilde { \chi } _ { m n } = \chi _ { m n } + \omega _ { m n } ^ { T } ( \chi _ { m n } ^ { e q } - \chi _ { m n } ),$$

where ( mn ) = (00) (10) (01) (20) (02) (11) (21) (12), and (22), and , , , , , , , ω T mn is the relaxation parameter for moment of order ( m + n ). The relaxation parameters of the first order moments ω T 10 = ω T 01 = ω T are related to the thermal diffusivity α = k/ ρc ( p ) via α = c 2 sT ( 1 /ω T -1 2 / ) ∆ , and the rest of the relaxation parameters of higher cent tral moments are typically set to unity. The results of Eq. (38) are then grouped in ˜ χ .

- · Compute post-collision raw moments from post-collision central moments via ˜ χ ′ = F -1 ˜ χ (see Eq. (47) in Appendix F for F -1 )
- · Compute post-collision distribution functions from post-collision raw moments via ˜ h = P -1 ˜ χ ′ (see Eq. (48) in Appendix F for P -1 )
- · Perform streaming step via h α ( x , t +∆ ) = ˜ t h α ( x -e α ∆ ), where t α = 0 1 2 , , , ..., 8.
- · Update the temperature field T is obtained from

$$T = \sum _ { \alpha = 0 } ^ { 8 } h _ { \alpha }.$$

## E LBM for surfactant concentration equation

The surfactant concentration equation (Eq. (4)), together with Eq. (5), is solved using a central moment LB technique by evolving a fourth distribution function q α , where α = 0 1 2 , , , . . . , 8, on the D2Q9 lattice. As before, we begin by defining the distribution function q α and its equilibrium q eq α , as well as the following raw and central moments, respectively:

$$\left ( \begin{array} { c } \Upsilon _ { m n } ^ { \prime } \\ \Upsilon _ { m n } ^ { \prime \, e q } \end{array} \right ) = \sum _ { \alpha = 0 } ^ { 8 } \left ( \begin{array} { c } q _ { \alpha } \\ q _ { \alpha } ^ { e q } \end{array} \right ) e _ { \alpha x } ^ { m } e _ { \alpha y } ^ { n },$$

$$\left ( \begin{array} { c } \Upsilon _ { m n } \\ \Upsilon _ { m n } ^ { e q } \end{array} \right ) = \sum _ { \alpha = 0 } ^ { 8 } \left ( \begin{array} { c } q _ { \alpha } \\ q _ { \alpha } ^ { e q } \end{array} \right ) ( e _ { \alpha x } - u _ { x } ) ^ { m } ( e _ { \alpha y } - u _ { y } ) ^ { n }.$$

For convenience, we list the components of the distribution function and its equilibrium, respectively, using q = ( q , q 0 1 , q 2 , . . . , q 8 ) † and q eq = ( q eq 0 , q eq 1 , q eq 2 , . . . , q eq 8 ) † , and analogously for the raw moments and central moments via

$$\mathfrak { Y } ^ { \prime } & = ( \mathfrak { Y } ^ { \prime } _ { 0 0 }, \mathfrak { Y } ^ { \prime } _ { 1 0 }, \mathfrak { Y } ^ { \prime } _ { 0 1 }, \mathfrak { Y } ^ { \prime } _ { 2 0 }, \mathfrak { Y } ^ { \prime } _ { 0 2 }, \mathfrak { Y } ^ { \prime } _ { 1 1 }, \mathfrak { Y } ^ { \prime } _ { 2 1 }, \mathfrak { Y } ^ { \prime } _ { 1 2 }, S ^ { \prime } _ { 2 2 } ),$$

$$\mathfrak { Y } & = ( \mathfrak { Y } _ { 0 0 }, \mathfrak { Y } _ { 1 0 }, \mathfrak { Y } _ { 0 1 }, \mathfrak { Y } _ { 2 0 }, \mathfrak { Y } _ { 0 2 }, \mathfrak { Y } _ { 1 1 }, \mathfrak { Y } _ { 2 1 }, \mathfrak { Y } _ { 1 2 }, \mathfrak { Y } _ { 2 2 } ).$$

Similar to Section. B, we obtain the discrete equilibrium central moments from the corresponding continuous counterpart of the Maxwellian by substituting the surfactant concentration ψ for the density ρ . This allows us to build a central moment-based collision model for solving the advection and diffusion parts of the surfactant concentration equation given in Eq. (4); moreover, we account for the fluxes due to R = ( R ,R x y ) via m ψ ψ (1 -ψ ) R appearing in Eq. (4) by means of further corrections to the first order equilibrium moments. The outcomes read as follows.

$$\Upsilon _ { 0 0 } ^ { e q } & = \psi, \quad \Upsilon _ { 1 0 } ^ { e q } = - M _ { \psi } \psi ( 1 - \psi ) R _ { x }, \\ \Upsilon _ { 0 1 } ^ { e q } & = - M _ { \psi } \psi ( 1 - \psi ) R _ { y }, \quad \Upsilon _ { 2 0 } ^ { e q } = c _ { s _ { \psi } } ^ { 2 } \psi, \quad \Upsilon _ { 0 2 } ^ { e q } = c _ { s _ { \psi } } ^ { 2 } \psi, \\ \Upsilon _ { 1 1 } ^ { e q } & = 0, \quad \Upsilon _ { 2 1 } ^ { e q } = 0, \quad \Upsilon _ { 2 2 } ^ { e q } = c _ { s _ { \psi } } ^ { 4 } \psi,$$

$$\Upsilon _ { 1 1 } ^ { e q } = 0, \quad \Upsilon _ { 2 1 } = 0, \quad \Upsilon _ { 1 2 } ^ { e q } = 0, \quad \Upsilon _ { 2 2 } ^ { e q } = c _ { s _ { \psi } } ^ { 4 } \psi,$$

where c 2 s ψ = 1 3 and / R = ( R ,R x y ) is defined in Eq. (5). Then, the computational procedure for solving the surfactant concentration equation for a time step ∆ t starting from q α = q α ( x , t ) can be summarized as follows:

- · Compute pre-collision raw moments from distribution functions via Υ ′ = P q (see Eq. (45) in Appendix F for P )
- · Compute pre-collision central moments from raw moments via Υ = F Υ ′ (see Eq. (46) in Appendix F for F )
- · Perform collision step via relaxation of central moments Υ mn to their equilibria Υ eq mn :

$$\tilde { \Upsilon } _ { m n } = \Upsilon _ { m n } + \omega _ { m n } ^ { \psi } ( \Upsilon _ { m n } ^ { e q } - \Upsilon _ { m n } ),$$

where ( mn ) = (00) (10) (01) (20) (02) (11) (21) (12), , , , , , , , and (22), and ω ψ mn is the relaxation parameter for moment of order ( m + n ). The relaxation parameters of the first order moments ω ψ 10 = ω ψ 01 = ω ψ are related to the local surfactant mobility M ψ = m ψ ψ (1 -ψ ) via m ψ = c 2 s ψ ( 1 /ω ψ -1 2 / ) ∆ , and the rest of the relaxation pat rameters of higher central moments are typically set to unity. The results of Eq. (43) are then grouped in ˜ Υ .

- · Compute post-collision raw moments from post-collision central moments via ˜ Υ ′ = F -1 ˜ Υ (see Eq. (47) in Appendix F for F -1 )
- · Compute post-collision distribution functions from post-collision raw moments via ˜ q = P -1 ˜ Υ ′ (see Eq. (48) in Appendix F for P -1 )
- · Perform streaming step via q α ( x , t +∆ ) = ˜ ( t q α x -e α ∆ ), where t α = 0 1 2 , , , ..., 8.
- · Finally, update the surfactant concentration field ψ which is obtained from

$$\psi = \sum _ { \alpha = 0 } ^ { 8 } q _ { \alpha }.$$

## F Transformation Matrices for Central Moment LB Schemes on a D2Q9 lattice

Here, we summarize the various mapping relations that are needed prior to and following the collision step, where different central moments are relaxed to their equilibria, in the central moment LB scheme on the D2Q9 lattice. See Appendix B for related discussion and the notations adopted.

The transformation matrix P mapping a vector of distribution functions f to a vector

## of raw moments κ ′ is given by

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00146v2-figure-aa74f8b6-6f0b-4c81-9c3d-8940a5cf76e3.jpg)

Next, the transformation matrix F mapping a vector of raw moments κ ′ to a vector of central moments κ reads as



1

0

0

0

0

0

0

0

0























































0

0

0

0

0









































0





0

-

u

x

-

u

y

u

2

x

u

2

y

u u

x

y

-

u u

2

x

-

y

u u

x

y

-

1

0

2

u

0

-

u

y

2

u u

x

y

u

y

0

1

0

-

2

u

-

u

y

x

u

2

x

2

u u

x

y

0

0

1

0

0

-

u

0

0

0

0

1

0

0

-

u

x

0

0

0

0

1

-

2

u

-

x

2

u

y

0

0

0

0

0

1

0

0

0

0

0

0

0

1







2



u u

x

2

y

-

u u

x

2

y

-

u u

2

x

y

u

2

y

u

2

x

4

u u

x

y

-

2

u

y

-

2

u

x

1

Then, the transformation matrix F -1 mapping a vector of (post-collision) central moments

2

2

F

=

x

y

(46)

˜ κ to a vector of (post-collision) raw moments ˜ κ ′ can be written as

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00146v2-figure-2f00b2f1-1986-4cb7-b366-5e82c694515e.jpg)

It may be noted that if F = F ( u , u x y ), then F -1 = F ( -u , x -u y

Finally, we express the transformation matrix P raw moments ˜ κ ′ to a vector of (post-collision) distribution functions ˜ f as

) (see [99]). -1 mapping a vector of (post-collision)

$$P ^ { - 1 } = \begin{bmatrix} 1 & 0 & 0 & - 1 & - 1 & 0 & 0 & 0 & 1 \\ 0 & \frac { 1 } { 2 } & 0 & \frac { 1 } { 2 } & 0 & 0 & 0 & - \frac { 1 } { 2 } & - \frac { 1 } { 2 } \\ 0 & 0 & \frac { 1 } { 2 } & 0 & \frac { 1 } { 2 } & 0 & - \frac { 1 } { 2 } \\ 0 & - \frac { 1 } { 2 } & 0 & \frac { 1 } { 2 } & 0 & 0 & \frac { 1 } { 2 } & - \frac { 1 } { 2 } \\ 0 & 0 & - \frac { 1 } { 2 } & 0 & \frac { 1 } { 2 } & 0 & - \frac { 1 } { 2 } & - \frac { 1 } { 2 } \\ 0 & 0 & 0 & 0 & 0 & \frac { 1 } { 4 } & \frac { 1 } { 4 } & \frac { 1 } { 4 } \\ 0 & 0 & 0 & 0 & 0 & - \frac { 1 } { 4 } & \frac { 1 } { 4 } & - \frac { 1 } { 4 } & \frac { 1 } { 4 } \\ 0 & 0 & 0 & 0 & 0 & \frac { 1 } { 4 } & - \frac { 1 } { 4 } & - \frac { 1 } { 4 } & \frac { 1 } { 4 } \\ 0 & 0 & 0 & 0 & 0 & - \frac { 1 } { 4 } & - \frac { 1 } { 4 } & \frac { 1 } { 4 } \end{bmatrix}$$

## References

| [1]   | Pierre-Gilles De Gennes, Fran¸oise c Brochard-Wyart, David Qu´r´, e e et al. Capillarity and wetting phenomena: drops, bubbles, pearls, waves , volume 315. Springer, 2004.                                                                                   |
|-------|---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| [2]   | LE Scriven and CV Sternling. The Marangoni effects. Nature , 187(4733):186-188, 1960.                                                                                                                                                                         |
| [3]   | Ronald F Probstein. Physicochemical hydrodynamics: an introduction . John Wiley & Sons, 2005.                                                                                                                                                                 |
| [4]   | NO Young, JS Goldstein, and Mi J Block. The motion of bubbles in a vertical temperature gradient. Journal of Fluid Mechanics , 6(3):350-356, 1959.                                                                                                            |
| [5]   | RShankar Subramanian. Slow migration of a gas bubble in a thermal gradient. AIChE Journal , 27(4):646-654, 1981.                                                                                                                                              |
| [6]   | R Shankar Subramanian. Thermocapillary migration of bubbles and droplets. Ad- vances in Space Research , 3(5):145-153, 1983.                                                                                                                                  |
| [7]   | Samuel WJ Welch. Transient thermocapillary migration of deformable bubbles. Jour- nal of colloid and interface science , 208(2):500-508, 1998.                                                                                                                |
| [8]   | Chen Ma and Dieter Bothe. Direct numerical simulation of thermocapillary flow based on the volume of fluid method. International Journal of Multiphase Flow , 37(9):1045- 1058, 2011.                                                                         |
| [9]   | Anton A Darhuber and Sandra M Troian. Principles of microfluidic actuation by modulation of surface stresses. Annual Review of Fluid Mechanics , 37:425-455, 2005.                                                                                            |
| [10]  | Alireza Karbalaei, Ranganathan Kumar, and Hyoung Jin Cho. Thermocapillarity in microfluidics-a review. Micromachines , 7(1):13, 2016.                                                                                                                         |
| [11]  | R Vochten and Georges Petre. Study of the heat of reversible adsorption at the air- solution interface. ii. experimental determination of the heat of reversible adsorption of some alcohols. Journal of Colloid and Interface Science , 42(2):320-327, 1973. |
| [12]  | Georges Petre and Maherzia Aza Azouni. Experimental evidence for the minimum of surface tension with temperature at aqueous alcohol solution/air interfaces. Journal of Colloid and Interface Science , 98(1):261-263, 1984.                                  |
| [13]  | MC Limbourg-Fontaine, Georges P´tr´, e e and Jean Claude Legros. Thermocapillary movements under microgravity at a minimum of surface tension. Naturwissenschaften , 73(7):360-362, 1986.                                                                     |
| [14]  | D Villers and JK Platten. Temperature dependence of the interfacial tension between water and long-chain alcohols. The Journal of Physical Chemistry , 92(14):4023-4024,                                                                                      |

| [15]   | Yoshiyuki Abe, Akira Iwasaki, and Kotaro Tanaka. Microgravity experiments on phase change of self-rewetting fluids. Annals of the New York Academy of Sciences , 1027(1):269-285, 2004.                                                                                                                                        |
|--------|--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| [16]   | William Schupbach, Bashir Elbousefi, and Kannan Premnath. Lattice Boltzmann simulations of thermocapillary convection in self-rewetting fluid layers and bubble dynamics. Bulletin of the American Physical Society , 2022.                                                                                                    |
| [17]   | Bashir Elbousefi, William Schupbach, and Kannan Premnath. Lattice Boltzmann simulations of self-rewetting drops impinging on nonuniformly heated fluid interfaces and solid surfaces. Bulletin of the American Physical Society , 2023.                                                                                        |
| [18]   | Bashir Mohamed Elbousefi. Lattice Boltzmann Simulations of Interfacial Dynamics in Self-Rewetting Fluid Layers, Bubbles, and Drops . PhD thesis, University of Colorado Denver, 2024.                                                                                                                                          |
| [19]   | Yoshiyuki Abe. Terrestrial and microgravity applications of self-rewetting fluids. Mi- crogravity Science and Technology , 19(3):11-12, 2007.                                                                                                                                                                                  |
| [20]   | Raffaele Savino, Anselmo Cecere, and Roberto Di Paola. Surface tension-driven flow in wickless heat pipes with self-rewetting fluids. International Journal of Heat and Fluid Flow , 30(2):380-388, 2009.                                                                                                                      |
| [21]   | Raffaele Savino, Anselmo Cecere, Stefan Van Vaerenbergh, Yoshiyuki Abe, G Pizzirusso, WTzevelecos, Mohamed Mojahed, and Quentin Galand. Some experi- mental progresses in the study of self-rewetting fluids for the SELENE experiment to be carried in the thermal platform 1 hardware. Acta Astronautica , 89:179-188, 2013. |
| [22]   | Yanxin Hu, Tengqing Liu, Xuanyou Li, and Shuangfeng Wang. Heat transfer enhance- ment of micro oscillating heat pipes with self-rewetting fluid. International Journal of Heat and Mass Transfer , 70:496-503, 2014.                                                                                                           |
| [23]   | Shen-Chun Wu, Tien-Ju Lee, Wei-Jhih Lin, and Yau-Ming Chen. Study of self- rewetting fluid applied to loop heat pipe with ptfe wick. Applied Thermal Engineering , 119:622-628, 2017.                                                                                                                                          |
| [24]   | Anselmo Cecere, Davide De Cristofaro, Raffaele Savino, Vincent Ayel, Thibaud Sole- Agostinelli, Marco Marengo, Cyril Romestant, and Yves Bertin. Experimental anal- ysis of a flat plate pulsating heat pipe with self-rewetting fluids during a parabolic flight campaign. Acta Astronautica , 147:454-461, 2018.             |
| [25]   | Yanxin Hu, Kaixin Huang, and Jin Huang. A review of boiling heat transfer and heat pipes behaviour with self-rewetting fluids. International Journal of Heat and Mass Transfer , 121:107-118, 2018.                                                                                                                            |
| [26]   | Minghan Zhu, Jin Huang, Mengjie Song, and Yanxin Hu. Thermal performance of a thin flat heat pipe with grooved porous structure. Applied Thermal Engineering , 173:115215, 2020.                                                                                                                                               |

| [27]   | Dongwei Zhang, Zhuantao He, Jian Guan, Songzhen Tang, and Chao Shen. Heat transfer and flow visualization of pulsating heat pipe with silica nanofluid: An exper- imental study. International Journal of Heat and Mass Transfer , 183:122100, 2022.                                  |
|--------|---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| [28]   | Jhih-Huang Gao, Shen-Chun Wu, Ya-Wei Lee, Ta-Li Chou, and Yan-Chun Chen. The study of novel self-rewetting fluid application to loop heat pipe. Applied Sciences , 12(6):3121, 2022.                                                                                                  |
| [29]   | Zipei Su, Yanxin Hu, Shaobin Zheng, Chao Dang, Kaizhao Liu, Jiapei Chen, and Jin Huang. Study on heat transfer characteristics of annular pulsating heat pipe with temperature oscillation characteristic. International Communications in Heat and Mass Transfer , 144:106779, 2023. |
| [30]   | Caihang Liang, Xiong Yao, Miaowei Liu, Juan Du, Yuxiang Hong, and Lei Zhao. Comparisons of thermal performances in a pulsating heat pipe by using a nanofluid and a self-rewetting nanofluid with carbon nanotubes. Case Studies in Thermal En- gineering , 41:102589, 2023.          |
| [31]   | Anze Sitar and Iztok Golobic. Heat transfer enhancement of self-rewetting aqueous n-butanol solutions boiling in microchannels. International Journal of Heat and Mass Transfer , 81:198-206, 2015.                                                                                   |
| [32]   | K Sefiane, X Yu, Gail Duursma, and J Xu. On heat and mass transfer using evaporating self-rewetting mixtures in microchannels. Applied Thermal Engineering , 179:115662, 2020.                                                                                                        |
| [33]   | Yanxin Hu, Suling Zhang, Xuanyou Li, and Shuangfeng Wang. Heat transfer en- hancement of subcooled pool boiling with self-rewetting fluid. International Journal of Heat and Mass Transfer , 83:64-68, 2015.                                                                          |
| [34]   | Yanxin Hu, Sixu Chen, Jin Huang, and Mengjie Song. Marangoni effect on pool boiling heat transfer enhancement of self-rewetting fluid. International Journal of Heat and Mass Transfer , 127:1263-1270, 2018.                                                                         |
| [35]   | Yanxin Hu, Hai Wang, Mengjie Song, and Jin Huang. Marangoni effect on microbub- bles emission boiling generation during pool boiling of self-rewetting fluid. Interna- tional Journal of Heat and Mass Transfer , 134:10-16, 2019.                                                    |
| [36]   | Matic Moˇe, z Viktor Vajc, Matevˇ z Zupanˇiˇ, c c Radek ˇ Sulc, and Iztok Golobiˇ. c Pool boil- ing performance of water and self-rewetting fluids on hybrid functionalized aluminum surfaces. Processes , 9(6):1058, 2021.                                                           |
| [37]   | Yong Il Kim, Boo-Hyoung Bang, Keunhee Jang, Seongpil An, Alexander L Yarin, and Sam S Yoon. Pool boiling enhancement via nanotexturing and self-propelled swing motion for bubble shedding. International Communications in Heat and Mass Transfer , 133:105934, 2022.                |
| [38]   | Ibrahim Zaaroura, Souad Harmand, Julien Carlier, Malika Toubal, Aur´lie e Fasquelle,                                                                                                                                                                                                  |

- [38] Ibrahim Zaaroura, Souad Harmand, Julien Carlier, Malika Toubal, Aur´lie Fasquelle, e and Bertrand Nongaillard. Thermal performance of self-rewetting gold nanofluids:

|      | Application to two-phase heat transfer devices. International Journal of Heat and Mass Transfer , 174:121322, 2021.                                                                                                                          |
|------|----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| [39] | Martin ER Shanahan and Khellil Sefiane. Recalcitrant bubbles. Scientific Reports , 4(1):1-9, 2014.                                                                                                                                           |
| [40] | Dimitrios Mamalis, Vasileios Koutsos, and Khellil Sefiane. Bubble rise in a non- isothermal self-rewetting fluid and the role of thermocapillarity. International Journal of Thermal Sciences , 117:146-162, 2017.                           |
| [41] | Milton J Rosen and Joy T Kunjappu. Surfactants and interfacial phenomena . John Wiley & Sons, 2012.                                                                                                                                          |
| [42] | Charles N Baroud, Francois Gallaire, and R´mi e Dangla. Dynamics of microfluidic droplets. Lab on a Chip , 10(16):2032-2045, 2010.                                                                                                           |
| [43] | Ralf Seemann, Martin Brinkmann, Thomas Pfohl, and Stephan Herminghaus. Droplet based microfluidics. Reports on progress in physics , 75(1):016601, 2011.                                                                                     |
| [44] | Nina M Kovalchuk, Marten Reichow, Thomas Frommweiler, Daniele Vigolo, and Mark JH Simmons. Mass transfer accompanying coalescence of surfactant-laden and surfactant-free drop in a microfluidic channel. Langmuir , 35(28):9184-9193, 2019. |
| [45] | Xiao Luo, Zheng Yuan Luo, and Bo Feng Bai. Effect of thermal convection on ther- mocapillary migration of a surfactant-laden droplet in a microchannel. Physics of Fluids , 32(9):092009, 2020.                                              |
| [46] | Nina MKovalchuk and Mark JH Simmons. Review of the role of surfactant dynamics in drop microfluidics. Advances in Colloid and Interface Science , page 102844, 2023.                                                                         |
| [47] | Shervin Issakhani, Omid Jadidi, Jafar Farhadi, and Vahid Bazargan. Geometrically- controlled evaporation-driven deposition of conductive carbon nanotube patterns on inclined surfaces. Soft Matter , 19(7):1393-1406, 2023.                 |
| [48] | Mohamed Qenawy, Yiqi Chen, Junfeng Wang, Jiameng Tian, and Bin Chen. Inter- mittent cryogen spray cooling coupled with cold air jet for heat transfer enhancement and cryogen saving of laser dermatology. Physics of Fluids , 36(2), 2024.  |
| [49] | Masoumeh Mirhoseini, Alireza Banaee, and Alireza Jalali. A novel phase-field lattice Boltzmann framework for diffusion-driven multiphase evaporation. Physics of Fluids , 36(8), 2024.                                                       |
| [50] | He Yan, Lei Wang, Jiangxu Huang, and Yuan Yu. Thermocapillary migration of a self-rewetting droplet on an inclined surface: A phase-field simulation. arXiv preprint arXiv:2407.12430 , 2024.                                                |
| [51] | SG Slavtchev and SP Miladinova. Thermocapillary flow in a liquid layer at minimum in surface tension. Acta Mechanica , 127(1):209-224, 1998.                                                                                                 |
| [52] | Alexander Oron and Philip Rosenau. On a nonlinear thermocapillary effect in thin liquid layers. Journal of Fluid Mechanics , 273:361-374, 1994.                                                                                              |

| [53]   | William Batson, Yehuda Agnon, and Alex Oron. Thermocapillary modulation of self-rewetting films. Journal of Fluid Mechanics , 819:562-591, 2017.                                                                                                                                   |
|--------|------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| [54]   | Z Yu. Thermocapillary instability of self-rewetting films on vertical fibers. Physics of Fluids , 30(8):082104, 2018.                                                                                                                                                              |
| [55]   | Mohammed Zubair and Rajagopal Vellingiri. Dynamics of thin self-rewetting liquid films on an inclined heated substrate. Physics of Fluids , 35(11), 2023.                                                                                                                          |
| [56]   | Manoj Kumar Tripathi, KC Sahu, G Karapetsas, K Sefiane, and OK Matar. Non- isothermal bubble rise: non-monotonic dependence of surface tension on temperature. Journal of Fluid Mechanics , 763:82-108, 2015.                                                                      |
| [57]   | Mounika Balla, Manoj Kumar Tripathi, Kirti Chandra Sahu, George Karapetsas, and Omar K Matar. Non-isothermal bubble rise dynamics in a self-rewetting fluid: three-dimensional effects. Journal of Fluid Mechanics , 858:689-713, 2019.                                            |
| [58]   | Mohammad Majidi, Reza Haghani-Hassan-Abadi, and Mohammad-Hassan Rahimian. Single recalcitrant bubble simulation using a hybrid lattice Boltzmann finite difference model. International Journal of Multiphase Flow , 127:103289, 2020.                                             |
| [59]   | TR Mitchell, MMajidi, MHRahimian, and CR Leonardi. Computational modeling of three-dimensional thermocapillary flow of recalcitrant bubbles using a coupled lattice Boltzmann-finite difference method. Physics of Fluids , 33(3):032108, 2021.                                    |
| [60]   | Ze-Lai Xu, Jun-Yuan Chen, Hao-Ran Liu, Kirti Chandra Sahu, and Hang Ding. Motion of self-rewetting drop on a substrate with a constant temperature gradient. Journal of Fluid Mechanics , 915:A116, 2021.                                                                          |
| [61]   | Bashir Elbousefi, William Schupbach, Kannan N Premnath, and Samuel WJ Welch. Thermocapillary convection in superimposed layers of self-rewetting fluids: Analytical and lattice Boltzmann computational study. International Journal of Heat and Mass Transfer , 208:124049, 2023. |
| [62]   | Haozhe Su, Chunxi Li, Dan Li, and Xuemin Ye. Enhanced spreading of surfactant- containing, self-rewetting fluids in pulmonary drug delivery. Physics of Fluids , 34(11), 2022.                                                                                                     |
| [63]   | V Sharanya, GP Sekhar, and Christian Rohde. Surfactant-induced migration of a spherical droplet in non-isothermal Stokes flow. Physics of Fluids , 31(1), 2019.                                                                                                                    |
| [64]   | Zhenlin Guo. Nonlinear coupling effects of the thermocapillarity and insoluble surfac- tants to droplet migration under poiseuille flow. Physical Review Fluids , 8(2):024001, 2023.                                                                                               |
| [65]   | Roberto Benzi, Sauro Succi, and Massimo Vergassola. The lattice Boltzmann equa- tion: theory and applications. Physics Reports , 222(3):145-197, 1992.                                                                                                                             |
| [66]   | Dazhi Yu, Renwei Mei, Li-Shi Luo, and Wei Shyy. Viscous flow computations with the method of lattice Boltzmann equation. Progress in Aerospace sciences , 39(5):329-367,                                                                                                           |

| [67]   | Pierre Lallemand, Li-shi Luo, Manfred Krafczyk, and Wen-An Yong. The lattice Boltzmann method for nearly incompressible flows. Journal of Computational Physics , 431:109713, 2021.                                                                                                                       |
|--------|-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| [68]   | Keerti Vardhan Sharma, Robert Straka, and Frederico Wanderley Tavares. Current status of lattice Boltzmann methods applied to aerodynamic, aeroacoustic, and ther- mal flows. Progress in Aerospace Sciences , 115:100616, 2020.                                                                          |
| [69]   | KJ Petersen and JR Brinkerhoff. On the lattice boltzmann method and its application to turbulent, multiphase flows of various fluids including cryogens: A review. Physics of Fluids , 33(4), 2021.                                                                                                       |
| [70]   | Runa Samanta, Himadri Chattopadhyay, and Chandan Guha. A review on the appli- cation of lattice Boltzmann method for melting and solidification problems. Compu- tational Materials Science , 206:111288, 2022.                                                                                           |
| [71]   | Seyed Ali Hosseini and Iliya V Karlin. Lattice Boltzmann for non-ideal fluids: Fun- damentals and practice. Physics Reports , 1030:1-137, 2023.                                                                                                                                                           |
| [72]   | Seyed Ali Hosseini, Pierre Boivin, Dominique Th´venin, e and Ilya Karlin. Lattice Boltzmann methods for combustion applications. Progress in Energy and Combustion Science , 102:101140, 2024.                                                                                                            |
| [73]   | Xiaoyi He, Shiyi Chen, and Raoyang Zhang. A lattice Boltzmann scheme for in- compressible multiphase flow and its application in simulation of Rayleigh-Taylor instability. Journal of computational physics , 152(2):642-663, 1999.                                                                      |
| [74]   | Xiaoyi He and Gary D Doolen. Thermodynamic foundations of kinetic theory and lattice Boltzmann models for multiphase flows. Journal of Statistical Physics , 107(1):309-328, 2002.                                                                                                                        |
| [75]   | Taehun Lee and Ching-Long Lin. A stable discretization of the lattice Boltzmann equation for simulation of incompressible two-phase flows at high density ratio. Jour- nal of Computational Physics , 206(1):16-47, 2005.                                                                                 |
| [76]   | Kannan N Premnath and John Abraham. Three-dimensional multi-relaxation time (MRT) lattice-Boltzmann models for multiphase flow. Journal of Computational Physics , 224(2):539-559, 2007.                                                                                                                  |
| [77]   | Farzaneh Hajabdollahi, Kannan NPremnath, and Samuel WJWelch. Central moment lattice Boltzmann method using a pressure-based formulation for multiphase flows at high density ratios and including effects of surface tension and Marangoni stresses. Journal of Computational Physics , 425:109893, 2021. |
| [78]   | Tahani Alsadik, Samuel WJ Welch, and Kannan N Premnath. Lattice Boltzmann simulations of quasi-steady film and axisymmetric nucleate boiling. Physics of Fluids ,                                                                                                                                         |

| [79]   | Haihu Liu, Yonghao Zhang, and Albert J Valocchi. Modeling and simulation of thermocapillary flows using lattice Boltzmann method. Journal of Computational Physics , 231(12):4433-4453, 2012.                                |
|--------|------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| [80]   | Robert Scherr, Matthias Markl, and Carolin K¨rner. o Volume of fluid based modeling of thermocapillary flow applied to a free surface lattice Boltzmann method. Journal of Computational Physics , 492:112441, 2023.         |
| [81]   | Martin Geier, Andreas Greiner, and Jan G Korvink. Cascaded digital lattice Boltz- mann automata for high Reynolds number flow. Physical Review E , 73(6):066705, 2006.                                                       |
| [82]   | Timm Kr¨ger, u Halim Kusumaatmaja, Alexandr Kuzmin, Orest Shardt, Goncalo Silva, and Erlend Magnus Viggen. The lattice Boltzmann method. Springer International Publishing , 10(978-3):4-15, 2017.                           |
| [83]   | RGM Van der Sman and S Van der Graaf. Diffuse interface model of surfactant adsorption onto flat and droplet interfaces. Rheologica acta , 46(1):3-11, 2006.                                                                 |
| [84]   | Haihu Liu and Yonghao Zhang. Phase-field modeling droplet dynamics with soluble surfactants. Journal of Computational Physics , 229(24):9166-9187, 2010.                                                                     |
| [85]   | RGM Van Der Sman and MBJ Meinders. Mesoscale models of dispersions stabilized by surfactants and colloids. Advances in colloid and interface science , 211:63-76, 2014.                                                      |
| [86]   | Ruud GM van der Sman and MBJ Meinders. Analysis of improved lattice Boltz- mann phase field method for soluble surfactants. Computer Physics Communications , 199:12-21, 2016.                                               |
| [87]   | Y Shi, GH Tang, LH Cheng, and HQ Shuang. An improved phase-field-based lattice Boltzmann model for droplet dynamics with soluble surfactant. Computers & Fluids , 179:508-520, 2019.                                         |
| [88]   | Stefan Engblom, Minh Do-Quang, Gustav Amberg, and Anna-Karin Tornberg. On diffuse interface modeling and simulation of surfactants in two-phase fluid flow. Com- munications in Computational Physics , 14(4):879-915, 2013. |
| [89]   | Giovanni Soligo, Alessio Roccon, and Alfredo Soldati. Coalescence of surfactant-laden drops by phase field method. Journal of Computational Physics , 376:1292-1311, 2019.                                                   |
| [90]   | Farzaneh Hajabdollahi. Cascaded Lattice Boltzmann Methods Based on Central Mo- ments for Thermal Convection, Multiphase Flows and Complex Fluids . PhD thesis, University of Colorado Denver, 2019.                          |
| [91]   | Pao-Hsiung Chiu and Yan-Ting Lin. A conservative phase field method for solving incompressible two-phase flows. Journal of Computational Physics , 230(1):185-204, 2011.                                                     |

| [92]   | Ying Sun and Christoph Beckermann. Sharp interface tracking using the phase-field equation. Journal of Computational Physics , 220(2):626-653, 2007.                                                                            |
|--------|---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| [93]   | R Folch, J Casademunt, A Hern´ndez-Machado, a and L Ramirez-Piscina. Phase-field model for hele-shaw flows with arbitrary viscosity contrast. i. theoretical approach. Physical Review E , 60(2):1724, 1999.                    |
| [94]   | Jeremiah U Brackbill, Douglas B Kothe, and Charles Zemach. A continuum method for modeling surface tension. Journal of computational physics , 100(2):335-354, 1992.                                                            |
| [95]   | Anand Kumar. Isotropic finite-differences. Journal of Computational Physics , 201(1):109-118, 2004.                                                                                                                             |
| [96]   | Hang Ding, Peter DM Spelt, and Chang Shu. Diffuse interface model for incompress- ible two-phase flows with large density ratios. Journal of Computational Physics , 226(2):2078-2095, 2007.                                    |
| [97]   | Kannan N Premnath and Sanjoy Banerjee. Incorporating forcing terms in cascaded lattice Boltzmann approach by method of central moments. Physical Review E , 80(3):036702, 2009.                                                 |
| [98]   | Kannan N Premnath and Sanjoy Banerjee. On the three-dimensional central moment lattice Boltzmann method. Journal of Statistical Physics , 143:747-794, 2011.                                                                    |
| [99]   | Eman Yahia and Kannan N Premnath. Central moment lattice Boltzmann method on a rectangular lattice. Physics of Fluids , 33(5):057110, 2021.                                                                                     |
| [100]  | Eman Yahia, William Schupbach, and Kannan N Premnath. Three-dimensional cen- tral moment lattice Boltzmann method on a cuboid lattice for anisotropic and inho- mogeneous flows. Fluids , 6(9):326, 2021.                       |
| [101]  | Roland Clift, John R Grace, and Martin E Weber. Bubbles, drops, and particles. 2005.                                                                                                                                            |
| [102]  | S Aland and A Voigt. Benchmark computations of diffuse interface models for two- dimensional bubble dynamics. International Journal for Numerical Methods in Fluids , 69(3):747-761, 2012.                                      |
| [103]  | Y Wang, Chang Shu, JY Shao, J Wu, and XD Niu. A mass-conserved diffuse interface method and its application for incompressible multiphase flows with large density ratio. Journal of Computational Physics , 290:336-351, 2015. |
| [104]  | Z Chen, C Shu, D Tan, XD Niu, and QZ3923148 Li. Simplified multiphase lattice Boltzmann method for simulating multiphase flows with large density ratios and com- plex interfaces. Physical Review E , 98(6):063314, 2018.      |
| [105]  | Kannan N Premnath, Farzaneh Hajabdollahi, and Samuel WJ Welch. Surfactant effects on interfacial flow and thermal transport processes during phase change in film boiling. Physics of Fluids , 30(4), 2018.                     |

[106] Martin Geier, Martin Sch¨nherr, Andrea Pasquali, and Manfred Krafczyk. o The cumulant lattice Boltzmann equation in three dimensions: Theory and validation. Computers &amp; Mathematics with Applications , 70(4):507-547, 2015.
Bashir Elbousefi William Schupbach Kannan N. Premnath Samuel W. J. Welch
Investigation of Surfactant-Laden Bubble Migration Dynamics in Self-Rewetting Fluids using Lattice Boltzmann Method
Self-rewetting fluids (SRFs) (e.g., aqueous solutions of long-chain alcohols)
show anomalous nonlinear (quadratic) variations of surface tension with
temperature involving a positive gradient, leading to different thermocapillary
convection compared to normal fluids (NFs). Moreover, surface-active materials
or surfactants can significantly alter interfacial dynamics by their adsorption
on fluid interfaces. The coupled effects of temperature- and surfactant-induced
Marangoni stresses on migration bubbles in SRFs remain unexplored. We use a
robust lattice Boltzmann (LB) method based on central moments to simulate the
two-fluid motions, capture interfaces, and compute the transport of energy and
surfactant concentration fields, and systematically study the surfactant-laden
bubble dynamics in SRFs. When compared to motion of bubbles in normal fluids,
in which they continuously migrate without a stationary behavior, our results
show that they exhibit dramatically different characteristics in SRFs. Not only
is the bubble motion directed towards the minimum temperature location in SRFs,
but, more importantly, the bubble attains an equilibrium position. In the
absence of surfactants, such an equilibrium position arises at the minimum
reference temperature occurring at the center of the domain. The addition of
surfactants moves the equilibrium location further upstream, which is
controlled by the magnitude of the Gibbs elasticity parameter that determines
the magnitude of the surface tension variation with surfactants. The quadratic
sensitivity coefficient of surface tension on temperature associated with the
SRF modulates this behavior. The lower this quantity, the greater is the role
of surfactants modifying the equilibrium position of the bubble in SRF. These
findings provide new means to potentially manipulate the bubble dynamics, and
especially to tune its equilibrium states.