Printer Friendly

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

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 ([[beta].sub.c]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-[[beta].sub.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 layer i are expressed as [15]

[mathematical expression not reproducible] (1a)

[mathematical expression not reproducible] (1b)

where T is the temperature; U is the displacement vector; t is time; [[lambda].sup.(i)] and [[mu].sup.(i)] are the Lame constants; [[rho].sup.(i)] is the density; [c.sup.(i)] is the specific heat at constant strain; [[beta].sup.(i)] = [[alpha].sup.(i)](3[[lambda].sup.(i)] + 2[[mu].sup.(i)]),where [[alpha].sup.(i)] is the coefficient of linear thermal expansion; [K.sup.(i)] is the thermal conductivity; [T.sub.0] is the reference temperature; [Q.sup.u] is the force vector in the displacement field; [Q.sup.t] is heat source; and [t.sup.(i).sub.0] is the thermal relaxation time. When [t.sup.(i).sub.0] = 0, (1a)-(1b) reduce to the classical coupled thermoelastic theory. When [[beta].sup.(i)] = 0, (1a) reduces to the non-Fourier hyperbolic heat conductive equation.

The corresponding dimensionless terms are defined as follows:

[mathematical expression not reproducible] (2)

in which

[mathematical expression not reproducible]. (3)

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

[mathematical expression not reproducible] (4a)

[mathematical expression not reproducible] (4b)

where [mathematical expression not reproducible].

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 [[GAMMA].sub.[theta]1] and [[GAMMA].sub.u1], respectively, then

[theta] = [[theta].sub.0](t) on [[GAMMA].sub.[theta]1] (5a)

u = [u.sub.0](t) on [[GAMMA].sub.u1], (5b)

where [[theta].sub.0] and [u.sub.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 [[GAMMA].sub.[theta]2] and [[GAMMA].sub.u2], the following boundary conditions should be satisfied:

[F.sup.[theta]] = [F.sup.[theta].sub.b](t) on [[GAMMA].sub.[theta]2] (5c)

[F.sup.u] = [F.sup.u.sub.b](t) on [[GAMMA].sub.u2], (5d)

where [F.sup.[theta].sub.b](t) and [F.sup.u.sub.b](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

[mathematical expression not reproducible], (6)

where [q.sub.0] and [t.sub.0] 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:

[mathematical expression not reproducible], (7)


[mathematical expression not reproducible] (8)

in which [N.sub.[theta]] and [N.sub.u] are the shape function for the temperature and displacement, respectively. And [B.sub.[theta]] and [B.sub.u] 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 < ... < [t.sub.n] < ... < [t.sub.N].

We assume that the unknown variable is not continuous at [t.sub.n] and is expressed as

[[[w.sub.n]]] = w([t.sup.+.sub.n]) - w([t.sup.-.sub.n]), (9)

where w([t.sup.+.sub.n]), w([t.sup.-.sub.n]) are discontinuities of variables at [t.sub.n] with [mathematical expression not reproducible].

For each time subdomain [I.sub.n] = ([t.sup.+.sub.n][t.sup.-.sub.n+1]), a time step length is defined [DELTA]t = [t.sup.n.sub.+1] - [t.sub.n], for n = 0,1, ..., N-1. For each time subdomain [I.sub.n], we interpolate the unknown variables, that is, the global nodal temperature-displacements vector [mathematical expression not reproducible], by using the third-order Hermite time shape functions as

[mathematical expression not reproducible], (10)

where [mathematical expression not reproducible],represent the global nodal values of the temperature-displacements vector and its time-derivative at times [t.sub.n], [t.sub.n+1], respectively. Among them, [mathematical expression not reproducible]. The global nodal values of temperature-displacements vector and its time-derivative, that is, [mathematical expression not reproducible], at time [t.sup.-.sub.n] can be obtained at the end of the previous time step. In addition,

[mathematical expression not reproducible] (11)

The global nodal values [mathematical expression not reproducible] (the temporal derivative of the temperature-displacements vector) at arbitrary time i [member of] [[t.sub.n], [t.sub.n+1]] are treated as an independent variable by linear time shape functions as

[mathematical expression not reproducible]. (12)

We define [mathematical expression not reproducible] as the basic independent variables; [mathematical expression not reproducible] are then the first-order and second-order temporal derivative, respectively, of [mathematical expression not reproducible].

By dropping the superscripts of the vectors [mathematical expression not reproducible], [V.sup.+.sub.n], [V.sup.-.sub.n+1] and the time variable t, we can rewrite (10) and (12) as

[mathematical expression not reproducible]. (13)

The generalised thermoelastic weak forms of the semidiscretised equation (9) and the constraint condition [mathematical expression not reproducible] - V = 0, along with the discontinuity conditions [mathematical expression not reproducible] and V on a typical time subdomain can be expressed as

[mathematical expression not reproducible]. (14)

Substituting (9) and (13) into (14), we obtain the following matrix equation of DGFEM for solving from independent variations of [delta][[??].sup.+.sub.n], [delta][[??].sup.-.sub.n+1], [delta][V.sup.+.sub.n], and [delta][V.sup.-.sub.n+1]

[[??].sub.n] = [[??].sup.-.sub.n] (15)

[mathematical expression not reproducible] (16)

[mathematical expression not reproducible], (17)

where [F.sup.e.sub.1] = ([DELTA]t/3)[F.sup.e.sub.n] + ([DELTA]t/6)[F.sup.e.sub.n+i], [F.sup.e.sub.2] = ([DELTA]t/6)[F.sup.e.sub.n] + ([DELTA]t/3)[F.sup.e.sub.n+i]

[mathematical expression not reproducible]. (18)

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.

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:

[[[C.sub.a]].sup.e] = [[beta].sup.(i).sub.c][[K].sup.e] (19)

[[beta].sup.(i).sub.c] = 2[[xi].sup.(i)]/[[omega].sub.n], (20)

where [[[C.sub.a]].sup.e] is the stiffness proportional matrix, [[beta].sup.(i).sub.c] is the damping constants, [[xi].sup.(i)] is the damping ratio estimated from the global behaviour of the system, and [[omega].sub.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 [[beta].sup.(i).sub.c] should be less than or equal to the time step.

Substituting the stiffness proportional artificial damping matrix equation (19) into (16), we obtain

[mathematical expression not reproducible], (21)

where [C.sub.w] = [C.sub.a] + 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 Figure 1, defined on the space-time domain [0,0.25] x [0,0.26]. The associated dimensionless material constants are taken as

[mathematical expression not reproducible] (22)

The dimensionless thermal relaxation time is chosen as [[tau].sup.(i)] = 0.1, i = 1,2,3 [15].

A time step size [DELTA]t = 1.0 x [10.sup.-3] and the artificial damping constants [[beta].sup.(i).sub.c] = 5.0 x [10.sup.-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

[mathematical expression not reproducible]. (23)

Boundary conditions are

[mathematical expression not reproducible]. (24)

To compare the reliability of the proposed DGFEM and traditional method, we evaluate the L-S model of [a.sup.(i)] = 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.

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 31 x 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: [[beta].sup.(1).sub.1] = 5.035x [10.sup.-3], [a.sup.(1)] = 0.926, [[tau].sup.(1)] = 0.1 [15], [[beta].sup.(2).sub.1] = 5.035x [10.sup.-3], [a.sup.(2)] = 0.926, [[tau].sup.(2)] = 0.2, [[beta].sup.(3).sub.1] = 5.035x[10.sup.-3], [a.sup.(3)] = 0.926, [[tau].sup.(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 [DELTA]t = 3.0 x [10.sup.-3], and an artificial damping constant [[beta].sup.(i).sub.c] = 1.5 x [10.sup.-3] is used.

Figures 5-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-[[beta].sub.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.

5. Conclusions

Multilayer materials are widely used in the development of microelectronics systems, photoelectric equipment, micro-sensors, 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-[[beta].sub.c] was demonstrated through comparison with traditional time-stepping algorithms such as the Newmark [gamma] 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).


[1] Y.-M. Yang, X.-Z. Wang, andW.-J. Zhang, "Thermoelastic stress analysis of multilayered films in a micro-thermoelectric cooling device," Acta Mechanica Sinica, vol. 28, no. 6, pp. 1644-1650, 2012.

[2] A. H. Akbarzadeh, M. H. Babaei, and Z. T. Chen, "Thermopiezoelectric analysis of a functionally graded piezoelectric medium," International Journal of Applied Mechanics, vol. 3, no. 1, pp. 47-68, 2011.

[3] J. W. Fu, A. H. Akbarzadeh, Z. T. Chen, L. F. Qian, and D. Pasini, "Non-Fourier heat conduction in a sandwich panel with a cracked foam core," International Journal of Thermal Sciences, vol. 102, pp. 263-273, 2016.

[4] J. W. Fu, Z. T. Chen, and L. F. Qian, "Coupled thermoelastic analysis of a multi-layered hollow cylinder based on the C-T theory and its application on functionally graded materials," Composite Structures, vol. 131, pp. 139-150, 2015.

[5] A. Singh and E. B. Tadmor, "Thermal parameter identification for non-Fourier heat transfer from molecular dynamics," Journal of Computational Physics, vol. 299, pp. 667-686, 2015.

[6] B. S. Yilbas and A. Y. Al-Dweik, "Closed form solutions for thermal stress field due to non-equilibrium heating during laser short-pulse irradiation," Physica B: Condensed Matter, vol. 407, no. 12, pp. 2169-2175, 2012.

[7] C. Catteneo, "A form of heat conduction equation which eliminates the paradox of instantaneous propagation," Computes Rendus, vol. 247, pp. 431-433, 1958.

[8] P. Vernotte, "Les paradoxes de la theorie continue de lequation de la chaleur," Computes Rendus, vol. 246, pp. 3154-3155, 1958.

[9] H. W. Lord and Y. Shulman, "A generalized dynamical theory of thermoelasticity," Journal of the Mechanics and Physics of Solids, vol. 15, no. 5, pp. 299-309, 1967.

[10] J. H. Prevost and D. Tao, "Finite element analysis of dynamic coupled thermoelasticity problems with relaxation times," Journal of Applied Mechanics, vol. 50, no. 4 a, pp. 817-822, 1983.

[11] K. K. Tamma and R. R. Namburu, "Computational approaches with applications to non-classical and classical thermomechanical problems," Applied Mechanics Reviews, vol. 50, no. 9, pp. 514-551, 1997.

[12] K. K. Tamma, X. Zhou, and R. Valasutean, "Computational algorithms for transient analysis: the burden of weight and consequences towards formalizing discrete numerically assigned [DNA] algorithmic markers: Wp-family," Computer Methods in Applied Mechanics and Engineering, vol. 149, no. 1-4, pp. 153-188, 1997.

[13] X. Li, D. Yao, and R. W. Lewis, "A discontinuous Galerkin finite element method for dynamic and wave propagation problems in non-linear solids and saturated porous media," International Journal for Numerical Methods in Engineering, vol. 57, no. 12, pp. 1775-1800, 2003.

[14] W. Wu and X. Li, "Application of the time discontinuous Galerkin finite element method to heat wave simulation," International Journal of Heat and Mass Transfer, vol. 49, no. 9-10, pp. 1679-1684, 2006.

[15] P. Guo, W.-H. Wu, and Z.-G. Wu, "A time discontinuous Galerkin finite element method for generalized thermo-elastic wave analysis, considering non-Fourier effects," Acta Mechanica, vol. 225, no. 1, pp. 299-307, 2014.

Pan Guo, (1) Wen-Hua Wu, (2) and Jun Zhao (3)

(1) National Center for International Joint Research of Micro-Nano Moulding Technology, School of Mechanics and Engineering Science of Zhengzhou University, Zhengzhou 450001, China

(2) The 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

(3) School of Mechanics and Engineering Science, Zhengzhou University, Zhengzhou 450001, China

Correspondence should be addressed to Wen-Hua Wu;

Received 2 May 2017; Accepted 6 July 2017; Published 11 September 2017

Academic Editor: Tai Thai

Caption: FIGURE 1: A uniform thermoelastic configuration with one end fixed and the other end subjected to a rectangular impulse heat load.

Caption: FIGURE 2: Comparison of the temperature profiles at t = 0.16, t = 0.20, t = 0.22, and t = 0.26 between the proposed DGFEM and the Newmark method.

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

Caption: FIGURE 4: Mesh, geometry, and boundary conditions of the two-dimensional model.

Caption: FIGURE 5: Comparison of the temperature distributions between the Newmark method and the DGFEM.

Caption: FIGURE 6: Comparison of displacement profiles between Newmark method and DGFEM.
COPYRIGHT 2017 Hindawi Limited
No portion of this article can be reproduced without the express written permission from the copyright holder.
Copyright 2017 Gale, Cengage Learning. All rights reserved.

Article Details
Printer friendly Cite/link Email Feedback
Title Annotation:Research Article
Author:Guo, Pan; Wu, Wen-Hua; Zhao, Jun
Publication:Shock and Vibration
Article Type:Report
Date:Jan 1, 2017
Previous Article:Dynamic Optimization Design of Cranes Based on Human-Crane-Rail System Dynamics and Annoyance Rate.
Next Article:Numerical Study on the In-Plane and Out-of-Plane Resistance of Brick Masonry Infill Panels in Steel Frames.

Terms of use | Privacy policy | Copyright © 2019 Farlex, Inc. | Feedback | For webmasters