This paper describes a time-discontinuous Galerkin finite element method (DGFEM-) for the generalised thermoelastic problem of multilayer materials subjected to a transient high-frequency heat source. The governing and constitutive relations are presented on the basis of the well-known Lord–Shulman (L–S) theory. A DGFEM- method is developed to allow the general temperature-displacement vector and its temporal gradient to be discontinuous at a fixed time . A stiffness proportional artificial damping term is added to the final DG discretisation form to filter out the spurious numerical oscillations in the wave-after stage and at adjacent-layer interfaces. The numerical results show that the present DGFEM- provides much more accurate solutions for generalised thermoelastic coupled behaviour of multilayer structures. Compared with widely used traditional numerical methods (e.g., the Newmark method), the present DGFEM- can effectively capture the discontinuities behaviours of impulsive waves in space in the simulation of high modes and sharp gradients.

1. 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 [36]. 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) [1012], 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 () 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.

2. 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]where is the temperature; is the displacement vector; is time; and are the Lame constants; is the density; is the specific heat at constant strain; , where is the coefficient of linear thermal expansion; is the thermal conductivity; is the reference temperature; is the force vector in the displacement field; is heat source; and is the thermal relaxation time. When , (1a)-(1b) reduce to the classical coupled thermoelastic theory. When , (1a) reduces to the non-Fourier hyperbolic heat conductive equation.

The corresponding dimensionless terms are defined as follows:in which

Using these nondimensional variables, we restate (1a) and (1b) (with asterisks omitted for convenience) aswhere , , , , , , , , .

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 and , respectively, thenwhere and 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 and , the following boundary conditions should be satisfied: where and are the given surface heat flux and the surface traction, respectively.

In the present paper, surface impulse heat source is applied with the formwhere and are the amplitude of thermal shock and the constant time.

3. 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: wherein 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.

3.2. 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, .

We assume that the unknown variable is not continuous at and is expressed aswhere , are discontinuities of variables at with .

For each time subdomain , a time step length is defined , for . For each time subdomain , we interpolate the unknown variables, that is, the global nodal temperature-displacements vector , by using the third-order Hermite time shape functions aswhere , , , and represent the global nodal values of the temperature-displacements vector and its time-derivative at times , , respectively. Among them, and . The global nodal values of temperature-displacements vector and its time-derivative, that is, , at time can be obtained at the end of the previous time step. In addition,

The global nodal values (the temporal derivative of the temperature-displacements vector) at arbitrary time are treated as an independent variable by linear time shape functions as

We define , as the basic independent variables; and are then the first-order and second-order temporal derivative, respectively, of .

By dropping the superscripts of the vectors , , , 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 , along with the discontinuity conditions and 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 , , , and where ,

This equation means that the continuity of nodal temperature-displacement vector at any time level is automatically ensured. The discontinuity is located at nodal derivative vector at discretised time levels.

3.3. 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 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 obtainwhere is the global damping term.

Equations (17), (23), and (19) form the final expressions for the present DGFEM.

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

4.1. 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 . The associated dimensionless material constants are taken as The dimensionless thermal relaxation time is chosen as [15].

A time step size and the artificial damping constants are used. Three layers are discretised into forty-five, ten, and forty-five elements along thickness, respectively.

Initial conditions are

Boundary conditions are

To compare the reliability of the proposed DGFEM and traditional method, we evaluate the L-S model of in this case. Figures 2 and 3 represent the temperature and stress distributions along the distance of -axis at different nondimensional time , , , and 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 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.

4.2. 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 four-node finite element elements under a plane strain assumption, as depicted in Figure 4.

The dimensionless material property data are defined as follows: , , [15], , , , , , . 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 , and an artificial damping constant is used.

Figures 5-6 show snapshots of the temperature and the stress distributions versus time for , 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.

5. 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 third-order Hermite time shape functions and linear time shape functions were applied to interpolate the temperature-displacement 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 wave-front 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.

Conflicts of Interest

The authors declare that they have no conflicts of interest.


The authors are pleased to acknowledge the support of this work by the National Natural Science Foundation of China through Contract Grant no. 15572072, National Key Basic Research and Development Program through Contract Grant nos. 2014CB046803 and 2016ZX05028-002-005, and Startup Research Fund of Zhengzhou University (1511327001).