Sintering process simulation of a solid oxide fuel cell anode and its predicted thermophysical properties

Accepted Manuscript Sintering Process Simulation of a Solid Oxide Fuel Cell Anode and Its Predicted Thermophysical Properties Pei Fu, Min Yan, Min Zeng, Qiuwang Wang PII: DOI: Reference:

S1359-4311(17)33975-3 http://dx.doi.org/10.1016/j.applthermaleng.2017.06.061 ATE 10580

To appear in:

Applied Thermal Engineering

Received Date: Revised Date: Accepted Date:

21 February 2016 7 April 2017 12 June 2017

Please cite this article as: P. Fu, M. Yan, M. Zeng, Q. Wang, Sintering Process Simulation of a Solid Oxide Fuel Cell Anode and Its Predicted Thermophysical Properties, Applied Thermal Engineering (2017), doi: http:// dx.doi.org/10.1016/j.applthermaleng.2017.06.061

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Revised Manuscript for Applied Thermal Engineering (ATE-2016-11703R1) Presented at the 5th Asian Symposium on Computational Heat Transfer and Fluid Flow (ASCHT2015), 2015.11.22-11.25, Busan, Republic of Korea. (Original paper title: “Simulation of the sintering process of porous anode for solid oxide fuel cells” and Paper number ASCHT15-Mon02-004)

Sintering Process Simulation of a Solid Oxide Fuel Cell Anode and Its Predicted Thermophysical Properties Pei Fu 1, Min Yan 2, Min Zeng 1 and Qiuwang Wang 1,*

1 Key Laboratory of Thermo-Fluid Science and Engineering, MOE, School of Energy and Power Engineering, Xi'an Jiaotong University, Xi'an, Shaanxi, 710049, P.R. China 2 School of Energy Engineering, Western Mine Exploitation and Hazard Prevention of the Ministry of Education, Xi’an University of Science & Technology, Xi’an, Shaanxi, 710054, P.R. China *Corresponding author: Tel and Fax: +86(29)82665539; E−mail: [email protected]

Highlights 1. A CG-MD method is improved to simulate the sintering process of SOFC anode. 2. Nanostructure and relevant thermal properties of the sintered anode are obtained. 3. Effects of sintering conditions and composition are systematically investigated. 4. Analyses and predictions provide a guide to obtain the desired properties.

1

Abstract Complex porous structure and properties of a solid oxide fuel cell (SOFC) anode significantly affect the SOFC performance. In this study, the sintering process of the typical NiO-YSZ anode is simulated by a coarse-grained molecular dynamics method. After the sintering process simulation, the porous nanostructure is obtained and the thermophysical properties of the sintered anode are evaluated. It is found that the volume shrinkage occurs along with the self-assembly phenomenon during the sintering process. Sintering at relatively low temperature and high pressure could contribute to the densification of the anode. Compared with the experimental structure, the sintered anode obtained from this study shows a common morphology. Sintering at the nanoscale is beneficial to distribute the functional components equally. Furthermore, the effects of the sintering temperature, sintering pressure and material composition on the volumetric heat capacity and thermal expansion coefficient of the sintered anode are systematically discussed. The predicted thermophysical properties of the anode agree well with the open literature data. Results in this study could provide a guide of the sintering conditions and composition during the experimental studies to obtain the desired properties of the SOFC anode.

Keywords: coarse-grained molecular dynamics (CG-MD); SOFC; anode; NiO-YSZ; sintering process; thermophysical properties

2

Nomenclature V

Interaction potential, kJ·mol-1; Volume, cm3

r

Distance between two particles, nm

Kb

Bond force constant, kJ·mol-1·nm-2

Kθ

Angle Force constant, kJ·mol-1·rad -2

b0

Equilibrium value of chemical bond, nm

q

Quantity of electric charge, C

F

Mechanical Force, kJ·mol-1·nm-1

m

Atomic mass, u

ri

Position of the particle i, nm

t

Time, ps

CV

Volumetric heat capacity, J·cm-3·K-1

E

Internal energy, J

T

Temperature, K

M

Molar mass, g·mol-1

R

Gas constant, 8.314 J·K-1·mol-1

L

Length, m

p

Pressure, bar

Greek symbols θ

Angle among three particles, deg

θ0

Equilibrium value of angle, deg

φ

Dihedral among four particles, deg

ε

The depth of the potential at the minimum, kJ·mol-1

σ

The distance where the potential is zero, nm

π

Constant

ε0

Vacuum dielectric constant, 8.85e-12 F·m-1

εr

Relative dielectric constant 3

ρ

Density, g·cm-3

αL

Linear thermal expansion coefficient, K-1

Subscripts i, j, k, l Particle index Abbreviation AA

All-atom

ASE

Atomistic Simulation Environment

CG

Coarse-grained

DFT

Density functional theory

FCC

Face-centered cubic

KMC

Kinetic Monte Carlo

LJ

Lenard-Jones

MD

Molecular dynamics

NPT

Isothermal-isobaric ensemble

NVT

Canonical ensemble

SBCG Shaped-based coarse-graining SC

Simple cubic

SHC

Specific heat capacity

SOFC

Solid oxide fuel cell

TEC

Thermal expansion coefficient

TEM

Transmission Electron Microscopy

TPB

Three-phase boundary

VHC

Volumetric heat capacity

VMD

Visual Molecular Dynamics

YSZ

Yttria-stabilized zirconia

8YSZ

8 mol% yttria stabilized zirconia

4

1. Introduction Compared with conventional power generation systems, a solid oxide fuel cell (SOFC) can effectively solve the problems, such as the low utilization rate of energy and environmental pollution, due to its high efficiency, low emission and extensive fuel selection [1, 2]. As shown in Fig. 1, a SOFC mainly contains three essential components: an anode, electrolyte and cathode, which are mostly constructed from ceramic, metal or metal ceramic composites [3]. It is well known that the electrochemical reactions only occur at the three-phase boundary (TPB) in the anode [4]. Exothermically electrochemical reactions in the anode cause extremely high local temperature, which would lead to mechanical failures even under steady-state operations. Heat capacity is an important factor in designing cell components. It has a significant effect on the cell temperature distribution. Thermal expansion coefficient (TEC) is a vital criterion for the anode selection and application at high temperatures to avoid mechanical stresses. Therefore, the anode structure and thermophysical properties strongly affect the performance of a SOFC [5]. At present, Ni-YSZ cermet is the most widely adopted material for a SOFC anode [6]. The porous Ni-YSZ cermet is reduced from the NiO-YSZ composite which is sintered from the NiO and YSZ green body [7]. Several experimental studies [8-10] have shown that pellets made of nano-sized NiO and YSZ powders could achieve dense microstructures to extend the TPB length and exhibit high performances after the sintering process. Specifically, as nanocrystalline powders in the size ranging from 5 to 10 nm are becoming increasingly available, denser and more nanostructured materials could be achieved by exploiting the sintering technology [11]. Furthermore, the thermophysical properties of an anode vary with the synthesis route and sintering conditions due to the resultant diverse local morphologies and microstructures. However, experimental studies at the extreme small scales tend to be costly and difficult to be carried out. Due to the limitations of the experimental research, simulation work is more 5

attractive to investigate the nanostructures and properties of the porous anode. The microscopic approaches are an effective way to model the ideal, theoretical or proposed microstructure under clearly defined conditions for anode properties prediction. Xu et al. [12, 13]adopted the molecular dynamics (MD) method to investigate the effect of the porous structure properties, such as the porosity and framework structure, on the sintering process. The results revealed that the porous structure played a critical role in suppression of the sintering. However, as the length scale increases, the computational cost of the MD method becomes relatively high. The Kinetic Monte Carlo (KMC) method, which is at larger length scales and longer time scales, has been used to investigate the reaction mechanism and the oxide ion diffusion process at the anode-electrolyte interface [14, 15]. In literature [16, 17], the density functional theory (DFT) method has been used to investigate the adsorption, dissociation, and diffusion of fuels at the interface. Nevertheless, these approaches failed to simulate the sintering technology and model the microstructure at the nanoscale. Wang et al. [18] firstly adopted the coarse-grained molecular dynamics (CG-MD) method to model the nanostructures and evaluate the thermal properties of the SOFC anode. The results showed that the CG-MD method was an effective tool to simulate the sintering process and assess the thermophysical properties of the anode. However, the effect of the sintering conditions on the sintering performance, the anode nanostructure, and the properties was not clarified in this paper. The present study mainly focused on this. In this study, a CG-MD method is adopted as a modelling technique to reflect the sintering process of the NiO-YSZ porous anode for investigating the sintering performance of the controlled nanostructures. The heat capacity and TEC of the anodic nanostructures sintered at various compositions and conditions are also evaluated and compared with the experimental results. Our study may provide deep insights into the sintering performance at the nanoscale and the local nanostructure of the NiO-YSZ anode. The evaluation of the anode properties provides a rapid and accurate prediction of various compositions and sintering conditions. 6

2. Model systems and computational method The CG-MD method based on the Martini force field can not only keep most properties of the MD simulation, but also greatly increase the length scales and the time scales [19]. The method has been widely used in the numerical research of the biomacromolecules, such as the phospholipid, the cholesterol, the amino acid and the sugar. Researchers also adopted the method to describe the interaction between CG beads of some other macromolecules, such as carbon slab, Pt nano-cluster, NiO and YSZ nanocrystals [18, 20]. The precise details of the method are presented in the following parts. 2.1 All-atom (AA) models AA models are used as the target for the coarse-grained modelization [19]. In the classical AA-MD method, every atom is explicitly simulated and fine details are provided at the atom level. However, this method is limited to smaller systems and shorter time due to the detailed description of the electronic behavior. In a representative Ni-YSZ anode, Ni is the functional component which provides electronic conductivity and catalytic activity at the TPB when the SOFC works. Before reduced in the humidified hydrogen atmosphere with 95% hydrogen and 5% water, Ni exists as its oxide NiO for the sintering process. YSZ mainly provides ionic conductivity to extend the reaction zone. In general, the chemical formula of YSZ is written as (Y2O3)x(ZrO2)(1−2x) with x ranging from 2% to 10%, which means YSZ is composed of the yttria and zirconia. In this study, the AA models are established to generate NiO and YSZ nanocrystals as sintering powders. As the nanocrystalline ceramics have been extensively investigated, it is quite possible to obtain the lattice structure of the NiO and YSZ nanocrystals. The initial atomistic configurations are accomplished in the Atomistic Simulation Environment (ASE) software which is written in the Python programming language with the aim of setting up, steering, and analyzing atomistic simulations [20]. NiO employs the rock salt (NaCl) structure with two nested face-centered cubic 7

(FCC) lattices of nickel and oxide. The NiO nanocrystal is modeled by 5×5×5 lattices, and each lattice constant is 4.195 Å [21]. Fig. 2 illustrates the AA model of the NiO nanocrystal which contains 500 nickel atoms and 500 oxygen atoms. Pure zirconia has three basic polymorphs: a monoclinic phase at room temperature, a tetragonal (distorted fluorite) phase above 1175 °C and below 2370 °C, and a cubic phase (fluorite) above 2370 °C. The tetragonal and cubic phases are rare in nature. However, the addition of yttria to zirconia stabilizes the tetragonal and cubic phases down to room temperature in the regimes of 2-4 mol% and 6-10 mol% of Y2O3, respectively. In this study, 8 mol% yttria stabilized zirconia (8YSZ) is adopted because the oxygen ionic conductivity shows a maximum at ~8.0 mol% Y2O3 [22]. Each cubic ZrO2 lattice exhibits a fluorite structure with the FCC lattice of zirconium atoms and the simple cubic (SC) lattice of oxygen atoms. The lattice parameters of the cubic ZrO2 is a=b=c=5.08 Å, and α=β=γ=90° [23]. The cubic ZrO2 nanocrystal is modeled by 5×5×5 lattices to obtain the AA model with the size of 2.54 nm×2.54 nm×2.54 nm. The addition of Y2O3 to ZrO2 improves the ionic conductivity due to the formation of oxygen vacancies when the host Zr4+ cations are replaced by lower-valent Y3+ cations. A specific number of Zr4+ are substituted by Y3+, which is determined by the composition, at randomly chosen cation sites. At the same time, for every two Y3+, an O2- is removed from a randomly selected anion site. Consequently, the AA model of the YSZ nanocrystal contains 426 zirconium atoms, 74 yttrium atoms and 963 oxygen atoms, as shown in Fig. 3. 2.2 Coarse-grained (CG) models and parameterization To model the SOFC porous anode accurately, tens of thousands of atoms need to be considered. Instead of explicitly describing every atom in the system, the CG-MD method uses CG beads to represent groups of atoms for reducing and simplifying the AA models. The CG-MD method is employed to model the nanostructure of the SOFC anode in two major steps: firstly, transform the AA models to CG models with specific sphere beads at a predefined length scale, which is implemented in the CG builder of 8

the Visual Molecular Dynamics (VMD) package; secondly, set the parameters of the potential interaction between the CG beads with the Martini force field [16], which is a mature force field for CG-MD simulations of the macromolecules. 2.2.1 CG models The AA representation is converted to the CG representation of NiO and YSZ nanocrystals by using a shape-based coarse graining (SBCG) tool in VMD [24]. SBCG adopts a neural network learning algorithm to determine the best location of CG beads. The placement is adjusted (usually slightly) based on the centers of mass of the atomic clusters comprising each bead. The learning algorithm is stochastic, i.e., each time it is used to create a CG model of the same AA structure, the results will be slightly different. However, the overall shape of the nanocrystals is maintained all the same. By this method, the NiO nanocrystals are represented by 20 polar CG beads to ensure the CG beads with the radius of 0.43 nm [20]. Analogously, the YSZ nanoparticles are replaced by 38 beads of radius 0.43 nm. Compared with CG beads of Wang et al. [18] without considering the bead radius, CG models in the present work are more reasonable for the simulations. Fig. 4 shows the AA and CG representations of the two nanocrystals. 2.2.2 Parameterization After building the CG representations, the parameters of the potential interaction between the CG beads are specified. The chemically connected CG beads are described by the bonded interactions and separate CG beads interact through the nonbonded interactions based on the Martini force field. The interaction potential between sites i and j is defined as follows:

V (rij )=Vbonded (rij ) +Vnonbonded (rij )

(1)

where the interaction potential is divided into bonded and nonbonded terms, which are respectively defined in Eqs. (2) and (5).

Vbonded (rij )=Vbond (rij ) +Vangle (θijk ) +Vdihedral (ϕijkl ) 9

(2)

where the bond and angle contributions are specified as harmonic potentials in the following equations. The dihedral term is quite different from the first two potentials because it doesn’t contain the vibrational properties. When Cartesian coordinates are used and the minimum energy structure is employed for determining the vibrational frequencies, the dihedral term can be omitted [25]. 1 b K ij ( rij − bij0 ) 2 2 1 θ (θ ijk − θ ijk0 ) 2 Vangle (θ ijk ) = K ijk 2

Vbond ( rij ) =

(3) (4)

where Kb and Kθ are the bond and angle force constants, respectively; the rij is the distance between bead i and bead j; θijk is the angle among bead i, j and k; b0 and θ0 represent the equilibrium values for the individual terms. In this study, the force constants Kb, Kθ are specified as 1250 kJ·mol-1·nm-2, 25 kJ·mol-1·rad-2 [26], respectively. The involved equilibrium values are determined by the CG models transformed from the AA representation implemented in the SBCG tool of the VMD package. In the bonded interaction potential, a given bead interacts with just a few of its neighbors. However, it interacts with every other bead in the system when the nonbonded interaction potential is taken into account. Therefore, the nonbonded terms are much more computationally costly.

Vnonbonded (rij ) = VLJ (rij ) +VCoulomb (rij )

(5)

where the Lenard-Jones (LJ) term represents the van der Waals interaction between sites i and j, and the Coulomb term represents the electrostatic interaction between charged groups. The two terms are described as follows:

 σ 12  σ 6  VLJ (rij ) = 4ε ij  ij  −  ij    rij   rij     VCoulomb ( rij ) =

10

1

qi q j

4πε 0ε r rij

(6)

(7)

where εij is the depth of the potential at the minimum; σij is the finite distance at which the interaction potential is zero; ε0 is the vacuum dielectric constant (8.85e-12 F·m-1); εr is the relative dielectric constant; qi and qj are the charges of the beads i and j, respectively. In order to calculate the LJ potential, the bead types should be defined firstly. In the Martini force field [26], four main types of interaction sites are taken into account: polar (P), nonpolar (N), apolar (C) and charged (Q). The polar sites represent neutral groups of beads which can be easily dissolved into water, and the apolar sites represent hydrophobic beads, and the nonpolar sites represent the mixed groups which are partly polar and partly apolar. The charged sites refer to the cations and anions. Each particle type has a number of subtypes to allow for a more accurate representation. For the type N and Q, subtypes are defined by a letter which represents the denoting hydrogen-bonding capacities: d=donor, a=acceptor, da=both, o=none. For the main type P and C, subtypes are defined by a number which represents the degree of polarity: from low polarity 1 to high polarity 5. In this study, the SOFC porous anode is composed of two functional components, NiO and YSZ. According to their chemical properties, CG bead of NiO nanocrystal is defined as P4 because ionic bond can be seen as an extreme example of polar covalent bond, and CG bead of YSZ nanocrystal is defined as P5 because of its higher polarity [27]. The levels of LJ interaction between the two CG sites are shown in Table 1. Levels of interaction indicate the well depth in the LJ potential: level O means super attractive, ε=5.6 kJ·mol-1; level I means attractive, ε=5.0 kJ·mol-1. The LJ parameter σ=0.47 nm is adopted to these two interaction levels. In the CG-MD simulation, the LJ parameters can be calculated on the basis of the level of the interaction as shown in Table 2. In this study, the switching function is used for adjusting the LJ potential with the starting point of 0.9 nm and a cut off distance of 1.2 nm. Without charged particles in the system, the Coulomb force is not taken into account.

11

2.3 Coarse-grained molecular dynamics simulations The essence of Molecular dynamics is solving Newton's second law to simulate the movement of the particles. Newton’s equation of motion is described as

Fi = mi

( )

∂V rij ∂ri = − ∂ri ∂t 2

(8)

where mi is the mass of the particle i; ri is the position of the particle i; V(rij) is the interaction potential between particle i and j which is previously described in Eq. (1). On the basis of the above CG representations, the MD method is carried out to simulate the sintering process with the sintering powders NiO and YSZ for modelling the nanostructure of the SOFC anode. Firstly, a target box with the size of 50 nm×50 nm×50 nm which represents the interface region between NiO and YSZ phase is generated. The volume fraction of NiO is regulated by controlling the number of NiO and YSZ beads which are randomly distributed in the target box. 46 vol%, 55 vol% and 67 vol% of NiO for different starting compositions are obtained respectively in this study. Secondly, energy minimization is conducted for the initial structure to achieve the deepest point of the potential energy. The randomly mixed beads are optimized using time steps of 20 fs to accomplish a step of integration procedure. Steep integrator is used because it is robust and easy to implement. After this short energy minimization, the overlapped beads of the system are displaced perfectly. Thirdly, the isothermal-isobaric ensemble (NPT) [28] simulation is performed for the system to model the nanostructure by simulating the sintering process. Moles (N), pressure (P) and temperature (T) remain constant during the simulation. The temperature coupling is controlled by the V-rescale algorithm, which mimics a weak coupling to an external heat bath with a given temperature T0 representing the sintering temperature. The weak coupling algorithm is applied separately for each component (NiO and YSZ) with a time constant of 1.0 ps. The pressure coupling is controlled by the Berendsen algorithm with the isotropic couple type. The coupling algorithm is applied with a time constant of 2.0 ps and a reference pressure representing the sintering pressure. This equilibrating 12

process is carried out for 30 ns using 20 fs time steps. Finally, the thermophysical properties of the sintered nanostructure are achieved by the annealing procedure. During the simulation, the results including the positions, the velocities and the forces are saved every 5000 steps (100 ps) for later analysis. The energy information is saved every 2000 steps (40 ps) to the log file and the energy file. All these MD simulations are performed in the command-line interface of the GROMACS simulation package by operating a number of recognized files [29]. All graphical images presented in this study are prepared using VMD [30].

3. Results and discussion In this section, the results of our MD simulation of the system at different components and various conditions are presented. The sintering behavior of the controlled

nanostructures,

the

sintered

nanostructures

and

the

predicted

thermophysical properties (volumetric heat capacity and linear thermal expansion coefficient) will be considered separately in the following sections. 3.1 Sintering performance In general, the sintering technique is traditionally implemented at high temperatures but below the melting points of the materials. To understand the sintering mechanism of NiO-YSZ anode and the effect of sintering conditions on the sintering process, the porous anode is modelled at varying sintering temperatures ranging from 1500 K to 2000 K [31-33]. The sintering pressure from 10 to 40 bar adopted mainly depends on the sintering simulation studies of Xu et al. [12] and Wang et al. [18, 27]. We also admit that, the sintering pressure from 10 to 40 bar is chosen as an artifact in this simulation, and the effect of the sintering pressure could provide a rapid prediction for the nanoscale sintering at such conditions. Figs. 5 and 6 reflect the sintering process of NiO-YSZ anode within 30 ns at various sintering temperatures and pressures, respectively. The red domains represent the NiO phase, and the blue parts represent the YSZ phase. In addition, the pores represent the gas phase. In view 13

of the two-dimensional (2D) x-y cross sections, the phenomenon of self-assembly could be easily observed. The NiO beads aggregate together to form the NiO phase for conducting electrons. The YSZ beads aggregate together to form the YSZ phase for delivering oxygen ions. Along with the self-assembly phenomenon during the sintering process, pores form to diffuse the reactant gas and water, and the volume shrinkage of the box occurs. The initial box size of cubic 50 nm containing 46 vol% NiO is reduced to a smaller one at different sintering conditions. The box size decreases with the simulation time and eventually achieves the equilibrium state. This is mainly because sintering is the process where interparticle pores in a granular material are eliminated by atomic diffusion driven by capillary forces [11]. By comparing Figs. 5 (a), 5 (b) and 5 (c), the box size at 30 ns increases as the sintering temperature increases. The main reason is the capillary pressure which drives the volume shrinkage decreases with an increase in temperature [34]. Figs. 6 (a), 6 (b) and 6 (c) depict that an increase in the sintering pressure leads to a decrease in the final box size. These results indicate that sintering at relatively low temperature and high pressure could contribute to the densification of the sintered nanostructure. It is possible to control the sintering temperature and pressure to improve the anode properties. 3.2 Nanostructure properties After the MD running for simulating the sintering process, the final nanostructures of NiO-YSZ anode are obtained. Fig. 7 (a) illustrates the nanostructure of the NiO0.30YSZ modelled at 1623 K in this study. As mentioned previously, the red domains, blue parts and pores represent the NiO, YSZ and gas phases, respectively. Fig. 7 (b) shows the Transmission Electron Microscopy (TEM) micrograph of Ni0.30YSZ cermet sintered at the same temperature [35]. The black and grey grains represent the YSZ and Ni phases which are in the range of 5-12 nm. In addition, the white parts represent the gas phase. Figs. 7 (a) and (b) show the common morphology of the porous anode of SOFC which depends on the sintering technology. Both of 14

them contain three phases and the length scales in these two figures are quite consistent, which indicates the sintering process simulation at the nanoscale is rather meaningful. Compared with the experimental result, the nanostructure obtained from our study could supply more details of the TPB regions and the self-assembly phenomenon in the anode. Furthermore, it is easily observed from Fig. 7 (a) that the sintering at the nanoscale is beneficial to distribute the functional components equally. This is meaningful for increasing the TPB lengths to improve the electrochemical reactions. However, the TEM micrograph has a small difference from our modelling nanostructure. Fig. 7 (b) shows the composite which is reduced under flowing hydrogen at 1073 K. The porosity of the reduced structure is higher than that of the sintered structure. In the present study, sintered nanostructures (40 bar-1600 K) with NiO contents of 46, 55 and 67 vol% are shown in Fig. 8. These systems are further applied to evaluate the properties of the anode. To assess properties of the modelling nanostructures, densities and porosities of the sintered anodes at various sintering conditions (Figs. 9 and 10) are calculated. Sintering at nanoscale in this study has obtained nanostructures with reasonable density and porosity. As shown in Fig. 9, it can be observed that the lower the sintering temperature, the higher the sintered density and the lower the porosity. It can also be noted that the NiO content of the modelling nanostructures has a slight effect on the sintered density and porosity with comparable sintering temperature. In analogy to the effect of sintering temperature on the sintered density and porosity, the sintering pressure has an opposite effect (Fig. 10). Sintered density increases and porosity decreases with an increase in sintering pressure. These could be concluded that sintering at relatively low temperature and high pressure is beneficial to improve sinterability of the sintered nanostructure. Many properties of the SOFC anode, such as mechanical strength, length of TPB and pore percolation benefit from relatively high density and appropriate porosity. The control of the sintering temperature and 15

pressure could adjust the nanostructures to achieve the required density and porosity. 3.3 Predicted thermophysical properties 3.3.1 Volumetric heat capacity Heat capacity is an important factor in designing cell components. As the operation conditions change, a high heat capacity of the solid phase brings about a reduction of the anode temperature change. The effect of the sintering temperature, pressure and NiO content on the anode heat capacity could supply a rapid prediction for the experimental research to obtain the desired heat capacity. In this study, the volumetric heat capacity (VHC), which describes the ability of a given volume of a substance to store internal energy while undergoing a given temperature change without undergoing a phase change, is evaluated. The VHC is described as follows: CV =

∇E ∇T ⋅V

(9)

where CV represents the VHC (J·cm-3·K-1); E is the internal energy (J); T represents the temperature (K); V is the volume (cm3). In order to compute CV, the system is expanded by gradually increasing the temperature from 300 K to 1600 K in a canonical ensemble (NVT). Fig. 11 illustrates the predicted VHCs of the sintered anodes with different NiO contents at various sintering temperatures. The values of VHC with a specific composition show a negative correlation to the sintering temperatures. The VHC of NiO0.67YSZ is obviously higher than that of NiO0.46YSZ and NiO0.55YSZ with a comparable sintering temperature. Fig. 12 depicts the predicted VHCs at various sintering pressures. The values of VHC show a positive correlation to the sintering pressures and the NiO contents. At a relatively low sintering pressure, the VHCs of NiO0.46YSZ and NiO0.55YSZ are very close and slightly lower than those of NiO0.67YSZ. As the sintering pressure increases, the VHC of NiO0.67YSZ increases much higher than that of NiO0.46YSZ and NiO0.55YSZ. In this work, the VHCs of the SOFC anodes modelled at 1500 K to 2000 K and 10 bar to 40 bar are in the range of 0.61-1.55 J·cm-3·K-1. 16

Dulong-Petit law describes the heat capacity per weight which is close to a constant value, after it has been multiplied by a number representing the presumed relative atomic weight of the solid elements, as shown in Eq. (10). The description of the specific heat capacity (SHC) is accurate for most of the crystals under high temperature. The VHC is directly proportional to the bulk density and SHC of the materials. Furthermore, the VHC is subject to the molar mass and bulk density at the same time (shown in Eq. 11). Therefore, the increasing values of VHC are due to the densification of the system with the decreasing sintering temperature and increasing sintering pressure. (10)

CM = 3R CV = ρ C = ρ

3R M

(11)

where C represents the SHC (J·g-1·K-1); M is the molar mass (g·mol-1); R is the gas constant (J·K-1·mol-1); ρ is the bulk density (g·cm-3). Table 3 illustrates the predicted VHCs calculated by linear combination and CG-MD method. VHCs of NiO and 8YSZ are evaluated by the Dulong-Petit law separately. For the same sintered composites, VHCs by the CG-MD method at various sintering conditions are obviously lower than those by simple linear combination of NiO and YSZ VHCs. The main reason causing the above situation is as follows: after the sintering process, pores existing in the sintered composite could affect the VHC of the sintered composite, which can not be considered by the simple linear combination method. Furthermore, various sintering conditions strongly affect the porous microstructure of the sintered composite, which determines the VHC of the sintered composite, as shown in Figs. 11 and 12. CG-MD method is an effective way to evaluate the composite VHC at various sintering conditions, which could supply a prediction for the experimental research. In Table 4, the values of VHC for the SOFC anodes are compared to the open literature data [36-38]. The simulation results in the range of 0.61-1.55 J·cm-3·K-1 agree well with the open literature data. As the heat capacity on a volumetric basis in 17

solid materials varies more widely, the deviation is small and acceptable. In addition, compared with the results of Wang et al. [18], the values of VHC in this study are more reasonable. Modification of the CG models and annealing procedure for calculating VHCs could contribute to the differences between these two cases. Therefore, the CG-MD method applied in this study has been validated, and it would be an effective way to simulate the sintering process and predict relevant properties. As expected, the results in this study could supply a rapid and economical way for the experimental studies to decide the sintering conditions and components for obtaining nanostructures with better performances. 3.3.2 Thermal expansion coefficient To accurately estimate thermal expansion of the NiO-YSZ cermet-anode during the SOFC works at a wide range from room temperature to above 1000 K, the linear TECs of the modelling nanostructures are calculated. For solids, linear TEC describes change in one dimension in response to change in temperature at a constant pressure, through heat transfer. Therefore, the TEC is described as follows:

αL =

1  ∂L    L  ∂T  p

(12)

where αL represents the TEC (K-1); L is a particular length measurement (m); ∂L/∂T represents the rate of change of that linear dimension per unit change in temperature (m·K-1). An annealing procedure is adopted in this study to compute the TEC by gradually increasing the temperature in a NPT ensemble. Fig. 13 shows the box sizes of NiO0.46YSZ, NiO0.55YSZ and NiO0.67YSZ composites modelled at 10 bar-1600 K with the variation of temperature. Box size increases with temperature, which means thermal expansion occurs during the heating process. Slopes of box size vs. temperature fitted curves change from 300 K to 1100 K. According to the slopes, the TECs over the same temperature range could be obtained. As shown in Fig. 14, the TECs of the composites with a specific NiO content tend to increase with temperature because of volume expansion due to the 18

formation of point defects at high temperature [39]. In addition, as the NiO content increases, the TECs of the composites decrease. This can be explained that microstructure of the initial configuration of NiO and YSZ particles in a composite strongly affects the TECs with NiO of 46, 55 and 67 vol% [39]. Nanoparticles sintering with the particular NiO content could lead to decrease in TECs of the composites as NiO increases. Fig. 15 depicts the negative effect of sintering temperature and the positive effect of sintering pressure on TECs at 300 K of NiO0.46YSZ nanostructures obtained from relevant sintering conditions. Compared with TECs of 8YSZ (5.0-11.5×10−6·K−1) and NiO (11.5-18.5×10−6·K−1), which are reported in literature [39], the TECs of composites calculated in this study (4.54-20.4×10 −6·K−1) could be considered in a reasonable range. To avoid cell cracking and delamination during cell fabrication and operation, the TEC of the anode is required to match with that of the electrolyte (8YSZ). Since the TEC of the electrolyte is relatively low, adopting a higher sintering temperature or an appropriate sintering pressure for modelling the anode is beneficial for TEC match between anode and electrolyte.

4. Conclusions In this study, the sintering process of the typical NiO-YSZ anode is investigated and thermophysical properties of the sintered nanostructures are evaluated through CG-MD simulations. The effects of the sintering conditions and composition on the sintering performance and nanostructure properties are also discussed. The main conclusions can be summarized as follows: (1) During the sintering process of the NiO-YSZ anode, volume shrinkage occurs along with the self-assembly phenomenon driven by capillary forces. The volume shrinkage increases with the decreasing sintering temperature and the increasing sintering pressure. Sintering at relatively low temperature and high pressure could contribute to the densification of the sintered nanostructure. 19

(2) Compared with the experimental sintered nanostructure, the SOFC anode obtained from this study shows the common morphology. Sintering at the nanoscale is beneficial to distribute the functional components equally, which is meaningful for increasing the TPB lengths to improve the electrochemical reactions. In addition, the control of the sintering temperature, pressure and NiO content might adjust the nanostructures to obtain desired sintered density and porosity. (3) Sintering temperature has a negative effect and sintering pressure has a positive effect on the VHC and TEC of the SOFC anode. Nanoparticles sintering with the particular NiO content lead to increase in VHC and decrease in TEC as NiO increases. The predicted VHCs in the range of 0.61-1.55 J·cm-3·K-1 agree well with the open literature data. TECs of composites calculated in this study (4.54-20.4×10−6·K−1) could be considered in a reasonable range. Compared with expensive experimental studies at extreme small scales, results in this study provide a rapid and accurate prediction of various compositions and sintering conditions for obtaining the optimal thermophysical properties of the SOFC anode.

Acknowledgements This work is financially supported by the National Science Foundation of China (Grant No. 51536007) and the National Science Foundation of China for International Cooperation and Exchange (Grant No. 51120165002).

References [1] H. Wen, J.C. Ordonez, J.V.C. Vargas, Optimization of single SOFC structural design for maximum power, Appl. Therm. Eng. 50 (2013) 12-25. [2] Q.Y. Chen, M. Zeng, J. Zhang, Q.W. Wang, Optimal design of bi-layer interconnector for SOFC based on CFD-Taguchi method, Int. J. Hydrogen Energy 35 (2010) 4292-4300. [3] M. Liu, P.V. Aravind, The fate of tars under solid oxide fuel cell conditions: A review, Appl. Therm. Eng. 70 (2014) 687-693. 20

[4] M. Yan, M. Zeng, Q.Y. Chen, Q.W. Wang, Numerical study on carbon deposition of SOFC with unsteady state variation of porosity, Appl. Energy 97 (2012) 754-762. [5] M. Radovic, E. Lara-Curzio, R.M. Trejo, H. Wang, W.D. Porter, Thermophysical properties of YSZ and Ni-YSZ as a function of temperature and porosity, in: A. Wereszczak，E. Lara-Curzio，N.P. Bansal, Advances in Solid Oxide Fuel Cells II: Ceramic Engineering and Science Proceedings, John Wiley & Sons, Inc., Hoboken, NJ, USA, 2009, pp. 79. [6] J.R. Wilson, W. Kobsiriphat, R. Mendoza, H.Y. Chen, J.M. Hiller, D.J. Miller, K. Thornton, P.W. Voorhees, S.B. Adler, S.A. Barnett, Three-dimensional reconstruction of a solid-oxide fuel-cell anode, Nat. Mater. 5 (2006) 541-544. [7] N.Q. Minh, Ceramic fuel cells, J. Am. Ceram. Soc. 76 (1993) 563-588. [8] S.D. Kim, H. Moon, S.H. Hyun, J. Moon, J. Kim, H.W. Lee, Ni-YSZ cermet anode fabricated from NiO-YSZ composite powder for high-performance and durability of solid oxide fuel cells, Solid State Ionics 178 (2007) 1304-1309. [9] Y.J. Leng, S.H. Chan, K.A. Khor, S.P. Jiang, P. Cheang, Effect of characteristics of Y2O3/ZrO2 powders on fabrication of anode-supported solid oxide fuel cells, J. Power Sources 117 (2003) 26-34. [10] M. Marinšek, K. Zupan, Microstructure evaluation of sintered combustion-derived fine powder NiO-YSZ, Ceram. Int. 36 (2010) 1075-1082. [11] I.W. Chen, X.H. Wang, Sintering dense nanocrystalline ceramics without final-stage grain growth, Nature 404 (2000) 168-171. [12] J.X. Xu, S. Bai, Y. Higuchi, N. Ozawa, K. Sato, T. Hashidac, M. Kubo, Multi-nanoparticle model simulations of the porosity effect on sintering processes in Ni/YSZ and Ni/ScSZ by the molecular dynamics method, J. Mater. Chem. A. 3 (2015) 21518-21527. [13] J.X. Xu, R. Sakanoi, Y. Higuchi, N. Ozawa, K. Sato, T. Hashidac, M. Kubo, Molecular dynamics simulation of Ni nanoparticles sintering process in Ni/YSZ 21

multi-nanoparticle system, J. Phys. Chem. C. 117 (2013) 9663-9672. [14] X. Wang, K.C. Lau, C.H. Turner, B.I. Dunlap, Kinetic Monte Carlo simulation of the elementary electrochemistry in a hydrogen-powered solid oxide fuel cell, J. Power Sources 195 (2010) 4177-4184. [15] R. Pornprasertsuk, J. Cheng, H. Huang, F.B. Prinz, Electrochemical impedance analysis of solid oxide fuel cell electrolyte using kinetic Monte Carlo technique, Solid State Ionics 178 (2007) 195-205. [16] D.B. Ingram, S. Linic, First-principles analysis of the activity of transition and noble metals in the direct utilization of hydrocarbon fuels at solid oxide fuel cell operating conditions, J. Electrochem. Soc. 156 (2009) B1457-B1465. [17] M.H. Weng, H.T. Chen, Y.C. Wang, S.P. Ju, J.G. Chang, M.C. Lin, Kinetics and mechanisms for the adsorption, dissociation, and diffusion of hydrogen in Ni and Ni/YSZ slabs: A DFT study, Langmuir 28 (2012) 5596-5605. [18] Y.F. Wang, J.L. Yuan, B. Sundén, Y.L. Hu, Coarse-grained molecular dynamics investigation of nanostructures and thermal properties of porous anode for solid oxide fuel cell, J. Power Sources 254 (2014) 209-217. [19] S. Kakac, A. Pramuanjaroenkij, X.Y. Zhou, A review of numerical modeling of solid oxide fuel cells, Int. J. Hydrogen Energy 32 (2007) 761-786. [20] Y. Xiao, M.L. Dou, J.L. Yuan, M. Hou, W. Song, B. Sundén, Fabrication process simulation of a PEM fuel cell catalyst layer and its microscopic structure characteristics, J. Electrochem. Soc. 159 (2012) B308-B314. [21] Y. Du, W.N. Wang, X.W. Li, J. Zhao, J.M. Ma, Y.P. Liu, G.Y. Lu, Preparation of NiO nanoparticles in microemulsion and its gas sensing performance, Mater. Lett. 68 (2012) 168-170. [22] K.C. Lau, B.I. Dunlap, Molecular dynamics simulation of yttria-stabilized zirconia (YSZ) crystalline and amorphous solids, J. Phys.: Condens. Matter. 23 (2011) 035401. [23] K.C. Lau, B.I. Dunlap, Lattice dielectric and thermodynamic properties of yttria 22

stabilized zirconia solids, J. Phys.: Condens. Matter. 21 (2009) 145402. [24] http://www.ks.uiuc.edu/Research/vmd/plugins/cgtools/ [25] A.D. MacKerell, D. Bashford, M. Bellott, R. Dunbrack, J.D. Evanseck, M.J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, All-atom empirical potential for molecular modeling and dynamics studies of proteins, J. Phys. Chem. B. 102 (1998) 3586-3616. [26] S.J. Marrink, H.J. Risselada, S. Yefimov, D.P. Tieleman, A.H. De Vries, The MARTINI force field: coarse grained model for biomolecular simulations, J. Phys. Chem. B. 111 (2007) 7812-7824. [27] Y.F. Wang , J.L. Yuan, B. Sundén, Y.L. Hu, Reconstruction of Anode Nanostructures for Solid Oxide Fuel Cells, Proceedings of the ASME 2013 11th Fuel Cell Science, Engineering, and Technology Conference, ESFuelCell2013, Minneapolis, Minnesota, USA, July 14-19, 2013, p. V001T02A005. [28] A.

Das,

H.C.

Andersen,

The

multiscale

coarse-graining method.

V.

Isothermal-isobaric ensemble. J. Chem. Phys. 132 (2010) 164106. [29] D. van der Spoel, E. Lindahl, B. Hess, and the GROMACS development team. GROMACS User Manual version 4.6.5. www.gromacs.org (2013). [30] http://www.ks.uiuc.edu/ [31] Y.X. Zhang, C.R. Xia, M. Ni, Simulation of sintering kinetics and microstructure evolution of composite solid oxide fuel cells electrodes, Int. J. Hydrogen Energy 37 (2012) 3392-3402. [32] Z.J. Jiao, N. Shikazono, Simulation of nickel morphological and crystal structures evolution in solid oxide fuel cell anode using phase field method, J. Electrochem. Soc. 161 (2014) F577-F582. [33] X.X. Liu, C.L. Martin, G. Delette, J. Laurencin, D Bouvard, T. Delahaye, Microstructure of porous composite electrodes generated by the discrete element method, J. Power Sources 196 (2011) 2046-2054. [34] S.A. Grant, A. Salehzadeh, Calculation of temperature effects on wetting 23

coefficients of porous solids and their capillary pressure functions, Water Resour. Res. 32 (1996) 261-270. [35] S.T. Aruna, M. Muthuraman, K.C. Patil, Synthesis and properties of Ni-YSZ cermet: anode material for solid oxide fuel cells, Solid State Ionics 111 (1998) 45-51. [36] A.C. Burt, I.B. Celik, R.S. Gemmen, A.V. Smirnov, A numerical study of cell-to-cell variations in a SOFC stack, J. Power Sources 126 (2004) 76-87. [37] V.M. Janardhanan, O. Deutschmann, Numerical study of mass and heat transport in solid-oxide fuel cells running on humidified methane, Chem. Eng. Sci. 62 (2007) 5473-5486. [38] C. Schluckner, V. Subotić, V. Lawlor, C. Hochenauer, Three-dimensional numerical and experimental investigation of an industrial-sized SOFC fueled by diesel reformat-Part II: Detailed reforming chemistry and carbon deposition analysis, Int. J. Hydrogen Energy 40 (2015) 10943-10959. [39] M. Mori, T. Yamamoto, H. Itoh, Thermal Expansion of Nickel-Zirconia Anodes in Solid Oxide Fuel Cells during Fabrication and Operation, J. Electrochem. Soc. 145 (1998) 1374-1381.

24

Figure captions Fig. 1 Schematic of SOFC components Fig. 2 AA model for the NiO nanocrystal (a) along the (001) direction (b) with the size of 2.0975 nm×2.0975 nm×2.0975 nm Fig. 3 AA model for the YSZ nanocrystal (a) along the (001) direction (b) with the size of 2.54 nm×2.54 nm×2.54 nm Fig. 4 AA and CG representations of the (a) NiO nanocrystal (b) YSZ nanocrystal Fig. 5 Time evolutions of sintering processes at 40 bar and (a) 1500 K (b) 1600 K (c) 1700K Fig. 6 Time evolutions of sintering processes at 1600 K and (a) 20 bar (b) 30 bar (c) 40 bar Fig. 7 (a) Snapshot of NiO0.30YSZ nanostructure sintered at 1623 K. (b) TEM micrograph of Ni0.30YSZ cermet sintered at 1623 K [35] Fig. 8 Sintered nanostructures with NiO contents of (a) 46 vol% (b) 55 vol% (c) 67 vol% Fig. 9 Sintered densities and porosities of the porous anodes modelled at varying sintering temperatures Fig. 10 Sintered densities and porosities of the porous anodes modelled at varying sintering pressures Fig. 11 Predicted VHCs of the porous anodes modelled at varying sintering temperatures Fig. 12 Predicted VHCs of the porous anodes modelled at varying sintering pressures Fig. 13 Box sizes of NiO0.46YSZ, NiO0.55YSZ and NiO0.67YSZ anodes with the variation of temperature Fig. 14 TECs of NiO0.46YSZ, NiO0.55YSZ and NiO0.67YSZ anodes with the variation of temperature Fig. 15 Predicted TECs at 300 K of the porous anodes modelled at varying sintering temperatures and pressures Table captions Table 1 Levels of LJ interaction between NiO and YSZ Table 2 LJ interaction parameters used in the CG-MD simulation (C6=4 εijσij6 and C12=4εijσij12) Table 3 Predicted VHCs: linear combination versus CG-MD Table 4 Comparison of predicted VHCs with the open literature data 25

Fig. 1 Schematic of SOFC components

26

(a)

(b)

Fig. 2 AA model for the NiO nanocrystal (a) along the (001) direction (b) with the size of 2.0975 nm×2.0975 nm×2.0975 nm

27

(a)

(b)

Fig. 3 AA model for the YSZ nanocrystal (a) along the (001) direction (b) with the size of 2.54 nm×2.54 nm×2.54 nm

28

(a)

(b)

Fig. 4 AA and CG representations of the (a) NiO nanocrystal (b) YSZ nanocrystal

29

t=0.1 ns

t=0.5 ns

t=2.5 ns

t=12.5 ns

t=30 ns

(a)

(b)

(c) Fig. 5 Time evolutions of sintering processes at 40 bar and (a) 1500 K (b) 1600 K (c) 1700K

30

t=0.1 ns

t=0.5 ns

t=2.5 ns

t=12.5 ns

t=30 ns

(a)

(b)

(c) Fig. 6 Time evolutions of sintering processes at 1600 K and (a) 20 bar (b) 30 bar (c) 40 bar

31

(a)

(b)

Fig. 7 (a) Snapshot of NiO0.30YSZ nanostructure sintered at 1623 K. (b) TEM micrograph of Ni0.30YSZ cermet sintered at 1623 K [35]

32

(a)

(b)

(c)

Fig. 8 Sintered nanostructures with NiO contents of (a) 46 vol% (b) 55 vol% (c) 67 vol%

33

100

NiO0.46YSZ-40 bar NiO0.55YSZ-40 bar

3

Density (kg·m-3)

5.0x10

NiO0.67YSZ-40 bar

90 80 70

4.5x103

60 3

4.0x10

50 40

3

3.5x10

Porosity (%)

5.5x103

30 3.0x103

20

3

2.5x10 1400

1500

1600

1700

1800

1900

2000

10 2100

Sintering Temperature (K)

Fig. 9 Sintered densities and porosities of the porous anodes modelled at varying sintering temperatures

34

6x103

100

NiO0.46YSZ-1600 K

80

3

4x10

70

3x103

60 50

2x103

Porosity (%)

Density (kg·m-3)

90

NiO0.55YSZ-1600 K NiO0.67YSZ-1600 K

5x103

40 1x103

30

0

20 10

15

20

25

30

35

40

Sintering Pressure (bar)

Fig. 10 Sintered densities and porosities of the porous anodes modelled at varying sintering pressures

35

VHC (J·cm-3·K-1)

1.7 1.6

NiO0.46YSZ-40 bar

1.5

NiO0.55YSZ-40 bar NiO0.67YSZ-40 bar

1.4 1.3 1.2 1.1 1.0 0.9 1500

1600

1700

1800

1900

2000

Sintering Temperature (K)

Fig. 11 Predicted VHCs of the porous anodes modelled at varying sintering temperatures

36

)

-1

NiO0.46YSZ-1600K

1.2

NiO0.55YSZ-1600K NiO0.67YSZ-1600K

-3

VHC (J·cm ·K

1.4

1.0

0.8

0.6 10

15

20

25

30

35

40

Sintering Pressure (bar)

Fig. 12 Predicted VHCs of the porous anodes modelled at varying sintering pressures

37

Box Size (nm)

43.55

NiO0.46YSZ

46.95

43.50

46.90

43.45

46.85

43.40 43.35

NiO0.55YSZ

52.70

NiO0.67YSZ

52.65 52.60

46.80

52.55

46.75 52.50 46.70

43.30

52.45 46.65

43.25 300 600 900 1200

Temperature (K)

52.40 300 600 900 1200 300 600 900 1200

Temperature (K)

Temperature (K)

Fig. 13 Box sizes of NiO0.46YSZ, NiO0.55YSZ and NiO0.67YSZ anodes with the variation of temperature

38

TEC (10-6·K-1)

14 12

NiO0.46YSZ NiO0.55YSZ

10

NiO0.67YSZ

8 6 4 2 0 200 300 400 500 600 700 800 900 1000 1100 1200

Temperature (K)

Fig. 14 TECs of NiO0.46YSZ, NiO0.55YSZ and NiO0.67YSZ anodes with the variation of temperature

39

Sintering Temperature (K) 1400 22

1500

1600

1700

1800

1900

30

35

2000

20 18

TEC (10-6·K-1)

16 14 12 10 8 6 4 2 10

15

20

25

40

Sintering Pressure (bar)

Fig. 15 Predicted TECs at 300 K of the porous anodes modelled at varying sintering temperatures and pressures

40

Table 1 Levels of LJ interaction between NiO and YSZ Bead type

NiO (P4)

YSZ (P5)

NiO (P4)

I

O

YSZ (P5)

O

O

41

Table 2 LJ interaction parameters used in the CG-MD simulation (C6=4εijσij6 and C12=4εijσij12) Level

C6

C12

O

0.241

0.0026

I

0.216

0.0023

42

Table 3 Predicted VHCs: linear combination versus CG-MD Component

Molar mass/ g·mol-1

Density/ g·cm-1

VHC/ J·cm-3·K-1

NiO

74.69

6.67

2.23

8YSZ

131.42

5.90

1.12

Predicted VHC/ J·cm-3·K-1

Predicted VHC/ J·cm-3·K-1

(linear combination)

(CG-MD)

NiO0.46YSZ

1.63

0.63-1.43

NiO0.55YSZ

1.73

0.61-1.43

NiO0.67YSZ

1.86

0.65-1.55

Composite

43

Table 4 Comparison of predicted VHCs with the open literature data Literature

Predicted VHC/ J·cm-3·K-1

Burt et al. [36]

1.20

Janardhanan et al. [37]

1.49

Schluckner et al. [38]

1.80

Wang et al. [18]

0.4-5

Present study

0.61-1.55

44