# The Inverse Problem of Identification of Hydrogen Permeability Model

Authors: Yury V. Zaika; Natalia I. Rodchenkova; Ekaterina K. Kostikova
Journal: Advances in Mathematical Physics
 (2018)

---

## Abstract

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)).


---

## Body

## 1. 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–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–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–23]. This work is a continuation of [24–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 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–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-factor-at-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.


## 2. 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 (l is the plate thickness and S is its area). The temperature T 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:(1)∂tct,x=DT∂x2ct,x,t>0,c0,x=0,x∈0,l,where c(t,x) is the concentration of diffusing hydrogen. The diffusion coefficient D depends on the sample temperature T in an Arrhenius way with preexponential factor D0 and activation energy ED:D=D0exp⁡{-ED/[RT(t)]}.Initial data are determined by the fact that the sample had been preliminarily degassed:c(0,x)=0, z(0,x)=0, x∈[0,l].Nonlinear boundary conditions are derived from the material flux balance:(2)-dQindt=2sTμTp0t-bTc02tS=-SD∂xct,0,-dQoutdt=2sTμTplt-bTcl2tS=SD∂xct,l.Here, Qin(t), Qout(t) are the amounts of hydrogen atoms in the input chamber of volume Vin and output chamber of volume Vout, c0(t)≡c(t,0), cl(t)≡c(t,l). 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 Jp is related to the pressure p by the Hertz Knudsen formula: Jp=p/2πmkT (k is the Boltzmann constant, m is the mass of a hydrogen molecule). In the context of the experiment it is convenient to choose the following measurement units [l]=cm, [p]=Torr. Then we numerically obtain the dependence Jp=μp, μ(T)≈2.474·1022/T ([μ]=1H2/(Torrcm2s), [T]=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 2s. One can write s (as a parameter of the model) instead of 2s, but it is more convenient to interpret the dimensionless probability factor s as the fraction of absorbed H atoms within the 2s notation. Thus, 2sμp 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,J0,l=bc0,l2 are the densities of the desorption flux from the sample (deviation from the square desorption law is significant only at extreme temperatures) and b is the desorption coefficient. We also assume that s and b depend on the temperature in an Arrhenius way. Formally, the activation energy Es 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 ps=const is maintained at a constant temperature T=const on both sides of the membrane, the equilibrium concentration c¯ of the dissolved atomic diffusionally mobile hydrogen is finally established. By equating the derivatives in model (2) to zero, we get Sieverts’ law c¯∝ps:c¯=Γps, Γ≡2sμ/b.Let us clarify the experimental conditions. The volumesVin,out comprise several liters, the thickness of the membrane is l<mm, the area S is about 1cm2, and the inflow pressure p0(0) is within several hundreds Torr.It now remains to find the magnitudes ofQin, Qout. Within the time of transfer through the membrane the gas is in the thermodynamical quasi-equilibrium, wherefore we use the formula N=pV/(kT). Here, N is the number of gas particles occupying the volume V at the temperature T and the pressure p (in the SI system [p]=Pa, [V]=m3, [k]=J/K). Taking into account the relationships Torr=133.322Pa and Pa=J/m3 (formally), we get the following expressions for the corresponding pressures and volumes in the boundary conditions (2): Q=2N=αpV/T, α≈1.931·1019. Here, p, V, T are the numerical values in the selected units (Torr, cm3, 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μ(T). The gas inside the volumes Vin,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 Q (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 T 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 V.


### 2.2. 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:(3)Jt=-D∂xc=Dl-1c0t-clt≈JRt=DΓl-1p0t-plt.Let us formulate the problem of modelling the concentrations c0,l using the pressures p0,l (the problem is also of interest per se) without the quasi-equilibrium simplification c(t)=Γp(t). The quasi-stationary state is achieved within a time t0, which is short compared to the total experiment time (∂xc=-[c0(t)-cl(t)]/l). So, the original model (1)-(2) can be simplified (taking into account the formula Q=αpV/T):(4)p˙0,lt=∓β0,lc0t-clt,β0,l≡SDαVin,outl-1T,(5)2sμp0,lt-bc0,l2t=±Dl-1c0t-clt,t≥t0>0.Since by virtue of the “inlet-outlet” balance the equalities hold true, (6)p˙lt=-VinVout-1p˙0t⟹plt=plt0+VinVout-1p0t0-p0t,it is sufficient to express the concentrations c0,l(t)=c0,l(p0(t)) 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 l). The dimensionless normalized variables are convenient for numerical simulation:(7)X0,lt=1+2lc0,ltbD-1,a0,lt=4l2Γ2p0,ltb2D-2-1.In addition, the system of equations (5) is compactly written in the symmetric form a0+2Xl=X02, al+2X0=Xl2. For the variable X≡Xl we obtain the incomplete quartic equation [X2-al]2=4[2X+a0], 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 p˙0=f(p0). Therefore, we shall aim to derive differential equations for X0,l, since information about the dynamics of the boundary concentrations c0,l is also important.Differentiate (5) with respect to time and substitute the pressure derivatives from (4). For the variables X0,l we get the system(8)X˙0t=-2sM0X0-Xl·Xl-VinVout-1X0Xl-1,M0≡μSTαVin,(9)X˙lt=2sMlX0-Xl·X0-Vin-1VoutX0Xl-1,Ml≡μSTαVout.The change of variables in (7) determines the concentrations c0,l(t) from which the model pressures p0,l(t), t≥t0 are calculated using (5).Let us formulate step-by-step the numerical algorithm of modelling the pressuresp0,l(t) for the current values of D, b, s coefficients (the authors used the Scilab freeware). We target at the “normal” experimental conditions [11, 12, 29–31], including the values of p, T, l, V, S.(1) We fixt=t0: 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 X≡Xl we choose the root of the biquadratic polynomial [X2-al(t0)]2-4[2X+a0(t0)]. From physical considerations it follows that cl(t0)>c¯l(t0) and thus X>1+1+al=1+2lΓpl(t0)bD-1.(2) The system of equationsa0+2Xl=X02, al+2X0=Xl2 (t=t0) yields the missing value of X0(t0). Formally, one equation is enough, but we take into account averaging procedures including determination of the values of p0,l(t0).(3) We numerically integrate the ODE system (8), (9) with the obtained initial data. The change of variables in (7) defines the concentrations c0,l(t), which are used to calculate the model pressures p0,l(t), t≥t0 from (5).Computational experiments show that the model curves almost coincide (att≥t0) 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Φ=DΓ. All the variable parameters of the original model that influence the permeability, namely, D, b, s, are important when running the above algorithm. Thus, the fast hydrogen permeability model does not lose in informativeness concerning the considered transfer parameters.


## 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 (l is the plate thickness and S is its area). The temperature T 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:(1)∂tct,x=DT∂x2ct,x,t>0,c0,x=0,x∈0,l,where c(t,x) is the concentration of diffusing hydrogen. The diffusion coefficient D depends on the sample temperature T in an Arrhenius way with preexponential factor D0 and activation energy ED:D=D0exp⁡{-ED/[RT(t)]}.Initial data are determined by the fact that the sample had been preliminarily degassed:c(0,x)=0, z(0,x)=0, x∈[0,l].Nonlinear boundary conditions are derived from the material flux balance:(2)-dQindt=2sTμTp0t-bTc02tS=-SD∂xct,0,-dQoutdt=2sTμTplt-bTcl2tS=SD∂xct,l.Here, Qin(t), Qout(t) are the amounts of hydrogen atoms in the input chamber of volume Vin and output chamber of volume Vout, c0(t)≡c(t,0), cl(t)≡c(t,l). 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 Jp is related to the pressure p by the Hertz Knudsen formula: Jp=p/2πmkT (k is the Boltzmann constant, m is the mass of a hydrogen molecule). In the context of the experiment it is convenient to choose the following measurement units [l]=cm, [p]=Torr. Then we numerically obtain the dependence Jp=μp, μ(T)≈2.474·1022/T ([μ]=1H2/(Torrcm2s), [T]=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 2s. One can write s (as a parameter of the model) instead of 2s, but it is more convenient to interpret the dimensionless probability factor s as the fraction of absorbed H atoms within the 2s notation. Thus, 2sμp 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,J0,l=bc0,l2 are the densities of the desorption flux from the sample (deviation from the square desorption law is significant only at extreme temperatures) and b is the desorption coefficient. We also assume that s and b depend on the temperature in an Arrhenius way. Formally, the activation energy Es 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 ps=const is maintained at a constant temperature T=const on both sides of the membrane, the equilibrium concentration c¯ of the dissolved atomic diffusionally mobile hydrogen is finally established. By equating the derivatives in model (2) to zero, we get Sieverts’ law c¯∝ps:c¯=Γps, Γ≡2sμ/b.Let us clarify the experimental conditions. The volumesVin,out comprise several liters, the thickness of the membrane is l<mm, the area S is about 1cm2, and the inflow pressure p0(0) is within several hundreds Torr.It now remains to find the magnitudes ofQin, Qout. Within the time of transfer through the membrane the gas is in the thermodynamical quasi-equilibrium, wherefore we use the formula N=pV/(kT). Here, N is the number of gas particles occupying the volume V at the temperature T and the pressure p (in the SI system [p]=Pa, [V]=m3, [k]=J/K). Taking into account the relationships Torr=133.322Pa and Pa=J/m3 (formally), we get the following expressions for the corresponding pressures and volumes in the boundary conditions (2): Q=2N=αpV/T, α≈1.931·1019. Here, p, V, T are the numerical values in the selected units (Torr, cm3, 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μ(T). The gas inside the volumes Vin,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 Q (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 T 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 V.


## 2.2. 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:(3)Jt=-D∂xc=Dl-1c0t-clt≈JRt=DΓl-1p0t-plt.Let us formulate the problem of modelling the concentrations c0,l using the pressures p0,l (the problem is also of interest per se) without the quasi-equilibrium simplification c(t)=Γp(t). The quasi-stationary state is achieved within a time t0, which is short compared to the total experiment time (∂xc=-[c0(t)-cl(t)]/l). So, the original model (1)-(2) can be simplified (taking into account the formula Q=αpV/T):(4)p˙0,lt=∓β0,lc0t-clt,β0,l≡SDαVin,outl-1T,(5)2sμp0,lt-bc0,l2t=±Dl-1c0t-clt,t≥t0>0.Since by virtue of the “inlet-outlet” balance the equalities hold true, (6)p˙lt=-VinVout-1p˙0t⟹plt=plt0+VinVout-1p0t0-p0t,it is sufficient to express the concentrations c0,l(t)=c0,l(p0(t)) 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 l). The dimensionless normalized variables are convenient for numerical simulation:(7)X0,lt=1+2lc0,ltbD-1,a0,lt=4l2Γ2p0,ltb2D-2-1.In addition, the system of equations (5) is compactly written in the symmetric form a0+2Xl=X02, al+2X0=Xl2. For the variable X≡Xl we obtain the incomplete quartic equation [X2-al]2=4[2X+a0], 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 p˙0=f(p0). Therefore, we shall aim to derive differential equations for X0,l, since information about the dynamics of the boundary concentrations c0,l is also important.Differentiate (5) with respect to time and substitute the pressure derivatives from (4). For the variables X0,l we get the system(8)X˙0t=-2sM0X0-Xl·Xl-VinVout-1X0Xl-1,M0≡μSTαVin,(9)X˙lt=2sMlX0-Xl·X0-Vin-1VoutX0Xl-1,Ml≡μSTαVout.The change of variables in (7) determines the concentrations c0,l(t) from which the model pressures p0,l(t), t≥t0 are calculated using (5).Let us formulate step-by-step the numerical algorithm of modelling the pressuresp0,l(t) for the current values of D, b, s coefficients (the authors used the Scilab freeware). We target at the “normal” experimental conditions [11, 12, 29–31], including the values of p, T, l, V, S.(1) We fixt=t0: 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 X≡Xl we choose the root of the biquadratic polynomial [X2-al(t0)]2-4[2X+a0(t0)]. From physical considerations it follows that cl(t0)>c¯l(t0) and thus X>1+1+al=1+2lΓpl(t0)bD-1.(2) The system of equationsa0+2Xl=X02, al+2X0=Xl2 (t=t0) yields the missing value of X0(t0). Formally, one equation is enough, but we take into account averaging procedures including determination of the values of p0,l(t0).(3) We numerically integrate the ODE system (8), (9) with the obtained initial data. The change of variables in (7) defines the concentrations c0,l(t), which are used to calculate the model pressures p0,l(t), t≥t0 from (5).Computational experiments show that the model curves almost coincide (att≥t0) 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Φ=DΓ. All the variable parameters of the original model that influence the permeability, namely, D, b, s, are important when running the above algorithm. Thus, the fast hydrogen permeability model does not lose in informativeness concerning the considered transfer parameters.


## 3. Modelling of Hydrogen Permeability

### 3.1. 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–13, 29–33]. We fix T=673K, l=0.05cm, μ(T)=2.474×1022/T1H2/(Torrcm2s), D=2×10-5cm2/s, Γ=2×10201H/(cm3Torr), Φ=DΓ=4×10151H/(cmsTorr). We set the value s=1.2×10-4 and calculate the corresponding desorption coefficient b=2μs/Γ2=5.7×10-24cm4/s. The input pressures p¯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 300K. 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~<c¯. Within the model we determine c¯i and c~i by the formulas (10)c¯i=2μsp¯ib-1,2μsp¯-bc~2=J¯=Dc~l-1⟹c~i=-a+a2+Γ2p¯i,a≡D2bl-1.For the given values of the parameters we have c~1=1.06×1021<c¯1=1.09×1021, c~2=1.38×1021<c¯2=1.41×1021, c~3=1.64×1021<c¯3=1.67×1021.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 (11)∂tc=D∂x2c,0<x<l,ct,0=c~1,ct,l=0,c0,x=0,x∈0,l.The penetrating flux is J1=-Dcxl=Dc~1l-1[1+2∑n=1∞‍(-1)nexp⁡{-n2π2Dt/l2}].From the computational point of view it is convenient to introduce the dimensionless timet′=Dt/l2 which is oriented at the characteristic diffusion time l2/D. As small t→0 a singularity appears if we directly use the partial sum of the expression(12)ft′=1+2∑n=1∞-1nexp⁡-n2π2t′.Let us provide another expression for f using the properties of the Jacobi theta function. More precisely, we are interested in the function (13)θ3t,x=1+2∑n=1∞exp⁡-n2π2tcos⁡2nπx,t>0.We have an alternative presentation for x=0 [34]:(14)θ3t,0=1+2∑n=1∞exp⁡-n2π2t=1πt∑-∞∞exp⁡-n2t.The series on the left is rapidly converging for large t. But the series on the right is rapidly converging for small t. If we define θ(t)=∑‍exp⁡{-πn2t} (n∈Z,t>0), then we get θ(1/t)=tθ(t), or t∑‍exp⁡{-πn2t}=∑‍exp⁡{-πn2t-1}(n∈Z), which is known as the functional equation for the theta function [35]. After some auxiliary transformations for 0<t′≪1 we get the appropriate representation(15)ft′=2πt′∑m=2n-1exp⁡-m24t′n∈N.The functionf(t′) has an S-shaped form of the saturation curve (see the inset in Figure 1). Stage I ends with c(t∗,x)=c~1(l-x)l-1.Figure 1
The establishment of fluxesJ(t).Stage II. t∗→t0=0 is the new t time zero: (16)ct=Dcxx,0<x<l,ct,0=c~2,ct,l=0,c0,x=c~1l-xl-1.Stage II ends with c(t∗,x)=c~2(l-x)l-1. The penetrating flux is (17)J2t=Dl-1c~1+c~2-c~11+2∑n=1∞-1nexp⁡-n2π2Dtl2.Stage III. The formulas are similar to the cyclic interchange c~1→c~2, c~2→c~3.The result of connecting the stages into a single “experimental” curve of the penetrating fluxJ(t) (conventionally, t=(t1,t1∗+t2,t2∗+t3), J=(J1,J2,J3)) is shown in Figure 1.


### 3.2. Converting Pressure into Flux

During the real experiment, the gas pressure inside the outlet volumeV=Vout 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 ν([ν]=m3/s). We take as the framework the ideal gas state equation: pV=NkT. Here, [p]=Pa, [V]=m3, N is the number of particles (H2 molecules), and k is the Boltzmann constant. Differentiating on t we get p˙V=N˙kT. Let us calculate the particle balance over the time Δt: (18)ΔN=0.5JlSΔt-νΔtNV-1,Jl=-SD∂xct,l,where [N/V]=1H2/m3 (the concentration), [Jl]=1H/(m2s). The factor 0.5 is due to the fact that the diffusion flux is atomic, and the particle in the volume V is the H2 molecule. Let us divide the equation of the material balance by Δt and direct Δt→0. Finally, we get the differential equation: (19)N˙t=0.5JltS-νV-1Nt⟹p˙V=N˙kT=kT0.5JlS-νV-1N=kTSJlt2-νVkTNt︸ptV⟹p˙t=kTS2VJlt-νVpt,p0=0.Denoting θ1=kTS/(2V), θ0=V/ν, we get p˙=θ1Jl-p/θ0 (p0=0), wherefrom (20)pt=θ1∫0texp⁡τ-tθ0Jlτdτ,Jlt=2VkTSp˙t+2νkTSpt.The measured function p(t) is noised, so we first apply a smoothing procedure and only then calculate the derivative p˙(t). 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 Jl≈2νp(t)/(kTS). Asymptotically, a steady state value is established: p˙→0⇒J¯=2νp¯/(kTS).Only the operation of integration is needed to calculate the lag time (see below):(21)∫0tJlτdτ=2VkTSpt+2νkTS∫0tpτdτ,so a preliminary approximation of the derivative p˙(t) is unnecessary. It suffices to use the composite Simpson formula. One should remember that we use the following measurement units within this subsection [kT]=J=N×m, [L]=m, [p]=Pa=N/m2 (Pa≈7.5×10-3Torr), [J]=1H/(m2s). In the following, we return to the measurement units accepted in this article.


### 3.3. 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 (p0,l(t) are measured). Considering the new time zero we also have t≥0. We are so far talking about the direct problem of modelling hydrogen pressures inside the volumes Vin,out. We specify the values: S=0.5cm2, Vin=1500cm3, Vout=2200cm3, p0(0)=p¯0=p¯3, c¯=Γp¯0. We target here the experimental conditions and the data for the V85Ni15 alloy [11, 33].We numerically solve the initial boundary problem of hydrogen permeability:(22)∂tc=DT∂x2c,c0,x=l-xc~3l-1,2sTμTp0,lt-bTc0,l2t=∓D∂xcx=0,l,dQin,outdt=-2sμp0,lt-bc0,l2tS,Qin,out=αp0,ltVin,outT-1.The membrane temperature is taken in the dependences ofD, b, s, μ on T, and the temperature of the gas inside the chambers is taken in the expressions for Qin,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 Vin,out into the expressions for M0,l). To this end one should skip the initial time t0 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”.Figure 2
Dynamics of pressuresp0,l(t).


## 3.1. 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–13, 29–33]. We fix T=673K, l=0.05cm, μ(T)=2.474×1022/T1H2/(Torrcm2s), D=2×10-5cm2/s, Γ=2×10201H/(cm3Torr), Φ=DΓ=4×10151H/(cmsTorr). We set the value s=1.2×10-4 and calculate the corresponding desorption coefficient b=2μs/Γ2=5.7×10-24cm4/s. The input pressures p¯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 300K. 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~<c¯. Within the model we determine c¯i and c~i by the formulas (10)c¯i=2μsp¯ib-1,2μsp¯-bc~2=J¯=Dc~l-1⟹c~i=-a+a2+Γ2p¯i,a≡D2bl-1.For the given values of the parameters we have c~1=1.06×1021<c¯1=1.09×1021, c~2=1.38×1021<c¯2=1.41×1021, c~3=1.64×1021<c¯3=1.67×1021.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 (11)∂tc=D∂x2c,0<x<l,ct,0=c~1,ct,l=0,c0,x=0,x∈0,l.The penetrating flux is J1=-Dcxl=Dc~1l-1[1+2∑n=1∞‍(-1)nexp⁡{-n2π2Dt/l2}].From the computational point of view it is convenient to introduce the dimensionless timet′=Dt/l2 which is oriented at the characteristic diffusion time l2/D. As small t→0 a singularity appears if we directly use the partial sum of the expression(12)ft′=1+2∑n=1∞-1nexp⁡-n2π2t′.Let us provide another expression for f using the properties of the Jacobi theta function. More precisely, we are interested in the function (13)θ3t,x=1+2∑n=1∞exp⁡-n2π2tcos⁡2nπx,t>0.We have an alternative presentation for x=0 [34]:(14)θ3t,0=1+2∑n=1∞exp⁡-n2π2t=1πt∑-∞∞exp⁡-n2t.The series on the left is rapidly converging for large t. But the series on the right is rapidly converging for small t. If we define θ(t)=∑‍exp⁡{-πn2t} (n∈Z,t>0), then we get θ(1/t)=tθ(t), or t∑‍exp⁡{-πn2t}=∑‍exp⁡{-πn2t-1}(n∈Z), which is known as the functional equation for the theta function [35]. After some auxiliary transformations for 0<t′≪1 we get the appropriate representation(15)ft′=2πt′∑m=2n-1exp⁡-m24t′n∈N.The functionf(t′) has an S-shaped form of the saturation curve (see the inset in Figure 1). Stage I ends with c(t∗,x)=c~1(l-x)l-1.Figure 1
The establishment of fluxesJ(t).Stage II. t∗→t0=0 is the new t time zero: (16)ct=Dcxx,0<x<l,ct,0=c~2,ct,l=0,c0,x=c~1l-xl-1.Stage II ends with c(t∗,x)=c~2(l-x)l-1. The penetrating flux is (17)J2t=Dl-1c~1+c~2-c~11+2∑n=1∞-1nexp⁡-n2π2Dtl2.Stage III. The formulas are similar to the cyclic interchange c~1→c~2, c~2→c~3.The result of connecting the stages into a single “experimental” curve of the penetrating fluxJ(t) (conventionally, t=(t1,t1∗+t2,t2∗+t3), J=(J1,J2,J3)) is shown in Figure 1.


## 3.2. Converting Pressure into Flux

During the real experiment, the gas pressure inside the outlet volumeV=Vout 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 ν([ν]=m3/s). We take as the framework the ideal gas state equation: pV=NkT. Here, [p]=Pa, [V]=m3, N is the number of particles (H2 molecules), and k is the Boltzmann constant. Differentiating on t we get p˙V=N˙kT. Let us calculate the particle balance over the time Δt: (18)ΔN=0.5JlSΔt-νΔtNV-1,Jl=-SD∂xct,l,where [N/V]=1H2/m3 (the concentration), [Jl]=1H/(m2s). The factor 0.5 is due to the fact that the diffusion flux is atomic, and the particle in the volume V is the H2 molecule. Let us divide the equation of the material balance by Δt and direct Δt→0. Finally, we get the differential equation: (19)N˙t=0.5JltS-νV-1Nt⟹p˙V=N˙kT=kT0.5JlS-νV-1N=kTSJlt2-νVkTNt︸ptV⟹p˙t=kTS2VJlt-νVpt,p0=0.Denoting θ1=kTS/(2V), θ0=V/ν, we get p˙=θ1Jl-p/θ0 (p0=0), wherefrom (20)pt=θ1∫0texp⁡τ-tθ0Jlτdτ,Jlt=2VkTSp˙t+2νkTSpt.The measured function p(t) is noised, so we first apply a smoothing procedure and only then calculate the derivative p˙(t). 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 Jl≈2νp(t)/(kTS). Asymptotically, a steady state value is established: p˙→0⇒J¯=2νp¯/(kTS).Only the operation of integration is needed to calculate the lag time (see below):(21)∫0tJlτdτ=2VkTSpt+2νkTS∫0tpτdτ,so a preliminary approximation of the derivative p˙(t) is unnecessary. It suffices to use the composite Simpson formula. One should remember that we use the following measurement units within this subsection [kT]=J=N×m, [L]=m, [p]=Pa=N/m2 (Pa≈7.5×10-3Torr), [J]=1H/(m2s). In the following, we return to the measurement units accepted in this article.


## 3.3. 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 (p0,l(t) are measured). Considering the new time zero we also have t≥0. We are so far talking about the direct problem of modelling hydrogen pressures inside the volumes Vin,out. We specify the values: S=0.5cm2, Vin=1500cm3, Vout=2200cm3, p0(0)=p¯0=p¯3, c¯=Γp¯0. We target here the experimental conditions and the data for the V85Ni15 alloy [11, 33].We numerically solve the initial boundary problem of hydrogen permeability:(22)∂tc=DT∂x2c,c0,x=l-xc~3l-1,2sTμTp0,lt-bTc0,l2t=∓D∂xcx=0,l,dQin,outdt=-2sμp0,lt-bc0,l2tS,Qin,out=αp0,ltVin,outT-1.The membrane temperature is taken in the dependences ofD, b, s, μ on T, and the temperature of the gas inside the chambers is taken in the expressions for Qin,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 Vin,out into the expressions for M0,l). To this end one should skip the initial time t0 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”.Figure 2
Dynamics of pressuresp0,l(t).


## 4. General Identification Technique

### 4.1. Determination of the Lag Time

For completeness, we briefly describe the method of estimating the diffusion coefficient proposed by Daines-Berrer. The curve of the fluxJ1(t) asymptotically moves to the stationary value J¯1=Dc~1/l. Hence, Q(t)=∫0t‍J(τ)dτ≈J¯1[t-l2/(6D)] (t≥t∗=t1∗). The intersection of the asymptote with the t axis gives the so-called lag time τ0=l2/(6D), which allows estimating the diffusion coefficient. Analytically, (23)τ0≈t∗-∫0t∗J1τJ¯1-1dτ,J1≈J¯1,t≥t∗.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 (J¯1=sup⁡J1(t)). In addition, the value τ0 does not depend on c~1. It is usually assumed that the locally equilibrium concentration c¯1=Γp¯1 is quickly established at the inlet, so one can additionally estimate the solubility Γ=2sμ/b and the permeability Φ=DΓ using the corresponding value J¯1=Dc¯1/l. 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 c0(t)≈c~1<c¯1 (t≥ε,ε≪1) is quickly established given that cl(t)≈0. The value of c~1 as such is yet to be clarified. Thus, only the estimate of the diffusion coefficient D is considered reliable at this stage.With a new zero time reference, integrating the expressionJ2(t) we get (24)∫0tJ2τ-J¯1dτ≈J¯2-J¯1·t-l26D-1t∗⟶t0=0,where J¯i=Dc~i/l, t≥t∗=t2∗. Formally, we obtain the same expressions for the lag time and the estimate of D if we change both the zero time and the flux baseline (J¯1 value overstatement). There is no additional information here (about the target values of the surface parameters b and s), but the triple penetration experiment allows refining D. For the model numerical experiment we have τ01≈τ02≈τ03≈20.82s.


### 4.2. Isotherm: Initial Estimates ofb, s

In experimental practices it is common to draw and analyze the isotherm, i.e., the curve of the steady state penetrating fluxJ¯ dependence on the inlet pressure p¯, while vacuum pumping is performed at the outlet side. If one targets the Sieverts’ law and the (quasi)equilibrium concentration c¯=Γp¯ at the inlet side (hence J¯=Dc¯/l), then it is natural to plot the dependence J¯=J¯(p¯).Let us analyze the steady state flux balance equation:a≡D(2lb)-1, (25)2sμp¯-bc~2=J¯=Dc~l-1⟹c~=a-1+1+a-2Γ2p¯,J¯=Dc~l-1=Dl-1a-1+1+a-2Γ2p¯.Asymptotic analysis shows that the dependence J¯(p¯) has a parabolic shape (J¯∝p¯) at low inlet pressures p¯: (26)a-2Γ2p¯≡x2x∝p¯⟹c~=a-1+1+x2=a0.5x2+⋯,and a straight line form at relatively high pressures p¯: (27)a-2Γ2p¯≫1⟹J¯=-D22l2b-1+DΓl-1p¯.Using the straight-line segment of the isotherm we find DΓ/l (the slope of the straight line), and knowing the estimate of D we determine the initial approximation of the solubility coefficient Γ. From the intersection of the straight line with the ordinate axis we find b. Knowing the values of Γ=2sμ/b and b we compute s and Φ=DΓ.A graphic illustration (using a minimal required set of the calculated modelJ¯1,2,3 values) is presented in Figure 3.Figure 3
Extrapolation of the isotherm ofJ¯.Note that the obtained initial approximations are in good agreement with the original “forgotten” parameters as shown in Table1.Table 1 Original value Approximate value D 2.00 × 1 0 - 5 2.0012 × 1 0 - 5 b 5.72×10-24 5.8371 × 1 0 - 24 s 1.20 × 1 0 - 4 1.2221 × 1 0 - 4 Γ 2.00 × 1 0 20 1.9984 × 1 0 20 Φ 4.00 × 1 0 15 3.9993 × 1 0 15


### 4.3. The Final Stage of Isothermic Identification

Table1 reflects 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. 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 theD, b, s 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 “fine-tuned” by means of local variation of the preliminary values of D, b, s in the fast hydrogen permeability model (the ODE system).


## 4.1. Determination of the Lag Time

For completeness, we briefly describe the method of estimating the diffusion coefficient proposed by Daines-Berrer. The curve of the fluxJ1(t) asymptotically moves to the stationary value J¯1=Dc~1/l. Hence, Q(t)=∫0t‍J(τ)dτ≈J¯1[t-l2/(6D)] (t≥t∗=t1∗). The intersection of the asymptote with the t axis gives the so-called lag time τ0=l2/(6D), which allows estimating the diffusion coefficient. Analytically, (23)τ0≈t∗-∫0t∗J1τJ¯1-1dτ,J1≈J¯1,t≥t∗.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 (J¯1=sup⁡J1(t)). In addition, the value τ0 does not depend on c~1. It is usually assumed that the locally equilibrium concentration c¯1=Γp¯1 is quickly established at the inlet, so one can additionally estimate the solubility Γ=2sμ/b and the permeability Φ=DΓ using the corresponding value J¯1=Dc¯1/l. 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 c0(t)≈c~1<c¯1 (t≥ε,ε≪1) is quickly established given that cl(t)≈0. The value of c~1 as such is yet to be clarified. Thus, only the estimate of the diffusion coefficient D is considered reliable at this stage.With a new zero time reference, integrating the expressionJ2(t) we get (24)∫0tJ2τ-J¯1dτ≈J¯2-J¯1·t-l26D-1t∗⟶t0=0,where J¯i=Dc~i/l, t≥t∗=t2∗. Formally, we obtain the same expressions for the lag time and the estimate of D if we change both the zero time and the flux baseline (J¯1 value overstatement). There is no additional information here (about the target values of the surface parameters b and s), but the triple penetration experiment allows refining D. For the model numerical experiment we have τ01≈τ02≈τ03≈20.82s.


## 4.2. Isotherm: Initial Estimates ofb, s

In experimental practices it is common to draw and analyze the isotherm, i.e., the curve of the steady state penetrating fluxJ¯ dependence on the inlet pressure p¯, while vacuum pumping is performed at the outlet side. If one targets the Sieverts’ law and the (quasi)equilibrium concentration c¯=Γp¯ at the inlet side (hence J¯=Dc¯/l), then it is natural to plot the dependence J¯=J¯(p¯).Let us analyze the steady state flux balance equation:a≡D(2lb)-1, (25)2sμp¯-bc~2=J¯=Dc~l-1⟹c~=a-1+1+a-2Γ2p¯,J¯=Dc~l-1=Dl-1a-1+1+a-2Γ2p¯.Asymptotic analysis shows that the dependence J¯(p¯) has a parabolic shape (J¯∝p¯) at low inlet pressures p¯: (26)a-2Γ2p¯≡x2x∝p¯⟹c~=a-1+1+x2=a0.5x2+⋯,and a straight line form at relatively high pressures p¯: (27)a-2Γ2p¯≫1⟹J¯=-D22l2b-1+DΓl-1p¯.Using the straight-line segment of the isotherm we find DΓ/l (the slope of the straight line), and knowing the estimate of D we determine the initial approximation of the solubility coefficient Γ. From the intersection of the straight line with the ordinate axis we find b. Knowing the values of Γ=2sμ/b and b we compute s and Φ=DΓ.A graphic illustration (using a minimal required set of the calculated modelJ¯1,2,3 values) is presented in Figure 3.Figure 3
Extrapolation of the isotherm ofJ¯.Note that the obtained initial approximations are in good agreement with the original “forgotten” parameters as shown in Table1.Table 1 Original value Approximate value D 2.00 × 1 0 - 5 2.0012 × 1 0 - 5 b 5.72×10-24 5.8371 × 1 0 - 24 s 1.20 × 1 0 - 4 1.2221 × 1 0 - 4 Γ 2.00 × 1 0 20 1.9984 × 1 0 20 Φ 4.00 × 1 0 15 3.9993 × 1 0 15


## 4.3. The Final Stage of Isothermic Identification

Table1 reflects 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. 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 theD, b, s 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 “fine-tuned” by means of local variation of the preliminary values of D, b, s in the fast hydrogen permeability model (the ODE system).


## 5. Thermal Desorption Spectroscopy (TDS)

Pressure values inside the chambers tend to equalize with time. The uniform equilibrium saturationc=c¯=Γp¯, Γ≡2μs/b, T=T¯=const corresponds to this pressure p=p¯. 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 dependenceD(T) 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 s(T) 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 b(T). This is the bulk desorption coefficient (the coefficient of effective recombination of atoms into H2 molecules): J0,l=bc0,l2, b=bvol.The heatingT(t) is usually linear T(t)=T0+vt. 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: T(t)=Tmax. We refine the boundary conditions taking into account that hydrogen atoms can accumulate at the surface during slow monotone heating and relatively low temperatures.


### 5.1. Dynamical Boundary Conditions

The surface is a significant potential energy barrier (see [2, p. 177–206]). Hence, the boundary conditions are modeled as follows: c0(t)≡c(t,0), cl(t)≡c(t,l),(28)c0,x=c¯,x∈0,l,t∈0,t∗,(29)c0t=gTq0t,clt=gTqlt,(30)q˙0t=2μTsTp0t-bTq02t+D∂xct,0,(31)q˙lt=2μTsTplt-bTql2t-D∂xct,l.Here, q0,l(t) are the surface concentrations (x=0,l,q˙≡dq/dt); g(T) is the parameter of local equilibrium between the surface and the subsurface bulk (coefficient of quick dissolution); b(T)=b0exp⁡{-Eb/[RT]} is the surface desorption coefficient: (32)J0,lt=bTtq0,l2t=bTtg-2Ttc0,l2t⟹b=bsurface=g2bvolume.The well-posedness of these dynamical boundary conditions was proved in [36].


### 5.2. “Surface–Bulk” Fluxes Balance

Model (29) of quick dissolution (local equilibrium) on the surface is derived from the more general ratios:(33)k-T1-c0,ltcmax-1q0,lt-k+T1-q0,ltqmax-1c0,lt=∓Dcx0,l.Coefficients k-, k+ are the descriptors of the rate of dissolution in the bulk and transfer to the surface. When concentrations are nowhere near maximum and Dcx≈0 on the relative scale, we obtain (29), where g=k-/k+. If the surface is isotropic (in terms of Ek-≈Ek+), then the parameter g(T) is a little dependent on temperature. The density of the H atoms adsorption flux can be modeled by the term 2μsp[1-θ]2 for balance equations (30), (31). For the ranges of weak concentrations and sufficiently high working temperatures the degree of surface coverage satisfies θ(t)=q(t)/qmax≪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”.


### 5.3. Symmetry Conditions

We shall hereinafter use a contracted notation for simplicity(34)bt≡bTt,Dt≡DTt,st≡sTt.The generalized quick dissolution coefficient g of local “surface-bulk” equilibrium is assumed to be constant (Ek-≈Ek+) in the given heating range. The TDS method fulfills the symmetry conditions:(35)p=p0,l,q=q0,l,c0=cl,D∂xc0=-D∂xcl,c¯x=c¯=const.The information capacities of the initial and final stages of the TDS experiment are low. So it is sufficient to limit ourselves to t≤t∗, 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 J(t) or TDS spectra J(T) (T(t)↔t) for different saturation conditions and heating rates. The conversion p(·)↦J(·) 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-7Torr). For this reason, the component P≡2μsp (P≫1) is the key control for the saturation stage, but resorption is neglected for the degassing stage (P≪1).


### 5.4. 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 coefficientD 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:(36)∂tc=D∂x2c-∑ν=1maν-1-Zνct,x-aν+zνt,x,∂tzν=aν-T1-Zνct,x-aν+Tzνt,x,where zν(t,x) is the concentration of hydrogen atoms captured by defects of different types; aν∓ are the coefficients of H capture and release by traps; Zν≡zν(t,x)/zmaxν is the defects saturation degree (zmaxν=max⁡zν). 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 aj-(T) is identically zero and aj+(T) value is positive only if the critical temperature is reached: T(t)≥Tcrit. It is easy to simulate the required number of TDS-peaks using different binding energies (coefficients Ea). Numerical algorithms based on difference schemes and modeling results were presented in [24, 26]. In this paper we use only (1) (D=Deff). For thin membranes of quickly permeable material used in experiments this approximation is usually sufficient.


### 5.5. 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 thep value.(1) On the surface the following ratio is applied:(37)2μsp1-θ12-bq2=0,θ1≡qqmax,ν2≡2μsbqmax2.Here, θ1 is the degree of surface coverage (next θi have a similar meaning). Let us formulate θ1: (38)pν21-θ12=θ12⟹θ1=νp1+νp,νp≪1⟹θ1∝p.(2) In the “surface–bulk” equilibrium we have (θ2≡c/cmax, k1≡k-qmax, k2≡k+cmax, γ≡k1/k2), then from Formula (33) we obtain (39)k1θ11-θ2=k2θ21-θ1⟹θ2=γθ11+γ-1θ1.(3) For definiteness, we take into account only one type of traps (m=1). We use (36). For symmetry, we add the saturation factor [1-θ2] for z=zν. Here we have θ3≡z/zmax, a≡a+zmax/[a-cmax]; then we obtain (40)θ21-θ3=θ31-θ2a⟹θ3=θ2a+1-aθ2.The “saturation-degassing” experiment provides information about a general average concentrationc~ in the bulk V=Sl (end surfaces are neglected): (41)Vc~=2θ1qmaxS+θ2cmaxV+θ3zmaxV.Normalize c~ by cmax and consider the dependence(42)Θp≡c~cmax=2qmaxlcmaxθ1+θ2+θ3zmaxcmax.Let us focus on the curves in the axes (p,Θ). Numerical results are presented in Figure 5 (a¯=a+/a-, cmax=1018, l=0.1cm). 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Θ=f(p) 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].


### 5.6. Richardson Approximation

When experimentally estimating the hydrogen permeability of structural materials, Richardson approximation is often used for the penetration flux density:(43)Jt=Φpint-pouttl-1,pout<pin.Membranes are usually very thin and their permeability is relatively high. We can accept for the concentration gradient cx=[cout-cin]/l. In accordance with the Fick’s law, the formula J=-D[cout-cin]/l is “exact”. Here, the function J(t) 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 c is substituted with the equilibrium concentration c¯=Γp in accordance with Sieverts’ law (Γ is the solubility coefficient, or solubility for short). On account of the inequalities c¯in>cin, c¯out<cout, this substitution overrates the second factor in the formula J=D[cin-cout]/l:c¯in-c¯out>cin-cout. Hence, the value of D 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 Φ<DΓ. The Φ and Γ values from the (typically used) formula Φ=DΓ determine the lower bound of the diffusion coefficient D. 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 c¯=Γp, Γ=g2μs/b. The “saturation-degassing” experiment yields the total value c~>c¯ and an overstated (for the permeability problem) value Γmax:c~=Γmaxp.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 D 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 Φ<DΓ<DΓmax. If the material has a high level of trapping by defects, then thecalculated permeability can be an order of magnitude higher than the real permeability. The permeability coefficient Φ (as a parameter of (43)) has an S-form (Arrhenius-like) of the saturation curve based on the order of pressures. It is only for relatively high pressures (where c0,l are near Sieverts’ concentrations) that we get Φ≈DΓ.


## 5.1. Dynamical Boundary Conditions

The surface is a significant potential energy barrier (see [2, p. 177–206]). Hence, the boundary conditions are modeled as follows: c0(t)≡c(t,0), cl(t)≡c(t,l),(28)c0,x=c¯,x∈0,l,t∈0,t∗,(29)c0t=gTq0t,clt=gTqlt,(30)q˙0t=2μTsTp0t-bTq02t+D∂xct,0,(31)q˙lt=2μTsTplt-bTql2t-D∂xct,l.Here, q0,l(t) are the surface concentrations (x=0,l,q˙≡dq/dt); g(T) is the parameter of local equilibrium between the surface and the subsurface bulk (coefficient of quick dissolution); b(T)=b0exp⁡{-Eb/[RT]} is the surface desorption coefficient: (32)J0,lt=bTtq0,l2t=bTtg-2Ttc0,l2t⟹b=bsurface=g2bvolume.The well-posedness of these dynamical boundary conditions was proved in [36].


## 5.2. “Surface–Bulk” Fluxes Balance

Model (29) of quick dissolution (local equilibrium) on the surface is derived from the more general ratios:(33)k-T1-c0,ltcmax-1q0,lt-k+T1-q0,ltqmax-1c0,lt=∓Dcx0,l.Coefficients k-, k+ are the descriptors of the rate of dissolution in the bulk and transfer to the surface. When concentrations are nowhere near maximum and Dcx≈0 on the relative scale, we obtain (29), where g=k-/k+. If the surface is isotropic (in terms of Ek-≈Ek+), then the parameter g(T) is a little dependent on temperature. The density of the H atoms adsorption flux can be modeled by the term 2μsp[1-θ]2 for balance equations (30), (31). For the ranges of weak concentrations and sufficiently high working temperatures the degree of surface coverage satisfies θ(t)=q(t)/qmax≪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”.


## 5.3. Symmetry Conditions

We shall hereinafter use a contracted notation for simplicity(34)bt≡bTt,Dt≡DTt,st≡sTt.The generalized quick dissolution coefficient g of local “surface-bulk” equilibrium is assumed to be constant (Ek-≈Ek+) in the given heating range. The TDS method fulfills the symmetry conditions:(35)p=p0,l,q=q0,l,c0=cl,D∂xc0=-D∂xcl,c¯x=c¯=const.The information capacities of the initial and final stages of the TDS experiment are low. So it is sufficient to limit ourselves to t≤t∗, 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 J(t) or TDS spectra J(T) (T(t)↔t) for different saturation conditions and heating rates. The conversion p(·)↦J(·) 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-7Torr). For this reason, the component P≡2μsp (P≫1) is the key control for the saturation stage, but resorption is neglected for the degassing stage (P≪1).


## 5.4. 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 coefficientD 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:(36)∂tc=D∂x2c-∑ν=1maν-1-Zνct,x-aν+zνt,x,∂tzν=aν-T1-Zνct,x-aν+Tzνt,x,where zν(t,x) is the concentration of hydrogen atoms captured by defects of different types; aν∓ are the coefficients of H capture and release by traps; Zν≡zν(t,x)/zmaxν is the defects saturation degree (zmaxν=max⁡zν). 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 aj-(T) is identically zero and aj+(T) value is positive only if the critical temperature is reached: T(t)≥Tcrit. It is easy to simulate the required number of TDS-peaks using different binding energies (coefficients Ea). Numerical algorithms based on difference schemes and modeling results were presented in [24, 26]. In this paper we use only (1) (D=Deff). For thin membranes of quickly permeable material used in experiments this approximation is usually sufficient.


## 5.5. 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 thep value.(1) On the surface the following ratio is applied:(37)2μsp1-θ12-bq2=0,θ1≡qqmax,ν2≡2μsbqmax2.Here, θ1 is the degree of surface coverage (next θi have a similar meaning). Let us formulate θ1: (38)pν21-θ12=θ12⟹θ1=νp1+νp,νp≪1⟹θ1∝p.(2) In the “surface–bulk” equilibrium we have (θ2≡c/cmax, k1≡k-qmax, k2≡k+cmax, γ≡k1/k2), then from Formula (33) we obtain (39)k1θ11-θ2=k2θ21-θ1⟹θ2=γθ11+γ-1θ1.(3) For definiteness, we take into account only one type of traps (m=1). We use (36). For symmetry, we add the saturation factor [1-θ2] for z=zν. Here we have θ3≡z/zmax, a≡a+zmax/[a-cmax]; then we obtain (40)θ21-θ3=θ31-θ2a⟹θ3=θ2a+1-aθ2.The “saturation-degassing” experiment provides information about a general average concentrationc~ in the bulk V=Sl (end surfaces are neglected): (41)Vc~=2θ1qmaxS+θ2cmaxV+θ3zmaxV.Normalize c~ by cmax and consider the dependence(42)Θp≡c~cmax=2qmaxlcmaxθ1+θ2+θ3zmaxcmax.Let us focus on the curves in the axes (p,Θ). Numerical results are presented in Figure 5 (a¯=a+/a-, cmax=1018, l=0.1cm). 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Θ=f(p) 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].


## 5.6. Richardson Approximation

When experimentally estimating the hydrogen permeability of structural materials, Richardson approximation is often used for the penetration flux density:(43)Jt=Φpint-pouttl-1,pout<pin.Membranes are usually very thin and their permeability is relatively high. We can accept for the concentration gradient cx=[cout-cin]/l. In accordance with the Fick’s law, the formula J=-D[cout-cin]/l is “exact”. Here, the function J(t) 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 c is substituted with the equilibrium concentration c¯=Γp in accordance with Sieverts’ law (Γ is the solubility coefficient, or solubility for short). On account of the inequalities c¯in>cin, c¯out<cout, this substitution overrates the second factor in the formula J=D[cin-cout]/l:c¯in-c¯out>cin-cout. Hence, the value of D 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 Φ<DΓ. The Φ and Γ values from the (typically used) formula Φ=DΓ determine the lower bound of the diffusion coefficient D. 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 c¯=Γp, Γ=g2μs/b. The “saturation-degassing” experiment yields the total value c~>c¯ and an overstated (for the permeability problem) value Γmax:c~=Γmaxp.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 D 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 Φ<DΓ<DΓmax. If the material has a high level of trapping by defects, then thecalculated permeability can be an order of magnitude higher than the real permeability. The permeability coefficient Φ (as a parameter of (43)) has an S-form (Arrhenius-like) of the saturation curve based on the order of pressures. It is only for relatively high pressures (where c0,l are near Sieverts’ concentrations) that we get Φ≈DΓ.


## 6. 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 (l usually is fractions of mm). Model (36) gives the possibility to get any number of peaks using traps with different parameters aν±. 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 g=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 (J=bq2). It is reasonable to try to avoid iterative solution of the boundary value problem for interim approximations of the model parameters D0, ED, b0, Eb, s0, Es, g. To this end, we will run the transformations to reduce the problem to the integration of a low order ODE system.


### 6.1. Derivation of Riccati-Type Equation

The accepted TDS degassing model is∂tc=D(t)∂x2c, c(0,x)=c¯, (44)c0,lt=gqt,q˙t=-btq2t+Dt∂xct,0,Jt≡btq2t.Let us replace the time t′=∫0t‍Ddτ (keeping the former notation t):(45)∂tct,x=∂x2ct,x,c0,x=c¯,c0,l=gqt,(46)∂xc0=-∂xcl=q˙t+JtD-1t.Here q(t) 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: (47)c^=ct,x-gqt,∂tc^t,x=∂x2c^t,x+ft,ft=-gq˙t,c^0,x=φ^x=0,c^0,l=0.Let us write the solution of the linear boundary value problem using the source function (Green’s function): (48)c^t,x=∫0t∫0lG1x,ξ,t-τfτdξdτ,G1x,ξ,t≡2l∑n=1∞exp⁡-n2π2l2tsin⁡nπxlsin⁡nπξl.Boundary conditions contain the derivative c^x(t,0): (49)∂xc^0=-4gl∫0tq˙τ∑′exp⁡n2π2l2τ-tdτ,∑′≡∑n=1,3,5,….For τ=t we have divergent series. So, term-by-term integration is implied. For the original time t we get ∂xc(t,0)=∂xc^(t,0)=∂xc(t,l)=∂xc^(t,l), (50)∂xct,0=-4gl∑′∫0tq˙τexp⁡-n2π2l2∫τtDsdsdτ.Finally, the dynamic boundary condition is written in the integrodifferential form:(51)q˙t=-btq2t-4gDl∑′∫0tq˙τexp⁡-n2π2l2∫τtDsdsdτ.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 q(t) uniquely determines the solution c(t,x). The analogy with the functional differential equation of the form x˙(t)=F[t,x(t),x˙(t),x[0,t],x˙[0,t]] 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 [t1,t2]⊂(0,t∗), which corresponds to the TDS peak. Measurement for t≈0,t∗ yields little information. There is a voluminous body of literature on Riccati-type equations (including matrix equations for the optimal control theory).


### 6.2. Dimensionless Form of the Problem

For more comfortable modeling we turn to dimensionless variables using substitution rules:t′=∫0t‍D(τ)dτ/l2, x′=x/l, v=q/q¯ (c¯=gq¯). Retaining the notation t, we obtain(52)v˙t=-b~tv2t-4gl∑′∫0tv˙τexp⁡-n2π2t-τdτ,v0=1,where b~(t)≡q¯b(t)l2/D(t) is a dimensionless parameter of quadratic “desorption”.


### 6.3. Initial Data

We specify the factorb~(t) for quadratic nonlinearity. Initial saturation is conducted under relatively high temperature T=T¯=const and pressure p¯=const. After the steady state of saturation is attained we get (53)2μT¯sT¯p¯=bT¯q¯2,c¯=gq¯⟹c¯=gb0-1/22μsp¯T¯exp⁡Eb2RT¯-1,b~t≡q¯btl2D-1=c¯btl2gDt-1≡b~0exp⁡-Eb~RTt-1,Eb~≡Eb-ED,b~0≡l2D0-1b02μsp¯T¯exp⁡Eb2RT¯-1.Hereafter, let us be guided by a maximum limit of surface concentration at aroundqmax~1015–1016 (monolayer in the context of geometric statics). If during initial saturation the surface concentration for the given model parameters is q¯>1014, the degree of surface coverage must be taken into account. Then, to calculate the initial value of q¯ we take the ratio 2μsp¯[1-q¯qmax-1]2=bq¯2 instead of 2μsp¯=bq¯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 thresholdqmax may be higher, if it is specified what meaning is attached to the term “bubbling” surface layer (at a quite high temperature). The real q¯ 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” q¯ concentration. High precision requirements are noncritical here.


## 6.1. Derivation of Riccati-Type Equation

The accepted TDS degassing model is∂tc=D(t)∂x2c, c(0,x)=c¯, (44)c0,lt=gqt,q˙t=-btq2t+Dt∂xct,0,Jt≡btq2t.Let us replace the time t′=∫0t‍Ddτ (keeping the former notation t):(45)∂tct,x=∂x2ct,x,c0,x=c¯,c0,l=gqt,(46)∂xc0=-∂xcl=q˙t+JtD-1t.Here q(t) 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: (47)c^=ct,x-gqt,∂tc^t,x=∂x2c^t,x+ft,ft=-gq˙t,c^0,x=φ^x=0,c^0,l=0.Let us write the solution of the linear boundary value problem using the source function (Green’s function): (48)c^t,x=∫0t∫0lG1x,ξ,t-τfτdξdτ,G1x,ξ,t≡2l∑n=1∞exp⁡-n2π2l2tsin⁡nπxlsin⁡nπξl.Boundary conditions contain the derivative c^x(t,0): (49)∂xc^0=-4gl∫0tq˙τ∑′exp⁡n2π2l2τ-tdτ,∑′≡∑n=1,3,5,….For τ=t we have divergent series. So, term-by-term integration is implied. For the original time t we get ∂xc(t,0)=∂xc^(t,0)=∂xc(t,l)=∂xc^(t,l), (50)∂xct,0=-4gl∑′∫0tq˙τexp⁡-n2π2l2∫τtDsdsdτ.Finally, the dynamic boundary condition is written in the integrodifferential form:(51)q˙t=-btq2t-4gDl∑′∫0tq˙τexp⁡-n2π2l2∫τtDsdsdτ.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 q(t) uniquely determines the solution c(t,x). The analogy with the functional differential equation of the form x˙(t)=F[t,x(t),x˙(t),x[0,t],x˙[0,t]] 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 [t1,t2]⊂(0,t∗), which corresponds to the TDS peak. Measurement for t≈0,t∗ yields little information. There is a voluminous body of literature on Riccati-type equations (including matrix equations for the optimal control theory).


## 6.2. Dimensionless Form of the Problem

For more comfortable modeling we turn to dimensionless variables using substitution rules:t′=∫0t‍D(τ)dτ/l2, x′=x/l, v=q/q¯ (c¯=gq¯). Retaining the notation t, we obtain(52)v˙t=-b~tv2t-4gl∑′∫0tv˙τexp⁡-n2π2t-τdτ,v0=1,where b~(t)≡q¯b(t)l2/D(t) is a dimensionless parameter of quadratic “desorption”.


## 6.3. Initial Data

We specify the factorb~(t) for quadratic nonlinearity. Initial saturation is conducted under relatively high temperature T=T¯=const and pressure p¯=const. After the steady state of saturation is attained we get (53)2μT¯sT¯p¯=bT¯q¯2,c¯=gq¯⟹c¯=gb0-1/22μsp¯T¯exp⁡Eb2RT¯-1,b~t≡q¯btl2D-1=c¯btl2gDt-1≡b~0exp⁡-Eb~RTt-1,Eb~≡Eb-ED,b~0≡l2D0-1b02μsp¯T¯exp⁡Eb2RT¯-1.Hereafter, let us be guided by a maximum limit of surface concentration at aroundqmax~1015–1016 (monolayer in the context of geometric statics). If during initial saturation the surface concentration for the given model parameters is q¯>1014, the degree of surface coverage must be taken into account. Then, to calculate the initial value of q¯ we take the ratio 2μsp¯[1-q¯qmax-1]2=bq¯2 instead of 2μsp¯=bq¯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 thresholdqmax may be higher, if it is specified what meaning is attached to the term “bubbling” surface layer (at a quite high temperature). The real q¯ 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” q¯ concentration. High precision requirements are noncritical here.


## 7. Extracting Integrable Singularity

The functional differential equation of thermal desorption (52) has a singularity hindering its numerical solution. The function, (54)Θs=4∑′exp⁡-n2π2s,∑′≡∑n=1,3,5,…,takes finite values for s>0. The series is rapidly converging for large s. Formally, substituting s=0 (the integration variable τ reaches the upper limit t) we obtain a divergent series. This can be “fixed” by term-by-term integration. The order of terms in the series is O(n-2) (n=2i-1,i≥1), wherefore convergence is slow. Let us formulate the problem. The equation has to be approximated by a low order ODE system to enable application of standard software packages (for example, Scilab).Let us run some transformations using the theory of Jacobi theta functions to explicitly extract the integrable singularity. We consider the presentation (14): (55)θ3t,0=1+2∑n=1∞exp⁡-n2π2t=1πt∑-∞∞exp⁡-n2t.The series on the left is rapidly converging for large t. But the series on the right is rapidly converging for small t. Let us run auxiliary transformations: (56)1πt∑-∞∞exp⁡-n2t=1+2∑k=1∞exp⁡-k2π24t+2∑′=2∑′+14πt∑-∞∞exp⁡-n24t.Proceeding from here, the last series is subtracted from the first series and the result is doubled. The following expression is obtained for s>0: (57)Θs≡4∑′exp⁡-n2π2s=1-Qπs,Qs≡-2∑n=1∞qn2,q≡-exp⁡{-1/(4s)}. The series Q is very rapidly converging for small s. As s→+0 we have Q→0 and the integrable singularity Θ≈1/πs. The graph for Q(s) has an S-shaped (“Arrhenius-like”) saturation curve form and S(1)≈0.9996. The function Q/πs increases monotonically to a maximum (≈0.828 for s≈0.334) and then decreases monotonically. To represent the series Q(s) it is reasonable to use only a small number (5–8) of series terms. As an alternative to (52), using the above-described presentation of Θ(s) we obtain the equation(58)v˙t=-b~tv2t-ϰ∫0t1-Qt-τπt-τv˙τdτt≥0,v0=1,where the fraction (weak singularity under the integral) rapidly decreases from infinity (τ=t) to near zero (τ=t-1). For t>1, the lower limit of the integral can be replaced with t-1. Thus, the neglectable background can be easily identified in the original physical time using the relation 1=t2′-t1′=∫t1t2‍D(T(τ))dτ/l2. The compact functional differential (v˙≡dv/dt) TDS equation (58) with initial data v(0)=1 replaces the original nonlinear boundary value problem (from Section 6.1) with dynamic boundary conditions in the sense that formally only the dynamics of the surface concentration (desorption) is required for TDS spectrum construction. Note that equation (58) contains a Caputo fractional derivative d1/2v(t)/dt1/2=∫0t‍v˙(τ)[π(t-τ)]-1/2dτ. The theory of integrodifferential equations (the “residual” integral containing the function Q can be transformed by parts) is a surging field of modern mathematics and its applications.


## 8. Numerical Method and Computer Simulation

To be specific in the paper, we use data for nickel and steel (12Cr18Ni10Ti) [6], tungsten [5], and beryllium [39, 40]. Estimates depend substantially on the experimental conditions and sample pretreatment. So the values are perceived as a model for numerical illustrations. The parameters common for all the materials are l=0.1cm, T0=300K, [E]=J/mol. The assumed values of model parameters are the following: (steel) b0=1.28×10-9 cm2/s, Eb=80×103, D0=3.09×10-4 cm2/s, ED=31×103, g=50cm-1, Eg=-5×103, s0=0.6, Es=110×103, c¯=6×1017 at.H/cm3, T¯=1170K, T˙=1K/s; (Ni) b0=1.53×10-14, Eb=43.2×103, D0=7.5×10-3, ED=40×103, g=100, Eg=0, s0=1.8×10-2, Es=61.4×103, p¯=37.4 Torr, T¯=770, T˙=0.5; (W) b0=6·10-4, Eb=69.559×103, D0=4.1·10-3, ED=37.629×103, g0=104, Eg=15×103, s0=9·10-3, Es=15, p¯=670, T¯=1300, T˙=5; (Be) b0=3.08·10-9, Eb=57.43×103, D0=3·10-3, ED=28×103, g0=200, Eg=1.5×103, s0=1.44·10-4, Es=1.82×103, p¯=760, T¯=1150, T˙=5.The main role in the degassing dynamics belongs to quadratic desorption. We therefore approximate the integral term in (58). The after-effect horizon here is h<1 (τ∈[t,t-h] for dimensionless time, Q(1)≈0.999). Let us fix h~0.3–0.4 due to the smooth function Q(s)/s graph (Figure 4) and the trapezoid rule for numerical integration. The replacement of t′=∫0t‍Ddτ/l2 is naturally targeted at the diffusion time scale l2/D so that the step h corresponds to a significant segment of the experiment physical time. Then, TDS equation (58) can be approximated step-by-step by a low order ODE system on segments of dimensionless time of length h.Figure 4
FunctionsQ(s) and Q(s)/πs.Figure 5
Saturation degree.
(a)
Effect of defects (b)
Effect of the surfaceThe greatest contribution to the integral is made by the valuev˙(τ), τ≈t due to nonlimited (but integrable) singularity. Thus, the quadratic approximation (59)v˙τ≈v˙t+At-τ+Bt-τ2is used due to function v˙ concavity. Let us consider the current segment of time t∈[kh,(k+1)h], k≥0. The conditions (60)v˙τkh=v˙kh,∫khtv˙τdτ=vt-vkhvt≈vkh+v˙tt-khdetermine the values of A(t), B(t) (constants on τ): (61)v˙τ≈v˙t+2v˙t-v˙kht-kht-τ-3v˙t-v˙kht-kh2t-τ2,τ∈[kh,t](t∈[kh,(k+1)h]). Represent the integral from Formula (58) as a sum τ∈[0,kh], τ∈[kh,t]. For the second additive, the singularity v˙(τ)/t-τ is explicitly integrated by substituting (9). The trapezoid rule and mean value theorem are used for the integral without singularity (Q(+0)/+0=0)(62)∫khtQt-τπt-τv˙τdτ=v˙ξ∫khtQt-τπt-τdτ≈v˙kh+v˙t2π·12·Qt-kht-kh.To approximate the first integral (63)∫0kh1-Qt-τπt-τv˙τdτ=4∑′∫0khv˙τe-n2π2t-τdτwe use only several terms (for definiteness n=1,3) on the right-hand side taking into account the factor n2. The negative additives (v˙<0) thereby drop out. Let us compensate for that by replacing kh with t and introducing the additional variables w1,2(t): (64)wit≡∫0texp⁡-2i-12π2t-τv˙τdτ,w˙it=-2i-12π2wit+v˙t,w1,20=0.As a result we obtain an ODE system instead of (58):(65)v˙t=-b^tv2t+ϰ^π215+Qt-kh4t-khv˙kh-4ϰ^w1t+w2t,v0=1,w˙1t=⋯-4ϰ^+π2w1t-4ϰ^w2t,w10=0,w˙2t=⋯-4ϰ^w1t-4ϰ^+9π2w2t,w20=0,b^≡b~1+ϰ/π32/15-Qt-kh/4t-kh,ϰ^≡ϰ1+ϰ/π32/15-Qt-kh/4t-kh.Ellipsis stands for the right-hand side of the first line of the first equation (-b^⋯v˙(kh)). The series for Q(s) converges very rapidly. It is sufficient to use 5–8 terms for software.The sought functionv(t)=q(t)/q¯ on the current segment of dimensionless time (in the system (65) t=t′, t′∈[kh,(k+1)h], k≥0) is computed by the Runge-Kutta 4th-order method for integration of ODE systems (the authors used Scilab software). If the TDS spectrum has two peaks, then w3,4 should be additionally used to improve the modeling accuracy. Returning to the physical time, we get the model spectrum J(T)=b(T)q2(t) (t′↔t↔T(t)).Qualitatively, thermal desorption spectra of metal structural materials have a typical form. Figure6 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 W=lbvolc¯/D=lbc¯/(Dg2) 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 W<10-2 and DLR corresponds to W>104. 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.Figure 6
TDS-spectra.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.Figure7 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.Figure 7
Dynamics of surface processes and diffusion.
(a)
Surface parameters (b)
Diffusion afflux to the surfaceLet us briefly present the numerical modelingalgorithm.(i)
Set the parameterss0, Es, b0, Eb, D0, ED, g. Determine the values of c¯, b~(t′), ϰ=gl under the saturation conditions p¯, T¯ as described inInitial saturation stage. If the equilibrium bulk concentration c¯ of dissolved diffusing hydrogen is known, then presetting of s(T¯) is not needed. Other scenarios of initial saturation are possible (see discussion inInitial saturation stage). The required adjustment of initial data is determined by the specifics of the actual TDS experiment.(ii)
The low order ODE system of type (65) is numerically integrated step by step in dimensionless time (on an h length interval). Other ODE approximations were presented in [27]. This procedure is standard for every modern software.(iii)
The graphJ(t)=b(T(t))q2(t) or spectrum J(T) is plotted in physical time using the function v(t′)=q(t′)/q¯.


## 9. 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 withl increase). This statistical information 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 D(T) is known (Daines-Berrer method is usually used for the DLR of permeability). Let us formulate the problem: to estimate the surface parameters b0, Eb, g using the J(t) (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 (bvol=b/g2) for the situation where accumulation on the surface can be neglected (q˙≈0) was presented in [25].Since the functionD(t)≡D(T(t)) is known, it is reasonable to move directly to dimensionless time t′, which is oriented at the characteristic time of the diffusion transfer l2/D. The old notation t is retained not to complicate the formulas. Let us present the nonlinear additive term in the system (65) in the following form: (66)b^tv2t=αϰ·lDtc¯·Jt,Jt≡btq2t.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 function J(t) (t=t′) is known from measurements. After substitution into (65) we obtain a system of three linear ODE. It is reasonable to add the variables w3,4 to improve the calculation accuracy. The measurements are usually noisy, but the function J(t) comprised of the right-hand side of the ODE and is smoothed by integration.Note that the right-hand sides of the equations now depend only on one estimated dissolution parameterϰ=gl. The identification algorithm can thus be decomposed. By setting the value of g we can eventually calculate the dependence q(t)=q(t;g) in the original physical time t. Be reminded now that(67)J=bq2⟹b=Jq-2⟹ln⁡b0-EbR·103·103T=ln⁡ψ,where ψ(t)≡J/q2, T=T(t). Since heating is monotonic t↔T(t), we can introduce the coordinates X≡103/T, Y≡ln⁡ψ. The parametric curve X(t), Y(t) is obtained on the plane (X,Y). Judging by the ratio (67), this curve has to be a straight line with a negative slope. Hence, we have the criterion for the choice of the “correct” value for g: the curve X(t), Y(t)=ln⁡{J(t)/q2(t;g)} has to be a straight line segment on the plane (X,Y). This subproblem is scalar, where only g varies.Formally, we prolong the straight line segment until it crosses the coordinate axes. The crossing with they-axis (X=0) yields the ln⁡b0 value. The crossing with the abscissa x¯=103×Rln⁡b0/Eb determines the Eb value.The algorithm is laborious due to iterative use of the procedure of numerical solving of the initial problem for the ODE system. It requires some effort and familiarity with mathematical packages, but it is much easier and quicker to use the standard built-in operation than to perform iterative solution of the original nonlinear boundary value problem with three varying parameters.The presented linearization method is illustrated in Figure8. Model spectra were plotted to test the algorithm. A system of five equations (for the variables v(t), w1–4(t)) 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.Figure 8
Estimation of desorption and dissolution parameters.The initial value ofg yields the mass balance(68)c¯l+2q¯=c¯l+2g-2=2∫0t∗Jτdτ.The experiments have demonstrated that the accuracy of such estimation decreases for g≫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(t) 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.


## 10. Conclusions

It is advisable to aggregate the thermal desorption (TDS) and penetration experiments (with and without vacuum pumping) to make the measurements substantially more informative for further estimation of the hydrogen permeability parameters and to improve the accuracy of parameter identification. This paper suggests a cascade experiment technique and the corresponding mathematical software.The identification algorithm uses only integral operators thus ensuring the noise resistance of experimental data treatment. The penetration method with vacuum pumping 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 molecular hydrogen pressures are measured over a long time, is characterized by a much higher accuracy of measurements.The first stage of the aggregated experiment is perceived as preliminary estimation of the diffusionD, desorption b=bvol, and absorption s 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 “fine-tuned” by means of local variation of the preliminary values of D, b, s in the ODE-model of fast hydrogen permeability.The final stage deals with the problem of identification of the spectra of hydrogen thermal desorption. This problem is of high relevance for the nuclear power industry. Qualitatively, the identification consists of revealing the causes of desorption peaks. It is usually assumed that TDS peaks appear due to hydrogen release from traps with different binding energies. Here, it was demonstrated using a diffusion model for a homogeneous material that if surface processes are taken into account, two-peak spectra can be obtained even for very thin experimental samples. The tendency to resort to the “theory of different traps” as the only explanation is understandable, but the volume of our samples was near zero for the trapping capacity to manifest itself.The nonlinear boundary value problem (standard diffusion equation with dynamical boundary conditions) is reduced to the functional differential equation for the surface concentration, because nothing but desorption dynamics is required to plot a TDS spectrum. An effective numerical algorithm oriented to the use of mathematical packages (including freeware) is proposed. The main final output of this part of the paper is a geometrically transparent method for solving the inverse problem of surface parameter identification where desorption and diffusion in the bulk are dynamically interrelated.



---

*Source: 4628346-2018-06-07.xml*