Front. Mater. Frontiers in Materials Front. Mater. 2296-8016 Frontiers Media S.A. 963149 10.3389/fmats.2022.963149 Materials Original Research Forced vibration analysis of inhomogeneous quasicrystal coating in a thermal environment Feng et al. 10.3389/fmats.2022.963149 Feng Xin 1 2 Zhang Liangliang 1 Zhu Zhaowei 1 Li Yang 3 Gao Yang 1 * 1 College of Science, China Agricultural University, Beijing, China 2 College of Engineering, China Agricultural University, Beijing, China 3 Department of Aerospace Science and Technology, Space Engineering University, Beijing, China

Edited by: Chunli Zhang, Zhejiang University, China

Reviewed by: Xiangyu Li, Southwest Jiaotong University, China

Junhong Guo, Inner Mongolia University of Technology, China

 *Correspondence: Yang Gao, gaoyangg@gmail.com

This article was submitted to Mechanics of Materials, a section of the journal Frontiers in Materials

 26 09 2022 2022 9 963149 07 06 2022 18 07 2022 Copyright © 2022 Feng, Zhang, Zhu, Li and Gao. 2022 Feng, Zhang, Zhu, Li and Gao

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

The tremendous attention of researchers has been attracted to the unusual properties of quasicrystals in coatings. In this paper, the exact solutions of the functionally graded multilayered two-dimensional quasicrystal coating structures in a thermal environment are derived for advanced boundary-value problems with mixed boundary conditions. The state space method is formulated to the thermal coupling with quasicrystal linear elastic theory that derives the state equations for functionally graded quasicrystal coating structures along the thickness direction. The mixed supported boundary conditions in the x-direction and the simply supported boundary conditions in the y-direction are subjected to time-harmonic temperature loadings, which are represented by means of the differential quadrature technique and Fourier series expansions, respectively. Traction on both the bottom and top surfaces is free, and perfect thermal and mechanical contacts between constituents are incorporated at the internal interfaces. A global propagator matrix, which connects the field variables at the top interface to those at the bottom interface for the whole coating structure, is further completed by joint coupling matrices to overcome the numerical instabilities. Finally, three application examples are proposed to throw light on various effects of the power law index, frequency, and different boundary conditions on the field variables in three-layer coating structures. The present solution can serve as a benchmark for the modeling of functionally graded quasicrystal coating structures based on various numerical methods.

 forced vibration analysis quasicrystal coating structures functionally graded materials mixed boundary conditions thermoelastic analysis 11972365 12102458 12102481 1101-240001 National Natural Science Foundation of China10.13039/501100001809 China Agricultural University Education Foundation10.13039/501100012454 Introduction

Multilayered quasicrystal (QC) composites possess huge potential applications in the aeronautical and aerospace industries, manufacturing of semiconductor and electrical devices as well as some modern high-technology sectors due to their unique quasi-periodic structure and excellent physical properties (Fan, 2011; Guo et al., 2020). In addition, due to the high hardness, good chemical stability, high temperature resistance, and corrosion resistance of QC materials, QC coating technology has been tremendously developed (Beardsley, 2008; Mora et al., 2020), which proposed new ways to improve the aeroengine performance. Compared with monolithic samples, artificially designed layered QC coating structures can provide better combinations of strength and ductility, strength-to-weight ratios, and stiffness-to-weight ratios. Meanwhile, the temperature deformation and thermal stress of the plates and coatings are also important factors that must be considered in the engineering structure design. Thus, the responses of the QC composite material structure in the thermal environment have always been a research hotspot, including static response [residual stress caused by high temperature (Polishchuk et al., 2012)], steady state response [stress problem in the heat insulation process (Li and Liu, 2018)], and dynamic response (instantaneous thermal stress in the impact state).

In addition, QC materials as multilayered structures have some weaknesses, such as they may crack at low temperature or peel off at high temperature, which will reduce the life-time and reliability of the devices. For instance, QC composite thermoelectric converters are exposed to thermal environments, including severe temperature gradients, which can generate interlayer thermal stresses that affect the structure in practice (Maciá, 2004). Thus, in order to improve thermal and mechanical material properties, various kinds of models have been proposed to overcome these cark and problems. Li (2013) and Li et al. (2015) obtained the fundamental solutions for the one-dimensional (1D) QC medium subjected to two identical thermal loadings on the upper and lower crack lips. Under surface heat loadings, the interlayer stress for 1D hexagonal piezoelectric QC simply supported nanoplate (Zhang et al., 2019) decreases continuously with the increase of the functional gradient factors. Based on the QC nonlocal theory, Yang et al. (2018) obtained the thermo-elastic analysis solution for the two-dimensional (2D) QC simply supported nanoplates subjected to a temperature change. According to the pseudo-Stroh formalism, the deformation of a conductive elliptic hole embedded in the 2D decagonal QC plane (Guo et al., 2016) was derived under a remotely uniform heat flow. In addition, the case of the thermal conductivity measurement along with periodic and quasi-periodic directions in decagonal QCs (Matsukawa et al., 1999; Dolinsek, 2012) and the thermo-elastic response of 2D QC plates (Yang et al., 2014; Huang et al., 2021) were investigated.

For many years there has been considerable interest in the static and dynamic response in multilayered crystal and QC media due to its many applications in composite materials, nondestructive evaluation, geophysics, surface acoustic wave devices, etc. (Bak, 1985a). The propagator matrix method is one of the most effective methods to analyze these media (Huang et al., 2019). The exact solutions of vibration frequency (Sun et al., 2021) and buckling load (Zhang et al., 2022) for 2D decagonal QC multilayered nanoplates can be derived by using the propagator matrix method, which connects the field variables at the upper and lower interfaces of each layer. However, propagator matrix method implementation may suffer from numerical instabilities (Chen et al., 2003; Lü et al., 2006) when the frequency is high, the thickness is large, and the discrete point number exceeds a certain value for the differential quadrature methods (DQM). Various techniques and methods have been proposed to overcome this problem of the propagator matrix method. Dunkin (1965) derived the delta operator technique to alleviate the problem for isotropic materials. Although this method has been developed for the anisotropic material (Castaings and Hosten, 1994), it requires the computation of large-order delta matrices. In addition, Mal (1988) presented the direct global matrix method by reducing the transfer distance. But this method occupied substantial memory storage and computation time. After that, an improved recursive algorithm of the stiffness matrix method (Tan, 2006; Tan, 2011) is presented for computing both total and surface stiffness matrices for a general multilayered anisotropic media. Different from the stiffness matrix method, the traditional propagation matrix method was firstly used to continuously connect the interface variables for multilayered plates with perfect/imperfect interfaces, which is further completed by the dual variable and position technique (Li et al., 2021; Vattré et al., 2021; Vattré and Pan, 2021) to improve the stability of frequency calculation. According to the continuity conditions, one stable approach that involves joint coupling matrices (Lü et al., 2008) is used to deal with the perfect and imperfect interfaces.

As reviewed above, previous analyses for the thermoelastic QC multilayered structure have been conducted under the assumption of static deformation, within which the exact solutions have been evaluated for models with simply supported boundary conditions only. In this paper, forced vibration analysis is carried out on the thermoelastic, anisotropic, multilayered, and inhomogeneous QC coating structures with mixed supported boundary conditions. DQM is utilized to simulate the simply and/or clamped supported boundary conditions along the x-direction, and the Fourier series expansion as the exact closed-form solutions are used to satisfy the boundary conditions in the y-direction. In addition, the basic equations of QC material are converted into the state equations composed of partial differential by utilizing the state-space method. A global propagator matrix is derived by joint coupling matrices to overcome the numerical instabilities. Finally, the typical numerical examples are presented to verify the effectiveness of this method, and the results are very useful for the design and understanding of the characterization of the functionally graded (FG) QC materials in their applications to multilayered systems.

 Description of thermoelastic problems of quasicrystal structures

Figure 1A presents the three-dimensionless multilayered QC thermal barrier coating structure model that is composed of M-bonded orthotropic, linearly thermoelastic, and FG layers. The relationship between the global Cartesian coordinate system and the local material coordinate system of this model is assumed to be (x, y, z) = (x 1, x 2, x 3). The arbitrary, finite, and inhomogeneous pth layer (p = 1, 2, 3, …, M) is surrounded by z p and z p-1 interfaces, respectively, and the in-plane x - and y -directions are aligned with the horizontal edges of the mode boundaries. In addition, the corresponding thickness is h p = z p − z p−1. It also follows that the bottom and top surfaces of the laminate are z 0 = 0 and z M = H, respectively.

(A) An M-layered FG 2D QC structure. (B) Material property variation in terms of volume fraction of exterior material.

Figure 1B describes the variation of the material properties in terms of the volume fraction of external material along the z-direction for pth layer. The pth layer FG material properties vary along the thickness direction, and they can be a combination of 2D QC materials. According to the multilayered structure theory, the pth layer can be divided into arbitrary mathematical layers along the thickness direction. For the material property distribution, the corresponding models are taken as (Saviz, 2017) F p ( z ) = ( F ( z p ) − F ( z p − 1 ) ) ( z − z p − 1 z p − z p − 1 ) n + F p − 1 , where F(z) is the material property; n is the non-negative power-law index.

Basic equation

In this part, the material local coordinate system (x 1, x 2, x 3) is utilized to describe the basic equations of 2D decagonal QC materials. To simplify the formula, the non-negative power-law index is taken as n = 1 in the following equations. QCs exhibit symmetries that are contradictory to the classical crystalline law of symmetry (Fan, 2011). Two different elementary excitations are associated with atomic motion in QCs (Hu et al., 2000; Fan, 2011): phonon and phason. Phonon displacements u i (i = 1, 2, 3) are related to the translation of atoms whereas phason displacements w k (k = 1, 2) are related to atomic rearrangements along the quasiperiodic direction. Both displacement fields are needed in the analysis of QCs and are actually coupled with each other. Based on the fundamentals of QC linear elastic theory (Hu et al., 2000), the generalized relationship of strain-displacement for 2D QCs is ε i j = ( u i , j + u j , i ) / 2 ,   w k j = w k , j , where j = 1, 2, 3. ε ij and w kj are the strains in the phonon and phason fields, respectively. The subscript comma denotes the partial differentiation with respect to the axis.

Based on Newton’s second law, the first equation of motion in Eq. 3 is related to phonon equilibrium equations. According to the structural fluctuations attributed to phonon-phason coupling, the second equation in Eq. 3 follows Bak’s model (Bak, 1985a; Bak, 1985b), as follows σ i j , j = ρ ∂ 2 u i ∂ t 2 ,   H k j , j = ρ ∂ 2 w k ∂ t 2 , where σ ij denote the phonon stresses and H kj are the phason stresses; ρ is the mass density.

According to the conventions in crystallography, all the elastic properties and thermal expansion properties possess intrinsic centrosymmetry (Hu et al., 2000). Thus, the constitutive equations with nonlocal and thermal effects for 2D decagonal QCs with point groups 10 mm, 1,022, 10 ¯ 2 m, and 10/mmm can be expressed as follows (Li et al., 2018) σ 11 = C 11 ε 11 + C 12 ε 22 + C 13 ε 33 + R 1 ( w 11 + w 22 ) − β 1 T , σ 22 = C 12 ε 11 + C 11 ε 22 + C 13 ε 33 − R 1 ( w 11 + w 22 ) − β 2 T , σ 33 = C 13 ε 11 + C 13 ε 22 + C 33 ε 33 − β 3 T , σ 23 = σ 32 = 2 C 44 ε 23 , σ 31 = σ 13 = 2 C 44 ε 13 , σ 12 = σ 21 = 2 C 66 ε 12 − R 1 w 12 + R 1 w 21 , H 11 = R 1 ( ε 11 − ε 22 ) + K 1 w 11 + K 2 w 22 , H 22 = R 1 ( ε 11 − ε 22 ) + K 1 w 22 + K 2 w 11 , H 23 = K 4 w 23 , H 12 = − 2 R 1 ε 12 + K 1 w 12 − K 2 w 21 , H 13 = K 4 w 13 , H 21 = 2 R 1 ε 12 − K 2 w 12 + K 1 w 21 , where C 11, C 12, C 13, C 33, C 44 are the elastic constants with the relationship 2C 66 = C 11 − C 12 in the phonon field; K 1, K 2, K 4, and R 1 represent the phason elastic constants, and phonon-phason coupling elastic constants, respectively; β i are the thermal constants with β 1 = β 2; T represents the variation of the temperature.

According to the heat conduction theory, the relationships of the heat fluxes q i and the temperature T (Yang et al., 2014) are q 1 = − k 11 T , 1 , q 2 = − k 22 T , 2 , q 3 = − k 33 T , 3 , where k 11, k 22, and k 33 represent the thermal conductivity coefficients.

Based on the Green–Naghdi generalized thermoelastic theory (Al-Qahtani and Datta, 2004), the extended dynamic equations for QC thermoelastic plate in the absence of heat generations (Yang et al., 2018) are q i , i − T 0 β i u i , i ∂ 2 u i ∂ t 2 = ρ C e ∂ 2 T ∂ t 2 , where T 0 is the environmental temperatures and C e denotes the specific heat capacity at constant strain.

 Formulations of state-space method

The Cartesian coordinate system (x, y, z) is utilized to describe the dynamic response of multilayered QC structures. Substituting Eq. 5 into Eq. 6, the use of the principle of conservation of energy obtains, as follows − k 11 ∂ 2 T ∂ x 2 − k 22 ∂ 2 T ∂ y 2 + ∂ q z ∂ z − T 0 ( β 1 ∂ 3 u x ∂ x ∂ t 2 + β 2 ∂ 3 u y ∂ y ∂ t 2 + β 3 ∂ 3 u z ∂ z ∂ t 2 ) = ρ C e ∂ 2 T ∂ t 2 ,    q z = − k 33 ∂ T ∂ z .

According to the state-space method and multilayered structure theory (Huang et al., 2019), the temperature T and heat flux q z in Eq. 7 can be taken as the state variables. Thus, substituting Eqs 2, 4 into Eq. 3 and combined with Eq. 7, the state equations can be derived in matrix forms as ∂ ∂ z θ 1 ( x , y , z , t ) = A θ 1 ( x , y , z , t ) , where θ 1(x, y, z, t) = [u x , u y , w x , w y , σ zz , T, σ xz , σ yz , H xz , H yz , u z , q z ]T is called the basic variables, in which the superscript “T” denotes transpose, and the state transition matrix A is A= [ 0 A 1 A 2 0 ] .

The submatrices A 1 and A 2 in Eq. 9 are A 1 = [ a 2 0 0 0 − ∂ ∂ x 0 0 a 2 0 0 − ∂ ∂ y 0 0 0 b 3 0 0 0 0 0 0 b 3 0 0 − ∂ ∂ x − ∂ ∂ y 0 0 c 8 ∂ 2 ∂ t 2 0 0 0 0 0 0 c 1 ] , A 2 = [ − a 1 ∂ 2 ∂ x 2 − a 6 ∂ 2 ∂ y 2 + c 8 ∂ 2 ∂ t 2 − ( a 3 + a 6 ) ∂ 2 ∂ x ∂ y − b 1 ∂ 2 ∂ x 2 + b 1 ∂ 2 ∂ y 2 − ( a 3 + a 6 ) ∂ 2 ∂ x ∂ y − a 6 ∂ 2 ∂ x 2 − a 1 ∂ 2 ∂ y 2 + c 8 ∂ 2 ∂ t 2 2 b 1 ∂ 2 ∂ x ∂ y − b 1 ∂ 2 ∂ x 2 + b 1 ∂ 2 ∂ y 2 2 b 1 ∂ 2 ∂ x ∂ y − b 2 ( ∂ 2 ∂ x 2 + ∂ 2 ∂ y 2 ) + c 8 ∂ 2 ∂ t 2 − 2 b 1 ∂ 2 ∂ x ∂ y − b 1 ∂ 2 ∂ x 2 + b 1 ∂ 2 ∂ y 2 0 a 4 ∂ ∂ x a 4 ∂ ∂ y 0 − c 2 c 6 ∂ 3 ∂ x ∂ t 2 − c 2 c 6 ∂ 3 ∂ y ∂ t 2 0 − 2 b 1 ∂ 2 ∂ x ∂ y a 4 ∂ ∂ x c 2 ∂ ∂ x − b 1 ∂ 2 ∂ x 2 + b 1 ∂ 2 ∂ y 2 a 4 ∂ ∂ y c 2 ∂ ∂ y 0 0 0 − b 2 ( ∂ 2 ∂ x 2 + ∂ 2 ∂ y 2 ) + c 8 ∂ 2 ∂ t 2 0 0 0 a 5 c 3 0 − c 3 c 6 ∂ 2 ∂ t 2 c 4 ∂ 2 ∂ x 2 + c 5 ∂ 2 ∂ y 2 − c 7 ∂ 2 ∂ t 2 ] , where the coefficients in Eq. 10 are written as follows a 1 = C 11 − C 13 2 C 33 ,   a 2 = 1 C 44 ,   a 3 = C 12 − C 13 2 C 33 ,   a 4 = − C 13 C 33 ,   a 5 = 1 C 33 , a 6 = C 66 ,   b 1 = R 1 ,   b 2 = K 1 ,   b 3 = 1 K 4 ,   c 1 = − 1 k 33 , c 2 = β 1 − C 13 β 3 C 33 , c 3 = β 3 C 33 , c 4 = k 11 , c 5 = k 22 , c 6 = T 0 , c 7 = ρ C e + T 0 β 3 2 C 33 , c 8 = ρ .

The eliminating in-plane variables θ 2(x, y, z, t) = [σ xx , σ yy , σ xy , H xx , H yy , H xy , H yx , q x , q y ]T can be expressed as σ x x = a 1 ∂ u x ∂ x + a 3 ∂ u y ∂ y + b 1 ∂ w x ∂ x + b 1 ∂ w y ∂ y − a 4 σ z z − c 2 T , σ y y = a 3 ∂ u x ∂ x + a 1 ∂ u y ∂ y − b 1 ∂ w x ∂ x − b 1 ∂ w y ∂ y − a 4 σ z z − c 2 T , σ x y = a 6 ∂ u x ∂ y + a 6 ∂ u y ∂ x − b 1 ∂ w x ∂ y + b 1 ∂ w y ∂ x , H x x = b 1 ∂ u x ∂ x − b 1 ∂ u y ∂ y + b 2 ∂ w x ∂ x + b 4 ∂ w y ∂ y , H y y = b 1 ∂ u x ∂ x − b 1 ∂ u y ∂ y + b 4 ∂ w x ∂ x + b 2 ∂ w y ∂ y , H x y = − b 1 ∂ u x ∂ y − b 1 ∂ u y ∂ x + b 2 ∂ w x ∂ y − b 4 ∂ w y ∂ x , H y x = b 1 ∂ u x ∂ y + b 1 ∂ u y ∂ x − b 4 ∂ w x ∂ y + b 2 ∂ w y ∂ x , q x = − c 4 ∂ T ∂ x , q y = − c 5 ∂ T ∂ y , with b 4 = K 2.

For the rectangular composite materials with ideal constraints at the edges, Simply supported (S), Clamped supported (C), and zero temperature boundary conditions are listed in the x- and y-directions as follows { S : u y = u z = w y = σ x x = H x x = 0 ,   T = 0 C : u x = u y = u z = w x = w y = 0 ,   T = 0    at  x = 0  and  x = L x ,  S : u x = u z = w x = σ y y = H y y = 0 ,   T = 0     at  y = 0  and  y = L y , where “S” indicates the simply supported boundary condition and “C” is the clamped supported boundary condition. For example, “CSCS” denotes the model with the clamped supported boundary conditions at x = 0 and x = L x , and others are simply supported boundary conditions at y = 0 and y = L y , respectively.

For the simply supported boundary conditions of the FG 2D QC coating structures in the y-direction, the Fourier series is sufficiently general to obtain the solutions, as follows [ u x ( x , y , z , t ) u y ( x , y , z , t ) w x ( x , y , z , t ) w y ( x , y , z , t ) σ z z ( x , y , z , t ) q z ( x , y , z , t ) σ x z ( x , y , z , t ) σ y z ( x , y , z , t ) H x z ( x , y , z , t ) H y z ( x , y , z , t ) u z ( x , y , z , t ) T ( x , y , z , t ) ] = ∑ l = 1 ∞ [ u ˜ x ( x , z ) sin ( q y ) u ˜ y ( x , z ) cos ( q y ) w ˜ x ( x , z ) sin ( q y ) w ˜ y ( x , z ) cos ( q y ) σ ˜ z z ( x , z ) sin ( q y ) q ˜ z ( x , z ) sin ( q y ) σ ˜ x z ( x , z ) sin ( q y ) σ ˜ y z ( x , z ) cos ( q y ) H ˜ x z ( x , z ) sin ( q y ) H ˜ y z ( x , z ) cos ( q y ) u ˜ z ( x , z ) sin ( q y ) T ˜ ( x , z ) sin ( q y ) ] e i ω t ,   [ σ x x ( x , y , z , t ) σ y y ( x , y , z , t ) σ x y ( x , y , z , t ) H x x ( x , y , z , t ) H y y ( x , y , z , t ) H x y ( x , y , z , t ) H y x ( x , y , z , t ) q x ( x , y , z , t ) q y ( x , y , z , t ) ] = ∑ l = 1 ∞ [ σ ˜ x x ( x , z ) sin ( q y ) σ ˜ y y ( x , z ) sin ( q y ) σ ˜ x y ( x , z ) cos ( q y ) H ˜ x x ( x , z ) sin ( q y ) H ˜ y y ( x , z ) sin ( q y ) H ˜ x y ( x , z ) cos ( q y ) H ˜ y x ( x , z ) cos ( q y ) q ˜ x ( x , z ) sin ( q y ) q ˜ y ( x , z ) cos ( q y ) ] e i ω t , where ω is the angular frequency, and imaginary i = − 1 ; the half-wave numbers are taken as q = lπ/L y with l being the positive integers so that all Fourier series expansion coefficients are related to the summations for l.

Incorporating Eq. 14 with Eqs 8, 12 yields ∂ ∂ z θ ˜ 1 ( x , z ) = A ˜ θ ˜ 1 ( x , z ) ,    θ ˜ 2 ( x , z ) = G ˜ θ ˜ 2 ( x , z ) , where θ ˜ 1 ( x , z ) = [ u ˜ x ,   u ˜ y ,   w ˜ x ,   w ˜ y ,   σ ˜ z z ,   T ˜ ,   σ ˜ x z ,   σ ˜ y z ,   H ˜ x z ,   H ˜ y z ,   u ˜ z ,   q ˜ z ] T ; A ˜ and G ˜ are both the coefficient matrixes; θ ˜ 2 ( x , z ) = [ σ ˜ x x ,   σ ˜ y y ,   σ ˜ x y ,   H ˜ x x ,   H ˜ y y ,   H ˜ y x ,   H ˜ x y ,   q ˜ x ,   q ˜ x ] T .

Although the Fourier series can be used to satisfy the simply supported boundary conditions in the y-direction, the coating structures with the clamped supported and mixed boundary condition cannot be solved in the x-direction. Here, DQM is used to deal with this problem. DQM is a numerical method for solving ordinary differential, partial differential, integral, and integrodifferential equations. And DQM has been approved as highly efficient for the rapid solution of various typical boundary value and initial value problems (Lü et al., 2008).

In this paper, Eq. 15 are dispersed using DQM to separate the dependence on them to the x-direction and convert the partial differentials Eq. 15 into ordinary differential equations. Here, the Chebyshev-Gauss-Lobatto grid space model (Lü et al., 2006) is adopted in the in-plane discrete direction x r = L x 2 [ 1 − cos ( r − 1 N − 1 π ) ]    ( r = 1,2,3 , ⋯ N ) , where N is the number of sampling points.

Based on the DQM, Eq. 15 can be rewritten as d u ˜ x r d z = a 2 σ ˜ x z r − ∑ k = 1 N X r k ( 1 ) u ˜ z k ,   d u ˜ y r d z = a 2 σ ˜ y z r − q u ˜ z r ,   d w ˜ x r d z = b 3 H ˜ x z r ,   d w ˜ y r d z = b 3 H ˜ y z r ,   d σ ˜ z z r d z = − c 8 ω 2 u z r − ∑ k = 1 N X r k ( 1 ) σ ˜ x z r + q σ ˜ y z r , d T ˜ r d z = c 1 q ˜ z r , d σ ˜ x z r d z = − c 8 ω 2 u x r − a 1 ∑ k = 1 N X r k ( 2 ) u ˜ x k + a 6 q 2 u ˜ x r + ( a 3 + a 6 ) q ∑ k = 1 N X r k ( 1 ) u ˜ y k − b 1 ∑ k = 1 N X r k ( 2 ) w ˜ x k − b 1 q 2 w ˜ x r + 2 b 1 q ∑ k = 1 N X r k ( 1 ) w ˜ y k + a 4 ∑ k = 1 N X r k ( 1 ) σ ˜ z z k + c 2 ∑ k = 1 N X r k ( 1 ) T ˜ k , d σ ˜ y z r d z = − ( a 3 + a 6 ) q ∑ k = 1 N X r k ( 1 ) u ˜ x k − c 8 ω 2 u y r − a 6 ∑ k = 1 N X r k ( 2 ) u ˜ y k + a 1 q 2 u ˜ y r + 2 b 1 q ∑ k = 1 N X r k ( 1 ) w ˜ x k − b 1 ∑ k = 1 N X r k ( 2 ) w ˜ y k − b 1 q 2 w ˜ y r + a 4 q σ ˜ z z r + c 2 q T ˜ r , d H ˜ x z r d z = − b 1 ∑ k = 1 N X r k ( 2 ) u ˜ x k − b 1 q 2 u ˜ x r − 2 b 1 q ∑ k = 1 N X r k ( 1 ) u ˜ y k − c 8 ω 2 w x r − b 2 ∑ k = 1 N X r k ( 2 ) w ˜ x k + b 2 q 2 w ˜ x r , d H ˜ y z r d z = − 2 b 1 q ∑ k = 1 N X r k ( 1 ) u ˜ x k − b 1 ∑ k = 1 N X r k ( 2 ) u ˜ y j − b 1 q 2 u ˜ y r − c 8 ω 2 w y r − b 2 ∑ k = 1 N X r k ( 2 ) w ˜ y k + b 2 q 2 w ˜ y r , d u ˜ z r d z = a 4 ∑ k = 1 N X r k ( 1 ) u ˜ x k − a 4 q u ˜ y r + a 5 σ ˜ z z r + c 3 T ˜ r , d q ˜ z r d z = c 2 c 6 ω 2 ∑ k = 1 N X r k ( 1 ) u ˜ x k + c 2 c 6 ω 2 q u ˜ y r + c 3 c 6 ω 2 σ ˜ z z r + c 4 ∑ k = 1 N x X i k ( 2 ) T ˜ k j − c 5 q 2 T ˜ r + c 7 ω 2 T ˜ r , σ ˜ x x r = a 1 ∑ k = 1 N X r k ( 1 ) u ˜ x k − a 3 q u ˜ y r + b 1 ∑ k = 1 N X r k ( 1 ) w ˜ x k − b 1 q w ˜ y r − a 4 σ ˜ z z r − c 2 T ˜ r , σ ˜ y y r = a 3 ∑ k = 1 N X r k ( 1 ) u ˜ x k − a 1 q u ˜ y r − b 1 ∑ k = 1 N X r k ( 1 ) w ˜ x k + b 1 q w ˜ y r − a 4 σ ˜ z z r − c 2 T ˜ r , σ ˜ x y r = a 6 q u ˜ x r + a 6 ∑ k = 1 N X r k ( 1 ) u ˜ y k − b 1 q w ˜ x r + b 1 ∑ k = 1 N X r k ( 1 ) w ˜ y k , H ˜ x x r = b 1 ∑ k = 1 N X r k ( 1 ) u ˜ x k + b 1 q u ˜ y r + b 2 ∑ k = 1 N X r k ( 1 ) w ˜ x k − b 4 q w ˜ y r , H ˜ y y r = b 1 ∑ k = 1 N X r k ( 1 ) u ˜ x k + b 1 q u ˜ y r + b 4 ∑ k = 1 N X r k ( 1 ) w ˜ x k − b 2 q w ˜ y r , H ˜ x y r = − b q u ˜ x r − b 1 ∑ k = 1 N X r k ( 1 ) u ˜ y k + b 2 q w ˜ x r − b 4 ∑ k = 1 N X r k ( 1 ) w ˜ y k , H ˜ y x r = b 1 q u ˜ x r + b 1 ∑ k = 1 N X r k ( 1 ) u ˜ y k − b 4 q w ˜ x r + b 2 ∑ k = 1 N X r k ( 1 ) w ˜ y k , q ˜ x r = − c 4 ∑ k = 1 N X r k ( 1 ) T ˜ k , q ˜ y r = − c 5 q T ˜ r , where X r k ( m ) (1≤ m ≤ N − 1) are the differential quadrature weight coefficients.

For a special problem, the boundary conditions should also be incorporated into the state Eq. 17 so as to obtain the solution. Two typical supporting conditions at x = 0 and x = L x can be rewritten as Simply supported  ( S ) :  u ˜ y d = u ˜ z d = w ˜ y d = σ ˜ x x d = H ˜ x x d = 0 ,    T ˜ = 0 , Clamped supported  ( C ) :  u ˜ x d = u ˜ y d = u ˜ z d = w ˜ x d = w ˜ y d = 0 ,    T ˜ = 0 , where d = 1 or N.

However, the phonon and phason stress boundary conditions σ ˜ x x d = H ˜ x x d = 0 cannot be introduced into Eq. 17 directly and should be considered by virtue of the expressions of state variables using Eq. 18. In this paper, three boundary conditions of SSSS, CSCS, and CSSS are considered, and the updated state equations for the QC multilayered structure with boundary conditions CSSS are listed in Eq. A1 of Appendix A.

 General solutions for 2D quasicrystal laminate

For the pth layer, the state Eq. A1 that satisfies the corresponding boundary conditions can be written as the following unified matrix form d dz δ ( p ) = T ( p ) δ ( p ) , where δ ( p ) = [ u ˜ xr ,  u ˜ yr ,  w ˜ xr ,  w ˜ yr ,  σ ˜ zzr ,  T ˜ r ,  σ ˜ xzr ,  σ ˜ yzr ,  H ˜ xzr ,  H ˜ yzr ,  u ˜ zr ,  q ˜ zr ] T ; T (p) is the coefficient matrix.

According to the theory of ordinary differential equation, the solutions of Eq. 21 are derived at the upper and lower interfaces of the pth layer δ ( p ) ( z ) =exp [ T ( p ) ( z- z p-1 ) ] δ ( p ) ( z p-1 )     ( z p-1 ≤ z ≤ z p ) .

Let z = z p in Eq. 22, we find that δ ( p ) = M ( p ) δ ( p ) , where M ( p ) =exp [ ( z p - z p-1 ) T ( p ) ] =exp [ h p T ( p ) ] .

Similarly, we get δ ( p+1 ) = M ( p+1 ) δ ( p+1 ) .

Assume that the interface z p is formed by bonding the upper surface z p + of the pth layer and the lower surface z p − of the (p+1)-th layer, and the nearby layers are perfectly linked.

The field variables are assumed to be continuous, so the interface conditions are as follows u ˜ i r ( p ) + = u ˜ i r ( p ) − ,  w ˜ k r ( p ) + = w ˜ k r ( p ) − , T ˜ r ( p ) + = T ˜ r ( p ) − , σ ˜ i z r ( p ) + = σ ˜ i z r ( p ) − , H ˜ k z r ( p ) + = H ˜ k z r ( p ) − , q ˜ z r ( p ) + = q ˜ z r ( p ) − , where the superscripts “+” and “−” are defined as variables on the upper surface of the pth layer and variables on the lower surface of the (p + 1)-th layer, respectively.

However, if the thickness h p (z − z p ) in Eq. 22 is large, or the coefficient matrix T (p) exists maxima or minima values, numerical instability will be encountered. In addition, when it comes to the cases of high-order frequencies and large discrete point numbers (Chen et al., 2003; Lü et al., 2006), the numerical instabilities will also arise by the conventional propagator matrix M (p) (h p ) in Eq. 23. Several methods have been developed to deal with this numerical difficulty. In this section, based on the joint coupling matrices (Lü et al., 2008), a new transfer relation is established to resolve this problem.

According to interface conditions, field variables can be rewritten as J c ( p ) [ τ ( p ) + τ ( p ) - ] = 0     ( p=1,2,3,…,M-1 ) , where J c = [I, −I] is the joint coupling matrix at the interface z p , and τ = [ u ˜ x r ,   u ˜ y r ,   w ˜ x r ,   w ˜ y r ,   u ˜ z r ,   q ˜ z r ,   σ ˜ x z r ,   σ ˜ y z r ,   σ ˜ z z r ,   H ˜ x z r ,   H ˜ y z r ,   T ˜ r ] T .

Similarly, we can obtain the relationship of the top and bottom surfaces for the whole QC model, as follows J t τ ( M ) − = t ( H ) ,   J w τ ( 0 ) + = t ( 0 ) , where J t = J w = [I, 0], and t = [ σ ˜ x z r , σ ˜ y z r , σ ˜ z z r , H ˜ x z r , H ˜ y z r , T ˜ r ] T .

It is assumed that the mechanical boundary conditions on the top and bottom surfaces of the model are zero, and a top surface harmonic heat temperature is considered, i.e., z = H : σ ˜ x z r = σ ˜ y z r = σ ˜ z z r = H ˜ x z r = H ˜ y z r = 0 ,   T ˜ r = T ˜ 0 e i ω t ⁡ sin ⁡ q y , z = 0 : σ ˜ x z r = σ ˜ y z r = σ ˜ z z r = H ˜ x z r = H ˜ y z r = T ˜ r = 0 , where T ˜ 0 = [ T s ⁡ sin ( π x 2 / L x ) ,   ⋯ ,   T s ⁡ sin ( π x r / L x ) ,   ⋯ ,   T s ⁡ sin ( π x N − 1 / L x ) ] T      ( 2 ≤ r ≤ N − 1 ) , with temperature amplitude T s .

From Eqs 26, 28, one gets Jτ=f, where J = diag [ J t ,   J c ,   ⋯ ,   J c ,   J w ] ,   f = [ t ( H ) ,   0 ,   ⋯ ,   0 ,   t ( 0 ) ] T , τ = [ ( τ ( M ) − ) T ,   ( τ ( M − 1 ) + τ ( M − 1 ) − ) T ,   … , ( τ ( 1 ) + τ ( 1 ) − ) T ,   ( τ ( 0 ) + ) T ] T .

τ can be rewritten as τ = [ ( τ ( M ) − τ ( M − 1 ) + ) T ,   ( τ ( M − 1 ) − τ ( M − 2 ) + ) T ,   … ,   ( τ ( 2 ) − τ ( 1 ) + ) T ,   ( τ ( 1 ) − τ ( 0 ) + ) T ] T .

The whole structure can be defined as τ=N τ + , where N = diag [ ( M ( M ) I ) ,   ( M ( M − 1 ) I ) ,   ⋯ ,   ( M ( 2 ) I ) ,   ( M ( 1 ) I ) ] , τ − = [ ( τ ( M − 1 ) + ) T ,   ( τ ( M − 2 ) + ) T ,   ⋯ ,   ( τ ( 1 ) + ) T ,   ( τ ( 0 ) + ) T ] T .

Substituting Eq. 34 into Eq. 31 results in JNτ=f .

Eq. 36 can be rewritten as τ= [ JN ] -1 f , where the superscript “−1” represents the inverse of the matrix.

Using Eqs 37, 18, the forced vibration solutions of the whole QC coating structure are obtained.

 Numerical examples

Typical illustrative examples of the fully coupled thermoelasticity theory in multilayered QC structures are provided for forced vibration analysis. The multilayered FG QC coating model is composed of three single layers of the same thickness, and the dimensions of this laminate are L x × L y × H = 1 mm × 1 mm × 0.3 mm. Two kinds of QC materials (called QC1 and QC2) are considered for every single layer, and the thermoelasticity material properties of QC1 (top) and QC2 (bottom) are listed in Table 1. And two kinds of material parameters completely match the elastic deformation energy density of the QC (Hwu, 2010; Fan, 2011), so they can be used to calculate the relative changes of the dynamic analysis for the 2D FG QC structures. In addition, the temperature amplitude T s is taken as 1 K. The following discussions are based on calculations with l = 1, and the discrete points are taken as N = 25. For the convenience of expression, the natural frequencies of the layered structures are considered to be dimensionless as Ω = ωL max/ C m a x / ρ m a x , where L max is the maximum lateral dimension of the structure model; C max and ρ max represent the maximum elastic constant and density in this model, respectively.

Material properties (C, K, and R in 109 N/m2, β in 106 N K/m2, k in W K/m, T 0 in K, C e in J kg/K, ρ in kg/m3).

C 11 C 12 C 13 C 33 C 44 K 1 K 2 K 4 R 1
QC1 234.33 57.41 66.63 232.22 70.19 122 24 12 8.846
QC2 166 77 78 162 43 65 13 6 5.417
β 1 β 3 k 11 k 22 k 33 T 0 C e ρ
QC1 1.383 1.798 5.3 5.3 6.89 293 1.187 4,186
QC2 1.644 2.226 3.672 3.672 4.972 293 1.228 5,800

It is assumed that the tractive force and harmonic heat temperature on the top and bottom surfaces of the model are zero. The first-order dimensionless natural frequencies Ω for the QC1/QC1/QC1 plate with boundary conditions SSSS is Ω = 0.33471. In order to the resonance of the QC laminates, the dimensionless frequencies for force vibration response are taken as Ω = 0.05, 0.1, 0.15, and 0.2. It is noted that the static closed-form solutions of this model can be obtained at Ω = 0.

One numerical example of the QC plate with the boundary condition SSSS to verify the validity and accuracy of the proposed method and the numerical solution is presented in Figure 2. The material parameters, the shape and size of the plate, the quasicrystal quasi-periodic direction, and the mechanical quantities used are consistent with those in Ref. (Yang et al., 2018). The phonon and phason displacements and stresses u y , w y , σ yy , and H yy are normalized by their maximum values among these four variables along the thickness direction. The results obtained by the pseudo-Stroh formalism (Yang et al., 2018) and the presented method are consistent. This feature indicates that the method in this paper has high precision and good convergence and the solutions are stable.

Normalized displacements and stresses for the QC plate.

 Functionally graded quasicrystal coating structures with boundary condition CSCS

In this part, the solutions of the FG QC coating structures with boundary conditions CSCS are presented. The power law index is set as n = 0.2, 0.6, 1, 2, and 6. The frequency is taken as Ω = 0.1. The horizontal coordinate is fixed at (x, y) = (0.25L x , 0.75L y ).

Figure 3 illustrates the effects of the power law index n on the phonon and phason displacement components of the structures along the thickness direction. Regarding the displacement solutions induced by the time harmonic temperature loadings only, where a complete reversal variation of u y and w y occurs in Figures 3A,D with increasing n, both displacement components are sensitive to n. The distribution of u z (Figure 3B) internal model is nonlinear rather than functional, which is consistent with the elastic theory of models. w x (Figure 3C) are more sensitive to n than w y , and the value of w x is greater than that of w y .

Variation of the phonon and phason displacements for QC coating structures: (A) u x , (B) u z , (C) w x , and (D) w y .

Figure 4 illustrates the effects of the power law index n on the phonon and phason stress components of the QC coating structures along the thickness direction. The solution strictly satisfies the basic equation of 2D QC material, and the distributions of field variables at any position of the multilayered structure can also be presented. Since each layer of the coating structure is inhomogeneous, the distributions of σ yz and σ zz along the z-direction are not symmetrical in Figures 4A,B. σ xy , H yx , and H xy (Figures 4C,E,F) are similar between layers, which are discontinuous between internal interfaces. Based on the classic elasticity theory, only the local stress of the basic equation of the QC is considered. But in fact, the stress state also includes the strong interlayer stress between the interfaces. The high-interlayer stress is considered to be one of the special failure mechanisms of composite materials in engineering applications. In addition, if the material properties at the interfaces are the same, these stresses are continuous for these structure models. The distribution of H yz (Figure 4D) is nonlinear, and the value is zero at the external surfaces. For the above physical variables, if DQM will be also utilized to discretize domains along the y-direction, the superposition l will not appear in the formulations of the exact solution.

Variation of the phonon and phason stresses for QC coating structures: (A) σ yz , (B) σ zz , (C) σ xy , (D) H yz , (E) H yx , and (F) H xy .

Figure 5 illustrates the effects of the power law index n on the temperature and heat flux of the QC coating structures along the thickness direction. The vertical heat flux q z (Figure 5B) for multilayers is more sensitive to n than the temperature T (Figure 5A). This feature is attributed to q z which is affected by a linear material constant k 33 in Eq. 5. In addition, n has a significant effect on q z and does not change its magnitude and direction at the top and bottom surfaces of the QC coating structures.

Variation of the temperature and heat flux for QC coating structures: (A) T and (B) q z .

 Functionally graded quasicrystal coating structures with boundary condition CSSS

In this part, the solutions of the FG QC coating structures with boundary conditions CSSS along the thickness direction are presented. The power law index is set as n = 1. The frequencies are taken as Ω = 0.05, 0.1, 0.15, and 0.2.

Figures 6–8 present the contour plots of u z , w x , and w y for the coating structures on the x-z plane with y = 0.75L y , respectively. The maximum magnitudes of u z (Figures 6A–D), w x (Figures 7A–D), and w y (Figures 8A–D) exist near the top surface and increase with the increase of Ω. And Ω has a slight effect on u z , w x , and w y . Due to the gradient distribution of the material within the coating structures, the variation of these displacements of this model is not symmetric at x = 0.5 mm. Different from u z and w y , the value of w x is not zero at x = L x . This feature is consistent with the linear elastic theory of models. In addition, the variation of u z , w x , and w y can also prove that the values of any position in the model can be obtained.

Effect of Ω on the phonon displacement u z (10−7 mm): (A) Ω = 0.05, (B) Ω = 0.1, (C) Ω = 0.15, and (D) Ω = 0.2.

Effect of Ω on the phason displacement w x (10−9 mm): (A) Ω = 0.05, (B) Ω = 0.1, (C) Ω = 0.15, and (D) Ω = 0.2.

Effect of Ω on the phonon displacement w y (10−9 mm): (A) Ω = 0.05, (B) Ω = 0.1, (C) Ω = 0.15, and (D) Ω = 0.2.

Figures 9–11 present the contour plots of T, q x , and q y for the coating structures on the x-z plane with y = 0.75L y , respectively. The distributions of T, q x , and q y are symmetric at x = 0.5 mm, and these variables are not affected by boundary conditions. T (Figures 9A–D) and q x (Figures 10A–D) obviously change with the different Ω, and they are slightly sensitive to the values of Ω. q x and q y (Figures 11A–D) are discontinuous at the internal interfaces, and interlaminar heat fluxes are little affected by Ω. Interlaminar thermal stress and interlaminar heat fluxes are attributed to the thermophysical and mechanical properties mismatching between the sublayers of the QC coating structures under the time-harmonic temperature loadings, which would inevitably involve interlayer failure of the laminate and reduce the carrying capacity.

Effect of Ω on the temperature T (K): (A) Ω = 0.05, (B) Ω = 0.1, (C) Ω = 0.15, and (D) Ω = 0.2.

Effect of Ω on the heat flux q x (10 W/m2): (A) Ω = 0.05, (B) Ω = 0.1, (C) Ω = 0.15, and (D) Ω = 0.2.

Effect of Ω on the heat flux q y (10 W/m2): (A) Ω = 0.05, (B) Ω = 0.1, (C) Ω = 0.15, and (D) Ω = 0.2.

 Functionally graded quasicrystal coating structures with boundary conditions SSSS, CSCS, and CSSS

In this part, the solution of the FG QC coating structures with boundary conditions SSSS, CSCS, and CSSS along the thickness direction is presented. The power law index is set as n = 1. The frequency is taken as Ω = 0.1. The horizontal coordinate is fixed at (x, y) = (0.25L x , 0.75L y ).

Figure 12 illustrates the effects of the different boundary conditions on the phonon and phason displacement components of the coating structures along the thickness direction. It can be observed that the different boundary conditions have an obvious influence on u y , u z , w x , and w y (Figures 12A–D). The maximum magnitudes of u z exist at the top surface of the coating structure. And the values of u z decrease with the increase of clamped-supported boundary conditions. This feature indicates that the clamped-supported boundary conditions have a great effect on u z , and the rigidity of the coating structure is constantly getting stronger.

Variation of the phonon and phason displacements for QC coating structures: (A) u x , (B) u z , (C) w x , and (D) w y .

Figure 13 illustrates the effects of the different boundary conditions on the phonon and phason stress components of the coating structures along the thickness direction. The distributions of σ xz , σ zz , σ xx , H xz , H xx , and H xy (Figures 13A–F) are considerably altered by the clamped supported boundary conditions. It is noted that the values of σ xz , σ zz , and H xz on the bottom and top surfaces are zero (i.e., free tractions for thermal loadings), which validates the present formulation and the numerical results by the surface boundary conditions. The interlaminar stress values of σ xx and H xx with boundary conditions SSSS are greater than that of σ xx and H xx with CSCS and CSSS, while σ xy and H xy are the opposite.

Variation of the phonon and phason stresses for QC coating structures: (A) σ xz , (B) σ zz , (C) σ xx , (D) H xz , (E) H xx , and (F) H xy .

Figure 14 illustrates the effects of the different boundary conditions on T and q z of the coating structures along the thickness direction. The values and distributions of T and q z (Figures 14A,B) are insensitive to different boundary conditions, and the heat fluxes q x and q y possess the same feature as q z .

Variation of the temperature and heat flux for QC coating structures: (A) T and (B) q z .

 Conclusion

The exact solutions for time-harmonic temperature loadings in multilayered anisotropic FG QC coating structures with mixed supported boundary conditions are derived in a thermal environment. While different boundary conditions for each inhomogeneous layer are satisfied using the DQM and the Fourier series to obtain the general field solutions, the multilayered feature is globally handed by improving the traditional propagation matrix method with joint coupling matrices. The present framework admits different mechanical and thermal boundary conditions to derive the exact solutions and can be applied to thick and thin laminated structures, without producing numerical instability issues from high frequency and/or large discrete point numbers. Application illustrations are proposed to throw light on various effects of the power-law index n, frequency Ω, and different boundary conditions on the thermoelastic fields in three multilayered structures, namely:

1) Different from the conventional propagator matrix method, the global propagator relation is reestablished to resolve numerical instabilities in the case of large discrete point numbers and high frequencies for QC coating structures. Numerical examples prove that this method has high precision and good convergence at N = 25.

2) Depending on the input frequency amplitude, severe oscillating displacements and stresses take place on each layer and interface that can endanger the safety-related strength and stiffness of the three-layer coatings. The remarkable change in the magnitude and distributions of the induced thermal stresses σ xx , σ yy , σ xy , H xx , H xy , H yy , H yx and heat fluxes q x , q y with increasing boundary conditions and frequency amplitude is exhibited.

3) The internal displacement, stress, and temperature fields due to time-harmonically forced vibration for QC coating structures are described using the power law index. The field variables are sensitive to the power law index except for temperature T, and their distributions are different with the change of n.

Overall, the methods and numerical results in this paper can be utilized to validate the accuracy of other numerical methods and serve for the analysis and design of multilayered QC samples with desired time-harmonic thermoelastic responses.

 Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

 Author contributions

XF wrote the first draft of the manuscript, and contributed to conception and design of the study. LZ performed the statistical analysis. ZZ and YL processed data and pictures. YG revised manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

 Funding

Funding: This work was supported by the National Natural Science Foundation of China (grant numbers 11972365, 12102458, and 12102481) and the China Agricultural University Education Foundation (grant numbers 1101-240001).

 Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

 Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

 References Al-Qahtani H. Datta S. K. (2004). Thermoelastic waves in an anisotropic infinite plate. J. Appl. Phys. 96 (7), 3645–3658. 10.1063/1.1776323 Bak P. (1985). Phenomenological theory of icosahedral incommensurate (quasiperiodic) order in Mn-Al alloys. Phys. Rev. Lett. 54 (14), 1517–1519. 10.1103/physrevlett.54.1517 Bak P. (1985). Symmetry, stability, and elastic properties of icosahedral incommensurate crystals. Phys. Rev. B 32 (9), 5764–5772. 10.1103/physrevb.32.5764 Beardsley M. B. (2008). Potential use of quasicrystalline materials as thermal barrier coatings for diesel engine components[D]. Iowa State University. Castaings M. Hosten B. (1994). Delta operator technique to improve the Thomson–Haskell‐method stability for propagation in multilayered anisotropic absorbing plates. J. Acoust. Soc. Am. 95 (4), 1931–1941. 10.1121/1.408707 Chen W. Q. Lü C. F. Bian Z. G. (2003). Elasticity solution for free vibration of laminated beams. Compos. Struct. 62 (1), 75–82. 10.1016/s0263-8223(03)00086-2 Dolinsek J. (2012). Electrical and thermal transport properties of icosahedral and decagonal quasicrystals. Chem. Soc. Rev. 41 (20), 6730–6744. 10.1039/c2cs35036j Dunkin J. W. (1965). Computation of modal solutions in layered, elastic media at high frequencies. Bull. Seismol. Soc. Am. 55 (2), 335–358. 10.1785/bssa0550020335 Fan T. Y. (2011). Mathematical theory of elasticity of quasicrystals and its applications[M]. Heidelberg: Springer. Guo J. H. Yu J. Xing Y. M. Pan E. Li L. (2016). Thermoelastic analysis of a two-dimensional decagonal quasicrystal with a conductive elliptic hole. Acta Mech. 227 (9), 2595–2607. 10.1007/s00707-016-1657-7 Guo J. H. Zhang M. Chen W. Q. Zhang X. (2020). Free and forced vibration of layered one-dimensional quasicrystal nanoplates with modified couple-stress effect. Sci. China Phys. Mech. Astron. 63 (7), 274621–274622. 10.1007/s11433-020-1547-3 Hu C. Z. Wang R. H. Ding D. H. (2000). Symmetry groups, physical property tensors, elasticity and dislocations in quasicrystals. Rep. Prog. Phys. 63 (1), 1–39. 10.1088/0034-4885/63/1/201 Huang Y. Z. Li Y. Yang L. Z. Gao Y. (2019). Static response of functionally graded multilayered one-dimensional hexagonal piezoelectric quasicrystal plates using the state vector approach. J. Zhejiang Univ. Sci. A 20 (2), 133–147. 10.1631/jzus.a1800472 Huang Y. Z. Zhao M. Feng M. L. (2021). Electric–elastic analysis of multilayered two-dimensional decagonal quasicrystal circular plates with simply supported or clamped boundary conditions. Math. Mech. Solids 26 (9), 1337–1353. 10.1177/1081286520981618 Hwu C. (2010). Anisotropic elastic plates[M]. Boston: Springer US. Li L. H. Liu G. T. (2018). Study on a straight dislocation in an icosahedral quasicrystal with piezoelectric effects. Appl. Math. Mech. 39 (9), 1259–1266. 10.1007/s10483-018-2363-9 Li X. Y. Wang T. Zheng R. F. Kang G. (2015). Fundamental thermo-electro-elastic solutions for 1D hexagonal QC. Z. Angew. Math. Mech. 95 (5), 457–468. 10.1002/zamm.201300095 Li X. Y. (2013). Fundamental solutions of penny-shaped and half-infinite plane cracks embedded in an infinite space of one-dimensional hexagonal quasi-crystal under thermal loading. Proc. R. Soc. A 469 (2154), 20130023. 10.1098/rspa.2013.0023 Li Y. Yang L. Z. Gao Y. (2018). Thermo-elastic analysis of functionally graded multilayered two-dimensional decagonal quasicrystal plates. Z. Angew. Math. Mech. 98 (9), 1585–1602. 10.1002/zamm.201700371 Li Y. Yang L. Z. Zhang L. L. Gao Y. (2021). Nonlocal free and forced vibration of multilayered two-dimensional quasicrystal nanoplates. Mech. Adv. Mater. Struct. 28 (12), 1216–1226. 10.1080/15376494.2019.1655687 Lü C. F. Chen W. Q. Shao J. W. (2008). Semi-analytical three-dimensional elasticity solutions for generally laminated composite plates. Eur. J. Mech. - A/Solids 27 (5), 899–917. 10.1016/j.euromechsol.2007.12.002 Lü C. Lee Y. Lim C. Chen W. (2006). Free vibration of long-span continuous rectangular Kirchhoff plates with internal rigid line supports. J. Sound. Vib. 297 (1-2), 351–364. 10.1016/j.jsv.2006.04.007 Maciá E. (2004). Compatibility factor of segmented thermoelectric generators based on quasicrystalline alloys. Phys. Rev. B 70 (10), 100201. 10.1103/physrevb.70.100201 Mal A. K. (1988). Wave propagation in layered composite laminates under periodic surface loads. Wave motion 10 (3), 257–266. 10.1016/0165-2125(88)90022-4 Matsukawa M. Yoshizawa M. Noto K. Yokoyama Y. Inoue A. (1999). Anisotropic thermal transport of 2D quasicrystals of decagonal Al–Ni–Co system. Phys. B Condens. Matter 263 (1-4), 146–148. 10.1016/s0921-4526(98)01469-0 Mora J. García P. Muelas R. Aguero A. (2020). Hard quasicrystalline coatings deposited by HVOF thermal spray to reduce ice accretion in aero-structures components. Coatings (Basel). 10 (3), 290. 10.3390/coatings10030290 Polishchuk S. Boulet P. Mézin A. de Weerd M. C. Weber S. Ledieu J. (2012). Residual stress in as-deposited Al–Cu–Fe–B quasicrystalline thin films. J. Mat. Res. 27 (5), 837–844. 10.1557/jmr.2011.415 Saviz M. R. (2017). Electro-elasto-dynamic analysis of functionally graded cylindrical shell with piezoelectric rings using differential quadrature method. Acta Mech. 228 (5), 1645–1670. 10.1007/s00707-016-1746-7 Sun T. Y. Guo J. H. Pan E. (2021). Nonlocal vibration and buckling of two-dimensional layered quasicrystal nanoplates embedded in an elastic medium. Appl. Math. Mech. 42 (8), 1077–1094. 10.1007/s10483-021-2743-6 Tan E. L. (2006). Hybrid compliance-stiffness matrix method for stable analysis of elastic wave propagation in multilayered anisotropic media. J. Acoust. Soc. Am. 119 (1), 45–53. 10.1121/1.2139617 Tan E. L. (2011). Recursive asymptotic hybrid matrix method for acoustic waves in multilayered piezoelectric media. Open J. Acoust. 1 (2), 27–33. 10.4236/oja.2011.12004 Vattré A. Pan E. Chiaruttini V. (2021). Free vibration of fully coupled thermoelastic multilayered composites with imperfect interfaces. Compos. Struct. 259, 113203. 10.1016/j.compstruct.2020.113203 Vattré A. Pan E. (2021). Thermoelasticity of multilayered plates with imperfect interfaces. Int. J. Eng. Sci. 158, 103409. 10.1016/j.ijengsci.2020.103409 Yang L. Z. Li Y. Gao Y. Pan E. (2018). Three-dimensional exact thermo-elastic analysis of multilayered two-dimensional quasicrystal nanoplates. Appl. Math. Model. 63 (11), 203–218. 10.1016/j.apm.2018.06.050 Yang L. Z. Zhang L. L. Song F. Gao Y. (2014). General solutions for three-dimensional thermoelasticity of two-dimensional hexagonal quasicrystals and an application. J. Therm. Stresses 37 (3), 363–379. 10.1080/01495739.2013.869149 Zhang L. Guo J. H. Xing Y. M. (2019). Nonlocal analytical solution of functionally graded multilayered one-dimensional hexagonal piezoelectric quasicrystal nanoplates. Acta Mech. 230 (5), 1781–1810. 10.1007/s00707-018-2344-7 Zhang M. Guo J. H. Li Y. S. (2022). Bending and vibration of two-dimensional decagonal quasicrystal nanoplates via modified couple-stress theory. Appl. Math. Mech. 43 (3), 371–388. 10.1007/s10483-022-2818-6 Appendix A: This appendix contains some equations cited in the text of the present paper.

The state equations for CSSS are d u ˜ x r d z = a 2 σ ˜ x z r − ∑ k = 2 N − 1 X r k ( 1 ) u ˜ z k   ( 2 ≤ r ≤ N ) ,   d u ˜ y r d z = a 2 σ ˜ y z r − q u ˜ z r   ( 2 ≤ r ≤ N − 1 ) , d w ˜ x r d z = b 3 H ˜ x z r   ( 2 ≤ r ≤ N ) ,   d w ˜ y r d z = b 3 H ˜ y z r   ( 2 ≤ r ≤ N − 1 ) ,   d σ ˜ z z r d z = − c 8 ω 2 u z r − ∑ k = 2 N X r k ( 1 ) σ ˜ x z k + q σ ˜ y z r − 1 a 2 ∑ k = 2 N − 1 f 1 r k u ˜ z k   ( 2 ≤ r ≤ N − 1 ) ,   d T ˜ r d z = c 1 q ˜ z r   ( 2 ≤ r ≤ N − 1 ) , d σ ˜ x z r d z = − c 8 ω 2 u x r − a 1 ∑ k = 2 N X r k ( 2 ) u ˜ x k − a 4 2 a 5 ∑ k = 2 N f 1 r k u ˜ x k + a 6 q 2 u ˜ x r + a 1 ∑ k = 2 N f N r k u ˜ x k + ( a 3 + a 6 ) q ∑ k = 2 N − 1 X r k ( 1 ) u ˜ y k − b 1 ∑ k = 2 N − 1 X r k ( 2 ) w ˜ x k − b 1 q 2 w ˜ x r + b 1 ∑ k = 2 N f N r k w ˜ x k + 2 b 1 q ∑ k = 2 N X r k ( 1 ) w ˜ y k + a 4 ∑ k = 2 N − 1 X r k ( 1 ) σ ˜ z z k + c 2 ∑ k = 2 N − 1 X r k ( 1 ) T ˜ k   ( 2 ≤ r ≤ N ) , d σ ˜ y z r d z = − ( a 3 + a 6 ) q ∑ k = 2 N X r k ( 1 ) u ˜ x k − c 8 ω 2 u y r − a 6 ∑ k = 2 N − 1 X r k ( 2 ) u ˜ y k + a 1 q 2 u ˜ y r + 2 b 1 q ∑ k = 2 N X r k ( 1 ) w ˜ x k − b 1 ∑ k = 2 N − 1 X r k ( 2 ) w ˜ y k − b 1 q 2 w ˜ y r + a 4 q σ ˜ z z r + c 2 q T ˜ r   ( 2 ≤ r ≤ N − 1 ) , d H ˜ x z r d z = − b 1 ∑ k = 2 N X r k ( 2 ) u ˜ x k − b 1 q 2 u ˜ x r − 2 b 1 q ∑ k = 2 N − 1 X r k ( 1 ) u ˜ y k − c 8 ω 2 w x r − b 2 ∑ k = 2 N X r k ( 2 ) w ˜ x k + b 2 q 2 w ˜ x r   ( 2 ≤ r ≤ N ) , d H ˜ y z r d z = − 2 b 1 q ∑ k = 2 N X r k ( 1 ) u ˜ x k − b 1 ∑ k = 2 N − 1 X r k ( 2 ) u ˜ y k − b 1 q 2 u ˜ y r − c 8 ω 2 w y r − b 2 ∑ k = 2 N − 1 X r k ( 2 ) w ˜ y k + b 2 q 2 w ˜ y r   ( 2 ≤ r ≤ N − 1 ) , d u ˜ z r d z = a 4 ∑ k = 2 N X r k ( 1 ) u ˜ x k − a 4 q u ˜ y r + a 5 σ ˜ z z r + c 3 T ˜ r   ( 2 ≤ r ≤ N − 1 ) , d q ˜ z r d z = c 2 c 6 ω 2 ∑ k = 2 N X r k ( 1 ) u ˜ x k + c 2 c 6 ω 2 q u ˜ y r + c 3 c 6 ω 2 σ ˜ z z r + c 4 ∑ k = 2 N − 1 X i k ( 2 ) T ˜ k j − c 5 q 2 T ˜ r + c 7 ω 2 T ˜ r   ( 2 ≤ r ≤ N − 1 ) , with f 1 r k = X r 1 ( 1 ) X 1 k ( 1 ) ,   f N r k = X r N ( 1 ) X N r ( 1 ) .