A new probabilistic tool for the determination and optimization of multiphase equation of state parameters : Application to tin

A thermodynamically consistent Equation Of State (EOS) was developed to predict and analyse the behaviour of multiphase metals under shock wave loading. Assuming the Mie-Gruneisen hypothesis together with the Birch (for example) formulation, the EOS gives the relation between pressure P, temperature T and atomic volume V. Experimental data (P,V,T) for each phase are provided mainly by X-ray diffraction measurements with diamond anvil cells. In this work, mathematical tools are designed to optimize the determination of the EOS parameters and evaluate uncertainty. The general EOS form is y = fθ(x) where y = P, x = (V,T) and θ is the parameter vector to calibrate. Using experimental data (xi, yi), the least square (non-linear) regression provides an optimal value θ∗ for the fit parameters. The measurement errors on y and x give biased estimation of θ∗ with the standard method. Assuming centered and known variance laws for the errors, a statistical procedure is proposed to estimate θ∗ and determine confidence intervals. Thanks to a Bayesian approach it is possible to introduce physical interval knowledge of the parameters in this procedure. Moreover, various EOS fθ∗ formulations are evaluated with a chi-squared type statistical test. The present method is applied on experimental data for multi phase tin (β and γ phases and liquid state) in order to provide an optimized multi-phase model. Furthermore, the method is used to design further experimental campaign and to evaluate the gain of new experimental data with the corresponding estimated errors.


Introduction
In shock waves experiments on multiphase materials, the diagnostics generally measure only the free surface velocities.The pressure P, specific volume V and temperature T need being inferred from hypotheses and calculations.In order to reproduce these velocities and understand the behaviour of metals under these specific loading conditions, hydrocodes and models (in particular Equation of State -EOS) are developed.In previous works, thermodynamically based EOS (Fig. 1) were developed and validated for tin [1][2][3].However, input EOS parameters gathered in these works are not always consistent as they originate from various experimental works and / or are determined from different EOS formalisms.
In situ x-ray diffraction experiments in Diamond Anvil Cells (DAC) provide Pressure-Volume-Temperature (P, V, T ) data sets that can lead, through an adjustment process, to the determination of an equilibrium EOS [4][5][6][7][8][9].This adjustment is classically obtained by a least square regression on the (P, V, T ) data with a P = f (V, T ) model.With the Mie-Grüneisen hypothesis for each phase the pressure model has the form: The first term P k (V) describes a reference isothermal (temperature T 0 ) pressure, for instance with the Birch a e-mail: emmanuel.fraizier@cea.frformulation [10]: where G 0 = Γ 0 C v 0 with C v 0 is the specific heat constant, Γ 0 is the Grüneisen coefficient, K 0 is the isothermal modulus and N 0 is its first pressure derivative.Concerning the units, T is in Kelvin, P is in GPa, V is in mm 3 kg −1 , and G 0 is in Jkg −1 K −1 .N 0 does not have any unit.(T 0 , P 0 ) defines the reference condition.It must be noted that experiments under extreme conditions of pressure and temperature are often difficult to perform and therefore data on the material of interest is often scarce, if not lacking.It is interesting, before planning a new experiment in order to increase the accuracy of a multiphase EOS, to evaluate the potential benefit of this experiment, taking into account the foreseen data accuracy and the amount of data to be recorded.For this purpose, some new mathematical tools are developed.These tools should give new information on Birch, Murnaghan [11], Vinet [12] and some others different [13][14][15] formulations for each phase.Furthermore, taking into account the phase transitions determined from experimental data and from the thermodynamically based EOS we plan to optimize each parameter (for each phase).This paper deals with the preliminary steps of this work.The methodology and the mathematical tools principles are detailed.The first results from literature data on tin [5,6] present the determination of the EOS parameters for EPJ Web of Conferences β and γ phases.The maximum uncertainties are used for every (P,V,T) data [5,6] but, later, the uncertainties should be defined for each point.

Methodology description
This work presents mathematical tools designed for optimizing the determination of the EOS parameters and evaluate the associated errors.This paragraph gives a first overview of the mathematical methodology.The complete mathematical methodology is postponed at the end of the document (see part 5).The general EOS form is Y = f ϑ (X) where the output is y = P, the input vector is x = (V, T ), and ϑ is the vector of parameters to be adjusted to best fit the experimental data (x i , y i ), i = 1, . . ., n.This adjustment is classically obtained by a least square regression when the input vector x has no measurement error and the output quantity y may have measurement error.Indeed, one can prove that the parameters vector vartheta * n which minimizes the quadratic misfit, which is the sum of the squared residuals (y i − f ϑ (x i )) 2 , is a consistent estimator of ϑ, which means that ϑ * n tends to the real ϑ when n tends to infinity.The problem is more difficult when the input vector x has measurement error.In this paper, we propose an estimate of ϑ that takes into account the measurement errors of both x and y and that has a reduced bias compared to the standard least square estimator.This estimator ϑ n is the ϑ which minimizes K(ϑ) (see part 5).
From this expression, we can study the estimation error distribution and we can deduce a confidence interval for ϑ, which is a direct effect of the measurement errors on x and y.Finally, we study the possibility of taking an a priori knowledge into account.Then we propose a new estimator of ϑ , whose expression is given in part 5.
Moreover, a goodness of fit test is proposed to answer to the question whether the EOS formulation is acceptable or not.A statistical P-value is computed.This quantity follows a uniform law on [0, 1] if the EOS formulation is right and is close to 0 if it is not the case.
These mathematical tools are applied to calibrate the studied EOS for tin from available experimental data [5,6].

Application
We study here the EOS given by the Birch model assuming the Mie-Gruneisen hypothesis defined by ( 1) and ( 2).The parameters (K 0 , N 0 , G 0 ) should be estimated given some set of experimental data.A R program has been developed (www.r-project.org) to estimate the parameters using the proposed methodology.The experimental data correspond to non isotherm experiments for β [5], and γ phases [6].

β phase
In [5], the β phase is characterized by n = 31 observations where the temperature ranges form 298 to 529 K and the pressure ranges form 0 to 8.5 GPa.The reference state is defined by (P 0 , T 0 ) = (0, 300).The standard deviations for the measurement errors are σ P = 0.13, σ V = 125, and σ T = 2.5.The a priori knowledge [1] on each parameter is introduced using a Gaussian law N(µ, σ 2 ) defined in Table 1.
The set of parameters is estimated using the raw data without introducing any a priori knowledge on the parameters and then using the a priori.The estimation results are presented in Table 2.
Based on the mathematical analysis, the error distribution on the parameters is computed with the a priori (cf.paragraph 5.4).Table 3 shows the approximate correlation matrix of the estimation errors, using (4) of section 5.5.We observe strong correlations for some parameters as for (K 0 , N 0 ) and (K 0 , V 0 ).
We have conducted a Monte Carlo study to estimate the correlation matrix in parallel (Table 4).
It should be stressed that both approaches bring similar results.The approximation given by 4 of section 5.5 seems .00 -0.28 -0.15 -0.66 K 0 -0.28 1.00 0.50 0.38 N 0 -0.15 0.50 1.00 -0.26 G 0 -0.66 0.38 -0.26 1.00 to be justified.The P-value to test the physical model equals 0.016.The goodness of the EOS formulation is not rejected at level 1%.

γ phase
In [5], the γ phase is characterized by n = 38 available observations where the temperature ranges form 298 K to 713 K and the pressure ranges form 4.49 to 16.15 GPa.The reference state is defined by (P 0 , T 0 ) = (9.40,300).The standard deviations for the measurement errors are σ P = 0.24, σ V = 116, and σ T = 4.2.
The a priori knowledge on each parameter is introduced using a Gaussian law N(µ, σ 2 ) defined in Table 5.The set of parameters is estimated using the raw data without introducing any a priori knowledge on the parameters and then using the a priori.The estimation results are presented in Table 6.
Table 7 shows the estimated correlation matrix.The respective standard deviation are S V 0 = 93.39,S K 0 = 1.61, S N 0 = 0.54, S G 0 = 37.24.As for the β phase, we have conducted a Monte Carlo study to estimate the correlation matrix in parallel (Table 8).
We observe similar correlation matrices computed with the mathematical approach or with Monte-Carlo.The Pvalue to test the physical model equals 0.027.The goodness of the EOS formulation is not rejected at level 1%.

Isothermal Birch regressions for β and γ phases
The Table 9 summarizes the parameters from literature [1] previously used for β phase and the present results.Then, the isothermal Birch regressions for data tin at 298 K and 498 K are computed and represented on Fig. 2 with the corresponding experimental data from [1].The Table 10 summarizes the parameters from literature [5] previously used for γ phase and the present results.
Then, the isothermal Birch regressions for γ phase at 298 K and 628 K are computed and represented on Fig. 3 with the corresponding experimental data from [6].
The present work regression fits the experimental data with a better accuracy than previously achieved with literature parameters, in particular for phase.It may be explained by our dedicated optimization for these particular 04024-p.3This work 623K Fig. 3. Isothermal Birch regression for γ tin and experimental data [6].

EPJ Web of Conferences
data [5,6] and Birch formulation whereas the literature parameters integrate several experimental data.

Conclusion, benefits and perspectives
A first step of a new tool development to determine the multiphase EOS parameters is reached.The resulting parameters enable to compute isothermal Birch regression in good agreement with experimental data.We can now integrate more experimental data and test different P k (V) formulations.Furthermore, the literature data can be used as an a priori on data.This tool will be used to determine the benefits of new experimental campaign or new diagnostics including the uncertainties.The matrix correlation of the V 0 , K 0 , N 0 and G 0 parameters is estimated and gives further information.A Monte Carlo study, conducted in parallel, shows similar results for this matrix and validates our approach.Furthermore, a P-value quantifies the goodness of physical model.Then, this new tool gives the opportunity to test different P k (V) formulations.The last perspective is to optimize every parameters of each phase, taking into account the experimental transitions and our thermodynamically based EOS model.

Statistical model
Consider two physical variables x = (x 1 , . . ., x d ) ∈ R d and y ∈ R linked by the following equation y = φ(x) where φ is a function from R d to R.An equation of state formulation is defined by a family of functions on R d : where ϑ = (ϑ 1 , ϑ 2 , . . ., ϑ k ) corresponds to the physical parameters.The EOS formulation is relevant if ϕ is mostly equal to f ϑ * for a given value ϑ * ∈ Θ. Considering n physical points (x 1 , y 1 ), . . ., (x n , y n ), an approximation of ϑ * is given by ϑ If the EOS formulation is relevant one can expect ϑ * n = ϑ * for n ≥ k, as in the Birch formulation for instance.
The physical points (x i , y i ) are only known through experimental measurements (X i , Y i ).The measurement errors are modeled as follows: ) is a Gaussian random vector with mean x i and covariance matrix σ 2 X (typically a diagonal matrix).
-Y i ∼ N(y i , σ 2 Y ).-The random vectors X i , Y i , 1 ≤ i ≤ n are supposed to be independent.

Estimation error distribution
Notation: For ϑ = (ϑ 1 , . . ., ϑ k ) and x = (x 1 , . . ., x d ), define the column vectors 1≤i, j≤k and the d × k matrix If M is a matrix, M denotes its transpose.We introduce the assumption: A 2 : A 1 holds true and ϑ * n belongs to the interior of Θ. Moreover the Hessian matrix of κ at ϑ * n is positive definite.Note that the Hessian matrix of κ at ϑ * n is equal to where: σ 2 X denotes the largest eigenvalue of σ 2

X
-Z is a k-dimensional random vector whose distribution converges to N(0, I k ) where I k denotes the k×k identity matrix.

04024-p.4
DYMAT 2012 When the second term of ( 3) is negligible, the first term provides the probability distribution of the estimation error, namely N(0, Σ/n) with Σ = Λ −1 ΓΛ −1 , which order of magnitude is less than n −1/2 (||σ 2 X || 1/2 + σ Y ).Remark 1 For a large number n of observations and when σ 2 X is known, the standard least square contrast K(ϑ) can be successfully replaced by which provides an estimation of ϑ * n with a reduced bias, still of order O(||σ 2 X ||), when f ϑ (x) is a bilinear function of ϑ and x.

Confidence interval
We construct confidence intervals for the component ϑ * n,i of ϑ * n and more generally for u ϑ * n where u is a vector of Moreover u Σu can be approximated by u Σu where Σ is the empirical version of Σ, replacing (x i , y i , ϑ * n ) by (X i , Y i , θ).Finally we obtain that This gives directly confidence intervals for u θn .

Use of a priori knowledge
We propose here a method which takes benefit of a previous estimation of ϑ * n .Assume that we have some information on the first p components of ϑ , that is we observe a pdimensional random vector τ, independent of θ, following the law N(τ * n , Σ) where τ * n = (ϑ * n,1 , ϑ * n,2 , . . ., ϑ * n,p ) and Σ is a given p × p covariance matrix (positive definite).
The new estimator θn we can compute is the "best" linear combination of θ and τ: Note that θ is always as good as θ since (n Σ−1 + Γ) −1 ≤ Σ/n.

Feedback and test on the physical model
Can we detect whether the EOS formulation is good or not?In other words are the y i − f ϑ * n (x i ), 1 ≤ i ≤ n, significantly not equal to zero?
Introduce the log-likelihood contrast: for ϑ ∈ Θ and n vectors ξ 1 , . . ., ξ n ∈ R d set Therefore, if t denotes the computed value of D, the probability P(χ 2 n−k > t) is an approximate P-value for testing if the EOS formulation is relevant (the null hypothesis here).

Table 2 .
β phase, parameter estimation using raw data and then with a priori knowledge.

Table 6 .
γ phase, parameter estimation using raw data and with a priori knowledge.

Table 9 .
Literature and current work parameters for β tin.

Table 10 .
Literature and current work parameters for γ tin.