Vibration of a Two-Layer “Metal+PZT” Plate Contacting with Viscous Fluid

1  Introduction

In the classical sense, the study of the dynamics of a piecewise homogeneous acoustic medium, such as “plate+fluid” systems in the cases where the plate is made of conventional metal or polymer materials, is necessary for various fields of modern industry, e.g., fluid transport, geophysical investigations, biomedical device construction, submarine construction and sound insulation. The results of these studies provide a theoretical basis for understanding the acoustic phenomena occurring in the piecewise homogeneous acoustic medium. Consequently, these results provide opportunities for the control and management of these phenomena.

In modern literature, the work of Lamb [1], completed a hundred years ago, is considered the beginning of these studies. Since then, a lot of further research has been done. For a detailed review, as well as some recent results, see the works of Amabili [2], Akbarov [3,4], Akbarov et al. [5–7], Guz [8], Guz et al. [9], Paimushin et al. [10], Paimushin et al. [11], Shuaib et al. [12], Sorokin et al. [13], Zamanov et al. [14] and many others listed therein.

Note that these investigations can be classified according to various aspects, such as the theories applied to describe the motion of the plate (approximate plate theories, which are described, for example in [1,2,10–12] and many other theories discussed in the work of Akbarov [3]); the exact equations of elastodynamics discussed for example in [5–9]; and modeling of the plate material as elastic, for example, in the work of Akbarov et al. [6] and many other works listed therein, or viscoelastic, for example, in the works [7,14]. Moreover, these works can be classified according to the type of dynamic problems considered (wave dispersion, free vibration, and forced vibration), the modeling of the fluid (inviscid fluid, viscous fluid, compressible and incompressible fluid), and so on.

The aforementioned works serve not only their direct purpose but also provide the theoretical basis for the study of the corresponding problems related to the cases in which the plate material is intelligent, such as piezoelectric or piezomagnetic. It should be noted that the results of such investigations have great importance for energy harvesting procedures (see, for example, references [15–18] and others listed therein), and in hydroacoustic transducers that receive (or generate) acoustic waves (see, for example references [19–21] and others listed therein).

Nevertheless, there has not been enough fundamental theoretical research in this field, to which the present research is dedicated, in particular, to the mechanical forced vibration of the hydro-piezoelectric system consisting of a two-layer “metal+PZT” plate, a compressible viscous fluid, and a rigid wall. The research presented here is carried out in the framework of the piecewise-homogeneous body model using the exact equations and relations of elastodynamics and electrodynamics in the plane-strain state to describe the motion of the “metal+piezoelectric” plate and using the linearized Navier-Stokes equations for the plane flow of the compressible viscous fluid.

To show the relevance and importance of the present work, let us consider a brief review of related research and begin with the work of Belkourchia et al. [16], in which the algorithm is developed for the numerical solution of the problem related to determination of the pressure exerted by ocean motion waves on a vertical cantilever beam with piezoelectric patches where the motion of the cantilever is described within the framework of the Euler-Bernoulli beam theory, but the flow of the fluid is described by the Navier-Stokes equations for the incompressible viscous fluid. This does not take into account the piezoelectricity of the patches, and it is assumed that determination of the fluid pressures acting on the left and right sides of the cantilever allows for calculation of the electric potentials occurring in the piezoelectric patches. In this way, extraction of the energy originating from the fluid motion can be controlled. The work of Akaydin et al. [22] and others listed therein also use a similar approach, but the material of the cantilever is taken as a piezoelectric one. In addition, the Euler-Bernoulli beam model is used in many studies listed in the conference book [23] on this topic.

The piezoelectric composite cantilever beam model is also used in the work of He et al. [24] to determine the optimal thickness of the PZT layer in the piezoelectric energy harvesters (PEH) to obtain the maximum electrical energy under static mechanical and external electrical excitation. FEM (finite element method) is used in the framework of the exact relations of linear piezoelectricity to obtain concrete numerical results on the influence of the volume fraction of the PZT material on electrical energy. However, the approach developed in the work of He et al. [24] is not only related to the cantilever beam model, but also to more general PEH models. In the work of Chen et al. [25], the effect of geometric nonlinearity and flexoelectricity on the uniform and bimorph PEH under harmonic mechanical excitation is analyzed. The model of a cantilever beam is also considered and the motion of these beams is described in the framework of the Euler-Bernoulli beam theory. Hamilton’s principle is used to obtain the equations of motion and these equations are solved using Galerkin’s method.

In contrast to the aforementioned work, the work by Amini et al. [15] attempts to use the three-dimensional exact equations of motion of elasto-electrodynamics to describe the motion of the piezoelectric cantilever. There are also many studies of the dynamics of intelligent systems that relate to linear and nonlinear vibrations, and dispersion of the waves propagating therein. Examples of such recent investigations can be found in the papers [26–28] and many others presented therein. Let us consider briefly the subjects of these works and begin with the work of Jiang et al. [26] which presents a nonlinear piezoelectric energy harvester with two degrees of freedom to improve the efficiency of energy harvesting at low frequencies. An L-shaped piezoelectric cantilever beam is considered, and this beam is assumed to be excited by nonlinear attractive forces of magnets attached to the end of the cantilever and to the base. The dynamic equation for the PZT beam is derived from the magnetic attraction force mentioned above. Under the harmonic excitation, this equation is solved numerically and the range of change of the problem parameters is determined, which improves the harvesting properties of the system.

The work of Nie et al. [27] investigates the dispersion of the propagating Rayleigh-type waves in the system consisting of the piezoelectric layer and a semi-infinite dielectric substrate. The material of the piezoelectric layer is assumed to be a cubic crystal and the contact conditions on the interface of the constituents of the system are assumed to be mechanically perfect but electrically imperfect. Dispersion relations for electrically open or short-circuited boundary conditions are obtained and analyzed.

The paper by Qiao et al. [28] presents a theoretical analysis of the technology for vibration control of wind turbine blades by using piezoelectric smart structures. Since the design of the blade structure made of piezoelectric material is similar to a shallow shell structure, the equations of motion are derived based on the theory of piezoelectric laminated shallow shells. Based on the solution of these equations, it is found that the blade vibration under wind load can be reduced by applying control voltage.

In all the above works dealing with hydro-piezoelectric systems, the fluid flow is described in the framework of the incompressible viscous or inviscid Newtonian fluid model using the linear Navier-Stokes equations. Moreover, in these works, the motion (or vibration) of the piezoelectric cantilevers is considered with finite sizes, so such a model cannot discover the local electromechanical and hydro-electromechanical coupling effects without the influence of the boundary conditions. Such an effect related to the subject of the present work can be discovered in the framework of infinite two-layer “metal+piezoelectric” plate and fluid systems.

Note also that the above work does not investigate how the electromechanical coupling effect of the piezoelectric material can affect the pressure at the interface between the fluid and the plate. Since the energy generation in the energy harvesting piezoelectric systems comes from this pressure, the above coupling effects must be considered in theoretical studies.

In order to address the above issues, Ekicioglu Kuzeci [29] carried out some investigations in this field and studied the mechanically forced vibration of the hydro-piezoelectric system consisting of the infinite piezoelectric plate and the compressible (barotropic) Newtonian viscous fluid with finite depth. The flow of the fluid is described by the linearized Navier-Stokes equations, while the motion of the plate is described by the exact equations of the theory of electroelasticity for piezoelectric materials.

In the present paper, an attempt is made to develop these investigations for the system consisting of a two-layer “metal+piezoelectric” plate, a compressible viscous fluid, and a rigid wall. The motion of the plate is described within the piecewise-homogeneous body model using the exact equations of elasto-and electrodynamics for piezoelectric materials, and the flow of the fluid is described using the linearized Navier-Stokes equations for a compressible viscous fluid.

The present work can also be considered as further development of the second author’s work from the last decade dealing with the forced vibration of single-layer purely elastic or viscoelastic plates in contact with a compressible viscous fluid (see, for example, the works [5–7]). A review of these works can be found in the paper by Akbarov [3]. Note that a number of studies have also been carried out on the problem of the dispersion of waves in hydro-elastic systems consisting of a single-layer purely elastic plate and a compressible fluid. Their data and considerations are contained in the works of Guz et al. [9] and Guz [8].

Thus, it is clear that the investigations on the subject of the present work were carried out for the case where the plate is single-layered. Thus, the main difference of the present work from the previous related works of the second author and his collaborators, as well as from the investigations described in the works [8,9], is the two-layer nature of the plate and the piezoelectricity of the material of one of the layers.

Note that in the formulation of the problem studied in this paper, we use the linearized Navier-Stokes equations obtained from the corresponding nonlinear Navier-Stokes equations for compressible barotropic viscous fluids. Linearization is one of the methods for solving nonlinear problems. Many other methods for solving nonlinear physical-mechanical problems have been developed in recent years. We briefly review some here and begin with the work of Mahdy et al. [30], in which the reduced differential transform method is used to solve the nonlinear fractional model tumor-immune. The same method is also used in the work of Gepreel et al. [31] to solve nonlinear biomathematical problems related to the spread of viruses in a computer network and to the problem of susceptible infected-recovered (SIR) childhood diseases.

Moreover, in a series of studies by Bazighifan et al. [32], Santra et al. [33], Bazighifan et al. [34], Bazighifan et al. [35], Moaaz et al. [36], and many others listed therein, the nonlinear differential equations of even order are studied with p-Laplacian-like operators, the results of which shed light on the solution of the corresponding nonlinear problems in continuum mechanics. In the work of Bazighifan [37], a new type of oscillation criterion is established for a neutral differential equation of even order. In the work of Santra et al. [33], sufficient conditions for oscillation of all solutions of a second-order functional differential equation are found. New results on the oscillation of solutions of a class of advanced even-order differential equations with p-Laplacian-like operators are obtained in the work of Bazighifan et al. [35]. The asymptotic peculiarities for a class of fourth order neutral differential equations are studied in the work of Bazighifan et al. [34]. New sufficient conditions for the oscillation of fourth order neutral differential equations are found in the work of Moaaz et al. [36]. The work of Moaaz et al. [38] deals with the oscillation of a class of third order nonlinear delay differential equations with mean term and a new description of the oscillation of third order linear differential equation without damping is proposed. Finally, in the work of Moaaz et al. [39], the general class of differential equations is considered (the particular cases of which have been studied previously by the authors of the related works mentioned above) and the corresponding general theorems for the asymptotic behavior of the solution of these equations are established.

In this context, since the present study refers to coupling field problems, we refer to the recent studies of Mahdy et al. [40], Mohammed et al. [41], Khamis et al. [42], Mahdy et al. [43] and many others within these works, which also refer to coupling field problems. The novel model of photothermoelasticity theory is studied in the work [40]. A novel mathematical model is proposed and applied in [41] to study the micro stretch properties of an elastic semiconductor medium. In papers [42,43], novel mathematical models are developed for the piezoelectric elastic semiconductor medium and photothermoelasticity theory, respectively.

The remainder of this paper is organized as follows. The mathematical formulation of the problem, the equation of motion of the metal/PZT plate and the equations of flow, as well as the corresponding boundary, contact, and compatibility conditions are given in Section 2. The solution method for the formulated boundary value problem and the algorithm for obtaining the numerical results are explained in Section 3. The numerical results are presented and discussed in Section 4, and finally, the results are summarized in Section 5.

2  Formulation of the Problem and Governing Field Equations and Relations

We consider the hydro-elastic-piezoelectric system shown in Fig. 1, which consists of a two-layer “metal+piezoelectric” plate, a compressible barotropic viscous fluid, and a rigid wall that limits the fluid depth. Associate the Cartesian coordinate system Ox1x2x3 with the upper surface plane of the metal layer of the plate, according to which the metal and piezoelectric layers of the plate occupy the regions {−∞<x1<+∞,−hm<x2<0,−∞<x3<+∞} and {−∞<x1<+∞,−hm−hp<x2<−hm,−∞<x3<+∞}, where h=(hm+hp) is the plate thickness, hm (hp) is the thickness of the metal (piezoelectric) layer of the plate, and the fluid occupies the region {−∞<x1<+∞,−h−hd<x2<−h,−∞<x3<+∞}, where hd is the fluid depth, i.e., the distance between the lower surface plane of the plate and the rigid wall. Let us assume that the direction of the Ox3 axis is perpendicular to the plane of Fig. 1. Since we consider the state of plane strain in the plate and the two-dimensional flow of the fluid on the plane Ox1x2, this axis is not seen in Fig. 1. We also assume that under −∞<x3<+∞,x1=0 and x2=0, uniformly distributed time-harmonic line forces with intensity P0 act on the upper plane of the plate.

Figure 1: The sketch of the hydro-elastic-piezoelectric system

In this framework, it is required to determine the mechanical and electrical fields in the above hydro-elastic-piezoelectric system by using the appropriate exact field equations and relations. For this purpose, we write these field equations and relations for each component of the system, using the upper indices (m) and (p) to indicate the affiliation of the quantities to the metallic and piezoelectric layers, respectively. According to the monograph by Yang [44], the equations of motion for the metal and piezoelectric layers can be written as follows:

∂σ11(α)∂x1+∂σ12(α)∂x2=ρ(α)∂2u1(α)∂t2,∂σ12(α)∂x1+∂σ22(α)∂x2=ρ(α)∂2u2(α)∂t2,∂D1(p)∂x1+∂D2(p)∂x2=0,(1)

where (α) in (1) is (m) for the metal layer and (α) is (p) for the piezoelectric layer.

In addition, the following relations are:

Mechanical relations for the metal layer material:

σ11(m)=λ(m)(γ11(m)+γ22(m))+2μ(m)γ11(m),σ22(m)=λ(m)(γ11(m)+γ22(m))+2μ(m)γ22(m),σ12(m)=2μ(m)γ12(m).(2)

Electro–mechanical relations for the piezoelectric layer material when it is polarized along the Ox2 axis:

σ11(p)=c11(p)γ11(p)+c13(p)γ22(p)−e31(p)E2(p),σ22(p)=c13(p)γ11(p)+c33(p)γ22(p)−e33(p)E2(p),σ12(p)=c44(p)(γ12(p)+γ21(p))−e15(p)E1(p),D1(p)=e15(p)(γ12(p)+γ21(p))+ε11(p)E1(p),D2(p)=e31(p)γ11(p)+e33(p)γ22(p)+ε33(p)E2(p),E1(p)=−∂φ(p)∂x1,E2(p)=−∂φ(p)∂x2.(3)

In (2) and (3), γ11(α), γ22(α), and γ12(α)(=γ21(α)) (where (α) is (m) for the metal layer, and (α) is (p) for the piezoelectric layer) are the components of the mechanical strain tensor and these components are determined through the components of the displacement vectors by the following relations:

γ11(α)=∂u1(α)∂x1,γ22(α)=∂u2(α)∂x2,γ21(α)=γ12(α)=12(∂u1(α)∂x2+∂u2(α)∂x1).(4)

In Eqs. (1)–(4), σ11(α), σ12(α), and σ22(α) are the components of the mechanical stress tensor in the plane-strain state; u1(α) and u2(1) are the components of the mechanical displacement vector; D1(p) and D2(p) are the components of the electrical displacement vector; φ(p) is an electric potential; c11(p), c33(p), c13(p), and c44(p) are elastic constants; e31(p), e33(p) and e15(p) are piezoelectric constants; ε11(p) and ε33(p) are dielectric constants; and λ(m) and μ(m) are Lame constants of the metal layer material.

This completes the field equations and relations describing the motion of the two-layer “metal+piezoelectric” plate. Now we try to write the field equations and relations for the compressible (barotropic) viscous fluid flow. According to the monograph of Guz [8], the linearized equations of the fluid flow for the considered case can be represented as follows:

∂ρ(1)∂t+ρ0(1)(∂V1∂x1+∂V2∂x2)=0,ρ0(1)∂V1∂t−μ(1)(∂2V1∂x12+∂2V1∂x22)+∂p(1)∂x1−(μ(1)+λ(1))∂θ∂x1=0,ρ0(1)∂V2∂t−μ(1)(∂2V2∂x12+∂2V2∂x22)+∂p(1)∂x2−(μ(1)+λ(1))∂θ∂x2=0.(5)

where the first equation of (5) is the continuity equation, but the last two equations of (5) are the linearized Navier-Stokes equations, which must be provided with the linearized rheological relations and the linearized equation of state given below:

T11=−p(1)+λ(1)θ+2μ(1)e11,T22=−p(1)+λ(1)θ+2μ(1)e22,T33=−p(1)+λ(1)θ,T12=2μ(1)e12,p(1)=a02ρ(1),a02=∂p0(1)∂ρ0(1)(6)

where

e11=∂V1∂x1,e22=∂V2∂x2,e12=12(∂V1∂x2+∂V2∂x1),θ=∂V1∂x1+∂V2∂x2.(7)

In (5)–(7), the following notation is used: ρ0(1) is the fluid density before perturbation; p0(1) is the hydrostatic pressure before perturbation; ρ(1) is the perturbation of the fluid density; p(1) is the perturbation of the hydrostatic pressure; V1 and V2 are the components of the fluid flow velocity vector in the directions of the Ox1 and Ox2 axes, respectively; Tij and eij(ij=11;22;33) are the components of the stress and strain velocity tensors in the fluid; a0 is the sound velocity in the fluid; and λ(1) and μ(1) are the coefficients of the fluid viscosity.

According to the monograph by Guz [8], we also add the following representations for the velocities V1 and V2, and the pressure p(1) to the previous equations and relations:

V1=∂φ(1)∂x1+∂ψ(1)∂x2,V2=∂φ(1)∂x2−∂ψ(1)∂x1,p(1)=ρ0(1)(λ(1)+2μ(1)ρ0(1)Δ−∂∂t)φ(1)(8)

where the potentials φ(1) and ψ(1) satisfy the following equations:

[(1+λ(1)+2μ(1)a02ρ0(1)∂∂t)Δ−1a02∂2∂t2]φ(1)=0,(ν(1)Δ−∂∂t)ψ(1)=0,Δ=∂2∂x12+∂2∂x22(9)

where ν(1) is the kinematic viscosity, i.e., ν(1)=μ(1)/ρ0(1). Supposing that p(1)=−(T11+T22+T33)/3, then from the constitutive relations in (6), λ(1)=−2μ(1)/3. Thus, the field equations and relations for describing the fluid flow are complete.

Now we try to formulate the corresponding boundary conditions at the x2=0 level; the contact conditions between the metal and piezoelectric layers at the x2=−hm level; the compatibility conditions at the x2=−h interface level between the plate and the fluid; and finally the impermeability conditions on the rigid wall at x2=−h−hd.

The boundary conditions on the upper face plane of the plate are:

σ21(m)|x2=0=0,σ22(m)|x2=0=−P0δ(x1)eiωt,(10)

where δ(x1) is the Dirac delta function.

Perfect contact conditions on the plane between the metal and piezoelectric layers are:

σ21(m)|x2=−hm=σ21(p)|x2=−hm,σ22(m)|x2=−hm=σ22(p)|x2=−hm,u1(m)|x2=−hm=u1(p)|x2=−hm,u2(m)|x2=−hm=u2(p)|x2=−hm.(11)

Compatibility conditions on the interface plane between the fluid and plate are:

∂u1(p)∂t|x2=−h=V1|x2=−h,∂u2(p)∂t|x2=−h=V2|x2=−h,σ21(p)|x2=−h=T21|x2=−h,σ22(p)|x2=−h=T22|x2=−h.(12)

The impermeability conditions on the rigid wall are:

V1|x2=−h−hd=0,V2|x2=−h−hd=0.(13)

We also add the boundary conditions on the planes x2=−hm and x2=−h with respect to the electrical quantities of the piezoelectric layer to the foregoing ones. In the present work, we assume that the following open-circuit conditions take place:

D2(p)|x2=−hm=0,D2(p)|x2=−h=0.(14)

This completes the formulation of the problems investigated in the present paper.

3  Method of Solution

The presence of piezoelectric components in composite materials complicates the analytical solution of the corresponding problems. Therefore, it sometimes seems necessary to apply numerical methods, as in the works of Aylikci et al. [45], Ray et al. [46] and in many other works listed therein, to solve the corresponding problems. However, in the present case, the exponential Fourier transform is successfully used to solve the corresponding boundary value problem in obtaining an analytical expression for the transforms of the sought values whose originals are found numerically.

According to the definition of the time-harmonic oscillation, let us represent all the sought values of the considered problem as (x1,x2,t)=g¯(x1,x2)eiωt. According to this statement, replacing the derivatives ∂(⋅)/∂t and ∂2(⋅)/∂t2 with iω(⋅¯) and −ω2(⋅¯), respectively in the foregoing equations, the corresponding equations are obtained with the boundary, contact, compatibility and impermeability conditions with respect to the amplitudes of the sought quantities. We omit the over-bar on the amplitudes and we employ the exponential Fourier transform for the solution to these equations.

fF(s,x2)=∫−∞+∞f(x1,x2)eisx1dx1.(15)

According to which, the amplitudes can be presented as follows:

u1(α)=12π∫−∞∞u1F(α)(s,x2)e−isx1ds,u2(α)=12π∫−∞∞u2F(α)(s,x2)e−isx1ds,σ11(α)=12π∫−∞∞σ11F(α)(s,x2)e−isx1ds,σ22(α)=12π∫−∞∞σ22F(α)(s,x2)e−isx1ds,σ12(α)=12π∫−∞∞σ12F(α)(s,x2)e−isx1ds,D1(p)=12π∫−∞∞D1F(p)(s,x2)e−isx1ds,D2(p)=12π∫−∞∞D2F(p)(s,x2)e−isx1ds,φ(p)=12π∫−∞∞φF(p)(s,x2)e−isx1ds,φ(1)=12π∫−∞∞φF(1)(s,x2)e−isx1ds,ψ(1)=12π∫−∞∞ψF(1)(s,x2)e−isx1ds,V1=12π∫−∞∞V1F(s,x2)e−isx1ds,V2=12π∫−∞∞V2F(s,x2)e−isx1ds,T11=12π∫−∞∞T11F(s,x2)e−isx1ds,T22=12π∫−∞∞T22F(s,x2)e−isx1ds,T12=12π∫−∞∞T12F(s,x2)e−isx1ds.(16)

Now, we consider separately the determination of the Fourier transforms of the quantities related to the metal and piezoelectric layers of the plate and related to the fluid.

3.1 Determination of the Fourier Transforms of the Quantities Related to the Elastic-Metal Plate

Substituting the expressions related to the metal layer given in (16) into the Eqs. (2) and (4), we obtain the following equations of motion for the metal layer with respect to the displacement terms:

Au1F(m)+Bdu2F(m)dx2+d2u1F(m)dx22=0,Du2F(m)+Bdu1F(m)dx2+Gd2u2F(m)dx22=0(17)

where

A=X2−s2(λ/μ+2),B=−is(λ/μ+1),

D=X2−s2,G=λ/μ+2,X2=ω2h2/c22,c2=μ/ρ(18)

and in (17) under x2, it must be understood that x2/h. Introducing the notation:

A0=AG−B2+DG,B0=ADG,k1=−A02+A024−B0,k2=−A02−A024−B0(19)

we can write the solution of the Eq. (17) as follows:

u2F(m)=Z1(m)ek1x2+Z2(m)e−k1x2+Z3(m)ek2x2+Z4(m)e−k2x2,

u1F(m)=Z1(m)a1ek1x2+Z2(m)a2e−k1x2+Z3(m)a3ek2x2+Z4(m)a4e−k2x2.(20)

where

a1=−D−Gk12Bk1,a2=−a1,a3=−D−Gk22Bk2,a4=−a3.(21)

Using the relations (2), (4) and (20), we obtain the expressions for the Fourier transforms of the stresses σ11F(m), σ12F(m), and σ22F(m), and to reduce the volume of the paper we do not give these expressions here.

3.2 Determination of the Fourier Transforms of the Quantities Related to the Piezoelectric Layer

Substituting the expressions in (16) into the foregoing equations and relations in (1), (3), and (4) and doing the corresponding mathematical manipulations, we obtain the following system of the ordinary differential equations for u1F(p), u2F(p) and φF(p):

((c(p))2s2−c~11(p))u1F(p)+d2u1F(p)d(sx2)2−i(1+c~13(p))du2F(p)d(sx2)−i(1+e~31(p))dφ~F(p)d(sx2)=0,−i(1+c~13(p))du1F(p)d(sx2)+((c(p))2s2−1)u2F(p)+c~33(p)d2u2F(p)d(sx2)2−φ~F(p)+e~33(p)d2φ~F(p)d(sx2)2=0,−i(1+e~31(p))du1F(p)d(sx2)+e~33(p)d2u2F(p)d(sx2)2−u2F(p)+ε~11(p)φ~F(p)−ε~33(p)dφ~F(p)d(sx2)=0.(22)

where

(c(p))2=ω2h2c44(p)/ρ(p),c~11(p)=c11(p)c44(p),c~13(p)=c13(p)c44(p),e~31(p)=e31(p)e15(p),c~33(p)=c33(p)c44(p),e~33(p)=e33(p)c44(p),φ~F(p)=e15(p)c44(p)φF(p),ε~11(p)=ε11(p)c44(p)(e15(p))2,ε~33(p)=ε33(p)c44(p)(e15(p))2.(23)

We recall that in Eq. (22) under x2, it must be understood that x2/h. To find the solution to the system of the ordinary differential equations in (22), we employ the well-known Euler method according to which, the particular solution to this system can be represented as follows:

u1Fpar(p)=iA(p)ebsx2,u2Fpar(p)=B(p)ebsx2,φ~Fpar(p)=C(p)ebsx2.(24)

In (24), the lower index par indicates the “particular” solution. Moreover, the constants A(p), B(p) and C(p) are unknown constants, and in (24) through the symbol b, the unknown constant is also noted and these unknowns will be determined in the solution procedure described below.

Thus, substituting the particular solutions in (24) into the system of equations in (22) and employing the usual procedures, we obtain the following system of homogeneous linear algebraic equations with respect to the unknown constants A(p), B(p) and C(p).

(b2+((c(p))2s2−c~11(p)))A(p)−(1+c~13(p))bB(p)−(1+e~31(p))bC(p)=0,(1+c~13(p))bA(p)+((c(p))2s2−1+c~33(p)b2)B(p)+(e~33(p)b2−1)C(p)=0,(1+e~31(p))bA(p)+(e~33(p)b2−1)B(p)+(ε~11(p)−ε~33(p)b2)C(p)=0.(25)

According to the requirement of the existence of the non-trivial solution of the system of equations in (25), we equate to zero the determinant of the coefficient matrices of this system and we obtain the following bi-cubic equation for determination of the unknown constant b.

b13+a4b12+a2b1+a0=0,(26)

where

b1=b2,a4=[c~33(p)ε~11(p)−ε~33(p)α22−ε~33(p)α11c~33(p)−e~33(p)(α21α13+α13α31)+α13α31c~33(p)−ε~33(p)α12α21−α11(e~33(p))2+2e~33(p)](−ε~33(p)c~33(p)−(e~33(p))2)−1,

a2=[ε~11(p)α22+ε~11(p)α11ε~33(p)−ε~33(p)α11α22+α21α13+α12α31+α13α31α22+ε~11(p)α12α21−2α11e~33(p)−1](−ε~33(p)c~33(p)−(e~33(p))2)−1,

a0=[ε~11(p)α11α22−α11](−ε~33(p)c~33(p)−(e~33(p))2)−1.(27)

In (27), the notation given below is used.

α11=(c(p))2s2−c~11(p),α12=1+c~13(p),α13=1+e~31(p),α21=1+c~13(p),α22=(c(p))2s2−1,α31=1+e~31(p).(28)

We employ Vieta’s trigonometric formula for cubic equations in order to find the roots of the Eq. (26), according to which, in the first step, the values of the following expressions must be calculated.

Q=(a42−3a2)/9,R=(2a43−9a4a2+27a0)/54.(29)

In the next step, the values of the expression in (30) are calculated.

S=Q3−R2(30)

According to the sign of S, the cases are distinguished. So if S>0, then the roots of Eq. (26) are determined through the expressions:

b11=−2(Q)12cosφ−a43,b12=−2(Q)12cos(φ+2π3)−a43,b13=−2(Q)12cos(φ−2π/3)−a43,φ=arccos(RQ3/2)/3.(31)

However, if S<0, then these roots are determined through the expressions:

b11=−2sgn(R)|Q|12coshφ−a43,b12=sgn(R)|Q|12coshφ−a43+i(3|Q|)12sinhφ,b13=sgn(R)|Q|12coshφ−a43−i(3|Q|)12sinhφ,φ=arccos⁡h(|R||Q|3/2)/3.(32)

In this way, we determine the roots of the characteristic Eq. (26) and we can write the six roots of the bi-cubic equation as follows:

b1=b11,b2=−b1,b3=b12,b4=−b3,b5=b13,b6=−b5.(33)

According to the expressions in (24) and to the solution technique used, the general solution to the system of ordinary differential equations in (22) is:

u1F=iA1eb1sx2+iA2eb2sx2+iA3eb3sx2+iA4eb4sx2+iA5eb5sx2+iA6eb6sx2,u2F=A1Y1eb1sx2+A2Y2eb2sx2+A3Y3eb3sx2+A4Y4eb4sx2+A5Y5eb5sx2+A6Y6eb6sx2,φF=A1Z1eb1sx2+A2Z2eb2sx2+A3Z3eb3sx2+A4Z4eb4sx2+A5Z5eb5sx2+A6Z6eb6sx2.(34)