Frictional contact problem of one-dimensional hexagonal piezoelectric quasicrystals layer

Based on three-dimensional (3D) general solutions for one-dimensional (1D) hexagonal piezoelectric quasicrystals (PEQCs), this paper studied the frictional contact problem of 1D hexagonal PEQCs layer. The frequency response functions for 1D hexagonal PEQCs layer are analytically derived by applying double Fourier integral transforms to the general solutions and boundary conditions, which are consequently converted to the corresponding influence coefficients. The conjugate gradient method is used to obtain the unknown pressure distribution, while the discrete convolution–fast Fourier transform technique is applied to calculate the displacements and stresses of phonon and phason, electric potentials and electric displacements. Numerical results are given to reveal the influences of layer thickness, material parameters and loading conditions on the contact behavior. The obtained 3D contact solutions are not only helpful for further analysis and understanding of the coupling characteristics of phonon, phason and electric field, but also provide a reference basis for experimental analysis and material development.

Piezoelectricity is an important property of QCs. By virtue of the piezoelectric effect, QCs have better sensor and actuator functions in the design of intelligent structures and systems. Fujiwara et al. [4] discovered the piezoelectric effect of QCs for the first time and discussed the electronic structure and electron transport of QCs. Hu et al. [5] proposed the elastic theory framework of PEQCs and pointed out that the piezoelectric effect of QCs is excited by phonon and phason field simultaneously. Rao et al. [6] obtained the number of independent nonvanishing piezoelectric coefficients. With the development of preparation technology of piezoelectric quasicrystals (PEQCs) and actual requirement to resist contact damage, the relevant researches on contact problem becomes more and more important. Scanning probe/piezoelectric response force microscope (PFM) is a very important technology, which can be used to test and analyze the properties of materials. These techniques all depend on the relevant 3D contact theory analysis of materials. In conclusion, the 3D contact analysis of PEQCs layer has significantly theoretical values, which also attracts many researchers' attention. Li and Liu [7] derived the character formulae of representation matrices of the piezoelectric property tensors for the 1D-QCs. Altay and Dokmeci [8] obtained 3D fundamental equations of PEQCs. The 3D general solution of 1D hexagonal PEQCs according to the rigorous operator theory and Almansi theorem is given by Chen and Li [9,10]. By using the same method, Gao et al. [11][12][13] derived the 3D elastic general solution of 2D-PEQCs. Based on the general solution of QCs, Li et al. [14] obtained the fundamental solution of 1D hexagonal PEQCs in 2014. Based on the general solutions of QCs, Fan et al. [15] obtained the analytical solution of the contact problem of 2D-QCs by using the complex function method. Zhao and Li [16,17] studied contact problems of QCs by using analytical function theory and dual integral equation theory. Based on generalized potential function theory, the frictionless contact problem of 1D and 2D-QCs half-spaces is considered by Wu and Li [18,19]. Using the trial and error method, Hou et al. [20,21] studied the 3D fundamental solution of 1D hexagonal QCs layer. Li and Zhou [22] studied the fundamental solutions and frictionless contact problem in a semi-infinite space of 2D hexagonal PEQCs.
The study on 3D contacts would provide more general solutions and more accurate understanding to real-world problems [23,24]. However, little has been done on 3D contact problems of PEQCs layers. The semi-analytical models (SAMs), supported by the conjugate gradient method (CGM) and the discrete convolution-fast Fourier transform ( DC-FFT) [25][26][27], appear to be efficient in solving 3D contact problems involving multifield coupling [28][29][30][31]. This paper reports the development of the frictional contact of 1D hexagonal PEQCs layer, based on the SAMs, by the first mathematical derivation of the closed-form FRFs for 1D hexagonal PEQCs layer and creative extension of the numerical DC-FFT framework. The contact responses of displacements and stresses of phonon and phason, electric potentials and electric displacements are solved. The first set of numerical results is given to reveal the influences of layer thickness, material parameters and loading conditions on the contact behavior.

Basic equations
1D hexagonal PEQCs are transversely isotropic materials in which the phonon, phason and electric field are coupled. Phonon is the quantum of elastic vibration. Phason is the local rearrangement of atoms in the cell based on Penrose. In the Cartesian coordinate system, it is assumed that the atoms are arranged periodically in the x-y plane and quasiperiodic arrangement along z-axis. Therefore, the constitutive equations for 1D hexagonal PEQCs are [14] where σ i j , ε i j and C i j are the stresses, strains and elastic stiffness constants in the phonon field, respectively; H i j , w i j and K i are the corresponding quantities in the phason field, respectively; D i , E i and ξ i j are the electric displacements, electric fields and dielectric coefficients, respectively; R i is the phonon-phason coupling elastic constant; e i j and e i j are the piezoelectric coefficients.
In the phason and electric fields, there are similar geometric relations between the strain and displacement component. The generalized strain-displacement relations (geometric equations) are where u i , w z and φ are phonon displacement, phason displacement and electric potential, respectively. In the absence of body forces and free charges, the equilibrium equations can be expressed as By substituting Eqs. (1) and (2) into Eq. (3), the following governing equations can be obtained where ∂ 2 ∂ x 2 + ∂ 2 ∂ y 2 is the plane Laplacian operator.

General solutions
Using the rigorous operator theory and generalized Almansi theorem [32], the 3D general solutions are [14] where the potential function ψ i satisfies the quasi-Laplace equation where λ ki , γ ki and ρ i are given in "Appendix A" as well.
2.3 Frictional contact of 1D hexagonal PEQCs layer Figure 1 illustrates the contact problem considered in which the 1D hexagonal PEQCs indenter (a sphere of radius R) slides on the 1D hexagonal PEQC layer along the x direction. The 1D hexagonal PEQCs layer with finite thickness of h was deposited on the rigid grounding substrate. The normal force P z , tangential force P x , normal phason force H and electric charge Q are applied on the indenter, which lead to pressure p z , shear traction p x , phason force h z and electric charge density q z at the contact interface S (contact radius a). Obviously, s p z dxdy P z , s p x dxdy P x , s h z dxdy H , s q z dxdy Q. Since there is no experimental study on how to apply the phason force on the QCs surface, but there are some artificial methods that can cause the disorder of QCs phase field [33,34]. To show the responses under the phason force from the perspective of theoretical study, we assume the indenter can induce not only the phason force but also the phason displacement, and the applied phason force is equivalent inside the contact region. Furthermore, in the sliding contact, friction follows the Coulomb friction law, so that p x μp z , with μ for friction coefficient. The boundary conditions at the upper surface of 1D hexagonal PEQCs layer can be written as [28] The lower surface of 1D hexagonal PEQCs layer is perfectly bonded to a rigid grounded substrate, so that the boundary conditions at the lower surface are In order to determine the potential functions ψ i , the following direct double Fourier transforms and its inverse aref where j √ −1 without specification elsewhere; Hat "≈" denotes the double Fourier transform; m and n represent the transform variables with respect to the x and y directions, respectively.
Applying Fourier transforms to Eq. (6) results in the following exponential form of potential functions ψ ĩ where α √ m 2 + n 2 ; A i (m,n) and A i (m,n) are ten unknowns to be determined from the boundary conditions. Applying Fourier transforms to Eqs. (5) and (7), combined with Eq. (11), the transformed displacements and stresses of phonon and phason field, electric potential and electric displacement are expressed as Table 1 The material constants of 1D hexagonal PEQCs [14]. C i j , K i and R i are in GPa, e i j and e i j are in Cm −2 and ξ i j are in Applying Fourier transforms to Eqs. (8) and (9), we have The following results can be obtained by solving the above equations where ki are given in "Appendix B." Accordingly, the Fourier transform solutions for the corresponding displacements and stresses of phonon and phason field, electric potential and electric displacement are established by substituting the coefficients The transformed solutionũ z (m, n, z) of the layer is function of the transformed pressurep z , phason forcẽ h z and electric charge densityq z , which can be further expressed as whereG p q denotes the frequency response functions (FRFs) response of p excited by load q. Take surface normal displacement u z as an example,G u z qz represents the FRFs of u z caused by q z . The transformed solution is analytically expressed in closed forms in the Fourier-frequency domain, as shown in Eq. (15), which should be further integrated to be converted into the space domain from the Fourierfrequency domain. Applying the inverse Fourier transforms to Eq. (15) with respect to the frequency numbers m and n leads to It is difficult to solve the complicated integration in Eq. (16) analytically; however, this solution is for the space-domain convolutions between surface excitations and mechanical responses. Therefore, DC-FFT is employed to solve the convolutions in the present contact problems.

DC-FFT algorithm
In the Fourier domain, the FRFs represent the material response of the Dirac delta functions (a generalized unit point load at the coordinate origin). Therefore, by lettingp z 1(h z 0,q z 0),h z 1(p z 0,q z 0), q z 1(p z 0,h z 0) and H zz ≡H z , respectively, in (12), the components of the FRFs can be obtained analytically as Once the FRFs are given, the components of continuous Fourier transforms of influence coefficients (ICs) can be solved by multiplying the FRFs and the shape functions of Fourier transforms [26] where Y is the shape function given as follows The where x and y are the mesh sizes along the x and y directions, respectively; AL represents the level of the aliasing control; M e 2 γ M and N e 2 γ N denote the refinement mesh numbers of the original mesh numbers M and N ;K is the discrete Fourier transform form of ICs. Then, the DC-FFT algorithm is implemented to process pressure, phason force, electric charge by means of the corresponding components of the ICs. The normal displacements and stresses of phonon and phason, electric potential and electric displacement are given by the following inverse fast Fourier transform (IFFT), respectively.
The above formulation shows the DC-FFT algorithm for the 1D hexagonal PEQCs layer, which can also be applied to the indenter by letting the layer thickness to infinite.

Contact conditions and numerical modeling
Now the contact problem of the indenter sliding on the 1D hexagonal PEQCs layer can be solved, as shown in Fig. 1. The mechanical loadings at the contact interface S can be represented by pressure p z , shear traction p x with p x μp z , phason force h z and electric charge density q z . Therefore, the following contact conditions apply to such a contact problem where u z1 and u z2 are the surface normal displacement for the indenter and layer, respectively; g is the gap between the indenter and layer; h 0 is the profile of indenter. Pressure p z in Eq. (22) can be solved with the conjugate gradient method (CGM) described in Ref. [25].

Results and discussion
In the numerical process for solving the problem shown in Fig. 1, the fixed computational domain of 4 a 0 × 4 a 0 × h with 256 × 256 × 64 grids in the x, y and z directions is employed, where a 0 is the maximum contact radius, which will be discussed later. The material properties are given in Table 1.

Model verification
The correctness of the model is verified by comparison the 1D hexagonal PEQCs layer contact solutions of h → ∞, e i j , e i j , ξ i j → 0 with that of Wu et al. [18]. For the problem of 1D hexagonal PEQCs layer indented by a rigid insulating sphere with radius R 1 mm under P z 12 N,h 10a 0 , e i j , e i j , ξ i j 10 −30 , the  subsurface phonon normal stress, shear stress and phason stress are identical to those obtained by Wu et al. [18], as shown in Fig. 2 (a,b,c). The model is further verified through degenerating our solutions from "PEQCs layer" to "PE half-space" by taking "h → ∞, K i , R i → 0." For the problem of the PZT-4 half-space indented by a rigid insulating sphere with radius R 1 mm under normal force P z 12, maximum contact radius a 0 0.9086 3 √ cP z R contact pressure p 0 P z (πa 2 0 ), electric potential φ 0 P z · 10 −2 , according to Ref. [37], we can easy to calculate c 4.06 × 10 −12 , the surface normal stress and electric potential are identical to those obtained by Ding et al. [37], as shown in Fig. 3.
The relationship between indentation depth d u z (0, 0, z) and indentation load P z , P x , H , Q is shown as in Fig. 4, where radius R 1 μm (or 200 μm), friction coefficient μ 0 (or 0.5), layer thickness h a 0 (or 10 a 0 ), normal force P z ranging from 0 to 10 mN, phason force H 0 (or 2 mN) and electric charge Q 0(or 1 mC). It can be seen that in order to achieve the same indentation depth, the sphere with larger radius needs more force. The 1D hexagonal PEQCs layer h a 0 needs more force than h 10 a 0 . The indentation force of rigid ball is larger than that of the corresponding 1D hexagonal PQCS ball. Figure 4 also reveals that increasing phason force have a negligible influence on the indentation depth. For 1D hexagonal PEQCs layer h 10 a 0 , increasing friction coefficient μ and electric charge Q will increase the indentation depth. However, the opposite tendency is true for 1D hexagonal PEQCs layer h a 0 . In addition, when both the indenter and layer are 1D hexagonal PEQCs, the friction force has a greater influence on the indentation depth. These curves are in accordance with the law of pressure variation.

Effect of parameters
The problem of the sliding contact between the 1D hexagonal PEQCs layer and the 1D hexagonal PEQCs sphere is investigated, where the applied force P z 0.2 mN, phason force H 2 mN, sphere radius R 1 μm, layer thickness h 0.5 a 0 (or 1.0 a 0 or 2.0 a 0 ), the friction coefficient μ 0 (or 0.3 or 0.5 or 0.8) and electric charge Q 0(or 1 mC or 4 mC or 8 mC). For the convenience of display, the numerical results are normalized by the referenced parameters of maximum contact pressure p 0 , contact radius a 0 , phason displacement w 0 and electric potential φ 0 . For 1D hexagonal PEQCs half-space, according to Ref. [38], we can calculate the maximum indentation depth δ 0 a 2 0 R, a 0 3 3πξ 11 R P z 4, p 0 2a 0 (π 2 ξ 11 R), w 0 ξ 12 δ 0 ξ 11 ,φ 0 ξ 13 δ ξ 11 , ξ 11 3.3779 × 10 −12 m 2 N, ξ 12 1.1381 × 10 −11 m 2 N and ξ 13 2.4084 × 10 −3 m 2 C. Figures 5, 6, 7, 8, 9 and 10 show the dimensionless surface (subsurface) distributions for the pressure, phonon stress, phason displacement and electric potential, respectively. It can be found that friction coefficient does not affect the contact area, but the shape of the pressure distribution is changed. The maximum values of pressure, phonon stress, phason displacement and electric potential are shifted away from the sliding direction. In addition, when the electric charge Q is larger, then the pressure and electric potential become larger, and the phason displacement becomes smaller. The increase in Q value leads to larger value of the maximum pressure and electric potential, smaller value of the maximum phason displacement, which reveals that the contact interface appears "harder" when a larger electric charge Q is applied to the material. The peak phason displacement increases with an increasing of layer thickness h. A thicker layer leads to a smaller peak value of the contact pressure and electric potential. The increase in layer thickness h pushes the peak pressure, electric potential and phason displacement closer to the peak value in half-space case. Moreover, the electric potential and phason displacement extend out of the contact region and do not decrease to zero even at the edges of the contact regions of both surfaces. The surface pressure, electric potential and phason displacement decrease with an increasing of the distance far away from the contact surface. Figures 11 and 12 plot the diagrams of the subsurface phason stresses and electric displacements under different loads. The larger friction coefficient μ is, the smaller subsurface phason stresses are, and the greater electric displacements are. Besides, with the increase in Q value, the phason stresses increase gradually, and electric displacements decrease. In general, subsurface phason stresses and electric displacements decrease with the increase in layer thickness h. Meanwhile, the position of maximum phason stresses and electric displacements deviates slightly from the contact surface. It can be seen from Fig. 13 that the von Mises stresses increase with an increasing of the friction coefficient μ. The increase in μ causes the location of the maximum von Mises stress to move toward the contact surface. It is clear that friction has significant effect on the von Mises stress, but influence of electric charge Q can be negligible. In addition, the von Mises stresses becomes smaller with an increasing of h. The stresses show obvious distortion due to friction. This distortion is somewhat alleviated by an increasing of layer thickness h. Figure 14 reveals that the contact behaviors affected by layer thickness h a 0 . For sufficiently large h a 0 , the maximum dimensionless u z a 0 , p z p 0 , and σ s √ 3 p 0 do not change their value along with the increasing layer thickness, and the layer problem becomes a half-space problem. However, w z w 0 and φ φ 0 approach certain asymptotic values for larger h a 0 , the phason and electric responses penetrate far deep from the contact surface, therefore, the influence of substrate cannot be neglected. The sensitivity function suggested by Makagon et al. [39] was employed to evaluate the coupling effects of material properties on contact behavior V. The dimensionless sensitivities below are numerically calculated with step length 0.01 f 0 : S( f ) where f is a selected material constant with f 0 being its reference value. Figure 15 reveals the strong influence of the elastic stiffnesses C 11 , C 13 , C 33 , C 44 on u z , p z , w z and φ, while the influence of C 12 is weak. The phonon-phason coupling elastic constants R i , phason elastic stiffness constants K i , dielectric coefficient ξ 33 , piezoelectric coefficients e 33

Model extension
Now, we extend the model to the contact behavior involving indenters of different shapes. The sliding contact between the 1D hexagonal PEQCs layer with an ellipsoidal 1D hexagonal PEQC punch is investigated, where the applied force P z 0.2 mN, major and minor radii R x 2 μm and R y 1 μm, layer thickness h 2 a 0 , the friction coefficient μ 0.5, phason force H 2 mN and electric charge Q 1 mC. From Figs. 16, 17, 18 and 19, we draw the following conclusions: Compared with the spherical indenter, it is more difficult for the ellipsoidal indenter to press in 1D hexagonal PEQCs layer. The location of the maximum contact pressure, phonon stress, phason displacement, electric potential, von Mises stresses, phason stress and electric displacement is closer to the surface. The electric charge, friction coefficient and layer thickness still strongly influence contact behavior. For example, the value of von Mises stress decreases with the increase in layer thickness h. The location of the maximum von Mises stress is somewhat affected by the surface friction.

Conclusions
In this paper, the frictional contact problem of 1D hexagonal PEQCs layer has been solved. The 3D general solutions for the problem of 1D hexagonal PEQCs are obtained by using rigorous operator theory and Almansi theorem. The first analytical FRFs of displacements and stresses of phonon and phason, electric potentials and electric displacements are obtained by applying double Fourier integral transforms. The influences of layer thickness, material properties, surface pressure, phason force, electric charge and friction are numerically investigated by DC-FFT.
The major advantage of the current algorithm is high calculation efficiency, which can solve complicated contact problems in a few minutes. Because the influence coefficients are independent of indenter shapes, the proposed semi-analytical model can be conveniently implemented to various indenters, such as flat-ended rectangular, circular and elliptical punches, spherical and ellipsoidal indenters, cylinders of a finite length, multi-connection contacts and multi-domain contacts. It is not only convenient for surface geometry treatment, but also can deal with complicated surface topographic features such as rough surfaces. Therefore, the model has a wide application prospect in the design of multi-functional engineering components.

Appendix B
The intermediate variables ki , a