**1. Introduction**

Foams exist in several chemical and industrial processes, and they often delay and obstruct the production. Foam is a complex multiphase multicomponent system, and it often consists of one or more coexisting gas components (air and vapor) in the bubble and liquid phase represented in the lamella (water). The lamellar interfacial length can reach orders of nanometers depending on whether the stabilizing agents are surfactants or proteins [1]. The lamella of protein foams normally have higher length scales than surfactants' ones [2]. The repulsive forces between the surfactants or protein molecules in the lamellar interface often induce high positive values of disjoining pressure, which is the reason for stable foams [3]. Quantification of these disjoining pressure isotherms is often undertaken by thin film pressure balance technique, where the thickness of the thin film is controlled and measured with the film pressure until reaching the Black Newton film and, hence, rupture [4].

Thermal effects and gradients add more physical complexities to the foam system, such as evaporation–condensation, viscosity variations and the local variation of surface tension, leading to the Marangoni effect [5]. These complex phenomena interact with each other. For example, one of the most foam-dynamics influencing phenomena is the liquid drainage through the lamella [2]. However, temperature variation can highly affect the surface tension, leading to an induced flow, which may be in an unfavorable direction opposing the drainage one, hence leading to more stable foam [6]. There also exists the coarsening or Ostwald ripening in foams, where gas diffuses from smaller bubbles to

**Citation:** Mobarak, M.; Gatternig, B.; Delgado, A. On the Simulations of Thermal Liquid Foams Using Lattice Boltzmann Method. *Energies* **2023**, *16*, 195. https://doi.org/10.3390/ en16010195

Academic Editors: Artur Bartosik and Dariusz Asendrych

Received: 14 November 2022 Revised: 19 December 2022 Accepted: 20 December 2022 Published: 24 December 2022

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

larger bubbles [1]. This process is normally more important on the larger time scales, when compared to the drainage and convective mass flow.

Modelling and simulation of these interacting processes are not straightforward, although some aspects can be tackled using the mesoscopic lattice Boltzmann method (LBM), which has proven its capability of modelling complex multiphase and thermal flows [7]. The idea of the current work is to propose an approach to simulate some of these phenomena and discuss the current limitations based on the literature and current experience.

For the multiphase models in LBM, the two most famous methods are the freeenergy [8,9] and the pseudopotential Shan-Chen (SC) [10,11] approaches. The limitation of the high density ratio and viscosity ratios are always issues for both methods. The adopted approach in this work is the SC model, since it is possible to extend it to approach the simulation of liquid foams by reaching positive values of disjoining pressure. The original SC model can also reach higher stable density ratios up to 100 and can be extended to reach ratios up to 1000, though with a modified equations of state (EOS). In the work by Yuan and Schaefer [12], they studied several equations of state, where Peng-Robinson (P-R) EOS [13] yielded the highest stable density ratio and the least spurious currents, which is a well-known issue in multiphase models. The higher viscosity ratios issue can generally be stabilized by switching from the standard LBM single relaxation time (SRT) collision operator to either the two or multiple relaxation time (TRT or MRT) [14]. SC model can also reach comparable and even slightly higher surface tensions than the free-energy approach [14].

There also exists the free-surface (FS) LBM approach, where the dense component is considered while the other one is neglected; this somehow resembles the Volume of Fluid (VOF) approach in LBM. In the work by Körner et al. [15] and later by Ataei et al. [16], much success has been gained with reaching stable foam by introducing linearly variable repulsive forces as a function of the lamellar distance. Interfacial curvature estimation [17] is required for the FS-LBM in general, in addition to ray tracing algorithms for the lamellar distance evaluation. Nevertheless, the presented work is mainly focused on the SC LBM approach, where neither the tracking nor the curvature estimation is required, not to mention the lamellar distance evaluation.

The SC LBM is based on modelling the molecular interaction between phases and components to reach the phase separation "bottom-up approach" [18]. The original SC can be extended to model both multiphase and multicomponent system, while being able to tune the density ratio and surface tension separately [19]. In order to simulate foams, the surfactants effects shall be implemented into these interaction forces. This can be carried out by using the mid-range interaction (repulsion) in addition to the short-range one [20,21]. This benefit emerges from the nature of the SC LBM, modeling the molecular interactions themselves. This can lead to a positive disjoining pressure, hence delaying bubbles coalescence and obtaining foam! It is worth mentioning that even using a mesoscale method such as LBM, it is practically impossible to resolve the length scales of surfactants molecules; that is why only their forces' effects are modelled.

In the literature, there also exists the model presented by Chen et al. [22] and Nekovee et al. [23], which is based on introducing two extra sets of LBM distribution functions in addition to the SC LBM system, representing the surfactants species (position) and their own dipolar orientation, although, in the current work, we preferred avoiding this approach for foam simulations for three reasons. First, the model adds much more complexity to the already saturated system of LBM distribution functions. Second, relating the relaxation time, hence the diffusivity of the surfactants position and orientation, to the physical ones could not be seen to be possible. On the other hand, in the mid-range interaction model, the resulting disjoining pressure can be evaluated and, hence, there might be a way to relate it to the experimental isotherms. Third, it is reported in the work from Mukherjee et al. [24] that the model by itself could not inhibit the bubbles coalescence.

Thermal phase change from liquid to gas and vice versa using LBM has been proposed by Gong and Cheng [25] and further enhanced by Li et al. [26]. The model is based on coupling the temperature field by including a source term responsible for the enthalpy of phase change with the SC multiphase model, which is modified by the temperature dependent P-R EOS. However, in addition to the stability issues which are encountered, maintaining both the thermodynamic consistency for phase change and obtaining positive disjoining pressure, which is an essential property of foam formation and stability, has to our knowledge not been developed in the literature.

In the literature, the Enthalpy-based LBM to model phase change for liquid to solid states also exists [27,28], which might be a compromise to investigate the phase change in extremely stable foam by only tracking the interface front between liquid and gas phases accordingly with the temperature field.

The current work aims to propose a recipe for model combination in the scope of SC LBM, in order to simulate the thermal effects in liquid foam without taking into account the phases change. The convective diffusive heat transfer and thermal effects consideration while simulating the flow field in liquid foam is the key innovation in this work. Coarsening or Ostwald ripening will not be discussed in this work, since controlling the mass diffusivity of gas in liquids using the SC model separately is numerically impossible.

The paper is arranged as follows: the LBM and SC model are briefly proposed, showing the multiphase–multicomponent simulation parameters, the implemented forcing scheme and the TRT. The mid-range interaction model is proposed and explained. The temperature field and how it is implemented is shown. Units and non-dimensional groups to change from lattice to physical units are discussed. For the following section, a model verification is carried out by relating the interaction numerical parameters to the physical ones through the Young–Laplace test, where surface tension and density ratio (using SC EOS) are tuned and evaluated for different bubble sizes. Rayleigh–Taylor instability (RTI) is shown as a validation case for the implemented multiphase–multicomponent model and compared with numerical test cases from the literature. Simulation samples from the bubble rise diagram by Clift et al. [29] are presented in brief as a further qualitative validation for different Reynolds and bonds numbers, showing the model capability of capturing the interfacial behavior for different bubble regimes. In the application section, the developed model is applied for a case study of one-way coupled heat transfer in dynamic liquid foam while using the mid-range interaction model. Finally, a new technique to consider the variable surface tension due to temperature gradients is proposed, hence suggesting a technique to account for the Marangoni effect in thermal foams using SC LBM. A summary and conclusion section is shown to finalize the whole presented model, discuss the achievement fairly and suggest future work in order to reach the final aim of simulating the foam in industrial scale processes.
