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.
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. Bahan [12] implemented for the numerical solution of nonlinear Kawahara equation via the Crank-Nicolson-Differential Quadrature Method based on modified cubic B-splines and Bahan et al. [13] applied DQM based on modified cubic B-splines for numerical solution of the complex modified Korteweg-de Vries equation. Also, Bahan [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 length , cross-section area , Young’s modulus , moment of area , and density excited by 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] where in which is the transverse displacement of the beam, is the axial coordinate, is the time, and are the width and thickness of the beam, respectively. is the number of the piezoelectric actuators. is the force moment induced by the piezoelectric actuator , and , , and are the piezoelectric strain constant, Young’s modulus, thickness of the -th actuator, respectively. is the derivative of the Dirac-delta function, , are the coordinates of the two ends of -th actuator and is the voltage applied on the -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 as follows
For a very small , the function becomes identical to the functional delta when the degree of polynomial 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 .
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.
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 function is sufficiently smooth over the interval . In this interval, 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 and are the weighting coefficients of the first and second-order derivatives with respect to , respectively. The coefficients and are given as follows [29].
Similarly, we can obtain higher-order derivative formulas by using the higher weighted coefficients, which are expressed in 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]
3.2. Numerical Accuracy of the DQM
Consider a function , which is approximated by the Lagrange interpolation polynomial of degree . The error for the order derivative approximation of this function at point can be obtained as [29]. where
As a result, form equation (13), the accuracy of the DQM can be directly proportional to . This means that a result can be achieved with very high accuracy even if the number of grid points is small.
3.3. DQM Analogues
For numerical solution, grid points with coordinates in the -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 and are the resulting mass and stiffness matrices of the beam, and are the displacement and acceleration vectors, and is the charge vector applied by the piezoelectric actuators.
These terms are given by where is the identity matrix of order and is the DQM weighting coefficient matrix of the fourth-order derivative.
The discrete classical boundary conditions of the beam at and , using the DQ method, can be written as following: where and can be chosen from 1, 2, or 3. The selection of the values and can have the following traditional boundary conditions [6]: (i)simply supported: ; (ii)clamped-clamped: ; (iii)clamped-simply supported: ; (iv)clamped-free: ; (v)free-free: ;
Similarly, the boundary conditions (19) can also be expressed in a matrix form where the subscripts and 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 and are and , respectively.
Implementing the boundary conditions into Eq. (16) leads to the following system of ordinary differential equations wherein
and ; where , , and are matrices of order , , and , 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
Substituting this approximate expression acceleration into (22), one gets:
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 voltage 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 is taken as 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 length , Young’s modulus , moment of area under one PZT actuator of thickness , the piezoelectric strain constant , and Young’s modulus with a constant voltage . 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 is given by Eq. (6) in which . By using the DOM procedure, Eq. (25) is reduced to the following algebraic equation where , , and 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 and , respectively. The error variations used in numerical solutions is defined by
In accordance with the number of grid points (), the error variations for different locations of the actuator are shown in Figures 2 and 3.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
Figure 2 shows the variance of the percent error in the numerical results for 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 . Therefore, relying on the case does not allow a good estimate of the original derivative of the Dirac-delta function.
The results for 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 Figure 4, the error percentage is presented according to the number of grid points for different values of and the location of the actuator . 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 .
(a)
(b)
(c)
(d)
Figure 5 presents the convergence of the dimensionless deflection of the beam normalized by the static deflection of a simply supported beam under a constant voltage 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 .
4.2. Multiharmonic Piezoelectric Actuation
Consider a beam with length equal to cm, width cm, thickness cm, Young’s modulus Pa, mass density kg/m3. In the first case, the beam is considered simply supported and excited by one piezoelectric actuator with harmonic excitation . 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, , is evaluated for different locations of the piezoelectric actuator and normalized by the static deflection . The time step used to solve the resulting dynamic equations is , 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 , the associated governing differential equation of the dynamic response of the beam can be rewritten from Eq. (1) as follows where
is the harmonic associated to the -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 . 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.
On the other hand, Figure 8 shows the convergence of solutions in accordance with the number of points in the grid (). 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 , , and for and . 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 for S-S and C-C beam normalized by the static deflection 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.
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 subsection 4.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 function is expressed as is the time duration and is the excitation frequency. The fundamental period is , and the time duration is too small (). Figure 14 shows the normalized displacement at the middle of the span of a simply supported beam with respect to the normalized time for , , the parameter of regularization , and the location of actuator fixed at . These results are obtained by the present method and analytical solution giving in the appendix B where .
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 () and the response amplitude . 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.
4.3.2. Rectangular Piezoelectric Impulse
In this case, the piezoelectric actuator excites the beam by a rectangular impulse and the excitation function is expressed as
Figure 16 shows the normalized displacement at the middle of the span of a simply supported beam with respect to the normalized time for , parameter of regularization , and the location of the actuator is fixed at the time and space deflection is shown in Figure 17.
4.3.3. Symmetric Triangular Piezoelectric Impulse
In this case, the forcing voltage is expressed as
These three steps functions can also be written as a single function
where is the Heaviside function: for and for . The numerical results for this case are presented in Figure 18 with respect to the normalized time for parameter of regularization and location of the actuator is fixed at . Again, the time-space response is shown in Figure 19.
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.
Appendix
A.Analytical Solution for S-S Beam Exited by One Piezoelectric Harmonic Actuator
Under the piezoelectric actuator, the motion governing equation can be expressed as where
For a simply supported beam, the solution can be expressed by the following Fourier series:
Substituting the solution into Equation (A.5), multiplying both sides
with and using the mode orthogonality, namely,
The following algebraic equation is resulted.
Letting one obtains as where and .
Thus, the displacement response can be obtained by substituting Equation (A.5) into Equation (A.2) as follows:
B.Vibration of S-S Beam Exited by a Sinusoidal Piezoelectric Impulse
The motion governing equation of this problem can be expressed as where
The analysis is organized in two phases.
B.1.Phase I
In this phase the excitation is harmonic and the response includes the transient and permanent solutions, and is given by equation (A.6)
B.2.Phase II
During this phase, the system starting at , the system is in free vibration with initial conditions and at the end of phase I.
Assuming the solution in the form of mode superposition, the transverse deflection of a free simply supported beam can be written as
Substituting the solution into equation (B.10) for phase II, we get where
The general solution of equation (B.12) can be written as:
The constants of integration and can be determined knowing the initial displacement and the initial velocity at time . Substituting these initial conditions into equation (B.7) and equation (B.8) one gets
Equation (B.7) can then be written as
Thus, the free vibration displacement of beam for can be obtained by substituting Equation (B.10) into Equation (B.4): where
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.
Acknowledgments
The authors would like to acknowledge the financial support of the CNRST and the Moroccan Ministry of Higher Education and Scientific Research with the project PPR2/06/2016, as well as to the DSR at the King Abdulaziz University, Jeddah, Saudi Arabia.