- Email: [email protected]
dehydriding reactions: An OED study
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 1 ( 2 0 1 6 ) 5 1 7 6 e5 1 8 7
Available online at www.sciencedirect.com
ScienceDirect journal homepage: www.elsevier.com/locate/he
Improving the estimation quality of parameters in kinetic models for hydriding/dehydriding reactions: An OED study Fusheng Yang a,*, Francesco Ciucci b, Zhen Wu a, Zaoxiao Zhang a,c, Yuqi Wang d a
School of Chemical Engineering and Technology, Xi'an Jiaotong University, Xi'an, 710049, China Department of Aerospace and Mechanical Engineering & Department of Chemical and Biomolecular Engineering, Hong Kong University of Science and Technology, Hong Kong c State Key Laboratory of Multiphase Flow in Power Engineering, Xi'an Jiaotong University, Xi'an, 710049, China d School of Chemical Engineering, Northwest University, Xi'an, 710069, China b
article info
abstract
Article history:
Among the hydrogen storage properties of solid materials, e.g., metal hydrides, the reaction
Received 31 October 2015
kinetics is of great importance and is often investigated by fitting reaction rate curves to given
Received in revised form
models. Thanks to such models, key parameters such as activation energies can be found.
14 January 2016
However, it is noteworthy that the parameters thus obtained often show great disparities
Accepted 15 January 2016
even for the same material. In addition to experimental inaccuracies, these disparities may
Available online 6 February 2016
also be attributed to the improper processing of the data and to suboptimal experimental design. In this paper, we present some theoretical tools based on optimal experimental
Keywords:
design (OED) to analyze and optimize the uncertainty on the estimated parameters. The
Metal hydride
uncertainty is derived from asymptotic statistics and it is shown to be closely linked to the
Hydriding/dehydriding kinetics
sensitivity of the experiment to the parameters. We found that by applying a new OED-based
Optimal experimental design
procedure, the estimate of the kinetic parameters can be effectively improved. Copyright © 2016, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved.
Introduction Due to the unique properties of reversible reaction with large amount of hydrogen, metal hydride (MH) materials have shown great application potential and attracted wide attention in recent years. The possible applications of MH include hydrogen storage [1e4], static hydrogen compression [5,6], gas mixture separation/purification [7e9], thermal energy storage [10,11], and refrigeration/heat pumps [12,13].
To realize such practical applications, we need to have an in-depth understanding on the fundamental physical and chemical properties of the MH material, including the hydriding/dehydriding kinetic properties. Like many other gasesolid reactions, the hydriding reaction is heterogeneous and generally thought to proceed in the following steps [14e16]: (1) Physisorption of H2 molecules on the material surface; (2) Dissociation of physisorbed H2 molecules and chemisorption; (3) Penetration of H atoms to the subsurface;
* Corresponding author. Tel.: þ86 29 82668566; fax: þ86 29 82660689. E-mail address: [email protected] (F. Yang). http://dx.doi.org/10.1016/j.ijhydene.2016.01.086 0360-3199/Copyright © 2016, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved.
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 1 ( 2 0 1 6 ) 5 1 7 6 e5 1 8 7
(4) Diffusion of H atoms in the hydride layer either by an interstitial or a vacancy mechanism; (5) Hydride formation at the a/b phase (i.e. solid solution and hydride) interface. On the other hand, the dehydriding process will follow the abovementioned steps in the reverse order. As is well known, if one of the steps is significantly slower than any other steps, it will be considered as rate-controlling step dictating the actual reaction rate. From the reaction mechanism and ratecontrolling step, various kinetic models can be formulated accordingly [14]. By careful processing the experimental data, relevant physical parameters are obtained. Generally, data fitting is conducted for several possible model candidates, and the best one according to a certain criterion will be the “true” model, which is thought to precisely represent the intrinsic reaction behavior of the material under discussion. After that the unknown parameters in the selected model can be found. Obviously, the abovementioned procedure provides insight into the reaction kinetics pertaining to the target material and facilitates the modeling of the corresponding MH-based application system [2,5,13]. For example, Bershadsky et al. [17] observed the kinetic behavior of TiFe0.8Ni0.2 alloy under strictly controlled conditions, and fitted the data against 3 models to obtain the activation energy and pre-exponential factor. However, no single model seemed to represent the data perfectly. Similarly, Oh et al. [18] conducted stepwise volumetric experiments to study the hydriding kinetics of LaNi4.5Al0.5, and found the rate-controlling step to be dissociative chemisorption of hydrogen molecules. The activation energy of the reaction was determined to be 39.8 kJ/mol H2. Inomata et al. [19] carried out an experimental study on the hydriding/dehydriding kinetics of LaNi5, and concluded that the rate-controlling step of hydriding reaction during the initial period is the nucleation and growth of hydride phase, and then changed to hydrogen diffusion. Through isothermal experiments on Mg1.9Al0.1Ni, Li et al. [20] discovered that 3-D diffusion dictates the hydriding/dehydriding kinetics of the material. They also estimated that the activation energies are 52 ± 2 and 48 ± 1 kJ/mol H2 for both reactions, respectively. Yang et al. [21] investigated the reaction kinetics of NaAlH4 in certain decomposition steps. Their results justified a moving boundary reaction model for hydriding kinetics, while the dehydriding reaction showed quite complicated behavior and could not be well fitted to the models they discussed. Muthukumar et al. [22] measured the hydriding kinetics of several La-based alloys including LaNi5, and their results indicated that diffusion of hydrogen atom is the ratecontrolling step for all three alloys in the plateau region. The activation energies of the alloys are mostly around 30 kJ/mol H2. Paya et al. [23] utilized gravimetric apparatus to obtain hydriding kinetics data of LaNi5 alloy. A Johnson-Mehl-Avrami (JMA) model suggesting nucleation and growth mechanism was found able to describe the reaction behavior in the plateau region satisfactorily. Cao et al. [24] investigated the hydriding kinetics of LaNi5 in repeated cycles using volumetric method and fit the data with JMA model. The activation energy of the material was determined to be 19.4 kJ/mol H2. From the literature review, we can see that even for the well-known material LaNi5, the estimated values of the key
5177
parameters in the hydriding/dehydriding kinetics model, e.g., activation energies, show remarkable disparities [25]. For example, quite a few research groups fit the hydrogenation rate curve successfully to a JMA model suggesting nucleationgrowth mechanism, while the activation energy they determined ranged from 19.5 to 36.8 kJ/mol H2. Such disparity can be partly attributed to the poor experimental settings and varied instrumentation/procedure, on which quite a few previous studies were focused. In an early investigation, Wang and Suda [26] summarized the factors affecting the accuracy of hydriding/dehydriding kinetics measurement by the volumetric method, i.e., thermal effects, relationship between gas holder volume and sample mass, history and number of hydriding/dehydriding cycles, surface state of the particles, particle size, and experimental conditions. Some proposals were made for precise kinetic determinations accordingly. Later, Blach and Gray [27] analyzed the error of hydrogen intake measurement by a Sieverts apparatus, and the importance of a large apparatus volume was highlighted. Moreover, some rules-of-thumb to limit the measurement uncertainty were also given by the authors. Broom [28] presented an extensive review on the hydrogen sorption measurements by various methods, i.e., volumetric method, gravimetric method, and temperature-programmed desorption. A number of practical issues related to the measurement accuracy were discussed. Nyamsi et al. [25] derived the correlation between entropy generation of the experimental system and the error on activation energy, by which the effects of various design and operating conditions were analyzed. They found the conclusion qualitatively consistent with those reported in the literature, and the predicted relative error also fell into the known range. Recently, a series of Round Robin Test on the hydrogen storage properties of typical materials including MH were organized by the Institute of Energy and Transportation, Joint Research Center of the European Union [29,30]. It was found that the inter laboratory deviation of the kinetic data, say t90, was as high as 60% [30], although the materials tested came from exactly the same batch and identical measurement protocols were used by the participating laboratories. The authors of the work also concluded that thermal effects and improper sample mass/system volume might be the causes of data dispersion. It should be noted however, that besides the quality of the experimental data, the parameter estimation also depends on the theoretical method used. Unfortunately, most works have overlooked this important issue. In this paper, by applying theoretical tools of frequentist analysis and optimal experimental design (OED) for parameter estimation, we obtained more accurate and consistent estimation of the parameters in hydriding/dehydriding kinetics model. We propose a practical procedure involving OED, and we show that this is useful for the analysis of hydriding/dehydriding kinetic measurements.
Theoretical methods Kinetics models and parameter estimation The kinetic models can be classified into two types [31], one is empirical and focuses on the dependence of differential
5178
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 1 ( 2 0 1 6 ) 5 1 7 6 e5 1 8 7
reaction rate dX/dt on the phenomenological factors like reacted fraction X, temperature T and pressure P, while the other is directly derived from the single particle analysis (SPA) based on gasesolid reaction theory. Since the former is favorable for kinetic data fitting, our discussion will be focused on this type of models, which can be written as: dX=dt ¼ F1 ðXÞf2 ðPÞf3 ðTÞ
(1)
Suppose that 1/F1(X) can be analytically integrated with respect to X to get f1(X), then after simple manipulation we are able to obtain the following integral kinetic equation, f1 ðXÞ ¼ f2 ðPÞ$f3 ðTÞ$t
(2)
For the reaction under isothermal and isobaric conditions, the product of f2 and f3 is actually the reaction rate constant. Denote this constant by k, then the relationship is expressed as: k ¼ f2 ðPÞ$f3 ðTÞ
(4)
For the commonly used functions with respect to X and P, i.e. f1 and f2, please refer to Li et al. [20] and Forde et al. [32]. On the other hand, f3(T) generally takes the form of Arrhenius equation as follows: Ea f3 ðTÞ ¼ C$exp Rg $T
of time points ti (i=1, 2,…, I) with a fixed interval Δt Xmeas(ti) Fit Xmeas to Eq. (4) by least squares method to estimate the reaction rate constant k and other parameters (if any) under Pj and Tk k(Pi, Tk) Keep temperature unchanged and vary pressure Pj (j=1, 2,…, J) in a group of experiments, repeat the above two steps and fit the estimated k to Eq.(3) to obtain the constant f3(T) (since T is kept the same) f3(Tk)
(3)
Substitute Equation (3) to Equation (2), we can get the model equation applicable to the reaction rate curve from single kinetic experiment at a given pressure P and temperature T, k$t ¼ f1 ðXÞ
Perform single experiment at a given pressure Pj and temperature Tk to measure reacted fraction Xmeas at a series
(5)
where Rg is the general gas constant, and C and Ea are the preexponential factor and activation energy to be estimated, respectively. Considering the abovementioned effects of different factors on the reaction rate, a general sequential procedure for estimation of the unknown parameters C and Ea can be summarized in the flow diagram shown in Fig. 1. As mentioned in Introduction, for a practical study of hydriding/dehydriding kinetics, model selection is also an important issue, i.e., determining which model among various candidates best describes the physics reflected by the data. For example, we may have several f1 models representing different controlling mechanisms, and one of them should be selected as the “best” model before getting the estimated k. The similar situation applies to f2 models, too. Model selection need to be made according to a well-justified criterion about what is the “best” model. The criterion can be as simple as the commonly-used R-square, however, some more rigorous criterions based on information theory [33], e.g., Akaike information criterion (AIC), Bayesian information criterion (BIC) are frequently used. The procedure shown in Fig. 1 makes such model selection possible because the dependence of reaction rate on different factors (X, P and T) is identified step-by-step. However, the inclusion of model selection will greatly complicates our discussion, hence in this investigation we focus on the parameter estimation itself and assume that we already know the “true” model. In this case, we are able to
Vary temperature Tk (k=1, 2,…, K) in groups of experiments, repeat the above three steps and fit the estimated
f3(T)
to
Eq.(5)
to
determine
the
pre-exponential factor C & activation energy Ea
Fig. 1 e Flow diagram of the conventional procedure for parameter estimation in the hydriding/dehydriding kinetics model.
combine the last two steps in Fig. 1 and fit the estimated k with respect to varied pressure P and temperature T simultaneously using Equation (3).
Parameter estimation and optimal experimental design According to the frequentist optimal experimental design (OED) theory, generally we can have the following expression for the experimental measurement [34], ymeas ¼ gðq0 ; uÞ þ ε
(6)
where ymeas is the experimental data, which can be interpreted as system response to the input variables u, g(.,.) is a parametric model dependent on the “true” value of parameters q0 and ε denotes the measurement error. q may be a vector with multiple elements: q ¼ [q1, q2, …, qN]T. The data ymeas are collected and matched to a model: ymodel ¼ g b q; u
(7)
where b q is an unknown parameter vector. Due to the error ε, the estimated value of the parameters, b q, will inevitably deviate from the true value q0 . Suppose all ε’s are identically and independently distributed, following a Gaussian distribution, ε N 0; s2
(8)
Then asymptotic statistics can be applied. If b q and q0 are reasonably close, the least squares estimate (LSE) of parameter b q is unbiased and follows a Gaussian distribution,
5179
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 1 ( 2 0 1 6 ) 5 1 7 6 e5 1 8 7
b q Nðq0 ; Vq Þ
(9)
q. If b q is a scalar, Vq Where Vq denotes the covariance matrix of b will only be a scalar representing the variance. Obviously, “larger” Vq suggests higher uncertainty, hence should be avoided. The uncertainty of the estimated parameters depends on two factors. One is the error ε on the experimental data, characterized by the standard deviation sε, and the other is the sensitivity of the model with respect to the estimated parameters. Since the latter is dependent on the variables u, the aim of OED is to minimize Vq by properly selecting u without increasing experimental effort. The implementation of OED actually follows an optimization procedure. Firstly, the “tunable” experimental variables u are specified and the constraints of the variables are given. Thereafter, for known forms of model and error structure of the data, the covariance matrix of b q is derived in an asymptotic manner as the inverse of Fisher information matrix F. See the details in Faller et al. [34] and Ciucci et al. [35e37]. Vq zðFÞ1 ¼
!1 M X 1 T ½V gðq ; u Þ ½V gðq ; u Þ q 0 q 0 i i s2 i¼1 ε;i
(10)
where sε,i is the standard deviation of data collected in the ith experiment, characterized by variables ui. M is the number of experiments conducted. Finally, certain measure of the asymptotic covariance matrix Vq by a function 4 is minimized to seek optimal experimental conditions: q u ¼ arg min f Vq u; b u2U
several or even tens of times of hydriding/dehydriding cycles are needed for the hydrogen storage properties to reach equilibrium. Therefore, a single experiment takes a considerably long time (several hours or even a couple of days) and significant effort. As a result, only a limited number of experiments with varied temperatures or pressures (the total number is J K in Fig. 1, and a number larger than 20 is rarely reported in the literature) are allowed for the second problem. This typically implies that the data set for kb is small, making the optimal selection of P and T conditions important for accurate estimation of C and Ea. On the other hand, for the first problem the continuous collection of Xmeas is conducted with small time interval Dt (say 5s, comparable to the apparatus in our laboratory [24]) over a single experiment, thus providing a potentially large data set and facilitating accurate estimation of k. However, uncertainty quantification following OED is still beneficial for the understanding of experimental settings and subsequent estimation of C and Ea, as will be discussed later.
Results and discussions To understand the effects of OED on hydriding/dehydriding kinetics measurement, we used the model and parameters presented by Forde et al. [32] from experiments on La0.83Ce0.10Pr0.04 Nd0.03Ni4.40Al0.60 material as a reference. The authors concluded that the hydriding/dehydriding kinetics can be modeled by the following relationship 1=n
(11)
The most frequently used functions 4 of Vq include the determinant, trace and the maximum eigenvalue of the matrix, termed D-optimality, A-optimality and E-optimality, respectively. Combining the above statement on OED and the flow diagram in Fig. 1, we can see that two OED problems are involved in parameter estimation for hydriding/dehydriding kinetics models. The first one is for a single experiment at constant pressure and temperature, in which Xmeas serves as “y” (6) and experimental time t is the key variable to be optimized, namely “u” (6). The reaction rate constant k and other parameters (if any) form the parameter vector “q” to be estimated. The second one is for groups of experiments under varied pressures and temperatures, in which the estimated reaction rate constant kb is taken as “y”, and the set of experimental conditions (including pressure P and temperature T) is optimized as “u”. The parameter vector for estimation comprises the pre-exponential factor C and the activation energy Ea. It is important to note that the second problem uses estimated parameter rather than true experimental data as “y”, however, the analytic methods of OED can still be borrowed considering the Gaussian error structure of the estimated parameters. If the parametric models g(.,.) (Equations (4) and (3) for the two problems) are known, we can easily find the mathematical expressions for the asymptotic covariance matrix from (10) and use them for the OED analysis. As is well known, the kinetic measurements of MH are generally conducted after activation of the material, i.e.,
k$t ¼ ðlnð1 XÞÞ
(12)
The form of f2(P) and key parameters of the model are shown in Table 1. The Van't Hoff parameters of DH and DS for the hydriding processes, which are needed for the calculation of equilibrium pressures Peq, were also taken from the original paper [32] as 27 kJ/mol and 103.9 J mol1 K1, respectively. DH DS þ ln Peq ¼ Rg T Rg
(13)
Here the following box constraints are given for the experimental variables: 0 < tend 10000s 293K T 333K Pe ðTÞ þ 0:1MPa P 2MPa
(14)
Please note that the difference of reaction pressure and the equilibrium pressure, which is often termed as the driving force of the reaction, is set greater than 0.1 MPa in this investigation so that a reasonable reaction rate can be achieved. Considering the multiple factors affecting the accuracy
Table 1 e Kinetic model and fitted parameters as reference. Hydriding Pressure dependence Reaction order n, Pre-exponential factor C, s1 Activation energy Ea, kJ mol1
f2(P) ¼ ln(P/Peq) 2 798 26
5180
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 1 ( 2 0 1 6 ) 5 1 7 6 e5 1 8 7
of hydrogen content measurement, it is plausible to assume that the error on an individual data follows normal distribution with zero mean. The standard deviation sε for the hydrogen content [H/M] measurement is approximately 5% according to Muthukumar et al. [22] and Qajar et al. [38] under different pressure and temperature conditions. Hence the measured reacted fraction Xmeas, which is defined as the ratio of instantaneous hydrogen content [H/M] to saturated content [H/M]s, has a standard deviation of about 7%.
Uncertainty of estimate for k and n From (12) and the statement in Theoretical methods, we can see parameters k and n are estimated from fitting Xmeas at a series of successive time points. For this problem, the time duration of the kinetic measurement, denoted here by td is the tunable experimental variable, hence was taken for discussion on the effects for parameter estimation. Moreover, one should note that the parameter k depends on T and P (n is pretty fixed). Therefore, the parameter estimation quality was also evaluated for different k0 (the maximum value in the range of pressures and temperatures under discussion is about 0.03s1). Firstly, the objective functions 4 with regard to Vq were calculated under different td and k0, as shown in Fig. 2. Although the magnitudes of the different functions 4D, 4A and 4E (corresponding to D-optimality, A-optimality and E-optimality) vary, their trends are quite consistent. At the beginning the function values decrease as td extends, until a critical value of td dependent on k is reached and no further change will occur, suggesting a similar tendency for parameter estimate uncertainty. This result coincides with our intuition that with the data collection interval Dt fixed, experiment with longer duration generates more data points, which are beneficial to the accurate estimation of parameters. While after a certain threshold, say the time for 99% reaction completion (t99), the reacted fraction hardly changes and can no longer provide useful information. We can conclude that no significant variation in 4(Vq) occur after t99. The variable t99 take the values of 72, 715, 7150s for k0 ¼ 0.03, 0.003, 0.0003s1, respectively. In order to examine the uncertainty of k and n estimate in a more intuitive way, we also made a prediction for the socalled asymptotic confidence region based on Fisher information matrix, in which the parameter values obtained from LSE are believed to fall with more than 99% confidence. For comparison, 1000 Monte Carlo simulations for synthetic experiments [35e37] mimicking the corresponding LSE processes were also conducted to find the real statistical distribution of the parameters. The results in the case of different k0 and td are shown in Figs. 3e5, where dn is defined b . The b n0 , hence dn/n0 reflects the relative error of n as n similar expression also applies for the other estimated parameters, i.e. k, C and Ea.
Fig. 2 e The variation of characteristic function 4(V) for k and n estimation with respect to td when: (a) k0 ¼ 0.03s¡1; (b) k0 ¼ 0.003s¡1; (c) k0 ¼ 0.0003s¡1.
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 1 ( 2 0 1 6 ) 5 1 7 6 e5 1 8 7
5181
Fig. 3 e The asymptotic confidence region and Monte Carlo simulation results for synthetic experiments on estimating k and n when k ¼ 0.03s¡1. The experimental settings are: (a) td ¼ 50s, Dt ¼ 5s (b) td ¼ 100s, Dt ¼ 5s (c) td ¼ 50s, Dt ¼ 5s, conducted twice (d) td ¼ 50s, Dt ¼ 2s.
From the good agreement of the asymptotic confidence region (the ellipse) and the simulation results (the scatter points), it is apparent that the asymptotic covariate well represents the confidence region and is a powerful tool for quantifying parametric uncertainty. In the case of k0 ¼ 0.03s1, the increase of td from 50 to 100s or longer only leads to marginal improvement in the quality of parameter estimation. Under such situation, we may have to repeat the same experiments or increase the time resolution (smaller Dt) so that sufficient confidence can be achieved, as indicated in Fig. 3(c) and (d). When k0 decreases down to 0.003s1, the asymptotic confidence region shows a long and narrow shape for a short experimental duration td of 50s (see Fig. 4(a)), suggesting poor identifiability from the relatively large condition number of Vq [34]. In this case, the uncertainty in k estimate may be larger than 10%, which is generally not acceptable. By extending td to 500s, the uncertainty in both parameters are reduced to less than 2%. Similar to the case of k0 ¼ 0.03s1, extension of td beyond the critical value (see Fig. 4(c)) is meaningless while repeating the experiments contributes to better estimation of k and n (see Fig. 4(d)). When k0 takes a smaller value of 0.0003s1, the general trend with varied td is the same as those in the previous discussion. However, as k0 decreases from 0.03 to 0.0003s1, the parameter estimation quality achieved when td approaches to the b critical value, turns out to be better. The uncertainty in kb and n are reduced by 90%, which should be attributed to longer
reaction time and more useful data points. Another interesting observation is that estimate of k and n show different sensitivities to the variations of both td and k0. Considering that n generally take integer or half integer values [32], the estimation of n, whose uncertainty varies in a narrow range from 0.3% to 7% with different td and k0, is relatively easy to implement. On the other hand, more attention should be paid to the estimation of k since it is more sensitive to k0 and td.
OED for estimation of C, Ea In the conventional setting for hydriding/dehydriding kinetic experiments, uniformly spaced temperatures are often chosen in the concerning range, and regularly spaced pressures are assigned for each temperature [18,39]. Although such setting is useful to test the pressure dependence and the general model applicability, it may not be advisable for highquality parameter estimation. For simplicity, we assume that the relative uncertainty in kb measured by sε, is constant (10%, 5%, 2%, 1%) for different temperatures and pressures, thus the comparison between the conventional and optimal setting can be made. In the conventional setting, the 16 experiments were organized as follows: 4 temperatures within the range of (14) (293, 303, 313, 323 K) are chosen, and then for each temperature, 4 pressures (1.4, 1.6, 1.8, 2.0 MPa), providing sufficient driving force of reaction, are adopted. On the other hand, the following optimization
5182
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 1 ( 2 0 1 6 ) 5 1 7 6 e5 1 8 7
Fig. 4 e The asymptotic confidence region and Monte Carlo simulation results for synthetic experiments on estimating k and n when k ¼ 0.003s¡1. The experimental settings are: (a) td ¼ 50s, Dt ¼ 5s (b) td ¼ 500s, Dt ¼ 5s (c) td ¼ 2000s, Dt ¼ 5s (d) td ¼ 500s, Dt ¼ 5s, conducted twice.
problem (based on D-optimality) was formulated to get the optimal settings of P and T for the estimation of C and Ea, fP; Tgopt ¼ arg min detðVq ðC0 ; Ea0 ; P; TÞÞ fP;Tg
(15)
The solution of the problem was found to be 8 repeated experiments at 293 K, 2 MPa and 8 repeated experiments at 333 K, 2 MPa. The asymptotic confidence region and the results from 1000 Monte Carlo simulations are shown in Fig. 6 for the conventional setting. Due to the strong nonlinearity of the b fitting function, in case of large uncertainty of k(10% or 5%), the predicted confidence region does not cover the scattered points in the figure well. This is because Vq is derived from asymptotic expansion around the “true” value q0, hence it fails to characterize the real error propagation accurately in the region far from this value. As sε decreases, the agreement of asymptotic confidence region and simulation results gets better. Moreover, the estimation quality of both C and Ea also improves. It is noted that the magnitude of the uncertainty for estimating C and Ea can grow to a prohibitively high level, for example, a reasonable sε of 5% results in an uncertainty of 100% and 10% for C and Ea estimation, respectively. For comparison, we present the corresponding results in Fig. 7 for the optimal setting. As is expected, in each case of sε (10%, 5%, 2%, 1%) the optimal setting facilitates more accurate estimation of C and Ea, achieving a reduction of uncertainty by
about 1/3. Undoubtedly such comparison, together with the discussion in next section, validates our point of view that large disparity in the parameter values reported in literature may be partly caused by improper experimental design and data post-processing. It is worth noting that even in the best b is still larger than 10% and, therecase, the uncertainty in C fore, challenges experimental consistency.
Sequential procedure integrating OED methods We should keep in mind that for a nonlinear model g(u, q) with respect to q, the evaluation of asymptotic covariance matrix and implementation of OED procedure depends on the “true” values of the parameters, q0, which can hardly be known. This can be resolved by robust experimental design [36] dealing with a minemax problem, or a sequential procedure doing parameter estimation and OED alternatively. Here we just give a discussion on the latter means. As is well known, least squares estimation of parameters is based on the solution of following problem: b q ¼ arg min q
m X
2 wi ymodel ðq; ui Þ ymeas ðq0 ; ui Þ
(16)
i¼1
Where wi is the weight for ith experimental data. The so-called maximum likelihood estimate (MLE), in which wi is calculated as the reciprocal of the squared error for the
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 1 ( 2 0 1 6 ) 5 1 7 6 e5 1 8 7
5183
Fig. 5 e The asymptotic confidence region and Monte Carlo simulation results for synthetic experiments on estimating k and n when k ¼ 0.0003s¡1. The experimental settings are: (a) td ¼ 50s, Dt ¼ 5s (b) td ¼ 5000s, Dt ¼ 5s (c) td ¼ 8000s, Dt ¼ 5s (d) td ¼ 5000s, Dt ¼ 5s, conducted twice.
corresponding data [34,35], weighs favorably the measurements with small error and proves beneficial to the parameter estimation [35]. Therefore, MLE should be adopted if possible. For the first OED problem, the relative error of each Xmeas data is assumed to be constant (7%), hence the parameter estimation using Equation (16) can be easily implemented. However, for the second OED problem the error on kb can not be obtained unless the uncertainty analysis in Uncertainty of estimate for k and n is conducted. In order to show the respective effects of uncertainty analysis and optimal experimental settings by OED methods, we present 3 different procedures for comparison. Procedure 1: 1. The experimental conditions (P and T) are specified beforehand following conventional setting, i.e. combination of uniformly spaced pressures and temperatures, and all the experiments (the maximum number is given) are performed at a time to collect data Xmeas. Each experiment proceeds until t99 elapses to generate enough data for accurate estimation of parameters k and n. 2. The MLE of k and n under each experimental condition is found by fitting Xmeas with respect to t using (4) and (16). 3. Simple least squares, in which wi in Equation (16) is unitary, are applied to the fitting of kb with respect to P
and T using (3). Thus, the parameters C and Ea are estimated. Procedure 2: 1 & 2. Identical to the 1st and 2nd step in Procedure 1. 3. Uncertainty analysis as described in Uncertainty of estimate for k and n is carried out for each kb under the specified conditions. Subsequently, the weight wi is calculated accordingly. 4. The MLE of C and Ea are obtained by fitting kb with respect to P and T using (3) and (16). Procedure 3: 1e4. Identical to the 1st to 4th steps of procedure 2 except that fewer experiments are conducted at the beginning, and the rest ones are to be completed after their conditions are optimized 5. OED as described in OED for estimation of C, Ea is b and E ba implemented (note that the estimated parameters C from preceding steps need to be used at the place of their “true” values in (15)) to determine the settings of the rest experiments. 6. The rest experiments are carried out under the optimized conditions.
5184
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 1 ( 2 0 1 6 ) 5 1 7 6 e5 1 8 7
Fig. 6 e The asymptotic confidence region and Monte Carlo simulation results for synthetic experiments on estimating C and Ea through conventional setting. Suppose the relative uncertainty in k estimate satisfies: (a) sk ¼ 10% (b) sk ¼ 5% (c) sk ¼ 2% (d)sk ¼ 1%.
7. Steps 2e4 are repeated to the data collected from the additional experiments and the estimation of C and Ea is refined. We should stress that each of the abovementioned 3 procedures can be implemented seamlessly in a practical kinetic study, since the “true” parameters are not needed. For this discussion, the total number of experiments is set to be 16. In Procedure 1 and 2 we firstly chose 4 equally-spaced temperatures covering the concerning range. For each temperature, 4 equally-spaced pressures from the equilibrium value to the maximum value were then used. In Procedure 3, 9 and 7 experiments were assigned before and after optimization of the experimental settings. For the 9 preliminary experiments, a conventional setting similar to that in Procedure 1, except that the numbers of both temperature and pressure levels were 3, was applied. Through 1000 Monte Carlo simulations, the estimation results of C and Ea in synthetic experiments following Procedure 1, 2 and 3 were obtained, as shown in Figs. 8e10. The maximum relative errors in estimated C and Ea are found to be about 40% and 4% for Procedure 1, respectively. On the other hand, for Procedure 2 the corresponding errors are less than 30% and 3%, while the addition of OED in Procedure 3 will further reduce the maximum errors in C and Ea to 20% and 2%. We can conclude from the above comparison that both uncertainty analysis and OED improve the quality of parameter estimation. While model generality
should also be considered when tuning the experimental settings, the uncertainty analysis which enhances the parameter estimation confidence by simply additional data processing, is particularly valuable in the practice of kinetics study. The settings of the 7 additional experiments from OED in Procedure 3 are given in Equation (17). 9 8 297 5:7709 > > > > > > > 297 5:7709 > > > > > > > > > 297 5:7709 = < (17) fP; Tgopt 297 5:7709 > > > > > > > 297 5:7709 > > > > > > 297 5:7709 > > > ; : 303 6:9243 As it is evident from this study, the repetition of experiments under carefully selected conditions as reported in the literature [35e37,40], is beneficial because it is linked to higher sensitivity of data to the estimated parameters. It should also be noted that the optimal pressure barely provides the minimum driving force of reaction (specified as 0.1 MPa here), leading to a decrease in reaction rate constant k0 and a smaller b as revealed in Uncertainty of relative uncertainty for k, estimate for k and n. A similar conclusion was also obtained by Muthukumar et al. [22] regarding isothermal conditions. The optimal experimental setting once again highlights the importance of keeping the uncertainty of the k estimate as low as possible in order to gain lower overall uncertainty.
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 1 ( 2 0 1 6 ) 5 1 7 6 e5 1 8 7
5185
Fig.7 e The asymptotic confidence region and Monte Carlo simulation results for synthetic experiments on estimating C and Ea through optimal setting. Suppose the relative uncertainty in k estimate satisfies: (a) sk ¼ 10% (b) sk ¼ 5% (c) sk ¼ 2% (d) sk ¼ 1%.
Fig. 8 e Predicted relative error in C and Ea estimate by Procedure 1 from 1000 Monte Carlo simulations.
Fig. 9 e Predicted relative error in C and Ea estimate by Procedure 2 from 1000 Monte Carlo simulations.
5186
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 1 ( 2 0 1 6 ) 5 1 7 6 e5 1 8 7
Acknowledgments This work was financially supported by the National Natural Science Foundation of China under grant No. 51506174 and No. 21276209.
references
Fig. 10 e Predicted relative error in C and Ea estimate by Procedure 3 from 1000 Monte Carlo simulations.
Conclusion In this paper, the uncertainty of estimated parameters in a hydriding/dehydriding kinetic model was analyzed with theoretical tools of OED. From the case study based on kinetic data reported by Forde et al. [32], we were able to draw the following conclusion: (1) There exists a critical value for the experiment duration td under isothermal and isobaric conditions. Extending td beyond the value leads to negligible improvement on the estimate of k and n. This threshold can be approximated by 99% completion time t99. (2) The estimate of different parameters (k vs n, C vs Ea) shows quite different trends as a function of the experimental conditions. For the present investigation, the estimate of k is more sensitive to the variation of td and k0 than n. While Ea can be estimated with sufficient confidence using a suitable procedure, the error in C is often too large to permit any meaningful comparison among the results from different research groups. (3) Prohibitively large uncertainty in the estimate of C and Ea may arise if the uncertainty of k estimate is not controlled carefully. This observation also supports our view that the disparity in the parameter values reported in the literature may also be caused by the estimation methods used. (4) The proposed sequential procedure integrating OED is shown to be effective for improving the parameter estimation quality (in our case, the area of asymptotic confidence region shrinks by more than 70%). The results also show that both uncertainty analysis and OED steps can reduce the estimation uncertainty of both C and Ea.
[1] Sakintuna B, Lamari-Darkrim F, Hirscher M. Metal hydride materials for solid hydrogen storage: a review,. Int J Hydrogen Energy 2007;32(9):1121e40. [2] Botzung M, Chaudourne S, Gillia O, Perret C, Latroche M, Percheron-Guegan A, et al. Simulation and experimental validation of a hydrogen storage tank with metal hydrides. Int J Hydrogen Energy 2008;33(1):98e104. [3] Visaria M, Mudawar I, Pourpoint T. Enhanced heat exchanger design for hydrogen storage using high-pressure metal hydride: Part 1. Design methodology and computational results. Int J Heat Mass Transf 2011;54(1e3):413e23. [4] Yang FS, Cao XX, Zhang ZX, Bao ZW, Wu Z, Nyamsi SN. Assessment on the long-term performance of a LaNi5 based hydrogen storage system. Energy Procedia 2012;29:720e30. [5] Wang YQ, Yang FS, Meng XY, Guo QF, Zhang ZX, Park IS, et al. Simulation study on the reaction process based single stage metal hydride thermal compressor. Int J Hydrogen Energy 2010;35:321e8. [6] Hopkins RR, Kim KJ. Hydrogen compression characteristics of a dual stage thermal compressor system utilizing LaNi5 and Ca0.6Mm0.4Ni5 as the working metal hydrides. Int J Hydrogen Energy 2010;35:5693e702. [7] Miura S, Fujisawa A, Ishida M. A hydrogen purification and storage system using metal hydride. Int J Hydrogen Energy 2012;37:2794e9. [8] Dunikov D, Borzenko V, Malyshenko S. Influence of impurities on hydrogen absorption in a metal hydride reactor. Int J Hydrogen Energy 2012;37:13843e8. [9] Chen XY, Wei LX, Deng L, Yang FS, Zhang ZX. A review on the metal hydride based hydrogen purification and separation technology. Appl Mech Mater 2014;448e453:3027e36. B. High temperature metal [10] Felderhoff M, Bogdanovic hydrides as heat storage materials for solar and related applications. Int J Mol Sci 2009;10:325e44. [11] Bao ZW, Yang FS, Wu Z, Nyamsi SN, Zhang ZX. Optimal design of metal hydride reactors based on CFD-Taguchi combined method. Energy Convers Manag 2013;65:322e30. [12] Muthukumar P, Groll M. Metal hydride based heating and cooling systems: a review. Int J Hydrogen Energy 2010:3817e31. [13] Yang FS, Zhang ZX, Wang GX, Bao ZW, Diniz da Costa JC, Rudolph V. Numerical study of a metal hydride heat transformer for low-grade heat recovery. Appl Therm Eng 2011;31:2749e956. [14] Martin M, Gommel C, Borkhart C, Fromm E. Absorption and desorption kinetics of hydrogen storage alloys. J Alloys Compd 1996;238:193e201. [15] Schweppe F, Martin M, Fromm E. Model on hydride formation describing surface control, diffusion control and transition regions. J Alloys Compd 1997;261:254e8. rube V, Radtke G, Dresselhaus M, Chen G. Size effects on [16] Be the hydrogen storage properties of nanostructured metal hydrides: a review. Int J Energy Res 2007;31(6e7):637e63. [17] Bershadsky E, Klyuch A, Ron M. Hydrogen absorption and desorption kinetics of TiFe0.8Ni0.2H. Int J Hydrogen Energy 1995;20(1):29e33.
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 1 ( 2 0 1 6 ) 5 1 7 6 e5 1 8 7
[18] Oh JW, Kim CY, Nahm KS, Sim KS. The hydriding kinetics of LaNi4.5Al0.5 with hydrogen,. J Alloys Compd 1998;278:270e6. [19] Inomata A, Aoki H, Miura T. Measurement and modeling of hydriding and dehydriding kinetics. J Alloys Compd 1998;278:103e9. [20] Li Q, Lin Q, Chou KC, Jiang LJ. A study on the hydridingdehydriding kinetics of Mg1.9Al0.1Ni. J Mater Sci 2004;39:61e5. [21] Yang H, Ojo A, Ogaro P, Goudy AJ. Hydriding and dehydriding kinetics of Sodium Alanate at constant pressure thermodynamic driving forces. J Phys Chem C 2009;113(32):14512e7. [22] Muthukumar P, Satheesh A, Linder M, Mertz R, Groll M. Studies on hydriding kinetics of some La-based metal hydride alloys. Int J Hydrogen Energy 2009;34:7253e62. J, Freni A, Corbera n JM, Compan ~ V. Hydriding kinetics of [23] Paya LaNi5 using nucleation-growth and diffusion models. J New Mater Electrochem Syst 2012;15(4):293e300. [24] Cao XX, Yang FS, Wu Z, Wang YQ, Zhang ZX. Experimental study on the hydrogen storage properties of LaNi5 alloy in repeated hydriding/dehydriding cycles. Adv Mater Res 2013;815:25e30. [25] Nyamsi SN, Yang FS, Wu Z, Bao ZW, Zhang ZX. Assessment of errors on the kinetic data by entropy generation analysis. Int J Hydrogen Energy 2012;37:12365e74. [26] Wang XL, Suda S. Consistent determination of the intrinsic kinetic properties between hydrogen and hydriding alloys. J Alloys Compd 1995;231(1):660e5. [27] Blach TP, Gray EM. Sieverts apparatus and methodology for accurate determination of hydrogen uptake by light-atom hosts. J Alloys Compd 2007;446e447:692e7. [28] Broom DP. The accuracy of hydrogen sorption measurements on potential storage materials. Int J Hydrogen Energy 2007;32:4871e88. [29] Zlotea C, Moretto P, Steriotis T. A Round Robin characterization of the hydrogen sorption properties of a carbon based material. Int J Hydrogen Energy 2009;34:3044e57.
5187
[30] Moretto P, Zlotea C, Dolci F, Amieiro A, Bobet JL, Borgschulte A, et al. A Round Robin Test exercise on hydrogen absorption/desorption properties of a magnesium hydride based material. Int J Hydrogen Energy 2013;38:6704e17. [31] Hua L, Yang FS, Meng XY, Bao ZW, Zhang ZX. Progress in kinetic models for hydrogenation/dehydrogenation of metal hydrides. Chem Ind Eng Progress 2010;29(3):413e8. [32] Forde T, Maehlen JP, Yartys VA, Lototsky MV, Uchida H. Influence of intrinsic hydrogenation/dehydrogenation kinetics on the dynamic behavior of metal hydrides: a semiempirical model and its verification. Int J Hydrogen Energy 2007;32:1041e9. [33] Burnham KP, Anderson DR. Multimodel inference: understanding AIC and BIC in model selection. Sociol Methods Res 2004;33(2):261e304. [34] Faller D, Klingmu¨ller U, Timmer J. Simulation methods for optimal experimental design in system biology. Simulation 2003;79(12):717e25. [35] Ciucci F. Revisiting parameter identification in electrochemical impedance spectroscopy: weighted least squares and optimal experimental design. Electrochimica Acta 2013;87:532e45. [36] Ciucci F, Carraro T, Chueh W, Lai W. Reducing error and measurement time in impedance spectroscopy using model based optimal experimental design. Electrochimica Acta 2011;56:5416e34. [37] Ciucci F. A statistical perspective on oxygen diffusion and surface exchange experiments: sensitivity analysis, parameter estimation and robust optimal experimental design,. Solid State Ion 2013;232:97e105. [38] Qajar A, Peer M, Rajagopalan R, Foley HC. High pressure hydrogen adsorption apparatus: design and error analysis. Int J Hydrogen Energy 2012;37:9123e36. [39] Jemni A, Ben Nasrallah S, Lamloumi J. Experimental and theoretical study of a metal-hydrogen reactor. Int J Hydrogen Energy 1999;24:631e44. [40] Pronzato L. Optimal experimental design and some related control problems. Automatica 2008;44:303e25.
Available online at www.sciencedirect.com
ScienceDirect journal homepage: www.elsevier.com/locate/he
Improving the estimation quality of parameters in kinetic models for hydriding/dehydriding reactions: An OED study Fusheng Yang a,*, Francesco Ciucci b, Zhen Wu a, Zaoxiao Zhang a,c, Yuqi Wang d a
School of Chemical Engineering and Technology, Xi'an Jiaotong University, Xi'an, 710049, China Department of Aerospace and Mechanical Engineering & Department of Chemical and Biomolecular Engineering, Hong Kong University of Science and Technology, Hong Kong c State Key Laboratory of Multiphase Flow in Power Engineering, Xi'an Jiaotong University, Xi'an, 710049, China d School of Chemical Engineering, Northwest University, Xi'an, 710069, China b
article info
abstract
Article history:
Among the hydrogen storage properties of solid materials, e.g., metal hydrides, the reaction
Received 31 October 2015
kinetics is of great importance and is often investigated by fitting reaction rate curves to given
Received in revised form
models. Thanks to such models, key parameters such as activation energies can be found.
14 January 2016
However, it is noteworthy that the parameters thus obtained often show great disparities
Accepted 15 January 2016
even for the same material. In addition to experimental inaccuracies, these disparities may
Available online 6 February 2016
also be attributed to the improper processing of the data and to suboptimal experimental design. In this paper, we present some theoretical tools based on optimal experimental
Keywords:
design (OED) to analyze and optimize the uncertainty on the estimated parameters. The
Metal hydride
uncertainty is derived from asymptotic statistics and it is shown to be closely linked to the
Hydriding/dehydriding kinetics
sensitivity of the experiment to the parameters. We found that by applying a new OED-based
Optimal experimental design
procedure, the estimate of the kinetic parameters can be effectively improved. Copyright © 2016, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved.
Introduction Due to the unique properties of reversible reaction with large amount of hydrogen, metal hydride (MH) materials have shown great application potential and attracted wide attention in recent years. The possible applications of MH include hydrogen storage [1e4], static hydrogen compression [5,6], gas mixture separation/purification [7e9], thermal energy storage [10,11], and refrigeration/heat pumps [12,13].
To realize such practical applications, we need to have an in-depth understanding on the fundamental physical and chemical properties of the MH material, including the hydriding/dehydriding kinetic properties. Like many other gasesolid reactions, the hydriding reaction is heterogeneous and generally thought to proceed in the following steps [14e16]: (1) Physisorption of H2 molecules on the material surface; (2) Dissociation of physisorbed H2 molecules and chemisorption; (3) Penetration of H atoms to the subsurface;
* Corresponding author. Tel.: þ86 29 82668566; fax: þ86 29 82660689. E-mail address: [email protected] (F. Yang). http://dx.doi.org/10.1016/j.ijhydene.2016.01.086 0360-3199/Copyright © 2016, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved.
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 1 ( 2 0 1 6 ) 5 1 7 6 e5 1 8 7
(4) Diffusion of H atoms in the hydride layer either by an interstitial or a vacancy mechanism; (5) Hydride formation at the a/b phase (i.e. solid solution and hydride) interface. On the other hand, the dehydriding process will follow the abovementioned steps in the reverse order. As is well known, if one of the steps is significantly slower than any other steps, it will be considered as rate-controlling step dictating the actual reaction rate. From the reaction mechanism and ratecontrolling step, various kinetic models can be formulated accordingly [14]. By careful processing the experimental data, relevant physical parameters are obtained. Generally, data fitting is conducted for several possible model candidates, and the best one according to a certain criterion will be the “true” model, which is thought to precisely represent the intrinsic reaction behavior of the material under discussion. After that the unknown parameters in the selected model can be found. Obviously, the abovementioned procedure provides insight into the reaction kinetics pertaining to the target material and facilitates the modeling of the corresponding MH-based application system [2,5,13]. For example, Bershadsky et al. [17] observed the kinetic behavior of TiFe0.8Ni0.2 alloy under strictly controlled conditions, and fitted the data against 3 models to obtain the activation energy and pre-exponential factor. However, no single model seemed to represent the data perfectly. Similarly, Oh et al. [18] conducted stepwise volumetric experiments to study the hydriding kinetics of LaNi4.5Al0.5, and found the rate-controlling step to be dissociative chemisorption of hydrogen molecules. The activation energy of the reaction was determined to be 39.8 kJ/mol H2. Inomata et al. [19] carried out an experimental study on the hydriding/dehydriding kinetics of LaNi5, and concluded that the rate-controlling step of hydriding reaction during the initial period is the nucleation and growth of hydride phase, and then changed to hydrogen diffusion. Through isothermal experiments on Mg1.9Al0.1Ni, Li et al. [20] discovered that 3-D diffusion dictates the hydriding/dehydriding kinetics of the material. They also estimated that the activation energies are 52 ± 2 and 48 ± 1 kJ/mol H2 for both reactions, respectively. Yang et al. [21] investigated the reaction kinetics of NaAlH4 in certain decomposition steps. Their results justified a moving boundary reaction model for hydriding kinetics, while the dehydriding reaction showed quite complicated behavior and could not be well fitted to the models they discussed. Muthukumar et al. [22] measured the hydriding kinetics of several La-based alloys including LaNi5, and their results indicated that diffusion of hydrogen atom is the ratecontrolling step for all three alloys in the plateau region. The activation energies of the alloys are mostly around 30 kJ/mol H2. Paya et al. [23] utilized gravimetric apparatus to obtain hydriding kinetics data of LaNi5 alloy. A Johnson-Mehl-Avrami (JMA) model suggesting nucleation and growth mechanism was found able to describe the reaction behavior in the plateau region satisfactorily. Cao et al. [24] investigated the hydriding kinetics of LaNi5 in repeated cycles using volumetric method and fit the data with JMA model. The activation energy of the material was determined to be 19.4 kJ/mol H2. From the literature review, we can see that even for the well-known material LaNi5, the estimated values of the key
5177
parameters in the hydriding/dehydriding kinetics model, e.g., activation energies, show remarkable disparities [25]. For example, quite a few research groups fit the hydrogenation rate curve successfully to a JMA model suggesting nucleationgrowth mechanism, while the activation energy they determined ranged from 19.5 to 36.8 kJ/mol H2. Such disparity can be partly attributed to the poor experimental settings and varied instrumentation/procedure, on which quite a few previous studies were focused. In an early investigation, Wang and Suda [26] summarized the factors affecting the accuracy of hydriding/dehydriding kinetics measurement by the volumetric method, i.e., thermal effects, relationship between gas holder volume and sample mass, history and number of hydriding/dehydriding cycles, surface state of the particles, particle size, and experimental conditions. Some proposals were made for precise kinetic determinations accordingly. Later, Blach and Gray [27] analyzed the error of hydrogen intake measurement by a Sieverts apparatus, and the importance of a large apparatus volume was highlighted. Moreover, some rules-of-thumb to limit the measurement uncertainty were also given by the authors. Broom [28] presented an extensive review on the hydrogen sorption measurements by various methods, i.e., volumetric method, gravimetric method, and temperature-programmed desorption. A number of practical issues related to the measurement accuracy were discussed. Nyamsi et al. [25] derived the correlation between entropy generation of the experimental system and the error on activation energy, by which the effects of various design and operating conditions were analyzed. They found the conclusion qualitatively consistent with those reported in the literature, and the predicted relative error also fell into the known range. Recently, a series of Round Robin Test on the hydrogen storage properties of typical materials including MH were organized by the Institute of Energy and Transportation, Joint Research Center of the European Union [29,30]. It was found that the inter laboratory deviation of the kinetic data, say t90, was as high as 60% [30], although the materials tested came from exactly the same batch and identical measurement protocols were used by the participating laboratories. The authors of the work also concluded that thermal effects and improper sample mass/system volume might be the causes of data dispersion. It should be noted however, that besides the quality of the experimental data, the parameter estimation also depends on the theoretical method used. Unfortunately, most works have overlooked this important issue. In this paper, by applying theoretical tools of frequentist analysis and optimal experimental design (OED) for parameter estimation, we obtained more accurate and consistent estimation of the parameters in hydriding/dehydriding kinetics model. We propose a practical procedure involving OED, and we show that this is useful for the analysis of hydriding/dehydriding kinetic measurements.
Theoretical methods Kinetics models and parameter estimation The kinetic models can be classified into two types [31], one is empirical and focuses on the dependence of differential
5178
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 1 ( 2 0 1 6 ) 5 1 7 6 e5 1 8 7
reaction rate dX/dt on the phenomenological factors like reacted fraction X, temperature T and pressure P, while the other is directly derived from the single particle analysis (SPA) based on gasesolid reaction theory. Since the former is favorable for kinetic data fitting, our discussion will be focused on this type of models, which can be written as: dX=dt ¼ F1 ðXÞf2 ðPÞf3 ðTÞ
(1)
Suppose that 1/F1(X) can be analytically integrated with respect to X to get f1(X), then after simple manipulation we are able to obtain the following integral kinetic equation, f1 ðXÞ ¼ f2 ðPÞ$f3 ðTÞ$t
(2)
For the reaction under isothermal and isobaric conditions, the product of f2 and f3 is actually the reaction rate constant. Denote this constant by k, then the relationship is expressed as: k ¼ f2 ðPÞ$f3 ðTÞ
(4)
For the commonly used functions with respect to X and P, i.e. f1 and f2, please refer to Li et al. [20] and Forde et al. [32]. On the other hand, f3(T) generally takes the form of Arrhenius equation as follows: Ea f3 ðTÞ ¼ C$exp Rg $T
of time points ti (i=1, 2,…, I) with a fixed interval Δt Xmeas(ti) Fit Xmeas to Eq. (4) by least squares method to estimate the reaction rate constant k and other parameters (if any) under Pj and Tk k(Pi, Tk) Keep temperature unchanged and vary pressure Pj (j=1, 2,…, J) in a group of experiments, repeat the above two steps and fit the estimated k to Eq.(3) to obtain the constant f3(T) (since T is kept the same) f3(Tk)
(3)
Substitute Equation (3) to Equation (2), we can get the model equation applicable to the reaction rate curve from single kinetic experiment at a given pressure P and temperature T, k$t ¼ f1 ðXÞ
Perform single experiment at a given pressure Pj and temperature Tk to measure reacted fraction Xmeas at a series
(5)
where Rg is the general gas constant, and C and Ea are the preexponential factor and activation energy to be estimated, respectively. Considering the abovementioned effects of different factors on the reaction rate, a general sequential procedure for estimation of the unknown parameters C and Ea can be summarized in the flow diagram shown in Fig. 1. As mentioned in Introduction, for a practical study of hydriding/dehydriding kinetics, model selection is also an important issue, i.e., determining which model among various candidates best describes the physics reflected by the data. For example, we may have several f1 models representing different controlling mechanisms, and one of them should be selected as the “best” model before getting the estimated k. The similar situation applies to f2 models, too. Model selection need to be made according to a well-justified criterion about what is the “best” model. The criterion can be as simple as the commonly-used R-square, however, some more rigorous criterions based on information theory [33], e.g., Akaike information criterion (AIC), Bayesian information criterion (BIC) are frequently used. The procedure shown in Fig. 1 makes such model selection possible because the dependence of reaction rate on different factors (X, P and T) is identified step-by-step. However, the inclusion of model selection will greatly complicates our discussion, hence in this investigation we focus on the parameter estimation itself and assume that we already know the “true” model. In this case, we are able to
Vary temperature Tk (k=1, 2,…, K) in groups of experiments, repeat the above three steps and fit the estimated
f3(T)
to
Eq.(5)
to
determine
the
pre-exponential factor C & activation energy Ea
Fig. 1 e Flow diagram of the conventional procedure for parameter estimation in the hydriding/dehydriding kinetics model.
combine the last two steps in Fig. 1 and fit the estimated k with respect to varied pressure P and temperature T simultaneously using Equation (3).
Parameter estimation and optimal experimental design According to the frequentist optimal experimental design (OED) theory, generally we can have the following expression for the experimental measurement [34], ymeas ¼ gðq0 ; uÞ þ ε
(6)
where ymeas is the experimental data, which can be interpreted as system response to the input variables u, g(.,.) is a parametric model dependent on the “true” value of parameters q0 and ε denotes the measurement error. q may be a vector with multiple elements: q ¼ [q1, q2, …, qN]T. The data ymeas are collected and matched to a model: ymodel ¼ g b q; u
(7)
where b q is an unknown parameter vector. Due to the error ε, the estimated value of the parameters, b q, will inevitably deviate from the true value q0 . Suppose all ε’s are identically and independently distributed, following a Gaussian distribution, ε N 0; s2
(8)
Then asymptotic statistics can be applied. If b q and q0 are reasonably close, the least squares estimate (LSE) of parameter b q is unbiased and follows a Gaussian distribution,
5179
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 1 ( 2 0 1 6 ) 5 1 7 6 e5 1 8 7
b q Nðq0 ; Vq Þ
(9)
q. If b q is a scalar, Vq Where Vq denotes the covariance matrix of b will only be a scalar representing the variance. Obviously, “larger” Vq suggests higher uncertainty, hence should be avoided. The uncertainty of the estimated parameters depends on two factors. One is the error ε on the experimental data, characterized by the standard deviation sε, and the other is the sensitivity of the model with respect to the estimated parameters. Since the latter is dependent on the variables u, the aim of OED is to minimize Vq by properly selecting u without increasing experimental effort. The implementation of OED actually follows an optimization procedure. Firstly, the “tunable” experimental variables u are specified and the constraints of the variables are given. Thereafter, for known forms of model and error structure of the data, the covariance matrix of b q is derived in an asymptotic manner as the inverse of Fisher information matrix F. See the details in Faller et al. [34] and Ciucci et al. [35e37]. Vq zðFÞ1 ¼
!1 M X 1 T ½V gðq ; u Þ ½V gðq ; u Þ q 0 q 0 i i s2 i¼1 ε;i
(10)
where sε,i is the standard deviation of data collected in the ith experiment, characterized by variables ui. M is the number of experiments conducted. Finally, certain measure of the asymptotic covariance matrix Vq by a function 4 is minimized to seek optimal experimental conditions: q u ¼ arg min f Vq u; b u2U
several or even tens of times of hydriding/dehydriding cycles are needed for the hydrogen storage properties to reach equilibrium. Therefore, a single experiment takes a considerably long time (several hours or even a couple of days) and significant effort. As a result, only a limited number of experiments with varied temperatures or pressures (the total number is J K in Fig. 1, and a number larger than 20 is rarely reported in the literature) are allowed for the second problem. This typically implies that the data set for kb is small, making the optimal selection of P and T conditions important for accurate estimation of C and Ea. On the other hand, for the first problem the continuous collection of Xmeas is conducted with small time interval Dt (say 5s, comparable to the apparatus in our laboratory [24]) over a single experiment, thus providing a potentially large data set and facilitating accurate estimation of k. However, uncertainty quantification following OED is still beneficial for the understanding of experimental settings and subsequent estimation of C and Ea, as will be discussed later.
Results and discussions To understand the effects of OED on hydriding/dehydriding kinetics measurement, we used the model and parameters presented by Forde et al. [32] from experiments on La0.83Ce0.10Pr0.04 Nd0.03Ni4.40Al0.60 material as a reference. The authors concluded that the hydriding/dehydriding kinetics can be modeled by the following relationship 1=n
(11)
The most frequently used functions 4 of Vq include the determinant, trace and the maximum eigenvalue of the matrix, termed D-optimality, A-optimality and E-optimality, respectively. Combining the above statement on OED and the flow diagram in Fig. 1, we can see that two OED problems are involved in parameter estimation for hydriding/dehydriding kinetics models. The first one is for a single experiment at constant pressure and temperature, in which Xmeas serves as “y” (6) and experimental time t is the key variable to be optimized, namely “u” (6). The reaction rate constant k and other parameters (if any) form the parameter vector “q” to be estimated. The second one is for groups of experiments under varied pressures and temperatures, in which the estimated reaction rate constant kb is taken as “y”, and the set of experimental conditions (including pressure P and temperature T) is optimized as “u”. The parameter vector for estimation comprises the pre-exponential factor C and the activation energy Ea. It is important to note that the second problem uses estimated parameter rather than true experimental data as “y”, however, the analytic methods of OED can still be borrowed considering the Gaussian error structure of the estimated parameters. If the parametric models g(.,.) (Equations (4) and (3) for the two problems) are known, we can easily find the mathematical expressions for the asymptotic covariance matrix from (10) and use them for the OED analysis. As is well known, the kinetic measurements of MH are generally conducted after activation of the material, i.e.,
k$t ¼ ðlnð1 XÞÞ
(12)
The form of f2(P) and key parameters of the model are shown in Table 1. The Van't Hoff parameters of DH and DS for the hydriding processes, which are needed for the calculation of equilibrium pressures Peq, were also taken from the original paper [32] as 27 kJ/mol and 103.9 J mol1 K1, respectively. DH DS þ ln Peq ¼ Rg T Rg
(13)
Here the following box constraints are given for the experimental variables: 0 < tend 10000s 293K T 333K Pe ðTÞ þ 0:1MPa P 2MPa
(14)
Please note that the difference of reaction pressure and the equilibrium pressure, which is often termed as the driving force of the reaction, is set greater than 0.1 MPa in this investigation so that a reasonable reaction rate can be achieved. Considering the multiple factors affecting the accuracy
Table 1 e Kinetic model and fitted parameters as reference. Hydriding Pressure dependence Reaction order n, Pre-exponential factor C, s1 Activation energy Ea, kJ mol1
f2(P) ¼ ln(P/Peq) 2 798 26
5180
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 1 ( 2 0 1 6 ) 5 1 7 6 e5 1 8 7
of hydrogen content measurement, it is plausible to assume that the error on an individual data follows normal distribution with zero mean. The standard deviation sε for the hydrogen content [H/M] measurement is approximately 5% according to Muthukumar et al. [22] and Qajar et al. [38] under different pressure and temperature conditions. Hence the measured reacted fraction Xmeas, which is defined as the ratio of instantaneous hydrogen content [H/M] to saturated content [H/M]s, has a standard deviation of about 7%.
Uncertainty of estimate for k and n From (12) and the statement in Theoretical methods, we can see parameters k and n are estimated from fitting Xmeas at a series of successive time points. For this problem, the time duration of the kinetic measurement, denoted here by td is the tunable experimental variable, hence was taken for discussion on the effects for parameter estimation. Moreover, one should note that the parameter k depends on T and P (n is pretty fixed). Therefore, the parameter estimation quality was also evaluated for different k0 (the maximum value in the range of pressures and temperatures under discussion is about 0.03s1). Firstly, the objective functions 4 with regard to Vq were calculated under different td and k0, as shown in Fig. 2. Although the magnitudes of the different functions 4D, 4A and 4E (corresponding to D-optimality, A-optimality and E-optimality) vary, their trends are quite consistent. At the beginning the function values decrease as td extends, until a critical value of td dependent on k is reached and no further change will occur, suggesting a similar tendency for parameter estimate uncertainty. This result coincides with our intuition that with the data collection interval Dt fixed, experiment with longer duration generates more data points, which are beneficial to the accurate estimation of parameters. While after a certain threshold, say the time for 99% reaction completion (t99), the reacted fraction hardly changes and can no longer provide useful information. We can conclude that no significant variation in 4(Vq) occur after t99. The variable t99 take the values of 72, 715, 7150s for k0 ¼ 0.03, 0.003, 0.0003s1, respectively. In order to examine the uncertainty of k and n estimate in a more intuitive way, we also made a prediction for the socalled asymptotic confidence region based on Fisher information matrix, in which the parameter values obtained from LSE are believed to fall with more than 99% confidence. For comparison, 1000 Monte Carlo simulations for synthetic experiments [35e37] mimicking the corresponding LSE processes were also conducted to find the real statistical distribution of the parameters. The results in the case of different k0 and td are shown in Figs. 3e5, where dn is defined b . The b n0 , hence dn/n0 reflects the relative error of n as n similar expression also applies for the other estimated parameters, i.e. k, C and Ea.
Fig. 2 e The variation of characteristic function 4(V) for k and n estimation with respect to td when: (a) k0 ¼ 0.03s¡1; (b) k0 ¼ 0.003s¡1; (c) k0 ¼ 0.0003s¡1.
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 1 ( 2 0 1 6 ) 5 1 7 6 e5 1 8 7
5181
Fig. 3 e The asymptotic confidence region and Monte Carlo simulation results for synthetic experiments on estimating k and n when k ¼ 0.03s¡1. The experimental settings are: (a) td ¼ 50s, Dt ¼ 5s (b) td ¼ 100s, Dt ¼ 5s (c) td ¼ 50s, Dt ¼ 5s, conducted twice (d) td ¼ 50s, Dt ¼ 2s.
From the good agreement of the asymptotic confidence region (the ellipse) and the simulation results (the scatter points), it is apparent that the asymptotic covariate well represents the confidence region and is a powerful tool for quantifying parametric uncertainty. In the case of k0 ¼ 0.03s1, the increase of td from 50 to 100s or longer only leads to marginal improvement in the quality of parameter estimation. Under such situation, we may have to repeat the same experiments or increase the time resolution (smaller Dt) so that sufficient confidence can be achieved, as indicated in Fig. 3(c) and (d). When k0 decreases down to 0.003s1, the asymptotic confidence region shows a long and narrow shape for a short experimental duration td of 50s (see Fig. 4(a)), suggesting poor identifiability from the relatively large condition number of Vq [34]. In this case, the uncertainty in k estimate may be larger than 10%, which is generally not acceptable. By extending td to 500s, the uncertainty in both parameters are reduced to less than 2%. Similar to the case of k0 ¼ 0.03s1, extension of td beyond the critical value (see Fig. 4(c)) is meaningless while repeating the experiments contributes to better estimation of k and n (see Fig. 4(d)). When k0 takes a smaller value of 0.0003s1, the general trend with varied td is the same as those in the previous discussion. However, as k0 decreases from 0.03 to 0.0003s1, the parameter estimation quality achieved when td approaches to the b critical value, turns out to be better. The uncertainty in kb and n are reduced by 90%, which should be attributed to longer
reaction time and more useful data points. Another interesting observation is that estimate of k and n show different sensitivities to the variations of both td and k0. Considering that n generally take integer or half integer values [32], the estimation of n, whose uncertainty varies in a narrow range from 0.3% to 7% with different td and k0, is relatively easy to implement. On the other hand, more attention should be paid to the estimation of k since it is more sensitive to k0 and td.
OED for estimation of C, Ea In the conventional setting for hydriding/dehydriding kinetic experiments, uniformly spaced temperatures are often chosen in the concerning range, and regularly spaced pressures are assigned for each temperature [18,39]. Although such setting is useful to test the pressure dependence and the general model applicability, it may not be advisable for highquality parameter estimation. For simplicity, we assume that the relative uncertainty in kb measured by sε, is constant (10%, 5%, 2%, 1%) for different temperatures and pressures, thus the comparison between the conventional and optimal setting can be made. In the conventional setting, the 16 experiments were organized as follows: 4 temperatures within the range of (14) (293, 303, 313, 323 K) are chosen, and then for each temperature, 4 pressures (1.4, 1.6, 1.8, 2.0 MPa), providing sufficient driving force of reaction, are adopted. On the other hand, the following optimization
5182
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 1 ( 2 0 1 6 ) 5 1 7 6 e5 1 8 7
Fig. 4 e The asymptotic confidence region and Monte Carlo simulation results for synthetic experiments on estimating k and n when k ¼ 0.003s¡1. The experimental settings are: (a) td ¼ 50s, Dt ¼ 5s (b) td ¼ 500s, Dt ¼ 5s (c) td ¼ 2000s, Dt ¼ 5s (d) td ¼ 500s, Dt ¼ 5s, conducted twice.
problem (based on D-optimality) was formulated to get the optimal settings of P and T for the estimation of C and Ea, fP; Tgopt ¼ arg min detðVq ðC0 ; Ea0 ; P; TÞÞ fP;Tg
(15)
The solution of the problem was found to be 8 repeated experiments at 293 K, 2 MPa and 8 repeated experiments at 333 K, 2 MPa. The asymptotic confidence region and the results from 1000 Monte Carlo simulations are shown in Fig. 6 for the conventional setting. Due to the strong nonlinearity of the b fitting function, in case of large uncertainty of k(10% or 5%), the predicted confidence region does not cover the scattered points in the figure well. This is because Vq is derived from asymptotic expansion around the “true” value q0, hence it fails to characterize the real error propagation accurately in the region far from this value. As sε decreases, the agreement of asymptotic confidence region and simulation results gets better. Moreover, the estimation quality of both C and Ea also improves. It is noted that the magnitude of the uncertainty for estimating C and Ea can grow to a prohibitively high level, for example, a reasonable sε of 5% results in an uncertainty of 100% and 10% for C and Ea estimation, respectively. For comparison, we present the corresponding results in Fig. 7 for the optimal setting. As is expected, in each case of sε (10%, 5%, 2%, 1%) the optimal setting facilitates more accurate estimation of C and Ea, achieving a reduction of uncertainty by
about 1/3. Undoubtedly such comparison, together with the discussion in next section, validates our point of view that large disparity in the parameter values reported in literature may be partly caused by improper experimental design and data post-processing. It is worth noting that even in the best b is still larger than 10% and, therecase, the uncertainty in C fore, challenges experimental consistency.
Sequential procedure integrating OED methods We should keep in mind that for a nonlinear model g(u, q) with respect to q, the evaluation of asymptotic covariance matrix and implementation of OED procedure depends on the “true” values of the parameters, q0, which can hardly be known. This can be resolved by robust experimental design [36] dealing with a minemax problem, or a sequential procedure doing parameter estimation and OED alternatively. Here we just give a discussion on the latter means. As is well known, least squares estimation of parameters is based on the solution of following problem: b q ¼ arg min q
m X
2 wi ymodel ðq; ui Þ ymeas ðq0 ; ui Þ
(16)
i¼1
Where wi is the weight for ith experimental data. The so-called maximum likelihood estimate (MLE), in which wi is calculated as the reciprocal of the squared error for the
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 1 ( 2 0 1 6 ) 5 1 7 6 e5 1 8 7
5183
Fig. 5 e The asymptotic confidence region and Monte Carlo simulation results for synthetic experiments on estimating k and n when k ¼ 0.0003s¡1. The experimental settings are: (a) td ¼ 50s, Dt ¼ 5s (b) td ¼ 5000s, Dt ¼ 5s (c) td ¼ 8000s, Dt ¼ 5s (d) td ¼ 5000s, Dt ¼ 5s, conducted twice.
corresponding data [34,35], weighs favorably the measurements with small error and proves beneficial to the parameter estimation [35]. Therefore, MLE should be adopted if possible. For the first OED problem, the relative error of each Xmeas data is assumed to be constant (7%), hence the parameter estimation using Equation (16) can be easily implemented. However, for the second OED problem the error on kb can not be obtained unless the uncertainty analysis in Uncertainty of estimate for k and n is conducted. In order to show the respective effects of uncertainty analysis and optimal experimental settings by OED methods, we present 3 different procedures for comparison. Procedure 1: 1. The experimental conditions (P and T) are specified beforehand following conventional setting, i.e. combination of uniformly spaced pressures and temperatures, and all the experiments (the maximum number is given) are performed at a time to collect data Xmeas. Each experiment proceeds until t99 elapses to generate enough data for accurate estimation of parameters k and n. 2. The MLE of k and n under each experimental condition is found by fitting Xmeas with respect to t using (4) and (16). 3. Simple least squares, in which wi in Equation (16) is unitary, are applied to the fitting of kb with respect to P
and T using (3). Thus, the parameters C and Ea are estimated. Procedure 2: 1 & 2. Identical to the 1st and 2nd step in Procedure 1. 3. Uncertainty analysis as described in Uncertainty of estimate for k and n is carried out for each kb under the specified conditions. Subsequently, the weight wi is calculated accordingly. 4. The MLE of C and Ea are obtained by fitting kb with respect to P and T using (3) and (16). Procedure 3: 1e4. Identical to the 1st to 4th steps of procedure 2 except that fewer experiments are conducted at the beginning, and the rest ones are to be completed after their conditions are optimized 5. OED as described in OED for estimation of C, Ea is b and E ba implemented (note that the estimated parameters C from preceding steps need to be used at the place of their “true” values in (15)) to determine the settings of the rest experiments. 6. The rest experiments are carried out under the optimized conditions.
5184
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 1 ( 2 0 1 6 ) 5 1 7 6 e5 1 8 7
Fig. 6 e The asymptotic confidence region and Monte Carlo simulation results for synthetic experiments on estimating C and Ea through conventional setting. Suppose the relative uncertainty in k estimate satisfies: (a) sk ¼ 10% (b) sk ¼ 5% (c) sk ¼ 2% (d)sk ¼ 1%.
7. Steps 2e4 are repeated to the data collected from the additional experiments and the estimation of C and Ea is refined. We should stress that each of the abovementioned 3 procedures can be implemented seamlessly in a practical kinetic study, since the “true” parameters are not needed. For this discussion, the total number of experiments is set to be 16. In Procedure 1 and 2 we firstly chose 4 equally-spaced temperatures covering the concerning range. For each temperature, 4 equally-spaced pressures from the equilibrium value to the maximum value were then used. In Procedure 3, 9 and 7 experiments were assigned before and after optimization of the experimental settings. For the 9 preliminary experiments, a conventional setting similar to that in Procedure 1, except that the numbers of both temperature and pressure levels were 3, was applied. Through 1000 Monte Carlo simulations, the estimation results of C and Ea in synthetic experiments following Procedure 1, 2 and 3 were obtained, as shown in Figs. 8e10. The maximum relative errors in estimated C and Ea are found to be about 40% and 4% for Procedure 1, respectively. On the other hand, for Procedure 2 the corresponding errors are less than 30% and 3%, while the addition of OED in Procedure 3 will further reduce the maximum errors in C and Ea to 20% and 2%. We can conclude from the above comparison that both uncertainty analysis and OED improve the quality of parameter estimation. While model generality
should also be considered when tuning the experimental settings, the uncertainty analysis which enhances the parameter estimation confidence by simply additional data processing, is particularly valuable in the practice of kinetics study. The settings of the 7 additional experiments from OED in Procedure 3 are given in Equation (17). 9 8 297 5:7709 > > > > > > > 297 5:7709 > > > > > > > > > 297 5:7709 = < (17) fP; Tgopt 297 5:7709 > > > > > > > 297 5:7709 > > > > > > 297 5:7709 > > > ; : 303 6:9243 As it is evident from this study, the repetition of experiments under carefully selected conditions as reported in the literature [35e37,40], is beneficial because it is linked to higher sensitivity of data to the estimated parameters. It should also be noted that the optimal pressure barely provides the minimum driving force of reaction (specified as 0.1 MPa here), leading to a decrease in reaction rate constant k0 and a smaller b as revealed in Uncertainty of relative uncertainty for k, estimate for k and n. A similar conclusion was also obtained by Muthukumar et al. [22] regarding isothermal conditions. The optimal experimental setting once again highlights the importance of keeping the uncertainty of the k estimate as low as possible in order to gain lower overall uncertainty.
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 1 ( 2 0 1 6 ) 5 1 7 6 e5 1 8 7
5185
Fig.7 e The asymptotic confidence region and Monte Carlo simulation results for synthetic experiments on estimating C and Ea through optimal setting. Suppose the relative uncertainty in k estimate satisfies: (a) sk ¼ 10% (b) sk ¼ 5% (c) sk ¼ 2% (d) sk ¼ 1%.
Fig. 8 e Predicted relative error in C and Ea estimate by Procedure 1 from 1000 Monte Carlo simulations.
Fig. 9 e Predicted relative error in C and Ea estimate by Procedure 2 from 1000 Monte Carlo simulations.
5186
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 1 ( 2 0 1 6 ) 5 1 7 6 e5 1 8 7
Acknowledgments This work was financially supported by the National Natural Science Foundation of China under grant No. 51506174 and No. 21276209.
references
Fig. 10 e Predicted relative error in C and Ea estimate by Procedure 3 from 1000 Monte Carlo simulations.
Conclusion In this paper, the uncertainty of estimated parameters in a hydriding/dehydriding kinetic model was analyzed with theoretical tools of OED. From the case study based on kinetic data reported by Forde et al. [32], we were able to draw the following conclusion: (1) There exists a critical value for the experiment duration td under isothermal and isobaric conditions. Extending td beyond the value leads to negligible improvement on the estimate of k and n. This threshold can be approximated by 99% completion time t99. (2) The estimate of different parameters (k vs n, C vs Ea) shows quite different trends as a function of the experimental conditions. For the present investigation, the estimate of k is more sensitive to the variation of td and k0 than n. While Ea can be estimated with sufficient confidence using a suitable procedure, the error in C is often too large to permit any meaningful comparison among the results from different research groups. (3) Prohibitively large uncertainty in the estimate of C and Ea may arise if the uncertainty of k estimate is not controlled carefully. This observation also supports our view that the disparity in the parameter values reported in the literature may also be caused by the estimation methods used. (4) The proposed sequential procedure integrating OED is shown to be effective for improving the parameter estimation quality (in our case, the area of asymptotic confidence region shrinks by more than 70%). The results also show that both uncertainty analysis and OED steps can reduce the estimation uncertainty of both C and Ea.
[1] Sakintuna B, Lamari-Darkrim F, Hirscher M. Metal hydride materials for solid hydrogen storage: a review,. Int J Hydrogen Energy 2007;32(9):1121e40. [2] Botzung M, Chaudourne S, Gillia O, Perret C, Latroche M, Percheron-Guegan A, et al. Simulation and experimental validation of a hydrogen storage tank with metal hydrides. Int J Hydrogen Energy 2008;33(1):98e104. [3] Visaria M, Mudawar I, Pourpoint T. Enhanced heat exchanger design for hydrogen storage using high-pressure metal hydride: Part 1. Design methodology and computational results. Int J Heat Mass Transf 2011;54(1e3):413e23. [4] Yang FS, Cao XX, Zhang ZX, Bao ZW, Wu Z, Nyamsi SN. Assessment on the long-term performance of a LaNi5 based hydrogen storage system. Energy Procedia 2012;29:720e30. [5] Wang YQ, Yang FS, Meng XY, Guo QF, Zhang ZX, Park IS, et al. Simulation study on the reaction process based single stage metal hydride thermal compressor. Int J Hydrogen Energy 2010;35:321e8. [6] Hopkins RR, Kim KJ. Hydrogen compression characteristics of a dual stage thermal compressor system utilizing LaNi5 and Ca0.6Mm0.4Ni5 as the working metal hydrides. Int J Hydrogen Energy 2010;35:5693e702. [7] Miura S, Fujisawa A, Ishida M. A hydrogen purification and storage system using metal hydride. Int J Hydrogen Energy 2012;37:2794e9. [8] Dunikov D, Borzenko V, Malyshenko S. Influence of impurities on hydrogen absorption in a metal hydride reactor. Int J Hydrogen Energy 2012;37:13843e8. [9] Chen XY, Wei LX, Deng L, Yang FS, Zhang ZX. A review on the metal hydride based hydrogen purification and separation technology. Appl Mech Mater 2014;448e453:3027e36. B. High temperature metal [10] Felderhoff M, Bogdanovic hydrides as heat storage materials for solar and related applications. Int J Mol Sci 2009;10:325e44. [11] Bao ZW, Yang FS, Wu Z, Nyamsi SN, Zhang ZX. Optimal design of metal hydride reactors based on CFD-Taguchi combined method. Energy Convers Manag 2013;65:322e30. [12] Muthukumar P, Groll M. Metal hydride based heating and cooling systems: a review. Int J Hydrogen Energy 2010:3817e31. [13] Yang FS, Zhang ZX, Wang GX, Bao ZW, Diniz da Costa JC, Rudolph V. Numerical study of a metal hydride heat transformer for low-grade heat recovery. Appl Therm Eng 2011;31:2749e956. [14] Martin M, Gommel C, Borkhart C, Fromm E. Absorption and desorption kinetics of hydrogen storage alloys. J Alloys Compd 1996;238:193e201. [15] Schweppe F, Martin M, Fromm E. Model on hydride formation describing surface control, diffusion control and transition regions. J Alloys Compd 1997;261:254e8. rube V, Radtke G, Dresselhaus M, Chen G. Size effects on [16] Be the hydrogen storage properties of nanostructured metal hydrides: a review. Int J Energy Res 2007;31(6e7):637e63. [17] Bershadsky E, Klyuch A, Ron M. Hydrogen absorption and desorption kinetics of TiFe0.8Ni0.2H. Int J Hydrogen Energy 1995;20(1):29e33.
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 1 ( 2 0 1 6 ) 5 1 7 6 e5 1 8 7
[18] Oh JW, Kim CY, Nahm KS, Sim KS. The hydriding kinetics of LaNi4.5Al0.5 with hydrogen,. J Alloys Compd 1998;278:270e6. [19] Inomata A, Aoki H, Miura T. Measurement and modeling of hydriding and dehydriding kinetics. J Alloys Compd 1998;278:103e9. [20] Li Q, Lin Q, Chou KC, Jiang LJ. A study on the hydridingdehydriding kinetics of Mg1.9Al0.1Ni. J Mater Sci 2004;39:61e5. [21] Yang H, Ojo A, Ogaro P, Goudy AJ. Hydriding and dehydriding kinetics of Sodium Alanate at constant pressure thermodynamic driving forces. J Phys Chem C 2009;113(32):14512e7. [22] Muthukumar P, Satheesh A, Linder M, Mertz R, Groll M. Studies on hydriding kinetics of some La-based metal hydride alloys. Int J Hydrogen Energy 2009;34:7253e62. J, Freni A, Corbera n JM, Compan ~ V. Hydriding kinetics of [23] Paya LaNi5 using nucleation-growth and diffusion models. J New Mater Electrochem Syst 2012;15(4):293e300. [24] Cao XX, Yang FS, Wu Z, Wang YQ, Zhang ZX. Experimental study on the hydrogen storage properties of LaNi5 alloy in repeated hydriding/dehydriding cycles. Adv Mater Res 2013;815:25e30. [25] Nyamsi SN, Yang FS, Wu Z, Bao ZW, Zhang ZX. Assessment of errors on the kinetic data by entropy generation analysis. Int J Hydrogen Energy 2012;37:12365e74. [26] Wang XL, Suda S. Consistent determination of the intrinsic kinetic properties between hydrogen and hydriding alloys. J Alloys Compd 1995;231(1):660e5. [27] Blach TP, Gray EM. Sieverts apparatus and methodology for accurate determination of hydrogen uptake by light-atom hosts. J Alloys Compd 2007;446e447:692e7. [28] Broom DP. The accuracy of hydrogen sorption measurements on potential storage materials. Int J Hydrogen Energy 2007;32:4871e88. [29] Zlotea C, Moretto P, Steriotis T. A Round Robin characterization of the hydrogen sorption properties of a carbon based material. Int J Hydrogen Energy 2009;34:3044e57.
5187
[30] Moretto P, Zlotea C, Dolci F, Amieiro A, Bobet JL, Borgschulte A, et al. A Round Robin Test exercise on hydrogen absorption/desorption properties of a magnesium hydride based material. Int J Hydrogen Energy 2013;38:6704e17. [31] Hua L, Yang FS, Meng XY, Bao ZW, Zhang ZX. Progress in kinetic models for hydrogenation/dehydrogenation of metal hydrides. Chem Ind Eng Progress 2010;29(3):413e8. [32] Forde T, Maehlen JP, Yartys VA, Lototsky MV, Uchida H. Influence of intrinsic hydrogenation/dehydrogenation kinetics on the dynamic behavior of metal hydrides: a semiempirical model and its verification. Int J Hydrogen Energy 2007;32:1041e9. [33] Burnham KP, Anderson DR. Multimodel inference: understanding AIC and BIC in model selection. Sociol Methods Res 2004;33(2):261e304. [34] Faller D, Klingmu¨ller U, Timmer J. Simulation methods for optimal experimental design in system biology. Simulation 2003;79(12):717e25. [35] Ciucci F. Revisiting parameter identification in electrochemical impedance spectroscopy: weighted least squares and optimal experimental design. Electrochimica Acta 2013;87:532e45. [36] Ciucci F, Carraro T, Chueh W, Lai W. Reducing error and measurement time in impedance spectroscopy using model based optimal experimental design. Electrochimica Acta 2011;56:5416e34. [37] Ciucci F. A statistical perspective on oxygen diffusion and surface exchange experiments: sensitivity analysis, parameter estimation and robust optimal experimental design,. Solid State Ion 2013;232:97e105. [38] Qajar A, Peer M, Rajagopalan R, Foley HC. High pressure hydrogen adsorption apparatus: design and error analysis. Int J Hydrogen Energy 2012;37:9123e36. [39] Jemni A, Ben Nasrallah S, Lamloumi J. Experimental and theoretical study of a metal-hydrogen reactor. Int J Hydrogen Energy 1999;24:631e44. [40] Pronzato L. Optimal experimental design and some related control problems. Automatica 2008;44:303e25.