A multiphase model for exploring electrochemical Marangoni flow

Electrochemistry Communications 155 (2023) 107567 Available online 1 September 2023 1388-2481/© 2023 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). A multiphase model for exploring electrochemical Marangoni flow E. Karimi-Sibaki a , A. Vakhrushev a , A. Kadylnykova a , M. Wu b , A. Ludwig b , J. Bohacek c , A. Kharicha a , b , * a Christian-Doppler Laboratory for Metallurgical Applications of Magnetohydrodynamics, Montanuniversitaet of Leoben, Franz-Josef-Str. 18, A-8700 Leoben, Austria b Chair of Simulation and Modeling of Metallurgical Processes, Montanuniversitaet of Leoben, Franz-Josef-Str. 18, A-8700 Leoben, Austria c Heat Transfer and Fluid Flow Laboratory, Faculty of Mechanical Engineering, Brno University of Technology, Technicka 2896/2, 616 69 Brno, Czech Republic ARTICLE INFO Keywords: Volume of fluid (VOF) Electro-Marangoni flow Eutectic gallium indium alloy (EGaIn) Droplet elongation Interfacial tension modulation ABSTRACT A multiphase numerical model based on the volume of fluid (VOF) method is proposed to simulate the transient, electrochemically-generated Marangoni flow in a system comprising a NaOH electrolyte and a eutectic gallium–indium (EGaIn) metal droplet. The model incorporates appropriate equations to accurately represent the transport phenomena, including flow, electric potential, and electric current density, within the entire system. The model includes the transient variation in the interfacial tension as a function of electric current density at the interface, leading to the generation of Marangoni flow and enabling the tracking of droplet shape evolution. Notably, the model successfully captures the elongation of the droplet towards the cathode, which is validated through comparison with available experimental data. 1. Introduction Within the realms of fluid mechanics, magnetohydrodynamics (MHD), and electrochemistry, the understanding of fluid instabilities and interfacial dynamics in electrochemical systems involving liquid metals and electrolytes, subject to the influence of electric and magnetic fields, is of utmost importance. These phenomena have significant implications in various processes such as drop casting [1], liquid metal batteries [2], electroslag remelting [3], aluminium electrolysis [4], welding [5], and, in general, electrolytic processes involving hydrogen bubbles [6]. Active control of interfacial dynamics with diverse applications in prosthetics, soft robotics, tunable electronics, and electrochemical sensors has prompted extensive interest in studying and harnessing the unique characteristics of the eutectic gallium–indium alloy (EGaIn), which possesses high surface tension, high electrical and thermal conductivities, low toxicity, and a melting point below room temperature [7–10]. Specifically, the substantial modulation of interfacial tension in EGaIn through the application of electric fields provides a fertile ground for studying and exploring a wide range of fascinating physical phenomena, including electrohydrodynamic [11] and Marangoni effects [12–14]. In this study, we focus on the phenomenon of Marangoni flow, which emerges due to the gradient in the interfacial tension resulting from the electrochemical oxidation of gallium at the surface of an EGaIn droplet [15]. This gradient in interfacial tension exerts stress on the metal- –electrolyte interface, thereby driving the movement of the electrolyte (e.g., NaOH) from regions of lower tension to regions of higher tension [11,16,17]. Herein, we propose a novel multiphase framework based on the volume of fluid (VOF) method [18] to accurately simulate the transport phenomena in the EGaIn droplet and NaOH electrolyte system. The framework utilizes the finite-volume method (FVM) for discretization [19], ensuring the conservation of mass, momentum, and current density. The VOF-based model effectively represents a sharp interface, enabling the capture of the complex shape evolution of the droplet across different process sizes and geometries. By investigating the effect of applied potential on the system, our aim is to gain fundamental insights into the electrochemically-generated Marangoni flow. These insights can be leveraged by industry professionals to optimize their practices and improve operational efficiency. 2. Modeling All symbols are defined in Appendix A. The electrolyte and the metal are considered as two incompressible and immiscible liquid phases. All governing equations are listed in Table 1. The volume of vluid (VOF) method [18] is used to recognize the two * Corresponding author. E-mail address: [email protected] (A. Kharicha). Contents lists available at ScienceDirect Electrochemistry Communications journal homepage: www.elsevier.com/locate/elecom https://doi.org/10.1016/j.elecom.2023.107567 Received 5 July 2023; Received in revised form 28 August 2023; Accepted 31 August 2023 Electrochemistry Communications 155 (2023) 107567 2 phases by incorporating a marker function (β) that depends on spatial and temporal factors, as stated in Eq. (1). To model the interface transport, an advection equation is solved for β on a fixed Eulerian grid, as described in Eq. (2). In the context of the continuum surface force (CSF) model [20], the unit normal vector ( n →) at the interface is computed using Eq. (3), while the interface curvature is determined using Eq. (4). The flow field is determined by applying continuity, Eq. (5), and momentum, Eq. (6), conservation principles. The density, ρ =β ρ m+ (1−β) ρ e, and viscosity, μ =β μ m+ (1−β) μ e, of the mixture are calculated as weighted averages of the respective properties of the two phases within each computational cell. The force of surface tension, Eq. (7), depends on the interfacial tension [20]. δ is one at the interface and zero elsewhere. Herein, we investigate the electrochemical oxidation of eutectic gallium–indium (EGaIn), which results in a modification of its surface tension. The anodic droplet undergoes faradaic reactions, leading to the generation of Ga 3+ ions [15]. The electrochemical oxidation of gallium creates a gradient of electric potential and interfacial tension across the droplet’s surface that induces Marangoni flow, where the electrolyte flows from low to high interfacial tension regions [15]. Despite the electrochemical oxidation process, the alteration in interfacial tension exhibits complete reversibility, accompanied by minimal hysteresis effects [16]. Therefore, it becomes feasible to characterize the interfacial tension as a function of electric potential or electric current density [15–17], with the objective of modeling the Marangoni flow induced by electrochemical reactions. The study conducted by Hillaire et al. [17] presents the observed changes in interfacial tension as a function of electric current density at the metal–electrolyte interface, as shown in Fig. 1(a). By employing a curve-fitting approach, we effectively model this relationship, which is expressed by Eq. (8). Considering primary current distribution to determine the electric potential field of the electrolyte [21,22], the Laplace equation, Eq. (9), must be solved. The electrical conductivity of the metal surpasses that of the electrolyte by several orders of magnitude. As a result, to accurately model the minimal electric field within the metal, a sink term for electric potential (Sφ= − β Kφ) is implemented solely within computational cells fully occupied by the metal [23,24]. This modification transforms the original Laplace equation into an elliptic Helmholtz-type equation. By using a small value for the permeability parameter (e.g., K=10−10), the electric potential within the metal effectively diminishes to zero [23]. It should be noted that employing an excessively small value (e.g., 10−20) could result in numerical instability and convergence issues. The electric current density is subsequently calculated utilizing Ohm’s law, as described in Eq. (10). The models were configured based on the study conducted by Li et al. [25]. Considering a 2D symmetrical model in this study, the computational domain and boundaries are schematically shown in Fig. 1 (b). The boundary conditions for the flow encompass symmetrical conditions and no-slip conditions at the anode, cathode, and insulating wall. Zero flux of potential is enforced at the insulating and symmetry walls. At the anode, the electric potential is set to zero, while at the cathode, the potential is either −0.5 V or −5 V, depending on the operating parameters of the system. The finite-volume method (FVM) was employed to discretize the governing equations listed in Table 1. The FVM employs a first-order implicit temporal discretization scheme along with a third-order MUSCL spatial discretization scheme [19]. The Geo-Reconstruct-based explicit scheme was used to discretize the volume fraction equation [26]. A set of user-defined functions (UDFs) was implemented into the ANSYS FLUENT v.17.1 commercial software. These UDFs were meticulously designed to account for critical aspects related to surface tension force and electric potential. Transient simulations were performed to dynamically capture and monitor the temporal deformation of the droplet. The computational domain was filled using a total of 5500 quadrilateral mesh elements. To achieve a converged solution, the CFL number, (uΔt a), is kept below 1 using a time-step size of ca. 10 -3 s. The computation time on a single processor Intel Core i7-4790 K CPU 4.00 Table 1 All governing equations. Volume fraction β(x →,t) = Vm V=⎧ ⎨ ⎩ 1,x →∈Metal 0<β<1,x →∈Interface 0,x →∈Electrolyte (1) ∂ β ∂ t+ ∇ • (βu) = 0 (2) n →= (nx,ny) = ∇β ‖∇β‖ (3) κ= ∇ • ( n →)(4) Flow field ∂ρ ∂ t+ ∇ • ( ρ u) = 0 (5) ∂ρ u ∂ t+ ∇ • ( ρ uu) = − ∇p+ ∇ • [ μ (∇u+ ∇uT)] +Fγ (6) Fγ=γκ n →δ+δ[∇γ−n →(n →• ∇γ)] (7) γ=max(0.031,0.245e−0.00246|j →•n →|)(8) Electric current density ∇ • ( − σ ∇φ) = − β Kφ (9) j →= − σ ∇φ (10) Fig. 1. (a) Interfacial tension dependence on the electric current density; (b) schematic presentation of the computational domain and boundaries. Experimental results are reproduced from ref. [17]. E. Karimi-Sibaki et al. Electrochemistry Communications 155 (2023) 107567 3 GHz ranges from one to two hours. All of the parameters used in the calculations are listed in Table 2, where material properties are obtained from ref. [27]. 3. Results and discussions We adopt the same set of geometrical and operational parameters as Li et al. [25] to investigate the temporal evolution of droplet shape. Due to the inherently transient nature of the system, it is challenging to effectively demonstrate field structures. As a result, transient simulation results are presented as Supplementary Material, namely “Droplet oscillation.avi” and “Droplet elongation.avi”. These animations offer a dynamic representation of the evolving droplet behavior. Readers are encouraged to utilize these resources to gain a deeper understanding of the complex interplay between field structures. Additionally, we present selected snapshots at different times, shown in Figs. 2 and 3, to provide insights into the dynamics of the system. In the study of Li et al. [25], it was observed that when a very low applied DC voltage (e.g., 0.5 V) was utilized, the droplet failed to reach the cathode and instead flattened and ceased motion after a certain period of time. In line with these findings, our simulation results, including volume fraction, velocity, electric potential, and electric current density, as shown in Fig. 2(a),(b),(c), exhibit similar behavior. The magnitude of vortical flow generated within the droplet, resulting from Marangoni effects, is insufficient to induce substantial deformation or sustain persistent deformation of the droplet. On the other hand, Li et al. [25] reported an interesting observation indicating that, under the application of a relatively moderate DC voltage of 5 V, the droplet exhibited a smooth and continuous flow from the anode, situated at the center of the chamber, towards the cathodes positioned at the side walls. As shown in Fig. 3(a),(b),(c), the elongation of the liquid metal protrusions towards the cathodes can be ascribed to the spatial distribution of the electric current density field, which induces significant variations in the interfacial tension. Consequently, a vigorous vortical flow emerges throughout the interface within the droplet, exerting a continuous propulsive force that drives the liquid metal towards the cathodes. To facilitate a more detailed analysis of the findings, the contour plot of the normalized current density (j*=j jmax) at the metal–electrolyte interface is illustrated in Fig. 4(a) and represents a valuable tool for tracking and visualizing the temporal evolution of the system. jmax represents the maximum current density (ca. 850 A m−2), beyond which the surface tension remains consistently minimal and invariant (see Fig. 1 (a)). At very low applied potential (0.5 V), the change in normalized current density and, consequently, the interfacial tension is less than 5%, as shown in Fig. 4(a). This indicates that the surface tension is sufficiently high to maintain the overall shape of the droplet, irrespective of any existing gradients. Conversely, as the normalized current density and, consequently, interfacial tension undergo variations ranging from 1 to 100% along the interface, a substantial gradient in interfacial tension force emerges, under an applied voltage of 5 V, as shown in Fig. 4(a). This gradient exerts a significant influence on the droplet’s behavior, affecting its deformation and overall dynamics. Fig. 4(b) illustrates temporal evolution of the normalized maximum width of the metal (W*=W 2Rd). Under a low applied voltage (0.5 V), we observe small fluctuations initially, followed by a relatively constant W*, indicating minimal variation in the maximum width of the metal. However, under a high applied voltage (5 V), the increase in W* initially progresses at a slow rate, but later undergoes a significant acceleration, resulting in a substantial growth in the droplet spreading rate. Our simulation result concurs with the experimental observation of Table 2 Parameters used in our calculations. Parameter Ra[mm]0.5 Rc[mm]1 Rv[mm]25 Rd[mm]8 Vapp[V]0.5 and 5 σ e[Sm−1]2.0 σ m[Sm−1]3.4 ×10 6 ρ e[kgm−3]2130 ρ m[kgm−3]6250 μ e[kgm−1s−1]0.004 μ m[kgm−1s−1]0.002 Fig. 2. With an applied voltage of 0.5 V, snapshots of the field structures at different times showing volume fraction, velocity, electric potential, and electric current density. Vectors of velocity, isolines of electric potential, and current density streamlines are provided. The arrow indicates time advancement. E. Karimi-Sibaki et al. Electrochemistry Communications 155 (2023) 107567 4 Li et al. [25], as shown in Fig. 4(c), which clearly illustrates the elongation of the liquid metal protrusion towards the cathodes. The 2D symmetrical nature of our simulation model may introduce differences compared to the referenced experiment [25]. Neglecting gravitational effects, thermal effects [28,29], droplet flattening, and droplet wetting in our model could impact the droplet’s shape and dynamics during elongation. Additionally, the experimental setup in [25] has some nonsymmetrical aspects, such as the droplet not being initially perfectly centered, contributing to the observed discrepancies. The current investigation sheds light on the transport phenomena involved in the shape evolution of liquid metal under the influence of Marangoni effects. Through the analysis of simulation results, valuable insights can be gained regarding the underlying mechanisms driving the electrolytic flow, particularly in relation to the interfacial tension gradient. The proposed VOF-based model holds significant potential for assisting industries to obtain a deeper understanding of the impact of operational parameters on the overall performance of the process. 4. Summary In this study, we propose a comprehensive multiphase model utilizing the volume of fluid (VOF) method to simulate the electrochemically-induced Marangoni flow phenomenon. We thoroughly investigated the transport phenomena, including flow, electric potential, and electric current density, within a NaOH electrolyte and a eutectic gallium–indium (EGaIn) metal droplet. The electrochemical oxidation of the anodic EGaIn droplet generates a gradient of interfacial tension along the metal–electrolyte interface, thereby inducing a flow. To accurately capture the transient variation of interfacial tension as a function of electric current density at the interface, we propose a novel algorithm. The model is then applied to investigate the influence of applied potential on the system behavior. At very low applied potential, where the interfacial tension gradient is minimal, the droplet exhibits slight oscillations while maintaining its original shape. However, with the application of a moderate potential, the droplet undergoes Fig. 3. With an applied voltage of 5 V, snapshots of the field structures at different times showing volume fraction, velocity, electric potential, and electric current density. Vectors of velocity, isolines of electric potential, and current density streamlines are provided. The arrow indicates time advancement. Fig. 4. (a) The contour plot of normalized current density captures the evolution of the droplet shape for applied voltages of 0.5 and 5 V; (b) temporal evolution of the normalized maximum width of the metal; (c) calculated droplet elongation is compared with experimental result reproduced from ref. [25], which is available under Creative Commons By Attribution 4.0 (CC BY 4.0) license. E. Karimi-Sibaki et al. Electrochemistry Communications 155 (2023) 107567 5 elongation towards the cathodes. To complement our findings, transient results are visualized in the form of animations in the supplementary material (Appendix B). CRediT authorship contribution statement E. Karimi-Sibaki: Conceptualization, Methodology, Investigation, Software, Writing – original draft. A. Vakhrushev: Visualization, Writing – review & editing. A. Kadylnykova: Investigation. M. Wu: Supervision. A. Ludwig: Supervision. J. Bohacek: Methodology, Software. A. Kharicha: Conceptualization, Methodology, Supervision, Writing – review & editing, Project administration. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Data availability The datasets used and/or analyzed during the current study are available from the corresponding author upon reasonable request. Acknowledgements The authors acknowledge financial support from the Austrian Federal Ministry of Economy, Family and Youth and the National Foundation for Research, Technology and Development within the framework of the Christian-Doppler Laboratory for Metallurgical Applications of Magnetohydrodynamics. Appendix Appendix A:. Nomenclature a Mesh element size, m Fγ Surface tension force, N m−3 j →Electric current density,A m−2 jmax Maximum electric current density,A m−2 K Permeability n →(nx,ny)Unit normal vector p Pressure,Pa Ra Radius anode, m Rc Radius cathode, m Rv Radius vessel, m Rd Radius droplet, m Sφ Electric potential sink term,A m−3s−1 t Time,s u →Velocity vector,m s−1 V Volume of computational cell,m3 Vapp Applied voltage,V Vm Metal volume in computational cell,m3 W Maximum width of metal, m x →Vector coordinates, m β Volume fraction of deposit γ Interfacial tension,N m−1 δ Dirac Delta function,m−1 κ Curvature,m−1 ρ Density,kg m−3 ρ m Metal density,kg m−3 ρ e Electrolyte density,kg m−3 μ Dynamic viscosity,kg s−1m−1 μ m Metal viscosity,kg s−1m−1 μ e Electrolyte viscosity,kg s−1m−1 φ Electric potential,V σ Electrical conductivity, S m −1 σ e Electrical conductivity of electrolyte, S m −1 σ m Electrical conductivity of metal, S m −1 Δt Time step size, s ()* A normalized parameter Appendix B. Supplementary material Supplementary data to this article can be found online at https://doi.org/10.1016/j.elecom.2023.107567. References [1] A.K. Selva Kumar, Y. Zhang, D. Li, R.G. Compton, Electrochem. Commun. 121 (2020), 106867. [2] D.H. Kelley, T. Weier, Appl. Mech. Rev. 70 (2018), 020801. [3] A. Kharicha, E. Karimi-Sibaki, M. Wu, A. Ludwig, J. Bohacek, Steel Res. Int. 89 (2018) 1700100. [4] T. Utigard, J.M. Toguri, Metall. Trans. B 17 (1986) 547–552. [5] K.C. Mills, B.J. Keene, R.F. Brooks, A. Shirali, Philosophical Transactions of the Royal Society of London Series A: Mathematical, Physical and Engineering Sciences 356 (1998) 911–925. E. Karimi-Sibaki et al. Electrochemistry Communications 155 (2023) 107567 6 [6] X. Yang, D. Baczyzmalski, C. Cierpka, G. Mutschke, K. Eckert, PCCP 20 (2018) 11542–11548. [7] Y. Lin, J. Genzer, M.D. Dickey, Adv. Sci. 7 (2020) 2000192. [8] M.S. Anwar, F. Ehlers, A. Bangert, Electron. Lett 58 (2022) 617–619. [9] G. Li, J. Du, A. Zhang, D. Lee, J. Appl. Phys. 126 (2019), 084505. [10] M. Song, K. Kartawira, K.D. Hillaire, C. Li, C.B. Eaker, A. Kiani, K.E. Daniels, M. D. Dickey, PNAS 117 (2020) 19026–19032. [11] C.-L. Song, Y. Tao, W.-Y. Liu, Y.-C. Chen, R. Xue, T.-Y. Jiang, B. Li, H.-Y. Jiang, Y.- K. Ren, Phys. Rev. E 105 (2022), 025102. [12] X. Li, L. Cao, B. Xiao, F. Li, J. Yang, J. Hu, T. Cole, Y. Zhang, M. Zhang, J. Zheng, S. Zhang, W. Li, L. Sun, X. Chen, S. Tang, Adv. Sci. 9 (2022) 2105289. [13] C.B. Eaker, D.C. Hight, J.D. O’Regan, M.D. Dickey, K.E. Daniels, Phys. Rev. Lett. 119 (2017), 174502. [14] H.J. Lee, D.J. Fermın, R.M. Corn, H.H. Girault, Electrochem. Commun. 1 (1999) 190–193. [15] M. Song, K.E. Daniels, A. Kiani, S. Rashid-Nadimi, M.D. Dickey, Advanced Intelligent Systems 3 (2021) 2100024. [16] M.R. Khan, C.B. Eaker, E.F. Bowden, M.D. Dickey, PNAS 111 (2014) 14047–14051. [17] K.D. Hillaire, M.D. Dickey, K.E. Daniels, ArXiv Preprint ArXiv:2011.10016 (2020). [18] C.W. Hirt, B.D. Nichols, J. Comput. Phys. 39 (1981) 201–225. [19] H.K. Versteeg, W. Malalasekera, An Introduction to Computational Fluid Dynamics: The Finite Volume Method, 2nd ed, Pearson Education Ltd, Harlow, England; New York, 2007. [20] J.U. Brackbill, D.B. Kothe, C. Zemach, J. Comput. Phys. 100 (1992) 335–354. [21] E.J.F. Dickinson, H. Ekstr¨ om, E. Fontes, Electrochem. Commun. 40 (2014) 71–74. [22] E. Karimi-Sibaki, A. Kharicha, M. Wu, A. Ludwig, J. Bohacek, J. Electrochem. Soc. 165 (2018) E604–E615. [23] E. Karimi-Sibaki, A. Kharicha, A. Vakhrushev, M. Wu, A. Ludwig, J. Bohacek, Electrochem. Commun. 112 (2020), 106675. [24] E. Karimi-Sibaki, A. Vakhrushev, M. Wu, A. Ludwig, J. Bohacek, A. Kharicha, J. Electrochem. Soc. 170 (2023), 072501. [25] M. Li, H.M.C.M. Anver, Y. Zhang, S.-Y. Tang, W. Li, Micromachines 10 (2019) 209. [26] D.J. Benson, Appl. Mech. Rev. 55 (2002) 151–165. [27] Z. Zhao, S. Soni, T. Lee, C.A. Nijhuis, D. Xiang, Adv. Mater. 35 (2023) 2203391. [28] J. Massing, G. Mutschke, D. Baczyzmalski, S.S. Hossain, X. Yang, K. Eckert, C. Cierpka, Electrochim. Acta 297 (2019) 929–940. [29] A.M. Meulenbroek, A.W. Vreman, N.G. Deen, Electrochim. Acta 385 (2021), 138298. E. Karimi-Sibaki et al.