# Three-Dimensional Semianalytical Solutions for Piezoelectric Laminates Subjected to Underwater Shocks

Authors: Xu Liang; Wenbin Lu; Ronghua Zhu; Changpeng Ye; Guohua Liu
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

In this study, a piezoelectric laminate is analyzed by applying the Laplace transform and its numerical inversion, Fourier transform, differential quadrature method (DQM), and state space method. Based on the modified variation principle for the piezoelectric laminate, the fundamental equations for dynamic problems are derived. The solutions for the displacement, stress, electric potential, and dielectric displacement are obtained using the proposed method. Durbin’s inversion method for the Laplace transform is adopted. Four boundary conditions are discussed through the DQM. The proposed method is validated by comparing the results with those of the finite element method (FEM). Moreover, this semianalytical method is further extended to describe the dynamic response of piezoelectric laminated plates subjected to underwater shocks by introducing Taylor’s fluid-structure interaction algorithm. Both air-backed and water-backed laminated plates are investigated, and the effect of the fluid is examined. In the time domain, the electric potential and displacements of sample points are calculated under four boundary conditions. The present method is shown to be accurate and can be a useful method to calculate the dynamic response of underwater laminated sensors.


---

## Body

## 1. Introduction

Piezoelectric materials are extensively used as sensing elements in various engineering fields, such as aerospace engineering, ocean engineering, civil engineering, and mechanical engineering. In addition, piezoelectric materials are one of the most important components in smart structures which can sense and drive changes in external circumstances and then be modified in response by computers and processors [1]. Since the Curie brothers first discovered the direct piezoelectric effect in single-crystal quartz in 1880 and Gabriel Lippmann discovered the converse piezoelectric effect in 1881, various types of piezoelectric materials have been found, including both natural materials and composite materials [2, 3].Composite laminates are a common smart structure. As shown in Figure1, the structure is composed of stacked layers of different materials; the individual layers are generally orthotropic or transversely isotropic. Owing to the various characteristics of the individual layers, a hybrid laminate may exhibit multiple functions. Piezoelectric laminates are usually piled as sensors and actuators that will react to or induce a displacement or electric potential [4, 5].Figure 1
Piezoelectric laminate model.Pioneering research on piezoelectric laminates has been conducted primarily using finite element models (FEMs) and the space state method. A FEM has stable convergence and its solution is effective; however, it requires complex grids and long computation times to produce results with high accuracy [6, 7]. Mukherjee and Saha [8] focused on large deformations of piezoelectric beams using a FEM based on the first-order shear deformation theory and the Newton–Raphson method. Garcia et al. [9] developed a layer-wise FEM for a piezoelectric plate using the Reissner mixed variational principle. Dash and Singh [10] employed a C0 isoparametric FEM to solve the nonlinear free vibration problem. Wankhade and Bajoria [11] considered the higher-order shear deformation theory to carry out a FEM of a piezoelectric beam under static and dynamic excitations.On the other hand, the state space method is also an efficient and accurate method. In view of the relationships among transverse vectors, it is suitable for laminated structures [12]. Different approaches have been investigated in the state space method for piezoelectric laminates. D’Ottavio and Kröplin [13] combined the Reissner variational statement and constitutive equations to simulate piezoelectric laminates. Tiersten [14] demonstrated Hamilton’s principle for linear piezoelectric media under different boundary conditions and discussed both homogeneous and inhomogeneous problems. In addition, Qing et al. [15] proposed a modified mixed Hellinger–Reissner (H–R) variational principle based on the work of Steele and Kim [16].The state space method is the primary focus of this study. The variational principle for the state space method is a complex differential equation including several variables, and thus some transformations should be employed. The Laplace transform and differential quadrature method (DQM) are the two main transformations used in this study.The Laplace transform is usually used to change a function in the time domain to one in a complex frequency field for dynamic problems [17]. For this reason, in recent studies, the Laplace transform has generally been adopted to analyze both statics and dynamics problems, and the numerical inversion of the Laplace transform has been studied. Zhao [18] proposed two algorithms that divide the integration interval into small subspaces with different integration steps. Durbin [19] combined a trigonometric series with an inversion series and presented rather smaller results. Qing et al. [15] subsequently made comparisons to the inverse algorithms proposed by both Zhao and Durbin and improved Zhao’s algorithms by using Subbotin splines. Wang et al. [6] improved the calculation efficiency with two novel algorithms and additionally compensated for the deviation node in Qing’s algorithms. Durbin’s inversion method is used in this study owing to its high efficiency, and the results are sufficiently precise.First proposed in 1971 by Bellman et al. [20] based on integral quadrature, DQM is another method for dealing with partial differential factors. Moreover, with DQM, complex boundary conditions can be analyzed. Extensive research has focused on the choice of grid points, boundary conditions, weight coefficients, and the convergence behavior. Several researchers [21–23] have eliminated the unreliable results from uniform grid points and offered various distribution forms that could satisfy higher orders and different simulations. Chen et al. [24] employed the Chebyshev polynomial to obtain the weighting coefficients of the matrix form. To date, DQM has been widely applied for composite structures. Du and Shu [23] chose DQM to deal with the clamped and simply supported boundary conditions of beams and plates. Pradhan and Murmu [25] used DQM to analyze the vibration of functionally graded material beams.Piezoelectric laminates are employed for numerous applications in the sea. It is important for underwater structures to consider the interaction between the water and the structure. There are few studies regarding underwater piezoelectric laminates, and thus the fluid-structure interaction (FSI) is also considered in this study. Taylor [26] proposed a one-dimensional model to analyze the FSI. It has been further verified under various conditions with different materials, fluid media, and structures. Schiffer and Tagarielli [27] examined the FSI response with various materials and structural combinations. Blom [28] investigated a piston problem and combined the one-dimensional equation and the Euler equation to propose an algorithm. Deshpande and Fleck [29, 30] studied sandwich beams subjected to underwater waves using a lumped parameter model and FEM.In this study, we discuss three-dimensional semianalytical solutions for piezoelectric laminates subjected to an underwater shock. The state space method is established based on the constitutive equations for piezoelectric materials. Both air-backed and water-backed laminated plates are investigated based on the effect of the FSI. Combined with the Fourier transform, DQM, and Laplace transform, the present model is analyzed under four different boundary conditions. The results are verified using finite element analysis software.


## 2. Mathematical Description

### 2.1. Fundamental Equations

A piezoelectric laminate consists of a series of layers of different materials, including piezoelectric materials. Consider a laminated plate comprising three layers of PZT/Gr-Epoxy/PZT materials, where PZT means lead zirconate titanate piezoelectric ceramics. As shown in Figure1, the laminate is subject to an underwater shock. Thus, the upper layer is always exposed to the water, while the lower layer may be exposed to either water or air.Briefly, the piezoelectric effect is demonstrated to change polarization when a mechanical stress is applied or to create a mechanical deformation with the introduction of additional voltage. In view of basic theory, the constitutive equation for piezoelectric materials can be written simply as the following form:(1)σp=CpqESq+epiEi,Di=eipSp+εijSEj,where σ, S, E, and D are the stress, strain, electric field, and electric displacement, respectively, and C, e, and ε are the elastic, piezoelectric, and dielectric coefficients, respectively. The subscripts p, q, i, and j denote different directions in the material coordinates (i, j = 1, 2, 3; p, q = 1, 2, 3, 4, 5, 6). The superscript E indicates parameters measured at a constant electric field, and the superscript S denotes those measured at a constant or zero strain field. The superscripts will be omitted in the subsequent derivation process.Assuming that the piezoelectric materials are orthotropic or have orthotropic symmetry relative to thex-y coordinate plane, constitutive equation (1) can be extended to the following matrix form:(2)σ1σ2σ3σ4σ5σ6D1D2D3=C11C12C1300000e31C12C22C2300000e32C13C23C3300000e33000C44000e2400000C550e150000000C660000000e150−ε1100000e24000−ε220e31e32e3300000−ε33S1S2S3S4S5S6−E1−E2−E3,in the Cartesian coordinate system; p, q = x, y, z, yz, xz, xy, and i, j = x, y, z. For the model of a piezoelectric laminate shown in Figure 1, the subsequent equations and derivations will be described in the Cartesian coordinate system.In terms of geometric relationships, the relationships between the strains and displacements can be described in the form of gradient equations:(3)Sx=αu,Sy=βv,Sz=γw,Syz=γv+βw,Sxz=γu+αw,Sxy=βu+αv,where u, v, and w denote the displacements in the three directions of the Cartesian axes; and α, β, and γ are partial differential operators of the three Cartesian axes with respect to the following variable: α=∂/∂x,β=∂/∂y,δ=∂/∂z.In light of the quasi-static Maxwell equations, the electric fields,Ei, and electric potential, ϕ, have similar relationships as the strains and displacements:(4)Ex=−αϕ,Ey=−βϕ,Ez=−γϕ.Two general designations are defined asP=σxz,σyz,σz,DzT and Q=u,v,w,ϕT for the two groups including four dual vectors, while the remaining vectors of the variables are given as P2=σx,σy,σxy,Dx,DyT, D1=Sxz,Syz,Sz,−EzT, and D2=Sx,Sy,Sxy,−Ex,−EyT. Thus, the subsequent derivations can be simplified with these five designations. Constitutive equation (2) can then be rewritten as follows:(5)D1P2=Φ11Φ12Φ21Φ22PD2,where Φ11, Φ12, Φ21, and Φ22 represent the transformations of the 9 × 9 matrix in equation (2) and Φ21=−Φ12T.Hence, equations (3) and (4) become the two following matrix forms:(6)D1=Q˙+G1Q,D2=G2Q,where Q˙=dQ/dz,G1=0000α0β000000000,andG2=α0000β00βα00000α000β.According to the Reissner mixed variational theorem, Reissner’s energy density function is given by the above matrices:(7)LMR=PTD1+P2TD2−∫D1TdP−∫D2TdP2.Owing to the special form of orthotropic constitutive equations, the secondary vectors in equation (7) are replaced, and the integral parts are derived to dot-product forms:(8)LMR=PTQ˙+PTG1+Φ21TG2Q+12QG2TΦ22G2−ΩQ−12PTΦ11P,where Ω=−ρ∂2/∂t20000ρ∂2/∂t20000ρ∂2/∂t200000 for dynamic problems in the time domain.In terms of the Reissner variational principle, the energy equation with an internal force is given by the following equation:(9)δΠ=δ∫0h∬ΓLMRdΩdz,where h indicates the thickness of the layer and Γ is the area of the layer.Therefore, based on equations (8) and (9), the state space method form can be deduced as follows:(10)ddzPQ=A11A12A21A22PQ,where A11=G1T+G2TΦ21, A12=G2TΦ22G2−Ω, A21=Φ11, and A22=−A11T.In addition, the in-plane vectors and rest vectors can be evaluated by the expression in the following equation:(11)σx,σy,τxy,Dx,Dy=Φ21P+Φ22G2Q.The partial differential operator ofz in equation (10) can be eliminated by transforming the coefficient matrix into an exponential matrix as follows:(12)PzQz=expA11A12A21A22z−zkPzkQzk,zk<z<zk+1,k=1,2,…,where k indicates the kth layer and zk represents the value of the thickness at the kth layer.Therefore, the state vectors of the piezoelectric laminate atz can be written as follows:(13)PzQz=expHz−zk∏j=k1expHzj−zj−1P0Q0,where(14)H=A11A12A21A22.


### 2.2. Fluid-Structure Interaction (FSI)

As shown in Figure1, the piezoelectric laminate is subject to an underwater explosion, resulting in shock waves on both the upper and lower surfaces of the laminate. In accordance with the one-dimensional FSI model proposed by Taylor [26] in 1963, the effect of the shock waves on the laminate can be calculated.The shock wave is simply described as a function of the time,t, and the vertical distance to the plane, z, which in this paper is equal to zero:(15)p0z,t=pmexp−t−z/cwtd,where cw denotes the wave velocity in water as a constant and pm and td are empirical formulas for the explosive charge weight, W (kg), and perpendicular distance from the upper surface of the laminate, S (m), respectively. TNT is used for the explosion in this study:(16)pm=52.16W1/3S1.13MPa,td=0.0965W1/3W1/3S−0.22ms.Under these conditions, the planes influenced by the shock wave include both the upper and lower surfaces of the laminate [29, 30]. The wave pressure on one plane consists mainly of the primary shock wave p0, reflected wave, pr1, and rarefaction wave, pr2:(17)pr1z,t=pmexp−t+z/cwtd,pr2z,t=−ρwcww˙z,t,where ρw is the density of water and w˙z,t=∂wz,t/∂t.Consequently, when the laminate is completely submerged underwater, the wave pressures on the upper and lower surfaces are given by the following:(18)px,y,h,t=p00,t+pr10,t+pr20,t=2pmexp−ttd−ρwcww˙x,y,h,t,px,y,0,t=pr20,t=−ρwcww˙x,y,0,t.In this paper, the interactions between the fluid and structure are divided into two conditions. As the wave pressure on the upper surface is caused directly by the underwater shock, the lower surface is exposed to both water and air, which here are referred to as water-backed and air-backed laminated plates. Because the density of air is far smaller than the density of water, the wave pressure on the air-backed surface is assumed to be zero:px,y,0,t=0.Hence, on the upper and lower surfaces of the laminate, the boundary conditions are confirmed as follows according to the FSI:(19)σxz,σyz,σz,Dzx,y,h,t=0,0,px,y,h,t,0.(20)σxz,σyz,σz,Dzx,y,0,t=0,0,px,y,0,t,0,where px,y,0,t=−ρwcww˙x,y,0,t for the water-backed plate and px,y,0,t=0 for the air-backed plate.


## 2.1. Fundamental Equations

A piezoelectric laminate consists of a series of layers of different materials, including piezoelectric materials. Consider a laminated plate comprising three layers of PZT/Gr-Epoxy/PZT materials, where PZT means lead zirconate titanate piezoelectric ceramics. As shown in Figure1, the laminate is subject to an underwater shock. Thus, the upper layer is always exposed to the water, while the lower layer may be exposed to either water or air.Briefly, the piezoelectric effect is demonstrated to change polarization when a mechanical stress is applied or to create a mechanical deformation with the introduction of additional voltage. In view of basic theory, the constitutive equation for piezoelectric materials can be written simply as the following form:(1)σp=CpqESq+epiEi,Di=eipSp+εijSEj,where σ, S, E, and D are the stress, strain, electric field, and electric displacement, respectively, and C, e, and ε are the elastic, piezoelectric, and dielectric coefficients, respectively. The subscripts p, q, i, and j denote different directions in the material coordinates (i, j = 1, 2, 3; p, q = 1, 2, 3, 4, 5, 6). The superscript E indicates parameters measured at a constant electric field, and the superscript S denotes those measured at a constant or zero strain field. The superscripts will be omitted in the subsequent derivation process.Assuming that the piezoelectric materials are orthotropic or have orthotropic symmetry relative to thex-y coordinate plane, constitutive equation (1) can be extended to the following matrix form:(2)σ1σ2σ3σ4σ5σ6D1D2D3=C11C12C1300000e31C12C22C2300000e32C13C23C3300000e33000C44000e2400000C550e150000000C660000000e150−ε1100000e24000−ε220e31e32e3300000−ε33S1S2S3S4S5S6−E1−E2−E3,in the Cartesian coordinate system; p, q = x, y, z, yz, xz, xy, and i, j = x, y, z. For the model of a piezoelectric laminate shown in Figure 1, the subsequent equations and derivations will be described in the Cartesian coordinate system.In terms of geometric relationships, the relationships between the strains and displacements can be described in the form of gradient equations:(3)Sx=αu,Sy=βv,Sz=γw,Syz=γv+βw,Sxz=γu+αw,Sxy=βu+αv,where u, v, and w denote the displacements in the three directions of the Cartesian axes; and α, β, and γ are partial differential operators of the three Cartesian axes with respect to the following variable: α=∂/∂x,β=∂/∂y,δ=∂/∂z.In light of the quasi-static Maxwell equations, the electric fields,Ei, and electric potential, ϕ, have similar relationships as the strains and displacements:(4)Ex=−αϕ,Ey=−βϕ,Ez=−γϕ.Two general designations are defined asP=σxz,σyz,σz,DzT and Q=u,v,w,ϕT for the two groups including four dual vectors, while the remaining vectors of the variables are given as P2=σx,σy,σxy,Dx,DyT, D1=Sxz,Syz,Sz,−EzT, and D2=Sx,Sy,Sxy,−Ex,−EyT. Thus, the subsequent derivations can be simplified with these five designations. Constitutive equation (2) can then be rewritten as follows:(5)D1P2=Φ11Φ12Φ21Φ22PD2,where Φ11, Φ12, Φ21, and Φ22 represent the transformations of the 9 × 9 matrix in equation (2) and Φ21=−Φ12T.Hence, equations (3) and (4) become the two following matrix forms:(6)D1=Q˙+G1Q,D2=G2Q,where Q˙=dQ/dz,G1=0000α0β000000000,andG2=α0000β00βα00000α000β.According to the Reissner mixed variational theorem, Reissner’s energy density function is given by the above matrices:(7)LMR=PTD1+P2TD2−∫D1TdP−∫D2TdP2.Owing to the special form of orthotropic constitutive equations, the secondary vectors in equation (7) are replaced, and the integral parts are derived to dot-product forms:(8)LMR=PTQ˙+PTG1+Φ21TG2Q+12QG2TΦ22G2−ΩQ−12PTΦ11P,where Ω=−ρ∂2/∂t20000ρ∂2/∂t20000ρ∂2/∂t200000 for dynamic problems in the time domain.In terms of the Reissner variational principle, the energy equation with an internal force is given by the following equation:(9)δΠ=δ∫0h∬ΓLMRdΩdz,where h indicates the thickness of the layer and Γ is the area of the layer.Therefore, based on equations (8) and (9), the state space method form can be deduced as follows:(10)ddzPQ=A11A12A21A22PQ,where A11=G1T+G2TΦ21, A12=G2TΦ22G2−Ω, A21=Φ11, and A22=−A11T.In addition, the in-plane vectors and rest vectors can be evaluated by the expression in the following equation:(11)σx,σy,τxy,Dx,Dy=Φ21P+Φ22G2Q.The partial differential operator ofz in equation (10) can be eliminated by transforming the coefficient matrix into an exponential matrix as follows:(12)PzQz=expA11A12A21A22z−zkPzkQzk,zk<z<zk+1,k=1,2,…,where k indicates the kth layer and zk represents the value of the thickness at the kth layer.Therefore, the state vectors of the piezoelectric laminate atz can be written as follows:(13)PzQz=expHz−zk∏j=k1expHzj−zj−1P0Q0,where(14)H=A11A12A21A22.


## 2.2. Fluid-Structure Interaction (FSI)

As shown in Figure1, the piezoelectric laminate is subject to an underwater explosion, resulting in shock waves on both the upper and lower surfaces of the laminate. In accordance with the one-dimensional FSI model proposed by Taylor [26] in 1963, the effect of the shock waves on the laminate can be calculated.The shock wave is simply described as a function of the time,t, and the vertical distance to the plane, z, which in this paper is equal to zero:(15)p0z,t=pmexp−t−z/cwtd,where cw denotes the wave velocity in water as a constant and pm and td are empirical formulas for the explosive charge weight, W (kg), and perpendicular distance from the upper surface of the laminate, S (m), respectively. TNT is used for the explosion in this study:(16)pm=52.16W1/3S1.13MPa,td=0.0965W1/3W1/3S−0.22ms.Under these conditions, the planes influenced by the shock wave include both the upper and lower surfaces of the laminate [29, 30]. The wave pressure on one plane consists mainly of the primary shock wave p0, reflected wave, pr1, and rarefaction wave, pr2:(17)pr1z,t=pmexp−t+z/cwtd,pr2z,t=−ρwcww˙z,t,where ρw is the density of water and w˙z,t=∂wz,t/∂t.Consequently, when the laminate is completely submerged underwater, the wave pressures on the upper and lower surfaces are given by the following:(18)px,y,h,t=p00,t+pr10,t+pr20,t=2pmexp−ttd−ρwcww˙x,y,h,t,px,y,0,t=pr20,t=−ρwcww˙x,y,0,t.In this paper, the interactions between the fluid and structure are divided into two conditions. As the wave pressure on the upper surface is caused directly by the underwater shock, the lower surface is exposed to both water and air, which here are referred to as water-backed and air-backed laminated plates. Because the density of air is far smaller than the density of water, the wave pressure on the air-backed surface is assumed to be zero:px,y,0,t=0.Hence, on the upper and lower surfaces of the laminate, the boundary conditions are confirmed as follows according to the FSI:(19)σxz,σyz,σz,Dzx,y,h,t=0,0,px,y,h,t,0.(20)σxz,σyz,σz,Dzx,y,0,t=0,0,px,y,0,t,0,where px,y,0,t=−ρwcww˙x,y,0,t for the water-backed plate and px,y,0,t=0 for the air-backed plate.


## 3. Boundary Conditions and Transform Methods

### 3.1. Laplace Transform Methods

The Laplace transform method is adopted to obtain the response in the time domain, which is described in the two following equations:(21)f˜s=ℒft=∫0∞ftexp−stdt,(22)ft=ℒ−1f˜s=12πi∫ξ−i∞ξ+i∞f˜sexpstds.In this paper, Durbin’s inversion method given by equation (23) is adopted. The trigonometric series is obtained, and the error of the results is independent of t, and thus the analytical model can yield good results with high efficiency:(23)ft=2eαtTf12+∑k=1KRefk+1coskπtT,where sk=α+kπ/Ti,k=0,1,…,K, fk+1=fsk, and α=5/T is the real part of imaginary number s, and T=2T0, where T0 denotes the duration of the model in the time domain. Here, the value of K is 100.


### 3.2. Fourier Transform Methods

To represent different means of support, the boundary conditions atx = 0, a are expressed in four ways, as given by equations (24)–(32).Simply supported (x = 0)-simply supported (x = a) (S-S):(24)atx=0,σx=v=w=0,(25)atx=a,σx=v=w=0.Clamped (x = 0)-clamped (x = a) (C-C):(26)atx=0,u=v=w=0,(27)atx=a,u=v=w=0.Clamped (x = 0)-free (x = a) (C-F):(28)atx=0,u=v=w=0,(29)atx=a,σx=σxy=σxz=0.Clamped (x = 0)-simply supported (x = a) (C-S):(30)atx=0,u=v=w=0,(31)atx=a,σx=v=w=0.The Fourier transform method is employed here to change the form of∂/∂y. According to the coefficient matrix in equation (10), the Fourier transform indicates the following:(32)σxzσyzσzDzuvwϕ=∑n=1∞σ¯xzsinnπy/bσ¯yzcosnπy/bσ¯zsinnπy/bD¯zsinnπy/bu¯sinnπy/bv¯cosnπy/bw¯sinnπy/bϕ¯sinnπy/b,n=1,2,3,….In addition, the Fourier transform of the pressure on the surfaces yields(33)p=∑n=1∞p¯sinnyπb,p¯=2b∫0bpsinnyπbdy=2η2pb,where(34)η2=−1+−1nnπ.


### 3.3. Differential Quadrature Method (DQM)

While the Fourier transform method is used to change the form of∂/∂y, DQM is used to parse the form of ∂/∂x.By selecting some grid points in the domain ofx, the differential form of a continuous function, fx, can be described by an approximate function that is calculated with the values of all grid points and a weighting coefficient matrix [7, 31, 32]:(35)∂ifx∂xi=∑n=1MBmnifxn,m=1,2,…,M;n=1,2,…,M−1,where the superscript i indicates the ith-order derivative, M is the total number of sampling points, and Bmni denotes the weighting coefficients.The Chebyshev–Gauss–Lobatto points are chosen for their efficient speed of convergence and high computational accuracy [7]. Therefore, the sample points of the x-coordinate can be chosen as follows:(36)xi=121−cosi−1πN−1,i=1,2,…,N.


### 3.4. Normalization of the Variables

The variables need to be normalized to prevent nonconvergence of the computation caused by large gaps in magnitude between different parameters.According to the constitutive equation in equation (1), the variables can be normalized in the following ways:(37)σ⟶pq=σpqC33,C⟶pq=CpqC33,e⟶pi=epiC33ε33,ε⟶ij=εijε33,D⟶i=DiC33ε33,ϕ⟶=ϕaε33C33,t⟶=tca,s⟶=shc,u⟶i=uia,ρ⟶=ρρ0,x⟶,y⟶,z⟶=x,y,za,where c=C33/ρ0.For FSI, the density of water,ρw, and wave velocity underwater, cw, are normalized as follows:(38)ρ⟶w=ρwρ0,c⟶w=cwc.The⟶ of the normalized variables will be omitted in the following equations.In view of the semianalytical transformations mentioned above, for the boundary conditions of C-C, S-S, C-F, and C-S, the state space method form in equation (10) is further deduced as follows:(39)dσ˜¯xzidz=k8∑j=1MAij1σ˜¯zj+k10∑j=1MAij1D˜¯zj−k12∑j=1MAij2u˜¯j+ρs2+k15ζ2u˜¯i+k13+k15ζ∑j=1MAij2w˜¯j,dσ˜¯yzidz=k9ζσ˜¯zi+k11ζD˜¯zmi−k13+k15ζ∑j=1MAij1u˜¯j+ρs2+k14ζ2v˜¯i−k15∑j=1MAij2v˜¯j,dσ˜¯zidz=−∑j=1MAij1σ˜¯xzj+ζσ˜¯yzi+ρs2w˜¯i,dD˜¯zidz=k6∑j=1MAij1σ˜¯xzj−k7ζσ˜¯yzi−k16∑j=1MAij2ϕ˜¯j+k17ζ2ϕ˜¯i,du˜¯idz=k1σ˜¯xzi−∑j=1MAij1w˜¯j+k6∑j=1MAij1ϕ˜¯j,dv˜¯idz=k2σ˜¯yzi−ζw˜¯i+k7ζϕ˜¯i,dw˜¯idz=k3σ˜¯zi+k4D˜¯zi+k8∑j=1MAij1u˜¯j−k9ζv˜¯i,dϕ˜¯idz=k4σ˜¯zi+k5D˜¯zi+k10∑j=1MAij1u˜¯j−k11ζv˜¯i,where the subscript i indicates the ith grid point in the DQM (i = 1, 2, …, M), Aij1=ηBij1, Aij2=η2Bij2, and ζ=nhπ/b,η=h/a.Due to DQM, the vectors of the four support means as equations (26)–(33) can be rewritten as(40)σx1,σxy1,σxz1,u1,v1,w1,atx=0,σxM,σxyM,σxzM,uM,vM,wM,atx=a.Based on equation (11) and the boundary conditions as stated above, equation (39) should be deduced to avoid singular matrices in it, which means the coefficient matrix H=A11A12A21A22 in equation (10) should be deduced into different forms for the four boundary conditions:(41)C−C:H=00k8Amn1H14H15H160000k9ςI2H24H25H2600−Amn1ςI20000H37H38H41H420000H47H48k1I200000−Amn1H580k2I20000−ςI2H6800k3I2H74kgAmn1−k9ςI20000H83H84H85H8600,where m, n = 2∼M − 1 in the matrix and(42)H41=k6Amn1,m=1∼M,n=2∼M−1,H42=k7ζ0M−2×1,I2,0M−2×1T,H83=k40M−2×1,I2,0M−2×1T,H84=k5IM−k42k3EM1+EM−1,H14=k10Amn1−k4k8k3Am11,0M−2×M−2,k4k8k3AmM1,m=2∼M−1,n=1∼M,H24=k11ζ0M−2×1,I2,0M−2×1,H74=k40M−2×1,I2,0M−2×1,H15=ρs2+k15ζ2I2−k12Amn2−k82k3A1n1Am11+AMn1AmM1,m,n=2∼M−1,H25=−k13+k15ζAmn1,m,n=2∼M−1,H58=k6Amn1,m=2∼M−1,n=1∼M,H85=k10Amn1−k4k8k3An11,0M−1×M−2,k4k8k3AnM1T,m=1∼M,n=2∼M−1,H85=k10Amn1−k4k8k3An11,0M−1×M−2,k4k8k3AnM1T,m=1∼M,n=2∼M−1,H16=k13+k15ζAmn1,m,n=2∼M−1,H26=ρs2+ζ2I2−k15Amn2,m,n=2∼M−1,H86=−k11ζ0M−2×1,I2,0M−2×1T,H37=ρs2I2−A1n1Am11k1−AMn1AmM1k1,m,n=2∼M−1,H47=k6k1A1n1Am11+AMn1AmM1,m=1∼M,n=2∼M−1,H68=k7ζ0M−2×1,I2,0M−2×1,H38=k6k1A1n1Am11+AMn1AmM1,m=2∼M−1,n=1∼M,H48=k17ζ2IM−k16Amn2−k62k1A1n1Am11+AMn1AmM1+k72k2EM1+EM−1,m,n=1∼M,C−F:H=00k8Amn1k10Amn1H15H160000H23H24H25H2600H31H320000H37H38H41H420000H47H48k1I200000H57H580k2I10000−ςI1H6800H73H74H75H760000H83H84H85H8600,where m = 2∼M − 1, n = 1∼M − 1 in the matrix, and(43)H31=−Amn1,m=1∼M−1,n=2∼M−1,H41=k6Amn1,m=1∼M,n=2∼M−1,H32=ζ0M−2×1,I2T,0M−1×1,H42=−k7ζ0M−1×1,I1T,H23=k9ζ0M−1×1,I2,0M−2×1T,H73=k30M−1×1,I2,0M−2×1T,H24=k11ζ0M−1×1,I1−k10k9ζk80M−1×1,EM−1−1,H83=k4I1,0M−1×1T,H74=k40M−1×1,I1−k10k3k80M−1×1,EM−1−1,H84=k5IM−k10k4k8EM−1,H15=k12AmM1AMn1−k12Amn2+ρs2+k15ζ2I2,m,n=2∼M−1,H25=−k13+k15ζAmn1,m,n=2∼M−1,H75=k8Amn1+0M−2×M−2,k12k3AnM1k8T,m=2∼M,n=2∼M−1,H85=k10Amn1+0M−2×M−1,k12k4AnM1k8T,m=1∼M,n=2∼M−1,H16=k15+k13ζAmn1+k13ζ0M−2×M−2,AmM1−k12AmM1AMn1AMM1−1,m=2∼M−1,n=2∼M,H26=ρs2+k14ζ2I1−k15Amn2+k13+k15AmM1AMn1−0M−1×M−2,k12k9AnM1AMM1k8T,m,n=2∼M,H76=k13k3ζk8EM−1−1−k8ζAmM1AMn1−k9ζI1−0M−1×M−2,k12k3AnM1AMM1k8ζT,m,n=2∼M−1,H86=k13k4ζk80M−1×1,EM−1−1T−k11ζ0M−1×1,IM−1T−k10ζAmM1AMn1−0M−1×M−1,k12k4AnM1AMM1k8ζT,m=1∼M,n=2∼M,H37=ρs20M−2×1,I2T,0M−1×1−Am11A1n1k1,m=1∼M−1,n=2∼M,H47=k6k1Am11A1n1,m=1∼M,n=2∼M,H57=−Amn1,m=2∼M−1,n=2∼M,H38=k6k1Am11A1n1−k7ζ2k2EM−11,0M−1×1,m=1∼M−1,n=1∼M,H48=k17ζ2IM−k16Amn2−k62k1Am11A1n1+k72ζ2k2EM1,m,n=1∼M,H58=k6Amn1,m=2∼M−1,n=1∼M,H68=k7ζ0M−1×1,I1,C−S:H=00H13H14H15H160000k9ςI2H24H25H2600H31ςI20000H37H38k6Amn1H420000H47H48k1I100000H57H580k2I20000−ςI2H6800k3I2H74H75−k9ςI20000H83H84H85H8600,where m, n = 1∼M in the matrix and(44)H31=−Amn1,m=2∼M−1,n=2∼M,H42=−k7ζ0M−2×1,I2,0M−2×1T,H13=k8Amn1,m=2∼M,n=2∼M−1,H83=k40M−2×1,I2,0M−2×1T,H14=k10Amn1,0M−1×1−k4k8k3Am11,0M−1×M−1,m=2∼M,n=1∼M−1,H24=k11ζ0M−2×1,I1,0M−2×1T,H74=k40M−2×1,I2,0M−2×1,H84=k5IM−k42k3EM1−k4k10k8EM−1,H25=−k13+k15ζAmn1,m=2∼M−1,n=2∼M,H15=−k12Amn2+k12AmM1AMn1−k82k3Am11A1n1+ρs2+k15ζ2I1,m,n=2∼M,H75=k8Amn1,m=2∼M−1,n=2∼M,H16=k13+k15ζAmn1,m=2∼M,n=2∼M−1,H85=k10Amn1−k4k8k3An11,0M−1×M−2,k4k12k8AnM1T,m=1∼M,n=2∼M,H26=−k15Amn2+ρs2+k14ζ2I2,m,n=2∼M−1,H86=−k11ζ0M−2×1,I2,0M−2×1T,H37=ρs2I2−A1n1Am11k1,m,n=2∼M−1,H47=k6k1A1n1Am11,m=1∼M,n=2∼M−1,H57=−Amn1,m=2∼M,n=2∼M−1,H38=k6k1A1n1Am11,m=2∼M−1,n=1∼M,H48=k17ζ2IM−k16Amn2−k62k1A1n1Am11+k72k2EM1+EM−1,m,n=1∼M,H58=k6Amn1,m=2∼M,n=1∼M,H68=k7ζ0M−2×1,I2,0M−2×1,S−S:H=00k8Amn1H14H15H160000k9ςI2H24H25H2600H31ςI20000ρs2I20k6Amn1H4200000H48k1IM00000−Amn1H580H620000−ςI2H6800k3I2H74H75−k9ςI20000H83H84H85H8600,where(45)H31=−Amn1,m=2∼M−1,n=1∼M,H42=k7ζ0M−2×1,I2,0M−2×1T,H62=k20M−2×1,I2,0M−2×1,H83=k40M−2×1,I2,0M−2×1T,H14=0M×1,k10Amn1,0M×1,m=2∼M−1,n=1∼M,H24=k11ζ0M−2×1,I2,0M−2×1,H74=k40M−2×1,I2,0M−2×1,H84=k5IM−k10k4k8EM1+EM−1,H86=−k11ζ0M−2×1,I2,0M−2×1T,H15=k12Am11A1n1−Amn2+AmM1AMM1+ρs2+k15ζ2IM,m,n=1∼M,H25=−k13+k15ζAmn1,m=2∼M−1,n=1∼M,H75=k8Amn1,m=2∼M−1,n=1∼M,H16=k15+k13ζAmn1,m=1∼M,n=2∼M−1,H85=k10Amn1+k12k4k3An11,0M×M−2,k4k12k8AnM1T,m,n=1∼M,H26=ρs2+k14ζ2I2−k15Amn2,m,n=2∼M−1,H48=k17ζ20M−1×1,I1T−k16Amn2+k72ζ2k2EM−1,m=1∼M,n=2∼M,H58=k6Amn1,m=1∼M,n=2∼M,H68=k7ζI2,0M−2×1,where(46)IM=10…001…0⋮⋮⋱⋮00…1M×M,I1=10…001…0⋮⋮⋱⋮00…1M−1×M−1,I2=10…001…0⋮⋮⋱⋮00…1M−2×M−2,I=1,0m×n=0⋯0⋮⋱⋮0⋯0m×n,Em1=10…000…0⋮⋮⋱⋮00…0m×m,Em−1=0…00⋮⋱⋮⋮0…000…01m×m.


## 3.1. Laplace Transform Methods

The Laplace transform method is adopted to obtain the response in the time domain, which is described in the two following equations:(21)f˜s=ℒft=∫0∞ftexp−stdt,(22)ft=ℒ−1f˜s=12πi∫ξ−i∞ξ+i∞f˜sexpstds.In this paper, Durbin’s inversion method given by equation (23) is adopted. The trigonometric series is obtained, and the error of the results is independent of t, and thus the analytical model can yield good results with high efficiency:(23)ft=2eαtTf12+∑k=1KRefk+1coskπtT,where sk=α+kπ/Ti,k=0,1,…,K, fk+1=fsk, and α=5/T is the real part of imaginary number s, and T=2T0, where T0 denotes the duration of the model in the time domain. Here, the value of K is 100.


## 3.2. Fourier Transform Methods

To represent different means of support, the boundary conditions atx = 0, a are expressed in four ways, as given by equations (24)–(32).Simply supported (x = 0)-simply supported (x = a) (S-S):(24)atx=0,σx=v=w=0,(25)atx=a,σx=v=w=0.Clamped (x = 0)-clamped (x = a) (C-C):(26)atx=0,u=v=w=0,(27)atx=a,u=v=w=0.Clamped (x = 0)-free (x = a) (C-F):(28)atx=0,u=v=w=0,(29)atx=a,σx=σxy=σxz=0.Clamped (x = 0)-simply supported (x = a) (C-S):(30)atx=0,u=v=w=0,(31)atx=a,σx=v=w=0.The Fourier transform method is employed here to change the form of∂/∂y. According to the coefficient matrix in equation (10), the Fourier transform indicates the following:(32)σxzσyzσzDzuvwϕ=∑n=1∞σ¯xzsinnπy/bσ¯yzcosnπy/bσ¯zsinnπy/bD¯zsinnπy/bu¯sinnπy/bv¯cosnπy/bw¯sinnπy/bϕ¯sinnπy/b,n=1,2,3,….In addition, the Fourier transform of the pressure on the surfaces yields(33)p=∑n=1∞p¯sinnyπb,p¯=2b∫0bpsinnyπbdy=2η2pb,where(34)η2=−1+−1nnπ.


## 3.3. Differential Quadrature Method (DQM)

While the Fourier transform method is used to change the form of∂/∂y, DQM is used to parse the form of ∂/∂x.By selecting some grid points in the domain ofx, the differential form of a continuous function, fx, can be described by an approximate function that is calculated with the values of all grid points and a weighting coefficient matrix [7, 31, 32]:(35)∂ifx∂xi=∑n=1MBmnifxn,m=1,2,…,M;n=1,2,…,M−1,where the superscript i indicates the ith-order derivative, M is the total number of sampling points, and Bmni denotes the weighting coefficients.The Chebyshev–Gauss–Lobatto points are chosen for their efficient speed of convergence and high computational accuracy [7]. Therefore, the sample points of the x-coordinate can be chosen as follows:(36)xi=121−cosi−1πN−1,i=1,2,…,N.


## 3.4. Normalization of the Variables

The variables need to be normalized to prevent nonconvergence of the computation caused by large gaps in magnitude between different parameters.According to the constitutive equation in equation (1), the variables can be normalized in the following ways:(37)σ⟶pq=σpqC33,C⟶pq=CpqC33,e⟶pi=epiC33ε33,ε⟶ij=εijε33,D⟶i=DiC33ε33,ϕ⟶=ϕaε33C33,t⟶=tca,s⟶=shc,u⟶i=uia,ρ⟶=ρρ0,x⟶,y⟶,z⟶=x,y,za,where c=C33/ρ0.For FSI, the density of water,ρw, and wave velocity underwater, cw, are normalized as follows:(38)ρ⟶w=ρwρ0,c⟶w=cwc.The⟶ of the normalized variables will be omitted in the following equations.In view of the semianalytical transformations mentioned above, for the boundary conditions of C-C, S-S, C-F, and C-S, the state space method form in equation (10) is further deduced as follows:(39)dσ˜¯xzidz=k8∑j=1MAij1σ˜¯zj+k10∑j=1MAij1D˜¯zj−k12∑j=1MAij2u˜¯j+ρs2+k15ζ2u˜¯i+k13+k15ζ∑j=1MAij2w˜¯j,dσ˜¯yzidz=k9ζσ˜¯zi+k11ζD˜¯zmi−k13+k15ζ∑j=1MAij1u˜¯j+ρs2+k14ζ2v˜¯i−k15∑j=1MAij2v˜¯j,dσ˜¯zidz=−∑j=1MAij1σ˜¯xzj+ζσ˜¯yzi+ρs2w˜¯i,dD˜¯zidz=k6∑j=1MAij1σ˜¯xzj−k7ζσ˜¯yzi−k16∑j=1MAij2ϕ˜¯j+k17ζ2ϕ˜¯i,du˜¯idz=k1σ˜¯xzi−∑j=1MAij1w˜¯j+k6∑j=1MAij1ϕ˜¯j,dv˜¯idz=k2σ˜¯yzi−ζw˜¯i+k7ζϕ˜¯i,dw˜¯idz=k3σ˜¯zi+k4D˜¯zi+k8∑j=1MAij1u˜¯j−k9ζv˜¯i,dϕ˜¯idz=k4σ˜¯zi+k5D˜¯zi+k10∑j=1MAij1u˜¯j−k11ζv˜¯i,where the subscript i indicates the ith grid point in the DQM (i = 1, 2, …, M), Aij1=ηBij1, Aij2=η2Bij2, and ζ=nhπ/b,η=h/a.Due to DQM, the vectors of the four support means as equations (26)–(33) can be rewritten as(40)σx1,σxy1,σxz1,u1,v1,w1,atx=0,σxM,σxyM,σxzM,uM,vM,wM,atx=a.Based on equation (11) and the boundary conditions as stated above, equation (39) should be deduced to avoid singular matrices in it, which means the coefficient matrix H=A11A12A21A22 in equation (10) should be deduced into different forms for the four boundary conditions:(41)C−C:H=00k8Amn1H14H15H160000k9ςI2H24H25H2600−Amn1ςI20000H37H38H41H420000H47H48k1I200000−Amn1H580k2I20000−ςI2H6800k3I2H74kgAmn1−k9ςI20000H83H84H85H8600,where m, n = 2∼M − 1 in the matrix and(42)H41=k6Amn1,m=1∼M,n=2∼M−1,H42=k7ζ0M−2×1,I2,0M−2×1T,H83=k40M−2×1,I2,0M−2×1T,H84=k5IM−k42k3EM1+EM−1,H14=k10Amn1−k4k8k3Am11,0M−2×M−2,k4k8k3AmM1,m=2∼M−1,n=1∼M,H24=k11ζ0M−2×1,I2,0M−2×1,H74=k40M−2×1,I2,0M−2×1,H15=ρs2+k15ζ2I2−k12Amn2−k82k3A1n1Am11+AMn1AmM1,m,n=2∼M−1,H25=−k13+k15ζAmn1,m,n=2∼M−1,H58=k6Amn1,m=2∼M−1,n=1∼M,H85=k10Amn1−k4k8k3An11,0M−1×M−2,k4k8k3AnM1T,m=1∼M,n=2∼M−1,H85=k10Amn1−k4k8k3An11,0M−1×M−2,k4k8k3AnM1T,m=1∼M,n=2∼M−1,H16=k13+k15ζAmn1,m,n=2∼M−1,H26=ρs2+ζ2I2−k15Amn2,m,n=2∼M−1,H86=−k11ζ0M−2×1,I2,0M−2×1T,H37=ρs2I2−A1n1Am11k1−AMn1AmM1k1,m,n=2∼M−1,H47=k6k1A1n1Am11+AMn1AmM1,m=1∼M,n=2∼M−1,H68=k7ζ0M−2×1,I2,0M−2×1,H38=k6k1A1n1Am11+AMn1AmM1,m=2∼M−1,n=1∼M,H48=k17ζ2IM−k16Amn2−k62k1A1n1Am11+AMn1AmM1+k72k2EM1+EM−1,m,n=1∼M,C−F:H=00k8Amn1k10Amn1H15H160000H23H24H25H2600H31H320000H37H38H41H420000H47H48k1I200000H57H580k2I10000−ςI1H6800H73H74H75H760000H83H84H85H8600,where m = 2∼M − 1, n = 1∼M − 1 in the matrix, and(43)H31=−Amn1,m=1∼M−1,n=2∼M−1,H41=k6Amn1,m=1∼M,n=2∼M−1,H32=ζ0M−2×1,I2T,0M−1×1,H42=−k7ζ0M−1×1,I1T,H23=k9ζ0M−1×1,I2,0M−2×1T,H73=k30M−1×1,I2,0M−2×1T,H24=k11ζ0M−1×1,I1−k10k9ζk80M−1×1,EM−1−1,H83=k4I1,0M−1×1T,H74=k40M−1×1,I1−k10k3k80M−1×1,EM−1−1,H84=k5IM−k10k4k8EM−1,H15=k12AmM1AMn1−k12Amn2+ρs2+k15ζ2I2,m,n=2∼M−1,H25=−k13+k15ζAmn1,m,n=2∼M−1,H75=k8Amn1+0M−2×M−2,k12k3AnM1k8T,m=2∼M,n=2∼M−1,H85=k10Amn1+0M−2×M−1,k12k4AnM1k8T,m=1∼M,n=2∼M−1,H16=k15+k13ζAmn1+k13ζ0M−2×M−2,AmM1−k12AmM1AMn1AMM1−1,m=2∼M−1,n=2∼M,H26=ρs2+k14ζ2I1−k15Amn2+k13+k15AmM1AMn1−0M−1×M−2,k12k9AnM1AMM1k8T,m,n=2∼M,H76=k13k3ζk8EM−1−1−k8ζAmM1AMn1−k9ζI1−0M−1×M−2,k12k3AnM1AMM1k8ζT,m,n=2∼M−1,H86=k13k4ζk80M−1×1,EM−1−1T−k11ζ0M−1×1,IM−1T−k10ζAmM1AMn1−0M−1×M−1,k12k4AnM1AMM1k8ζT,m=1∼M,n=2∼M,H37=ρs20M−2×1,I2T,0M−1×1−Am11A1n1k1,m=1∼M−1,n=2∼M,H47=k6k1Am11A1n1,m=1∼M,n=2∼M,H57=−Amn1,m=2∼M−1,n=2∼M,H38=k6k1Am11A1n1−k7ζ2k2EM−11,0M−1×1,m=1∼M−1,n=1∼M,H48=k17ζ2IM−k16Amn2−k62k1Am11A1n1+k72ζ2k2EM1,m,n=1∼M,H58=k6Amn1,m=2∼M−1,n=1∼M,H68=k7ζ0M−1×1,I1,C−S:H=00H13H14H15H160000k9ςI2H24H25H2600H31ςI20000H37H38k6Amn1H420000H47H48k1I100000H57H580k2I20000−ςI2H6800k3I2H74H75−k9ςI20000H83H84H85H8600,where m, n = 1∼M in the matrix and(44)H31=−Amn1,m=2∼M−1,n=2∼M,H42=−k7ζ0M−2×1,I2,0M−2×1T,H13=k8Amn1,m=2∼M,n=2∼M−1,H83=k40M−2×1,I2,0M−2×1T,H14=k10Amn1,0M−1×1−k4k8k3Am11,0M−1×M−1,m=2∼M,n=1∼M−1,H24=k11ζ0M−2×1,I1,0M−2×1T,H74=k40M−2×1,I2,0M−2×1,H84=k5IM−k42k3EM1−k4k10k8EM−1,H25=−k13+k15ζAmn1,m=2∼M−1,n=2∼M,H15=−k12Amn2+k12AmM1AMn1−k82k3Am11A1n1+ρs2+k15ζ2I1,m,n=2∼M,H75=k8Amn1,m=2∼M−1,n=2∼M,H16=k13+k15ζAmn1,m=2∼M,n=2∼M−1,H85=k10Amn1−k4k8k3An11,0M−1×M−2,k4k12k8AnM1T,m=1∼M,n=2∼M,H26=−k15Amn2+ρs2+k14ζ2I2,m,n=2∼M−1,H86=−k11ζ0M−2×1,I2,0M−2×1T,H37=ρs2I2−A1n1Am11k1,m,n=2∼M−1,H47=k6k1A1n1Am11,m=1∼M,n=2∼M−1,H57=−Amn1,m=2∼M,n=2∼M−1,H38=k6k1A1n1Am11,m=2∼M−1,n=1∼M,H48=k17ζ2IM−k16Amn2−k62k1A1n1Am11+k72k2EM1+EM−1,m,n=1∼M,H58=k6Amn1,m=2∼M,n=1∼M,H68=k7ζ0M−2×1,I2,0M−2×1,S−S:H=00k8Amn1H14H15H160000k9ςI2H24H25H2600H31ςI20000ρs2I20k6Amn1H4200000H48k1IM00000−Amn1H580H620000−ςI2H6800k3I2H74H75−k9ςI20000H83H84H85H8600,where(45)H31=−Amn1,m=2∼M−1,n=1∼M,H42=k7ζ0M−2×1,I2,0M−2×1T,H62=k20M−2×1,I2,0M−2×1,H83=k40M−2×1,I2,0M−2×1T,H14=0M×1,k10Amn1,0M×1,m=2∼M−1,n=1∼M,H24=k11ζ0M−2×1,I2,0M−2×1,H74=k40M−2×1,I2,0M−2×1,H84=k5IM−k10k4k8EM1+EM−1,H86=−k11ζ0M−2×1,I2,0M−2×1T,H15=k12Am11A1n1−Amn2+AmM1AMM1+ρs2+k15ζ2IM,m,n=1∼M,H25=−k13+k15ζAmn1,m=2∼M−1,n=1∼M,H75=k8Amn1,m=2∼M−1,n=1∼M,H16=k15+k13ζAmn1,m=1∼M,n=2∼M−1,H85=k10Amn1+k12k4k3An11,0M×M−2,k4k12k8AnM1T,m,n=1∼M,H26=ρs2+k14ζ2I2−k15Amn2,m,n=2∼M−1,H48=k17ζ20M−1×1,I1T−k16Amn2+k72ζ2k2EM−1,m=1∼M,n=2∼M,H58=k6Amn1,m=1∼M,n=2∼M,H68=k7ζI2,0M−2×1,where(46)IM=10…001…0⋮⋮⋱⋮00…1M×M,I1=10…001…0⋮⋮⋱⋮00…1M−1×M−1,I2=10…001…0⋮⋮⋱⋮00…1M−2×M−2,I=1,0m×n=0⋯0⋮⋱⋮0⋯0m×n,Em1=10…000…0⋮⋮⋱⋮00…0m×m,Em−1=0…00⋮⋱⋮⋮0…000…01m×m.


## 4. Numerical Results

The model in Figure1 has two components: the first and third layers are composed of PZT, and the second layer is Gr/Epoxy. All the material parameters are listed in Table 1, including the elastic modulus, density, piezoelectric constant, and dielectric constant. The length, width, and thickness of the model are represented as a, b, and h, respectively, and each layer has a thickness of h/3. The values of a, b, and h are provided in Table 2 (see Tables 3 and 4).Table 1
Properties of PZT and Gr/Epoxy.
Elastic modulus (GPa)C11C12C13C22C23C33C44C55C66PZT2383.982.1923.61.9210.65.434.42.15Gr/Epoxy1562.2962.298.521.98.524.554.551.7Density (kg/m3)ρPZT7500Gr/Epoxy1600Piezoelectric constant (C/m2)e31e32e33e24e15PZT−0.13−0.14−0.28−0.01−0.01Gr/Epoxy00000Dielectric constantε11/ε0ε22/ε0ε33/ε0PZT1697.51697.51468.3Gr/Epoxy333Note. ε0=8.854×10−12F/m.Table 2
Geometric sizes of the model.
Symbols (m)abhh/34.1220.020.02/34.2110.010.01/3Table 3
Properties of PZT-2 and PZT-4.
Elastic modulus (GPa)C11C12C13C22C23C33C44C55C66PZT-2134.967.968.1134.968.1113.322.222.233.4PZT-4139.077.874.3139.074.3115.425.625.630.6Density (kg/m3)ρPZT-27500PZT-47600Piezoelectric constant (C/m2)e31e32e33e24e15PZT-2−1.82−1.829.059.789.78PZT-4−5.20−5.2015.0812.7212.72Dielectric constantε11/ε0ε22/ε0ε33/ε0PZT-2990990450PZT-4147514751300Table 4
Geometric sizes of the model.
Symbols (m)abhh/34.3110.010.01/3The semianalytical underwater laminated model based on the state space method incorporating the Laplace transform, Fourier transform, and DQM is programmed using the Mathematica platform.


### 4.1. Comparison of the Present Method and FEM without FSI

First, the piezoelectric laminate under mechanical pressure is analyzed without considering the FSI. A sinusoidal transverse load,q0, is applied to the upper surface of the model, where q0=sinπtsinπy/b. For the electrical boundary, the electric potential, ϕ, is set to zero at both the upper and lower surfaces. Meanwhile, the numerical results are compared with the simulation with the FEM.Taking the C-C boundary condition as an example, the results of both methods are shown in Figures2(a) and 2(b); the chosen sample points are x = a/2, y = b/2 in the x-y plane along the thickness, z. Figure 3 shows the distribution of the electric potential in the x-z plane at y = b/2. Figures 2 and 3 represent the results at t = 0.5 s, while Figure 4 displays the time-history curves of the vectors of the sample point at x = a/2, y = b/2, and z = h/2.Figure 2
Comparisons of the vectors obtained with the FEM and the present method atx = a/2 and y = b/2. (a) Relationship between the displacement in the (z)-direction, (w) (mm), and the thickness. (b) Relationship between the electric potential, ϕ (V), and the thickness.
(a)(b)Figure 3
Distributions of the electric potential,ϕ (V), on the x-z plane at y = b/2 obtained with the FEM and present method. (a) FEM. (b) Present method.
(a)(b)Figure 4
Comparisons of the vector-time curves obtained with the FEM and present method atx = a/2 and y = b/2. (a) Time-history curves of the displacement, (w) (mm), in the (z)-direction. (b) Time-history curves of the electric potential, ϕ (V).
(a)(b)As shown in Figure2, compared with the FEM results, the values of w (m) and ϕ (V) of the sample points obtained with the present method have similar trends along the thickness; the correlation coefficient in Figure 2(a) is 99.997%, while that in Figure 2(b) is 99.999%. The correlation coefficient can be calculated as in the following equation:(47)CorrelX=∑x−x′x−x′2.Figure3 shows that both the FEM and present method produce similar distributions, and Figure 4 shows almost the same trends over time obtained with the two methods. These similar trends, distributions, and correlation coefficients confirm the reliability of the present method.


### 4.2. Comparison of the Present Method and FEM with FSI

The dynamic responses of underwater piezoelectric laminated plates are verified by FEM. Taking the air-backed C-C boundary condition into consideration, the upper layer of the plate is exposed to the water. Therefore, a spring foundation is set on the upper layer in the FEM model, where the spring constant is zero, and the viscous damping here isCV=0.022ρWcw, which shows the underwater properties. Equation (48) denotes the relationship of the spring constant and viscous damping. Besides, the time step is set as 0.0001 s due to equation (16) of shock wave:(48)FA=−kAu−u0−dA∂u−u0∂t,where kA is spring constant and dA is viscous damping.The same sample point atx = a/2, y = b/2, and z = h/2 in the time domain is chosen. Figures 5 and 6 show the comparisons between two curves of both present method and FEM. The chosen period of FEM curves is from 0 s to 0.3 s.Figure 5
Time-history curves of different vectors for the C-C air-backed boundary condition. (a) Time history ofw (m). (b) Time history of ϕ (V).
(a)(b)Figure 6
Time-history curves of different vectors for the C-C water-backed boundary condition. (a) Time history ofw (m). (b) Time history of ϕ (V).
(a)(b)It is clear that the time history ofw and the peak values of ϕ are close, but the curves of ϕ show a sharp fluctuation during the first 0.03 s in both the present method and FEM. The FEM takes more steps at first 0.03 seconds to analyze the shock wave; besides the values are more unstable. While the present method owns a stable change during the first 0.003 s, the values show a slight fluctuation in the time domain due to the limitation of computational accuracy.


### 4.3. Influence of the FSI Based on the Present Method

The FSI is coupled to the piezoelectric laminate in this section. The vectors on the upper and lower surfaces are applied as equations (19) and (20). In addition, the plane stress on the surfaces is defined as zero to satisfy the one-dimensional FSI theory. To investigate the influence of the FSI, air-backed and water-backed laminated plates are simulated for each boundary condition. The electric vectors are defined as ϕ = 0 on the upper surface and the lower surface.On account of the equations for the shock wave given in equations (15) and (16), for pm and td, W = 1 kg and S = 500 m to establish the boundary conditions on the upper and lower surfaces. The results for the four boundary conditions of S-S, C-C, C-F, and C-S are shown in Figures 7–14. In Figures 7–10, t = 0.05 s, and sample points of x = a/2 and y = b/2 along the thickness, z, are chosen. The values on the boundaries of some figures are shown in white, which indicates that they are much larger or much smaller than the rest of the figure.Figure 7
Distribution ofϕ (V) at the x-z plane, y = b/2, for the C-C boundary condition. (a) Air-backed. (b) Water-backed.
(a)(b)Figure 8
Distribution ofϕ (V) at the x-z plane, y = b/2, for the C-F boundary condition. (a) Air-backed. (b) Water-backed.
(a)(b)Figure 9
Distribution ofϕ (V) at the x-z plane, y = b/2, for the C-S boundary condition. (a) Air-backed. (b) Water-backed.
(a)(b)Figure 10
Distribution ofϕ (V) at the x-z plane, y = b/2, for the S-S boundary condition. (a) Air-backed. (b) Water-backed.
(a)(b)Figure 11
Time-history curves of different vectors for the C-C boundary condition. (a) Time history ofw (m). (b) Time history of ϕ (V).
(a)(b)Figure 12
Time-history curves of different vectors for the C-F boundary condition. (a) Time history ofw (m). (b) Time history of ϕ (V).
(a)(b)Figure 13
Time-history curves of different vectors for the C-S boundary condition. (a) Time history ofw (m). (b) Time history of ϕ (V).
(a)(b)Figure 14
Time-history curves of different vectors for the S-S boundary condition. (a) Time history ofw (m). (b) Time history of ϕ (V).
(a)(b)For various boundary conditions, both the distributions and values of the electric potential at thex-z plane are different. For the C-C, C-F, and C-S boundary conditions, the distributions of the values are almost the same under both air-backed and water-backed conditions, and the FSI on the surfaces lessens the external force, which leads to a reduction in ϕ. For the S-S condition, the distribution of value shows an obvious difference; the maximum values shown in the figures are at different positions, which is because there are rotational degrees of freedom and the stress tends to disperse at 1/4 and 3/4 of the x-axis under the smaller external force. The same phenomenon can also be observed in the C-S boundary condition, where the large values disperse near the simply supported boundary. For the C-F condition, the stress concentrates on the clamped boundary and the values are obviously large, while the values close to the free boundary are almost 0.The response time history of each vector can also be calculated. Here, the sample point atx = a/2, y = b/2, and z = h/2 in the time domain is chosen for analysis.In Figures11–14, it is clear that the values of each vector at the chosen point are lessened by the rarefaction wave on the lower surface, and the values under water-backed condition are about half of the values under air-backed condition. The decreasing of water-backed value seems to be slower than that of air-backed value. The water-backed condition shows more buffering against the shock wave. The four kinds of boundary conditions indicate the close value of w, and ϕ under S-S condition seems to be the smallest due to the more degrees of freedom on the boundary.


## 4.1. Comparison of the Present Method and FEM without FSI

First, the piezoelectric laminate under mechanical pressure is analyzed without considering the FSI. A sinusoidal transverse load,q0, is applied to the upper surface of the model, where q0=sinπtsinπy/b. For the electrical boundary, the electric potential, ϕ, is set to zero at both the upper and lower surfaces. Meanwhile, the numerical results are compared with the simulation with the FEM.Taking the C-C boundary condition as an example, the results of both methods are shown in Figures2(a) and 2(b); the chosen sample points are x = a/2, y = b/2 in the x-y plane along the thickness, z. Figure 3 shows the distribution of the electric potential in the x-z plane at y = b/2. Figures 2 and 3 represent the results at t = 0.5 s, while Figure 4 displays the time-history curves of the vectors of the sample point at x = a/2, y = b/2, and z = h/2.Figure 2
Comparisons of the vectors obtained with the FEM and the present method atx = a/2 and y = b/2. (a) Relationship between the displacement in the (z)-direction, (w) (mm), and the thickness. (b) Relationship between the electric potential, ϕ (V), and the thickness.
(a)(b)Figure 3
Distributions of the electric potential,ϕ (V), on the x-z plane at y = b/2 obtained with the FEM and present method. (a) FEM. (b) Present method.
(a)(b)Figure 4
Comparisons of the vector-time curves obtained with the FEM and present method atx = a/2 and y = b/2. (a) Time-history curves of the displacement, (w) (mm), in the (z)-direction. (b) Time-history curves of the electric potential, ϕ (V).
(a)(b)As shown in Figure2, compared with the FEM results, the values of w (m) and ϕ (V) of the sample points obtained with the present method have similar trends along the thickness; the correlation coefficient in Figure 2(a) is 99.997%, while that in Figure 2(b) is 99.999%. The correlation coefficient can be calculated as in the following equation:(47)CorrelX=∑x−x′x−x′2.Figure3 shows that both the FEM and present method produce similar distributions, and Figure 4 shows almost the same trends over time obtained with the two methods. These similar trends, distributions, and correlation coefficients confirm the reliability of the present method.


## 4.2. Comparison of the Present Method and FEM with FSI

The dynamic responses of underwater piezoelectric laminated plates are verified by FEM. Taking the air-backed C-C boundary condition into consideration, the upper layer of the plate is exposed to the water. Therefore, a spring foundation is set on the upper layer in the FEM model, where the spring constant is zero, and the viscous damping here isCV=0.022ρWcw, which shows the underwater properties. Equation (48) denotes the relationship of the spring constant and viscous damping. Besides, the time step is set as 0.0001 s due to equation (16) of shock wave:(48)FA=−kAu−u0−dA∂u−u0∂t,where kA is spring constant and dA is viscous damping.The same sample point atx = a/2, y = b/2, and z = h/2 in the time domain is chosen. Figures 5 and 6 show the comparisons between two curves of both present method and FEM. The chosen period of FEM curves is from 0 s to 0.3 s.Figure 5
Time-history curves of different vectors for the C-C air-backed boundary condition. (a) Time history ofw (m). (b) Time history of ϕ (V).
(a)(b)Figure 6
Time-history curves of different vectors for the C-C water-backed boundary condition. (a) Time history ofw (m). (b) Time history of ϕ (V).
(a)(b)It is clear that the time history ofw and the peak values of ϕ are close, but the curves of ϕ show a sharp fluctuation during the first 0.03 s in both the present method and FEM. The FEM takes more steps at first 0.03 seconds to analyze the shock wave; besides the values are more unstable. While the present method owns a stable change during the first 0.003 s, the values show a slight fluctuation in the time domain due to the limitation of computational accuracy.


## 4.3. Influence of the FSI Based on the Present Method

The FSI is coupled to the piezoelectric laminate in this section. The vectors on the upper and lower surfaces are applied as equations (19) and (20). In addition, the plane stress on the surfaces is defined as zero to satisfy the one-dimensional FSI theory. To investigate the influence of the FSI, air-backed and water-backed laminated plates are simulated for each boundary condition. The electric vectors are defined as ϕ = 0 on the upper surface and the lower surface.On account of the equations for the shock wave given in equations (15) and (16), for pm and td, W = 1 kg and S = 500 m to establish the boundary conditions on the upper and lower surfaces. The results for the four boundary conditions of S-S, C-C, C-F, and C-S are shown in Figures 7–14. In Figures 7–10, t = 0.05 s, and sample points of x = a/2 and y = b/2 along the thickness, z, are chosen. The values on the boundaries of some figures are shown in white, which indicates that they are much larger or much smaller than the rest of the figure.Figure 7
Distribution ofϕ (V) at the x-z plane, y = b/2, for the C-C boundary condition. (a) Air-backed. (b) Water-backed.
(a)(b)Figure 8
Distribution ofϕ (V) at the x-z plane, y = b/2, for the C-F boundary condition. (a) Air-backed. (b) Water-backed.
(a)(b)Figure 9
Distribution ofϕ (V) at the x-z plane, y = b/2, for the C-S boundary condition. (a) Air-backed. (b) Water-backed.
(a)(b)Figure 10
Distribution ofϕ (V) at the x-z plane, y = b/2, for the S-S boundary condition. (a) Air-backed. (b) Water-backed.
(a)(b)Figure 11
Time-history curves of different vectors for the C-C boundary condition. (a) Time history ofw (m). (b) Time history of ϕ (V).
(a)(b)Figure 12
Time-history curves of different vectors for the C-F boundary condition. (a) Time history ofw (m). (b) Time history of ϕ (V).
(a)(b)Figure 13
Time-history curves of different vectors for the C-S boundary condition. (a) Time history ofw (m). (b) Time history of ϕ (V).
(a)(b)Figure 14
Time-history curves of different vectors for the S-S boundary condition. (a) Time history ofw (m). (b) Time history of ϕ (V).
(a)(b)For various boundary conditions, both the distributions and values of the electric potential at thex-z plane are different. For the C-C, C-F, and C-S boundary conditions, the distributions of the values are almost the same under both air-backed and water-backed conditions, and the FSI on the surfaces lessens the external force, which leads to a reduction in ϕ. For the S-S condition, the distribution of value shows an obvious difference; the maximum values shown in the figures are at different positions, which is because there are rotational degrees of freedom and the stress tends to disperse at 1/4 and 3/4 of the x-axis under the smaller external force. The same phenomenon can also be observed in the C-S boundary condition, where the large values disperse near the simply supported boundary. For the C-F condition, the stress concentrates on the clamped boundary and the values are obviously large, while the values close to the free boundary are almost 0.The response time history of each vector can also be calculated. Here, the sample point atx = a/2, y = b/2, and z = h/2 in the time domain is chosen for analysis.In Figures11–14, it is clear that the values of each vector at the chosen point are lessened by the rarefaction wave on the lower surface, and the values under water-backed condition are about half of the values under air-backed condition. The decreasing of water-backed value seems to be slower than that of air-backed value. The water-backed condition shows more buffering against the shock wave. The four kinds of boundary conditions indicate the close value of w, and ϕ under S-S condition seems to be the smallest due to the more degrees of freedom on the boundary.


## 5. Discussion and Conclusions

A novel and efficient method is proposed in this paper to solve the problem of the dynamic responses of underwater piezoelectric laminated plates. This method first combines the Laplace transform, Fourier transform, and DQM to change the differential forms ofx, y, and t. The state space method is then introduced to deduce the fundamental equations of the piezoelectric material. Finally, the FSI is employed to couple the piezoelectric laminate and influence of the water.In this study, the finite element method is chosen to verify the results of the present method for the piezoelectric problem. Although there is some deviation in the values of some vectors, the correlation coefficients of the results and overall trends indicate that the present method is accurate and reliable. For the comparison of FSI on the C-C condition, the FEM dealing with shock load reveals instability, and the present method has its advantage in the dynamic response.As the FEM has limits in coupling the piezoelectric-mechanical-fluid problem, the one-dimensional FSI is used to analyze this condition. The results show that the present method has good convergence in conveying external forces such as a shock wave in the form of the FSI. Additionally, the FSI has an obvious buffer effect on the external force.When the underwater shock was applied on the model under different boundary conditions, it can be easily found that the water on the surfaces has buffering to the flow. Different boundary conditions show their specific distributions of the vectors, which can be the basis of sensors under different supporting boundaries.The method proposed in this paper aims to analyze underwater piezoelectric laminates, which can be applied as sensors and actuators in marine engineering. Furthermore, it also provides a new semianalytical solution to multiphysics problems.



---

*Source: 2439521-2020-08-30.xml*