Advances in Mathematical Physics

Volume 2018, Article ID 4628346, 19 pages

https://doi.org/10.1155/2018/4628346

## The Inverse Problem of Identification of Hydrogen Permeability Model

Institute of Applied Mathematical Research of the Karelian Research Centre of the Russian Academy of Sciences, 11 Pushkinskaya St., Petrozavodsk 185910, Russia

Correspondence should be addressed to Yury V. Zaika; ur.ailerak.crk@akiaz

Received 7 December 2017; Accepted 2 May 2018; Published 7 June 2018

Academic Editor: Alexander Iomin

Copyright © 2018 Yury V. Zaika et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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

#### 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 ( is the plate thickness and is its area). The temperature is constant throughout the experiment. The concentration of dissolved diffusing hydrogen (in monatomic state) is sufficiently low and the diffusion flux can be considered proportional to the concentration gradient. The membrane is thin and the material has a sufficiently high hydrogen permeability coefficient, so we restrict ourselves to a standard diffusion equation:where is the concentration of diffusing hydrogen. The diffusion coefficient depends on the sample temperature in an Arrhenius way with preexponential factor and activation energy .

Initial data are determined by the fact that the sample had been preliminarily degassed: , , .

Nonlinear boundary conditions are derived from the material flux balance:Here, , are the amounts of hydrogen atoms in the input chamber of volume and output chamber of volume , , . The identity sign is frequently used here in the sense of equality by definition. Within the considered operating temperature range the gaseous hydrogen is in molecular form, but for consistency (considering that atomic hydrogen diffuses through the metal membrane) we use atoms as the unit. According to the kinetic gas theory, the incident particle flux density is related to the pressure by the Hertz Knudsen formula: ( is the Boltzmann constant, is the mass of a hydrogen molecule). In the context of the experiment it is convenient to choose the following measurement units , . Then we numerically obtain the dependence , (, ). 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 . One can write (as a parameter of the model) instead of , but it is more convenient to interpret the dimensionless probability factor as the fraction of absorbed H atoms within the notation. Thus, 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, are the densities of the desorption flux from the sample (deviation from the square desorption law is significant only at extreme temperatures) and is the desorption coefficient. We also assume that and depend on the temperature in an Arrhenius way. Formally, the activation energy in the exponent can as well be negative, being a linear combination of the activation energies and heats of the surface processes on the way “from gas to the solution”. If a constant saturation pressure of molecular hydrogen is maintained at a constant temperature on both sides of the membrane, the equilibrium concentration of the dissolved atomic diffusionally mobile hydrogen is finally established. By equating the derivatives in model (2) to zero, we get Sieverts’ law , .

Let us clarify the experimental conditions. The volumes comprise several liters, the thickness of the membrane is , the area is about , and the inflow pressure is within several hundreds .

It now remains to find the magnitudes of , . Within the time of transfer through the membrane the gas is in the thermodynamical quasi-equilibrium, wherefore we use the formula . Here, is the number of gas particles occupying the volume at the temperature and the pressure (in the SI system , , ). Taking into account the relationships and (formally), we get the following expressions for the corresponding pressures and volumes in the boundary conditions (2): , . Here, , , are the numerical values in the selected units (, , ).

Within the experimental unit the membrane is situated in a tube (which is heated to a predetermined temperature) between the inlet and the outlet chambers. The tube diameter is large enough to consider the equality of pressures as the criterion of thermodynamical quasi-equilibrium between the gas in the tube and in the chambers. The membrane temperature should be taken for the formula for the kinetic constant . The gas inside the volumes (whose massive walls are at room temperature) may get heated up. During the preliminary experiment it is recommended to fill the chambers with a practically impermeable metal membrane between them with “ambient” gas. Then we heat the tube and record the pressure rises. Within the framework of the ideal gas approximation (equation of state) this procedure enables estimation of the increments of the gas temperature inside the chambers. The corresponding gas temperatures are the ones to be used in the formula given for (and the subsequent ones, only excluding the value ). The need of such a refinement arises from the characteristics of this particular experimental unit. Such an adjustment of the values of should not cause difficulties in further calculations. Besides, this procedure has relatively little effect on the final calculation of the model pressures taking into account the measurement errors and relatively large volumes .

##### 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: Let us formulate the problem of modelling the concentrations using the pressures (the problem is also of interest per se) without the quasi-equilibrium simplification . The quasi-stationary state is achieved within a time , which is short compared to the total experiment time (). So, the original model (1)-(2) can be simplified (taking into account the formula ):Since by virtue of the “inlet-outlet” balance the equalities hold true, it is sufficient to express the concentrations from the boundary conditions (5) and substitute them into the first equation of (4) (the sign is chosen depending on whether the index is or ). The dimensionless normalized variables are convenient for numerical simulation:In addition, the system of equations (5) is compactly written in the symmetric form , . For the variable we obtain the incomplete quartic equation , which can be solved in radicals (for physical considerations we are interested in the positive root). However, the explicit expression is somewhat cumbersome and we will anyway have to numerically integrate the first equation of (4) in the form . Therefore, we shall aim to derive differential equations for , since information about the dynamics of the boundary concentrations is also important.

Differentiate (5) with respect to time and substitute the pressure derivatives from (4). For the variables we get the system

The change of variables in (7) determines the concentrations from which the model pressures , are calculated using (5).

Let us formulate step-by-step the numerical algorithm of modelling the pressures for the current values of , , coefficients (the authors used the Scilab freeware). We target at the “normal” experimental conditions [11, 12, 29–31], including the values of , , , , .

(1) We fix : omit fast transient processes (the duration of the transient processes is about tens of seconds on the hours-long experimental time scale). For the variable we choose the root of the biquadratic polynomial . From physical considerations it follows that and thus .

(2) The system of equations , () yields the missing value of . Formally, one equation is enough, but we take into account averaging procedures including determination of the values of .

(3) We numerically integrate the ODE system (8), (9) with the obtained initial data. The change of variables in (7) defines the concentrations , which are used to calculate the model pressures , from (5).

Computational experiments show that the model curves almost coincide (at ) with those generated by the originally proposed model, i.e., the nonlinear distributed initial boundary value problem.

Observe the difference from the quasi-equilibrium model (the Richardson approximation), where the only parameter for approximation of the experimental pressure is the complex . All the variable parameters of the original model that influence the permeability, namely, , , , are important when running the above algorithm. Thus, the fast hydrogen permeability model does not lose in informativeness concerning the considered transfer parameters.

#### 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 , , , , , . We set the value and calculate the corresponding desorption coefficient . The input pressures 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 . 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): . Within the model we determine and by the formulas For the given values of the parameters we have , , .

Next we express the solutions of the standard boundary value problems with Dirichlet boundary conditions corresponding to the jumps of the inlet pressure.

*Stage I*. The boundary value problem of the penetration method is The penetrating flux is .

From the computational point of view it is convenient to introduce the dimensionless time which is oriented at the characteristic diffusion time . As small a singularity appears if we directly use the partial sum of the expressionLet us provide another expression for using the properties of the Jacobi theta function. More precisely, we are interested in the function We have an alternative presentation for [34]:The series on the left is rapidly converging for large . But the series on the right is rapidly converging for small . If we define (), then we get , or , which is known as the functional equation for the theta function [35]. After some auxiliary transformations for we get the appropriate representation

The function has an -shaped form of the saturation curve (see the inset in Figure 1). Stage I ends with .