**Simulation of Solid Oxide Fuel Cell Anode in Aspen HYSYS—A Study on the E**ff**ect of Reforming Activity on Distributed Performance Profiles, Carbon Formation, and Anode Oxidation Risk**

#### **Khaliq Ahmed 1, Amirpiran Amiri 2,\* and Moses O. Tadé <sup>1</sup>**


Received: 11 November 2019; Accepted: 24 February 2020; Published: 27 February 2020

**Abstract:** A distributed variable model for solid oxide fuel cell (SOFC), with internal fuel reforming on the anode, has been developed in Aspen HYSYS. The proposed model accounts for the complex and interactive mechanisms involved in the SOFC operation through a mathematically viable and numerically fast modeling framework. The internal fuel reforming reaction calculations have been carried out in a plug flow reactor (PFR) module integrated with a spreadsheet module to interactively calculate the electrochemical process details. By interlinking the two modules within Aspen HYSYS flowsheeting environment, the highly nonlinear SOFC distributed profiles have been readily captured using empirical correlations and without the necessity of using an external coding platform, such as MATLAB or FORTRAN. Distributed variables including temperature, current density, and concentration profiles along the cell length, have been discussed for various reforming activity rates. Moreover, parametric estimation of anode oxidation risk and carbon formation potential against fuel reformation intensity have been demonstrated that contributes to the SOFC lifetime evaluation. Incrementally progressive catalyst activity has been proposed as a technically viable approach for attaining smooth profiles within the SOFC anode. The proposed modeling platform paves the way for SOFC system flowsheeting and optimization, particularly where the study of systems with stack distributed variables is of interest.

**Keywords:** SOFC; simulation; internal reforming; anode oxidation; carbon formation

#### **1. Introduction**

Fuel cells convert the chemical energy in fuel directly into electricity and heat, without combustion, leading to high efficiencies with low or even zero emissions. The SOFC is becoming a mature technology and can make the commercial breakthrough if cost targets can be met by achieving cost reductions through volume manufacturing, improved lifespan/performance, and lower cost materials [1–3]. Research and development in the last twenty years have led to significant advances in all areas of the technology including cell, seal, interconnect, and stack design, as well as peripheral components and the entire balance of plant (BoP) [4,5]. Manufacturing achievements have led to defect identification and minimization, quality control, and scale-up of stack components and the entire stack assembly manufacture.

The SOFC system is appropriate to operate on a pipeline fuel such as reticulated natural gas with its well-established supply infrastructure throughout the world. For such a fuel, minimal fuel processing is required, which includes desulphurization of the fuel to remove sulphur compounds that are naturally present in the hydrocarbon fuel source and those that are added as odorants to meet legislative requirements such as for natural gas, propane, and LPG for domestic applications. The preprocessing also includes a level of conversion of the hydrocarbon fuel, conventionally known as prereforming, which functions to convert the hydrocarbon feed such as natural gas to a hydrogen-rich mixture or a methane-rich mixture depending on the type of anode in the SOFC stack, i.e., noninternal reforming or internal reforming type. For an internal reforming type SOFC, where methane can be converted by steam reforming on the anode, it suffices to prereform the fuel to a level where all higher hydrocarbons (C1+) are converted leading to a mixture of methane, hydrogen, and carbon oxides with little or no conversion of methane [6,7]. For a noninternal reforming type SOFC, all hydrocarbon components including methane need to be fully converted to a mixture of hydrogen and carbon oxides. Owing to its high electrical efficiency, the SOFC technology results in reduced emissions of CO2 and is practically noise free. Furthermore, it is free of NOx emissions due to its relatively low operating temperatures. The SOFC system is particularly attractive as a combined heat and power generation (CHP) system, since the waste heat generated can be used to supply heat to a hot-water system which can be interfaced to the SOFC system [8].

A system-level flowsheet model of the SOFC system including the complete BoP is a useful platform for simulating the performance of the plant and for sizing of individual components of the BoP. Commercially available process simulation software such as Aspen Plus or Aspen HYSYS, PRO/II, etc., contains extensive thermodynamic and physical properties database and includes in-built modules for a number of components which are commonly used in a process plant, such as heat-exchangers of various types, reactors of various types, compressors, pumps, valves, separating columns, tanks, mixers, etc. It allows for energy optimization via heat and work integration of system components. However, it does not include a module for fuel cell reactions, i.e., it cannot directly account for the electrochemical reactions involving ions and electrons. There are two approaches for modeling SOFC-based systems with commercial process simulators. In one approach, the SOFC model is developed in a separate platform such as FORTRAN, VB, C++, MATLAB, etc. and then linked to the process simulator [9–15]. In another approach, the SOFC reactions are modeled using the equilibrium reactor module GIBBS [16,17]. Anderson et al. [17] modeled the SOFC as a combination of an isothermal plug flow reactor (PFR) module, to account for methane reforming kinetics on the anode, and a GIBBS reactor for the fuel cell reactions of hydrogen and CO oxidation. However, this was not a system-level model and focused on reactions and mass transport processes at the cell level. Using established theoretical and/or empirical correlations from literature, they tested the validity of their model by comparing their simulation results with those of others reported in the literature. Two main drawbacks of this work are the assumptions of isothermal conditions for internal reforming and use of a GIBBS reactor for the fuel cell reactions. In a real system the SOFC stack does not operate in an isothermal mode. There are two opposing contributions to stack temperature profiles in the case of an internal reforming anode. The endothermic steam reforming reaction absorbs heat from the gas stream which results in cooling of the stack and the fuel cell reaction(s) release heat which results in heating of the stack; the net effect is determined by the extents of these reactions.

A conversion reactor is more appropriate for representing the fuel oxidation reactions by setting the percent conversion equal to fuel utilization. An equilibrium approach using the GIBBS reactor does not allow setting of the reaction conversion to match fuel utilization. In this work, the internal reforming of methane via steam reforming and the accompanying water-gas shift (WGS) reaction is modeled via the PFR module with the kinetic expressions from literature [18,19] and the fuel cell reactions are modeled using the conversion reactor where the conversion is linked to the fuel utilization value calculated in a spreadsheet block. Another feature of the current work is that unlike the work of Anderson et al. [17], where the PFR is modeled as an isothermal reactor, in this work the energy stream of the PFR is linked to the cell in the spreadsheet block which calculates the heat generated by the fuel cell reaction and is available as reaction heat in a direct internal reforming SOFC. The axial temperature profile created in the PFR is therefore, representative of the temperature profile on an SOFC anode with direct internal reforming, as the coupling of the endothermic methane steam reforming (MSR) reaction and accompanying mildly exothermic WGS reaction with heat available from the fuel cell reaction is appropriately captured with this approach. The corresponding composition profile under current load cannot be generated within the PFR module as this module only works with kinetic schemes and it is not possible to add the fuel cell reactions to the PFR as conversion reactions based on fuel utilization. An option to include the reaction kinetics of the electrochemical oxidation of hydrogen is also not available in the software. Nevertheless, the reaction extents and accompanying heat exchanges can be calculated in the spreadsheet module and linked to the PFR module. Firstly, this allows generation of open-circuit composition profiles of the internally reformed gas which sets the boundary for the Nernst voltage profile under load, after accounting for the extents of the fuel cell reactions. Secondly, the current density and composition profiles can be calculated within the spreadsheet block using appropriate correlations.

Previous work [9–13,16] largely focused on the issues that can predict and improve the fuel cell operation in terms of current generation and voltage losses. For instance, the effect of air flow rate, steam to carbon ratio (S/C), current density, fuel utilization (Uf), inlet temperatures, or operating pressure have been extensively investigated. By contrast, in this work we have mainly focused on an analysis of processes that significantly affect anode performance and lifetime and consequently impact on the SOFC system as a whole. We analyzed anode performance for various levels of reforming activity. Three cases are considered: (i) Full reforming activity, (ii) 1/3rd reforming activity, and (iii) 1/6th reforming activity. Reduced reforming activity may be the result of engineered design [20] or may result from progressive degradation of the anode from poisoning or sintering due to nickel coarsening over the useful life of the stack [4], which extends the reaction zone and requires more of the anode segment from the leading edge to fully convert methane. Reforming kinetics reported by Ahmed and Föger [18] and WGS reaction kinetics reported by Tingey [19] were employed as the reaction rate details for PFR, leading to 1D pseudo-homogeneous results. The three different levels of activity were assigned by reducing the Arrhenius factor in the rate expression by the reduction factors 0.33 (~1/3rd) and 0.67 (~1/6th). In physical terms this signifies loss of reforming activity by poisoning or sintering. For these levels of activity, we assess the anode oxidation risk and carbon formation potential on the anode, both of which have severe life-limiting consequences on the anode [7].

This paper contributes to the SOFC research considering its two novel contents including; (i) the novel simulation methodology: The simulation approach proposed for complicated internal reforming SOFC process offers simplicity and calculation speed without compromising the internal operations details. This is of particular interest for SOFC system modeling and design where several operational concepts including heat/mass transfer and electrochemical and fuel reformation reactions interactively occur at wide time and length scales. (ii) The understanding of distributed reformation potentials in controlling SOFC performance profiles. The incremental reformation process is demonstrated to be a promising strategy to moderate the undesirable gradients of SOFC internal profiles. This is promising to achieve higher homogeneity in temperature and concertation profiles inside the SOFC stack that subsequently offers enhanced current and voltage profiles. This is crucially important for SOFC efficiency and durability. In this paper, we demonstrate internal fuel reformation as an opportunity not only for heat integration and external reformer cost reduction but also for thermal management goals that eventually results in fuel cell longevity.

#### **2. Simulation Methodology**

The SOFC stack is simulated in the following way: Internal reforming of methane on the SOFC anode is modeled in a PFR module, using the reforming kinetics of Ahmed and Föger [18] and the WGS reaction kinetics of Tingey [19]. The electrochemical conversion of hydrogen and carbon monoxide are modeled as chemical conversions in a conversion reactor module. The associated electrical aspects including cell voltage, air utilization, and fuel utilization at a given operating current are calculated in a spreadsheet tool of Aspen HYSYS. The spreadsheet block in Aspen HYSYS is essentially an Excel-based spreadsheet with features of exporting data to and importing data from other modules of the flowsheet. This provides facility for post-processing of the calculated or entered values of the process variables. Computations of cell voltage entailing calculations of electrical losses, both ohmic and overpotential using empirical correlations [21] and calculation of Nernst voltage from the concentrations of the reacted gas are carried out in the spreadsheet, with compositions imported from the PFR module. Similarly, the sum of electrical losses and entropy change, after allowing for losses to the surroundings from the stack, is calculated and this value is exported to the energy stream linked to the PFR, as heat available for direct internal reforming.

The average stack operating temperature is obtained by a trial-and-error method. An average stack temperature is assumed for calculating Nernst voltage, operating voltage, and electrical losses. The assumed temperature is then compared to the value returned by the PFR, which simulates internal reforming and utilizes heat generated from the fuel cell reaction to compute the temperature profile within the fuel cell and the composition profile of the internally reformed fuel. The amount of heat generated, i.e., the electrical and entropy losses are calculated at the assumed average temperature. These steps are iterated until agreement is reached between the assumed average cell temperature and the calculated average cell temperature.

The composition of utilized gas is obtained by matching the degree of conversion in the conversion reactor, with the fuel utilization level calculated in the spreadsheet, based on fuel flow rate and operating current. Since the electrochemical conversion of H2 and CO are modeled as combustion reactions, there is a temperature rise in the reactor due to the exothermicity of the combustion reactions. Since the exothermicity of the fuel cell reactions have already been accounted for by calculating the electrical losses and entropy change of the fuel cell reaction and entering this value as heat input to the internal reforming PFR, temperature rise in the conversion reactor is suppressed by use of the HYSYS object SET, to set the temperature of the stream leaving the conversion reactor to be the same as the temperature of the stream leaving the PFR.

Figure 1 and Table 1 show the integration of the modules representing the anode operation and the set of equations used in the calculation of all distributed variables within the stack—temperature, current density, and composition of all chemical species on the anode side.

**Figure 1.** Schematic of integration of the internal reforming reaction modules and the fuel cell reactor module.

The reactions taking place in the reaction modules, shown in Figure 1, are as follows. Two parallel reactions take place in the internal reforming module (PFR) including MSR and WGS, presented by Equations (1) and (2), respectively:

$$\text{CH}\_4 + \text{H}\_2\text{O} = \text{CO} + 3\text{H}\_2\tag{1}$$

$$\text{CO} + \text{H}\_2\text{O} = \text{CO}\_2 + \text{H}\_2\tag{2}$$

*Processes* **2020**, *8*, 268

Two parallel electrochemical reactions take place in the fuel cell reactor module (conversion reactor). The main reaction is hydrogen oxidation (Equation (3)) that may occur with the CO oxidation (Equation (4)), simultaneously.

$$\text{H}\_2 + 0.5\text{O}\_2 = \text{H}\_2\text{O} \tag{3}$$

$$\text{CO} + 0.5\text{O}\_2 = \text{CO}\_2 \tag{4}$$

**Table 1.** Modeling approach for solid oxide fuel cell (SOFC) performance approximation.


#### **3. Results and Discussion**

The base case is a system operating at 65% fuel utilization. Input parameters for the base case simulation and also estimated results returned by the flowsheet calculations are presented in Tables 2 and 3, respectively.


**Table 2.** Model input parameters for the base case simulation.


**Table 3.** Estimated SOFC performance results.

#### *3.1. Distributed Profiles*

A high degree of internal reforming of methane is desirable as it reduces the prereformer size/cost and contributes to the cooling of the cell and therefore, reduces the cooling air flow requirement. However, the internal fuel reformation process may undesirably cause higher gradients in the distributed variables profiles. Stack operation homogeneity is of crucial importance from both efficiency and lifetime viewpoints. Figure 2 shows when the kinetics of reforming is too rapid, the cell temperature drops sharply near the cell inlet preventing an even distribution of current density and temperature. The steep temperature gradient is caused by the cooling impact associated with the endothermic steam reforming of methane on the anode. Figure 3 shows the reforming in this case is so fast that it is completed within the first 40% of the cell. Clearly, the maximum temperature gradient can be reduced if the sharp drop in temperature can be avoided, improving all temperature-dependent profiles including species concentration, current density, and overpotentials distributions, all of which have interacting effects. An effective way of achieving this objective is by using catalysts with lower overall reforming activity in anode and/or designing an anode with progressively increasing local activity along the cell length. A reduction of the anode activity to 1/3rd of that of the fast reforming anode results in a relatively uniform temperature profile as shown in Figure 2.

**Figure 2.** Temperature profiles for various levels of reforming activity.

**Figure 3.** Methane concentration profiles for different reforming activity.

A key technical question is whether the fuel reformation should be fully accomplished in the first half of the cell length. If not, there are significant opportunities to design anode reforming activity targeting a more homogenous operation. Accordingly, the objective of achieving smoother temperature profiles in the operating cell must also be assessed from another viewpoint. The activities are compared to the basis case activity (fast reformation). For instance, 1/3rd activity means 66% less activity (lower rate) compared to the basis case. We have adjusted this via E (activation energy value) adjustment in the simulations. With activities lower than the basis case, the catalyst still may offer sufficient reformation as the methane consumption (Figure 3) and hydrogen generation profiles show. The speed of reformation, however, is declined while it is compensated by residence time (axial length). This is the distributed conversion in contrast to the sharp conversion seen in basis case. To further reduce the temperature differential, an anode with a reforming activity 1/6th of that of the base case fast reforming anode is considered. This results in a more uniform temperature profile, Figure 2, but comes at a cost of increasing methane slip as shown in Figure 3. While a more homogenous electrochemical reactor is achievable under distributed/progressive reforming conditions, an immediate concern relevant to the reduced activity is the possibility of reduction in the generated current due to less local H2 availability in the first 30%–40% of the cell length. In order to assess this rigorously, the overall current produced over the cell surface must be calculated via integration of the local currents. Since no variation of current in cell width direction is assumed, that is a reasonable assumption for co-flow cell, two-dimensional integration can be replaced by one dimensional integration over the cell length. In such a case, therefore, the modeling approach proposed in this paper suffices. The surface under current density profiles in Figure 4 compares the total current production under different reforming approaches. Even though current production considerably drops in the cell inlet region when fast reforming occurs, the overall current production variation for fast and slow reformation activities can be reasonably ignored, indicating that cell efficiency will not be compromised for homogeneity.

Methane reforming profile (Figure 3) gives an indication of how the anode exhaust stream might be post-processed. In a fast reforming method, methane is completely consumed inside the stack leaving anode tail gas mainly including, hydrogen, CO, and a significant amount of steam. In such a case anode gas recycling is an appropriate process strategy. For slow/distributed fuel reforming, anode tail gas might contain some methane and less hydrogen and CO compared to the fast reforming. This can be understood by interpretation of methane and current profiles simultaneously. As current

generations are equal, the H2 consumption in all cases are almost the same. For a given methane rate at inlet, therefore, higher methane fraction in cell outlet indicates lower amount of hydrogen and CO. This becomes even more considerable for an activity as low as 1/6th of the base case. Therefore, an after-burner can be designed in the system to achieve high quality heat from the tail gas. It may be concluded that the progressive fuel reforming offers some advantage for a CHP system in which both high quality heat and homogenous stack performance are desirable.

**Figure 4.** Current density profiles for different levels of reforming activity and Uf = 65%.

Stack exhaust temperature is reduced with progressive reforming activity compared to fast reforming, but the stack temperature profile for most of the anode length is higher as shown in Figure 2. While the latter is expected due to less steeper cooling along the cell, the reduction in stack exhaust temperature can be explained from a consideration of the associated current density profile. As shown in Figure 4, the current density profile is steepest for the fast reforming case resulting in higher joule heating effect. A lower joule heating effect results in lower increase in temperatures for the distributed reforming cases, to the point that the temperature at the end of the cell is lower for distributed reforming, even though the temperatures near the inlet are higher. A more uniform current density profile results in the case of lower reforming activity as expected which is beneficial from the point of stress reduction.

#### *3.2. Anode Oxidation*

The SOFC anode is susceptible to oxidation by steam in the reaction mixture according to Equation (5):

$$\text{Ni} + 2\text{H}\_2\text{O} = \text{Ni}(\text{OH})\_2 + \text{H}\_2 \tag{5}$$

In this paper, anode oxidation risk is determined on the basis of industrial experience [22] with nickel-based catalysts, where it reported finding that steam-to-hydrogen ratios greater than 6–8 increases the risk of nickel oxidation in nickel-based steam reforming catalysts. Analyses, therefore, are based on the local partial pressures as a characteristic indicator rather than estimation of rate for reaction 5. Distributed profiles show that the risk of anode oxidation is particularly high at high fuel utilization, where the partial pressure of steam in the reaction mixture is high, i.e., at high pH2O/pH2. Anode oxidation risk profile along the length of the cell is shown in Figure 5 for 65% and 75% fuel utilization with three levels of reforming activity. Anode oxidation risk is highest at the anode inlet for internal reforming anodes where enough hydrogen has not been generated (Figure 6). The risk

increases with slow reforming activity. At higher utilization, more hydrogen is utilized, further lowering the H2 content and increasing the H2O content, thereby increasing the risk of anode oxidation compared to lower utilization. In such a case, anode tail gas recycle may worsen the situation by introducing more steam upstream. This is an additional reason why tail gas in the slow internal reforming case is recommended to be burnt rather than being recycled. In general, anode gas recycling may change the pH2O and pH2 balance in favour of steam.

**Figure 5.** Anode oxidation risk profile for different reforming activity at 65% and 75% fuel utilization.

**Figure 6.** Open-circuit hydrogen concentration profiles for different reforming activity.

#### *3.3. Carbon Formation*

Carbon formation on the SOFC anode, as a challenging issue for SOFC degradation, might occur principally via two routes: Boudouard reaction and methane cracking. In this work, carbon formation activity is calculated based on the thermodynamical equilibrium estimated using the local temperature and compositions of the relevant species participating in the respective reactions. Carbon may form by cracking or dissociation of a methane molecule into carbon and hydrogen molecules according to the reaction in Equation (6):

$$CH\_4 = C + 2H\_2 \tag{6}$$

For this reaction, the carbon activity for dissociation of a methane (*a<sup>d</sup> <sup>c</sup>* ) can be calculated by using Equation (7):

$$a\_c^d = \mathcal{K}\_c^d \frac{p\_{CH\_4}}{\left(p\_{H\_2}\right)^2} \tag{7}$$

The Boudouard or carbon disproportionation reaction can be presented as Equation (8):

$$\text{'} \color{orange}{2} \color{orange}{CO} = \color{orange}{C} + \color{purple}{CO\_2} \tag{8}$$

Similar to the methane cracking, a carbon activity can be defined for Boudouard reaction (*ab <sup>c</sup>* ) by using the local concentration/partial pressure of the gases involved, as presented in Equation (9):

$$a\_c^b = K\_c^b \frac{\left(p\_{\text{CO}}\right)^2}{p\_{\text{CO}\_2}} \tag{9}$$

Simulated distribution of the carbon formation risks from Boudouard reaction and methane cracking are depicted in Figures 7 and 8, respectively. Figure 7 shows that under the operational conditions used in this study and fuel utilization ranging from 65% to 75%, the probability of carbon formation through Boudouard reaction is low and well below 1 regardless of the fuel reforming pattern, i.e., fast or distributed reformation patterns. This is primarily due to the fact that this reaction is not thermodynamically benefited by the elevated temperature specifically above 700 ◦C. The trend along the anode length can be explained with respect to the CO and CO2 concentrations profiles and the anode thermal behaviour. The risk is relatively high at lower fuel utilization levels due to relatively higher and lower levels of CO and CO2, respectively, in the anode gas mixture and also lower local temperature, all enhancing the Boudouard reaction chance to occur. Note that according to Equation (9), at any given temperature, the higher the ratio of (pCO) <sup>2</sup>/pCO2, the higher the carbon activity. Near the inlet, CO is low, as methane reforming has not progressed much. Further down, the ratio depends on how much CO is formed by reforming and how much CO2 is formed by WGS (Equation (2)). The WGS equilibrium is also affected by the current draw, as the WGS equilibrium shifts to the right as more H2 are consumed by the hydrogen electrochemical oxidation reaction (Equation (3)). It will also be affected by electrochemical oxidation of CO to CO2, but the extent of this reaction is generally small, as this reaction is not as fast as H2 oxidation. Moreover, due to equal stoichiometry, this effect is accounted for by the WGS equilibrium reaction (Equation (2)). As the WGS reaction generates one mole of hydrogen per mole of CO, the electrochemical and chemical balance is unaffected whether the CO conversion is modeled as a WGS or as electrochemical oxidation. In a recent work [23], the effect of different reaction kinetics and equilibrium of the methane steam reforming reaction and WGS reaction were shown to have significant effect on the concentration profiles along the cell length. Temperature change along the cell length also affects the carbon activity as it changes the value of the equilibrium constant.

Figure 8 shows a very high risk of carbon formation by methane cracking near the anode inlet where methane concentration is high. The risk increases for lower fuel utilization and for slower reforming activity. The carbon formation activity is calculated based on thermodynamic equilibrium; in practice, whether carbon will be formed will depend on the kinetics of the reactions involved.

**Figure 7.** Profile of carbon formation risk from Boudouard reaction for different reforming activity at 65% and 75% fuel utilization.

**Figure 8.** Profile of carbon formation risk from methane cracking for different reforming activity at 65% and 75% fuel utilization.

#### **4. Conclusions**

A model of a SOFC cell incorporating a 1D model of the anode has been developed in Aspen HYSYS. A salient feature of this model is its ability to predict simultaneous direct internal reforming on the anode and electrochemical reaction with these two reactions thermally integrated. In spite of being based on empirical correlations that makes modeling platform computationally efficient it successfully captures the distributed variables associated with stack level performance, particularly with respect to internal reforming kinetics and electrical performance. This has been achieved through interlinking

of the in-built PFR module with a spreadsheet block inside the Aspen HYSYS environment. Cases for various levels of reforming activity have been compared to demonstrate the effect and relative advantages and disadvantages in terms of temperature and current density profiles. Two technically challenging aspects of SOFC operation: Anode oxidation risk and carbon formation potential have been evaluated. The methodology proposed in this paper is flexible to deploy more detailed fundamental correlations such as explicit equations based on exchange current.

**Author Contributions:** Conceptualization, K.A., A.A., and M.O.T.; methodology, K.A. and A.A.; software, K.A. and A.A.; validation, K.A. and A.A.; formal analysis, K.A., A.A., and M.O.T.; investigation, K.A., A.A., and M.O.T.; resources, M.O.T.; data curation, K.A. and A.A.; writing—original draft preparation, K.A., M.O.T., and A.A.; writing—review and editing, M.O.T., K.A. and A.A.; visualization, K.A. and A.A.; supervision, M.O.T.; project administration, M.O.T. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Studies on Influence of Cell Temperature in Direct Methanol Fuel Cell Operation**

**R. Govindarasu 1,\* and S. Somasundaram <sup>2</sup>**


**\*** Correspondence: rgovind@svce.ac.in

Received: 2 January 2020; Accepted: 17 March 2020; Published: 19 March 2020

**Abstract:** Directmethanol fuel cells (DMFCs) offer one of the most promising alternatives for the replacement of fossil fuels. A DMFC that had an active Membrane Electrode Assembly (MEA) area of 45 cm2, a squoval-shaped manifold hole design, and a Pt-Ru/C catalyst combination at the anode was taken for analysis in simulation and real-time experimentation. A mathematical model was developed using dynamic equations of a DMFC. Simulation of a DMFC model using MATLAB software was carried out to identify the most influencing process variables, namely cell temperature, methanol flow rate and methanol concentration during a DMFC operation. Simulation results were recorded and analyzed. It was observed from the results that the cell temperature was the most influencing process variable in the DMFC operation, more so than the methanol flow rate and the methanol concentration. In the DMFC, real-time experimentation was carried out at different cell temperatures to find out the optimum temperature at which maximum power density was obtained. The results obtained in simulation and the experiment were compared and it was concluded that the temperature was the most influencing process variable and 333K was the optimum operating temperature required to achieve the most productive performance in power density of the DMFC.

**Keywords:** direct methanol fuel cell; methanol crossover; power density; catalyst; membrane electrode assembly

#### **1. Introduction**

In the present juncture of the energy crisis, it has become inevitable to find alternatives for the mainly exploited fossil fuels. Recent developments in the field of fuel cells have given encouraging results suggesting their possible use of replacing the conventional highly polluting, less efficient combustible engines [1,2]. In this context, the direct methanol fuel cell (DMFC) appears to be the most promising tool to provide power to portable electronic devices. This is an electrochemical cell that has advantages such as offering a simple and easy method to store fuel, a simple design and green emissions. DMFC is a subcategory of a proton exchange membrane fuel cell (PEMFC) in which methanol is used as fuel. The salient features of DMFCs are the ease of transport of methanol, low or zero emissions, reliability in operation and utilization of methanol directly as a fuel to convert chemical energy into electric power. These DMFCs are designed especially for portable applications, where energy and power density are more important than efficiency [3,4]. A simple DMFC consists of a methanol distributor, gasket, anode, membrane, cathode, and oxygen gas distributor.

Regarding the working principles of the DMFC system, methanol is oxidized to hydrogen ions (H<sup>+</sup>) and electrons (e- ) at the anode. The released electrons are transported from the anode to the cathode through an electrical circuit where power is withdrawn, and at the same time, the hydrogen ions travel to the cathode through the electrolyte membrane. At the cathode, both the electrons and hydrogen ions react with oxygen and produce water and heat [5].

Increasing demand for clean sources of energy and sustainable energy development has led to the exploration of alternative energy. One such activity was to develop a membrane-protected anode in a fuel cell [6]. This was conductedusingan anionic backbone of sulfonated polystyrene block–(ethylene–ran–butylene)–block polystyrene polymer on top of an anode which was used to increase the oxygen evolution. Another method developed to increase the enhancement of oxygen deposition was by selective oxidation using a cation-selective polymer material such as Nafion which improved the electrolysis and enhanced the oxygen evolution. The evolved oxygen can be used in a DMFC wherein electricity is generated [7,8]. In a direct methanol fuel cell, a solution of methanol is internally reformed. This is conducted with the help of a suitable catalyst, which is then oxidized at the anode and liberates electrons and protons. Currently, the Hilbert fractal curve is used to design a DMFC [9]. It is a continuous space-filling curve that can be applied to grids of power. These curves are used to study the current collectors of the direct methanol fuel cell.

Currently, nanostructured catalysts lead to enhanced efficiency, robustness, and reliability in energy conversion and storage systems due to their unique physicochemical and electrochemical properties [10]. Carbon nanofiber webs have beenused as a porous methanol barrier to reduce the effect of methanol crossover in DMFCs. They have enhanced the cell performance at low methanol concentrations due to their balanced effect on reactant and product management, with an increase in peak power density compared to the conventional DMFC [11]. A mathematical model was developed and validated in realtime, to study the transient temperature distribution across a DMFC. The model was used to study temperature distribution across passive and active DMFCs as a function of process parameters and performance parameters [12].

Advanced control strategies based on methanol concentration were designed and implemented to increase the efficiency of a DMFC. This method is more efficient than performing several modifications to system layouts and operating strategies [13]. In addition to the control of an integrated fuel processing system, the fuel cell system in a DMFC provides efficient fuel cell operation and sustainable power source for various utilities [14].

In this study, a small DMFC with a 45 cm2sizedmembrane electrode assembly (MEA) was fabricated, and the effects of cell temperature on the DMFC performance were studied in simulation and in realtime.

#### **2. Materials and Methods**

#### *2.1. Model-Based Simulation*

Design of fuel cells does not necessarily mean higher efficiency. Improvement of fuel cell efficiency is much dependent on the modelling and control operations [15,16]. A Laplace domain model of a DMFC was developed using six different first-order equations which represented the DMFC operation (Table 1, Table 2). These equations were obtained from the material balances, potential balance against the anode and the cathode, coverage of species on the catalyst surface, and cell current. Aqueous methanol in the feed channel was considered to be in direct contact with the catalytic layer. Using the mathematical modelling equations of a DMFC, the MATLAB program was written with the s-function technique.


**Table 1.** Equation for the variation of CCH3OH and theanode/cathode overpotential.

**Table 2.** Equation for coverage of adsorbed species at the anode.


This program is used to compute the DMFC output voltage based on the derived anode overpotential, cathode overpotential, and the cell current as given below:

$$
\mathcal{U} = \mathcal{U}\_0 - \eta\_A + \eta\_C - \frac{d^M}{k^M} i\_{\text{Cell}} \tag{1}
$$

The output voltage *U* has been influenced by three variables, namely, the flow rate of methanol, the concentration of methanol, and the cell temperature. Among these variables, it was necessary to investigate which variable had a larger influence on the output voltage [17]. The developed model was simulated using MATLAB software (Figure 1). The influence of the above three variables over the DMFC output voltage was recorded and analyzed.

#### *2.2. Real-Time Experimentation*

In this research work, a single DMFC system was designed with a 45 cm2 active cross-sectional area (Figure 2). The cell was molded into a frame of a fiber-reinforced Teflon-coated jacket and was kept between two graphite blocks. Flow grooves for methanol and oxygen flow were provided in this system. The flow field comprised a series of 25 parallel channels with a 2 mm depth, a 1 mm wide rib and a 1 mm groove [18]. A provision for the tapping current and voltage measurement was provided. Electrical plate type heaters were placed behind each of the graphite blocks to heat the cell for achieving the desired operating temperature. Methanol solution was circulated with the peristaltic pump. Oxygen was supplied with 2 atm pressure. Catalysts of Pt-Ru/C and Pt/C were loaded at concentrations of 2 mg/cm<sup>2</sup> at anode and cathode, respectively. Nafion 117, a polymer electrolyte membrane, was used in MEA preparations. It has a high protonic conductivity, zero electronic conductivity, excellent mechanical stability, and low resistance. This polymer electrolyte membrane also acts as a separator between the anode and the cathode. A graphite plate is commonly used as reactant distributor cum current collector [19,20]. The reactions involved in the DMFC operation were as follows:

> *Anode reaction* : *CH*3*OH* + *<sup>H</sup>*2*<sup>O</sup>* <sup>→</sup> *CO*<sup>2</sup> + <sup>6</sup>*H*<sup>+</sup> + <sup>6</sup>*e*<sup>−</sup> Cathodereaction : 3/2*O*<sup>2</sup> + <sup>6</sup>*H*<sup>+</sup> + <sup>6</sup>*e*<sup>+</sup> <sup>→</sup> <sup>3</sup>*H*2O Overallreaction : C*H*3*OH* + 3/2*O*<sup>2</sup> → *CO*<sup>2</sup> + 2*H*2*O*

**Figure 1.** Real-time experimental setup of the direct-methanol fuel cell (DMFC) system.

**Figure 2.** Schematic diagram of the direct methanol fuel cell system.

#### **3. Results and Discussion**

#### *3.1. Model-Based Simulation*

A program written in MATLAB with the mathematical modelling equations of the DMFC was simulated as discussed in Section 2.1. In order to find out the most influencing variable of the DMFC, the analysis was carried out for cell temperature, methanol flow rate, and methanol concentration. During the simulation, the effect of different operating temperatures between 303 K and 343 K, different methanol flow rates in the ranges of 0.25 CCM to 2 CCM and different methanol concentrations in the ranges of 0.25 M to 2 M were recorded. The results obtained are plotted in Figures 3–5.

It was also observed from Figures 3 and 4 that the DMFC output voltage increased the methanol flow rate by up to 1 CCM and increased the 1M methanol concentration, thereafter, (i.e., greater than 1CCM methanol flow rate and higher than 1M methanol concentration) DMFC output voltage did not vary significantly and remained constant due to the methanol crossover to cathode. It was also observed from Figure 5 that the output voltage of the DMFC increased linearly with an increase in temperature up to 333 K and decreased beyond that due to the dry-out effect in the cell. From the dynamic studies, it is clear that the temperature is the most influencing variable on the output voltage of DMFC. Hence, this working temperature was selected as a manipulated variable for further investigation.

Step response analysis was carried out in simulation to validate the developed mathematical model and its behavior. Step input changes were given at the operating temperatures of 313 K, 323 K, and 333 K (25%, 50%, 75% of the cell temperature) with ± 10% and ± 15%. Step response of the system at 25% of the cell temperature for a step change of ± 10%, and a step change of ± 15% of cell temperature were observed and are given in Figure 6. Similar runswere conducted at 50% of the cell temperature and 75% of cell temperature and are recorded in Figures 7 and 8, respectively. From the step response curves, it was concluded that the developed model was behaving linearly in these operating temperatures.

**Figure 4.** Effect of methanol concentration.

**Figure 5.** Effect of cell temperature.

**Figure 6.** Step response of DMFC at 25% ofthe operating temperature.

**Figure 7.** Step response of DMFC at 50% of the operating temperature.

**Figure 8.** Step response of DMFC at 75% of the operating temperature.

#### *3.2. Real-Time Experimentation*

Experimentation of the squoval-shaped manifold holedesign(SSMHD)-based DMFC was carried out with a Pt-Ru/C catalyst combination in order to determine the optimum operating temperature (**TCell**) of the DMFC with a methanol concentration of 1 M and a methanol flow rate of 1 mL/min. The output voltage of DMFC for a defined current was measured using an electronic load bank and the performance of DMFC was analyzed. Voltage and power density obtained during the run of the experiment against current density were plotted (I–V Curve and I–P Curve) and are shown in Figures 9 and 10, respectively.

**Figure 9.** CCH3OM: Steady-state current density–voltage curve.

**Figure 10.** CCH3OM: Steady-state current density–power density curve.

It was observed from Figures 9 and 10 that an increase in the DMFC temperature from 303 K (30◦C) to 333 K (60◦C) led to an increase in DMFC performance, butfor a further increase in temperature beyond its boiling point 333 K (60◦C), DMFC performance was reduced due to the dry-out effect. Hence, it is reported that 333 K (60◦C) is the optimum operating temperature at which a maximum power density of 57.6 mW/cm<sup>2</sup> was obtained for a squoval-shaped manifold design (SSMD)-incorporated direct methanol fuel cell. In addition to the studies on temperature, the durability of a DMFC at different load currents, namely 2A, 4A, 6A, and 10A, was also studied.Output voltages were recorded and are shown in Figure 11. It was confirmed by the consistent, steady output that the SSMD-based DMFC was durable at the optimized operating conditions, namely a DMFC temperature of 333 K, a methanol concentration of 1 M, and a methanol flow rate of 1 mL/min.

**Figure 11.** Durability test on a squoval-shaped manifold hole design-based DMFC.

#### **4. Conclusions**

A direct methanol fuel cell (DMFC) is a particular type of fuel cell that has its own merits. In this work, a direct methanol fuel cell with a 45 cm<sup>2</sup> activation area in the MEA, and a Pt-Ru/C catalyst at the anode was fabricated and tested with the squoval-shaped manifold hole design-incorporated DMFC work station available in the laboratory. Performance of the DMFC was evaluated with different operating parameters, namely cell temperature, methanol flow rate, and methanol concentration in simulation. From the simulation results, it was concluded that cell temperature was the most influencing process variable. From the results, the optimum temperature (TCell) of the DMFC operation was identified as 333 K. Real-time experimentation was carried out with different cell temperatures, and the results were recorded. It was observed from the experimental data that the maximum power density of 57.6 mW/cm<sup>2</sup> at 288 mA/cm2 was achieved with the said operating temperature of 333 K. With this operating temperature, the durability of the DMFC was also verified at different load currents and found to be durable.

**Author Contributions:** This paper is come to the final shape with the conceptualization of R. Govindarasu and recourse of S. Somasundaram. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors are grateful to P. Kanthabhaba, Former Professor and Head, Department of Chemical Engineering, Annamalai University, Chidambaram-608002, India.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Air-Forced Flow in Proton Exchange Membrane Fuel Cells: Calculation of Fan-Induced Friction in Open-Cathode Conduits with Virtual Roughness**

**Dejan Brki´c 1,2,\*,**† **and Pavel Praks 1,\*,**†


Received: 8 April 2020; Accepted: 3 June 2020; Published: 11 June 2020

**Abstract:** Measurements of pressure drop during experiments with fan-induced air flow in the open-cathode proton exchange membrane fuel cells (PEMFCs) show that flow friction in its open-cathode side follows logarithmic law similar to Colebrook's model for flow through pipes. The stable symbolic regression model for both laminar and turbulent flow presented in this article correlates air flow and pressure drop as a function of the variable flow friction factor which further depends on the Reynolds number and the virtual roughness. To follow the measured data, virtual inner roughness related to the mesh of conduits of fuel cell used in the mentioned experiment is 0.03086, whereas for pipes, real physical roughness of their inner pipe surface goes practically from 0 to 0.05. Numerical experiments indicate that the novel approximation of the Wright-ω function reduced the computational time from half of a minute to fragments of a second. The relative error of the estimated friction flow factor is less than 0.5%.

**Keywords:** Colebrook equation; fuel cells; flow friction factor; open-cathode; pressure drop; symbolic regression; numerically stabile solution; roughness

#### **1. Introduction**

Flow friction is a complicated physical phenomenon and it is not constant but depends on flow rate and pressure drop. Because of its complexity, equations which describe the flow friction are mostly empirical [1]. The most used empirical formulation for turbulent pipe flow is given by Colebrook's equation [2]. Here we will evaluate flow friction factor caused by air flow in the cathodic side of the observed PEMFCsand we will provide accurate and consistent solution.

#### *1.1. Colebrook Equation for Pipe Flow Friction*

The standard Colebrook's friction equation for turbulent pipe flow [2]; Equation (1), follows logarithmic law and is based on an experiment performed by Colebrook and White in the 1930s [3], while its graphical interpretation was given by Moody in 1944 [4].

$$\frac{1}{\sqrt{\lambda\_{T(p)}}} = -2 \cdot \log\_{10} \left( \frac{2.51}{\mathcal{Re}\_{(p)}} \cdot \frac{1}{\sqrt{\lambda\_{T(p)}}} + \frac{\varepsilon\_{(p)}}{3.71} \right) \tag{1}$$

where:

λ*T*(*p*)—turbulent Darcy flow friction factor for pipes (dimensionless)


In Moody's diagram, for the turbulent regime the Reynolds number *Re*(*p*) goes from around 2320 to 108 while the relative roughness of inner pipe surface ε(*p*) from 0 for smooth surfaces to about 0.05 for very rough surfaces (on the other hand, flow friction for laminar regime for pipe flow λ*L*(*p*) does not depend of roughness ε(*p*) while the Reynolds number *Re*(*p*) goes from 0 to around 2320; the formula for laminar Darcy flow friction factor for pipe flow, λ*L*(*p*) = 64/*Re*(*p*) is not empirical but theoretically founded).

Colebrook and White experimented with airflow through pipes of different roughness of inner surface [3]. They used a set of pipes from which one left with a smooth inner surface, while others were covered with sand of different size of grains. For each pipe, one uniform grain size with glue as adhesive material was used. Thus, the pipes in the experiment were gradually smooth to the fully rough. The experiment revealed that the turbulent flow friction depends on the Reynolds number *Re*(*p*) and on the relative roughness of inner pipe surface ε(*p*). As can be seen from the Moody diagram [4], for the same values of the Reynolds number *Re*(*p*) the turbulent flow friction factor λ*T*(*p*) is higher in the pipes with higher relative roughness ε(*p*), where subsequently for the same flow, the corresponding pressure drop is higher, too.

The Colebrook equation is implicitly given in respect to the unknown turbulent flow friction factor λ*T*(*p*) and it can be rearranged in explicit form only in terms of Lambert W-function [5] or its cognate Wright ω-function, and even then further evaluation can be only approximated. Very accurate approximate formulas for the Colebrook equation for pipe flow based on the Wright ω-function are available [6,7].

#### *1.2. Modified Colebrook Equation for Flow Friction*

The Colebrook equation is empirical [3] and therefore possible modifications based on different conditions of flow can be done. We will evaluate here a modification for fuel cells [8,9]. Further, for example, US Bureau of Mines published a report in 1956 [10] that introduced a modified form of the Colebrook equation for gas flow where coefficient 2.51 was replaced with 2.825. Also, very recent modifications of a variety of empirical equations for pipe flow are available [11]. Here we analyze one modification for air flow friction in the open-cathode side of the observed type of PEMFCs [8], and subsequently, we give a stable and computationally efficient explicit solution which is valid in this case. Analogous analyses can be done for air flow for cooling of electronic products from hand-held devices to supercomputers [12,13]. Here we discuss the only influence of hydraulic effects of flow through the open-cathode conduits of fuel cells while available literature [14–21] should be consulted for thermodynamic aspects.

The Colebrook equation can be rearranged following data obtained from the experiment by Barreras et al. with fuel cells [8]. Based on this data and following analogy with pipe flow, it is estimated that virtual roughness is fixed by value ε(*FC*) = 0.03086. Therefore, the Colebrook equation can be rewritten for the observed fuel cell as Equation (2). We will further analyze this equation to provide its stable numerical solution. The value of virtual roughness ε(*FC*) = 0.03086 is from [9].

$$\frac{1}{\sqrt{\lambda\_{T(FC)}}} = -10 \cdot \log\_{10} \left( \frac{65.6}{\text{Re}\_{(FC)}} \cdot \frac{1}{\sqrt{\lambda\_{T(FC)}}} + \frac{\varepsilon\_{(FC)}}{0.1415596} \right) \tag{2}$$

where:

λ*T*(*FC*)—turbulent Darcy flow friction factor for fuel cells (dimensionless)

*Re*(*FC*)—Reynolds number (dimensionless)—the same definition as for pipes

ε(*FC*)—virtual relative roughness of fuel cell (dimensionless) *log*10—logarithmic function with base 10 *FC*—index related to Fuel Cells

#### **2. Proposed Model**

Proton exchange membrane fuel cells (PEMFCs) transform chemical energy from electrochemical reaction of hydrogen and oxygen to electrical energy [22–25]. Here we analyze fan-induced air-forced flow, based on data from the experiment with pressure drop in the cathode side of air-forced open-cathode PEMFCs by Barreras et al. [8]. Their experiment with fuel cells can be compared with the experiment performed by Colebrook and White with air flow through pipes [3]. Barreras et al. [8] use three different cathode configurations with aspect ratios *h*/*H* from 0.83 to 2.5 to form a mesh of cathodic channels to supply the fuel cell with enough air for cooling and with oxygen for a chemical reaction (roughness is real physical characteristic of pipe surface [26], but not of cathodic channels of fuel cells in terms of hydraulic performances).

Value of the Reynolds number *Re*(*FC*) during the experiment was from 45 to 4000, while as a difference from pipes, during flow through the observed fuel cell, the transition from laminar to turbulent flow occurred around *Re*(*FC*) = 500. In our numerical experiments, we use *Re*(*FC*) between 50 and 4100.

As already mentioned, the original Colebrook equation is valid for turbulent flow of air, water, or natural gas through pipes. On the other hand, for laminar flow, λ*L*(*p*) = 64/*Re*(*p*) is used whereas the transition from laminar to turbulent flow is around *Re*(*p*) ≈ 2320. This transitional border at the Moody's plot [4] is sharp where the equivalent sharp transition from laminar to turbulent flow for the observed fuel cell starts at around *Re*(*FC*) ≈ 500, as explained in [8]. Therefore, for airflow through the cathode side of the observed fuel cells, the flow friction factor λ*L*(*FC*) consists of two clearly defined types of flow:


For solving implicitly given equations, instead of iterative procedures [27,28], appropriate explicit approximations which are accurate but also computationally efficient can be used (review of appropriate explicit approximations for pipe flow friction is available in [29]). A computationally efficient and stabile unified equation for the observed type of fuel cells which is valid both for laminar and turbulent regime will be given in this article [30].

#### *2.1. Turbulent Flow*

In case of turbulent airflow during experiments with open-cathode PEMFCs, measurements show that pressure drop during turbulent flow at its cathode side follows logarithmic law, which form is comparable to the Colebrook's flow friction equation for pipe flow, but with different numerical values [8]. The flow friction related to air flow is given by Equation (3):

$$\frac{1}{\sqrt{\lambda\_{T(FC)}}} = -10 \cdot \log\_{10} \left( \frac{65.6}{\text{Re}\_{(FC)}} \cdot \frac{1}{\sqrt{\lambda\_{T(FC)}}} + 0.218 \right) \tag{3}$$

where:

λ*T*(*FC*)—turbulent Darcy flow friction factor for fuel cells (dimensionless)

*Re*(*FC*)—Reynolds number (dimensionless) – the same definition as for pipes

*log*10—logarithmic function with base 10

*FC*—index related to Fuel Cells

During turbulent flow, numerical values for the flow friction factor in pipe and fuel cells are different and that difference can go up to 60%. To make a direct connection between Equation (1) for pipe flow and Equation (3) for the observed fuel cell, Equation (4) can be used [25]:

$$\frac{1}{\sqrt{\lambda\_{T(FC)}}} = -14.17 + 5 \cdot \left( -2 \cdot \log\_{10} \left( \frac{2.51}{\text{Re}\_{(FC)}} \cdot \frac{1}{\sqrt{\lambda\_{T(FC)}}} + \frac{\varepsilon\_{(FC)}}{3.71} \right) \right) \tag{4}$$

where:

λ*T*(*FC*)—turbulent Darcy flow friction factor for fuel cells (dimensionless)

*Re*(*FC*)—Reynolds number (dimensionless)—the same definition as for pipes

ε(*FC*)—virtual relative roughness of fuel cell (dimensionless)

*log*10—logarithmic function with base 10

*FC*—index related to Fuel Cells

For Equation (4), virtual roughness can be recalculated based on the Colebrook equation as ε(*FC*) = 0.03086 for the observed fuel cell in the experiment [8]. This fuel cell was tested with three different cathode configurations [8]. As noted in [31], this roughness ε(*FC*) is not a real measurable physical characteristic of the surface of the used material for conduits (on the contrary for pipes ε(*p*) can be measured or at least estimated accurately [32–37]).

Both, Equations (3) and (4) are numerically unstable for *Re*(*FC*) < 575, which can be a critical problem knowing that turbulent zone starts for *Re*(*FC*) > 500. However, the novel solution proposed in this article is numerically stable.

Generally, implicitly given equations can be transformed in explicit form through the Lambert W-function [38,39]. The Lambert W-function [5] is defined as the multivalued function W that satisfies *<sup>z</sup>* = *eW*(*z*)·*W*(*z*). However, such transformation for the Colebrook equation for pipes contains a fast-growing term *e<sup>x</sup>* and because of that, overflow error in computers is possible [40,41]. Happily, results with fuel cells show that the solution is not in the zone where *e<sup>x</sup>* is too big to be stored in registers of computers. The model for fuel cells is given in Equation (5), while the related model for pipe flow friction model can be seen in [42].

$$\begin{array}{c} \frac{1}{\sqrt{\lambda\_{\Gamma(IC)}}} = a\_{\langle FC \rangle} \cdot W(\boldsymbol{\varepsilon^{\{IC\}}}) - \frac{R\_{\{IC\}}}{b\_{\langle FC \rangle}} \cdot \frac{\boldsymbol{\varepsilon\_{\{IC\}}}}{\boldsymbol{\varepsilon\_{\{IC\}}}} \\ \boldsymbol{x\_{\{FC\}}} = \ln \left( \frac{R\_{\{IC\}}}{a\_{\langle IC \rangle} \cdot b\_{\langle IC \rangle}} \right) + \frac{R \boldsymbol{\varepsilon\_{\{IC\}}}}{a\_{\langle IC \rangle} \cdot b\_{\langle IC \rangle}} \cdot \frac{\boldsymbol{\varepsilon\_{\{IC\}}}}{\boldsymbol{\varepsilon\_{\{IC\}}}} \\ a\_{\langle FC \rangle} = \frac{10}{\ln(10)} \\ b\_{\langle FC \rangle} = 65.6 \\ \frac{\boldsymbol{\varepsilon\_{\{FC\}}}}{\boldsymbol{\varepsilon\_{\{FC\}}}} = 0.218 \to \boldsymbol{\varepsilon\_{\{FC\}}} = 0.1415596 \end{array} \tag{5}$$

where:

λ*T*(*FC*)—turbulent Darcy flow friction factor for fuel cells (dimensionless) *Re*(*FC*)—Reynolds number (dimensionless)—the same definition as for pipes ε(*FC*)—virtual relative roughness of fuel cell (dimensionless) *a*(*FC*), *b*(*FC*), *c*(*FC*)—constants *x*(*FC*)—variable *log*10—logarithmic function with base 10 *ln*—natural logarithm *e*—exponential function

*W*—Lambert function

*FC*—index related to Fuel Cells

The parameter *c*(*FC*) for fuel cell is *c*(*FC*) = 0.1415596.

After procedures from [6,7,43], the following form for fuel cells expressed through the Lambert W-function and its cognate Wright ω-function is given in Equation (6):

$$\begin{array}{c|c} \frac{1}{\sqrt{\mathbf{r}\_{\rm{(FC)}}}} = \frac{10}{\ln(10)} \cdot \left[ \ln \left( \Delta\_{(\rm{FC)})} \right) + \mathcal{W} \left( \mathbf{c}^{\mathbf{x}/\rm{FC}} \right) - \mathbf{x}\_{\rm{(FC)}} \right] \\ \frac{1}{\sqrt{\mathbf{r}\_{\rm{(FC)}}}} = \frac{1}{a\_{(\rm{(FC)})}} \cdot \left[ B\_{(\rm{FC)}} + \alpha \left( \mathbf{x}\_{\rm{(FC)}} \right) - \mathbf{x}\_{\rm{(FC)}} \right] \\ \mathbf{x}\_{\rm{(FC)}} = B\_{(\rm{FC)}} + 0.218 \cdot \Delta\_{(\rm{FC)}} \\ B\_{(\rm{(FC)})} = \ln \left( \Delta\_{(\rm{(FC)})} \right) = \ln \left( \text{Re}\_{(\rm{(FC)})} \right) - 5.652138 \\ \Delta\_{(\rm{FC)}} = \frac{\mathcal{R}\_{(\rm{(FC)})} \cdot \mathbf{a}\_{(\rm{(FC)})}}{65.6} = \frac{\mathcal{R}\_{(\rm{(FC)})}}{284.9} \\ \frac{1}{a\_{(\rm{(FC)})}} = \frac{10}{\ln(10)} \approx 4.343 \end{array} \tag{6}$$

where:

λ*T*(*FC*)—turbulent Darcy flow friction factor for fuel cells (dimensionless) *Re*(*FC*)—Reynolds number (dimensionless)—the same definition as for pipes *a*(*FC*)—constant *x*(*FC*), Δ(*FC*), *B*(*FC*)—variable *ln*—natural logarithm *W*—Lambert function ω—Wright function *FC*—index related to Fuel Cells

However, symbolic regression applied on the explicit formulation, Equation (6), which involves *W*(*e <sup>x</sup>*(*p*)) <sup>−</sup> *<sup>x</sup>*(*p*) <sup>=</sup> <sup>ω</sup> - *x*(*p*) − *x*(*p*) gives very simple, but still accurate results in case of pipe flow [6,7] ([44,45] confirm these results), but unfortunately these analytical formulas, which are optimized for pipes, cannot be directly applied on the fuel cell equation. Fortunately, symbolic regression gives also very promising results for fuel cells as given in Equation (7):

$$\mathcal{W}(\boldsymbol{\varepsilon}^{\mathrm{(FC)}}) - \mathbf{x}\_{\mathrm{(FC)}} = \omega(\mathbf{x}\_{\mathrm{(FC)}}) - \mathbf{x}\_{\mathrm{(FC)}} \approx \frac{26.723}{\ln\left(\frac{R\varepsilon\_{\mathrm{(FC)}}}{284.9}\right) + 0.218 \cdot \frac{R\varepsilon\_{\mathrm{(FC)}}}{284.9} + 6.2611}} - 3.6795 \tag{7}$$

where:

*Re*(*FC*)—Reynolds number (dimensionless)—the same definition as for pipes *x*(*FC*)—variable *ln*—natural logarithm *W*—Lambert function ω—Wright function *FC*—index related to Fuel Cells

To avoid repetitive computations, parameters Δ(*FC*) and *B*(*FC*) are introduced, in Equation (8). Both symbolic regression analyses were performed in Eureqa, a commercial software tool, which automates the process of model building and interpretation [46,47].

$$\begin{array}{c} \frac{1}{\sqrt{\lambda\_{\text{(FC)}}}} = 4.343 \cdot \left( B\_{\text{(FC)}} + \frac{26.723}{B\_{\text{(FC)}} + 0.218 \cdot \Delta\_{\text{(FC)}} + 6.2611} - 3.6795 \right) \\ \qquad \qquad \qquad \Delta\_{\text{(FC)}} = \frac{R\_{\text{(FC)}}}{\text{2B4.9}} \\ B\_{\text{(FC)}} = \ln \text{(} \Delta\_{\text{(FC)}} \end{array} \tag{8}$$

where:

λ*T*(*FC*)—turbulent Darcy flow friction factor for fuel cells (dimensionless)

*Re*(*FC*)—Reynolds number (dimensionless)—the same definition as for pipes Δ(*FC*), *B*(*FC*)—variable

*ln*—natural logarithm *FC*—index related to Fuel Cells

#### *2.2. Unified Model*

Although the expression for laminar flow through pipes is λ*L*(*p*) = 64/*Re*(*p*), for fuel cells it is different, as given in Equation (9) [8]:

$$
\lambda\_{L(FC)} = \frac{58.91 + 50.66 \cdot e^{-\frac{3.44}{H}}}{Re\_{(FC)}} \tag{9}
$$

where:

λ*L*(*FC*)—laminar Darcy flow friction factor for fuel cells (dimensionless)

*h <sup>H</sup>*—channel depth/channel width used only in laminar flow (dimensionless)

*e*—exponential function

*FC*—index related to Fuel Cells

Values of *h*/*H* are from 0.83 to 2.5.

The experiment [8] shows that air flow through the cathode side of air-forced open-cathode PEMFCs are (1) laminar for the lower values of the Reynolds number, *Re*(*FC*) < 500 and (2) turbulent for the higher values, 500 < *Re*(*FC*) < 4000, where the Reynolds number is in hydraulics a very well-known dimensionless parameter that is used as a criterion for foreseeing flow patterns in a fluid's behavior (defined in the same way for air flow through pipes and here discussed air flow through fuel cells). The dimensionless Darcy's unified flow friction factor λ(*FC*), is the function of the switching function *y*, the laminar flow friction λ*L*(*FC*), and the turbulent flow friction λ*T*(*FC*). The unified coherent flow friction model that covers both laminar and turbulent zones is set by Equation (10) [30]:

$$
\lambda\_{\rm (FC)} = y \cdot \lambda\_{T(FC)} + (1 - y) \cdot \lambda\_{L(FC)} \tag{10}
$$

where:

λ(*FC*)—unified Darcy flow friction factor for fuel cells (dimensionless)

λ*T*(*FC*)—turbulent Darcy flow friction factor for fuel cells (dimensionless)

λ*L*(*FC*)—laminar Darcy flow friction factor for fuel cells (dimensionless)

*y*—switching function

*FC*—index related to Fuel Cells

The novel switching function *y* is given in Equation (11):

$$y = \frac{\mathcal{Re}\_{(FC)}}{\mathcal{Re}\_{(FC)} + \epsilon^{\text{558} - \mathcal{Re}\_{(FC)}}} \tag{11}$$

where:

*Re*(*FC*)—Reynolds number (dimensionless)—the same definition as for pipes

*y*—switching function

*e*—exponential function

*FC*—index related to Fuel Cells

The switching function was obtained by symbolic regression using HeuristicLab [47] and it is given in Figure 1.

**Figure 1.** Switching function given in Equation (11).

The laminar flow friction λ*L*(*FC*) depends on the Reynolds number *Re*(*FC*), but also on geometry of conduits, while the turbulent flow friction λ*T*(*FC*) depends only on the Reynolds number *Re*(*FC*). In the case of fuel cells, both coefficients are empirical. In addition, the switching function *y* contains the exponential function, (the similar situation is for calculation of λ*L*(*FC*) as already explained).

To avoid numerical instability, it is recommended to use the explicit approximation which gives the following unified formula in Equation (12).

$$\begin{aligned} \lambda\_{\langle FC\rangle} &= y \cdot \lambda\_{T\langle FC\rangle} + (1 - y) \cdot \lambda\_{L\langle FC\rangle} \\ \lambda\_{L\langle FC\rangle} &= \frac{58.91 + 50.66 \cdot e^{-\frac{3.41}{H}}}{Re\_{\langle FC\rangle}} \\ \frac{1}{\sqrt{\lambda\_{T\langle FC\rangle}}} &= 4.343 \cdot \left( B\_{\langle FC\rangle} + \frac{26.723}{\lambda\_{\langle FC\rangle} + 6.2611} - 3.6795 \right) \\ &\qquad \Delta\_{\langle FC\rangle} &= \frac{Re\_{\langle FC\rangle}}{284.9} \\ B\_{\langle FC\rangle} &= \ln\left(\Delta\_{\langle FC\rangle}\right) \\ x\_{\langle FC\rangle} &= B\_{\langle FC\rangle} + 0.218 \cdot \Delta\_{\langle FC\rangle} \\ y &= \frac{R\_{\langle FC\rangle}}{Re\_{\langle FC\rangle} + e^{\frac{358 - R\_{\langle FC\rangle}}{\lambda\_{\langle FC\rangle}}}} \end{aligned} \tag{12}$$

where:

λ(*FC*)—unified Darcy flow friction factor for fuel cells (dimensionless) λ*T*(*FC*)—turbulent Darcy flow friction factor for fuel cells (dimensionless) λ*L*(*FC*)—laminar Darcy flow friction factor for fuel cells (dimensionless) *Re*(*FC*)—Reynolds number (dimensionless)—the same definition as for pipes *h <sup>H</sup>*—channel depth/channel width used only in laminar flow (dimensionless) *x*(*FC*), Δ(*FC*), *B*(*FC*)—variables *y*—switching function *e*—exponential function *ln*—natural logarithm

*FC*—index related to Fuel Cells

For 216 = 65536 Sobol Quasi Monte-Carlo pairs [48], which cover *Re*FC = 50–4100 and for *h*/*H* from 0.83 to 2.45, the maximal relative error of the final calculated flow friction factor λFC using Equation (12) is 0.46% compared with the original Equation (2). The accuracy and speed of execution are tested through the code given in the next section.

#### **3. Software Code and Measurement of Execution Speed**

The unified equation for laminar and turbulent fan-induced air flow through open-cathode side of the observed PEMFCs is given in Equation (12), where the algorithm from Figure 2 should be followed.

**Figure 2.** Algorithm for solving Equation (12).

The turbulent flow friction <sup>λ</sup>*T*(*FC*) <sup>∼</sup> *LT* (with the intermediate step through 1/ λ*T*(*FC*) ∼ *LT*) can be expressed using "wrightOmega" function as *LT* = (*B* + *wrightOmega*(*x*) − *x*)/*a*, but it can be executed faster using approximation as given in Equation (8), as *LT* = 4.343 ∗ (*b* + 26.723./(*b* + 0.218 ∗ *a* + 6.2611) − 3.6795). The MATLAB code also works in GNU Octave, but it can be easily translated in any programming language. The final unified fuel cell model given by Equation (12) is coded in MATLAB as:

$$\begin{aligned} \textit{function } \textit{L = P} & \textit{P} & \textit{F (R, h)}\\ \textit{a = log(10) / 10;}\\ d &= R\* \textit{a / 65.6};\\ \textit{B = log(d);}\\ \textit{x = B + 0.218\*d;}\\ LT &= \textit{4.343\*(B + 26.723./(\textit{x} + 6.2611) - 3.6795);}\\ LT &= \textit{1./LT2;}\\ LL &= \textit{(58.91 + 50.66\* \exp(-3.4\*h))./R;}\\ y &= \textit{R./(R + \exp(558 - \text{R}));}\\ L &= y.\*LT + (1 - y).\*LL; \end{aligned}$$

⎤

⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥

in this code:

Output parameters of the function:


Input parameters of the function:



In MATLAB, symbol log() denotes the natural logarithm.

Implicitly given equations can be solved using iterative procedures [49,50], but also using appropriate accurate but also computationally efficient explicit approximations especially developed for the observed purpose. Our numerical results show that the computationally efficient approximation does not contain the time-consuming evaluation of the Wright ω-function [6,7], but simple polynomials found by symbolic regression [51], which can be easily evaluated on computers. It is because simple functions such as +, −, ∗ and / are executed directly in the CPU and hence they are very fast [43].

Using 216 = 65536 Sobol Quasi Monte-Carlo pairs [48], which cover *Re*(*FC*) = 50–4100 and for *h*/*H* from 0.83 to 2.45, the evaluation of the MATLAB built-in "wrightOmega" function required 30.9 s, while our novel approximation required only 0.0022 s (our approximation is around fourteen thousand times faster). Consequently, numerical tests show that the novel approximation presented here is very suitable for modeling of fan-induced flow friction in a mesh with virtual roughness for air-forced flow in the open-cathode PEMFCs, as it is very fast and still very accurate.

#### **4. Conclusions**

This paper gives a novel numerically stable explicit solution for flow friction during airflow in cathode side of open-cathode PEMFCs. Symbolic regression is successfully used for obtaining a cheap but still accurate approximation of the Wright-ω function for fuel cells-based explicit friction model and also for approximations of the switching function, which is needed for the unified formulation of fuel cell friction model. Numerical experiments indicate that the novel approximation of the Wright-ω function reduced the computational time from half of a minute to fragments of a second. The relative error of the estimated friction flow factor is less than 0.5%.

In future research, further analyses and experiments are foreseen to show if and how this value of the virtual roughness ε(*FC*) of fan-induced friction in a mesh of conduits for air-forced flow in the open-cathode PEMFCs can be changed and how it depends on fuel cells parameters (type of fuel cells, size, geometry of channels, etc.). The study should be extended to cover other types of fuel cell, to include water and heat management, etc. [52–60].

**Author Contributions:** D.B. had the initial idea for this article and he dealt with physical interpretation of the problem. P.P. conducted detailed numerical experiments (included symbolic regression) following initial inputs of D.B. D.B. wrote a draft of this article. Both authors contributed equally to this article. All authors have read and agreed to the published version of the manuscript.

**Funding:** The work of D.B. has been supported by the Ministry of Education, Youth and Sports of the Czech Republic through the National Programme of Sustainability (NPS II) project "IT4Innovations excellence in science-LQ1602", while the work of P.P. has been partially supported by the Technology Agency of the Czech Republic under the project "Energy System for Grids" TK02030039. Additionally, the work of D.B. has been partially supported by the Ministry of Education, Science, and Technological Development of the Republic of Serbia.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Notations**

The following symbols are used in this article:



#### **References**


© 2020 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 (http://creativecommons.org/licenses/by/4.0/).

## *Article* **Environmental Sustainability Assessment of Typical Cathode Materials of Lithium-Ion Battery Based on Three LCA Approaches**

#### **Lei Wang 1, Haohui Wu 1, Yuchen Hu 1, Yajuan Yu 1,\* and Kai Huang 2,\***


Received: 2 January 2019; Accepted: 24 January 2019; Published: 7 February 2019

**Abstract:** With the rapid increase in production of lithium-ion batteries (LIBs) and environmental issues arising around the world, cathode materials, as the key component of all LIBs, especially need to be environmentally sustainable. However, a variety of life cycle assessment (LCA) methods increase the difficulty of environmental sustainability assessment. Three authoritative LCAs, IMPACT 2002+, Eco-indicator 99(EI-99), and ReCiPe, are used to assess three traditional marketization cathode materials, compared with a new cathode model, FeF3(H2O)3/C. They all show that four cathode models are ranked by a descending sequence of environmental sustainable potential: FeF3(H2O)3/C, LiFe0.98Mn0.02PO4/C, LiFePO4/C, and LiCoO2/C in total values. Human health is a common issue regarding these four cathode materials. Lithium is the main contributor to the environmental impact of the latter three cathode materials. At the midpoint level in different LCAs, the toxicity and land issues for LiCoO2/C, the non-renewable resource consumption for LiFePO4/C, the metal resource consumption for LiFe0.98Mn0.02PO4/C, and the mineral refinement for FeF3(H2O)3/C show relatively low environmental sustainability. Three LCAs have little influence on total endpoint and element contribution values. However, at the midpoint level, the indicator with the lowest environmental sustainability for the same cathode materials is different in different methodologies.

**Keywords:** LIBs; environmental sustainability; cathode material; LCA

#### **1. Introduction**

With the expansion of the LIBs market, new cathode materials are constantly being developed [1]. In terms of weight fraction and cost, the cathode part for LIB is the most significant sector [2]. However, long-standing effort has been devoted to the development of high energy density and capability cathode materials [3], meeting the demand of electric vehicles, power tools, and large electric power storage units [4]. In fact, with the increase in energy density and capacity, many trace LIBs have an increasing impact on the environment [5]. Meanwhile, modern society must overcome many difficulties, such as obtaining natural resources and protecting the natural environment [4]. Before we commercialize a new cathode model, its environmental cost should be considered. Andersson confirmed the feasibility of environmental sustainability assessment in LCA for product development [6]. Numerous studies have quantified the impact of different types of LIB on the production process in its lifecycle [7]. Slowly, LCA is becoming more commonly used as a standardized method to determine the impact of a product or service on the environment throughout its whole life [8].

Although many new cathode models are being researched, the introduction of LiCoO2/C has enabled the commercialization of the first LIB [2]. LIBs have been available on the market from Sony Corp. since the early 1990s, and LiCoO2/C has become the leading LIB system [9]. LiFePO4/C stands as a competitive candidate cathode material for the next generation of a green and sustainable LIB system due to its long life span, abundant resources, low toxicity, and high thermal stability [10]. Meanwhile, the LiFe0.98Mn0.02PO4/C, as an improved cathode material for LiFePO4/C [11], also shows promising development potential. As Zeng [12] confirmed, the electrochemical performance of LiFePO4/C was remarkably improved by a slight manganese substitution, creating the general formula LiFeXMn1-XPO4/C [13]. However, most commercial cathode materials are LiCoO2/C, whose actual capacity is 140–155 mAhg<sup>−</sup>1, while LiFePO4/C or LiFe0.98Mn0.02PO4/C's theoretical capacity is only 170 mAhg<sup>−</sup>1, so LIBs urgently need a new cathode model, like FeF3(H2O)3/C [14].

We can see that electrochemical properties [15,16] have been the driving force for the development of different cathode materials. Many new studies have focused on electrochemical properties. However, due to the increasing prevalence of different cathode materials in electronic equipment and vehicles, their impact on the environment also needs to be considered [17]. Peters summarized the reviewed studies, focusing on energy demand and warming gas emissions for LiCoO2/C and LiFePO4/C [7]. Wang [18] evaluated the LIB with lithium-rich materials used in an electric vehicle throughout the life cycle of the battery. Wang [19] measured the carbon footprints of three industry lithium-ion secondary battery chains, and came to the conclusion that electric energy consumption is the main factor of lithium ion battery production companies in generating a carbon footprint. Liang [20] directly adopted LCA to assess the greenhouse gas emissions of LIBs. Deng [21] used LCA to assess high capacity molybdenum disulfide LIB for electric vehicles. Zackrission [22] elevated lithium-air batteries by LCA to quantify its climate impact, abiotic resource depletion and toxicity. Gong [23] evaluated four cathode materials by family footprint including carbon footprint and water footprint, and both direct and indirect water footprint [24].

In this paper, we not only pay attention to the total environmental sustainability of these four cathode materials, as previous studies have done, but also find that different indicators damage their environmental sustainability in different ways. Indeed, except for these common problems, certain indicators play different roles in specific cathode materials. In addition, from the results, we find that the different LCAs we choose have impacts on our evaluation. Based on the analysis of the calculation in three kinds of LCA, some suggestions are put forward for different cathodes, and the influence of different LCA is studied. There are three LCAs we chose for our study. Firstly, IMPACT 2002+ life cycle impact assessment methodology proposes a feasible implementation of a combined midpoint/damage approach, linking all types of life cycle inventory results (elementary flows and other interventions) via 14 midpoint categories to four damage categories, especially to human toxicity and ecotoxicity [25]. Secondly Eco-indicator 99(EI-99) is also used, a damage-oriented method to assess the emissions, extractions, and land use in all processes, and the damage to human health, ecosystem quality, and resources is calculated [26] ReCiPe assesses 18 impact categories at midpoint level (ozone depletion, human toxicity, etc.), and three endpoint categories (human health, ecosystems, resources) at endpoint level [27]. Based on these three LCAs, we divide the results into four innovative parts: (1) according to the total value of four different cathode materials, we have a preliminary understanding of the environmental sustainability ranking of these four cathode materials; (2) according to the endpoint value, we find one common environmental problem of the four cathode materials; (3) from the midpoint value, indicators showing the best or worst environmental sustainability of four cathode materials are summarized; (4) from the perspective of element contribution, the key elements that have obvious influence on the environmental sustainability of these four cathode materials are found. In fact, LiCoO2/C does not always show the lowest environmental sustainability among all indicators. The new cathode model FeF3(H2O)3/C also shows the lowest environmental sustainability in some respects. Although the chemical composition of LiFe0.98Mn0.02PO4/C is similar to that of LiFePO4/C, the metal resource consumption for the former is far larger. Instead of focusing on the overall environmental sustainability, we care more about these specific indicators in different LCAs.

#### **2. Materials and Methods**

#### *2.1. LCAs and Cathode Materials*

At present, IMPACT 2002+ and EI-99 have been widely recognized in many studies used to quantify the environmental impact of products [28]. Moreover, ReCiPe, as an improved LCA for EI-99, can evaluate more detailed indicators [29]. These three LCA methods are selected as the research methods. Although the results of environmental impacts are similar for equivalent categories, the ReCiPe and IMPACT 2002+ methods provide more categories for evaluation and comparison than EI-99 [30]. Furthermore, the assessment categories among these three LCAs are sufficient for us to make a comprehensive comparison. For example, the water issue is only evaluated in IMPACT 2002+, and acidification issues are only considered in both IMPACT 2002+ and EI-99. Indeed, there are many differences among these three LCAs, like the characterization of the unit and categories [31]. In addition, three cathode materials of LIB industry, LiCoO2/C, LiFePO4/C, and LiFe0.98Mn0.02PO4/C, are selected as research objects. These three cathode models are marketed and mass produced every year. The evaluation of these major cathode materials is of great significance to the environmental sustainability development of LIBs. Furthermore, a new cathode model, FeF3(H2O)3/C, is also evaluated as a comparison. The discovery of its environmental advantages can provide some enlightenment for these three market-oriented cathode models, which will help us to manufacture environmentally friendly cathode materials. Finally, the impact of different LCAs on sustainability assessment is also considered.

#### *2.2. Experimental Designation*

#### 2.2.1. Scope and Function Unit

To compare the four cathode materials on the same basis, we should stipulate the scope and function unit of these four cathode materials. We only research the environmental impact of the cathode part. The original quality list comes from the existing literature and laboratory. After a normalized conversion, we make the quality list meet the functional unit (1 kg), which means that the mass of each cathode model we evaluated is equal.

#### 2.2.2. Experimental Devices

This study is conducted by global LCA software Simapro released in 1990, and more than 80 countries have recognized its authority. Simapro allows researchers to collect, analyzes, and monitor the sustainability performance of products and services, from extraction of raw materials to manufacturing, distribution, use, and disposal [32]. In general, this includes (a) a user interface for modeling the product system, (b) a life cycle unit process database, (c) an impact assessment database with data supporting several life cycle impact assessment methodologies, and (d) a calculator that combines numbers from the databases in accordance with the modeling of the product system in the user interface [33].

#### 2.2.3. Experimental Process

We researched these four cathodes at a deeper level. Of course, the total value of environmental sustainability for these four cathodes is shown in three LCAs. However, we also looked at endpoint values and midpoint values. In short, the lower the value is, the better the environmental sustainability. Finally, the element contribution proportion to endpoint and total values is calculated.

(1) Simulated Assembly

To meet the requirements of functional unit (1 kg), the original data were converted into a standard quality list. Then we assembled four complete cathode models in Simapro.

#### (2) Calculation

To compare these four cathode materials, we have processed on the raw calculated data included in Supplementary Materials (Tables S1–S3), like normalization and contribution analysis. In fact, the data in the Supplementary Materials is calculated firstly by Simapro and all values has their own units which undoubtedly make it more difficult to compare these cathode materials. The detailed calculations of these four cathodes are listed as follows:


Some indicators are simplified: the total value in three LCAs, Ti(i = 1, 2, 3); the endpoint indicators in different LCAs(IMPACT 2002+, Xm(m = 1, 2, 3, 4); EI-99, Yn(n = 1, 2, 3); ReCiPe, Zh(h = 1, 2, 3)). The research process is shown in Figure 1.

**Figure 1.** The research process.

#### **3. Results and Discussion**

#### *3.1. The Total Value of the Four Cathodes' Environmental Impact*

It is necessary to qualify the total environmental impact of these four cathode materials before a deeper investigation, because all midpoint and endpoint indicators are divided from the total environmental impact. Unlike Gong's research on evaluating environmental, economic, and electrochemical performance indicators by footprint family and Peters's summary focusing on energy demand and warming gas emissions for LiCoO2/C and LiFePO4/C, we make a pure and further environmental assessment for these four cathode materials, and do not consider the economic benefits, electrochemical properties, energy demand, etc. This study aims to make a comprehensive assessment of these four typical cathodes directly. The total values for these four cathode materials in three LCAs are shown in Figure 2.

**Figure 2.** Total values for four cathodes in three LCAs: (**A**) LiCoO2/C, (**B**) LiFePO4/C, (**C**) LiFe0.98Mn0.02PO4/C, (**D**) FeF3(H2O)3/C.

As Figure 2 shows, we sort the four cathode types in descending order of environmental sustainability potential: D, C, B, and A. Although different LCAs have their own calculation standards, the trend of these four cathodes is the same. As Gong [23] confirms, D has the best environmental performance compared with B and C. However, the study does not explain the differences in specific indicators between B, C, and D. More detailed indicators that are important for each of the three cathodes are found in this study. In addition, material A is added to the study as another major market cathode model. No matter what LCA we choose, the new cathode model D always presents the best potential for environmental sustainability, while material A performs the worst. The environmental sustainability of material B is close to that of material C. Methodological emphasis on environmental assessment is only reflected in quantitative values, rather than the qualitative environmental sustainable potential among these four cathodes.

#### *3.2. Endpoint Level*

To distinguish concrete environmental impacts of different cathode materials, we calculate all endpoint indicators in Table 1. In any LCA, material A always has the highest endpoint value among these four materials. D is the smallest. Except the resource consumption in ReCiPe, 0.449 Pt of C larger than 0.393 Pt of B, other values for material B are slight larger than those for material C. IMPACT 2002+ and EI-99 have less impact on these four cathode materials' ranking of environmental sustainable potential. The resource consumption value between B and C calculated by ReCiPe is different to that calculated by the other two LCAs.


**Table 1.** Endpoint values for four cathodes.

Note: refer all endpoints (X1-4; Y1-2; Z1-3) to Figure 1.

The total value consists of the values of three or four endpoint indicators. To make a more intuitive observation, we calculated the contribution rates of different endpoint indexes to the total value, as shown in Figure 3. In three LCAs, the maximum contribution proportion of these four cathodes comes from human health. For material A, the contribution rate of ecosystem quality is relatively large in IMPACT 2002+.

**Figure 3.** Endpoint indicators' contribution to the total value.

Obviously, in three LCAs, the environmental sustainability capacity for these four cathodes is related to human health, which means human health is a common problem for these four cathode materials to solve. In fact, Wanger [34] has confirmed that the effect of LIBs on human health is a common problem for LIBs. For the cathode, the effect on human health remains a major concern. Its existence may require a major technical improvement to overcome. For these individual problems in different cathodes, we can learn from the strengths and weaknesses of different cathodes. For example, material A always shows the largest environmental load in these four cathodes. Reducing its yield or finding alternative models, elements, or mechanisms is a feasible way to reduce its impact on ecosystem quality. As we know, IMPACT 2002+ and EI-99 have the same ranking for the environmental sustainability among these four cathode materials. Three LCAs show that the impact on human health is a common problem for these four cathode materials. However, in ReCiPe, material C consumes more resources compared with B. In IMPACT 2002+, material A's impact on ecosystem quality makes a relative contribution to its total environmental impact.

#### *3.3. Midpoint Level*

Similarly, each endpoint indicator can be divided into a number of midpoint indicators. In order to avoid interference from different units and magnitude, we choose all values of material B as the benchmark and normalize all values of the other three cathode materials. These indicators with an extreme value always show great disadvantages and advantages for different cathodes, and these normalized values, obviously larger or less than 1, are more meaningful for individual cathode improvement. In IMPACT 2002+, the difference values between all normalized values and B' normalized value (1) are shown in Figure 4. In order to consider the impact of different LCA, we still give the unit of each midpoint indicator, reflecting their evaluation criteria. As we can see, there are 14 midpoint indicators, each of which has a different unit of measurement. To some extent, IMPACT 2002+ is more suitable for characterization evaluation. For example, we can use unit kg PO4 p-lim to express the land use problem. Moreover, because of the presence of phosphorus, the data are meaningful for eutrophication.

**Figure 4.** The normalized midpoint value in IMPACT 2002+.

As Figure 4 shows, when we regarded all normalized values of B as the baseline, except for the value of ionized radiation, other normalized values of A are larger than those of B. In particular, the eco-toxicity values of water body and surface are much larger, about 10.722 and 12.824, respectively. The high toxicity of cobalt may be the main reason for this situation. In fact, water problems and surface problems are difficult to completely separate, e.g., the toxicity of surface water. For water toxicity of A, water footprint assessment [35] may be a great method to quantify its water problems and reflect its toxicity from another perspective. For material D, except for the value for mineral refinement, other values are less than for B. Though D shows the best environmental sustainability, a green mineral refinement process is needed for material D. Finally, all normalized values for material C are slightly less than those for B. Material C, as an improved cathode to B by slight manganese substitution, has similar environmental sustainability potentiality to material B. The close element composition between material B and C, as the common formula shows, LiFeXMn1-XPO4/C (x = 0.98 in this study), accounts for the similar environmental sustainability. IMPACT 2002+ shows a sensitive assessment for cathode A, especially on water and surface ecotoxicity.

In EI-99, normalized values are shown in Figure 5. These indicators have the same cells separated from the same endpoint indicators. Compared with the IMPACT 2002+, these indicators cannot reflect specific substances due to their common units.

**Figure 5.** The normalized midpoint value in EI-99.

As Figure 5a–c shows, except for the radiation value, the normalized value of A is greater than that of B, and the respiratory inorganic matter, land occupation, and mineral resource problems of A are obvious, at 2.286, 2.247, and 1.982, respectively. For D, its mineral resource value is far larger, about 0.938. The ecological toxicity value is slightly larger. In addition to the mineral problems noted in IMPACT 2002+, the ecological toxicity of material D in EI-99 also becomes a low environmental sustainability index. Finally, all values of C are close to 1, as IMPACT 2002+ shows.

In ReCiPe, all normalized values are shown in Figure 6. The midpoint indicators in ReCiPe are more detailed than EI-99. The unit for these indictors is different.

**Figure 6.** The normalized midpoint value in ReCiPe.

As Figure 6a–c shows, the normalized values of nature land transformation and urban land occupancy for material A are far larger, at 2.975 and, 2.014 respectively. Actually, material D performs very well on these two indicators, with −0.586 and −0.702, respectively. Land occupation is an important part of the assessment of ecosystem quality [36]. That is a key difference between material A and D. For material C, not all normalized values are close to B, especially metal resource consumption, which is as high as 1.752. The method ReCiPe concentrates more on the economic costs (\$) in resource consumption [37], which means the economic cost in the slight manganese substitution process for material B needs to be cut down. To reduce the consumption of metal resources in the manganese substitution process must be a key issue for the development of material C. Finally, material D had low environmental sustainability in human toxicity, freshwater toxicity, and metal consumption compared with material B.

In addition, the problem of resource consumption deserves our attention. This endpoint in three LCAs is divided into two common midpoints, renewable and non-renewable resource consumption. We calculated the ratio of non-renewable resources to renewable ones. Mineral refinement in IMPACT 2002+, the mineral resource in EI-99 and the metal resource in ReCiPe are regarded as the renewable resource consumption. Likewise, the non-renewable energy, the fossil fuels and the fuel exhaustion are divided into the non-renewable resource consumption. The ratio is the non-renewable resources consumption per unit renewable resources consumption.

As Figures 4b, 5d and 6d show, material B always has the highest ratio in three LCAs, 626.614 (MJ primary/MJ surplus) in IMPACT 2002+, 41.359 (MJ/MJ) in EI-99, 8.780 (\$/\$) in ReCiPe. That means every unit renewable resource consumption needs more non-renewable resource in the whole life cycle of material B. Material B has the lowest environmental sustainability among these four cathodes. More green processes with low non-renewable resource consumption are needed for material B. As worldwide concern about fossil fuels grows, efforts at non-renewable resource protection are urgently required [38]. These high ratios need to be reduced. Integrating the renewable resources in a

small isolated power system, an isolated and complete battery [39], and improving the capacity for cathodes [40] are promising directions to achieve this goal.

Among the four cathode materials, the emphasis in the three LCAs is different. For material A, three LCAs all think that iron radiation is not a serious issue. The main problem in IMPACT 2002+ is ecotoxicity. On the contrary, EI-99 and ReCiPe think that the land issue is a serious issue. For materials B and C, the values are mostly close to each other except for the metal resource consumption. For material D, three LCAs all show its low environmental sustainability in terms of mineral resource consumption, and its toxicity is noted in EI-99 and ReCiPe.

#### *3.4. The Element Contribution to Environmental Sustainability*

In this part, we use the elemental symbol to represent all elements in tables as follows: lithium (Li), cobalt (Co), manganese (Mn), iron (Fe), fluorine (F), nitrogen (N), phosphorus (P), and oxygen (O). The contribution of the element to the endpoint and the total value is calculated.

#### 3.4.1. The Element Contribution in Endpoint Level

In the three LCAs, the element contribution rate of material A is shown in Table 2. Cobalt is the largest contributor to ecosystem quality, followed by lithium. The other endpoint indexes were mainly affected by lithium, followed by cobalt. So cobalt impairs the environmental sustainability of material A in terms of ecosystem quality (endpoint) and land and toxicity issues (midpoint).


**Table 2.** The element contribution proportion for material A.

For material B, the contribution rate of elements is shown in Table 3. The highest contribution proportion of each endpoint value always comes from lithium (the largest value is Y1 (83.70%), and the smallest value is Y2 (60.50%)), followed by phosphorus. Oxygen, iron, and nitrogen contribute little to each endpoint value.

**Table 3.** The element contribution proportion for material B.


Table 4 shows the element contribution rate of material C. The highest contribution rate of element to each endpoint value is still lithium, followed by phosphorus. The rest of the elements have less weight.


**Table 4.** The element contribution proportion for material C.

Finally, the element contribution proportion for material D is showed in Table 5. As a non-lithium composition model (no lithium in quality list), fluorine has the highest contribution rate among all mid-point indicators, at almost 90.00%.


**Table 5.** The element contribution proportion for material D.

Some issues can be found by the element contribution. The ecosystem quality in material A is mainly affected by cobalt. To seek a substitute material or reduce the quality of cobalt in production process will be a method to improve its environmental sustainability. As one heavy metal element, cobalt has higher ecosystem toxicity and pollution capacity than other elements in the four cathode materials. This is why many efforts to recover A are concentrated not only on lithium but also on cobalt [41]. There are two examples of Co substitutes. Xiang [42] improved the electrochemical kinetics of lithium manganese phosphate via Co-substitution. Di Lecce [43] investigated a new Sn-C/LiFe0.1Co0.9PO4 full lithium-ion cell with ionic liquid-based electrolyte. The two studies demonstrated the feasibility of producing Co-substitutes and for improving the environmental sustainability for A.

Except the ecosystem quality for material A mainly affected by cobalt, others indicators' environmental sustainability for materials A, B and C is mainly affected by lithium. However, the phosphorus in material B and C also has a great impact on their environmental sustainability. Reducing the consumption of lithium is an ongoing aim in the development for LIBs. More research on cobalt and phosphorus is also needed. All three LCAs show the same result—that is, different LCAs have no influence on the determination of the main elements contributing to the endpoint index values.

#### 3.4.2. The Element Contribution to the Total Value

The element proportions in the total level are shown in Figure 7.

**Figure 7.** The element contribution proportion to the total values in three LCAs.

In all LCAs, the largest contribution to the total value in materials A, B and C comes from lithium, up to 80.00% in EI-99 for materials B and C. For material D, it comes from the element fluorine, at around 90.00% in three LCAs. Moreover, the second contribution element for material A is cobalt, at 38.78% in EI-99 and 48.67% in IMPACT 2002+. For materials B and C it is phosphorus, at about 20.00%. Other elements made just a small contribution to the total values, less than 5.00%. This result is consistent with Yang's [44] research that there is significant waste of valuable metallic resources in LIBs and the environmental load of lithium consumption is the largest among all elements. Furthermore, except for the contribution of lithium, phosphorus plays an important role in the environmental sustainability potential of materials B and C. For material A, the influence of cobalt also cannot be ignored.

#### *3.5. The Methodologies' Influence to Environmental Sustainability Assessment*

Based on the discussion at the above four levels, we find that these four cathodes have different environmental sustainability potential due to the different LCA we use.


#### **4. Conclusions**

Different LCAs show different quantitative results in these four cathode materials. Qualitative assessments of these three LCAs is similar, both in terms of contributing elements and the total value. At the endpoint level, except for the resource consumption for LiFe0.98Mn0.02PO4/C and LiFePO4/C in ReCiPe, the ranking of other indicators' values is consistent with the total values. Four cathode models are ranked in descending order of environmental sustainability potential: FeF3(H2O)3/C, LiFe0.98Mn0.02PO4/C, LiFePO4/C, and LiCoO2/C. At the midpoint level, most indicators are consistent with the ranking. However, the most serious problem is determined differently based on different methodologies.

**Supplementary Materials:** The following are available online at www.mdpi.com/2227-9717/7/2/83/s1, Table S1: Calculation for mid point indicators in IMPACT 2002+, Table S2: Calculation for midpoint indicators in EI-99, Table S3: Calculation for midpoint indicators in ReCiPe.

**Author Contributions:** Conceptualization, L.W. and H.W.; Methodology, L.W. and H.W.; Software, Y.Y.; Validation, Y.Y. and K.H.; Formal Analysis, H.W. and Y.Y.; Investigation, L.W. and H.W.; Resources, Y.Y.; Writing—Original Draft Preparation, L.W.; Writing—Review and Editing, L.W. and H.W.; Visualization, L.W. and H.W.; Supervision, Y.Y. and K.H.; Project Administration, Y.Y. and K.H.; Funding Acquisition, Y.Y. and K.H.

**Funding:** This research was funded by: The Fundamental Research Funds for the Central Universities of China (No. 2018BLCB-05), the National Natural Science Foundation of China (No. 51474033), and Beijing Natural Science Foundation (9172012).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 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 (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Processes* Editorial Office E-mail: processes@mdpi.com www.mdpi.com/journal/processes

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18