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

Authors: Pan Guo; Wen-Hua Wu; Jun Zhao
Journal: Shock and Vibration
 (2017)

---

## Abstract

This paper describes a time-discontinuous Galerkin finite element method (DGFEM-βc) 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-βc method is developed to allow the general temperature-displacement vector and its temporal gradient to be discontinuous at a fixed time t. 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-βc 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-βc can effectively capture the discontinuities behaviours of impulsive waves in space in the simulation of high modes and sharp gradients.


---

## Body

## 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 [3–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–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 (βcC) 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-βc 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 layeri are expressed as [15](1a)ρiciτ0i∂2T∂t2+ρici∂T∂t-T0βidiv∂U∂t+τ0i∂2U∂t2-KiΔT=QT(1b)ρi∂2U∂t2-βi∇T-μiΔU-λi+μi∇div⁡U=QU,where T is the temperature; U is the displacement vector; t is time; λi and μi are the Lame constants; ρi is the density; ci is the specific heat at constant strain; βi=αi(3λi+2μi), where αi is the coefficient of linear thermal expansion; Ki is the thermal conductivity; T0 is the reference temperature; QU is the force vector in the displacement field; QT is heat source; and τ0i is the thermal relaxation time. When τ0i=0, (1a)-(1b) reduce to the classical coupled thermoelastic theory. When βi=0, (1a) reduces to the non-Fourier hyperbolic heat conductive equation.The corresponding dimensionless terms are defined as follows:(2)u=UU-,θ=T-T-T-,x∗=xx-,t∗=tt-,fθ∗=QTQ-T,fU∗=QUQ-Uin which(3)U-=K1c1ρ1λ1+2μ1,T-=T0,x-=K1c1ρ1λ1+2μ1,t-=K1c1λ1+2μ1,Q-T=T0ρ1c12λ1+2μ1K1,Q-U=ρ1c1λ1+2μ12K1ρ1λ1+2μ1.Using these nondimensional variables, we restate (1a) and (1b) (with asterisks omitted for convenience) as(4a)RρiRciτi∂2θ∂t2+RρiRci∂θ∂t-aidiv∂u∂t+τi∂2u∂t2-RKjΔθ=fθ(4b)Rρi∂2u∂t2-β1i∇θ-RμjΔu-Rλ+μj∇divu=fu,where τi=τ0ic1(λ1+2μ1)/K1, β1i=βiT0/(λ1+2μ1), ai=βi/ρ1c1, Rρi=ρi/ρ1, Rci=ci/c1, RKj=Kj/K1, Rμj=μj/μ1, Rλ+μj=(λj+μj)/(λ1+μ1), Rλ+2μj=(λj+2μj)/(λ1+2μ1).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 Γu1, respectively, then(5a)θ=θ0tonΓθ1(5b)u=u0tonΓu1,where θ0 and u0 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 Γu2, the following boundary conditions should be satisfied: (5c)Fθ=FbθtonΓθ2(5d)Fu=Fbuton⁡Γu2,where Fbθ(t) and Fbu(t) are the given surface heat flux and the surface traction, respectively.In the present paper, surface impulse heat source is applied with the form(6)Fbθt=q0n,t∈0,t00,t>t0on⁡Γθ2,where q0 and t0 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:(7)MθθMθu0Muuθ¨u¨+CθθCθu00θ˙u˙+Kθθ0KuθKuuθu=FθFu,where(8)Mθθ=∫ΩRρiRciτiNθNθdΩ,Mθu=∫ΩτiaiNθBudΩ,Muu=∫ΩRρiNuNudΩ,Cθu=∫ΩaiNθBudΩ,Cθθ=∫ΩRρiRciNθNθdΩ,Kθθ=∫ΩRKjBθBθdΩ,Kuθ=∫Ωβ1iBuNθdΩ,Kuu=∫ΩRλ+2μjBuBudΩ,Fθ=∫ΩNθfθdΩ,Fu=∫ΩNufudΩin which Nθ and Nu are the shape function for the temperature and displacement, respectively. And Bθ and Bu 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, 0<⋯<tn<⋯<tN.We assume that the unknown variable is not continuous attn and is expressed as(9)wn=wtn+-wtn-,where w(tn+), w(tn-) are discontinuities of variables at tn with w(tn±)=limε→0±⁡(w+ε).For each time subdomainIn=tn+,tn+1-, a time step length is defined Δt=tn+1-tn, for n=0,1,…,N-1. For each time subdomain In, we interpolate the unknown variables, that is, the global nodal temperature-displacements vector {θu}, by using the third-order Hermite time shape functions as(10)θut=θun+N1t+θun+1-N2t+qvn+M1t+qvn+1-M2t,where θun+, θun+1-, qvn+, and qvn+1- represent the global nodal values of the temperature-displacements vector and its time-derivative at times tn, tn+1, respectively. Among them, q=θ˙ and v=u˙. The global nodal values of temperature-displacements vector and its time-derivative, that is, θun-,qvn-, at time tn- can be obtained at the end of the previous time step. In addition, (11)N1=N1t=λ12λ1+3λ2,N2=N2t=λ22λ2+3λ1M1=M1t=λ12λ2Δt,M2=M2t=-λ1λ22Δt.λ1=tn+1-tΔt,λ2=t-tnΔtThe global nodal values{qv} (the temporal derivative of the temperature-displacements vector) at arbitrary time t∈tn,tn+1 are treated as an independent variable by linear time shape functions as(12)qvt=qvn+λ1t+qvn+1-λ2t.We defineU~={θu}, V={qv} as the basic independent variables; U~˙ and U~¨ are then the first-order and second-order temporal derivative, respectively, of U~.By dropping the superscripts of the vectorsU~n+, U~n+1-, Vn+, Vn+1- and the time variable t, we can rewrite (10) and (12) as(13)U~=U~nN1+U~n+1N2+VnM1+Vn+1M2Vt=Vnλ1t+Vn+1λ2t.The generalised thermoelastic weak forms of the semidiscretised equation (9) and the constraint condition U~˙-V=0, along with the discontinuity conditions U~ and V on a typical time subdomain can be expressed as(14)∫InδVTMV˙+CV+KU~-Fdt+∫InδU~TKU~˙-Vdt+δU~nTKU~n+δVnTMVn=0.Substituting (9) and (13) into (14), we obtain the following matrix equation of DGFEM for solving from independent variations of δU~n+, δU~n+1-, δVn+, and δVn+1-(15)U~n=U~n-(16)M+Δt6C-Δt212K-Δt6C-Δt212KΔt2C+Δt23KM+Δt2C+Δt26KVnVn+1=F1e-F2e+MVn-F1e+F2e+MVn--ΔtKU~n-(17)U~n+1=U~n-+12ΔtVn+Vn+1,where F1e=Δt/3Fne+Δt/6Fn+1e, F2e=Δt/6Fne+Δt/3Fn+1e(18)M=MθθMθu0Muu,C=CθθCθu00,K=Kθθ0KuθKuu,F=FθFu.This equation means that the continuity of nodal temperature-displacement vectorU at any time level is automatically ensured. The discontinuity is located at nodal derivative vector V 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:(19)Cae=βciKe(20)βci=2ξiωn,where Cae is the stiffness proportional matrix, βci is the damping constants, ξi is the damping ratio estimated from the global behaviour of the system, and ωn 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 βci should be less than or equal to the time step.Substituting the stiffness proportional artificial damping matrix equation (19) into (16), we obtain(21)M+Δt6Cw-Δt212K-Δt6Cw-Δt212KΔt2Cw+Δt23KM+Δt2Cw+Δt26KVnVn+1=F1e-F2e+MVn-F1e+F2e+MVn--ΔtKU~n-,where Cw=Ca+C is the global damping term.Equations (17), (23), and (19) form the final expressions for the present DGFEM.


## 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:(7)MθθMθu0Muuθ¨u¨+CθθCθu00θ˙u˙+Kθθ0KuθKuuθu=FθFu,where(8)Mθθ=∫ΩRρiRciτiNθNθdΩ,Mθu=∫ΩτiaiNθBudΩ,Muu=∫ΩRρiNuNudΩ,Cθu=∫ΩaiNθBudΩ,Cθθ=∫ΩRρiRciNθNθdΩ,Kθθ=∫ΩRKjBθBθdΩ,Kuθ=∫Ωβ1iBuNθdΩ,Kuu=∫ΩRλ+2μjBuBudΩ,Fθ=∫ΩNθfθdΩ,Fu=∫ΩNufudΩin which Nθ and Nu are the shape function for the temperature and displacement, respectively. And Bθ and Bu 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, 0<⋯<tn<⋯<tN.We assume that the unknown variable is not continuous attn and is expressed as(9)wn=wtn+-wtn-,where w(tn+), w(tn-) are discontinuities of variables at tn with w(tn±)=limε→0±⁡(w+ε).For each time subdomainIn=tn+,tn+1-, a time step length is defined Δt=tn+1-tn, for n=0,1,…,N-1. For each time subdomain In, we interpolate the unknown variables, that is, the global nodal temperature-displacements vector {θu}, by using the third-order Hermite time shape functions as(10)θut=θun+N1t+θun+1-N2t+qvn+M1t+qvn+1-M2t,where θun+, θun+1-, qvn+, and qvn+1- represent the global nodal values of the temperature-displacements vector and its time-derivative at times tn, tn+1, respectively. Among them, q=θ˙ and v=u˙. The global nodal values of temperature-displacements vector and its time-derivative, that is, θun-,qvn-, at time tn- can be obtained at the end of the previous time step. In addition, (11)N1=N1t=λ12λ1+3λ2,N2=N2t=λ22λ2+3λ1M1=M1t=λ12λ2Δt,M2=M2t=-λ1λ22Δt.λ1=tn+1-tΔt,λ2=t-tnΔtThe global nodal values{qv} (the temporal derivative of the temperature-displacements vector) at arbitrary time t∈tn,tn+1 are treated as an independent variable by linear time shape functions as(12)qvt=qvn+λ1t+qvn+1-λ2t.We defineU~={θu}, V={qv} as the basic independent variables; U~˙ and U~¨ are then the first-order and second-order temporal derivative, respectively, of U~.By dropping the superscripts of the vectorsU~n+, U~n+1-, Vn+, Vn+1- and the time variable t, we can rewrite (10) and (12) as(13)U~=U~nN1+U~n+1N2+VnM1+Vn+1M2Vt=Vnλ1t+Vn+1λ2t.The generalised thermoelastic weak forms of the semidiscretised equation (9) and the constraint condition U~˙-V=0, along with the discontinuity conditions U~ and V on a typical time subdomain can be expressed as(14)∫InδVTMV˙+CV+KU~-Fdt+∫InδU~TKU~˙-Vdt+δU~nTKU~n+δVnTMVn=0.Substituting (9) and (13) into (14), we obtain the following matrix equation of DGFEM for solving from independent variations of δU~n+, δU~n+1-, δVn+, and δVn+1-(15)U~n=U~n-(16)M+Δt6C-Δt212K-Δt6C-Δt212KΔt2C+Δt23KM+Δt2C+Δt26KVnVn+1=F1e-F2e+MVn-F1e+F2e+MVn--ΔtKU~n-(17)U~n+1=U~n-+12ΔtVn+Vn+1,where F1e=Δt/3Fne+Δt/6Fn+1e, F2e=Δt/6Fne+Δt/3Fn+1e(18)M=MθθMθu0Muu,C=CθθCθu00,K=Kθθ0KuθKuu,F=FθFu.This equation means that the continuity of nodal temperature-displacement vectorU at any time level is automatically ensured. The discontinuity is located at nodal derivative vector V 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:(19)Cae=βciKe(20)βci=2ξiωn,where Cae is the stiffness proportional matrix, βci is the damping constants, ξi is the damping ratio estimated from the global behaviour of the system, and ωn 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 βci should be less than or equal to the time step.Substituting the stiffness proportional artificial damping matrix equation (19) into (16), we obtain(21)M+Δt6Cw-Δt212K-Δt6Cw-Δt212KΔt2Cw+Δt23KM+Δt2Cw+Δt26KVnVn+1=F1e-F2e+MVn-F1e+F2e+MVn--ΔtKU~n-,where Cw=Ca+C 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 Figure1, defined on the space-time domain 0,0.25×0,0.26. The associated dimensionless material constants are taken as (22)β1i=5.035×10-3,i=1,2,3Rρi=5,1,2,i=1,2,3ai=0.926,i=1,2,3.Rci=1.0,0.2,1.0,i=1,2,3RKi=1.0,0.2,1.0,i=1,2,3Rλ+2μi=2.0,1.0,5.0,i=1,2,3The dimensionless thermal relaxation time is chosen as τi=0.1,i=1,2,3 [15].Figure 1
A uniform thermoelastic configuration with one end fixed and the other end subjected to a rectangular impulse heat load.A time step sizeΔt=1.0×10-3 and the artificial damping constants βci=5.0×10-5,i=1,2,3 are used. Three layers are discretised into forty-five, ten, and forty-five elements along thickness, respectively.Initial conditions are(23)Tx,0=0,T˙x,0=0U~x,0=0,U~˙x,0=0.Boundary conditions are(24)TL,t=0,U~L,t=0Fθ0,t=100t∈0,5×10-30t>5×10-3.To compare the reliability of the proposed DGFEM and traditional method, we evaluate the L-S model ofai=0.926,i=1,2,3 in this case. Figures 2 and 3 represent the temperature and stress distributions along the distance of x-axis at different nondimensional time t=0.16, t=0.20, t=0.22, and t=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 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.Figure 2
Comparison of the temperature profiles att=0.16, t=0.20, t=0.22, and t=0.26 between the proposed DGFEM and the Newmark method.Figure 3
Comparison of the stress profiles att=0.16, t=0.20, and t=0.26 between the proposed DGFEM and the Newmark method.


### 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 with31×31 four-node finite element elements under a plane strain assumption, as depicted in Figure 4.Figure 4
Mesh, geometry, and boundary conditions of the two-dimensional model.The dimensionless material property data are defined as follows:β11=5.035×10-3, a1=0.926, τ1=0.1 [15], β12=5.035×10-3, a2=0.926, τ2=0.2, β13=5.035×10-3, a3=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 Δt=3.0×10-3, and an artificial damping constant βci=1.5×10-3 is used.Figures5-6 show snapshots of the temperature and the stress distributions versus time for t=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-βc 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.Figure 5
Comparison of the temperature distributions between the Newmark method and the DGFEM.
(a)
Newmark method att=0.06, t=0.12, t=0.18, and t=0.24
(b)
DGFEM att=0.06, t=0.12, t=0.18, and t=0.24Figure 6
Comparison of displacement profiles between Newmark method and DGFEM.
(a)
Newmark method att=0.06, t=0.12, t=0.18, and t=0.24
(b)
DGFEM att=0.06, t=0.12, t=0.18, and t=0.24


## 4.1. 1D Thermoelastic Wave Propagation Problem

The first example is a 1D thermoelastic wave propagation problem, as shown in Figure1, defined on the space-time domain 0,0.25×0,0.26. The associated dimensionless material constants are taken as (22)β1i=5.035×10-3,i=1,2,3Rρi=5,1,2,i=1,2,3ai=0.926,i=1,2,3.Rci=1.0,0.2,1.0,i=1,2,3RKi=1.0,0.2,1.0,i=1,2,3Rλ+2μi=2.0,1.0,5.0,i=1,2,3The dimensionless thermal relaxation time is chosen as τi=0.1,i=1,2,3 [15].Figure 1
A uniform thermoelastic configuration with one end fixed and the other end subjected to a rectangular impulse heat load.A time step sizeΔt=1.0×10-3 and the artificial damping constants βci=5.0×10-5,i=1,2,3 are used. Three layers are discretised into forty-five, ten, and forty-five elements along thickness, respectively.Initial conditions are(23)Tx,0=0,T˙x,0=0U~x,0=0,U~˙x,0=0.Boundary conditions are(24)TL,t=0,U~L,t=0Fθ0,t=100t∈0,5×10-30t>5×10-3.To compare the reliability of the proposed DGFEM and traditional method, we evaluate the L-S model ofai=0.926,i=1,2,3 in this case. Figures 2 and 3 represent the temperature and stress distributions along the distance of x-axis at different nondimensional time t=0.16, t=0.20, t=0.22, and t=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 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.Figure 2
Comparison of the temperature profiles att=0.16, t=0.20, t=0.22, and t=0.26 between the proposed DGFEM and the Newmark method.Figure 3
Comparison of the stress profiles att=0.16, t=0.20, and t=0.26 between the proposed DGFEM and the Newmark method.


## 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 with31×31 four-node finite element elements under a plane strain assumption, as depicted in Figure 4.Figure 4
Mesh, geometry, and boundary conditions of the two-dimensional model.The dimensionless material property data are defined as follows:β11=5.035×10-3, a1=0.926, τ1=0.1 [15], β12=5.035×10-3, a2=0.926, τ2=0.2, β13=5.035×10-3, a3=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 Δt=3.0×10-3, and an artificial damping constant βci=1.5×10-3 is used.Figures5-6 show snapshots of the temperature and the stress distributions versus time for t=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-βc 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.Figure 5
Comparison of the temperature distributions between the Newmark method and the DGFEM.
(a)
Newmark method att=0.06, t=0.12, t=0.18, and t=0.24
(b)
DGFEM att=0.06, t=0.12, t=0.18, and t=0.24Figure 6
Comparison of displacement profiles between Newmark method and DGFEM.
(a)
Newmark method att=0.06, t=0.12, t=0.18, and t=0.24
(b)
DGFEM att=0.06, t=0.12, t=0.18, and t=0.24


## 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-βc 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.



---

*Source: 8398673-2017-09-11.xml*