Anisotropy of multi-layered structure with sliding and bonded interlayer conditions

Abstract

A better understanding of the mechanical behavior of the multi-layered structure under external loading is the most important item for the structural design and the risk assessment. The objective of this study are to propose and develop an analytical solution for the mechanical behaviors of multi-layered structure generated by axisymmetric loading, and to investigate the impact of anisotropic layers and interlayer conditions on the multi-layered structure. To reach these objectives, first, according to the governing equations, the analytical solution for a single layer was formulated by adopting the spatial Hankel transform. Then the global matrix technique is applied to achieve the analytical solution of multi-layered structure in Hankel domain. The sliding and bonded interlayer conditions were considered in this process. Finally, the numerical inversion of integral transform was used to solve the components of displacement and stress in real domain. Gauss-Legendre quadrature is a key scheme in the numerical inversion process. Moreover, following by the verification of the proposed analytical solution, one typical three-layered flexible pavement was applied as the computing carrier of numerical analysis for the multi-layered structure. The results have shown that the anisotropic layers and the interlayer conditions significantly affect the mechanical behaviors of the proposed structure.

Keywords

multi-layered structure / Hankel transformation / anisotropic / transversely isotropic / interlayer condition / Gauss-Legendre quadrature

Introduction

The multi-layered structure is universal in nature ranging from macro-scale to micro-scale [¹–⁸]. For example, the flexible or rigid pavements constructed on soil-foundation for highway to sustain vehicle loading is a traditional macro-scaled multi-layered structure [¹,⁹–¹²], while the clay nano-composites is one type of micro-scaled multi-layered structure [¹³–¹⁶]. A better understanding of the mechanical behavior of the multi-layered pavement under external loading is the most important item for the structural design and the risk assessment (e.g., deformation and stress). Multiple models and calculation methods have been applied to simulate the mechanical behavior of the multi-layered structure [¹⁷–²⁰]. In respect of analytical methods, Lefeuve-Mesgouez et al. [²¹] applied the wave equation and the Fourier Transformation to study the dynamic response of an elastic structure subjected to a moving load. Hung and Yang et al. [²²] used the Fourier Transformation to analyze the response of a viscoelastic structure subjected to four types of moving loads, and studied the mechanism of wave propagation in this dynamic model. In addition to the analytical solution, the numerical methods such as finite element method (FEM) and boundary element method (BEM) were also extensively applied. Zhai and Song [²³] developed a three dimensional moving FEM in moving coordinates and analyzed the vibrations of a railway ground system under moving loads.

Almost all materials applied in the multi-layered structure exhibit apparent anisotropic properties [²⁴–²⁸]. The anisotropic property is a prerequisite of direction-dependent, which implies different performance in different directions. Generally, asphalt mixture and macadam are the main components of multi-layered flexible pavement [²⁹]. Asphalt mixture and macadam are composed of differently shaped granular materials, and its inner structure presented apparent anisotropic properties [³⁰–³³]. The simplest anisotropic model is of a transverse isotropy, which is an assumption symmetry that appraises the material to behave differently in one direction from the others [²,³⁴]. The simplification comes from this assumption is equivalent to stating that some coefficients of an anisotropic material are equal or linearly dependent. Thus, for the mechanical analysis of anisotropic media, the anisotropic properties of materials could be approximated as transverse isotropy, while transverse isotropy can be observed as a special case of anisotropy [³⁵]. Kim et al. [³⁶] developed a procedure to account for the effects of aggregate physical properties, i.e., gradation and shape characteristics, in predicting the transversely isotropic of unbound granular layers. The analysis results indicated that the aggregate physical properties significantly affect the transversely isotropic resilient behavior of unbound granular bases, and this transversely isotropic resilient behavior has substantial effect on the pavement response. Tutumluer et al. [³⁷] reported validated model for predicting field performance of flexible pavement considering isotropic and transversely isotropic properties of unbound aggregate base (UAB) materials. The validated modeling approach was accomplished by analyzing conventional flexible pavement test sections using GT-PAVE finite element program to predict responses to load in the UAB layer and comparing these predicted responses to measured data. Masad et al. [³⁸] analyzed the influence of characterizing pavement layers as isotropic and anisotropic on the predictions of permanent deformation and fatigue performance of flexible pavements. The results showed that the anisotropic behavior of pavement layers explained parts of shift and calibration factors used to relate laboratory measurements to field performance. Wang et al. [³⁹] studied the stress fields of transversely isotropic and isotropic flexible pavement under vehicle loading for the several cases. Mukhopadhyay [⁴⁰] analyzed the response of a transversely isotropic layer lying on a rigid foundation under a moving load. Yan et al. [⁴¹] reported an analytical solution considering the transversely isotropic properties of structure layers to analyze the dynamic response of flexible pavements under impact loads. You et al. [²,¹⁹,⁴²,⁴³] provided three analytical solutions by spectral element method (SEM) to analyze the mechanical behaviors of multi-layered flexible pavement subjected to static and dynamic loads. These results show that the influence of transversely isotropic properties is significant for the mechanical behavior of flexible pavements. Therefore, the anisotropic properties of the structural layer is important for the structure design and the risk assessment of multi-layered flexible pavement.

In addition, the multi-layered structure is proposed to effectively transfer and distribute loads to the bearing layers. Therefore, the interlayer condition properties between the layers also have a significant effect on the mechanical behavior of the multi-layered structure [⁴⁴–⁴⁷]. The appropriate interlayer condition will help the different layers to combine together to resist loading by providing the bond between layers, while the lack of proper interlayer conditions result in slippage and reduced shear strength between layers in the multi-layered structure [⁴⁸]. Stenzler and Goulbourne [⁴⁹] studied the mechanical behavior of several polyacrylate adhesives with varying microstructure. Chupin et al. [⁵⁰] calculated the mechanical behavior of a viscoelastic layered structure to a moving load by considering the interlayer slip. You et al. [¹⁹,⁴³,⁵¹,⁵²] employed the interlayer conditions, sliding and bonded, in an analysis for the multi-layered structure under different loads. The result from their study apprises the obvious differences of the mechanical behavior of structures with various different interlayer conditions. Therefore, during the mechanical behavior analysis of multi-layered flexible pavement, the interlayer conditions between the structure layers also should also be considered.

The objective of this study is to develop an analytical solution to evaluate the impact of the anisotropic layer and interlayer conditions on the mechanical behavior of multi-layered structures subjected to axisymmetric loading. However, in order to simplify the derivation processes, the anisotropy is approximated as transverse isotropy.

Model development and derivation

Governing equations

The relationships of strain and stress for the transversely isotropic elastic multi-layered structure can be presented in the transversely isotropic constitutive equation [²,⁵³,⁵⁴], as shown in Eq. (1):

(1)\left\{\begin{array}{l}{\sigma }_r(r,z)\\ {\sigma }_{\theta }(r,z)\\ {\sigma }_z(r,z)\\ {\tau }_{rz}(r,z)\end{array}\right\}=\left[\begin{array}{llll}{c}_{\mathrm{11}}& c_{\mathrm{12}}& c_{\mathrm{13}}& \mathrm{0}\\ c_{\mathrm{21}}& c_{\mathrm{22}}& c_{\mathrm{23}}& \mathrm{0}\\ c_{\mathrm{31}}& c_{\mathrm{32}}& c_{\mathrm{33}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& c_{\mathrm{44}}\end{array}\right]\left\{\begin{array}{l}{\epsilon }_r(r,z)\\ {\epsilon }_{\theta }(r,z)\\ {\epsilon }_z(r,z)\\ {\gamma }_{rz}(r,z)\end{array}\right\},

where c_ij is the elastic coefficient, and it can be rewritten as Eq. (2). r, q, and z are the radial, circumferential, and axial coordinates, respectively. σ_r, σ_θ, σ_z, and τ_rz are normal stress in the radial, normal stress in the circumferential, normal stress in the axial coordinate, and shear stress in the planes normal to the plane of transverse isotropy, respectively. ε_r, ε_θ, ε_z, and γ_rz are the components of strain corresponding to the above stresses.

(2)\begin{array}{l}{c}_{\mathrm{11}}=c_{\mathrm{22}}=\kappa n(\mathrm{1}-n{\mu }_{\mathrm{vh}}^{\mathrm{2}}),c_{\mathrm{12}}=c_{\mathrm{21}}=\kappa n({\mu }_{\mathrm{h}}+n{\mu }_{\mathrm{v}}^{\mathrm{2}}),\\ c_{\mathrm{13}}=c_{\mathrm{23}}=c_{\mathrm{31}}=c_{\mathrm{32}}=\kappa n{\mu }_{\mathrm{vh}}({\mu }_{\mathrm{h}}+\mathrm{1}),c_{\mathrm{33}}=\kappa (\mathrm{1}-{\mu }_{\mathrm{vh}}^{\mathrm{2}}),\\ c_{\mathrm{44}}=G_{\mathrm{v}},n=E_{\mathrm{h}}/E_{\mathrm{v}},\\ \kappa =E_{\mathrm{v}}/(\mathrm{1}+{\mu }_{\mathrm{vh}})(\mathrm{1}-{\mu }_{\mathrm{vh}}-\mathrm{2}n{\mu }_{\mathrm{v}}^2), \end{array}

where E_v, E_h, and G_v are the Young’s modulus in the vertical direction, Young’s modulus in the horizontal direction, and shear modulus in the planes normal to the plane of transverse isotropy, respectively. μ_v and μ_vh are the Poisson’s ratios in the vertical and horizontal direction, respectively. And the relationship between strain and displacement for the transversely isotropic elastic multi-layered medium are shown in Eq. (3), and the governing equations under the static load are rewritten in Eq. (4). Equation (5) was obtained by introducing Eqs. (1) and (3) into Eq. (4), where r is the material density of structure.

(3){\epsilon }_r={\displaystyle \frac{\partial u_r}{\partial r}},{\epsilon }_{\theta }={\displaystyle \frac{u_r}{r}},{\epsilon }_z={\displaystyle \frac{\partial u_z}{\partial z}},{\gamma }_{rz}={\displaystyle \frac{\partial u_r}{\partial z}}+{\displaystyle \frac{\partial u_z}{\partial r}},

(4)\begin{array}{l}{\displaystyle \frac{\partial {\sigma }_r}{\partial r}}+{\displaystyle \frac{\partial {\tau }_{rz}}{\partial z}}+{\displaystyle \frac{{\sigma }_r-{\sigma }_{\theta }}{r}}=\mathrm{0},\\ {\displaystyle \frac{\partial {\tau }_{rz}}{\partial r}}+{\displaystyle \frac{\partial {\sigma }_z}{\partial z}}+{\displaystyle \frac{{\tau }_{rz}}{r}}=\mathrm{0},\end{array},

(5)\begin{array}{l}\left(c_{\mathrm{11}}\left({\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial r^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{1}}{r^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial }{\partial r}}\right)+c_{\mathrm{44}}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial z^{\mathrm{2}}}}\right)u_r+(c_{\mathrm{13}}+c_{\mathrm{44}}){\displaystyle \frac{{\partial }^{\mathrm{2}}{u}_z}{\partial r\partial z}}=\mathrm{0},\\ (c_{\mathrm{13}}+c_{\mathrm{44}}){\displaystyle \frac{\partial }{\partial z}}\left({\displaystyle \frac{\partial }{\partial r}}+{\displaystyle \frac{\mathrm{1}}{r}}\right)u_r+\left(c_{\mathrm{44}}\left({\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial r^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial }{\partial r}}\right)+c_{\mathrm{33}}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial z^{\mathrm{2}}}}\right)u_z=\mathrm{0.}\end{array}.

Analytical solution for a single layer

To decrease the partial differential equations of normal differential equations, the spatial Hankel transform were adopted, as shown in Eq. (6). x is the Hankel transform parameters, while J_m(x_r) denotes the m-order Bessel function of the first kind. The Eq. (7) can be obtained by taking the Hankel transform of the variables of Eq. (5).

(6)\overline{f}(\xi )={\displaystyle {\int }_{\mathrm{0}}^{\infty }f(x)r{J}_m(\xi r)}\mathrm{d}r,

(7)\begin{array}{l}{\displaystyle \frac{c_{\mathrm{13}}+c_{\mathrm{44}}}{c_{\mathrm{44}}}}{\displaystyle \frac{\partial {\overline{u}}_r(\xi ,z)}{\partial z}}+\left(-{\xi }^{\mathrm{2}}+{\displaystyle \frac{c_{\mathrm{33}}}{c_{\mathrm{44}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial z^{\mathrm{2}}}}\right){\overline{u}}_z(\xi ,z)=\mathrm{0},\\ {\displaystyle \frac{c_{\mathrm{13}}+c_{\mathrm{44}}}{c_{\mathrm{44}}}}{\displaystyle \frac{\partial {\overline{u}}_r(\xi ,z)}{\partial z}}+\left(-{\xi }^{\mathrm{2}}+{\displaystyle \frac{c_{\mathrm{33}}}{c_{\mathrm{44}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial z^{\mathrm{2}}}}\right){\overline{u}}_z(\xi ,z)=\mathrm{0},\end{array},

where {\bar{u}}_z and {\bar{u}}_r are the Hankel transforms with respect to the displacement u_z and u_r. It is noteworthy that a linear hysteretic material damping was used with a damping constant value (0.001). ξ is the parameter of Hankel transform corresponding to Ref. [⁴¹]. In addition, Eq. (7) can be rewritten as Eq. (8), just one formula.

(8)\left({\displaystyle \frac{{\partial }^{\mathrm{4}}}{\partial z^{\mathrm{4}}}}-{\alpha }^{\mathrm{2}}{\xi }^{\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial z^{\mathrm{2}}}}+{\beta }^{\mathrm{4}}{\xi }^{\mathrm{4}}\right){\bar{u}}_z(\xi ,z)=\mathrm{0.}

To solve Eq. (8), we define the solution of {\bar{u}}_z shown in Eq. (9), therefore, the Eq. (8) could be rewritten as Eq. (10) by introducing the Eq. (9), where

b^{\uparrow }\mathrm{2}=I(c{{\displaystyle I}}_{\downarrow }\mathrm{11}{{\displaystyle c}}_{\downarrow }\mathrm{33}\_{{\displaystyle c}}_{\downarrow }\mathrm{13}{{\displaystyle c}}_{\downarrow }\mathrm{44}I{{\displaystyle c}}_{\downarrow }\mathrm{1}{\mathrm{3}}^{\uparrow }\mathrm{2})/{{\displaystyle c}}_{\downarrow }\mathrm{33}{{\displaystyle c}}_{\downarrow }\mathrm{44}\mathrm{and}{c}^{\uparrow }\mathrm{4}={\displaystyle \frac{{{\displaystyle c}}_{\mathrm{11}}}{{{\displaystyle c}}_{\mathrm{33}}}}.

(9){\bar{u}}_z(\xi ,z)=e^{\lambda z},

(10){\lambda }^{\mathrm{4}}-b^{\mathrm{2}}{\xi }^{\mathrm{2}}{\lambda }^{\mathrm{2}}+c^{\mathrm{4}}{\xi }^{\mathrm{4}}=\mathrm{0.}

In the derivation process of this study, the variables are assumed to be in complex field. Thus, according to the theory of ordinary differential equations, there are only two situations with different b² and c⁴ values.

1) When b^{\mathrm{4}}-\mathrm{4}{c}^{\mathrm{4}}=\mathrm{0}, \lambda =\pm b\xi /\sqrt{\mathrm{2}}, then {\bar{u}}_z(\xi ,z) can be written as:

(11){\bar{u}}_z(\xi ,z)=[A(\xi )+zB(\xi )]e^{-\lambda z}+[C(\xi )+zD(\xi )]e^{\lambda z},

where A(ξ), B(ξ), C(ξ), and D(ξ) are four coefficient vectors. To avoid the numerical overflow due to the positive exponentials in Eq. (11), e^{\lambda z} is represented as e^{-\lambda (h-z)} in the present coordinate system, the coefficient vectors is expressed as Eq. (12) for every layer of the proposed structure. Thus, {\bar{u}}_z(\xi ,z) can be rewritten in Eq. (13).

(12){A^{\prime }}_{}(\xi )=A(\xi ),{B^{\prime }}_{}(\xi )=B(\xi ),{C^{\prime }}_{}(\xi )=C(\xi )e^{\lambda h},{D^{\prime }}_{}(\xi )=D(\xi )e^{\lambda h},

(13){\bar{u}}_z(\xi ,z)=[A^{\prime }(\xi )+z{B}^{\prime }(\xi )]e^{-\lambda z}+[C^{\prime }(\xi )+z{D}^{\prime }(\xi )]e^{-\lambda (h-z)}.

Substituting Eq. (13) into Eqs. (1) and (3), the formula of {\bar{u}}_r(\xi ,z), {\bar{\tau }}_{rz}(\xi ,z), and {\bar{\sigma }}_z(\xi ,z) can be expressed as Eq. (14a), where T⁽¹⁾ as presented in Eq. (14b), the elements of the vector T⁽¹⁾ as shown in Appendix (A-1), e andbare shown in Eq. (14c).

(14a){\left[\begin{array}{llll}{\overline{u}}_z(\xi ,z)& {\overline{u}}_r(\xi ,z)& {\overline{\tau }}_{rz}(\xi ,z)& {\overline{\sigma }}_z(\xi ,z)\end{array}\right]}^{\mathrm{T}}={\mathit{T}}_{}^{(\mathrm{1})}\mathit{e}\mathit{\beta },

(14b){\mathit{T}}_{}^{(\mathrm{1})}=\left[\begin{array}{llll}{t}_{\mathrm{11}}^{(\mathrm{1})}& t_{\mathrm{12}}^{(\mathrm{1})}& t_{\mathrm{13}}^{(\mathrm{1})}& t_{\mathrm{14}}^{(\mathrm{1})}\\ t_{\mathrm{21}}^{(\mathrm{1})}& t_{\mathrm{22}}^{(\mathrm{1})}& t_{\mathrm{23}}^{(\mathrm{1})}& t_{\mathrm{24}}^{(\mathrm{1})}\\ t_{\mathrm{31}}^{(\mathrm{1})}& t_{\mathrm{32}}^{(\mathrm{1})}& t_{\mathrm{33}}^{(\mathrm{1})}& t_{\mathrm{34}}^{(\mathrm{1})}\\ t_{\mathrm{41}}^{(\mathrm{1})}& t_{\mathrm{42}}^{(\mathrm{1})}& t_{\mathrm{43}}^{(\mathrm{1})}& t_{\mathrm{44}}^{(\mathrm{1})}\end{array}\right],

(14c)\begin{array}{l}\mathit{\beta }={\left[\begin{array}{llll}{A}^{\prime }(\xi )& B^{\prime }(\xi )& C^{\prime }(\xi )& D^{\prime }(\xi )\end{array}\right]}^{\mathrm{T}}\hfill \\ diag(\mathit{e})={\left[\begin{array}{llll}{e}^{-\lambda z}& e^{-\lambda z}& e^{-\lambda (h-z)}& e^{-\lambda (h_i-z)}\end{array}\right]}^{\mathrm{T}}\hfill \end{array}.

2) When b^{\mathrm{4}}-\mathrm{4}{c}^{\mathrm{4}}\ne \mathrm{0}, \lambda ={\lambda }_{\mathrm{1},\mathrm{2}}=\pm \xi \sqrt{b^{\mathrm{2}}+\sqrt{b^{\mathrm{4}}-\mathrm{4}{c}^{\mathrm{4}}}/\mathrm{2}}, then {\bar{u}}_z(\xi ,z) yields:

(15){\overline{u}}_z(\xi ,z)=A\text{'}(\xi ){\mathrm{e}}^{-\lambda z}+B\text{'}(\xi ){\mathrm{e}}^{-\lambda z}+C\text{'}(\xi ){\mathrm{e}}^{-\lambda (h-z)}+D\text{'}(\xi ){\mathrm{e}}^{-\lambda (h-z)}.

Similarly, in this case, the formula of {\bar{u}}_r(\xi ,z), {\bar{\tau }}_{rz}(\xi ,z), and {\bar{\sigma }}_z(\xi ,z) can be expressed as Eq. (16a), where T⁽²⁾ as presented in Eq. (16b), and the elements of the vector T⁽²⁾ as shown in Appendix (A-2).

(16a){\left[\begin{array}{llll}{\overline{u}}_z(\xi ,z)& {\overline{u}}_r(\xi ,z)& {\overline{\tau }}_{rz}(\xi ,z)& {\overline{\sigma }}_z(\xi ,z)\end{array}\right]}^{\mathrm{T}}={\mathit{T}}_{}^{(\mathrm{2})}\mathit{e}\mathit{\beta }\mathfrak{\beta },

(16b){\mathit{T}}_{}^{(\mathrm{2})}=\left[\begin{array}{llll}{t}_{\mathrm{11}}^{(\mathrm{2})}& t_{\mathrm{12}}^{(\mathrm{2})}& t_{\mathrm{13}}^{(\mathrm{2})}& t_{\mathrm{14}}^{(\mathrm{2})}\\ t_{\mathrm{21}}^{(\mathrm{2})}& t_{\mathrm{22}}^{(\mathrm{2})}& t_{\mathrm{23}}^{(\mathrm{2})}& t_{\mathrm{24}}^{(\mathrm{2})}\\ t_{\mathrm{31}}^{(\mathrm{2})}& t_{\mathrm{32}}^{(\mathrm{2})}& t_{\mathrm{33}}^{(\mathrm{2})}& t_{\mathrm{34}}^{(\mathrm{2})}\\ t_{\mathrm{41}}^{(\mathrm{2})}& t_{\mathrm{42}}^{(\mathrm{2})}& t_{\mathrm{43}}^{(\mathrm{2})}& t_{\mathrm{44}}^{(\mathrm{2})}\end{array}\right].

Analytical solution for the multi-layered structure

The boundary condition is a place on the multi-layered structure where either the external force or the displacement need be defined at the beginning of the mathematical derivation for the analytical solution. In this study, the boundary conditions on the top-surface, i.e., free-surface, and at the bottom-side, i.e., infinity, of the multi-layered structure yield six equations. The boundary condition at the top-surface read as the Eq. (17).

(17)\bar{\sigma }(r,\mathrm{0})e_z=\bar{q},

(18){C^{\prime }}_n(\xi )={D^{\prime }}_n(\xi )=\mathrm{0},

where \bar{q}is the imposed external force vector. In the bottom-side of the structure, the radial condition is considered as infinity, and the displacement field scatter to infinity. In the displacement equations, i.e., Eqs. (14a) and (16a), the radial condition implies that the vector is equal to zero at the bottom-side of the multi-layered structure.

In addition, the interlayer condition, sliding and bonded, between two adjacent layers is considered as a key requirement to make up the global vector matrix in this study. Figure 1 illustrates the interlayer parameters between layer-i and layer-(i + 1), where z_iand r_iare the vertical coordinate and horizontal coordinate in the coordinate system, respectively. The bonded interlayer condition, i.e., continuity, is interlayer condition is commonly employed and read that the traction vectors, i.e., stresses, from both sides of an interlayer are equal. The bonded interlayer condition between the layer-i and layer-(i + 1) is demonstrated in Eq. (19).

(19){\left[\begin{array}{l}{u}_z(r,z_i)\\ u_r(r,z_i)\\ {\tau }_{rz}(r,z_i)\\ {\sigma }_z(r,z_i)\end{array}\right]}_i={\left[\begin{array}{l}{u}_z(r,z_i)\\ u_r(r,z_i)\\ {\tau }_{rz}(r,z_i)\\ {\sigma }_z(r,z_i)\end{array}\right]}_{i+\mathrm{1}},

(20){\left[\begin{array}{l}{u}_z(r,z_i)\\ {\tau }_{rz}(r,z_i)\\ \mathrm{0}\\ {\sigma }_z(r,z_i)\end{array}\right]}_i={\left[\begin{array}{l}{u}_z(r,z_i)\\ \mathrm{0}\\ {\tau }_{rz}(r,z_i)\\ {\sigma }_z(r,z_i)\end{array}\right]}_{i+\mathrm{1}}.

Furthermore, the sliding interlayer condition stipulates that the shear stresses between the adjacent layers are equal to zero at both sides of the interlayer, while the vertical displacements and the third component of the traction vector remain the same as Eq. (19). Therefore, the sliding case for the interlayer between the layer-i and layer-(i + 1) is illustrated in Eq. (20). The boundary condition at the top-surface and the bottom-side of the multi-layered structure plus the interlayer conditions lead to solve the unknowns [⁵⁰].

(21)({\mathit{M}}_i{\mathit{m}}_i){\mathit{\beta }}_i=({\mathit{N}}_{i+\mathrm{1}}{\mathit{n}}_{i+\mathrm{1}}){\mathit{\beta }}_{i+\mathrm{1}}.

The global matrix technique is applied to compute the mechanical behavior of multi-layered structure in this study, and this technique eases the manipulation of the combination of sliding and bonded interlayer conditions [⁵⁵]. The assembly of the global matrix relies on Eq. (21) that implies quantities of the interlayer between the layer-i and layer-(i + 1) in the multi-layered structure.

(22a)\begin{array}{ll}{\mathit{M}}_i^{}=\left[\begin{array}{llll}{t}_{\mathrm{11}}^{}& t_{\mathrm{12}}^{}& t_{\mathrm{13}}^{}& t_{\mathrm{14}}^{}\\ t_{\mathrm{21}}^{}& t_{\mathrm{22}}^{}& t_{\mathrm{23}}^{}& t_{\mathrm{24}}^{}\\ t_{\mathrm{31}}^{}& t_{\mathrm{32}}^{}& t_{\mathrm{33}}^{}& t_{\mathrm{34}}^{}\\ t_{\mathrm{41}}^{}& t_{\mathrm{42}}^{}& t_{\mathrm{43}}^{}& t_{\mathrm{44}}^{}\end{array}\right],& {\mathit{N}}_i^{}=\left[\begin{array}{llll}{t}_{\mathrm{11}}^{}& t_{\mathrm{12}}^{}& t_{\mathrm{13}}^{}& t_{\mathrm{14}}^{}\\ t_{\mathrm{21}}^{}& t_{\mathrm{22}}^{}& t_{\mathrm{23}}^{}& t_{\mathrm{24}}^{}\\ t_{\mathrm{31}}^{}& t_{\mathrm{32}}^{}& t_{\mathrm{33}}^{}& t_{\mathrm{34}}^{}\\ t_{\mathrm{41}}^{}& t_{\mathrm{42}}^{}& t_{\mathrm{43}}^{}& t_{\mathrm{44}}^{}\end{array}\right]\end{array},

(22b)\begin{array}{ll}{\mathit{M}}_i^{}=\left[\begin{array}{llll}{t}_{\mathrm{11}}^{}& t_{\mathrm{12}}^{}& t_{\mathrm{13}}^{}& t_{\mathrm{14}}^{}\\ t_{\mathrm{31}}^{}& t_{\mathrm{32}}^{}& t_{\mathrm{33}}^{}& t_{\mathrm{34}}^{}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ t_{\mathrm{41}}^{}& t_{\mathrm{42}}^{}& t_{\mathrm{43}}^{}& t_{\mathrm{44}}^{}\end{array}\right],& {\mathit{N}}_i^{}=\left[\begin{array}{llll}{t}_{\mathrm{11}}^{}& t_{\mathrm{12}}^{}& t_{\mathrm{13}}^{}& t_{\mathrm{14}}^{}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ t_{\mathrm{31}}^{}& t_{\mathrm{32}}^{}& t_{\mathrm{33}}^{}& t_{\mathrm{34}}^{}\\ t_{\mathrm{41}}^{}& t_{\mathrm{42}}^{}& t_{\mathrm{43}}^{}& t_{\mathrm{44}}^{}\end{array}\right]\end{array}\mathrm{.}

The matrices of M_i and N_i are different from bonded interlayer to sliding interlayer. The definition of matrices M_i and N_i for bonded interlayer are obtained by inserting Eq. (14a) (or Eq. (16a)) into Eq. (19), and substitution of Eq. (14a) (or Eq. (16a)) into Eq. (20) induces to the determination of matrices M_i and N_i for the sliding interlayer case. Both the cases of b⁴ = 4c⁴ and b⁴≠4c⁴ are applied to these interlayer conditions, where the matrices M_i and N_i in Eq. (22a) with respect to the bonded interlayer case, and the matrices M_i and N_i for the sliding interlayer case as presented in Eq. (22b). Note that the matrix element t_ij can be substituted by t_ij⁽¹⁾ and t_ij⁽²⁾ for the case of b⁴ = 4c⁴ and b⁴≠4c⁴, respectively. However, the exponential matrices m_i and n_i that associated with e_i, and it is a diagonal matrix, which is not related to the interlayer conditions. The expressions of m_i and n_i are presented in Eqs. (23) and (24) for the cases of b⁴ = 4c⁴ and b⁴≠4c⁴, respectively. Furthermore, because of every layer of the multi-layered structure uses an independent local coordinate system, the depth z takes h_i and zero for the lower interlayer side of layer-i and the upper interlayer side of layer-(i + 1), respectively, in the process of substituting these parameters into Eq. (14a) (or Eq. (16a)).

(23)\begin{array}{l}diag({\mathit{m}}_i^{(\mathrm{1})})={\left[\begin{array}{llll}{e}^{-\lambda h_i}& e^{-\lambda h_i}& \mathrm{1}& \mathrm{1}\end{array}\right]}^{\mathrm{T}}\\ diag({\mathit{n}}_i^{(\mathrm{1})})={\left[\begin{array}{llll}\mathrm{1}& \mathrm{1}& e^{-\lambda h_i}& e^{-\lambda h_i}\end{array}\right]}^{\mathrm{T}}\end{array},

(24)\begin{array}{l}diag({\mathit{m}}_i^{(\mathrm{2})})={\left[\begin{array}{llll}{e}^{-{\lambda }_{\mathrm{1}}{h}_i}& e^{-{\lambda }_{\mathrm{2}}{h}_i}& \mathrm{1}& \mathrm{1}\end{array}\right]}^{\mathrm{T}}\\ diag({\mathit{n}}_i^{(\mathrm{2})})={\left[\begin{array}{llll}\mathrm{1}& \mathrm{1}& e^{-{\lambda }_{\mathrm{1}}{h}_i}& e^{-{\lambda }_{\mathrm{2}}{h}_i}\end{array}\right]}^{\mathrm{T}}\end{array}.

To obtain the mechanical behavior, i.e., displacements and stresses, of the multi-layered structure, the linear system, as shown in Eq. (25), needs be solved. The filling procedure of the global matrix, the first block of linear system, based on the repeating of the Eq. (21). Equation (25) verified any one of the interlayer cases, i.e., sliding and bonded.

(25)\left[\begin{array}{llllllll}{\mathit{D}}_{\mathrm{0}}{\mathit{e}}_{\mathrm{1}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ {\mathit{{\rm M}}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}& -{\mathit{N}}_{\mathrm{2}}{\mathit{n}}_{\mathrm{2}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& {\mathit{M}}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}& -{\mathit{N}}_{\mathrm{3}}{\mathit{n}}_{\mathrm{3}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{...}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& {\mathit{M}}_i{\mathit{m}}_i& -{\mathit{N}}_{i+\mathrm{1}}{\mathit{n}}_{i+\mathrm{1}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{...}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& {\mathit{M}}_{n-\mathrm{1}}{\mathit{m}}_{n-\mathrm{1}}& -{\mathit{N}}_{n-\mathrm{1}}{\mathit{n}}_{n-\mathrm{1}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& {\mathit{M}}_{n-\mathrm{1}}{\mathit{m}}_{n-\mathrm{1}}& {\mathit{D}}_n{\mathit{m}}_n\end{array}\right]\left[\begin{array}{l}{\mathit{\beta }}_{\mathrm{1}}\\ {\mathit{\beta }}_{\mathrm{2}}\\ ⋮\\ ⋮\\ {\mathit{\beta }}_i\\ ⋮\\ {\mathit{\beta }}_{n-\mathrm{1}}\\ {\mathit{\beta }}_n^{\text{'}}\end{array}\right]=\left[\begin{array}{l}\mathit{f}\\ \mathrm{0}\\ ⋮\\ ⋮\\ ⋮\\ ⋮\\ ⋮\\ \mathrm{0}\end{array}\right],

where the first block of the global matrix D₀ reflects the boundary condition at the top-surface of the multi-layered structure. Meanwhile, the boundary condition at the bottom-side of the multi-layered structure affects the last block of the global matrix D_n [⁵¹,⁵⁶]. Equation (26) depicts the detailed boundary condition at the top-surface and the bottom-side for the bonded interlayer condition, while Eq. (27) illustrates those information for the sliding interlayer condition. It should be noted that both the cases of b⁴ = 4c⁴ and b⁴≠4c⁴ are also applied to these interlayer conditions for the determination of matrices of D₀ and D_n, which indicates that the matrix element t_ij can be replaced by t_ij⁽¹⁾ and t_ij⁽²⁾ for the case of b⁴ = 4c⁴ and b⁴≠4c⁴, respectively. The zero vector in the above-mentioned linear system represents a 4th order vector. The number of unknown coefficient vector is equal to the number of coefficient equations in the linear system. Therefore, the unknown coefficient vectors b_i and b_n’ will be obtained by solving Eq. (25).

(26)\begin{array}{ll}{\mathit{D}}_{\mathrm{0}}^{}=\left[\begin{array}{llll}{t}_{\mathrm{31}}^{}\hfill & t_{\mathrm{32}}^{}\hfill & t_{\mathrm{33}}^{}\hfill & t_{\mathrm{34}}^{}\hfill \\ t_{\mathrm{41}}^{}\hfill & t_{\mathrm{42}}^{}\hfill & t_{\mathrm{43}}^{}\hfill & t_{\mathrm{44}}^{}\hfill \end{array}\right]& ,{\mathit{D}}_n^{}={\left[\begin{array}{llll}{t}_{\mathrm{11}}^{}\hfill & t_{21}^{}& t_{\mathrm{31}}^{}\hfill & t_{\mathrm{41}}^{}\hfill \\ t_{\mathrm{12}}^{}\hfill & t_{22}^{}& t_{\mathrm{32}}^{}\hfill & t_{\mathrm{42}}^{}\hfill \end{array}\right]}_{}^{\mathrm{T}}\end{array}.

(27)\begin{array}{ll}{\mathit{D}}_{\mathrm{0}}^{}=\left[\begin{array}{llll}{t}_{\mathrm{31}}^{}\hfill & t_{\mathrm{32}}^{}\hfill & t_{\mathrm{33}}^{}\hfill & t_{\mathrm{34}}^{}\hfill \\ t_{\mathrm{41}}^{}\hfill & t_{\mathrm{42}}^{}\hfill & t_{\mathrm{43}}^{}\hfill & t_{\mathrm{44}}^{}\hfill \end{array}\right]& ,{\mathit{D}}_n^{}={\left[\begin{array}{llll}{t}_{\mathrm{11}}^{}\hfill & \mathrm{0}\hfill & t_{\mathrm{31}}^{}\hfill & t_{\mathrm{41}}^{}\hfill \\ t_{\mathrm{12}}^{}\hfill & \mathrm{0}\hfill & t_{\mathrm{32}}^{}\hfill & t_{\mathrm{42}}^{}\hfill \end{array}\right]}_{}^{\mathrm{T}}\end{array}.

Solving the linear system Eq. (25), and then substituting the solved coefficient vectors β_i and β_n’ into Eq. (14a) (or Eq. (16a)) generates the components of displacement and stress in each layer of the multi-layered structure in Hankel domain. Then, the numerical inversion of integral transform, i.e. inverse Hankel transform, was applied to solve the components of displacement and stress in real domain. Gauss-Legendre quadrature is a key scheme in the numerical inversion process of this study. The Gauss-Legendre quadrature formula is presented in Eq. (28), while the infinite integral formula in the right side of Eq. (28) can be split into a series of finite integral formula sums, as shown in Eq. (29), where ξ is the parameter of Hankel transform corresponding to the radius (r).

(28)f(r)=K_m^{-\mathrm{1}}[f(\xi )]={\displaystyle {\int }_{\mathrm{0}}^{\infty }\overline{f}(\xi )\xi J_m(\xi r)\mathrm{d}\xi },

(29)f(r)={\displaystyle {\int }_{\mathrm{0}}^{b_{\mathrm{1}}}\overline{f}(\xi )\xi J_m(\xi r)\mathrm{d}\xi }+{\displaystyle {\int }_{b_{\mathrm{1}}}^{b_{\mathrm{2}}}\overline{f}(\xi )\xi J_m(\xi r)\mathrm{d}\xi }+\cdots +{\displaystyle {\int }_{b_n}^{b_{n+\mathrm{1}}}\overline{f}(\xi )\xi J_m(\xi r)\mathrm{d}\xi }+\cdots .

The finite integral item on the right side of Eq. (29) can be approximately computed by using the numerical methods. The six-point Gauss-Legendre quadrature was employed in the calculation for the approximate items of Eq. (29) in this study, and the approximate calculation formula is as follows:

(30){\displaystyle {\int }_{b_n}^{b_{n+\mathrm{1}}}\overline{f}(\xi )J_m(\xi r)\xi \mathrm{d}\xi }={\displaystyle \frac{b_{n+\mathrm{1}}-b_n}{\mathrm{2}}}{\displaystyle \sum _{p=\mathrm{1}}^{\mathrm{6}}{w}_p}\overline{f}({\beta }_p)J_m(r{\beta }_p){\beta }_p,

(31){\beta }_p=\left({\displaystyle \frac{b_{n+\mathrm{1}}-b_n}{\mathrm{2}}}\right)x_p+\left({\displaystyle \frac{b_{n+\mathrm{1}}+b_n}{\mathrm{2}}}\right),

where {{\displaystyle \beta }}_p is shown in Eq. (31). x_p and w_p are the Gaussian node and Gaussian coefficient, respectively. b_n is the limits of the integration, which is selected arbitrarily. In addition, Cornille [⁵⁷] reported that the convergence can be accelerated if the continuous zero points of the first type Bessel-Function J_m(x) is used as the integration interval. After conducted a series of sensitivity analysis, Cornille [⁵⁷] also concluded that the interval formed by the continuous zero points (zero-interval, divided into ten equal-spaced cells) can make the numerical integration result with acceptable accuracy. In Eq. (30), the upper-integral limit of the left side formula is infinite, which means that a series of integral sums on the right-side formulas should also span an infinite range. However, Kim’s research [⁵⁸] shows that the infinite integration converges rapidly after numerical quadrature in the first five zero-intervals of the J_m(x). Therefore, in this study, the first five zero-intervals of J_m(x) are used in the case of integral summation near the loading center (r=0), while fewer zero-intervals are used in the area far from the loading center.

Model validations

The validation of the present analytical solution was performed by comparison with the Refs. [³⁷,⁵⁹]. and the software BISAR. The following three examples in Refs. [³⁷,⁵⁹]. and the software BISAR were used to verify the accuracy of the present analytical solutions.

First, Ref. [⁵⁹] and present analytical solution computed the vertical displacement on the surface of the structure at different transversely isotropic conditions (using same structure parameters). The multi-layered structure specimen is based on a multi-layered soil structure subjected to axisymmetric loading as used for that model validation. In Ref. [⁵⁹], the soil structure was transversely isotropic linear elastic with E_v = 40 MPa, μ_vh = μ_h = 0.3, G_v = 30 MPa, and the modulus ratio n = E_h/ E_v is 1.0 and 0.5. The axisymmetric load was idealized as a single vertical circular uniform load q = 0.71 MPa with a diameter of 0.213 m. Figure 2 illustrates the vertical displacements comparison between the computed results obtained from present analytical solution and Ref. [⁵⁹]. when the transversely isotropic properties of soil structure were considered. It can be seen that the curves obtained from present analytical solution, with different transversely isotropic conditions, were a good match with the ones stemming from Ref. [⁵⁹].

Secondly, BISAR verified the present analytical solution and calculated the radial stresses at different depths of double-layered structure. The interlayer conditions between the top and bottom layer, sliding and bonded, were also considered in the comparison. In this case of validation, the multi-layered structure sample is based on an isotropic linear elastic double flexible pavement structure subjected to axisymmetric loading. BISAR is a stable/popular program used in the mechanical behavior analysis of multi-layered structure, but the program cannot be suited for the anisotropic structures. Thus, in order to verify the accuracy of the present analytical solutions by using the BISAR, the anisotropic problem was degenerated into an isotropic problem. The present analytical solution in physical domain can be obtained by assuming that the Young’s modules are the same in both horizontal and vertical directions. The axisymmetric load was idealized as a single vertical circular uniform load q= 0.71 MPa, and 0.15 m diameter. The asphalt mixture layer (top-layer) and the structure layer were isotropic linear elastic with E = 100 MPa, μ = 0.25, ρ = 2000 kg/m³, h = 0.2 m; E = 50 MPa, μ = 0.25, ρ = 1800 kg/m³, h = infinite, respectively. Figure 3 demonstrated the comparison between the calculated results obtained from the present analytical solution and software BISAR when the interlayer condition between these two layers was considered. The comparison results depict that the curves from present analytical solutions are good match with the curves obtained from BISAR.

Finally, the calculated response variables based on the present analytical solution in this study were compared with those obtained by Tutumluer et al. [³⁷], which involves predicted and measured response variables for conventional linear elastic pavements. The predicted response variables were obtained by using GT-PAVE finite element program, while the measure response variables are collected from a full-scale pavement test study conducted at Georgia Tech. Table 1 summarizes the material properties and the model parameters used in this study according to the prediction and measurement by Tutumluer et al. [³⁷]. The layer-1 and layer-3 were treated as isotropic linear elastic materials, while the layer-2 was given transversely isotropic material properties. For the anisotropic layer-2, a 0.15 horizontal modulus reduction factor (E_h/E_v) and 0.33 shear modulus reduction factor (G_v/E_v) are used in the GT-PAVE finite element program. The axisymmetric load was idealized as a single vertical circular uniform load q= 0.69 MPa (100 psi) with a diameter of 0.116 m (4.55 in.). Table 2 demonstrates the response variable comparison between the computed results obtained from present analytical solution and Ref. [³⁷]. When the transversely isotropic properties of layer-2 were considered. In contrast to the GT-PAVE finite element program results, the response variables obtained from the present analytical solution are in reasonably good agreement with the observed behavior of measurement in field. This is especially true for the calculated values of vertical displacement ({{\displaystyle \delta }}_z) on the top layer-1 and vertical strain ({{\displaystyle \epsilon }}_z )/ stress (\sigma \downarrow Z) on the top layer-3. Therefore, the response variables from present analytical solution are good match with the field measurement data, and this proximity is slightly better than GT-PAVE finite element program. In summary, the analytical solution presented in this study is accurate to be used in the analysis of the influence of anisotropic layer/interlayer conditions on multi-layered structure.

Numerical results and discussions

The multi-layered structure model under consideration is a composite structure, i.e., three-layered flexible pavement, which consists of three layers that are defined in Table 3. The transversely isotropic properties can be introduced to layer-1 and layer-2, i.e., the asphalt mixture layer and the macadam base layer in the flexible pavement [³⁰,⁶⁰,⁶¹], while the layer-3, soil-foundation, is considered as a linear elastic body. The elastic modulus ratio E_h/E_v, i.e., the ratio between the modulus in horizontal direction and vertical direction, of the layer-1 and layer-2 were defined as n₁ and n₂, respectively. The elastic modulus ratio n_i (i = 1 and 2) was selected as 1.0 (isotropy), 0.75, 0.50, 0.25. In addition, the axisymmetric load was idealized as a single vertical circular uniform load q= 0.71MPa, and 0.15m diameter. The investigation of the effects of interlayer condition and transversely isotropic properties on multi-layered structure as shown in below two subsections.

Impact of interlayer condition on multi-layered flexible pavement

The structure layers of multi-layered flexible pavement are initially full bonded initially, but it will eventually fail due to the differential expansion and vehicle loads between the asphalt mixture layer (layer-1) and the macadam base layer (layer-2). This section mainly investigates the impact of the interlayer condition between the layer-1 and layer-2. Two cases of interlayer conditions were studied, i.e. sliding and bonded. In addition, the computed models were divided into several combination cases to consider, and the specific combination cases are present in Table 4.

The vertical displacement of the top-surface of the flexible pavement structure is critical to reflect the deformation resistance of the flexible pavement construction. Figure 4 shows the simulation results for the vertical displacements of the top-surface distributed in the radial direction. It can be seen that the vertical displacement decreased with the increase in radial distance from loading center, and the magnitude of the vertical displacement appeared in the loading center. The vertical displacement with sliding condition was greater than those of the bonded condition at the same transversely isotropic condition, i.e., K3>K1 and K4>K2. The magnitude value of the vertical displacement of K3 is larger than those of K1 by 37.0%, and the magnitude value of the vertical displacement of K4 is larger than those of K2 by 37.9%. Furthermore, in the same interlayer condition, the vertical displacement with the transversely isotropic layer-1 and layer-2 was greater than those with the isotropic layer-1 and layer-2, i.e., K2>K1 and K4>K3. The magnitude value of the vertical displacement of K2 is larger than those of K1 by 14.0%, and the magnitude value of the vertical displacement of K4 is larger than those of K3 by 14.2%. Therefore, the sliding interlayer between layer-1 and layer-2 and the transversely isotropic (n<1.0) layer-1 and layer-2 lead to a higher vertical displacement, and the influence of the former is more obvious.

The fatigue cracking in the flexible pavement is usually caused by the radial stress in the pavement structure being too large, and the horizontal shear stress is a key factor to cause the rutting or shear-deformation of flexible pavement. The simulation results of the radial stress and horizontal shear stress are illustrated in Figs. 5 and 6, respectively. Figure 5 depicts the radial stress as a function of structure depth. These vertical profiles were collected under the center of loading (r = 0 m). The radial stresses for the sliding interlayer cases were higher than for the bonded interlayer cases, i.e., K3>K1 and K4>K2. As expected, the sliding interlayer case, i.e., K3 and K4, enforced between layer-1 and layer-2 yields a discontinuity in the radial stress at the position of z= 0.18m. This discontinuity was noticeable in all curves of Fig. 5, the bonded interlayer cases, i.e., K1 and K2, enforced between layer-1 and layer-2 also yielded a discontinuity in the radial stresses at z= 0.18 m, but this gap was obviously smaller than that of the sliding interlayer cases. Note that the discontinuity and peak value was greater when considering the transverse isotropy of layer-1 and layer-2, i.e., K2>K1 and K4>K3.

In addition, Fig. 6 presents the relationship between the horizontal shear stress and the structure depth. These vertical profiles were also collected under the center of loading (r = 0 m). The horizontal shear stresses increased with the increase of structure depth and the peak values occur at the position of z = 0.12m, and then the horizontal shear stresses decreased with the increase of structure depth. The horizontal shear stress for the sliding interlayer cases were higher than for the bonded interlayer cases, i.e., K3>K1 and K4>K2. The peak value of horizontal shear stress for the sliding cases is higher than those for the bonded cases by 11.4% (K3>K1) and 11.6% (K4>K2). With the structure depth increasing, as expected, the sliding interlayer case, i.e., K3 and K4, enforced between layer-1 and layer-2 causes a discontinuity in the horizontal shear stresses at the position of z = 0.18 m. The horizontal shear stresses are zero at the interface of layer-1 and layer-2. This discontinuity was also noticeable for the bonded cases, i.e., K1 and K2, but this gap was obviously smaller than that of the sliding interlayer cases. Note that the discontinuity and peak value was greater when considering the transverse isotropy of layer-1 and layer-2. i.e., K2>K1 and K4>K3. Therefore, the sliding interlayer and the transversely isotropic layer-1 and layer-2 also lead to higher radial stress and horizontal shear stress, with the influence of the former is more obvious.

Impact of transversely isotropic properties on multi-layered flexible pavement

To investigate the impact of modulus ratio n = E_h/ E_v of layer-1 and layer-2 on the multi-layered flexible pavement, the simulation was carried out by using the same structural parameters used in last section with one exception. The modulus ratio were selected as 0.25, 0.50, 0.75, and 1.00. The simulation results are demonstrated in Figs. 7–9.

Figures 7(a) and 7(b) illustrate the relationships between the maximum displacements on the top-surface of layer-1 and the modulus ratio, for layer-1 and layer-2, respectively. Note that the interlayer condition between layer-1 and layer-2 effects are also to be analyzed here. The maximum displacement exhibits a decreasing trend with the increasing modulus ratio, which include transversely isotropic (n<1.0) of layer-1 and layer-2. Figure 7(a) depicted that the maximum displacements with the sliding interlayer were higher than those with a bonded interlayer, and the transversely isotropic properties (n<1.0) of layer-2 lead to higher displacement in the same modulus ratio of layer-1. Likewise, in the Fig. 7(b), the transversely isotropic properties (n<1.0) of layer-1 also lead to a higher displacement with the same modulus ratio of layer-2, and the maximum displacement with sliding interlayer were higher than those with a bonded interlayer.

Figures 8(a) and 8(b) show the relationships between the maximum radial stress and the modulus ratio, for layer-1 and layer-2, respectively. As the ratio of the layer modulus increased, radial stresses decreased. It is noteworthy that the maximum radial stresses with the sliding interlayer were higher than those with bonded interlayer, and the transversely isotropic properties (n<1.0) of layer-2 (or layer-1) lead to a higher maximum displacement in the same modulus ratio of layer-1 (or layer-2). Therefore, the transversely isotropic layers have a significant influence on the maximum displacement and radial stress of multi-layered structure.

In addition, the transversely isotropic properties layers also have a significant influence on the maximum horizontal shear stress of multi-layered flexible pavement, as shown in Figs. 9(a) and 9(b). Figure 9(a) noted that the maximum horizontal shear stress increased with the increase of modulus ratio of layer-1. The maximum horizontal shear stresses with a sliding interlayer is higher than those with a bonded interlayer, and the transversely isotropic properties (n<1.0) of layer-2 cause to a higher maximum horizontal shear stress in the same modulus ratio of layer-1. However, it can be seen in Fig. 9(b) that the maximum horizontal shear stress decreased with the increase of modulus ratio of layer-1. The transversely isotropic properties (n<1.0) of layer-1 cause to a higher maximum horizontal shear stress in the same modulus ratio of layer-2, and the sliding interlayer condition leads to the higher maximum horizontal shear stress as compared to the bonded interlayer condition between the layer-1 and layer-2. Therefore, the transversely isotropic properties of layers also significantly influence the maximum horizontal shear stress of multi-layered structure. The transversely isotropic properties (n<1.0) of layer-1 lead to lower horizontal shear stress, while the transversely isotropic properties (n<1.0) of layer-2 caused higher horizontal shear stress. Furthermore, it should be noted that the influence of the change in modulus ratio of layer-2 on multi-layered structure was slightly larger than those of layer-1. In other words, for the background of multi-layered flexible pavement structure, the transversely isotropic properties of macadam base layer (layer-2) has a larger impact on the multi-layered flexible pavement structure compare with asphalt mixture layer (layer-1).

Summary and conclusions

In this study, the impact of the anisotropic layers and the interlayer conditions on the mechanical behaviors of multi-layered structure was provided by an application dedicated to multi-layered flexible pavement structures. The computing results show that both the anisotropic properties of layers and the interlayer conditions between the structural layers have significant effect on the multi-layered structure. The detailing conclusions are as follow:

1) Displacement on the top-surface of the layer-1 is higher for the sliding interlayer condition than for that of bonded interlayer condition. The obvious discontinuity in the radial stress and horizontal shear stress is observed at the interface between the layers, which means that the expired internal bond between the layers will increase the risk of structural deformation and damage.

2) Maximum displacement and radial stress decreased with an increase in modulus ratio n = E_h/ E_v of layer-1 and layer-2. The transversely isotropic properties (n<1.0) of layer-1 lead to lower horizontal shear stress, while the transversely isotropic properties (n<1.0) of layer-2 caused the higher horizontal shear stress. The influence of the change in modulus ratio of layer-2 on multi-layered structure is slightly larger than those of layer-1.

Research prospects

The present analytical solution considered the impact of bonded and sliding interlayers on the mechanical behavior of anisotropic multi-layered structure under axisymmetric loading. However, the current analytical solution lacks the thorough understanding of bonding behaviors in the multi-layered structures. The actual interlayer state may between full bonded and sliding in the structure. Therefore, if possible, there is need to develop an analytical solution to accurately investigate the change in the multi-layered structure critical responses due to the actual interlayer conditions. In addition, as the mechanical behaviors of multi-layered structure are affected by multiple uncertain parameters such as the interlayer conditions, transversely isotropic coefficients, temperatures, and loading, a comprehensive study of the effects of these parameters on the mechanical behavior is required. According to the literature review [⁸,⁶²–⁶⁶], the sensitivity analysis method could be employed to quantify the influence of the parameters on the mechanical behaviors in the further researches.

References

[1] Tafreshi S M, Khalaj O, Dawson A. Repeated loading of soil containing granulated rubber and multiple geocell layers. Geotextiles and Geomembranes, 2014, 42(1): 25–38

[2] You L, Yan K, Hu Y, Zollinger D G. Spectral element solution for transversely isotropic elastic multi-layered structures subjected to axisymmetric loading. Computers and Geotechnics, 2016, 72: 67–73

[3] Cai Y, Pan E. Surface loading over a transversely isotropic and multilayered system with imperfect interfaces: Revisit enhanced by the dual-boundary strategy. International Journal of Geomechanics, 2018, 18(6): 04018032

[4] Xu Q, Prozzi J A. A time-domain finite element method for dynamic viscoelastic solution of layered-half-space responses under loading pulses. Computers & Structures, 2015, 160: 20–39

[5] Pandolfi A, Gizzi A, Vasta M. Visco-electro-elastic models of fiber-distributed active tissues. Meccanica, 2017, 52(14): 3399–3415

[6] Ai Z Y, Cheng Y C, Zeng W Z. Analytical layer-element solution to axisymmetric consolidation of multilayered soils. Computers and Geotechnics, 2011, 38(2): 227–232

[7] Zhang Y, Luo R, Lytton R L. Microstructure-based inherent anisotropy of asphalt mixtures. Journal of Materials in Civil Engineering, 2011, 23(10): 1473–1482

[8] Hamdia K M, Silani M, Zhuang X, He P, Rabczuk T. Stochastic analysis of the fracture toughness of polymeric nanoparticle composites using polynomial chaos expansions. International Journal of Fracture, 2017, 206(2): 215–227

[9] Chen G. Steady-state solutions of multilayered and cross-anisotropic poroelastic half-space due to a point sink. International Journal of Geomechanics, 2005, 5(1): 45–57

[10] Kim I, Tutumluer E. Field validation of airport pavement granular layer rutting predictions, Transportation Research Record. Transportation Research Record: Journal of the Transportation Research Board, 2006, 1952(1): 48–57

[11] Lv S, Wang S, Liu C, Zheng J, Li Y, Peng X. Synchronous testing method for tension and compression moduli of asphalt mixture under dynamic and static loading states. Journal of Materials in Civil Engineering, 2018, 30(10): 04018268

[12] Al-Qadi I L, Wang H, Yoo P J, Dessouky S H. Dynamic analysis and in situ validation of perpetual pavement response to vehicular loading. Transportation Research Record: Journal of the Transportation Research Board, 2008, 2087(1): 29–39

[13] Talebi H, Silani M, Bordas S P, Kerfriden P, Rabczuk T. A computational library for multiscale modeling of material failure. Computational Mechanics, 2014, 53(5): 1047–1071

[14] Msekh M A, Cuong N, Zi G, Areias P, Zhuang X, Rabczuk T. Fracture properties prediction of clay/epoxy nanocomposites with interphase zones using a phase field model. Engineering Fracture Mechanics, 2018, 188: 287–299

[15] Kim D W, Choi H S, Lee C, Blumstein A, Kang Y. Investigation on methanol permeability of Nafion modified by self-assembled clay-nanocomposite multilayers. Electrochimica Acta, 2004, 50(2–3): 659–662

[16] Malwitz M M, Lin-Gibson S, Hobbie E K, Butler P D, Schmidt G. Orientation of platelets in multilayered nanocomposite polymer films. Journal of Polymer Science. Part B, Polymer Physics, 2003, 41(24): 3237–3248

[17] Ai Z Y, Cao Z. Vibration isolation of row of piles embedded in transverse isotropic multi-layered soils. Computers and Geotechnics, 2018, 99: 115–129

[18] Ai Z Y, Ye Z, Zhao Z, Wu Q L, Wang L J. Time-dependent behavior of axisymmetric thermal consolidation for multilayered transversely isotropic poroelastic material. Applied Mathematical Modelling, 2018, 61: 216–236

[19] You L, Yan K, Liu N, Shi T, Lv S. Assessing the mechanical responses for anisotropic multi-layered medium under harmonic moving load by Spectral Element Method (SEM). Applied Mathematical Modelling, 2019, 67: 22–37

[20] Liu K, Li Y, Wang F, Xie H, Pang H, Bai H. Analytical and model studies on behavior of rigid polyurethane composite aggregate under compression. Journal of Materials in Civil Engineering, 2019, 31(3): 04019007

[21] Lefeuve-Mesgouez G, Le Houédec D, Peplow A. Ground vibration in the vicinity of a high-speed moving harmonic strip load. Journal of Sound and Vibration, 2000, 231(5): 1289–1309

[22] Hung H H, Yang Y B. Elastic waves in visco-elastic half-space generated by various vehicle loads. Soil Dynamics and Earthquake Engineering, 2001, 21(1): 1–17

[23] Zhai W, Song E. Three dimensional FEM of moving coordinates for the analysis of transient vibrations due to moving loads. Computers and Geotechnics, 2010, 37(1–2): 164–174

[24] Dehghan Manshadi S, Khojasteh A, Nategh S, Rahimian M. Interaction between annular crack and rigid disc inclusion in a transversely isotropic solid. Meccanica, 2018, 53(11–12): 2973–2997

[25] Biot M A. Theory of elasticity and consolidation for a porous anisotropic solid. Journal of Applied Physics, 1955, 26(2): 182–185

[26] Biot M A. Theory of deformation of a porous viscoelastic anisotropic solid. Journal of Applied Physics, 1956, 27(5): 459–467

[27] Hematiyan M, Khosravifard A, Shiah Y. A new stable inverse method for identification of the elastic constants of a three-dimensional generally anisotropic solid. International Journal of Solids and Structures, 2017, 106–107: 240–250

[28] Eskandari M, Samea P, Ahmadi S F. Axisymmetric time-harmonic response of a surface-stiffened transversely isotropic half-space. Meccanica, 2017, 52(1-2): 183–196

[29] Lv S, Liu C, Yao H, Zheng J. Comparisons of synchronous measurement methods on various moduli of asphalt mixtures. Construction & Building Materials, 2018, 158: 1035–1045

[30] Wang L, Hoyos L R, Mohammad L, Abadie C. Characterization of asphalt concrete by multi-stage true triaxial testing. Journal of ASTM International, 2005, 2(10): 12276

[31] Pan T, Tutumluer E, Anochie-Boateng J. Aggregate morphology affecting resilient behavior of unbound granular materials, Transportation Research Record. Transportation Research Record: Journal of the Transportation Research Board, 2006, 1952(1): 12–20

[32] Tutumluer E, Seyhan U. Laboratory determination of anisotropic aggregate resilient moduli using an innovative test device, Transportation Research Record. Transportation Research Record: Journal of the Transportation Research Board, 1999, 1687(1): 13–21

[33] Oda M. Initial fabrics and their relations to mechanical properties of granular material. Soils and foundations, 1972, 12(1): 17–36

[34] Alkhalifah T, Tsvankin I. Velocity analysis for transversely isotropic media. Geophysics, 1995, 60(5): 1550–1566

[35] You L, Yan K, Hu Y, Ma W. Impact of interlayer on the anisotropic multi-layered medium overlaying viscoelastic layer under axisymmetric loading. Applied Mathematical Modelling, 2018, 61: 726–743

[36] Kim S H, Little D N, Masad E, Lytton R L. Estimation of level of anisotropy in unbound granular layers considering aggregate physical properties. International Journal of Pavement Engineering, 2005, 6(4): 217–227

[37] Tutumluer E, Little D, Kim S H. Validated model for predicting field performance of aggregate base courses. Transportation Research Record: Journal of the Transportation Research Board, 2003, 1837(1): 41–49

[38] Masad S, Little D, Masad E. Analysis of flexible pavement response and performance using isotropic and anisotropic material properties. Journal of Transportation Engineering, 2006, 132(4): 342–349

[39] Wang L, Hoyos L R, Wang J, Voyiadjis G, Abadie C. Anisotropic properties of asphalt concrete: Characterization and implications for pavement design and analysis. Journal of Materials in Civil Engineering, 2005, 17(5): 535–543

[40] Mukhopadhyay A. Stresses produced by a normal load moving over a transversely isotropic layer of ice lying on a rigid foundation. Pure and Applied Geophysics, 1965, 60(1): 29–41

[41] Yan K, Xu H, You L. Analytical layer-element approach for wave propagation of transversely isotropic pavement. International Journal of Pavement Engineering, 2016, 17(3): 275–282

[42] You L, Yan K, Hu Y, Liu J, Ge D. Spectral element method for dynamic response of transversely isotropic asphalt pavement under impact load. Road Materials and Pavement Design, 2018, 19(1): 223–238

[43] Yan K, Shi T, You L, J Man. Spectral element method for dynamic response of multilayered half medium subjected to harmonic moving load. International Journal of Geomechanics, 2018, 18(12): 04018161

[44] You L, You Z, Dai Q, Xie X, Washko S, Gao J. Investigation of adhesion and interface bond strength for pavements underlying chip-seal: Effect of asphalt-aggregate combinations and freeze-thaw cycles on chip-seal. Construction & Building Materials, 2019, 203: 322–330

[45] Canestrari F, Ferrotti G, Lu X, Millien A, Partl M N, Petit C, Phelipot-Mardelé A, Piber H, Raab C. Mechanical testing of interlayer bonding in asphalt pavements. Advances in Interlaboratory Testing and Evaluation of Bituminous Materials, 2013, 9: 303–360

[46] Borawski B, Singh J, Todd J A, Wolfe D E. Multi-layer coating design architecture for optimum particulate erosion resistance. Wear, 2011, 271(11–12): 2782–2792

[47] Sui Y, Appenzeller J. Screening and interlayer coupling in multilayer graphene field-effect transistors. Nano Letters, 2009, 9(8): 2973–2977

[48] Chun S, Kim K, Greene J, Choubane B. Evaluation of interlayer bonding condition on structural response characteristics of asphalt pavement using finite element analysis and full-scale field tests. Construction & Building Materials, 2015, 96: 307–318

[49] Stenzler J S, Goulbourne N. Effect of Polyacrylate Interlayer microstructure on the impact response of multi-layered polymers. Time Dependent Constitutive Behavior and Fracture/Failure Processes, 2011, 3: 241–258

[50] Chupin O, Chabot A, Piau J M, Duhamel D. Influence of sliding interfaces on the response of a layered viscoelastic medium under a moving load. International Journal of Solids and Structures, 2010, 47(25–26): 3435–3446

[51] Liu N, Yan K, Shi C, You L. Influence of interface conditions on the response of transversely isotropic multi-layered medium by impact load. Journal of the Mechanical Behavior of Biomedical Materials, 2018, 77: 485–493

[52] You L, Yan K, Shi T, Man J, Liu N. Analytical solution for the effect of anisotropic layers/interlayers on an elastic multi-layered medium subjected to moving load. International Journal of Solids and Structures, 2019, 172: 10–20

[53] Jiang Y, Lee F C. Single-stage single-phase parallel power factor correction scheme. In: Power Electronics Specialists Conference, PESC’94 Record. IEEE, 1994, 1145–1151

[54] You L, You Z, Yan K. Effect of anisotropic characteristics on the mechanical behavior of asphalt concrete overlay. Frontiers of Structural and Civil Engineering, 2019, 13(1): 110–122

[55] Schmidt H, Tango G. Efficient global matrix approach to the computation of synthetic seismograms. Geophysical Journal International, 1986, 84(2): 331–359

[56] Zhang P, Liu J, Lin G, Wang W. Axisymmetric dynamic response of the multi-layered transversely isotropic medium. Soil Dynamics and Earthquake Engineering, 2015, 78: 1–18

[57] Cornille P. Computation of Hankel transforms. SIAM Review, 1972, 14(2): 278–285

[58] Kim J. General viscoelastic solutions for multilayered systems subjected to static and moving loads. Journal of Materials in Civil Engineering, 2011, 23(7): 1007–1016

[59] Ai Z, Cang N, Cheng C. Analytical layer-element method for axisymmetric problem of transversely isotropic multi-layered soils. Chinese Journal of Geotechnical Engineering, 2012, 34: 863–867

[60] Masad E, Tashman L, Somedavan N, Little D. Micromechanics-based analysis of stiffness anisotropy in asphalt mixtures. Journal of Materials in Civil Engineering, 2002, 14(5): 374–383

[61] Tutumluer E, Thompson M. Anisotropic modeling of granular bases in flexible pavements. Transportation Research Record: Journal of the Transportation Research Board, 1997, 1577(1): 18–26

[62] Vu-Bac N, Lahmer T, Keitel H, Zhao J, Zhuang X, Rabczuk T. Stochastic predictions of bulk properties of amorphous polyethylene based on molecular dynamics simulations. Mechanics of Materials, 2014, 68: 70–84

[63] Vu-Bac N, Lahmer T, Zhang Y, Zhuang X, Rabczuk T. Stochastic predictions of interfacial characteristic of polymeric nanocomposites (PNCs). Composites. Part B, Engineering, 2014, 59: 80–95

[64] Vu-Bac N, Lahmer T, Zhuang X, Nguyen-Thoi T, Rabczuk T. A software framework for probabilistic sensitivity analysis for computationally expensive models. Advances in Engineering Software, 2016, 100: 19–31

[65] Vu-Bac N, Rafiee R, Zhuang X, Lahmer T, Rabczuk T. Uncertainty quantification for multiscale modeling of polymer nanocomposites with correlated parameters. Composites. Part B, Engineering, 2015, 68: 446–464

[66] Vu-Bac N, Silani M, Lahmer T, Zhuang X, Rabczuk T. A unified framework for stochastic predictions of mechanical properties of polymeric nanocomposites. Computational Materials Science, 2015, 96: 520–535

RIGHTS & PERMISSIONS

Higher Education Press