The Inverse Problem of Identification of Hydrogen Permeability Model

One of the technological challenges for hydrogen materials science is the currently active search for structural materials with important applications (including the ITER project and gas-separation plants). One had to estimate the parameters of diffusion and sorption to numerically model the different scenarios and experimental conditions of the material usage (including extreme ones). The article presents boundary value problems of hydrogen permeability and thermal desorption with dynamical boundary conditions. A numerical method is developed for TDS spectrum simulation, where only integration of a nonlinear system of low order ordinary differential equations is required. The main final output of the article is a noise-resistant algorithm for solving the inverse problem of parametric identification for the aggregated experiment where desorption and diffusion are dynamically interrelated (without the artificial division of studies into the diffusion limited regime (DLR) and the surface limited regime (SLR)).


Introduction
Studies on the interaction of hydrogen isotopes with structural materials are mainly necessitated by problems in the energy industry, metal protection from hydrogen corrosion, and the design of chemical reactors [1][2][3][4][5][6][7][8][9][10].The limiting factors are diffusion processes as well as physical and chemical phenomena at the surface.The transfer parameters depend on the technological features of material batch production.It is therefore unreasonable to target at "tabular data".Instead, effective algorithms for solving the inverse problems of parametric identification of adequate mathematical models by experimental curves (data) are necessary.In this study, we consider the permeability model taking into account the main factors and the self-descriptiveness of the experiment.We shall focus on the methods of permeability and thermal desorption taking into account only the main limiting factors for the applied membrane filtering problem and the informative capabilities of the considered experimental methods.The mathematical research is based on the articles [11][12][13], which provide descriptions of the experimental techniques and experimental material on promising alloys for hydrogen separation and purification.
One had to estimate the parameters of diffusion and sorption to numerically model the different scenarios and experimental conditions of the material usage (including extreme ones) and identify the limiting factors.Some particular problems of the hydrogen materials science related to the topic of this study were presented and investigated in [14][15][16][17][18][19][20][21][22][23].This work is a continuation of [24][25][26][27], where the results of modelling hydrogen thermal desorption under various limiting factors are presented.This article deals with the inverse problem of parametric identification based on the suggested "cascade" experiment.
Experimental practices usually employ various modifications of the penetration method and TDS.The results of measurements depend both on the unit design features and on the procedure of preparing samples for hydrogen permeability testing.A successive use of various methods often causes, for example, impurities to appear on the sample surface, which significantly affects the reproducibility of the results.These data are the input for the inverse problems of parametric identification, which are sensitive to the level of different errors.It is therefore advisable to aggregate experiments to improve the accuracy and informative value of the 2 Advances in Mathematical Physics measurements.We suggest the following set-up of the "cascade" experiment.
A membrane heated to a fixed temperature served as the partition in the vacuum chamber.Degassing was performed in advance.A sufficiently high pressure of hydrogen gas was built up in steps at the inlet side.The penetrating flux was determined by mass-spectrometry in the vacuum maintained at the outlet side.This is the penetration method.Its advantage is a reliable determination of the diffusion coefficient by the Daines-Berrer method (based on the so-called lag time).It allows distinguishing between the bulk and the surface processes in the model, keeping in mind that surface parameters are significantly more difficult to estimate.When the steady state level of the penetrating flux is registered, we increase the inlet pressure and wait until a new steady state value is established.Using (at least) three pressure jumps at the inlet side we record the steady state flux values at the outlet side, thus determining "the degree of rectilinearity" of the isotherm.Then the pumping for vacuum building is stopped and the experiment proceeds as the "communicating vessels" method.When pressure values become nearly equal (the sample is almost uniformly saturated with hydrogen) it is possible to turn off the heating, create the vacuum at both sides of the membrane, and begin slowly reheating the sample (TDS-experiment).In addition, there is no depressurization of the diffusion cell and the sample surface remains uncontaminated by additional impurities.We will clarify the details of the aggregated experiment stages as we describe the method of solving the inverse problem of parametric identification.An important consideration is the uniqueness of the parameter estimates of the investigated model.Mathematicians are often reproached for "fascination with uniqueness theorems".But after all, in justifying the choice of, for example, structural materials for the ITER project, the results obtained on thin laboratory samples are extrapolated to "walls".Uniqueness allows for a correct recomputation.
Papers [18][19][20]27] were dedicated to interpretation of TDS peaks.Analysis of the causes of various TDS surges is crucial for the selection of reactor structural materials contacting with hydrogen isotopes.A sufficiently detailed review is presented in [19,20].Where modeling does not take the dynamics of surface processes into account, TDS peaks are inevitably interpreted as a result of capture of diffusing atomic hydrogen by the structural defects (traps) of the material with different binding energies and (or) as a manifestation of multichannel diffusion [28].The above listing of causes is not exhaustive.In this article, modelling shows that the appearance of a desorption peak can be caused by some combinations of the rates of surface and diffusion processes.This complicates the problem of TDS spectrum interpretation even more.
The main result of the paper is the method of parametric identification from experimental data.The difficulties of inverse problems solving in mathematical physics are well known.There is extensive mathematical literature and a number of specialized journals (inverse problems, ill-posed problems, etc.).In experimental practice, the inverse problem of multiparameter estimation is reduced to the one-factorat-a-time method for DLR and SLR.In real life, however, a material is used in the presence of a dynamic "surface-bulk" interplay.Thoroughly elaborated techniques are available for estimation of the diffusion coefficient.The determination of desorption and dissolution parameters is far more complex (unless the temperature is artificially lowered to "turn off" processes in the bulk).The paper presents an algorithm allowing the estimation of desorption and dissolution coefficients where diffusion and surface processes interact intensively.

Hydrogen Permeability Model
2.1.Distributed Transfer Model.Let us briefly describe the experiment.A sample of a structural material preheated to a fixed temperature acts as a vacuum chamber barrier.The sample degassing is performed in advance.At the initial time moment, pressure is built up at the inlet side by puffing of a portion of molecular hydrogen.The declining pressure in the input chamber and increasing pressure in the output chamber are measured.
Consider hydrogen transfer through the sample (ℓ is the plate thickness and  is its area).The temperature  is constant throughout the experiment.The concentration of dissolved diffusing hydrogen (in monatomic state) is sufficiently low and the diffusion flux can be considered proportional to the concentration gradient.The membrane is thin and the material has a sufficiently high hydrogen permeability coefficient, so we restrict ourselves to a standard diffusion equation: where (, ) is the concentration of diffusing hydrogen.The diffusion coefficient  depends on the sample temperature  in an Arrhenius way with preexponential factor  0 and activation energy   :  =  0 exp{−  /[()]}.Initial data are determined by the fact that the sample had been preliminarily degassed: (0, ) = 0, (0, ) = 0,  ∈ [0, ℓ].
Nonlinear boundary conditions are derived from the material flux balance: Here,  in (),  out () are the amounts of hydrogen atoms in the input chamber of volume  in and output chamber of volume  out ,  0 () ≡ (, 0),  ℓ () ≡ (, ℓ).The identity sign is frequently used here in the sense of equality by definition.Within the considered operating temperature range the gaseous hydrogen is in molecular form, but for consistency (considering that atomic hydrogen diffuses through the metal membrane) we use atoms as the unit.According to the kinetic gas theory, the incident particle flux density   is related to the pressure  by the Hertz Knudsen formula:   = / √ 2 ( is the Boltzmann constant,  is the mass of a hydrogen molecule).In the context of the experiment it is convenient to choose the following measurement units [ℓ] = cm, [] = Torr.Then we numerically obtain the dependence   = , () ≈ 2.474 ⋅ 10 22 / √  ([] = 1 H 2 /(Torr cm 2 s), [] = K).The processes of physical adsorption, chemisorption, dissociation of molecules into atoms, and dissolution take place on the surface.Only a small part of "incident" H atoms will, however, be absorbed into the membrane volume.This is taken into account by the factor 2.One can write  (as a parameter of the model) instead of 2, but it is more convenient to interpret the dimensionless probability factor  as the fraction of absorbed H atoms within the 2 notation.Thus, 2 is the resulting flux of atoms through the surface into the bulk without differentiation into more elementary stages.We will omit the word "density" assuming that the surface has unit area.Hereinafter,  0,ℓ =  2 0,ℓ are the densities of the desorption flux from the sample (deviation from the square desorption law is significant only at extreme temperatures) and  is the desorption coefficient.We also assume that  and  depend on the temperature in an Arrhenius way.Formally, the activation energy   in the exponent can as well be negative, being a linear combination of the activation energies and heats of the surface processes on the way "from gas to the solution".If a constant saturation pressure of molecular hydrogen  s = const is maintained at a constant temperature  = const on both sides of the membrane, the equilibrium concentration  of the dissolved atomic diffusionally mobile hydrogen is finally established.By equating the derivatives in model (2) to zero, we get Sieverts' law  ∝ √ s :  = Γ√ s , Γ ≡ √2/.
Let us clarify the experimental conditions.The volumes  in,out comprise several liters, the thickness of the membrane is ℓ < mm, the area  is about 1 cm 2 , and the inflow pressure  0 (0) is within several hundreds Torr.
It now remains to find the magnitudes of  in ,  out .Within the time of transfer through the membrane the gas is in the thermodynamical quasi-equilibrium, wherefore we use the formula  = /().Here,  is the number of gas particles occupying the volume  at the temperature  and the pressure  (in the SI system [] = Pa, [] = m 3 , [] = J/K).Taking into account the relationships Torr = 133.322Pa and Pa = J/m 3 (formally), we get the following expressions for the corresponding pressures and volumes in the boundary conditions (2):  = 2 = /,  ≈ 1.931 ⋅ 10 19 .Here, , ,  are the numerical values in the selected units (Torr, cm 3 , K).
Within the experimental unit the membrane is situated in a tube (which is heated to a predetermined temperature) between the inlet and the outlet chambers.The tube diameter is large enough to consider the equality of pressures as the criterion of thermodynamical quasi-equilibrium between the gas in the tube and in the chambers.The membrane temperature should be taken for the formula for the kinetic constant ().The gas inside the volumes  in,out (whose massive walls are at room temperature) may get heated up.During the preliminary experiment it is recommended to fill the chambers with a practically impermeable metal membrane between them with "ambient" gas.Then we heat the tube and record the pressure rises.Within the framework of the ideal gas approximation (equation of state) this procedure enables estimation of the increments of the gas temperature inside the chambers.The corresponding gas temperatures are the ones to be used in the formula given for  (and the subsequent ones, only excluding the value ).
The need of such a refinement arises from the characteristics of this particular experimental unit.Such an adjustment of the values of  should not cause difficulties in further calculations.Besides, this procedure has relatively little effect on the final calculation of the model pressures taking into account the measurement errors and relatively large volumes .

Fast Hydrogen Permeability Model.
It is clear from physical considerations that a quasi-stationary state is quickly established when the membrane is thin and the material has a sufficiently high hydrogen permeability coefficient: the diffusant concentration distribution is practically linear with respect to the thickness.In this sense, the results of numerical modelling based on the "general" model (the presented boundary-value problem) confirm its adequacy.Since near-to-surface concentrations cannot be measured, the Richardson approximation is usually used in practice to analyze the penetrating flux: Let us formulate the problem of modelling the concentrations  0,ℓ using the pressures  0,ℓ (the problem is also of interest per se) without the quasi-equilibrium simplification () = Γ√().The quasi-stationary state is achieved within a time  0 , which is short compared to the total experiment time (   = −[ 0 () −  ℓ ()]/ℓ).So, the original model ( 1)-( 2) can be simplified (taking into account the formula  = /): ≥  0 > 0. ( Since by virtue of the "inlet-outlet" balance the equalities hold true, it is sufficient to express the concentrations  0,ℓ () =  0,ℓ ( 0 ()) from the boundary conditions (5) and substitute them into the first equation of (4) (the sign is chosen depending on whether the index is 0 or ℓ).The dimensionless normalized variables are convenient for numerical simulation: In addition, the system of equations ( 5) is compactly written in the symmetric form  0 + 2 ℓ =  2 0 ,  ℓ + 2 0 =  2 ℓ .For the variable  ≡  ℓ we obtain the incomplete quartic equation which can be solved in radicals (for physical considerations we are interested in the positive root).However, the explicit expression is somewhat cumbersome and we will anyway have to numerically integrate the first equation of (4) in the form ṗ 0 = ( 0 ).Therefore, we shall aim to derive differential equations for  0,ℓ , since information about the dynamics of the boundary concentrations  0,ℓ is also important.
Differentiate (5) with respect to time and substitute the pressure derivatives from (4).For the variables  0,ℓ we get the system The change of variables in (7) determines the concentrations  0,ℓ () from which the model pressures  0,ℓ (),  ≥  0 are calculated using (5).
(1) We fix  =  0 : omit fast transient processes (the duration of the transient processes is about tens of seconds on the hours-long experimental time scale).For the variable  ≡  ℓ we choose the root of the biquadratic polynomial (2) The system of equations  0 + 2 ℓ =  2 0 ,  ℓ + 2 0 =  2 ℓ ( =  0 ) yields the missing value of  0 ( 0 ).Formally, one equation is enough, but we take into account averaging procedures including determination of the values of  0,ℓ ( 0 ).
Computational experiments show that the model curves almost coincide (at  ≥  0 ) with those generated by the originally proposed model, i.e., the nonlinear distributed initial boundary value problem.
Observe the difference from the quasi-equilibrium model (the Richardson approximation), where the only parameter for approximation of the experimental pressure is the complex Φ = Γ.All the variable parameters of the original model that influence the permeability, namely, , , , are important when running the above algorithm.Thus, the fast hydrogen permeability model does not lose in informativeness concerning the considered transfer parameters.

Numerical Modelling of the Penetration Experiment.
The proposed model is adapted to the experimental conditions and the data range for alloys based on V group metals with high hydrogen permeability, in particular, data for vanadium alloys which are presented in [11][12][13][29][30][31][32][33].We fix  = 673 K, ℓ = 0.05 cm, () = 2.474 × 10 22 / √  1 H 2 /(Torr cm 2 s),  = 2 × 10 −5 cm 2 /s, Γ = 2 × 10 20 1 H /(cm 3 √ Torr), Φ = Γ = 4 × 10 15 1 H /(cm s √ Torr).We set the value  = 1.2 × 10 −4 and calculate the corresponding desorption coefficient  = 2/Γ 2 = 5.7 × 10 −24 cm 4 /s.The input pressures  1,2,3 = {30, 50, 70} Torr were built up in steps at the inlet and maintained to achieve steady state fluxes at the outlet.Degassing of the membrane was done in advance and continuous vacuum pumping was performed at the outlet side.The gas temperature inside the inlet and outlet chambers (whose volumes are sufficiently high) is assumed to be equal to 300 K.This slight difference from the room temperature is due to the heating of the diffusion cell with the sample inside (specified by the characteristics of the equipment).The experimental conditions are such that the concentration at the membrane outlet side is near zero and at the inlet side a stationary concentration is quickly established (but it is lower than the equilibrium one): c < .Within the model we determine   and c by the formulas For the given values of the parameters we have c1 = 1.06 × 10 Next we express the solutions of the standard boundary value problems with Dirichlet boundary conditions corresponding to the jumps of the inlet pressure.Stage I.The boundary value problem of the penetration method is The penetrating flux is From the computational point of view it is convenient to introduce the dimensionless time   = /ℓ 2 which is oriented at the characteristic diffusion time ℓ 2 /.As small  → 0 a singularity appears if we directly use the partial sum of the expression Let us provide another expression for  using the properties of the Jacobi theta function.More precisely, we are interested in the function We have an alternative presentation for  = 0 [34]: The series on the left is rapidly converging for large .But the series on the right is rapidly converging for small .
, which is known as the functional equation for the theta function [35].After some auxiliary transformations for 0 <   ≪ 1 we get the appropriate representation The function (  ) has an -shaped form of the saturation curve (see the inset in Figure 1).Stage I ends with ( * , ) = c1 (ℓ − )ℓ −1 .

Converting Pressure into Flux.
During the real experiment, the gas pressure inside the outlet volume  =  out is measured, not the penetrating flux.Therefore, in this subsection we will provide the corresponding recalculation formula.Denote the pumping rate of the vacuum system by ] ([]] = m 3 /s).We take as the framework the ideal gas state equation:  = .Here, [] = Pa, [] = m 3 ,  is the number of particles (H 2 molecules), and  is the Boltzmann constant.Differentiating on  we get ṗ  = Ṅ.Let us calculate the particle balance over the time Δ: where ).The factor 0.5 is due to the fact that the diffusion flux is atomic, and the particle in the volume  is the H 2 molecule.Let us divide the equation of the material balance by Δ and direct Δ → 0. Finally, we get the differential equation: Denoting  1 = /(2),  0 = /], we get ṗ =  1  ℓ − / 0 ( 0 = 0), wherefrom The measured function () is noised, so we first apply a smoothing procedure and only then calculate the derivative ṗ ().Note that when the pumping system is powerful enough and permeability is relatively slow the first summand acts as a minor correction to the approximation  ℓ ≈ 2]()/().Asymptotically, a steady state value is established: Only the operation of integration is needed to calculate the lag time (see below): so a preliminary approximation of the derivative ṗ () is unnecessary.It suffices to use the composite Simpson formula.One should remember that we use the following measurement units within this subsection In the following, we return to the measurement units accepted in this article.

Boundary Value Problem of Hydrogen Permeability.
When the steady state permeability value is established during the penetration experiment, continuous pumping at the outlet and maintenance of constant pressure at the inlet are stopped.The aggregated experiment moves to the stage of "communicating vessels": inlet pressure declines and outlet pressure grows ( 0,ℓ () are measured).Considering the new time zero we also have  ≥ 0. We are so far talking about the direct problem of modelling hydrogen pressures inside the volumes  in,out .We specify the values:  = 0.5 cm 2 ,  in = 1500 cm 3 ,  out = 2200 cm 3 ,  0 (0) =  0 =  3 ,  = Γ√ 0 .We target here the experimental conditions and the data for the V 85 Ni 15 alloy [11,33].We numerically solve the initial boundary problem of hydrogen permeability: The membrane temperature is taken in the dependences of , , ,  on , and the temperature of the gas inside the chambers is taken in the expressions for  in,out (take into account the correction to room temperature due to the heating of the diffusion cell).The model curves of molecular hydrogen pressures are presented in Figure 2. If we use the ODE system (8), (9) instead of the "full" model, standard software packages will suffice (we substitute the values of the hydrogen temperature inside the volumes  in,out into the expressions for  0,ℓ ).To this end one should skip the initial time  0 within several minutes until a quasi-stationary (not quasi-equilibrium) mode is established.Then the above algorithm is applied to the fast hydrogen permeability model.Figure 2 shows that there is no significant additional error during the numerical simulation (the curves visually coincide).Model curves were numerically generated to test the following algorithm for solving the inverse problem of parametric identification.The parameters generating these curves were then "forgotten".The intersection of the asymptote with the  axis gives the so-called lag time  0 = ℓ 2 /(6), which allows estimating the diffusion coefficient.Analytically,

General Identification Technique
Note that the value under the integral sign is a relative magnitude which does not require absolute values of the penetrating flux in any measurement units ( 1 = sup  1 ()).In addition, the value  0 does not depend on c1 .It is usually assumed that the locally equilibrium concentration  1 = Γ√ 1 is quickly established at the inlet, so one can additionally estimate the solubility Γ = √2/ and the permeability Φ = Γ using the corresponding value  1 =  1 /ℓ.This assumption is not used in this article.We weaken it to increase the accuracy of further estimations.We assume that according to the experimental conditions the stationary inlet concentration  0 () ≈ c1 <  1 ( ≥ ,  ≪ 1) is quickly established given that  ℓ () ≈ 0.
The value of c1 as such is yet to be clarified.Thus, only the estimate of the diffusion coefficient  is considered reliable at this stage.
With a new zero time reference, integrating the expression  2 () we get where   = c  /ℓ,  ≥  * =  * 2 .Formally, we obtain the same expressions for the lag time and the estimate of  if we change both the zero time and the flux baseline ( 1 value overstatement).There is no additional information here (about the target values of the surface parameters  and ), but the triple penetration experiment allows refining .For the model numerical experiment we have  01 ≈  02 ≈  03 ≈ 20.82 s.Let us analyze the steady state flux balance equation:

Isotherm: Initial
Asymptotic analysis shows that the dependence (√) has a parabolic shape ( ∝ ) at low inlet pressures : and a straight line form at relatively high pressures : Using the straight-line segment of the isotherm we find Γ/ℓ (the slope of the straight line), and knowing the estimate of  we determine the initial approximation of the solubility coefficient Γ.From the intersection of the straight line with the ordinate axis we find .Knowing the values of Γ = √2/ and  we compute  and Φ = Γ.
A graphic illustration (using a minimal required set of the calculated model  1,2,3 values) is presented in Figure 3.
Note that the obtained initial approximations are in good agreement with the original "forgotten" parameters as shown in Table 1.The penetration method is characterized by a significant measurement error, and data on the penetrating flux are required (and this, in turn, requires a more accurate determination of the vacuum system characteristics).The model of the dissolved hydrogen concentration jump at the inlet side is not very precise either.We are brought to a conclusion that the "communicating vessels" stage, where hydrogen pressures are measured over a long time, offers a much higher accuracy of measurements.Thus, the first stage of the aggregated experiment is perceived as a preliminary estimation of the , ,  coefficients.It is essential that the solution of the inverse problem of parametric identification is unique, since the results obtained for thin laboratory membranes are extrapolated (recalculated) to the dimensions of real-life structures.The results are "finetuned" by means of local variation of the preliminary values of , ,  in the fast hydrogen permeability model (the ODE system).

Thermal Desorption Spectroscopy (TDS)
Pressure values inside the chambers tend to equalize with time.The uniform equilibrium saturation  =  = Γ, Γ ≡ √2/,  =  = const corresponds to this pressure  = .We turn off the heating.Next, we create the vacuum at both sides of the membrane and begin slowly reheating the sample from room temperature (the degassing is practically completed) to activate the desorption process.
The thermal desorption flux of hydrogen from the sample is measured by mass spectrometer.What additional information can be obtained from the TDS experiment?The dependence () is assumed to be known (for example, as a result of a series of experiments in the DLR mode at various temperatures).We obtain no new information about () under vacuum conditions since modern, quite powerful vacuum systems allow neglecting the resorption.It remains to more precisely define the two-parameter Arrhenius-like dependence ().This is the bulk desorption coefficient (the coefficient of effective recombination of atoms into H 2 molecules):  0,ℓ =  2 0,ℓ ,  =  vol .The heating () is usually linear () =  0 + V.The heating rate V is not too high (<K/s).When the temperature maximum is reached (if degassing is not yet completed), then heating is stopped: () =  max .We refine the boundary conditions taking into account that hydrogen atoms can accumulate at the surface during slow monotone heating and relatively low temperatures.

"
Surface-Bulk" Fluxes Balance.Model (29) of quick dissolution (local equilibrium) on the surface is derived from the more general ratios: Coefficients  − ,  + are the descriptors of the rate of dissolution in the bulk and transfer to the surface.When concentrations are nowhere near maximum and   ≈ 0 on the relative scale, we obtain (29), where  =  − / + .If the surface is isotropic (in terms of   − ≈   + ), then the parameter () is a little dependent on temperature.The density of the H atoms adsorption flux can be modeled by the term 2[1 − ] 2 for balance equations ( 30), (31).For the ranges of weak concentrations and sufficiently high working temperatures the degree of surface coverage satisfies () = ()/ max ≪ 1.These simplifications are in agreement with the limited information capacities of the TDS experiments.Experimental data are more easily approximated for a large number of parameters.But the uniqueness of the estimations is then failed, and thus essential errors may occur at extrapolation of the results "from thin plate to wall".
The generalized quick dissolution coefficient  of local "surface-bulk" equilibrium is assumed to be constant (  − ≈   + ) in the given heating range.The TDS method fulfills the symmetry conditions: The information capacities of the initial and final stages of the TDS experiment are low.So it is sufficient to limit ourselves to  ≤  * , when the flux from the sample has decreased by an order of magnitude compared to the maximum.The experimental data are the desorption density curves () or TDS spectra () (() ↔ ) for different saturation conditions and heating rates.The conversion (⋅)  → (⋅) is made more specific by taking into account the parameters of the experimental unit.Modern equipment provides a means of creating deep vacuum (10 −9 -10 −7 Torr).For this reason, the component  ≡ 2 ( ≫ 1) is the key control for the saturation stage, but resorption is neglected for the degassing stage ( ≪ 1).

Diffusion Model with Reversible
Capture.We can take into account different channels of diffusion, but the information content of the TDS experiment is limited.Therefore, the coefficient  is assumed to be an integral effective index.Seeking a write-up of the TDS peaks set, it is handier to use the following model: where  ] (, ) is the concentration of hydrogen atoms captured by defects of different types;  ∓ ] are the coefficients of H capture and release by traps;  ] ≡  ] (, )/ ] max is the defects saturation degree ( ] max = max  ] ).Capture is taken into account at its integral level for practical purposes.A more precise definition of the defects' geometry and distribution would add complexity to the model.If, for instance, the defect is not a microcavity but hydride phase inclusions, then at the degassing stage the corresponding coefficient  −  () is identically zero and  +  () value is positive only if the critical temperature is reached: () ≥  crit .It is easy to simulate the required number of TDS-peaks using different binding energies (coefficients   ).Numerical algorithms based on difference schemes and modeling results were presented in [24,26].In this paper we use only (1) ( =  eff ).For thin membranes of quickly permeable material used in experiments this approximation is usually sufficient.

Equilibrium Analysis.
Let us take a brief look at the equilibrium ratios at the accepted detail level of modeling.We assume that pressure and temperature are constant.Formally, equilibrium is characterized by all derivatives being equal to zero.Keeping in mind the extensively used Sieverts' law, we shall observe proportionality to the √ value.
(1) On the surface the following ratio is applied: Here,  1 is the degree of surface coverage (next   have a similar meaning).Let us formulate  1 : (2) In the "surface-bulk" equilibrium we have ( 2 ≡ / max ,  1 ≡  −  max ,  2 ≡  +  max ,  ≡  1 / 2 ), then from Formula (33) we obtain (3) For definiteness, we take into account only one type of traps ( = 1).We use (36).For symmetry, we add the saturation factor [1− 2 ] for  =  ] .Here we have  3 ≡ / max ,  ≡  +  max /[ −  max ]; then we obtain The "saturation-degassing" experiment provides information about a general average concentration c in the bulk  = ℓ (end surfaces are neglected): Normalize c by  max and consider the dependence Let us focus on the curves in the axes (√, Θ).Numerical results are presented in Figure 5 ( =  + / − ,  max = 10 18 , Advances in Mathematical Physics ℓ = 0.1 cm).Quite many parameters have an effect, but one can usually estimate the orders of magnitude to find "the distance to linearity" (to the Sieverts' law adequacy range).
The curves Θ = (√) have the "growing wave" form.It is noticeable in a wide pressure range only.For the given parameters, the inflection point is the most vivid on the curves marked with pentagrams.The analysis of each additive component in (42) for the total concentration shows that the position of the inflection point is determined by the moment when the first and third addends have turned to the saturation mode with the following prevalent rise of  2 .Note that the curves for a narrower pressure range are practically linear, in line with ranges of the Sieverts' law.In a wide range of pressures, the wave-like character of the graphs is observed experimentally [37].

Richardson Approximation.
When experimentally estimating the hydrogen permeability of structural materials, Richardson approximation is often used for the penetration flux density: Membranes are usually very thin and their permeability is relatively high.We can accept for the concentration gradient In accordance with the Fick's law, the formula  = −[ out −  in ]/ℓ is "exact".Here, the function () is the monoatomic hydrogen permeation flux density.The actual boundary volume concentrations cannot however be measured, and thus a quasi-equilibrium approximation is used, where the concentration  is substituted with the equilibrium concentration  = Γ√ in accordance with Sieverts' law (Γ is the solubility coefficient, or solubility for short).On account of the inequalities  in >  in ,  out <  out , this substitution overrates the second factor in the formula Hence, the value of  needs to be formally marked down to retain the equality.So, if the value of Φ (permeability) is found from experimental data by fitting as indicated in formula (43), then we obtain Φ < Γ.The Φ and Γ values from the (typically used) formula Φ = Γ determine the lower bound of the diffusion coefficient .Furthermore, dissolved diffusing hydrogen is mainly involved in the steady state (quasi-steady state) permeability mode.In the context of the model ( 29)-( 31), we obtain  = Γ√, Γ = √2/.The "saturation-degassing" experiment yields the total value c >  and an overstated (for the permeability problem) value Γ max : c = Γ max √.Hence, we can estimate the Φ value by fitting using formula (43).This information has practical value as a convenient coefficient for translation from pressures to flux.If, however, the  value is taken from one experiment (or paper) and the Γ value is taken from another source, then, strictly speaking, we get a ranking of three different numbers Φ < Γ < Γ max .If the material has a high level of trapping by defects, then the calculated permeability can be an order of magnitude higher than the real permeability.The permeability coefficient Φ (as a parameter of (43)) has an -form (Arrhenius-like) of the saturation curve based on the order of pressures.It is only for relatively high pressures (where  0,ℓ are near Sieverts' concentrations) that we get Φ ≈ Γ.

Functional Differential TDS Equation
Identification of TDS spectra is required not only to reveal the causes of different thermal desorption peaks, but also to enable numerical extrapolation and generalization of the results received for laboratory samples (ℓ usually is fractions of mm).Model (36) gives the possibility to get any number of peaks using traps with different parameters  ± ] .The question, however, is whether different peaks can occur when degassing an almost homogeneous material.To answer this question let us restrict ourselves to the basic diffusion equation ( 1), but retaining symmetric dynamical boundary conditions ( 29)-( 31) and (35).The surface is considered isotropic in terms of  = const over the heating range.The resorption during vacuum building is neglected.Thus, we are limited to a minimal number of parameters for the model which takes into account the dynamical interplay of surface processes and diffusion.In the following, let us operate at this approximation.
The comparison of simulated and experimental TDS spectra with a focus on parametric identification requires only the surface concentration ( =  2 ).It is reasonable to try to avoid iterative solution of the boundary value problem for interim approximations of the model parameters  0 ,   ,  0 ,   ,  0 ,   , .To this end, we will run the transformations to reduce the problem to the integration of a low order ODE system.

Derivation of Riccati
Let us replace the time   = ∫  0   (keeping the former notation ): Here () is considered as the parameter and ( 46) is an additional relation for the linear problem (45).We perform a replacement to get homogenous boundary conditions: Let us write the solution of the linear boundary value problem using the source function (Green's function): Boundary conditions contain the derivative ĉ (, 0): For  =  we have divergent series.So, term-by-term integration is implied.For the original time  we get Finally, the dynamic boundary condition is written in the integrodifferential form: The resultant equation (51) with quadratic nonlinearity will be classified as a functional differential Riccati equation of the neutral type.The equation is equivalent to the original boundary value problem in that the solution () uniquely determines the solution (, ).The analogy with the functional differential equation of the form ẋ () = F[, (), ẋ (),  [0,] , ẋ [0,] ] of the neutral type [38] is that it is impossible to eliminate the derivative q on the right side of the equation through integration by parts lest a divergent series arises.We are concerned with the time interval [ 1 ,  2 ] ⊂ (0,  * ), which corresponds to the TDS peak.Measurement for  ≈ 0,  * yields little information.There is a voluminous body of literature on Riccati-type equations (including matrix equations for the optimal control theory).
6.3.Initial Data.We specify the factor b() for quadratic nonlinearity.Initial saturation is conducted under relatively high temperature  =  = const and pressure  = const.
After the steady state of saturation is attained we get Hereafter, let us be guided by a maximum limit of surface concentration at around  max ∼ 10 15 -10 16 (monolayer in the context of geometric statics).If during initial saturation the surface concentration for the given model parameters is  > 10 14 , the degree of surface coverage must be taken into account.Then, to calculate the initial value of  we take the ratio 2[1 −  −1 max ] 2 =  2 instead of 2 =  2 .In the meantime, this a priori restriction (arising from the assumption that stationary "balls" are ordered geometrically on a plane) is highly questionable.For a dynamics model, it appears that the concentration threshold  max may be higher, if it is specified what meaning is attached to the term "bubbling" surface layer (at a quite high temperature).The real  value is somewhat conventional, since it is strongly influenced by the experiment pretreatment (building of vacuum before the start of heating).However, most of the hydrogen is in the bulk, the diffusion equation (and the diffusion process itself) has a smoothing effect, and we are interested in evident thermal desorption peaks, since the initial and final time intervals of the experiment offer little information.For initial calculations it is therefore sufficient to correctly estimate the magnitude of the "effective"  concentration.High precision requirements are noncritical here.(58) can be approximated step-by-step by a low order ODE system on segments of dimensionless time of length ℎ.
The greatest contribution to the integral is made by the value V (),  ≈  due to nonlimited (but integrable) singularity.Thus, the quadratic approximation is used due to function V concavity.Let us consider the current segment of time  ∈ [ℎ, ( + 1)ℎ],  ≥ 0. The conditions determine the values of (), () (constants on ): ∈ [ℎ, ] ( ∈ [ℎ, ( + 1)ℎ]).Represent the integral from Formula (58) as a sum  ∈ [0, ℎ],  ∈ [ℎ, ].For the second additive, the singularity V ()/√ −  is explicitly integrated by substituting (9).The trapezoid rule and mean value theorem are used for the integral without singularity ((+0)/ √ +0 = 0) To approximate the first integral we use only several terms (for definiteness  = 1, 3) on the right-hand side taking into account the factor  2 .The negative additives ( V < 0) thereby drop out.Let us compensate for that by replacing ℎ with  and introducing the additional variables  1,2 (): As a result we obtain an ODE system instead of (58):  Qualitatively, thermal desorption spectra of metal structural materials have a typical form.Figure 6 illustrates numerically simulated TDS spectra for the above-listed hydrogen permeability parameter values of the structural materials at different heating rates.Even when defects (traps) are not taken into account, the suggested model can yield different types of two-peak graphs (where the low temperature peak is more pronounced or where the peaks are comparable).For reference, the change of the so-called transport parameter  = ℓ vol / = ℓ/( 2 ) is shown in the top part of the graphs.This parameter is crucial for the study of membrane hydrogen permeability [5].We supposed that SLR corresponds to  < 10 −2 and DLR corresponds to  > 10 4 .For beryllium and steel, the low temperature TDS peak takes place in the range where diffusion and surface processes play commensurable roles.The high temperature TDS peak occurs in DLR.For nickel, surface processes and diffusion have similar effect throughout the experimental temperature range, and only a peak-like step appears in the TDS spectrum at a low temperature.
Let us focus on TDS spectra for tungsten because these spectra are qualitatively different from the previous ones.The aforementioned parameter values "for tungsten" are, of course, formal.Besides microimpurities, the parameters depend on how the sample surface was treated.This is especially true because in practice a plate is thin, the volume is small, and the effect of surface processes is more vivid.This is one of the reasons for such a high variation of the quantitative estimates of hydrogen permeability parameters.Another reason (but not the last) is the following.Different models are used for experimental data treatment (although coefficients formally have the same name).For the model data accepted here, a narrow splash followed by a lengthy movement to the second peak (less visible) was observed.In this model the sample has no high-capacity traps (standard diffusion equation), but more detailed consideration is given to the surface (see ( 29)-( 31), the plate is thin).At first, near-surface hydrogen is actively desorbed.Then, diffusion is slowly activated by heating.The concentration gradient is substantial, and pumping to the surface is growing.These effects are usually categorically explained by the presence of traps, not by the dynamics of "surface-bulk" interplay.An amazingly similar experimental curve (marked with black circles) is found in [41], Figure 1, although this paper discussed deuterium implantation on tungsten.The essence of our model results is that various TDS peaks may have other causes apart from the routinely blamed traps.
Figure 7 separately shows the dynamics of surface processes and diffusion.High-temperature peaks correspond to the significantly enlarged desorption coefficient and activation of the diffusion afflux from the bulk under heating.Peaks at a low temperature take place where the diffusion towards the surface is not high but near-surface concentrations are higher.
Let us briefly present the numerical modeling algorithm.(ii) The low order ODE system of type (65) is numerically integrated step by step in dimensionless time (on an ℎ length interval).Other ODE approximations were presented in [27].This procedure is standard for every modern software.

Inverse Problem of Parametric Identification
The presented algorithm of numerical modeling allows quick scanning of different scenarios and operating conditions of a material (including the heating law and extrapolation of the results with ℓ increase).This statistical is useful when designing the strategy of experimental research.New materials (various alloys) first have to be analyzed for their hydrogen permeability.In dealing with this task we encounter inverse mathematical physics problems, which are well-known for their difficulty.Let us suppose that the diffusion coefficient () is known (Daines-Berrer method is usually used for the DLR of permeability).Let us formulate the problem: to estimate the surface parameters  0 ,   ,  using the () (desorption flux) information.These conditions correspond to the real conditions of an experiment where desorption dynamics closely correlates with diffusion in the bulk.The problem of estimating the desorption coefficient ( vol = / 2 ) for the situation where accumulation on the surface can be neglected ( q ≈ 0) was presented in [25].
Since the function () ≡ (()) is known, it is reasonable to move directly to dimensionless time   , which is oriented at the characteristic time of the diffusion transfer ℓ 2 /.The old notation  is retained not to complicate the formulas.Let us present the nonlinear additive term in the system (65) in the following form: (66) Here,  is a function of the parameter .It is obtained through transformations using notations from (52) and (65).Elementary but somewhat lengthy formulas are omitted.The ＦＨ{ J(T (t)) q -2 (T(t) , g)} ＦＨ{ J(T (t)) q -2 (T(t) , g)} ＦＨ{ J(T (t)) q -2 (T(t) , g)} ＦＨ{ J(T (t)) q -2 (T(t) , g)} The presented linearization method is illustrated in Figure 8. Model spectra were plotted to test the algorithm.A system of five equations (for the variables V(),  1-4 ()) was numerically solved to obtain more accurate estimates for two-peak spectra.In addition, a series of experiments was performed for spectra with a random error not higher than 20% (we used standard Scilab uniform random number generator).The identification algorithm based on ODE system integration demonstrated the noise resistance of experimental data treatment.
The initial value of  yields the mass balance The experiments have demonstrated that the accuracy of such estimation decreases for  ≫ 1. Curves based on the "correct" numerical  value (the resultant curve is close to the straight line segment) and curves for 20% deviation from  are shown.The error of the identification algorithm testing is less than several percentages.This error appears due to the numerical error of direct and inverse problem solving.There is high sensitivity (appearing as vertical "beak" singularity) to situations where the "true"  value is exceeded.Mathematical singularity appears because, formally, the function V() changes sign, taking also negative values.We did not take efforts to avoid this "nonphysical" transition because it is a vivid sign of a wrong  value.The same numerical experiments were done for noisy spectra.Activation energies are determined more accurately because energy parameters more actively influence desorption during heating.From the mathematical point of view, there is no need to strive for visual agreement of simulated and experimental curves (because the experimental error is tens of percentages).ODE approximation is quite adequate.

Figure 3 :
Figure 3: Extrapolation of the isotherm of .

Figure 8 :
Figure 8: Estimation of desorption and dissolution parameters.

Table 1
The Final Stage of Isothermic Identification.Table1reflects only the computational errors arising when solving the direct and inverse problems.The real experimental data are noisy.The identification algorithm uses only integral operators thus ensuring the noise resistance of experimental data treatment.