# A New Numerical Procedure for Vibration Analysis of Beam under Impulse and Multiharmonics Piezoelectric Actuators

Authors: Yassin Belkourchia; Lahcen Azrar
Journal: Journal of Applied Mathematics
 (2020)

---

## Abstract

The dynamic behavior of structures with piezoelectric patches is governed by partial differential equations with strong singularities. To directly deal with these equations, well adapted numerical procedures are required. In this work, the differential quadrature method (DQM) combined with a regularization procedure for space and implicit scheme for time discretization is used. The DQM is a simple method that can be implemented with few grid points and can give results with a good accuracy. However, the DQM presents some difficulties when applied to partial differential equations involving strong singularities. This is due to the fact that the subsidiaries of the singular functions cannot be straightforwardly discretized by the DQM. A methodological approach based on the regularization procedure is used here to overcome this difficulty and the derivatives of the Dirac-delta function are replaced by regularized smooth functions. Thanks to this regularization, the resulting differential equations can be directly discretized using the DQM. The efficiency and applicability of the proposed approach are demonstrated in the computation of the dynamic behavior of beams for various boundary conditions and excited by impulse and Multiharmonics piezoelectric actuators. The obtained numerical results are well compared to the developed analytical solution.


---

## Body

## 1. Introduction

Many industrial and engineering problems can be generally modeled by partial differential equations. Currently, various analytical and numerical methods have been developed in the last decades to deal with these equations. Often, analytical methods are preferred because they give an exact solution allowing them to get useful information on the domain of the problem. However, analytical methods are almost available for simple engineering problems with simple geometries. To address these weaknesses, many researchers have resorted to numerical methods. The finite difference, finite volume, and finite element methods are the widely used numerical methods. As these methods are classically part of the low order methods, it is necessary to refine the mesh and this may require a relatively high computational effort.To overcome the abovementioned difficulties of low order methods, some researchers have thought of using high order methods. The Differential Quadrature Method (DQM) is one of the most widely used methods to solve this weakness due to many features such as high accuracy and performance. The DQM, initially presented by Bellman et al. [1, 2] in the mid-1970s, is an amazing method for the direct numerical solution of partial differential equations showing up in many fields in engineering, mathematics, and physics [3–5]. Most applications of this method involve static and dynamic analysis of structural components such as beams and plates [6–8]. In addition, Bert and Malik [9] presented an analysis of thin rectangular plates simply supported on two opposite ends, based on the classical thin plate theory. Du et al. [10] introduced a generalized DQM to examine the buckling problems of rectangular plates with internal support and variable bending stiffness. Wang et al. [11] investigated the vibration and buckling problems of thin rectangular plates with nonlinear distributed loads along two adjacent plate edges with nine boundary conditions by using the differential quadrature method. Bashan [12] implemented for the numerical solution of nonlinear Kawahara equation via the Crank-Nicolson-Differential Quadrature Method based on modified cubic B-splines and Bashan et al. [13] applied DQM based on modified cubic B-splines for numerical solution of the complex modified Korteweg-de Vries equation. Also, Bashan [14, 15] investigated numerical solutions of the system of coupled Korteweg-de Vries equation based on Finite difference method and DQM and numerical solutions of the coupled Korteweg-de Vries equation based on a combination of Crank-Nicolson method and quintic B-spline based differential quadrature method. Ucar et al. [16] proposed a numerical approximate solution of the nonlinear modified Burgers equation via the modified cubic B-spline and DQM. Also, Tang and Wang [17] used the DQM to solve the buckling analysis of symmetrically laminated rectangular plates under plane compression along two opposite edges. It has been observed from the evidence provided by various researchers that DQM is computationally productive and can give excellent results for the problems under investigation.Despite the abovementioned advantages of the DQM, it presents certain challenges when applied to partial differential equations containing singular functions such as the derivative of the Dirac-delta function. To overcome this difficulty, some authors have suggested coupling the DQM and the integral quadrature method (IQM) in which this type of problem can be handled [18]. This approach is applied to the manipulation of the Dirac-delta function applied to the problem of vibration of beams and rectangular plates subjected to a moving point load. On the other hand, Eftekhari [19] presented a methodological procedure for the numerical solution of the moving load problem, and Eftekhari and Jafari [20] combined the finite element method (FEM) and the DQM to study the free and forced vibration and buckling of rectangular plates. As previously mentioned in [20], this approach may have certain difficulties when applied to problems contains a singular time-dependent function. In general, it is more favorable to find the solution of these problems by the DQM itself and not by combining the DQM with other techniques.In addition, some researchers have also mentioned the difficulty described above when similar methods such as collocation and finite difference methods are employed. In [21], it is argued that the difficulty of the collocation method in dealing with such singular functions is due to the Gibbs phenomenon in which numerical solutions oscillate around. To solve this problem, a regularization of the singular function has been proposed in order to obtain a softer representation of the singular function and to stabilize the undesirable oscillatory behavior of solutions close to the singularities.Many regularization approaches have been developed by different researchers in this subject. Wei et al. [22] presented a computational approach of the regularized Dirac-delta function based on the discrete singular convolution algorithm for vibration analysis of rectangular plates with partial internal line or point supports. Burko and Khanna [23] suggested a regularization of the Dirac-delta function using the Gaussian function. Walden [24] solved a few differential equations with singular source terms by finite difference and finite element methods and showed that the full order of convergence can be achieved if the singularities are treated in the right way. Engquist et al. [25] presented two methods to construct consistent approximations for the regularization of the Dirac-delta function. Rivera et al. [26] proposed numerical resolution of the hyperbolic heat equation using smoothed mathematical functions for approximation of the Heaviside and Dirac-delta functions.In the present work, a numerical procedure based on the combination of the DQM with the regularization of the derivatives of the Dirac-delta function is elaborated for the numerical solution of the vibration response of the beam under the impulse and multiharmonic piezoelectric actuators. Based on this regularization, the DQM can be applied directly to discretize the resulting partial differential equations. To establish its applicability and reliability, the proposed approach is applied here to solve the static and dynamic analyses of beams excited by piezoelectric pulses and multiharmonic actuators, where the installation of piezoelectric actuators is characterized by derivatives of a Dirac-delta function. Various numerical results are presented and compared with the analytical solution developed herein. The presented numerical results demonstrated that the proposed methodology is simple, efficient, and accurate.


## 2. Mathematical Formulation

### 2.1. Equilibrium Equation

In this work, a beam with the lengthL, cross-section area A, Young’s modulus E, moment of area I, and density ρ excited by Na piezoelectric actuators is considered as shown in Figure 1. Based on Euler-Bernoulli beam model, the transverse motion is modeled by the following partial differential equation [27]
(1)EI∂4wx,t∂x4+ρA∂2wx,t∂t2=fx,t,where
(2)fx,t=b∑i=0Na∂2Mia∂x2=b∑i=1NaEpid31ihpi+h2Vitδ′x−xa1i−δ′x−xa2i,in which wx,t is the transverse displacement of the beam, x is the axial coordinate, t is the time, b and h are the width and thickness of the beam, respectively. Na is the number of the piezoelectric actuators. Mia is the force moment induced by the piezoelectric actuator i, and d31i, Epi, and hpi are the piezoelectric strain constant, Young’s modulus, thickness of the i-th actuator, respectively. δ′. is the derivative of the Dirac-delta function, xa1i, xa2i are the coordinates of the two ends of i-th actuator and Vit is the voltage applied on the i-th actuator.Figure 1
Beam under a PZT actuators.


### 2.2. Regularization Procedure

As mentioned above, with regard to the derivative of the Dirac-delta function, the direct implementation of this function by the DQM is not a simple matter due to the particular characterizations associated with it. One way to solve this challenge is to replace the derivative of the Dirac-delta function with a regular and soft function. In this context, different forms of regularized derivation of the Dirac-delta function have been proposed in the literature [22]. In this work, the Hermite polynomials are used [28]. An approximation of the distribution delta is constructed using the usual Hermite polynomials H2n as follows
(3)δMx−x0=1αexp−x−x022α2∑n=0M2−14n12πn!H2nx−x02α.For a very smallα, the function δMx−x0 becomes identical to the functional delta when the degree of polynomial M is fixed.
(4)limα→0δMx−x0=δx−x0In addition, the derivatives of the regularized formulation give an approximation to the derivatives of the Dirac-delta function and are given by(5)δMlx−x0=2−l/2αl+1exp−x−x022α2∑n=0M2−14n×−1l12πn!H2n+lx−x02α.Equations (3) and (5) mean that the differentiations have been transformed into an algebraic process in the approximate representation. This important feature of this approximate representation allows it to be a powerful computational tool for solving various partial and ordinary differential equations with singularities. The parameter α must be as small as possible with fixed M.In view of Eq. (5), the excitation force can be rewritten as
(6)fx,t=b∑i=0NaEpid31ihpi+h2Vit1α22exp−x−xa2i22α2∑m=0M/2−14m×12πm!H2mx−xa2i2α−1α22exp−x−xa1i22α2∑m=0M2−14m×12πm!H2mx−xa1i2α.It should be observed that decreasing of the regularization parameterα yields a more accurate representation. As the DQM is a higher-order method, this can considerably increase the cost of calculation, especially when the value α is too small. Consequently, it is necessary to find a suitable value of α for the numerical accuracy and efficiency of the calculation.


## 2.1. Equilibrium Equation

In this work, a beam with the lengthL, cross-section area A, Young’s modulus E, moment of area I, and density ρ excited by Na piezoelectric actuators is considered as shown in Figure 1. Based on Euler-Bernoulli beam model, the transverse motion is modeled by the following partial differential equation [27]
(1)EI∂4wx,t∂x4+ρA∂2wx,t∂t2=fx,t,where
(2)fx,t=b∑i=0Na∂2Mia∂x2=b∑i=1NaEpid31ihpi+h2Vitδ′x−xa1i−δ′x−xa2i,in which wx,t is the transverse displacement of the beam, x is the axial coordinate, t is the time, b and h are the width and thickness of the beam, respectively. Na is the number of the piezoelectric actuators. Mia is the force moment induced by the piezoelectric actuator i, and d31i, Epi, and hpi are the piezoelectric strain constant, Young’s modulus, thickness of the i-th actuator, respectively. δ′. is the derivative of the Dirac-delta function, xa1i, xa2i are the coordinates of the two ends of i-th actuator and Vit is the voltage applied on the i-th actuator.Figure 1
Beam under a PZT actuators.


## 2.2. Regularization Procedure

As mentioned above, with regard to the derivative of the Dirac-delta function, the direct implementation of this function by the DQM is not a simple matter due to the particular characterizations associated with it. One way to solve this challenge is to replace the derivative of the Dirac-delta function with a regular and soft function. In this context, different forms of regularized derivation of the Dirac-delta function have been proposed in the literature [22]. In this work, the Hermite polynomials are used [28]. An approximation of the distribution delta is constructed using the usual Hermite polynomials H2n as follows
(3)δMx−x0=1αexp−x−x022α2∑n=0M2−14n12πn!H2nx−x02α.For a very smallα, the function δMx−x0 becomes identical to the functional delta when the degree of polynomial M is fixed.
(4)limα→0δMx−x0=δx−x0In addition, the derivatives of the regularized formulation give an approximation to the derivatives of the Dirac-delta function and are given by(5)δMlx−x0=2−l/2αl+1exp−x−x022α2∑n=0M2−14n×−1l12πn!H2n+lx−x02α.Equations (3) and (5) mean that the differentiations have been transformed into an algebraic process in the approximate representation. This important feature of this approximate representation allows it to be a powerful computational tool for solving various partial and ordinary differential equations with singularities. The parameter α must be as small as possible with fixed M.In view of Eq. (5), the excitation force can be rewritten as
(6)fx,t=b∑i=0NaEpid31ihpi+h2Vit1α22exp−x−xa2i22α2∑m=0M/2−14m×12πm!H2mx−xa2i2α−1α22exp−x−xa1i22α2∑m=0M2−14m×12πm!H2mx−xa1i2α.It should be observed that decreasing of the regularization parameterα yields a more accurate representation. As the DQM is a higher-order method, this can considerably increase the cost of calculation, especially when the value α is too small. Consequently, it is necessary to find a suitable value of α for the numerical accuracy and efficiency of the calculation.


## 3. Solution Procedure

A numerical methodological approach based on the DQM combined with a regularization procedure is adapted to space and an implicit scheme for the time derivative.


### 3.1. Differential Quadrature Method

The derivative with respect to the spatial variable, it is discretized by applying the DQM. The principle of this method consists of approximating the derivative of a function at any location by a weighted linear sum of the values of the function at all points of the discretization of the domain. Suppose that the functionwx is sufficiently smooth over the interval x1,xN. In this interval, N distinct nodes are defined.
(7)x1<x2<⋯<xN.The function values on these points are assumed to be(8)w1,w2⋯,wN.According to the DQM, the first and second-order derivatives on each of these nodes are given by [29].
(9a)dwxidx≈∑j=1Naijwji=1,2,⋯,N,(9b)d2wxidx2≈∑j=1Nbijwji=1,2,⋯,N.The coefficientsaij and bij are the weighting coefficients of the first and second-order derivatives with respect to x, respectively. The coefficients aij and bij are given as follows [29].
(10a)aij=1xj−xi∏k=1,k≠i,jNxi−xkxj−xki,j=1,2,⋯,Ni≠j,(10b)aii=−∑j=1,j≠iNaiji=1,2,⋯,N,(10c)bij=2aijaii−aijxj−xii,j=1,2,⋯,Ni≠j,(10d)bij=−∑j=1,j≠iNbiji=1,2,⋯,N.Similarly, we can obtain higher-order derivative formulas by using the higher weighted coefficients, which are expressed ineijm to avoid confusion. They are characterized by induction [4].
(11a)eijm=maijeiim−1−eiim−1xi−xji,j=1,2,⋯,Ni≠j,m=2,3,⋯,N−1,(11b)eiim=−∑j=1,j≠iNeijmi=1,2,⋯,N.One of the key factors in the accuracy of DQ solutions is the choice of grid points. The zeros of some orthogonal polynomials are commonly adopted as grid points. In this work, the DQM grid points are taken nonuniformly spaced and are given by the following equations [6]
(12)xi=x1+121−cosi−1N−1πxN−x1.


### 3.2. Numerical Accuracy of the DQM

Consider a functionwx, which is approximated by the Lagrange interpolation polynomial of degree N−1. The error for the r−th order derivative approximation of this function at point xi can be obtained as [29].
(13)Erwxi=wNξMrxiN!,i=1,2,⋯,N,where
(14)Mx=∏j=1Nx−xj.As a result, form equation (13), the accuracy of the DQM can be directly proportional to N. This means that a result can be achieved with very high accuracy even if the number of grid points N is small.


### 3.3. DQM Analogues

For numerical solution,n grid points with coordinates x1,x2,⋯,xn in the x-direction are considered. Applying the quadrature rule, to Eq. (1), the following ordinary differential system is obtained
(15)EI∑k=1neik4wkt+ρA∂2wit∂t2=fit,i=1,⋯,n.This system can be rewritten in the following matrix from(16)KWt+MW¨t=Ft,where M and K are the resulting mass and stiffness matrices of the beam, Wt and W¨t are the displacement and acceleration vectors, and Ft is the charge vector applied by the piezoelectric actuators.These terms are given by(17)K=EIe4,M=ρAI,Wt=wx1,twx2,t⋯wxn,tT,W¨t=w¨x1,tw¨x2,t⋯w¨xn,tT,(18)Ft=b2α2∑i=0NaEpid31ihpi+h2Vitexp−x1−xa2i22α2∑m=0M2−14m×12πm!H2mx1−xa2i2α−exp−x1−xa1i22α2∑m=0M2−14m×12πm!H2mx1−xa1i2α⋯exp−xn−xa2i22α2∑m=0M2−14m×12πm!H2mxn−xa2i2α−exp−xn−xa1i22α2∑m=0M2−14m×12πm!H2mxn−xa1i2αT,where I is the identity matrix of order n×n and e4 is the DQM weighting coefficient matrix of the fourth-order derivative.The discrete classical boundary conditions of the beam atx=0 and x=1, using the DQ method, can be written as following:
(19)w1t=0,∑k=1ne1kn0wkt=0,(20)wnt=0,∑k=1nenkn1wkt=0,where n0 and nl can be chosen from 1, 2, or 3. The selection of the values n0 and nl can have the following traditional boundary conditions [6]:
(i)
simply supported:n0=2; n1=2(ii)
clamped-clamped:n0=1; n1=1(iii)
clamped-simply supported:n0=1; n1=2(iv)
clamped-free:n0=1; n1=3(v)
free-free:n0=2; n1=3Similarly, the boundary conditions (19) can also be expressed in a matrix form
(21)KBBWBt+KBIWIt=0,where the subscripts B and I indicate the grid points used for writing the quadrature analog of the boundary conditions and the governing differential equation, respectively. It is noted that the dimensions of matrices KBB and KBI are 4×4 and 4×n−4, respectively.Implementing the boundary conditions into Eq. (16) leads to the following system of ordinary differential equations
(22)K~WIt+M~W¨It=FIt,whereinK~=KII−KIBKBB−1KBI and M~=MII; where K~, M~, and KIB are matrices of order n−4×n−4, n−4×n−4, and n−4×4, respectively.Equation (22) can be solved by using various step-by-step time integration schemes. In this study, the time derivative is discretized using the centered finite difference scheme, then
(23)∂2WI∂t2=WIi+1−2WIi+WIi−1Δt2.Substituting this approximate expression acceleration into (22), one gets:
(24)K~WIi+1+M~WIi+1−2WIi+WIi−1Δt2=FIi+1+initialconditions.


## 3.1. Differential Quadrature Method

The derivative with respect to the spatial variable, it is discretized by applying the DQM. The principle of this method consists of approximating the derivative of a function at any location by a weighted linear sum of the values of the function at all points of the discretization of the domain. Suppose that the functionwx is sufficiently smooth over the interval x1,xN. In this interval, N distinct nodes are defined.
(7)x1<x2<⋯<xN.The function values on these points are assumed to be(8)w1,w2⋯,wN.According to the DQM, the first and second-order derivatives on each of these nodes are given by [29].
(9a)dwxidx≈∑j=1Naijwji=1,2,⋯,N,(9b)d2wxidx2≈∑j=1Nbijwji=1,2,⋯,N.The coefficientsaij and bij are the weighting coefficients of the first and second-order derivatives with respect to x, respectively. The coefficients aij and bij are given as follows [29].
(10a)aij=1xj−xi∏k=1,k≠i,jNxi−xkxj−xki,j=1,2,⋯,Ni≠j,(10b)aii=−∑j=1,j≠iNaiji=1,2,⋯,N,(10c)bij=2aijaii−aijxj−xii,j=1,2,⋯,Ni≠j,(10d)bij=−∑j=1,j≠iNbiji=1,2,⋯,N.Similarly, we can obtain higher-order derivative formulas by using the higher weighted coefficients, which are expressed ineijm to avoid confusion. They are characterized by induction [4].
(11a)eijm=maijeiim−1−eiim−1xi−xji,j=1,2,⋯,Ni≠j,m=2,3,⋯,N−1,(11b)eiim=−∑j=1,j≠iNeijmi=1,2,⋯,N.One of the key factors in the accuracy of DQ solutions is the choice of grid points. The zeros of some orthogonal polynomials are commonly adopted as grid points. In this work, the DQM grid points are taken nonuniformly spaced and are given by the following equations [6]
(12)xi=x1+121−cosi−1N−1πxN−x1.


## 3.2. Numerical Accuracy of the DQM

Consider a functionwx, which is approximated by the Lagrange interpolation polynomial of degree N−1. The error for the r−th order derivative approximation of this function at point xi can be obtained as [29].
(13)Erwxi=wNξMrxiN!,i=1,2,⋯,N,where
(14)Mx=∏j=1Nx−xj.As a result, form equation (13), the accuracy of the DQM can be directly proportional to N. This means that a result can be achieved with very high accuracy even if the number of grid points N is small.


## 3.3. DQM Analogues

For numerical solution,n grid points with coordinates x1,x2,⋯,xn in the x-direction are considered. Applying the quadrature rule, to Eq. (1), the following ordinary differential system is obtained
(15)EI∑k=1neik4wkt+ρA∂2wit∂t2=fit,i=1,⋯,n.This system can be rewritten in the following matrix from(16)KWt+MW¨t=Ft,where M and K are the resulting mass and stiffness matrices of the beam, Wt and W¨t are the displacement and acceleration vectors, and Ft is the charge vector applied by the piezoelectric actuators.These terms are given by(17)K=EIe4,M=ρAI,Wt=wx1,twx2,t⋯wxn,tT,W¨t=w¨x1,tw¨x2,t⋯w¨xn,tT,(18)Ft=b2α2∑i=0NaEpid31ihpi+h2Vitexp−x1−xa2i22α2∑m=0M2−14m×12πm!H2mx1−xa2i2α−exp−x1−xa1i22α2∑m=0M2−14m×12πm!H2mx1−xa1i2α⋯exp−xn−xa2i22α2∑m=0M2−14m×12πm!H2mxn−xa2i2α−exp−xn−xa1i22α2∑m=0M2−14m×12πm!H2mxn−xa1i2αT,where I is the identity matrix of order n×n and e4 is the DQM weighting coefficient matrix of the fourth-order derivative.The discrete classical boundary conditions of the beam atx=0 and x=1, using the DQ method, can be written as following:
(19)w1t=0,∑k=1ne1kn0wkt=0,(20)wnt=0,∑k=1nenkn1wkt=0,where n0 and nl can be chosen from 1, 2, or 3. The selection of the values n0 and nl can have the following traditional boundary conditions [6]:
(i)
simply supported:n0=2; n1=2(ii)
clamped-clamped:n0=1; n1=1(iii)
clamped-simply supported:n0=1; n1=2(iv)
clamped-free:n0=1; n1=3(v)
free-free:n0=2; n1=3Similarly, the boundary conditions (19) can also be expressed in a matrix form
(21)KBBWBt+KBIWIt=0,where the subscripts B and I indicate the grid points used for writing the quadrature analog of the boundary conditions and the governing differential equation, respectively. It is noted that the dimensions of matrices KBB and KBI are 4×4 and 4×n−4, respectively.Implementing the boundary conditions into Eq. (16) leads to the following system of ordinary differential equations
(22)K~WIt+M~W¨It=FIt,whereinK~=KII−KIBKBB−1KBI and M~=MII; where K~, M~, and KIB are matrices of order n−4×n−4, n−4×n−4, and n−4×4, respectively.Equation (22) can be solved by using various step-by-step time integration schemes. In this study, the time derivative is discretized using the centered finite difference scheme, then
(23)∂2WI∂t2=WIi+1−2WIi+WIi−1Δt2.Substituting this approximate expression acceleration into (22), one gets:
(24)K~WIi+1+M~WIi+1−2WIi+WIi−1Δt2=FIi+1+initialconditions.


## 4. Numerical Results and Discussion

In order to demonstrate the applicability of the proposed methodological approach and its numerical implementation, multiple computational examples are investigated. Firstly, the static analysis of submitted to a piezoelectric actuator with constant voltageV0 is considered to demonstrate the feasibility and effectiveness of the proposed procedure. This numerical example is introduced to further verify the exactitude and convergence of the proposed approach. The second example concerns the study of a beam exited by multiharmonic piezoelectric actuators. The third example is focused on the study of a beam excited by various types of impulse piezoelectric actuators. Three boundary conditions, namely, simply supported edges (S-S), clamped edges (C-C), and clamped free edges (C-F) conditions are considered in the second and third examples. In the present computation, the parameter M is taken as M=4 for all examples. The analytical solutions of simply supported beam exited by harmonic and impulse excitation using a piezoelectric actuator and for impulse excitation by a piezoelectric actuator is developed and given in the Appendices A and B. The accuracy of the proposed method can be verified by comparing the calculated results with those of the developed analytically ones.


### 4.1. Beams under a Piezoelectric Actuator with Constant Voltage

In this subsection, we consider a beam with lengthL, Young’s modulus E, moment of area I under one PZT actuator of thickness hp, the piezoelectric strain constant d31, and Young’s modulus Ep with a constant voltage V0. The PZT actuator is assumed to be perfectly bonded and its stiffness is neglected due to its limited contribution to the behavior of the beam, as shown in Figure 1.This problem can be modeled by the differential equation reduced from Eq. (1) as follows
(25)EI∂4wx∂x4=fx,where fx is given by Eq. (6) in which Vt=V0. By using the DOM procedure, Eq. (25) is reduced to the following algebraic equation
(26)KW=F,where K, W, and F are given by
(27)K=EIe4,W=wx1wx2⋯wxnT,F=bEpd31hp+h2V0χ,The vectorχ is given by:
(28)χ=12α2exp−x1−xa222α2∑m=0M2−14m×12πm!H2mx1−xa22α−exp−x1−xa122α2∑m=0M2−14m×12πm!H2mx1−xa12α⋯exp−xn−xa222α2∑m=0M2−14m×12πm!H2mxn−xa22α−exp−xn−xa122α2∑m=0M2−14m×12πm!H2mxn−xa12αT.To ensure the validity of the proposed methodology and its implementation, the problem of a beam simply supported excited by a concentrated actuator is addressed by applying the proposed methodology for various actuator locations and two different values of the regularization parameterα=0.25 and α=0.08, respectively. The error variations used in numerical solutions is defined by
(29)wDQM−wExactwExact×100.In accordance with the number of grid points (n), the error variations for different locations xa1,xa2 of the actuator are shown in Figures 2 and 3.Figure 2
Convergence and accuracy of the dimensionless deflection of a simply supported beam under a constant voltageV0 applied by actuator at various locations for (α=0.25) at x/L=1/2 (a), x/L=1/5 (b), x/L=2/3 (c), and x/L=3/4 (d).
(a)(b)(c)(d)Figure 3
Convergence and accuracy of the dimensionless deflection of a simply supported beam under a constant voltageV0 applied by actuator at different locations for (α=0.08) at x/L=1/2 (a), x/L=1/5 (b), x/L=2/3 (c), and x/L=3/4 (d).
(a)(b)(c)(d)Figure2 shows the variance of the percent error in the numerical results for α=0.25 and for different actuator locations. It is clearly noticed that the results obtained converge uniformly towards their final values. On the other hand, the numerical precision of the results is not very good because the maximum error in the numerical results is around 10%. Therefore, relying on the case α=0.25 does not allow a good estimate of the original derivative of the Dirac-delta function.The results forα=0.08 is illustrated in Figure 3. It can be clearly seen from Figure 3 that the error in this case is high and that its response is oscillating when the grid-point number is small. On the other hand, this trend in solutions is explained by the inaccurate representation of the regularized derivative of the Dirac-delta function in the discretized model for a few points on the grid. In addition, the percent error converges quickly to zero by augmenting the grid-point number. However, from the numerical results shown in Figure 3, we can also see that the error is influenced by the location of the actuator. We can also observe that the solutions at different points on the beam are of the same order of precision. On the contrary, when comparing the results of Figures 3 and 2, we notice that when higher values of α are utilized in the methodology approach, a more convergent tendency of solutions is achieved, especially for a few grid points. Therefore, when smaller values of α are utilized, a greater precision of the solution is obtained at a relatively larger number of grid points.In Figure4, the error percentage is presented according to the number of grid points for different values of α and the location of the actuator xa. It is observed that when the grid-point number increases, the error in the numerical outcomes tends to a constant value. This constant value is due to the error of the regularization procedure. Also, we can observe clearly that the value of this error can be decreased considerably by decreasing the value by α.Figure 4
Convergence and accuracy of the dimensionless deflection of a simply supported beam under a constant voltageV0 applied by actuator at x/L=1/4 for xa/L=0.3,0.6 (a), xa/L=0.54,0.84 (b), and at x/L=1/3 for xa/L=0.3,0.6 (c), xa/L0.54,0.84 (d) for different values of α.
(a)(b)(c)(d)Figure5 presents the convergence of the dimensionless deflection of the beam normalized by the static deflection wcs=bEbd13L3hp+hV0/2EI of a simply supported beam under a constant voltage V0 applied by an actuator for different locations and with different values of the regularization parameter α. By comparing the results with those of exact solutions, one can conclude that the solutions have good accuracy with small values of the regularization parameter α.Figure 5
Convergence and accuracy of the numerical results with decreasingα of the dimensional deflection of a simply supported beam under a constant voltage V0 applied by actuator for different locations n=51.


### 4.2. Multiharmonic Piezoelectric Actuation

Consider a beam with lengthL equal to 10.16 cm, width b=0.635 cm, thickness h=0.635 cm, Young’s modulus E=2.068×1011 Pa, mass density ρ=10686.9kg/m3. In the first case, the beam is considered simply supported and excited by one piezoelectric actuator with harmonic excitation Ω=9.8695. In the following cases, the beam is considered SS, C-C, and C-F excited by three piezoelectric actuators with different harmonic excitation and the same property. The dynamic response at the center of the beam, wcd, is evaluated for different locations of the piezoelectric actuator and normalized by the static deflection wcs=bEpL3d31hp+h/2EI. The time step used to solve the resulting dynamic equations is Δt=0.001, which is sufficient as the implicit time scheme is used to assure the stability of the algorithm.The results obtained numerically in this problem demonstrated that the value ofα at which the convergence is achieved depends significantly on the length value of the beam. In particular, for small beams, the numerical convergence was obtained with a small specific value of α.One approach to overcome the disadvantages noted above is to express the differential equation governing the beam in a dimensionless form. Afterward, to discretize the resulting nondimensional differential equation using DQM, a procedure similar to that outlined in Section2 is also used. Then, after this manipulation, larger values of α could be used in the proposed procedure. By inserting the dimensionless variable X=x/L, the associated governing differential equation of the dynamic response of the beam can be rewritten from Eq. (1) as follows
(30)∂4wX,t∂X4+η1∂2wX,t∂t2=∑i=1Naη2isinΩitδ′X−Xa1i−δ′X−Xa2i,where
(31)η1=ρhL4EI,η2i=L3bEpid31ihpi+h2EI,Xa1i=xa1iL,Xa2i=xa2iL.Ωi is the harmonic associated to the i-th actuator. The introduction of the regularized Dirac-delta derivative (4) in Eq. (30) leads to
(32)∂4wX,t∂X4+η1∂2wX,t∂t2=∑i=1Naη2isinΩit1α22exp−X−Xa2i22α2×∑n=0M2−14n×12πn!H2nX−Xa2i2α−exp−X−Xa1i22α2∑n=0M2−14n×12πn!H2nX−Xa1i2α.The influence of the regularization valueα on the precision and the consistency of the simulation results for two different actuator locations is illustrated in Figures 6 and 7. For the comparison, analytical solutions developed in the appendix A are used, when only one actuator is considered Na=1. The numerical results obtained are presented in two separate cases: (1) the case where the dimensional differential equation of the beam is taken into account (see Eq. (1)), and (2) the case where the dimensionless differential equation is examined (see Eq. (32)). Based on Figures 6 and 7, the results, of case (1) are not very satisfactory in terms of precision and convergence. It can be seen that the numerical results obtained in the case (2) are more precise than the one obtained in case (1), mainly in terms of numerical convergence. Hence, it should be noted that when the case (2) is taken into account, better accuracy is achieved.Figure 6
Convergence and accuracy of the numerical results for different values ofα for the case where the dimensional differential equation of motion of the beam is discretized using the DQM (n=31).Figure 7
Convergence and accuracy of the numerical results for different values ofα for the case where the dimensionless differential equation of motion of the beam is discretized using the DQM (n=31).On the other hand, Figure8 shows the convergence of solutions in accordance with the number of points in the grid (n). It is observed that the results of the proposed methodology converge uniformly and coincide with the analytical solutions.Figure 8
Convergence of numerical results for normalized central deflection of a simply supported beam excited by a harmonic piezoelectric actuator.Now, three actuators with the same parameters exciting the beam by various harmonic excitation are considered and the piezoelectric actuators are located atxa1/L=0.1,0.3, xa2/L=0.4,0.5, and xa3/L=0.7,0.9 for α=0.062 and n=31. As the presented DQM allows considering various boundary conditions, beams with simply supported (S-S), clamped (C-C), and clamped free (C-F) are investigated. The central displacement wcd for S-S and C-C beam normalized by the static deflection wcs is presented in Figure 9 and the obtained 3D plot is shown the Figure 10 for S-S beam. The time-space responses of C-C are plotted in Figures 11. Also the numerical result for the C-F beam, it is shown in Figures 12 and 13. From the numerical results presented in this section, it can be concluded that the proposed approach is very well suitable for the problem of partial differential equations involving singular functions like the Dirac-delta function and its derivatives. These numerical tests demonstrated that the proposed procedure is reliable and accurate for different boundary conditions.Figure 9
Time response for normalized central deflection of S-S and C-C beam excited by three actuators with different harmonics (Ω1,Ω2,Ω3).Figure 10
Time-space solution of normalized deflection of S-S beam excited by three actuators with different harmonics (Ω1,Ω2,Ω3).Figure 11
Numerical results 3D for normalized deflection of C-C beam excited by three actuators with excitation harmonic different.Figure 12
Numerical results for normalized central deflection of C-F beam excited by three actuators with excitation harmonic different.Figure 13
Numerical results 3D for normalized deflection of C-F beam excited by three actuators with excitation harmonic different.


### 4.3. Impulse Piezoelectric Actuation

This subsection focuses on the dynamic response of beams under various types of piezoelectric impulse excitations. The beam is under a short duration excitation and the maximum response is reached in a very short time. The transient response is considered as well as the permanent one. The used beam parameters are the same to those described in subsection4.2 and the used numerical procedure is based on the DQM and implicit time scheme.


#### 4.3.1. Sinusoidal Piezoelectric Impulse

The piezoelectric actuator excites the beam by an impulse that is described by a half-cycle sinusoidal load. The voltage functionVt is expressed as
(33)Vt=V0sinΩtfor0≤t≤td,0fort>td.td is the time duration and Ω is the excitation frequency. The fundamental period T is T=2πρh/EI, and the time duration td/T is too small (td<<T). Figure 14 shows the normalized displacement wcd/wcs at the middle of the span of a simply supported beam with respect to the normalized time t/T for Ω=π, td=1.5T, the parameter of regularization α=0.062, and the location of actuator fixed at xa/L=0.6,0.9. These results are obtained by the present method and analytical solution giving in the appendix B where wcs=V0bEpL3d31hp+h/2EI.Figure 14
Central displacement of a simply supported beam subjected to a sinusoidal piezoelectric impulse fortd=1.5T, Ω=π, xa/L=0.6,0.9.Figure14 demonstrates that the results of the proposed methodology converge uniformly and the results are identical to analytical ones. It is clearly shown that the maximum amplitude is reached in the time impulse duration (tmax≃0.72T) and the response amplitude wmax/wcs≃9wFree/wcs. This large amplitude may damage the structure and then can be used to design safe microbeams with piezoelectric patches. The associated time-space plot is given in Figure 15 and shows the time-space response behavior.Figure 15
Time-space response of S-S beam subjected to a sinusoidal impulse fortd=1.5T,xa/L=0.6,0.9.


#### 4.3.2. Rectangular Piezoelectric Impulse

In this case, the piezoelectric actuator excites the beam by a rectangular impulse and the excitation functionVt is expressed as
(34)Vt=V0for0≤t≤td,0fort>td.Figure16 shows the normalized displacement w/wcs at the middle of the span of a simply supported beam with respect to the normalized time t/T for Ω=π, parameter of regularization α=0.062, and the location of the actuator is fixed at xa/L=0.6,0.9 the time and space deflection is shown in Figure 17.Figure 16
Central displacement of a simply supported beam subjected to rectangular piezoelectric impulse voltage fortd=1.5T,xa/L=0.6,0.9.Figure 17
Time-space normalized deflection of S-S beam subjected to rectangular piezoelectric impulse voltage fortd=1.5T, xa/L=0.6,0.9.


#### 4.3.3. Symmetric Triangular Piezoelectric Impulse

In this case, the forcing voltageVt is expressed as
(35)Vt=2V0ttdfort≤0,2V0td−ttdfort>td2,0fort>td.These three steps functions can also be written as a single function(36)Vt=2V0ttd−4V0t−td/2tdHt−td2+2V0t−tdtdHt−td,whereHx−a is the Heaviside function: Hx−a=1 for x≤a and Hx−a=0 for x<a. The numerical results for this case are presented in Figure 18 with respect to the normalized time t/T for parameter of regularization α=0.062 and location of the actuator is fixed at xa/L=0.6,0.9. Again, the time-space response is shown in Figure 19.Figure 18
Central displacement of a simply supported beam subjected to symmetric triangular piezoelectric impulse fortd=1.5T, xa/L=0.6,0.9.Figure 19
Time-space normalized deflection of S-S beam subjected to Symmetric triangular piezoelectric impulse fortd=1.5T, xa/L=0.6,0.9.


## 4.1. Beams under a Piezoelectric Actuator with Constant Voltage

In this subsection, we consider a beam with lengthL, Young’s modulus E, moment of area I under one PZT actuator of thickness hp, the piezoelectric strain constant d31, and Young’s modulus Ep with a constant voltage V0. The PZT actuator is assumed to be perfectly bonded and its stiffness is neglected due to its limited contribution to the behavior of the beam, as shown in Figure 1.This problem can be modeled by the differential equation reduced from Eq. (1) as follows
(25)EI∂4wx∂x4=fx,where fx is given by Eq. (6) in which Vt=V0. By using the DOM procedure, Eq. (25) is reduced to the following algebraic equation
(26)KW=F,where K, W, and F are given by
(27)K=EIe4,W=wx1wx2⋯wxnT,F=bEpd31hp+h2V0χ,The vectorχ is given by:
(28)χ=12α2exp−x1−xa222α2∑m=0M2−14m×12πm!H2mx1−xa22α−exp−x1−xa122α2∑m=0M2−14m×12πm!H2mx1−xa12α⋯exp−xn−xa222α2∑m=0M2−14m×12πm!H2mxn−xa22α−exp−xn−xa122α2∑m=0M2−14m×12πm!H2mxn−xa12αT.To ensure the validity of the proposed methodology and its implementation, the problem of a beam simply supported excited by a concentrated actuator is addressed by applying the proposed methodology for various actuator locations and two different values of the regularization parameterα=0.25 and α=0.08, respectively. The error variations used in numerical solutions is defined by
(29)wDQM−wExactwExact×100.In accordance with the number of grid points (n), the error variations for different locations xa1,xa2 of the actuator are shown in Figures 2 and 3.Figure 2
Convergence and accuracy of the dimensionless deflection of a simply supported beam under a constant voltageV0 applied by actuator at various locations for (α=0.25) at x/L=1/2 (a), x/L=1/5 (b), x/L=2/3 (c), and x/L=3/4 (d).
(a)(b)(c)(d)Figure 3
Convergence and accuracy of the dimensionless deflection of a simply supported beam under a constant voltageV0 applied by actuator at different locations for (α=0.08) at x/L=1/2 (a), x/L=1/5 (b), x/L=2/3 (c), and x/L=3/4 (d).
(a)(b)(c)(d)Figure2 shows the variance of the percent error in the numerical results for α=0.25 and for different actuator locations. It is clearly noticed that the results obtained converge uniformly towards their final values. On the other hand, the numerical precision of the results is not very good because the maximum error in the numerical results is around 10%. Therefore, relying on the case α=0.25 does not allow a good estimate of the original derivative of the Dirac-delta function.The results forα=0.08 is illustrated in Figure 3. It can be clearly seen from Figure 3 that the error in this case is high and that its response is oscillating when the grid-point number is small. On the other hand, this trend in solutions is explained by the inaccurate representation of the regularized derivative of the Dirac-delta function in the discretized model for a few points on the grid. In addition, the percent error converges quickly to zero by augmenting the grid-point number. However, from the numerical results shown in Figure 3, we can also see that the error is influenced by the location of the actuator. We can also observe that the solutions at different points on the beam are of the same order of precision. On the contrary, when comparing the results of Figures 3 and 2, we notice that when higher values of α are utilized in the methodology approach, a more convergent tendency of solutions is achieved, especially for a few grid points. Therefore, when smaller values of α are utilized, a greater precision of the solution is obtained at a relatively larger number of grid points.In Figure4, the error percentage is presented according to the number of grid points for different values of α and the location of the actuator xa. It is observed that when the grid-point number increases, the error in the numerical outcomes tends to a constant value. This constant value is due to the error of the regularization procedure. Also, we can observe clearly that the value of this error can be decreased considerably by decreasing the value by α.Figure 4
Convergence and accuracy of the dimensionless deflection of a simply supported beam under a constant voltageV0 applied by actuator at x/L=1/4 for xa/L=0.3,0.6 (a), xa/L=0.54,0.84 (b), and at x/L=1/3 for xa/L=0.3,0.6 (c), xa/L0.54,0.84 (d) for different values of α.
(a)(b)(c)(d)Figure5 presents the convergence of the dimensionless deflection of the beam normalized by the static deflection wcs=bEbd13L3hp+hV0/2EI of a simply supported beam under a constant voltage V0 applied by an actuator for different locations and with different values of the regularization parameter α. By comparing the results with those of exact solutions, one can conclude that the solutions have good accuracy with small values of the regularization parameter α.Figure 5
Convergence and accuracy of the numerical results with decreasingα of the dimensional deflection of a simply supported beam under a constant voltage V0 applied by actuator for different locations n=51.


## 4.2. Multiharmonic Piezoelectric Actuation

Consider a beam with lengthL equal to 10.16 cm, width b=0.635 cm, thickness h=0.635 cm, Young’s modulus E=2.068×1011 Pa, mass density ρ=10686.9kg/m3. In the first case, the beam is considered simply supported and excited by one piezoelectric actuator with harmonic excitation Ω=9.8695. In the following cases, the beam is considered SS, C-C, and C-F excited by three piezoelectric actuators with different harmonic excitation and the same property. The dynamic response at the center of the beam, wcd, is evaluated for different locations of the piezoelectric actuator and normalized by the static deflection wcs=bEpL3d31hp+h/2EI. The time step used to solve the resulting dynamic equations is Δt=0.001, which is sufficient as the implicit time scheme is used to assure the stability of the algorithm.The results obtained numerically in this problem demonstrated that the value ofα at which the convergence is achieved depends significantly on the length value of the beam. In particular, for small beams, the numerical convergence was obtained with a small specific value of α.One approach to overcome the disadvantages noted above is to express the differential equation governing the beam in a dimensionless form. Afterward, to discretize the resulting nondimensional differential equation using DQM, a procedure similar to that outlined in Section2 is also used. Then, after this manipulation, larger values of α could be used in the proposed procedure. By inserting the dimensionless variable X=x/L, the associated governing differential equation of the dynamic response of the beam can be rewritten from Eq. (1) as follows
(30)∂4wX,t∂X4+η1∂2wX,t∂t2=∑i=1Naη2isinΩitδ′X−Xa1i−δ′X−Xa2i,where
(31)η1=ρhL4EI,η2i=L3bEpid31ihpi+h2EI,Xa1i=xa1iL,Xa2i=xa2iL.Ωi is the harmonic associated to the i-th actuator. The introduction of the regularized Dirac-delta derivative (4) in Eq. (30) leads to
(32)∂4wX,t∂X4+η1∂2wX,t∂t2=∑i=1Naη2isinΩit1α22exp−X−Xa2i22α2×∑n=0M2−14n×12πn!H2nX−Xa2i2α−exp−X−Xa1i22α2∑n=0M2−14n×12πn!H2nX−Xa1i2α.The influence of the regularization valueα on the precision and the consistency of the simulation results for two different actuator locations is illustrated in Figures 6 and 7. For the comparison, analytical solutions developed in the appendix A are used, when only one actuator is considered Na=1. The numerical results obtained are presented in two separate cases: (1) the case where the dimensional differential equation of the beam is taken into account (see Eq. (1)), and (2) the case where the dimensionless differential equation is examined (see Eq. (32)). Based on Figures 6 and 7, the results, of case (1) are not very satisfactory in terms of precision and convergence. It can be seen that the numerical results obtained in the case (2) are more precise than the one obtained in case (1), mainly in terms of numerical convergence. Hence, it should be noted that when the case (2) is taken into account, better accuracy is achieved.Figure 6
Convergence and accuracy of the numerical results for different values ofα for the case where the dimensional differential equation of motion of the beam is discretized using the DQM (n=31).Figure 7
Convergence and accuracy of the numerical results for different values ofα for the case where the dimensionless differential equation of motion of the beam is discretized using the DQM (n=31).On the other hand, Figure8 shows the convergence of solutions in accordance with the number of points in the grid (n). It is observed that the results of the proposed methodology converge uniformly and coincide with the analytical solutions.Figure 8
Convergence of numerical results for normalized central deflection of a simply supported beam excited by a harmonic piezoelectric actuator.Now, three actuators with the same parameters exciting the beam by various harmonic excitation are considered and the piezoelectric actuators are located atxa1/L=0.1,0.3, xa2/L=0.4,0.5, and xa3/L=0.7,0.9 for α=0.062 and n=31. As the presented DQM allows considering various boundary conditions, beams with simply supported (S-S), clamped (C-C), and clamped free (C-F) are investigated. The central displacement wcd for S-S and C-C beam normalized by the static deflection wcs is presented in Figure 9 and the obtained 3D plot is shown the Figure 10 for S-S beam. The time-space responses of C-C are plotted in Figures 11. Also the numerical result for the C-F beam, it is shown in Figures 12 and 13. From the numerical results presented in this section, it can be concluded that the proposed approach is very well suitable for the problem of partial differential equations involving singular functions like the Dirac-delta function and its derivatives. These numerical tests demonstrated that the proposed procedure is reliable and accurate for different boundary conditions.Figure 9
Time response for normalized central deflection of S-S and C-C beam excited by three actuators with different harmonics (Ω1,Ω2,Ω3).Figure 10
Time-space solution of normalized deflection of S-S beam excited by three actuators with different harmonics (Ω1,Ω2,Ω3).Figure 11
Numerical results 3D for normalized deflection of C-C beam excited by three actuators with excitation harmonic different.Figure 12
Numerical results for normalized central deflection of C-F beam excited by three actuators with excitation harmonic different.Figure 13
Numerical results 3D for normalized deflection of C-F beam excited by three actuators with excitation harmonic different.


## 4.3. Impulse Piezoelectric Actuation

This subsection focuses on the dynamic response of beams under various types of piezoelectric impulse excitations. The beam is under a short duration excitation and the maximum response is reached in a very short time. The transient response is considered as well as the permanent one. The used beam parameters are the same to those described in subsection4.2 and the used numerical procedure is based on the DQM and implicit time scheme.


### 4.3.1. Sinusoidal Piezoelectric Impulse

The piezoelectric actuator excites the beam by an impulse that is described by a half-cycle sinusoidal load. The voltage functionVt is expressed as
(33)Vt=V0sinΩtfor0≤t≤td,0fort>td.td is the time duration and Ω is the excitation frequency. The fundamental period T is T=2πρh/EI, and the time duration td/T is too small (td<<T). Figure 14 shows the normalized displacement wcd/wcs at the middle of the span of a simply supported beam with respect to the normalized time t/T for Ω=π, td=1.5T, the parameter of regularization α=0.062, and the location of actuator fixed at xa/L=0.6,0.9. These results are obtained by the present method and analytical solution giving in the appendix B where wcs=V0bEpL3d31hp+h/2EI.Figure 14
Central displacement of a simply supported beam subjected to a sinusoidal piezoelectric impulse fortd=1.5T, Ω=π, xa/L=0.6,0.9.Figure14 demonstrates that the results of the proposed methodology converge uniformly and the results are identical to analytical ones. It is clearly shown that the maximum amplitude is reached in the time impulse duration (tmax≃0.72T) and the response amplitude wmax/wcs≃9wFree/wcs. This large amplitude may damage the structure and then can be used to design safe microbeams with piezoelectric patches. The associated time-space plot is given in Figure 15 and shows the time-space response behavior.Figure 15
Time-space response of S-S beam subjected to a sinusoidal impulse fortd=1.5T,xa/L=0.6,0.9.


### 4.3.2. Rectangular Piezoelectric Impulse

In this case, the piezoelectric actuator excites the beam by a rectangular impulse and the excitation functionVt is expressed as
(34)Vt=V0for0≤t≤td,0fort>td.Figure16 shows the normalized displacement w/wcs at the middle of the span of a simply supported beam with respect to the normalized time t/T for Ω=π, parameter of regularization α=0.062, and the location of the actuator is fixed at xa/L=0.6,0.9 the time and space deflection is shown in Figure 17.Figure 16
Central displacement of a simply supported beam subjected to rectangular piezoelectric impulse voltage fortd=1.5T,xa/L=0.6,0.9.Figure 17
Time-space normalized deflection of S-S beam subjected to rectangular piezoelectric impulse voltage fortd=1.5T, xa/L=0.6,0.9.


### 4.3.3. Symmetric Triangular Piezoelectric Impulse

In this case, the forcing voltageVt is expressed as
(35)Vt=2V0ttdfort≤0,2V0td−ttdfort>td2,0fort>td.These three steps functions can also be written as a single function(36)Vt=2V0ttd−4V0t−td/2tdHt−td2+2V0t−tdtdHt−td,whereHx−a is the Heaviside function: Hx−a=1 for x≤a and Hx−a=0 for x<a. The numerical results for this case are presented in Figure 18 with respect to the normalized time t/T for parameter of regularization α=0.062 and location of the actuator is fixed at xa/L=0.6,0.9. Again, the time-space response is shown in Figure 19.Figure 18
Central displacement of a simply supported beam subjected to symmetric triangular piezoelectric impulse fortd=1.5T, xa/L=0.6,0.9.Figure 19
Time-space normalized deflection of S-S beam subjected to Symmetric triangular piezoelectric impulse fortd=1.5T, xa/L=0.6,0.9.


## 4.3.1. Sinusoidal Piezoelectric Impulse

The piezoelectric actuator excites the beam by an impulse that is described by a half-cycle sinusoidal load. The voltage functionVt is expressed as
(33)Vt=V0sinΩtfor0≤t≤td,0fort>td.td is the time duration and Ω is the excitation frequency. The fundamental period T is T=2πρh/EI, and the time duration td/T is too small (td<<T). Figure 14 shows the normalized displacement wcd/wcs at the middle of the span of a simply supported beam with respect to the normalized time t/T for Ω=π, td=1.5T, the parameter of regularization α=0.062, and the location of actuator fixed at xa/L=0.6,0.9. These results are obtained by the present method and analytical solution giving in the appendix B where wcs=V0bEpL3d31hp+h/2EI.Figure 14
Central displacement of a simply supported beam subjected to a sinusoidal piezoelectric impulse fortd=1.5T, Ω=π, xa/L=0.6,0.9.Figure14 demonstrates that the results of the proposed methodology converge uniformly and the results are identical to analytical ones. It is clearly shown that the maximum amplitude is reached in the time impulse duration (tmax≃0.72T) and the response amplitude wmax/wcs≃9wFree/wcs. This large amplitude may damage the structure and then can be used to design safe microbeams with piezoelectric patches. The associated time-space plot is given in Figure 15 and shows the time-space response behavior.Figure 15
Time-space response of S-S beam subjected to a sinusoidal impulse fortd=1.5T,xa/L=0.6,0.9.


## 4.3.2. Rectangular Piezoelectric Impulse

In this case, the piezoelectric actuator excites the beam by a rectangular impulse and the excitation functionVt is expressed as
(34)Vt=V0for0≤t≤td,0fort>td.Figure16 shows the normalized displacement w/wcs at the middle of the span of a simply supported beam with respect to the normalized time t/T for Ω=π, parameter of regularization α=0.062, and the location of the actuator is fixed at xa/L=0.6,0.9 the time and space deflection is shown in Figure 17.Figure 16
Central displacement of a simply supported beam subjected to rectangular piezoelectric impulse voltage fortd=1.5T,xa/L=0.6,0.9.Figure 17
Time-space normalized deflection of S-S beam subjected to rectangular piezoelectric impulse voltage fortd=1.5T, xa/L=0.6,0.9.


## 4.3.3. Symmetric Triangular Piezoelectric Impulse

In this case, the forcing voltageVt is expressed as
(35)Vt=2V0ttdfort≤0,2V0td−ttdfort>td2,0fort>td.These three steps functions can also be written as a single function(36)Vt=2V0ttd−4V0t−td/2tdHt−td2+2V0t−tdtdHt−td,whereHx−a is the Heaviside function: Hx−a=1 for x≤a and Hx−a=0 for x<a. The numerical results for this case are presented in Figure 18 with respect to the normalized time t/T for parameter of regularization α=0.062 and location of the actuator is fixed at xa/L=0.6,0.9. Again, the time-space response is shown in Figure 19.Figure 18
Central displacement of a simply supported beam subjected to symmetric triangular piezoelectric impulse fortd=1.5T, xa/L=0.6,0.9.Figure 19
Time-space normalized deflection of S-S beam subjected to Symmetric triangular piezoelectric impulse fortd=1.5T, xa/L=0.6,0.9.


## 5. Conclusion

There are many various fields of engineering and physics, whose governing partial differential equations with strong singularities like the derivative of the Dirac-delta function. For instance, the behavior of structures under piezoelectric patches can be modelled mathematically by the derivative of the function Dirac-delta. The direct discretization of the derivative of the Dirac-delta function using point discretization techniques like the DQM is not a facile task and special processing is required. In this work, the regularization procedure of the derivative of the Dirac-delta function by the distributed functionals using the Hermite polynomials combined with DQM and a time Implicit scheme was elaborated for the numerical solution of the dynamic behavior of beams with various boundary conditions excited by impulse and harmonic piezoelectric actuators. The DQM combined with regularization procedure was elaborated developed to the numerical solution for space discretization and implicit scheme for time discretization was used. Analytical solutions are also derived for these problems to validate the proposed approach. The location of the actuators is described by the derivatives of the Dirac-delta function that are regularized. The obtained numerical results proved that this proposed methodology is efficient and appropriate. Its main advantage is its simplicity ability to consider an arbitrary number of piezoelectric patches as well its high accuracy. Most importantly, the numerical results reveal that the methodological approach can be used as an efficient tool for many physical and mechanical phenomena modeled by partial differential equation with strong singular coefficients and excitation.



---

*Source: 7391848-2020-08-28.xml*