Modelling of Generalised Thermoelastic Wave Propagation of Multilayer Material under Thermal Shock Behaviour

1National Center for International Joint Research of Micro-Nano Moulding Technology, School of Mechanics and Engineering Science of Zhengzhou University, Zhengzhou 450001, China 2The State Key Laboratory for Structural Analysis of Industrial Equipment, Department of Engineering Mechanics, Faculty of Vehicle Engineering and Mechanics, Dalian University of Technology, Dalian 116024, China 3School of Mechanics and Engineering Science, Zhengzhou University, Zhengzhou 450001, China


Introduction
Multilayer composite material structures are widely used to satisfy global mechanical and thermal requirements in the development of microelectronics system, photoelectric equipment, microsensors, and so on [1,2].Thermomechanical interaction is one of the major mechanical behaviours of the aforementioned structures.Heat action, especially heat shock, can cause a rapid temperature increase with considerable stresses in each layer.
In general, classical Fourier heat conduction is the conventional theory used extensively to study heat conduction in multilayer materials subjected to a low-or intermediate-frequency heat source.However, one limitation of the Fourier theory is its accuracy, which is deficient for extreme low-temperature or ultra-short-pulse thermal heating in temporal-spatial scale [3][4][5][6].As a result, non-Fourier heat conduction schemes have been proposed, the simplest of which is Cattaneo-Vernotte (C-V) approach [7,8].
In 1967, an alternative explanation of thermal wave disturbances was proposed by Lord and Shulman [9].The generalised thermoelastic theory, which is based on a new heat conductive law to replace the classical Fourier law, has been successfully used to explore the thermal shock problem.The essential feature of the Lord-Shulman (L-S) theory is that it contains one constant that acts as a relaxation time under the supposition that a certain time delay exists between the heat flux and the temperature gradient.
Traditional time integration methods (e.g., the Newmark method) [10][11][12], despite their popularity, fail to properly capture discontinuities of impulsive waves in space and produce spurious numerical oscillations in the simulation of thermal shock problem of high modes and sharp gradients.Numerical analysis of multilayer structures under impact thermal shock loading has attracted the attention of many researchers.In 2003, Li et al. presented a novel time-discontinuous Galerkin finite element method (DGFEM) for solving dynamics and wave propagation in nonlinear solids and saturated porous media [13].Wu and Li (2006) successfully used the DGFEM to simulate heat wave propagation [14].One of the main characteristics of the DG method is that the basic variables (temperature) together with their temporal derivatives are assumed to be discontinuous at a defined time point and to be interpolated separately in a time step, respectively.Such relaxation of restrictions on the continuity of the basic variables provides a mechanism to filter the spurious oscillations in the wave-after stage.
In this article, an additional artificial damping discontinuous Galerkin numerical algorithm for the generalised thermoelastic problems of multilayer materials is developed on the basis of the generalised L-S and Wu and Li theories [13].An additional stiffness proportional damping scheme (  C) is brought into the final finite element discretised form to reduce the numerical oscillations appearing in the wave-front stage and at the interface between two adjacent layers.Numerical results demonstrate clearly that the present DGFEM-  method is valid for filtering out spurious oscillations in both wave-front and wave-after stages and at the adjacent-layer interface and for providing more accurate solutions than those obtained using traditional time integration methods (e.g., the Newmark method).The present modelling method could also be extended to the study of generalised piezothermoelastic and magnetothermoelastic problems under high-frequency loading.

Theory of Linear Thermoelastic of Multilayer Materials
This section summarises the basic governing equations of linearized thermoelastic theory of multilayer materials.The basic thermoelastic coupled governing equations of layer  are expressed as [15]  ()  ()  () 0 where  is the temperature; U is the displacement vector;  is time;  () and  () are the Lame constants;  () is the density;  () is the specific heat at constant strain;  () =  () (3 () +2 () ), where  () is the coefficient of linear thermal expansion;  () is the thermal conductivity;  0 is the reference temperature; Q  is the force vector in the displacement field; Q  is heat source; and  () 0 is the thermal relaxation time.When  () 0 = 0, (1a)-(1b) reduce to the classical coupled thermoelastic theory.When  () = 0, (1a) reduces to the non-Fourier hyperbolic heat conductive equation.
Notably, the appropriate boundary conditions associated with the governing equations (4a)-(4b) must be adopted.When the temperature and displacements are prescribed on the surfaces Γ 1 and Γ 1 , respectively, then where  0 and u 0 are the prescribed values of the temperature and the displacement vector.However, if surface heat flux and surface traction are applied to the corresponding surfaces Γ 2 and Γ 2 , the following boundary conditions should be satisfied: where F   () and F   () are the given surface heat flux and the surface traction, respectively.
In the present paper, surface impulse heat source is applied with the form where  0 and  0 are the amplitude of thermal shock and the constant time.

Finite Element Discretisation in Spatial-Temporal Domain
3.1.Spatial Discretisation.The solution of the Galerkin counterpart of the weak formulation is expressed in terms of the shape functions and gives rise to the following system of equations for L-S theory: where Ω, +2     Ω, in which   and   are the shape function for the temperature and displacement, respectively.And   and   are the first-order spatial derivatives of the respective shape function.

Discontinuous Galerkin Discretisation in Time Domain.
Time integration of governing equation ( 7) is performed using a discontinuous Galerkin scheme, which permits the discontinuity of basic variables at the discrete-time sequence, 0 < ⋅ ⋅ ⋅ <   < ⋅ ⋅ ⋅ <   .We assume that the unknown variable is not continuous at   and is expressed as where ( +  ), ( −  ) are discontinuities of variables at   with ( ±  ) = lim →0 ± ( + ).
For each time subdomain   = ( +  ,  − +1 ), a time step length is defined Δ =  +1 −  , for  = 0, 1, . . ., −1.For each time subdomain   , we interpolate the unknown variables, that is, the global nodal temperature-displacements vector {  u }, by using the third-order Hermite time shape functions as where , and { q k } − +1 represent the global nodal values of the temperature-displacements vector and its time-derivative at times   ,  +1 , respectively.Among them, q = θ and k = u .The global nodal values of temperature-displacements vector and its time-derivative, that is, {  u } −  , { q k } −  , at time  −  can be obtained at the end of the previous time step.In addition, The global nodal values { q k } (the temporal derivative of the temperature-displacements vector) at arbitrary time  ∈ [  ,  +1 ] are treated as an independent variable by linear time shape functions as We define { Ũ} = {  u }, {V} = { q k } as the basic independent variables; { U} and { Ü} are then the first-order and secondorder temporal derivative, respectively, of { Ũ}.
By dropping the superscripts of the vectors Ũ+  , Ũ− +1 , V +  , V − +1 and the time variable , we can rewrite (10) and (12) as The generalised thermoelastic weak forms of the semidiscretised equation ( 9) and the constraint condition U − V = 0, along with the discontinuity conditions Ũ and V on a typical time subdomain can be expressed as Substituting ( 9) and ( 13) into ( 14), we obtain the following matrix equation of DGFEM for solving from independent variations of  Ũ+  ,  Ũ− +1 , V +  , and where This equation means that the continuity of nodal temperature-displacement vector {U} at any time level is automatically ensured.The discontinuity is located at nodal derivative vector {V} at discretised time levels.

Artificial Damping Scheme for DGFEM.
As previously reported [15], the selection of stiffness proportional and a mass proportional damping coefficient is effective for high-frequency oscillations and low-frequency oscillations, respectively.To filter out the oscillations in the wave-front stage caused by a high-frequency impulse load, an artificial stiffness proportional Rayleigh-type damping scheme is introduced with the following form: where [C  ]  is the stiffness proportional matrix,  ()  is the damping constants,  () is the damping ratio estimated from the global behaviour of the system, and   is the basic frequency of the considered structure.As a result of a number of trial runs, we also determined that the effective stiffness proportional damping constants  ()   should be less than or equal to the time step.Substituting the stiffness proportional artificial damping matrix equation ( 19) into (16), we obtain where

Numerical Results and Discussions
In this section, two numerical examples are investigated to show the reliability of the proposed algorithm.The first 1D example is based on hyperbolic thermal wave equations and is presented to demonstrate the ability of the proposed DGFEM to filter out the numerical oscillations in wave after, wave-front, and the interface between two layers.The second example of nonclassical thermoelastic based on the L-S model subjected to impulse heat load is investigated.
Simulated results in 1D and 2D are presented to demonstrate that the proposed DGFEM formulations are capable of producing reliable results in the generalised thermoelastic wave problem in multilayer structures.

1D Thermoelastic Wave Propagation Problem.
The first example is a 1D thermoelastic wave propagation problem, as shown in Figure 1, defined on the space-time domain [0, 0.25] × [0, 0.26].The associated dimensionless material constants are taken as The dimensionless thermal relaxation time is chosen as  () = 0.1,  = 1, 2, 3 [15].( To compare the reliability of the proposed DGFEM and traditional method, we evaluate the L-S model of  () = 0.926,  = 1, 2, 3 in this case.Figures 2 and 3   unknown variables play the significant roles in filtering out the spurious numerical oscillations and in providing much more accurate solutions for generalised thermoelastic coupled problem subjected to a thermal impulse load.

2D Coupled Thermoelastic Wave Propagation Problem.
In this problem, we consider a 2D nondimensional form of the thermal wave propagation problem on the basis of the L-S model in a domain with 31 × 31 four-node finite element elements under a plane strain assumption, as depicted in Figure 4.
The dimensionless material property data are defined as follows:  (1)  1 = 5.035×10 −3 ,  (1) = 0.926,  (1) = 0.1 [15],  (2) 1 = 5.035×10 −3 ,  (2) = 0.926,  (2) = 0.2,  (3)  1 = 5.035×10 −3 ,  (3) = 0.926,  (3) = 0.3.The impulsive thermal source is the same as that in the first example and is applied along the curve AB, as shown in Figure 4. Time step is chosen as Δ = 3.0 × 10 −3 , and an artificial damping constant  ()  = 1.5 × 10 −3 is used.Figures 5-6 show snapshots of the temperature and the stress distributions versus time for  = 0.06, 0.12, 0.18, and 0.24, as obtained by the proposed DGFEM and by the Newmark method.Because the effects of high modes cannot be filtered out, the spurious numerical oscillation is observed at the free end after the wave tail passes in the calculation by the Newmark method.A comparison of the temperature and the stress distributions between the proposed DGFEM-  and the Newmark method reveals that the proposed DGFEM can successfully describe the behaviour of the generalised thermoelastic coupled problem subjected to a thermal impulse load in the 2D case.

Conclusions
Multilayer materials are widely used in the development of microelectronics systems, photoelectric equipment, microsensors, and so on.In this article, we presented a numerical treatment of the generalised non-Fourier thermoelastic theory of a multilayered material subjected to a high-frequency impact heat source.The temperature-displacements vector, together with its first-order temporal derivative, was treated as independent basic unknown variables.The thirdorder Hermite time shape functions and linear time shape functions were applied to interpolate the temperaturedisplacement and its temporal derivative, respectively.An additional artificial damping discontinuous Galerkin numerical algorithm was incorporated into the final finite element discretised form to reduce the numerical oscillations in the wave-front stage and at the interface between two adjacent layers.The advantage of the present DGFEM-  was demonstrated through comparison with traditional time-stepping algorithms such as the Newmark  family algorithm via a series of numerical examples.The numerical results demonstrated that the present modelling is valid to filter out the spurious numerical oscillations in the wavefront and wave-after stages and at the interface between the layers.

Additional Points
Highlights. (1) A time-discontinuous Galerkin finite element method is presented for the coupled Lord-Shulman theory of multilayer materials under a high-frequency heat source.
(2) An artificial damping scheme is introduced into the finite element equations to filter out the spurious numerical oscillations.(3) Numerical simulations show that the present modelling effectively filters out the spurious numerical oscillations in both the wave-front and wave-after stages and at adjacent-layer interfaces.

Figure 1 :
Figure 1: A uniform thermoelastic configuration with one end fixed and the other end subjected to a rectangular impulse heat load.
represent the temperature and stress distributions along the distance of axis at different nondimensional time  = 0.16,  = 0.20,  = 0.22, and  = 0.26 by the proposed DGFEM and Newmark method.Serious spurious numerical oscillations could be observed in both temperature and stress profiles at different time levels by Newmark method.By contrast, the smooth and continuous temperature and stress distributions are exhibited by using DGFEM.This indicates that artificial damping constant and discontinuous characteristic of basic Shock and Vibration 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0t = 0.16 Newmark t = 0.26 DGFEM t = 0.16 DGFEM t = 0.26

Figure 3 :Figure 4 :
Figure 3: Comparison of the stress profiles at  = 0.16,  = 0.20, and  = 0.26 between the proposed DGFEM and the Newmark method.