ANewNumerical Procedure for Vibration Analysis of Beam under Impulse and Multiharmonics Piezoelectric Actuators

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.


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][4][5]. Most applications of this method involve static and dynamic analysis of structural components such as beams and plates [6][7][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.

Equilibrium Equation.
In this work, a beam with the length L, cross-section area A, Young's modulus E, moment of area I, and density ρ excited by N a 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] EI where in which wðx, 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. N a is the number of the piezoelectric actuators. M a i is the force moment induced by the piezoelectric actuator i, and d i 31 , E i p , and h i p are the piezoelectric strain constant, Young's modulus, thickness of the i-th actuator, respectively. δ ′ ð:Þ is the derivative of the 2 Journal of Applied Mathematics Dirac-delta function, x i a1 , x i a2 are the coordinates of the two ends of i-th actuator and V i ðtÞ is the voltage applied on the i-th actuator.
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 H 2n as follows For a very small α, the function δ M ðx − x 0 Þ becomes identical to the functional delta when the degree of polynomial M is fixed.
In addition, the derivatives of the regularized formulation give an approximation to the derivatives of the Dirac-delta function and are given by 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 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.

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 function wðxÞ is sufficiently smooth over the interval ½x 1 , x N . In this interval, N distinct nodes are defined.
The function values on these points are assumed to be According to the DQM, the first and second-order derivatives on each of these nodes are given by [29].
The coefficients a ij and b ij are the weighting coefficients of the first and second-order derivatives with 3 Journal of Applied Mathematics respect to x, respectively. The coefficients a ij and b ij are given as follows [29].
Similarly, we can obtain higher-order derivative formulas by using the higher weighted coefficients, which are expressed in e ðmÞ ij to avoid confusion. They are characterized by induction [4].
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] x 3.2. Numerical Accuracy of the DQM. Consider a function wðxÞ, 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 x i can be obtained as [29]. where 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.

DQM Analogues.
For numerical solution, n grid points with coordinates x 1 , x 2 , ⋯, x n in the x-direction are considered. Applying the quadrature rule, to Eq. (1), the following ordinary differential system is obtained This system can be rewritten in the following matrix from where ½M and ½K are the resulting mass and stiffness matrices of the beam, fWðtÞg and f € WðtÞg are the displacement and acceleration vectors, and fFðtÞg is the charge vector applied by the piezoelectric actuators.
These terms are given by where ½I is the identity matrix of order n × n and ½e ð4Þ is the DQM weighting coefficient matrix of the fourth-order derivative.

Journal of Applied Mathematics
The discrete classical boundary conditions of the beam at x = 0 and x = 1, using the DQ method, can be written as following: where n 0 and n l can be chosen from 1, 2, or 3. The selection of the values n 0 and n l can have the following traditional boundary conditions [6]:  26 28 30 x a /L = (0.6,0.9)

Journal of Applied Mathematics
Similarly, the boundary conditions (19) can also be expressed in a matrix form 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 ½K BB and ½K BI are 4 × 4 and 4 × ðn − 4Þ, respectively.
Implementing the boundary conditions into Eq. (16) leads to the following system of ordinary differential equa-tionsK   26 28 30 x a /L = (0.6,0.9)  Journal of Applied Mathematics Equation (22) can be solved by using various step-bystep time integration schemes. In this study, the time derivative is discretized using the centered finite difference scheme, then Substituting this approximate expression acceleration into (22), one gets:K

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 voltage V 0 is considered to demonstrate the  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 where f ðxÞ is given by Eq. (6) in which VðtÞ = V 0 . By using the DOM procedure, Eq. (25) is reduced to the following algebraic equation where ½K, fWg, and fFg are given by The vector χ is given by: 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 In accordance with the number of grid points (n), the error variations for different locations ðx a1 , x a2 Þ of the actuator are shown in Figures 2 and 3. Figure 2 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 gridpoint 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.  In Figure 4, the error percentage is presented according to the number of grid points for different values of α and the location of the actuator x a . 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 5 presents the convergence of the dimensionless deflection of the beam normalized by the static deflection 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, w cd , is evaluated for different locations of the piezoelectric actuator and normalized by the static deflection w cs = ðbE p L 3 d 31 ðhp + 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 Section 2 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 where ð31Þ Ω i is the harmonic associated to the i-th actuator. The introduction of the regularized Dirac-delta derivative (4) in Eq. (30) leads to 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 ðN a = 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)).

Journal of Applied Mathematics
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.
On the other hand, Figure 8 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. Now, three actuators with the same parameters exciting the beam by various harmonic excitation are considered and the piezoelectric actuators are located at x 1 a /L = ð0:1,0:3Þ, x 2 a /L = ð0:4,0:5Þ, and x 3 a /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 w cd for S-S and C-C beam normalized by the static deflection w cs 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.
0 for t > t d : ( ð33Þ t d is the time duration and Ω is the excitation frequency. The fundamental period T is T = 2π ffiffiffiffiffiffiffiffiffiffiffi ffi ρh/EI p , and the time duration t d /T is too small (t d < <T). Figure 14 shows the normalized displacement w cd /w cs at the middle of the span of a simply supported beam with respect to the normalized time t/T for Ω = π, t d = 1:5T, the parameter of regularization α = 0:062, and the location of actuator fixed at x a /L = ð0:6,0:9Þ. These results are obtained by the present method and analytical solution giving in the appendix B where w cs = ðV 0 bE p L 3 d 31 ðhp + hÞÞ/ð2EIÞ. Figure 14 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 (t max ≃ 0:72T) and the response amplitude w max /w cs ≃ 9ðw Free /w cs Þ. This large amplitude may damage the structure and then can be used to design safe  14 Journal of Applied Mathematics microbeams with piezoelectric patches. The associated time-space plot is given in Figure 15 and shows the time-space response behavior.

Rectangular Piezoelectric Impulse.
In this case, the piezoelectric actuator excites the beam by a rectangular impulse and the excitation function VðtÞ is expressed as Figure 16 shows the normalized displacement w/w cs 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 x a /L = ð0:6,0:9Þ the time and space deflection is shown in Figure 17 0 for t > t d : These three steps functions can also be written as a single function where Hðx − aÞ is the Heaviside function: Hðx − aÞ = 1 for x ≤ a and Hðx − 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 x a /L = ð0:6,0:9Þ. Again, the time-space response is shown in Figure 19.

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. The following algebraic equation is resulted.

Data Availability
The code source data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.