Thermoelastic Diffusion Multicomponent Half-Space under the Effect of Surface and Bulk Unsteady Perturbations

This article presents an algorithm for solving the unsteady problem of one-dimensional coupled thermoelastic diffusion perturbations propagation in a multicomponent isotropic half-space, as a result of surface and bulk external effects. One-dimensional physico-mechanical processes, in a continuum, have been described by a local-equilibrium model, which included the coupled linear equations of an elastic medium motion, heat transfer, and mass transfer. The unknown functions of displacement, temperature, and concentration increments were sought in the integral form, which was a convolution of the surface and bulk Green’s functions and external effects functions. The Laplace transform on time and the Fourier sine and cosine transforms on the coordinate were used to find the Green’s functions. The obtained Green’s functions was analyzed. Test calculations were performed on the examples of some technological processes.


Introduction
Creation of materials with pre-programmed properties and modification of existing materials under new quality standards is one of the development of modern areas, in the construction materials processing technological field.The mathematical model construction of the various factors impact the material being processed and the obtained result analyses allow us to consider a greater number of modification options, while minimizing financial costs in the case of complex and high-tech processes.
There are a number of approaches to mathematical model constructions.The approach to coupled fields model construction is one that is favorable, is actively developing, and gives the opportunity to most accurately and analytically describe complex dynamic technological and physical processes.An example of such a model is a thermomechanical diffusion model, which describes the interaction of temperature, displacement, and concentration fields.
The relevance of this research direction is confirmed by the presence of the large number of works in this study area, worldwide.Among the main studies, we can identify [1][2][3][4][5] as the source of the problem formulation.
Most of the currently available work are devoted to solving static [6], quasistatic [7], and stationary [8] problems of thermomechanical diffusion.However, unsteady coupled one-dimensional [9] and two-dimensional problems [10,11] are the most interesting and difficult.In these articles, the solution of the unsteady problem has been reduced to the Laplace transform on time.Its inversion is associated Math.Comput.Appl.2019, 24, 26 2 of 14 with great mathematical difficulties.Numerical algorithms, ready-made computational mathematics and mechanics packages have been used most often, in this case, to get the originals of the unknown functions [12].In addition to these studies, it should be mentioned that a sufficiently detailed review on various issues of thermomechanical diffusion processes modeling for the 20th century, is available in [13].It is important to note that the very few studies devoted to the unsteady multicomponent thermoelastic diffusion continuum.
In the represented article, we consider the one-dimensional unsteady thermoelastic diffusion problem for homogeneous isotropic multicomponent half-space.The linearized local-equilibrium model of coupled thermoelastic diffusion has been used to describe the medium perturbations propagation with finite speed.The model includes the equations of elastic medium motion, heat transfer, and mass transfer [2,4,5,14].Initial conditions are zero.
The solution of the problem, similar to problems discussed in [15][16][17][18][19][20], is sought in the integral convolution form of Green's functions and functions of bulk or surface external perturbations.The Laplace transform on time and the Fourier sine and cosine transform on the coordinate are used to find the Green's functions.As a result of the transformations, the problem is reduced to a system of linear algebraic equations.Its solution is represented as rational fractions, which depend on the Laplace transform parameter.Their originals are found using well-known theorems and tables of operational calculus.Moreover, there is no need to develop new complex approaches for the numerical inversion of the Laplace transform, as it has been done in the above-mentioned works [9][10][11].The inversion of the Fourier sine and cosine transforms, in general, can be done numerically, using such well-known approaches as the Filon's method or the inverse Discrete Fourier Transform (DFT).Such an algorithm for the problem makes it possible to minimize the use of numerical algorithms and enables the analyses of the Green's functions.
The unsteady problems of thermomechanical diffusion (as opposed to steady ones) make it possible to consider rapid, complex, and coupled dynamic technological processing of metals and alloys, accompanied by an explosive growth of loads and intensive mass and heat transfer [4,14,21,22].Examples of such processes can be various pulse-periodic methods of surface treatment, such as laser [23] or electro-discharge in an oncoming liquid flow [24].In addition to this, it is necessary to add that one of the processes that can be quite fully described by the presented mathematical model is the process of ion implantation.This technology makes it possible to obtain modified coatings, with exceptional characteristics, and create many different stable solid solutions [25,26].

Problem Formulation
We consider an isotropic, homogeneous, N-component half-space bounded by a surface x 1 = 0 (Ox 1 x 2 x 3 is a rectangular Cartesian coordinate system, the axis Ox 1 is directed deep into the half-space).The dimensionless local equilibrium linear model of coupled thermoelastic diffusion is used to describe the medium perturbations propagation with finite speed and cross-diffusion effects [2,4,5,16,18,21].This model includes the equation of elastic medium motion, the equation of heat transfer, and N equations of mass transfer (prime and point, respectively denotes derivatives by the dimensionless spatial variable x and by the dimensionless time τ) [2,4,5] Math.Comput.Appl.2019, 24, 26 3 of 14 Displacements, thermal, and diffusion flows are set at the boundary of the half-space [5]: The limitation conditions are set at infinity (x → ∞) [5]: Initial conditions are zero: In Equations ( 1)-( 4) and further, dimensionless quantities are used (if a letter coincides, then the dimensional quantity is indicated with an asterisk; q = 1, N): Here t is the time variable; u 1 is the displacement vector component; L is the characteristic length of half-space; q is the component number; n (q) 0 and n (q) are initial and actual concentrations (mass fractions); t ϑ is the thermal relaxation time; t (q) is the diffusion relaxation time; C 1111 is the elastic constant; ρ is the mass density; b 11 is a temperature constant characterizing thermal deformations; 11 is a coefficient characterizing the medium volumetric change due to diffusion; D (qp) 11 is the diffusion coefficient; m (q) is the molar mass; R is the universal gas constant; T 0 and T are initial and actual temperatures; κ 11 is the coefficient of thermal conductivity; γ (q) is the activity coefficient; c 0 is the specific heat at constant concentration and deformation; F * 1 (x 1 , t) is the body force; F * 2 (x 1 , t) is the bulk heat inflow; and F * q+2 (x 1 , t) is the bulk mass inflow.

Integral Representation of the Solution
We present the solution of Equations ( 1)-( 4) in the form of convolutions [15][16][17][18][19][20]: Here G ik = G ik (x, ξ, τ) i = 1, N + 2 are the bulk and G ik = G ik (x, τ) are the surface Green's functions of Equations ( 1)-( 4) i = 1, N + 2 .The convolutions have the form: Bulk Green's functions are solutions of problems that include equations similar to system of Equation ( 1) [27]: Bulk Green's functions problem also includes zero initial conditions like Equation ( 4) and the following boundary conditions: Here and below, δ(x − ξ) and δ(τ) are the Dirac delta-function, δ kl is the Kronecker symbol.Surface Green's functions are solutions of problems that include homogeneous equations similar to system of Equation (1) [27]: Surface Green's functions problem also include zero initial conditions like Equation ( 4) and the following boundary conditions: Thus, the solution of Problem ( 1)-( 4) is reduced to finding the Green's functions.First of all, we consider an algorithm for finding the surface Green's functions.On the basis of this, we describe an algorithm for the bulk Green's functions.

Algorithm for Finding the Surface Green's Functions
We apply the Laplace transform on time to the homogeneous system of Equation ( 10) and the boundary conditions Equation (11), and get a system of differential equations [17][18][19][20] («s» is the transformation parameter, the index «L» denotes the Laplace transformant): Further, we apply to system of Equation ( 12) the Fourier sine and cosine transform («F» denotes the Fourier transform) [28]: Math.Comput.Appl.2019, 24, 26 6 of 14 To do this, we multiply the first equation in system of Equation ( 12) by sin(λx) and the rest by cos(λx).
Then, we integrate it along Ox in the interval from 0 to +∞ and take into account Equation (13).As a result, we obtain the system of linear algebraic equations for the Fourier-Laplace images [17,20]: where The solution of System (15) has the form of rational fractions: Here P(λ, s) is the determinant of the homogeneous System (15): P ik (λ, s) are the determinants obtained from P(λ, s) by replacing the i-th column with the right-side column of the System (15): The above is done following Cramer's rule.Below are the polynomials P ik (λ, s) for the two-component half-space that we will use for the calculation examples: Let s j = s j (λ) ∈ C j = 1, 2N + 4 be simple zeros of the polynomial P(λ, s).Based on the well-known residue theorems, the originals of the surface Green's functions are written as follows [16,19]: where Here prime denotes derivative by parameter «s».
For the final expression of the unknown functions, the Green's functions found in this way are substituted into the convolutions Equation (6).The originals in Equation ( 19) can be calculated using well-known numerical algorithms, such as inverse DFT-type methods or the Filon's method.The quadrature formula on an irregular sparsening grid presented in [17] can also be used.

Algorithm for Finding the Bulk Green's Functions
By analogy to the algorithm for finding the surface Green's functions, we apply the Laplace transform on time and then the sine and cosine transform on the spatial coordinates to the inhomogeneous system of Equation ( 8) and boundary conditions Equation (9).Then, we get the following system of linear algebraic equations: Its solution has the form similar to Formula ( 16): Here P(λ, s) ≡ P(λ, s) is the determinant of the homogeneous system in Equation ( 21).P ik (λ, ξ, s) are the determinants obtained from Formula ( 17) by replacing the i-th column with the right-side column of the System ( 21): (δ 1k sin(λξ), δ 2k cos(λξ), δ 3k cos(λξ), . . . ,δ N+2,k cos(λξ) The above is done following Cramer's rule.It is important to note that the polynomials P ik (λ, ξ, s) can be written through the P ik (λ, s) for i = 1, N + 2 and k = 2, N + 2, as follows: Originals of the bulk Green's functions can be found by a method similar to system of Equation (19) and Formula (20).

Analysis of Singularity
In [20] it was shown that some of the surface Green's functions originals had singularities in the form of Dirac delta functions and Heaviside functions on time, when s → ∞ .If considering the finite propagation speed, there remains only one singularity of the function G 21 (x, τ), in the Dirac delta function form.In the domain of Fourier-Laplace images it has the form: This singularity can be distinguished as an integer part as follows: where When we substitute Equation ( 26) into Equation ( 19), the function G L 21 (x, s) takes the form [20]: In this case, its original can be presented as follows, similar to Equation ( 19): where The analysis of such singularity can help avoid errors in numerical integration and simplify the calculation of convolutions.

The Calculation Example on the Ion Implantation Technology
Ion implantation is a method of introducing impurity atoms into the surface of a substrate, by bombarding one with a high-energy ion beam [25,26].As an assumption, the full cycle of ion implantation technology can be divided into several stages-the formation of the ion beam with specified characteristics, the fall of the beam on the substrate surface (target bombing process), the saturation of the target material with insertion elements and, the last stage, material annealing to form a layer of various chemical compounds and reduce internal stresses after bombing.The proposed model of thermoelastic diffusion, in a first approximation, is capable of describing the sequential stages of bombing and annealing, under the following assumptions: 1.
The infinite half-space is processed evenly over the entire area.

2.
Characteristics of the ion beam during the entire implantation process are constant.

3.
The possible inhomogeneity of the ion beam is neglected.

4.
The beam consists of only one chemical component.

5.
Chemical and phase transformations are not considered.
For example, we choose the two-component problem (N = 2) on the implantation of copper (Cu, q = 2) into an aluminum substrate (Al, q = 1) at the initial temperature T 0 = 700 K.At the initial time, copper already has a mass fraction in the substrate medium, equal to 0.05.All the rest (0.95) is aluminum.The characteristic length of a half-space L = 10 −2 m.Then, the following dimensionless values, obtained using Formula (5) [22,29,30], will correspond to the substrate medium: We will consider the whole technological process inseparably and accept the switch time τ o = 6.37 • 10 6 (which corresponds to 10 3 seconds) from the implantation mode (τ < τ o ) to the annealing mode (τ ≥ τ o ).The switching occurs instantly.To finalize a calculated example, we assume that the half-space boundary [25,26] • Is always rigidly fixed: f 1 (τ) ≡ 0; • At the implantation stage has no heat inflow, and at the annealing stage the heat flow is At the implantation stage is under inflow of Cu, but at the annealing stage mass transfer with the environment is absent: Bulk and volumetric external perturbations are absent.
Below are the distribution graphs of thermoelastic diffusion perturbations in the medium, over depth x of the half-space.The graphs were obtained by substituting the above boundary conditions into convolution Equation (6).All convolutions were found analytical.Only the roots s j of P(λ, s) and the inverse sine and cosine transforms were calculated numerically.For inverse transforms we used the quadrature formula on an irregular sparsening grid [17], with 10 5 integration points.This method is quite fast, but it has a slow convergence over the grid splitting.It should be used with good convergence of the integrals in Equation ( 19), over the transform parameter.The left column (Figures 1a, 2a, 3a and 4a) of the graphs shows implantation (τ < τ o ), the right one (Figures 1b, 2b, 3b and 4b) shows annealing (τ ≥ τ o ).

•
Is always rigidly fixed: ( ) • At the implantation stage has no heat inflow, and at the annealing stage the heat flow is ( ) ( ) • Bulk and volumetric external perturbations are absent.
Below are the distribution graphs of thermoelastic diffusion perturbations in the medium, over depth x of the half-space.The graphs were obtained by substituting the above boundary conditions into convolution Equation (6).All convolutions were found analytical.Only the roots j s of ( ) , P s λ and the inverse sine and cosine transforms were calculated numerically.For inverse transforms we used the quadrature formula on an irregular sparsening grid [17], with        In Figures 1-4, τ = 0.5τ o is indicated by the solid line; τ = 0.9τ o is indicated by the dotted line; τ = 2τ o is indicated by the dashed line; τ = 3τ o is indicated by the dash-dotted line.

The Calculation Example on the Pulse-Periodic Processes
Let us consider some general pulse-periodic method of surface treatment.They can include technologies such as laser modification [23] or electro-discharge modification in an oncoming liquid flow [24].They can be described as a first approximation of the above mathematical model, if we accept the following assumptions: 1.
The infinite half-space is processed, and the impact depth is significantly less than the surface being treated.

2.
A single impulse is considered.

3.
Possible non-uniformity of the material is neglected.

4.
Chemical and phase transformations, as well as electromagnetic fields and electron-phonon interaction are not taken into account.
For brevity and simplicity, we take the duralumin from the first calculation example with L = 10 −3 m and at T 0 = 600 K: To build the calculation example, we assume that deformations, thermal, and diffusion flows are given at the boundary of the half-space [23,24]: Of all the bulk effects, only the bulk heat inflow is set: Below in Figures 5 and 6 are distribution graphs of thermoelastic diffusion perturbations in the half-space.The graphs were obtained by substituting the boundary conditions above into convolution Equation ( 6).All convolutions were found analytically.Inverse transforms were done similar to the ion implantation example.It is necessary to add that the elastic problem solution obtained using the proposed algorithm, converges with one obtained analytically.

Conclusions
An algorithm is proposed for solving the one-dimensional unsteady thermoelastic diffusion problem for a multicomponent half-space.The main advantage of this approach is the ability to analytically find the Laplace transform originals of the surface and bulk Green's functions and analyze them.The originals of the Laplace transform were found analytically.It allowed us to reduce the final solution of the problem to the Fourier sine and cosine transform.The method effectiveness is demonstrated on the calculation examples technological processes.