Trefftz Method of Solving a 1D Coupled Thermoelasticity Problem for One- and Two-Layered Media

: This paper discusses a 1D one-dimensional mathematical model for the thermoelasticity problem in a two-layer plate. Basic equations in dimensionless form contain both temperature and displacement. General solutions of homogeneous equations (displacement and temperature equations) are assumed to be a linear combination of Trefftz functions. Particular solutions of these equations are then expressed with appropriately constructed sums of derivatives of general solutions. Next, the inverse operators to those appearing in homogeneous equations are deﬁned and applied to the right-hand sides of inhomogeneous equations. Thus, two systems of functions are obtained, satisfying strictly a fully coupled system of equations. To determine the unknown coefﬁcients of these linear combinations, a functional is constructed that describes the error of meeting the initial and boundary conditions by approximate solutions. The minimization of the functional leads to an approximate solution to the problem under consideration. The solutions for one layer and for a two-layer plate are graphically presented and analyzed, illustrating the possible application of the method. Our results show that increasing the number of Trefftz functions leads to the reduction of differences between successive approximations.


Introduction
The multilayered materials have several unique characteristics that are frequently required and essential in structural design problems. To design an efficient composite structure, a detailed analysis must be carried out to check the stresses at the layer interfaces, where high local stresses may occur.
Various approaches and methods have been presented in the literature for the stress analysis of a multilayered medium. The authors of [1,2] presented the transfer matrix method for the isothermal elasticity problem that has been developed. In the paper [3] this approach has been extended to thermoelasticity problems. The authors of [4] employ the flexibility and stiffness matrix methods. To obtain the complete solution of the entire layered medium, this method has also been adopted by the authors of [5]. A transfer matrix formulation has been developed by the authors of [6] for longitudinal wave component propagation in both continuum and layered structures. The authors of [7] coupled thermoelastic problems for multilayered spheres are considered. Using the Laplace transform with respect to time, the general solutions of the governing equations are obtained in the transform domain. The solution is obtained using the matrix similarity transformation and inverse Laplace transform. The behavior of thermoelastic waves at the interface of layered medium and distributions of these waves through the domain is examined by applying the direct finite element method by the authors of [8]. The authors of [9] discuss the thermoelastic effects in the deformation of plates with arbitrary changing elastic parameters and temperature through-thickness. Using the semi-inverse method, a simple analytic solution is obtained for a thermoelastic problem of a nonhomogeneous plate with arbitrary contour. An analysis of displacements and stresses of FGM thick spherical pressure vessels with exponential varying material properties using the semi-analytical solution has been considered by the authors of [10]. In the paper [11] coupled thermoelasticity problems for a "coating-substrate" system are solved with the use made of the variational principle of coupled thermoelasticity in the Laplace transforms space and FEM. The transform inversion is realized according to the Durbin method [12]. The authors of [13] show that the exact solution for the temperature field has been obtained by separating variables and splitting the problem into two parts homogeneous transient and nonhomogeneous steady state. The set of equations obtained are solved using the rigorous applications of analytic techniques with the help of the eigen value expansion method. However, the thermoelastic response is studied in the context of uncoupled thermoelasticity. The authors of [14] study the mechanical behavior of two linear isotropic thermoelastic solids, bonded together by a thin layer, constituting a linear isotropic thermoelastic material, using an asymptotic analysis based on a finite element approach.
This work aims to obtain an approximate solution for predicting the temperature, displacement, and stress fields inside one-and two-layered thermoelastic media based on the coupled thermoelasticity equations. General solutions of homogeneous equations (displacement and temperature equations) are assumed to be a linear combination of Trefftz functions (T-functions). Particular solutions of these equations are then expressed by appropriately constructed sums of derivatives of general solutions, which is a novelty in solving the coupled equations. Then, the inverse operators to those appearing in homogeneous equations are defined. Applying them to the right-hand sides of the equations makes it possible to obtain such a form of solutions that can be used to determine the coefficients for T-functions.
The Trefftz method is widely described in many papers and monographs, and articles [15][16][17][18][19][20][21][22][23][24]. The T-functions for a linear partial differential equation usually can be obtained using the method of variables separation. Expanding the solution(s), called generating function(s), into power series with respect to the parameter(s) resulting from the method of variable separation, one arrives at the T-functions. In 2000, two new methods of deriving the T-functions were published. The one based on expanding function in the Taylor series is particularly simple and effective; the second is based on the so-called inverse operators [25][26][27]. To determine the unknown coefficients of these linear combinations, a functional is constructed that describes the error of meeting the initial and boundary conditions by approximate solutions. The minimization of the functional leads to an approximate solution to the problem under consideration.

Heat Conduction Equation
Consider a perfectly conducting elastic infinite isotropic homogeneous medium without any heat sources or body forces. The relationship between the heat flux vector and the thermal disturbance in the thermal wave model is described as, ( * ) ,i denotes differentiation with respect to the Cartesian coordinate x i , i = 1, 2, 3 and u i , i = 1, 2, 3, are the components of the displacement vector. The thermal conductivity, k, will be considered as a constant. Finally, the heat equation takes the form, with a = k ρc , [m 2 /s], standing for thermal diffusivity.

Equations of Motion and Constitutive Equations
The equations of motion have the form, Define components of strain tensor as follows: Then the constitutive equations take the form, σ ij are the components of stress tensor and δ ij is the Kronecker's symbol. Moreover, ε kk ≡ e = ε 11 + ε 22 + ε 33 is dilatation.

Formulation of the Problem
The coordinate system is so chosen that the x-axis is taken perpendicularly to the considered layers, and the yand z-axes in parallel ( Figure 1). Assume that the y-axis is the symmetry axis of a layer above and below the y axis with the same boundary and initial conditions. Then, for x = 0, the temperature and displacement conditions will follow the assumed symmetry. In further considerations, the layers will be numbered as 1 and 2.

Equations of Motion and Constitutive Equations
The equations of motion have the form, Define components of strain tensor as follows: Then the constitutive equations take the form, are the components of stress tensor and is the Kronecker's symbol. Moreover, ≡ = + + is dilatation.

Formulation of the Problem
The coordinate system is so chosen that the x-axis is taken perpendicularly to the considered layers, and the y-and z-axes in parallel ( Figure 1). Assume that the y-axis is the symmetry axis of a layer above and below the y axis with the same boundary and initial conditions. Then, for x = 0, the temperature and displacement conditions will follow the assumed symmetry. In further considerations, the layers will be numbered as 1 and 2. The problem is described by a system of two coupled equations:  The problem is described by a system of two coupled equations: For 1D problem with u ≡ u 1 and e = ∂u ∂x . Moreover, σ ≡ σ 11 and the rest of the stress components are zero. Let us introduce the substitutions: Hence, Equations (8) and (9) take the form, The constitutive equation takes then the form, By virtue of (1.1), in 1D, the relationship between the heat flux and the thermal disturbance can be expressed as, In further consideration, the primes are dropped for convenience.

Case I: One Layer
Consider Equations (11) and (12) and denote the dimensionless length l 1 as l = Kl 1 . The following dimensionless initial and boundary conditions are assumed: Here, the dimensionless length l 1 is denoted as l = l 1 a λ+2µ ρ and BC(t) stands for the boundary condition for θ (l, t) (to be shown later).

Case II: Two Layers
Equations (11) and (12) hold for layer 1. The initial, boundary, and symmetry conditions have the same form in dimensional and dimensionless coordinates. The layers will be numbered as 1 and 2. The conditions on the border of the layers (i.e., for x = l 1 ) in dimensional form are as follows: In dimensionless form, we can find l 1 ≡ l 1 1 a λ+2µ ρ 1 = l 1 K 1 as an upper limit of layer 1 and l 1 ≡ l 1 K 2 as a lower limit of layer 2. Hence, l 1 = l 1 The length l 2 can be expressed dimensionless by the thermophysical parameters of layer 1, l 2 = l 2 K 1 and of layer 2, l 2 = l 2 K 2 . Then, l 2 = l 2 Thus, the conditions (18) in dimensionless form read (dropping the primes for convenience), with R = γ (λ+2µ) , and t ∈ (0, t end ). Finally, (19) and all relations below refer to the dimensionless layer 1 coordinates.
Initial conditions for x ∈ 0, l 1 (in layer 1) and x ∈ l 1 , l 1 + l 2 (in layer 2), i.e., x ∈ l 1 , take the form (with the primes dropped for convenience): Boundary conditions on the upper border of layer 2 in the dimensionless form read: with dimensionless t as for layer 1.
With the y axis anticipated as the symmetry axis, the following conditions are obtained for x = 0: Energies 2021, 14, 3637 6 of 16

Trefftz Functions for Cases I and II
In order to find an approximate solution of the problems presented as Cases I and II, Trefftz functions (also called T-functions) for homogeneous Equations (11) and (12), i.e., for the equations (24) can be easily derived from the harmonic functions putting in harmonic equation y = it, i = √ −1 being the imaginary unit. Hence, one obtains, The T-functions F 1n and F 2n , n = 0, 1, . . ., are then renumbered to form a string {F m } m=0,1,...,M , M standing for number of functions. The approximate general solution of Equation (24), u b , is a linear combination of T-functions: T-functions for Equation (25) are known [18], and read, The approximate general solution of Equation (25), θ b , is a linear combination of the T-functions: d m , m = 0, . . . , M and f r , r = 0, . . . , R stand for coefficients of linear combinations (27) and (29) (to be calculated). (24) and (25) The inverse operator L −1 U for Equation (24) can be determined as follows,

Inverse Operator for Operators
The inverse operator L −1 θ for Equation (25) can be determined as, It is easy to check that in both cases LL −1 x k t m = x k t m . For all inverse operators The T-functions and inverse operators for the case II are the same as for case I.

Solution of Nonhomogeneous System of Equations
The approximate solutions of Equations (11) and (12) will be expressed as a sum of the general solutions (27) and (29) and the particular solutions received by the operation of inverse operators, (30) and (31), on the right sides of Equations (11) and (12). Hence, with u p and θ p standing for particular solutions of Equations (11) and (12). Of course, Denote, But, Hence, Assuming that in the above presentation of u ap the number of products of the operators A, B, C, D is equal to n, the function u ap can be rewritten as, Similar operation applied to θ ap leads to the following result: Hence, In the case of two layers, we do the same, remembering about different parameters for layers I and II.
It should be clearly emphasized that the solution is given by Formulas (36) and (37) fulfill exactly the system of Equations (11) and (12). For example, taking three Trefftz functions in Formula (27) and four Trefftz functions in Formula (29), we get, It is easy to show that Functions (38) fulfill the system of Equations (11) and (12). The initial and boundary conditions are fulfilled approximately (see Section 5).

The Objective Functional for the Case I
The objective functional, having a form J x, t, c n f n , {d m } (in the case I denoted as J 1 ), which describes the mean square error of fulfilling the initial, the boundary, and the internal conditions, is defined as follows: dt for conds (16) + (u ap (0, t)) 2 dt for conds (17) As BC(t), we take for calculations, Finally, the objective functional J 1 reads, Figure 2 shows the boundary condition (39) for t ∈ (0, 0.4).
Energies 2021, 14, x FOR PEER REVIEW 8 of 17 In the case of two layers, we do the same, remembering about different parameters for layers I and II.
It should be clearly emphasized that the solution is given by Formulas (36) and (37) fulfill exactly the system of Equations (11) and (12). For example, taking three Trefftz functions in Formula (27) It is easy to show that Functions (38) fulfill the system of Equations (11) and (12). The initial and boundary conditions are fulfilled approximately (see Section 5).

The Objective Functional for the Case I
The objective functional, having a form , , , (in the case I denoted as ), which describes the mean square error of fulfilling the initial, the boundary, and the internal conditions, is defined as follows:

The Objective Functional for the Case II
The objective functional ( , , , , , ), which describes the mean square error of fulfilling the initial, the boundary, and the internal conditions, consists of

Numerical Results and Discussion
For numerical computations, the physical constants of stainless steel and Lead Zirconate titanate (PZT-5A) [28], are considered as in Table 1. The computations were carried out for dimensionless length l 1 = 1 and l 2 = 0.1 the numerical values of the temperature, stress, and displacement are represented graphically. The lengths l 1 and l 2 are expressed dimensionless by the thermophysical parameters of layer 1. The calculations for one layer with dimensionless length l 1 = 1 are shown in the graphs presented in Figure 3. They were performed with the approximation of the solutions by a 34degree polynomial (34 T-functions for Equation (25) and 67 T-functions for Equation (24)).

Results for Case I
The calculations for one layer with dimensionless length l1 = 1 are shown in the graphs presented in Figure 3. They were performed with the approximation of the solutions by a 34-degree polynomial (34 T-functions for Equation (25) and 67 T-functions for Equation (24)).
In order to investigate whether the process of calculating temperature, displacement, and stress is convergent, the maximum norm has been examined, for successive differences of solutions obtained for degree 10 (2) 34 (i.e., 10, 12, 14, …, 34) polynomials which approximate solutions. The maximum norm was investigated for differences in solutions in the space-time domain for = ( , ) ∈ = (0, 1) × (0, 0.4) ( Table  2) and for an approximation of the boundary condition for l = 1 for ∈ (0, 0.4) ( Table  3).   In order to investigate whether the process of calculating temperature, displacement, and stress is convergent, the maximum norm has been examined, for successive differences of solutions obtained for degree 10 (2) 34 (i.e., 10, 12, 14, . . . , 34) polynomials which approximate solutions. The maximum norm was investigated for differences in solutions in the space-time domain for z = (x, t) ∈ Z = (0, 1) × (0, 0.4) ( Table 2) and for an approximation of the boundary condition for l 1 = 1 for t ∈ (0, 0.4) ( Table 3). Table 2. The maximum norm for differences (upper row) for boundary condition BC and for temperature, displacement, and stress in the space-time domain.  Table 3. The maximum norm for the difference between the boundary condition BC for temperature for l 1 = 1 and its approximation. The data presented in Tables 2 and 3 show that the process of approximating solutions for temperature, displacement, and stress, as well as the process of approximating the boundary condition for temperature for l 1 = 1 are convergent. Figure 3 shows that the plot for the stress field is waving strongly. This is due to the differentiation of displacements, cf. (13). The displacement field shows slight wavings; they strengthen for stresses. This is a characteristic feature of solutions described by a linear combination of polynomials. This phenomenon disappears as the number of polynomials increases. However, increasing the number of polynomials significantly increases computation time. Since the calculations were performed on a PC, this indicates the need to use a much faster computer. It is worth mentioning that the T-functions for each of the equations constitute a complete system of functions [18,27].

Results for Case II
Calculations for two layers with dimensionless lengths l 1 = 1, l 2 = 0.1 are shown in the graphs in Figure 4. There are presented graphs for 11th and 25th degree polynomials. The comparison of the graphs shows how the degree of the used polynomial (combination of T-functions) affects the accuracy of the solution.
The graphs in Figure 4 show that in the case of temperature (top row), the number of polynomials slightly affects the shape of the graph. On the other hand, with higher degrees of the polynomial, both the displacement graphs (middle row), and stress graphs (bottom row) smooth out.
In order to investigate whether the process of calculating temperatures, displacements, and stresses in each of the two layers is convergent, the maximum norm was again tested for successive differences of solutions obtained for 9, 11, 13, . . . , 25 for layer 1 (for layer 2 it was 9, 13, 17, . . . , 25) T-functions for Equation (25) and a correspondingly doubled number of T-functions for Equation (24) in both layers. The maximum norm was investigated for the differences in solutions in the space-time domain for z = (x, t) ∈ Z = (0, 1) × (0, 0.4). The data presented in Tables 4 and 5 show that the process of approximating solutions for temperature, displacement, and stress in each layer separately are convergent processes.  In order to investigate whether the process of calculating temperatures, displacements, and stresses in each of the two layers is convergent, the maximum norm was again tested for successive differences of solutions obtained for 9, 11, 13, …, 25 for layer 1 (for layer 2 it was 9, 13, 17, …, 25) T-functions for Equation (25) and a correspondingly doubled number of T-functions for Equation (24) in both layers. The maximum norm was investigated for the differences in solutions in the space-time domain for = ( , ) ∈ =  The influence of the thermal conductivity of layer 2 and its thickness on the distribution of temperatures, displacements, and stresses in the layers was also investigated. The results are presented in the graphs in  (0, 1) × (0, 0.4). The data presented in Tables 4 and 5 show that the process of approximating solutions for temperature, displacement, and stress in each layer separately are convergent processes. The influence of the thermal conductivity of layer 2 and its thickness on the distribution of temperatures, displacements, and stresses in the layers was also investigated. The results are presented in the graphs in Figures 5-7.
The comparison of the temperature, displacement, and stress graphs in Figure 4 with the graphs in Figure 5 shows much higher temperature values in the first layer, resulting from the greater conductivity in the second layer, and higher values of displacements with reduced stress values.      Table 1, the thickness of layer 2 is two times smaller. Figure 6 shows a significant influence of the thickness of layer 2 and its thermal conductivity coefficient on the values of temperature and displacements in layer 2 at significantly reduced stresses in this layer. Figure 7 shows that reducing the thickness of layer 2 Figure 7. The thermal conductivity coefficient for layer 2 as in Table 1, the thickness of layer 2 is two times smaller. The comparison of the temperature, displacement, and stress graphs in Figure 4 with the graphs in Figure 5 shows much higher temperature values in the first layer, resulting from the greater conductivity in the second layer, and higher values of displacements with reduced stress values. Figure 6 shows a significant influence of the thickness of layer 2 and its thermal conductivity coefficient on the values of temperature and displacements in layer 2 at significantly reduced stresses in this layer. Figure 7 shows that reducing the thickness of layer 2 slightly increases the temperature in layer 1, slightly reduces the displacements on the outer edge of layer 1, and slightly disturbs the stresses.

Conclusions
The application of the Trefftz functions to the approximate solution of the thermoelastic problem in a single-and double-layered medium showed the effectiveness of the Trefftz method in connection with the method introduced in the article for determining particular solutions of coupled thermoelasticity equations. However, the solution procedure requires a powerful computer that will approximate solutions with high-order polynomials. The temperature and displacement fields in one and two layers could be obtained in a satisfactory form for a small number of T-functions (even for the solutions to be approximated by a 10-degree polynomial). The calculation of stresses requires a much higher degree of the approximating polynomial, due to the differentiation in Formula (13), which introduces significant undulations in the stress diagram. The applied method obtains an analytical solution with numerically determined coefficients at the Trefftz functions.
As composite materials are used on a larger scale, a rationale for research is the effect of sudden changes in their thermophysical properties [29,30]. The presented problem of thermo-structural analysis is part of this research area.
It should be emphasized that the approximate solution (in the form of a linear combination of the Trefftz function), which strictly fulfilled the fully coupled system of thermoelasticity equations is presented in the paper. It is probably the first time that such a solution has been obtained in this manner.