DOI: 10.3303/CET25117164 Paper Received: 29 December 2024; Revised: 9 March 2025; Accepted: 18 May 2025 Please cite this article as: Zhou F., Jiang Y., Shi Y., Gelfgat A., 2025, A Computational Model for Two-Phase Fluid Flow, Chemical Reactions, and Heat Transfer in LAA-SOFC Fuel Cell , Chemical Engineering Transactions, 117, 979-984 DOI:10.3303/CET25117164 CHEMICAL ENGINEERING TRANSACTIONS VOL. 117, 2025 A publication of The Italian Association of Chemical Engineering Online at www.cetjournal.it Guest Editors: Fabrizio Bezzo, Flavio Manenti, Gabriele Pannocchia, Almerinda di Benedetto Copyright © 2025, AIDIC Servizi S.r.l. ISBN 979-12-81206-17-5; ISSN 2283-9216 A Computational Model for Two-Phase Fluid Flow, Chemical Reactions, and Heat Transfer in LAA-SOFC Fuel Cell Fangzhe Zhoua, Yidong Jiangb, Yixiang Shia, Alexander Gelfgatc a Key Laboratory for Thermal Science and Power Engineering of the Ministry of Education, Department of Energy and Power Engineering, Tsinghua University, Beijing, 100084, People’s Republic of China b Institute for New Energy Materials & Engineering, School of Materials Science & Engineering, Fuzhou University, Fuzhou, Fujian Province, 350108, People’s Republic of China c School of Mechanical Engineering, Faculty of Engineering, Tel-Aviv University, Israel, 677980 gelfgat@tau.ac.il Preliminary results of computational modelling of a combined electrochemical, fluid dynamic, and heat transfer processes taking place in a liquid antimony anode (LAA) of a solid oxide fuel cell (SOFC) are reported. The model includes computation of two-phase fluid flow of Sb and Sb2O3, electric field and current, heat transfer that takes into account the Joule heating and the heat released or consumed by chemical reactions, and production of Sb and Sb2O3 by reduction and oxidation reactions, respectively. Following motion of interface separating the two phases, we monitor changes in the velocity, temperature and electric fields that take place during cell working process. 1. Introduction Solid oxide fuel cell (SOFC) can directly convert chemical energy into electrical energy. However, solid coke formation at the porous anode surface may block the reactive sites and thus hinder the anodic reaction. Liquid metal anodes have better transport properties and hold stable even with solid carbon. Compared with other liquid metals, such as Ga, Sn, Pb and In, antimony (Sb) and antimony oxide (Sb2O3) are both in the liquid state at the SOFC working temperature, avoiding solid oxide formation at the anode-electrolyte interface. This motivates present study of the liquid antimonide anode solid oxide fuel cell (LAA-SOFC). Several researchers have attempted to simulate the flow pattern in the LAA-SOFC. Jiang et al. (2024) used a mixture model to simulate the two-phase flow in the LAA and analyses the electrochemical performance of the fuel cell, but it cannot represent the real pattern and have limitations. In the present study we present a single computational model, which includes two-phase fluid flow of Sb and Sb2O3, electric field and current, heat transfer that takes into account the Joule heating and the heat released or consumed by chemical reactions, and production of Sb and Sb2O3 by reduction and oxidation reactions, respectively. To the best of our knowledge, such modelling was never conducted for the LAA-SOFC fuel cells. The corresponding time-dependent partial differential equations are discretized using the finite volume method, and then are integrated in time by the second-order projection method. Motion of interface separating the two phases is modelled using the volume of fluid (VOF) method. The computations start from an initial state, in which most of the anode volume is occupied by Sb, while a small portion of its upper part is filled by lighter Sb2O3. The oxidation reaction takes place at the electrolyte-anode boundary. The reduction reaction happens inside the upper part of the anode, where carbon fuel is assumed to be supplied. It is shown that during the time evolution, the amounts of Sb and Sb2O3 together with the interface arrive to the time-dependent asymptotic regime, in which production of Sb and Sb2O3 become balanced over the oscillation period. The velocity, temperature, concentration and electric potential field are monitored during the whole time evolution process. 979 mailto:gelfgat@tau.ac.il 2. Description of the model An axisymmetric model of the LAA-SOFC cell is sketched in Figure 1. The detailed description of the cell is given in Jiang et al. (2023) and references therein. The cell working cycle includes supply of air to the cathode, where oxygen ions are formed and transported across the electrolyte. The liquid metallic Sb in the anode is oxidized by oxygen ions at the anode-electrolyte interface, which is described by Eq(1). Sb + 3 2 O2− → 1 2 Sb2O3 + 3e− (1) Electrons are released to the external circuit to power the load. Antimony oxide Sb2O3 is then reduced by carbon fuel according to Eq(2). 2 3 Sb2O3 + C → 4 3 Sb + CO2 (2) This reaction takes place at the upper part of the cell, as is shown in Figure 1. Since the density of Sb2O3 is smaller than that of Sb, 𝜌Sb = 6420 kg/m3, 𝜌Sb2O3 = 5600 kg/m3, the antimony oxide always tends to be located above the pure antimony To model the process, we define the initial state as a layer of Sb2O3 above a layer of Sb, with the prescribed height ratio, and observe the time corresponding evolution in time. Figure 1. Sketch of the axisymmetric model of LAA-SOFC cell. 1 – cathode, 2 – solid electrolyte, 3 – interface between liquid Sb and Sb2O3. Table 1: Parameters of the simulation 3. Mathematical model The following assumptions are made for the mathematical model described below: • Only flow of Sb and Sb2O3 in the anode is considered. • The Boussinesq approximation is applied. • The flow in both phases is laminar. • Mass transfer of the oxygen and the carbon fuel is not included in the model. Given reaction rates are used to account for production and consumption of Sb and Sb2O3. • The boundaries describing electrolyte surface, bottom and sidewall are equipotential. Parameter Value Parameter Value Cell radius 0.04 m Sb density 6420 kg/m3 Cell height 0.1 m Sb2O3 density 5600 kg/m3 Electrolyte radius 0.01 m Sb viscosity 0.0012 Pa˖s Electrolyte height 0.05 m Sb2O3 viscosity 0.25 Pa˖s Carbon fuel thickness 0.004 m Sb thermal conductivity 25.9 W/m˖K Average current density 2000 A/m2 Sb2O3 thermal conductivity 1.72 W/m˖K Initial temperature 973 K Sb electric conductivity 1.1˖106 Ω-1m-1 Heat generation by oxidation -106.8 kW/mol Sb Sb2O3 electric conductivity 2.7 Ω-1m-1 Heat generation by reduction 32.2 kW/mol Sb 980 3.1 Fluid dynamics and heat transfer Within the Boussinesq approximation, we assume that the densities remain constant everywhere except the buoyancy terms. Choosing 𝝁𝟏 𝝆𝟏𝑹⁄ , 𝑹𝟐𝝆𝟏 𝝁𝟏⁄ , and 𝝁𝟏 𝟐 𝑹𝟐𝝆𝟏⁄ to be the scales of length, velocity, time, and pressure, respectively, rendering the temperature dimensionless by 𝜽 = (𝑻 − 𝑻𝒎𝒊𝒏) (𝑻𝒎𝒂𝒙 − 𝑻𝒎𝒊𝒏)⁄ , we arrive to the non-dimensional form of the following equations [ 𝜕𝑢𝑗 𝜕𝑡 + (𝑢𝑗 ∙ ∇)𝑢𝑗] = − 𝜌2,0 𝜌21𝜌𝑗,0 ∇𝑝𝑗 + 𝜇𝑗 𝜌2,0 𝜌21𝜌𝑗,0 𝜇21𝜇𝑗 𝜇2 𝐷𝑖𝑣(𝜏) + 𝐺𝑎𝑒𝑧 + 𝐺𝑟 𝛽21𝛽𝑗 𝛽2 𝜃𝑗𝑒𝑧 + 𝐶𝑎 𝜌2,0 𝜌21𝜌𝑗,0 (∇ ∙ 𝑛)δ(𝑉) (3) ∇ ∙ 𝑢𝑗 = 0 (4) 𝜕𝜃 𝜕𝑡 + (𝑢𝑗 ∙ ∇)𝜃𝑗 = 1 𝜌21𝑐𝑝,21 1 𝑃𝑟 [𝜅21∆𝜃𝑗 + 𝐽𝑜 𝜎21 𝑗2] + 𝑆𝑟𝑒𝑎𝑐𝑡𝑖𝑜𝑛𝑠 (5) where the dimensionless governing parameters are the density, viscosity, thermal expansion, and thermal diffusivity ratios 𝝆𝟐𝟏 = 𝝆𝟐 𝝆𝟏⁄ , 𝝁𝟐𝟏 = 𝝁𝟐 𝝁𝟏⁄ , 𝜷𝟐𝟏 = 𝜷𝟐 𝜷𝟏⁄ , 𝜿𝟐𝟏 = 𝜿𝟐 𝜿𝟏⁄ , and 𝒄𝒑,𝟐𝟏 = 𝒄𝒑,𝟐 𝒄𝒑,𝟏⁄ , the Prandtl number 𝑷𝒓 = 𝝁𝟏 𝜶𝟏𝝆𝟏⁄ , the Grashof number 𝑮𝒓 = 𝒈𝜷𝟏(𝑻𝒎𝒂𝒙 − 𝑻𝒎𝒊𝒏)𝑹𝟑𝝆𝟏,𝟎 𝟐 𝝁𝟏 𝟐⁄ , the capillary number 𝑪𝒂 = 𝜸𝟎𝑹𝝆/𝝁𝟐, the Galileo number is 𝑮𝒂 = 𝒈𝝆𝟏,𝟎 𝟐 𝑹𝟑/𝝁𝟏 𝟐, and the Joule number 𝑱𝒐 = 𝝈𝟏,𝒎𝝋𝟎 𝟐 𝜿𝟏.𝒎(𝑻𝒎𝒂𝒙 − 𝑻𝒎𝒊𝒏)⁄ . Here 𝝋𝟎 is the scale of electric potential. Note, that equations for the fluid 1 are the same as they would be for a single phase buoyancy convection flow, while ratios of thermophysical properties necessarily appear in the equations describing the flow in the fluid 2. The Dirac delta function 𝜹(𝑽) must be smoothed along with the Heaviside function (see below). Since the electric field is assumed to be irrotational, it is described by the electric potential 𝝓, and the electric current is defined by the Ohm’s law, 𝒋 = −𝝈𝛁𝝓. Assuming electro neutrality of both liquid phases, the electric potential is obtained from the equation ∇ ∙ (𝜎∇𝜙) = 0 (6) For the boundary conditions, we assume an overpotential V0 at the electrolyte-anode interface, and zero potential at the other boundaries. 3.2 Volume of fluid interface tracking For the details on the Volume-of-Fluid (VOF) method and numerical schemes applied, the reader is referred to Tryggvason et al. (2011) and Patel & Natarajan (2015). The interface between the two liquids is described by the volume of fluid (VOF) function, defined as 𝑉 = { 0 𝑖𝑛 𝑙𝑖𝑞𝑢𝑖𝑑 1 0.5 𝑜𝑛 𝑡ℎ𝑒 𝑖𝑛𝑡𝑒𝑟𝑓𝑎𝑐𝑒 1 𝑖𝑛 𝑙𝑖𝑞𝑢𝑖𝑑 2 (7) All liquids properties are described as, e.g., �̃� = 𝜌1 (1 − 𝐻 (𝑉 − 1 2 )) + 𝜌2𝐻(𝑉 − 1 2 ) (8) where 𝐻(𝑉) is the Heaviside step function. The VOF function 𝑉 is advected by the equation 𝜕𝑉 𝜕𝑡 + (𝒖 ∙ ∇)𝑉 = 𝜕𝑉 𝜕𝑡 + ∇ ∙ (𝒖𝑉) = 0 (9) With the addition of VOF approach, the two-phase flow is described as a “single fluid” with variable properties. The Eqs. (3) and (5) are replaced by (𝝉 is the viscous stress tensor) [ 𝜕𝒖 𝜕𝑡 + (𝒖 ∙ ∇)𝒖] = − 1 �̃� ∇𝑝 + 1 �̃� 𝐷𝑖𝑣(𝜇𝝉) + 𝐺𝑎𝒆𝑧 + 𝐺𝑟�̃�𝜃𝒆𝑧 + 𝐶𝑎 �̃� (∇ ∙ 𝒏)δ(𝑉) (10) 𝜕𝜃 𝜕𝑡 + (𝒖 ∙ ∇)𝜃𝑗 = 1 �̃��̃�𝑝 1 𝑃𝑟 [𝑑𝑖𝑣(𝜅 𝑔𝑟𝑎𝑑𝜃) + 𝐽𝑜 �̃� 𝒋2] + 𝑆𝑟𝑒𝑎𝑐𝑡𝑖𝑜𝑛𝑠 (11) while the continuity Eq. (4) remains unchanged. The parameters with a tilde are defined for the ratios of material parameters, e.g., �̃� = (1 − 𝐻 (𝑉 − 1 2 )) + 𝜌21𝐻(𝑉 − 1 2 ) (12) To finalize the formulation, we consider so-called “diffuse interface”, in which the interface has a finite thickness 𝜉. The Heaviside function is smeared over the finite thickness interface as proposed by Tryggvason et al. (2011). 981 The normal unit vector in Eq. (10) is then defined as 𝒏 = ∇𝑉/|∇𝑉|, so that its divergence ∇ ∙ 𝒏 is equal to the main curvature of the interface. We impose no-slip boundary conditions on all solid boundaries. At the interface of the electrolyte we impose the Beavers-Joseph condition 𝜕𝒗/𝜕𝑛 = 𝑁𝒗, where 𝑁 is a slip coefficient. At large 𝑁, this condition approaches the no-slip one. The heat generation and consumption by chemical reactions are described by two separate source terms defined as 𝑆𝑜𝑥𝑖𝑑𝑎𝑡𝑖𝑜𝑛 = −𝑆𝑜𝑗𝑛𝛿(𝜙), 𝑆𝑟𝑒𝑑𝑢𝑐𝑡𝑖𝑜𝑛 = 𝑆𝑟�̇�, 𝑆𝑜 = 𝐻1 3𝑧𝑂𝐹 𝜎𝑆𝑏𝜑0𝑅 𝑐𝑝,𝑆𝑏𝜇𝑆𝑏(𝑇𝑚𝑎𝑥−𝑇𝑚𝑖𝑛) , 𝑆𝑟 = 𝜌𝑆𝑏2𝑂3𝐻2 𝑀𝑆𝑏2𝑂3𝜌𝑆𝑏𝑐𝑝,𝑆𝑏(𝑇𝑚𝑎𝑥−𝑇𝑚𝑖𝑛) (13) where 𝐻1 and 𝐻2 are heat generation intensities listed in Table 1. 3.3 Chemical reactions Mass flux per unit area of the oxygen ions generated at the electrolyte surface is governed by the Faraday’s law in the form of Beale (2005) �̇�𝑂2− = − 𝑀𝑂 𝑧𝑂𝐹 𝑗𝑛d𝛤 (14) where 𝑀𝑂 is the oxygen atom weight in kg/mol, 𝑧𝑂 = −2 is the valency, 𝑗𝑛 is the current density, and 𝑗𝑛𝑑𝛤 is the total current through an area 𝑑𝛤. Total number of 𝑂2− ions produced per unit time per unit area is (𝑁𝐴 is the Avogadro number) �̇�𝑂 = 𝑁𝐴 �̇�𝑂2− 𝑀𝑂 = − 𝑁𝐴 𝑧𝑂𝐹 𝑗𝑛 (15) Assume that the volume of a grid cell adjacent to the electrolyte surface is 𝑑(𝑉𝑜𝑙) = 𝑑𝛤 ∙ 𝑑𝑙. Since the VOF function is the volume fraction of Sb2O3, d�̇� = ΔΩ̇𝑆𝑏2𝑂3 d(𝑉𝑜𝑙) = − 𝑀𝑆𝑏2𝑂3𝑆𝑏 𝑧𝑂ρ𝑆𝑏2𝑂3 𝐹 𝑗𝑛 d𝑙 (16) Change of the Sb2O3 volume fraction by the reduction reaction is described by �̇� = −𝜈𝑆𝑏2𝑂3 𝑉𝑘𝑓𝑢𝑒𝑙 , 𝑘𝑓𝑢𝑒𝑙 = 3.36 × 1015𝑒𝑥𝑝 (− 371441 𝑅𝑇 ) (17) where 𝜈𝑆𝑏2𝑂3 is the stoichiometric coefficient of Sb2O3, which depends of chemical formula of the fuel. In the reported computations, we assume that the reduction reaction takes place inside the upper layer of width 𝑑𝑑𝑢𝑒𝑙, which is taken to be 4 mm. Change of the VOF function during the time step ∆𝑡 is given by 𝑉(𝑡 + ∆𝑡) = 𝑚𝑎𝑥(𝑉(𝑡) − ∆𝑉, 0), ∆𝑉 = 𝜈𝑆𝑏2𝑂3 𝑘𝑓𝑢𝑒𝑙𝑉∆𝑡 (18) 4. Computational procedure The problem is discretized by the finite volume method defined on the staggered non-uniform grid that can be stretched near the boundaries. The finite volume schemes and the stretching function are the same as in Gelfgat (2007). The convective terms of velocity are discretized by a conservative scheme, the convective terms of the temperature by the second order UPWIND scheme, and the convective terms in (10) by the interface capturing scheme proposed in Patel & Natarajan (2019). The time derivatives are discretized by the second order three time levels backward scheme. The time integration is performed by the velocity and pressure splitting projection method. For solution of the systems of linear equations we apply BiCG(stab) method. In cases it saturates, the code switches to the GMRES method. The results presented are obtained on 𝟓𝟎 × 𝟏𝟎𝟎 r- and z- nodes grid. The whole numerical process is carried out using a specially designed Fortran code. 5. Results Parameters used for the first simulation are listed in Table 1. At this stage we assume that the average current through the cell, and therefore through the anode is constant and is equal 0.2 A/cm2. Due to relatively large electric conductivity of Sb (see Table 1), the voltage drop needed to support such an average current is about 7×10-5 V. Distributions of the electric potential and electric current are shown in Figure 2. Calculations are started when 10% of the upper part of the anode volume is occupied by Sb2O3, while all the rest is filled by Sb. The main result is presented as an animation of the whole process. Here we report several characteristic snapshots shown in Figure 3. Change of the blue and red colors show how amount and location of Sb and Sb2O3 change in time. Because the kinetics of the chemical reaction between fuel and Sb2O3 are better than the kinetics of the anodic electrochemical reaction, the proportion of Sb2O3 in the anode remains at 982 a lower level despite the restricted distribution of the fuel. The isotherms are shown by colored lines, and flow direction by arrows. Figure 2. Characteristic distribution of the electric potential (colour) and the electric current (arrows) in the cell. To arrive to the average current of 0.2 A/cm2 the electric potential difference between the electrolyte and the outer wall should be approximately 7×10-5 V. The black line shows position of the electrolyte-anode interface. Figure 3. Characteristic snapshots. Blue colour shows position of Sb, and red colour – position of Sb2O3.The isotherms are shown by coloured lines. Arrows show direction of flow. Area occupied by the cathode and electrolyte is filled by white colour. The isolines correspond to the levels of the dimensionless temperature 𝜃. Our present findings show that the Joule heating and the heat released or consumed by the chemical reactions do not produce any noticeable change in the temperature field due to the quick heat transfer and the low Joule heat of the liquid metal, while the temperature difference within a conventional SOFC single cell can reach 20 − 50℃. To allow for weak buoyancy convection we assume that the electrolyte and the outer boundary are at the temperature difference of 1℃. In fact, the main source of the flow, in the presented case, is the Sb and Sb2O3 density difference. The Sb2O3 phase always tends to ascend. At the top of the cell it reacts with the fuel, 983 producing Sb, which tends to descend. In very long time the system arrives to an asymptotic unsteady state, at which production of Sb2O3 by the oxidation reaction at the electrolyte interface becomes balanced with the production of Sb by the reduction reaction in the upper part of the cell. 6. Conclusions To the best of our knowledge, this is the first attempt to model two-phase flow in a Sb LAA-SOFC fuel cell using the VOF method for the interface tracking. The presented model takes into account heat transfer and chemical reactions and mimics how the cell arrives to an asymptotic state, which appears to by non-stationary. More work is planned to address the mass transfer, and to formulate more realistic boundary conditions. As the main sources of heat generation and voltage loss, the conduction of oxygen ions in the electrolyte and the overpotential of the electrochemical reactions will be included in the model calculation. Meanwhile, the heat absorption caused by the enthalpy change of the reaction between fuel and Sb2O3 will also be considered, which will lead to an increase of the temperature difference inside the anode. After completion of the latter consideration of 3D models is planned. Nomenclature u – velocity, m/s 𝝆𝟐𝟏– ratio of densities p – pressure, N/m2 𝝁𝟐𝟏– ratio of viscosities T – temperature, K 𝜿𝟐𝟏 – ratio of thermal conductivities θ – dimensionless temperature 𝒄𝒑,𝟐𝟏- ratio of heat capacities ρ – density kg/m3 𝝈𝟐𝟏 – ratio of thermal conductivities μ – viscosity, Pa˖s Ga – Galileo number β – thermal expansion coefficient Gr – Grashof number κ – thermal conductivity, W/m˖K Ca – Capillary number 𝒄𝒑 – heat capacity, J/kg˖K Pr – Prandtl number 𝝈 – electric conductivity, 1/Ω˖m Jo – Joule number j – electric current density, A/m2 NA – Avogadro number 𝜙 – electric potential, V F – Faraday constant V – VOF function �̇� – mass generation rate 𝜹 – Dirac delta function M – molar weight H – Heaviside function z - valency 𝑆𝑟𝑒𝑎𝑐𝑡𝑖𝑜𝑛𝑠 – heat generation by chemical reactions Acknowledgments This study was supported by Israel Science Foundation (ISF) grant 2979/23 (to A. Gelfgat) and National Natural Science Foundation of China (to Y. Shi). References Beale S., 2005, Numerical models for planar solid oxide fuel cells, Chapter In: Transport Phenomena in Fuel Cells, B. Sunden & M. Faghri (Ed.), WIT Press, UK, 43-82. Gelfgat A. 2007 Three-dimensional instability of axisymmetric flows, solution of benchmark problems by a low- order finite volume method, Int. J. Numer. Meths. Fluids, 54, 269-294. Jiang Y., Gu X., Shi J., Shi Y., Cai N., 2023, Co-generation of gas and electricity on liquid antimony anode solid oxide fuel cells for high efficiency, long-term kerosene power generation, Energy 263, 125758. Jiang Y., Liu C., Gu X., Shi Y., Yan W., Zhang J., 2024, Development of liquid antimony anode-based fuel cells: Effects of reaction-induced convection on mass transfer and electrochemical performance, Energy Convers. Manage., 319, 118874. Patel J. K., Natarajan G. T., 2015, A generic framework for design of interface capturing schemes for multi-fluid flows, Computers & Fluids, 106, 108-118. Tryggvason G., Scardovell, R., and Zalesky S., 2011, Direct numerical simulations of gas–liquid multiphase flows, Cambridge Univ Press, UK. 984 CET-vol117-b.pdf 288fangzhe.pdf A Computational Model for Two-Phase Fluid Flow, Chemical Reactions, and Heat Transfer in LAA-SOFC Fuel Cell