A Differential Quadrature Procedure with Regularization of the Dirac-delta Function for Numerical Solution of Moving Load Problem

The differential quadrature method (DQM) is one of the most elegant and efficient methods for the numerical solution of partial differential equations arising in engineering and applied sciences. It is simple to use and also straightforward to implement. However, the DQM is well-known to have some difficulty when applied to partial differential equations involving singular functions like the Dirac-delta function. This is caused by the fact that the Dirac-delta function cannot be directly discretized by the DQM. To overcome this difficulty, this paper presents a simple differential quadrature procedure in which the Dirac-delta function is replaced by regularized smooth functions. By regularizing the Dirac-delta function, such singular function is treated as non-singular functions and can be easily and directly discretized using the DQM. To demonstrate the applicability and reliability of the proposed method, it is applied here to solve some moving load problems of beams and rectangular plates, where the location of the moving load is described by a time-dependent Dirac-delta function. The results generated by the proposed method are compared with analytical and numerical results available in the literature. Numerical results reveal that the proposed method can be used as an efficient tool for dynamic analysis of beamand plate-type structures traversed by moving dynamic loads.


INTRODUCTION
Most practical problems in engineering are governed by partial differential equations with proper boundary conditions.In general, two kinds of methods, i.e., analytical and numerical methods, can be used to solve engineering problems.Analytical methods are often preferred since they allow better information of the solution over the problem domain.However, analytical methods are often limited to simple engineering problems (i.e., problems with simple geometry, simple boundary conditions, etc.).To overcome the limitations of the analytical methods, one may resort to numerical methods for finding approximately the required solution.The numerical methods, in general, can be very useful since they can be applied to more realistic engineering problems with complex boundary conditions and irregular geometries.Among them, the finite difference and finite element methods have been extensively used by many researchers to obtain approximate solutions for the engineering problems.However, these methods generally converge slowly with respect to mesh refinement and are extremely expensive for achieving high precision.
The high order methods can be adopted to do spatial discretization to tackle the difficulties of the low order finite difference and finite element methods.Among them, the differential quadrature method (DQM) has attracted many attentions because of many advantages such as high accuracy, simplicity and efficiency.The DQM, which is in the family of collocation and finite difference methods, was firstly introduced by Bellman and his associates (1971,1972) in the early 1970s for the numerical solution of initial and boundary value problems.Since then, the DQM has been successfully applied to a variety of engineering problems (Bert and Malik, 1996a;Shu, 2000, Zong andZhang, 2009).Most of these applications are related to static and dynamic analysis of structural components like beams, plates, and shells (Wang et al., 1993;Bert and Malik, 1996b;Shu, 1996;Civalek, 2004).Newer applications include the use of the DQM for approximation of time-derivative terms (Eftekhari and Jafari, 2012a), and solving fluid-structure interaction problems (Eftekhari and Jafari, 2014).The results of many researchers reveal that the DQM is computationally efficient and is applicable to a large class of initial and/or boundary value problems.
In spite of above-mentioned favorable features, the DQM has some difficulty in handling partial differential equations involving singular functions like the Dirac-delta function.This is mainly caused by the particular characteristics of the Dirac-delta function.For example, the one dimensional Dirac-delta function has the following properties: where f(ξ) may be any arbitrary function.It can be seen from Eqs.
(1) to (3) that two properties of the Dirac-delta function are in the form of integrals.It is well known that the weak-form-based methods such as the Ritz and finite element methods can easily handle the problems involving the Dirac-delta function, since they directly integrate the governing differential equation of the problem.However, the strong-form-based methods such as the DQM cannot easily handle such type of problems since they directly discretize the governing differential equation of the problem.A way for overcoming this issue is to combine the DQM with the integral quadrature (IQ) method (Eftekhari and Jafari, 2012b) so that the properties of the Dirac-delta function are approximately satisfied.However, as discussed earlier by Eftekhari and Jafari (2012c), this approach may encounter some difficulties when applied to the problems involving the time-dependent Dirac-delta function.
The above mentioned difficulty has also been addressed earlier by many researchers when similar methods such as the collocation and finite difference methods are used.In Gheorghiu (2007), a discussion is given that the difficulty of the collocation method in handling such singular functions is caused by the Gibbs phenomenon in which the numerical solutions become oscillatory around singularities.To obviate this difficulty, it has been suggested that the singular function (here, the Diracdelta function) should be regularized in order to achieve a smoother representation of the singular function and to stabilize the unwanted oscillatory behavior of solutions near the singularities.Various regularization techniques have been developed in Waldén (1999), Wei et al. (2002), Engquist et al. (2005), Burko and Khanna (2007), Sundararajan et al. (2007), Towers (2007), and Rivera et al. (2013).Waldén (1999) solved some elliptic and parabolic problems using the finite difference method and showed that the right regularization of the singular source terms can improve the convergence rate of solutions.In Wei et al. (2002), various forms of the regularized Dirac-delta function have been proposed in the context of the discrete singular convolution method.Engquist et al. (2005) proposed a level-set finite difference approach for regularization of the Dirac-delta function.Burko and Khanna (2007) proposed a very narrow Gaussian function for regularization of the Dirac-delta function.Recently, Sundararajan et al. (2007) proposed various advanced approaches for modeling the Dirac-delta function.More recently, Towers (2007) proposed two closely related finite difference methods for discretization of the Dirac-delta function.Most recently, Rivera et al. (2013) proposed various smooth functions for approximation of the Dirac-delta function.
Alternatively, Jung (2009), Jung and Don (2009) also Jung et al. (2009) made use of an alternative definition of the Dirac-delta function that the derivative of the Heaviside function, H(ξ), is the Dirac-delta function, δ(ξ), in the distribution sense, namely, dH(ξ)/dξ = δ(ξ).The method was referred to as the direct projection method.It was shown that although the direction projection of the Heaviside function yields highly oscillatory approximation of the Dirac-delta function, it can yield a non-oscillatory approximation of the solution and the error can also decay uniformly for certain types of differential equations.Most recently, Wang et al. (2014) investigated at length the properties of various ways for approximation of the Dirac-delta function.The techniques considered in their study were the regularization using the Gaussian smooth function, the points source approach, the direct projection technique, and the domain decomposition method.They discretized the differential equation using the finite difference and pseudo-spectral methods and concluded that the domain decomposition method is the best among these methods at highest accuracy for solving sine-Gordon equation involving the Dirac-delta function.
From the above literature survey, it is observed that a significant amount of research has been done on discretization of the Dirac-delta function using the collocation and finite difference methods.Only a limited number of research studies have addressed this subject when using the DQM.As pointed out earlier, the present author and his co-worker have recently made some attempts to solve such problems (Eftekhari and Jafari, 2012b;2012c).In Eftekhari and Jafari (2012c) the capability of the DQM for handling the time-dependent Dirac-delta function was discussed at length and it was concluded that there are some major limitations in the choice of time steps and mesh sizes.Most recently, Nikkhoo et al. (2012), Nikkhoo and Kananipour (2014) also Kananipour et al. (2014) have proposed a formulation based on the DQM for transient dynamic analysis of curved beams traversed by moving concentrated loads.But, no details were given how the time-dependent Diracdelta function was approximated or treated numerically in their study.Their numerical solutions also were not in a good agreement with the analytical and finite element solutions (Nikkhoo et al., 2012;Nikkhoo and Kananipour, 2014).The main source of error in their numerical results may be due to an incorrect approximate modeling of the Dirac-delta function.
In this study, a differential quadrature procedure based on the regularization of the Dirac-delta function is presented for the numerical solution of the moving load problem.Since the Dirac-delta function is replaced by a regularized smooth function, the DQM can be easily and directly applied to discretize the governing equation of the problem.To demonstrate its applicability and reliability, the method is applied here to solve some moving load problems of beams and rectangular plates, where the location of the moving load is described by a time-dependent Dirac-delta function.To the best of our knowledge, this is the first attempt in applying the DQM to moving load problems of straight beams and rectangular plates.Numerical results are presented and compared with analytical and numerical ones available in the literature.The numerical results prove that the proposed method is reliable and accurate and can be used as an efficient tool for handling the moving load problem.

DIFFERENTIAL QUADRATURE METHOD
Let w(ξ) be an arbitrary function and ξ 1 , ξ 2 , …, ξ n be a set of grid points in the ξ-direction.According to the DQM, the rth-order derivative of the function w(ξ) at any grid point can be approximated by the following formulation where w(ξ j ) represents the functional value at a grid point ξ j , w (r) (ξ i ) indicates the rth-order derivative of w(ξ) at a grid point ξ i , and ) (r ij A are the weighting coefficients of the rth-order derivative.The first-order weighting coefficients, i.e., ) 1 ( ij A , can be obtained from the following explicit formula (Quan and Chang, 1989) where The weighting coefficients of the higher-order derivatives (r > 1) can be computed from the following recurrence relationship (Shu, 2000) One of the key factors in the accuracy of the DQ solutions is the choice of grid points.The zeros of some orthogonal polynomials are commonly adopted as the grid points.In this work, the DQM grid points are taken nonuniformly spaced and are given by the following equations (Bert and Malik, 1996a) where Δ is a parameter that determines the closeness between the boundary points (ξ 1 and ξ n ) and their immediate adjacent points (ξ 2 and ξ n-1 ).In practice, the magnitude of Δ should be as small as possible (≤ 10 -3 ).In this study, the magnitude of parameter Δ is assumed to be Δ = 10 -3 .

Numerical Accuracy of the DQM
Let the nth-order derivative of the function w(ξ) be a constant (say C).It can be shown that the error estimated for Eq. ( 4) is (Shu, 2000) !
where ) (  is defined in Eq. ( 6).It may be seen from Eq. ( 9) that the accuracy of the DQM may be proportional to n or its powers.This implies that very high accuracy can be achieved using the DQM even if the number of grid points, n, is rather small.

REGULARIZATION OF THE DIRAC-DELTA FUNCTION
As pointed out earlier in introduction, due to particular characteristics of the Dirac-delta function, the direct discretization of this function using the DQ method is not an easy task.A way for overcoming this issue is to approximate (or regularize) the Dirac-delta function with simple mathematical functions.In this regard, various forms of the regularized Dirac-delta function have been proposed in the literature (for example, see Wei et al. (2002)).Among them, the following form of the regularized Dirac-delta function has been commonly used due to its excellent numerical properties such as smoothness and regularity: where α is the regularization parameter that describes the relationship between smoothness and the desired accuracy of approximation.To obtain an accurate representation of the Dirac-delta function using above regularized function, the parameter α should be as small as possible (α → 0). Figure 1 shows the effect of α-value on accuracy and smoothness of the approximate function.
It can be clearly seen from Figure 1 that the accuracy and smoothness of the approximate function can be easily and efficiently adjusted by changing the value of the parameter α.In other words, by decreasing the magnitude of the parameter α, the shape of approximate function becomes nonsmoother and approaches to that of the real Dirac-delta function.It is noted that although decreasing the value of the parameter α results in a better representation of the Dirac-delta function, it can introduce another difficulty when such function is approximated by the point discretization methods like the DQM.The reason for this is that a higher-order DQM should be used in these cases and this may considerably increase the computational cost especially when the α-value is chosen to be too small.Therefore, when a point discretization method like the DQM is employed, one should find a proper value for the α such that the numerical accuracy and computational efficiency of the method are balanced.

FORMULATION FOR FORCED VIBRATION OF BEAMS CARRYING MOVING LOADS
Consider an Euler-Bernoulli beam with length L, mass per unit length ρA, and flexural stiffness EI subjected to a concentrated load f moving at a constant velocity v.The governing differential equation of motion of the beam is given by (Meirovitch, 1967;Fryba, 1999) where w(x, t) is the lateral deflection of the beam, t is the time, and x f (t) is position coordinate of the moving load.In view of Eq. ( 10), Eq. ( 11) can be rewritten as Consider n grid points with coordinates x 1 , x 2 , …, x n in the x-direction.Satisfying Eq. ( 12) at any grid point x = x i and substituting the quadrature rule, given in Eq. ( 4), into results gives where ] [M and ] [K are the mass and stiffness matrices of the beam, are the displacement and acceleration vectors, and where [I] is an identity matrix of order n × n and [A ] (4) is the DQM weighting coefficient matrix of the fourth-order derivative.After applying the boundary conditions of the beam, Eq. ( 13) can be solved for the unknown grid-point values using various step-by-step time integration schemes.The details of implementation of the beam boundary conditions can be found in Bert and Malik (1996a).
In this study, the Newmark scheme (Bathe and Wilson, 1976) is employed to solve system (13).

FORMULATION FOR FORCED VIBRATION OF RECTANGULAR PLATES CARRYING MOVING LOADS
Consider an isotropic thin rectangular plate with length a, width b, mass per unit area ρh, and flexural rigidity D subjected to a moving concentrated load f.The governing differential equation of motion of the rectangular plate is given by (Rao, 2007) where w(x, y, t) is the lateral deflection of the plate, x f (t) and y f (t) are position coordinates of the moving load in the x-and y-directions, respectively, and v x and v y are the corresponding moving speeds.In view of Eq. ( 10), Eq. ( 18) can be rewritten as Consider n grid points with coordinates x 1 , x 2 , …, x n in the x-direction, and m grid points with coordinates y 1 , y 2 , …, y m in the y-direction.Satisfying Eq. ( 19) at any grid point x = x i and substituting the quadrature rule, given in Eq. ( 4), into results gives where [I] is an identity matrix of order n × n, [A ] (2) and [A ] (4) are the DQM weighting coefficient matrices of the second-and fourth-order x-derivatives, respectively.Furthermore, Satisfying Eq. ( 20) at any grid point y = y i and substituting the quadrature rule, given in Eq. ( 4), into results gives After applying the plate boundary conditions, one can solve the resulting system of ordinary differential equations for unknowns using various time marching schemes.The details of implementation of the plate boundary conditions can be found in Malik and Bert (1996) also in Shu andDu (1997a, 1997b).

NUMERICAL RESULTS AND DISCUSSION
To validate the proposed formulation and its implementation, a number of numerical examples are now presented.The first and second examples demonstrate the applications of the proposed method to static analysis of simply supported beams and rectangular plates subjected to concentrated point loads.These numerical examples are presented herein to better understand the main source of errors in the method and to better verify the accuracy and convergence of the method.The third example is of moving load problem of beams while the last one is that of rectangular plates.Two boundary conditions, namely, fully simply supported and fully clamped end conditions are considered in the third and last examples.For beams and plates with simply supported end conditions, analytical solutions can be found in the literature (Timoshenko and Woinowsky-Krieger, 1959;Fryba, 1999).Therefore, the accuracy of the proposed method can be verified by comparing the calculated results with those of analytical ones.

Deflection of a Simply Supported Beam Due to a Concentrated Point Load
Consider the bending problem of a simply supported beam with length L and flexural stiffness EI subjected to a concentrated transverse load f acting at a point x = x f , as shown in Figure 2. The governing differential equation for the bending of the beam can be written from Eqs. ( 11) and ( 12) as The procedure for discretizing Eq. ( 30) using the DQM is very similar to that described in Section 4. By applying the DQM to Eq. ( 30) we obtain where the stiffness matrix ] [K is defined already in Eq. ( 14), and It may be seen from Eqs. ( 33) and ( 34) that the elements of the load vector } {F depend directly on the coordinates of the grid points.This implies that the coordinates of the grid points can directly affect the values of the load vector.We will show in this section that the accuracy of the solutions may also be influenced by the grid points distribution.To demonstrate this, we computed the elements of the vector } {δ using Eqs.( 34) and ( 8) for different values of x f (location of the applied load), n (number of grid points) and α (regularization parameter), and compared the results with those obtained using Eq. ( 10).Figures 3 and 4 illustrate the results.These results are obtained using two different values of α, namely, α = 0.25 and α = 0.05.
When α = 0.25, Figure 3 compares the results of discrete model with those of the continuous model.It can be seen that, except for n = 9, the discrete model shows a good representation to the continuous model.Note that the peak point of the discrete model with n = 9 differs slightly from that of the continuous model for x f / L = 0.17 and x f / L = 0.34.This is because the location of the applied load does not coincide with the DQM grid points in these cases.
Figure 4 presents the results for α = 0.05.It can be seen from Figure 4 that when the number of grid points is small, the discretized regularized Dirac-delta function does not provide a good repre-x f x f sentation to the original continuous function.Most importantly, the nonzero zone of the original continuous function is not accurately represented in the discretized model when the grid number is small.However, by increasing the grid number, the nonzero zone is better represented in the discretized model.On the other hand, by comparing the results of Figure 4 with those of Figure 3, one may conclude that as the parameter α increases, the discretized model provide better representation to the original continuous model.34)) with the continuous regularized Dirac-delta function (see Eq. ( 10)) for different number of grid points and locations of applied load (α = 0.05).
As we discussed earlier in Section 3, the inappropriate choice of parameter α may be a source of error in the numerical results.This error, however, can be easily controlled by changing the value of the regularization parameter α.On the other hand, the numerical results presented in Figures 3 and  4 show that there is another source of error (say the main source of error) which may significantly influence the accuracy of the solutions.This error is mainly due to inappropriate representation of the regularized Dirac-delta function in the discretized form.
First, from Eq. ( 9) we note that the error of the DQM for approximation of w(ξ) and its derivatives should decrease exponentially to zero due to the presence of nth factorial in the denominator.Therefore, we expect that very high accuracy to be achieved when using a significantly small number of grid points.However, when the DQM is applied to the present problem, we see from Figures 3 and 4 that the regularized Dirac-delta function is not appropriately represented in the discretized model for small number of grid points and this may introduce another source of error in the numerical results.As it can be seen from Figures 3 and 4, this error is mainly caused by the coordinates of the grid points and is more significant for smaller values of α.A way for overcoming this difficulty is to change the coordinates of grid points such that more grid points are clustered in the nonzero zone, for instance by stretching the grid points toward the location of the applied load.However, our numerical experiments showed that this may results in ill-conditioning of the DQM resultant matrices (to save the space, these results are not presented herein).
To demonstrate the convergence and accuracy of the proposed approach, the bending problem of a simply supported beam under a concentrated point load is solved using the proposed approach for different locations of the applied load (x f ) and for two different values of the regularization parameter (α).Figures 5 and 6 present the variations of the percent error in numerical solutions (defined as 100 ) with n for different values of x f .
The numerical results for α = 0.25 are shown in Figure 5.As it can be seen, the obtained solutions converge to their final values uniformly.However, the numerical accuracy of the solutions is not very satisfactory as the maximum error in the numerical results is about 17 %.This is because the case α = 0.25 does not provide a good approximation to the original Dirac-delta function (see Figure 1 for more details).It can also be seen from Figure 5 that, in most cases, the error in numerical results increases as the value of x f increases.But, the error in different points of the beam (i.e., x/L = 1/5, 1/2, 2/3, and 3/4) is almost the same order.Figure 6 illustrates the results for α = 0.05.It can be seen from Figure 6 that when the number of grid points is too small, the error in quadrature solutions is very high and its behavior is oscillatory.Particularly, erroneous results are obtained when n = 9, 11, 13, and 15.As we discussed earlier, the reason for such trend of solutions is the inappropriate representation of the regularized Dirac-delta function in the discretized model for small number of grid points (see Figure 4 for details).However, the error in quadrature solutions approaches to zero rapidly by increasing the number of grid points.On the other hand, from the numerical results presented in Figure 6, one sees that the error in quadrature solutions is also influenced by the location of the applied load, as to be expected.From Figure 6, one also sees that the quadrature solutions at different points of the beam (i.e., x/L = 1/5, 1/2, 2/3, and 3/4) are of the same order of accuracy.On the other hand, by comparing the results of Figure 6 with those of Figure 5, one sees that when larger values of α are used in the method, then better converging trend of solutions is achieved particularly at small number of grid points.However, when smaller values of α are used in the method, then better accuracy of the solutions is achieved at rather large number of grid points.
In Figure 7, the percent error in numerical results is plotted against the number of grid points for different values of α and x f .It can be seen that the error in numerical results converges to a constant nonzero value by increasing the number of grid points.This nonzero constant value is, in fact, the error of the regularization function (see Eq. ( 10)).It is clear from Figure 7 that the value of this error can be decreased significantly by decreasing the value of α.
To better see the effects of regularization parameter on accuracy and convergence of solutions, the relative error in numerical results (defined as ) is plotted against α in Figure 8 for different values of x f .These results are obtained using a sufficiently large number of grid points (n = 45) and are given at different points of the beam.Note that the values of α in the plot are in the range 0.05 to 0.5.As it can be seen, depending on location of the applied load, the overall order of convergence of the method with α is in the range 1 to 2. This implies that the employed regularized Dirac-delta function is at most second-order accurate.On the other hand, the results of Figure 8 show that the order of convergence of the method (with respect to α) depends on the location of the applied load, and increases as x f → 0.5.

Deflection of a Simply Supported Rectangular Plate Due to a Concentrated Point Load
Consider the bending problem of a simply supported rectangular plate with length a, width b, thickness h, and bending stiffness D subjected to a concentrated transverse load f acting at a point (x, y) = (x f , y f ).The governing differential equation for the bending of the plate can be written from Eqs. ( 18) and ( 19) as The procedure for discretizing Eq. ( 35) using the DQM is very similar to that described in Section 5 and is not repeated here again.
Our numerical results for the present problem showed that the error in numerical results may have similar trend as that shown in Figures 5 and 6 for the beam problem.Most importantly, the convergence problem of the method for the case of small values of α and n has been found to be more critical than that of the beam problem due to the existence of the double Dirac-delta function in the governing differential equation of the plate.This can be clearly observed from the numerical results presented in Figure 9.The effects of regularization parameter value on accuracy and convergence of the solutions are studied in Figure 10.Again, one sees that the regularization function ( 10) is at most second-order accurate.

Vibration of a Beam Due to a Moving Point Load
Consider a simply supported beam with a span L of 10.16 cm, width b = 0.635 cm, thickness h = 0.635 cm, modulus of elasticity E = 2.068 × 10 11 Pa, mass density ρ = 10686.9kg/m 3 , subjected to a f = 4.45 N moving force (Rieker et al., 1996).The dynamic responses at the beam center, w cd , are evaluated for different values of v (velocity of moving load) and normalized by the static deflection w cs = fL 3 /(48EI) of a simply supported beam under a point load f at mid-span.The Newmark method with 1000 time intervals is used to solve the resulting dynamic equations.
Our numerical experiments for the present problem showed that the value of α at which the convergence is attained depends considerably on the value of the beam length.Particularly, for very short beams, the numerical convergence has been achieved using only too small values of α.This phenomenon can be easily justified from the results presented in Figure 1, wherein the variations of the regularized Dirac-delta function versus x are shown.From Figure 1, we see that the regularized Dirac-delta function has nonzero values at x < 0.2 for every value of α ≥ 0.0625.This implies that for the case of beams with L < 0.2, one should use very small values of α in the method to insure the convergence and accuracy of the solutions (α << 0.0625).
A way for overcoming the above-mentioned drawback is to express the governing differential equation of the beam in dimensionless form such that the domain of the problem is changed from 0 ≤ x ≤ L to 0 ≤ X ≤ 1.Then a similar procedure as that explained in Section 4 can be applied to discretize the resulting non-dimensional differential equation.By doing so, larger values of α can be used in the method, as we will show in the following (This can also be seen from the numerical results of the static analysis (see Figures 5-8)).By introducing the dimensionless variable X = x/L, the governing differential equation of motion of the beam can be written from Eq. ( 11) as Now, introducing the regularized Dirac-delta function (10) to Eq. ( 36) gives  12)); and (ii) the case where the dimensionless differential equation of the beam is considered (see Eq. ( 37)).As it can be seen from Figures 11 and 12, the results of case (i) are not very satisfactory in terms of both accuracy and convergence.The results of case (ii) are found to be more superior to those of case (i), particularly in terms of numerical convergence.Note that when case (ii) is considered, better accuracy is achieved using significantly larger vales of α, as compared with the case (i). Figure 12: Convergence and accuracy of the numerical results with decreasing α-value for the case where the dimensionless differential equation of motion of the beam is discretized using the DQM (n = 51).
On the other hand, from numerical results presented in Figure 12, one sees the rate of convergence of the method with α increases as the load speed increases.Therefore, larger values of α can be used in the method for the case of high speed of travelling load.Figure 13 shows the convergence of solutions with respect to the number of grid points (n).It can be seen from Figure 13 that the results of proposed method converge uniformly and agree well with the analytical solutions.It can also be clearly seen from this figure that the rate of convergence of the method with n increases as the load speed increases.This implies that smaller number of grid points can be used in the method for the case of high speed of travelling load.The numerical experiments are also carried out to obtain the dynamic magnification factors (DMFs) (maximum value of normalized deflection at the beam center) for various values of moving speed.The numerical results are given in Now, consider a moving point load on a beam with clamped-clamped boundary conditions.Although analytical solutions do not exist for this condition, accurate numerical solutions are available in the literature (Rieker et al., 1996).Therefore, the accuracy of the method can also be checked for this case.The dynamic responses at the beam center are evaluated for different travel speeds (v) and normalized by the static deflection w cs = fL 3 /(192EI) of a fully clamped beam under a point load f at mid-span.The numerical results are tabulated in Tables 2 and 3  element solutions obtained by Rieker et al. (1996).Note that these results are obtained using two different values of α.Excellent agreements are observed between the results of present analysis and those of Rieker et al. (1996).Besides, the convergence behavior of solutions is more uniform in the case α = 0.025 as compared with the case α = 0.02.However, the accuracy of the solutions in the case α = 0.02 is better than that in the case α = 0.025.From the numerical results presented in this section, it can be concluded that the proposed approach is very well suitable for handling the moving load problem on beam-like structures.

Vibration of a Rectangular Plate Due to a Moving Point Load
To demonstrate the applicability of the proposed method for the forced vibration analysis of rectangular plates carrying moving dynamic loads, application is made to a numerical example given by Eftekhari and Jafari (2012d).The parameters used in this numerical example are as follows: where v is the velocity of the moving load in the x-direction and t is the time.Other parameters are defined already in Section 5.As pointed out earlier, an analytical solution for the present problem is given in Fryba (1999) for plates with simply supported end conditions which will be used here to validate the proposed method.Figure 14 demonstrates the influence of the regularization parameter, α, on the accuracy of solutions and Figure 15 shows the convergence behavior of solutions with respect to the number of grid points.It can be seen from Figure 14 that the accuracy of solutions is improved significantly by decreasing the magnitude of α.On the other hand, by comparing these results with those of static analysis (see Figures 9 and 10), one sees that a smaller value of α should be used in the dynamic analysis of plates.For instance, as it can be seen from Figure 9, an acceptable accurate solutions for the static analysis of plates are achieved using α = 0.05.This value for the dynamic analysis of plates is α = 0.03 which is slightly smaller than that of the static analysis.The reason for this may be due to accumulation of the error in the time integration scheme, noting that the location of the applied load may affect the error of numerical results at each time step (see Figures 9 and 10 for more details).Now, consider a moving point load on a plate with fully clamped boundary conditions.The numerical results for this case are presented in Figure 16.The Ritz-DQM solutions of Eftekhari and Jafari (2012d) are also shown for comparison.It can be seen that the results of proposed method converge smoothly and agree well with those of the Ritz-DQM.

CONCLUSIONS
There are many problems in applied mechanics whose governing partial differential equations involve the Dirac-delta function.For instance, the behavior of structures under impulsive loading, the static and dynamic behavior of structures under impact loading, the thermoelastic dynamic behavior of structures under heat or pressure point sources can be mathematically described by means of the Dirac-delta function.In each of these problems, the accurate prediction of the structure response is very important and is a key factor in the design and analysis of such structures.In general, such type of problems can be solved by the help of analytical and/or numerical methods.Among the numerical methods, the point discretization techniques such as the collocation, finite difference, discrete singular convolution, and the differential quadrature methods have been successfully applied by many researchers to various problems encountered in applied mechanics.
The direct discretization of the Dirac-delta function using point discretization techniques like the DQM is not an easy task and special treatment is required.A way for overcoming this difficulty is to approximate the Dirac-delta function with simpler mathematical functions.Based on this idea, in this paper, a simple DQ procedure is presented for the numerical solution of moving load problem, where the location of the applied load is described by a time-dependent Dirac-delta function.To demonstrate its applicability, the proposed procedure is applied here to moving load problems of beams and rectangular plates.Numerical results prove that the proposed procedure is reliable and accurate.

Figure 1 :
Figure 1: Effect of α-value on accuracy and smoothness of the regularized Dirac-delta function.
26)where ij I are the elements of m × m identity matrix, the DQM weighting coefficients of the second-and fourth-order y-derivatives, respectively.Furthermore,

Figure 2 :
Figure 2: Simply-supported beam with a concentrated point load.

Figure 5 :Figure 6 :
Figure 5: Convergence and accuracy of the dimensionless deflection (w/β, β = fL 3 /EI) of a simply supported beam subjected to a concentrated point load for different locations of the applied load (α = 0.25)

Figure 7 :
Figure 7: Convergence and accuracy of the dimensionless deflection (w/β, β = fL 3 /EI) of a simply supported beam subjected to a concentrated point load for different values of α.

Figure 8 :
Figure 8: Variations of relative error in numerical results with α for different locations of the applied load (n = 45)

Figure 9 :
Figure 9: Convergence and accuracy of the dimensionless deflection (w/γ, γ = fa 2 /D) of a simply supported square plate subjected to a concentrated point load for different values of α.

Figure 10 :
Figure 10: Variations of relative error in numerical results with respect to α for different locations of the applied load (n = m = 45) Figures 11 and 12 show the effect of regularization parameter, α, on accuracy and convergence of the results for two different values of moving speed.The analytical solution values are also shown for comparison.Note that the numerical results are shown for two different cases: (i) the case where the dimensional differential equation of the beam is considered (see Eq. (12)); and (ii) the case where the dimensionless differential equation of the beam is considered (see Eq. (37)).As it can be seen from Figures11 and 12, the results of case (i) are not very satisfactory in terms of both accuracy and convergence.The results of case (ii) are found to be more superior to those of case (i), particularly in terms of numerical convergence.Note that when case (ii) is considered, better accuracy is achieved using significantly larger vales of α, as compared with the case (i).

Figure 11 :
Figure 11: Convergence and accuracy of the numerical results with decreasing α-value for the case where the dimensional differential equation of motion of the beam is discretized using the DQM (n = 51).

Figure 13 :
Figure 13: Convergence and accuracy of numerical results for normalized central deflection of a simply supported beam subjected to a moving load for different load speeds.

Figure 14 :
Figure 14: Effect of α-value on the accuracy of numerical results for central deflection of a simply supported square plate subjected to a moving load for different load speeds (n = m = 45).

Figure 15 :
Figure 15: Convergence and accuracy of numerical results for central deflection of a simply supported square plate subjected to a moving load for different load speeds (α = 0.03).

Figure 16 :
Figure 16: Convergence and accuracy of numerical results for central deflection of a clamped square plate subjected to a moving load for different load speeds (α = 0.03).

Table 1 .
These results are obtained using α = 0.018 and different values of n (number of grid points).The analytical solution values are also tabulated in this table for comparison purposes.It can be seen from Table1that the results generated by the proposed method converge rapidly and agree very closely with the analytical solutions.

Table 1 :
Convergence and accuracy of DMFs of a simply supported beam subjected to a moving concentrated load for different load speeds (α = 0.018).

Table 2 :
Convergence and accuracy of DMFs of a clamped beam subjected to a moving concentrated load for different load speeds (α = 0.025).

Table 3 :
Convergence and accuracy of DMFs of a clamped beam subjected to a moving concentrated load for different load speeds (α = 0.02).