A Differential Quadrature Procedure with Direct Projection of the Heaviside Function for Numerical Solution of Moving Load Problem

Owing to its particular characteristics, the direct discretization of the Dirac-delta function is not feasible when point discretization methods like the differential quadrature method (DQM) are applied. A way for overcoming this difficulty is to approximate (or regularize) the Dirac-delta function with simple mathematical 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. On the other hand, it is possible to combine the DQM with the integral quadrature method (IQM) to handle the Dirac-delta function. Alternatively, one may use another definition of the Dirac-delta function that the derivative of the Heaviside function, H(x), is the Dirac-delta function, δ(x), in the distribution sense, namely, dH(x)/dx = δ(x). This approach has been referred in the literature as the direct projection approach. It has been shown that although this approach yields highly oscillatory approximation of the Dirac-delta function, it can yield a non-oscillatory approximation of the solution. In this paper, we first present a modified direct projection approach that eliminates such difficulty (oscillatory approximation of the Dirac-delta function). We then demonstrate the applicability and reliability of the proposed method by applying it to some moving load problems of beams and rectangular plates.


INTRODUCTION
The differential quadrature method (DQM), which was firstly introduced by Bellman and his associates (1971,1972) in the early 1970s, is a powerful numerical method for the direct solution of partial differential equations that arise in various fields of engineering, mathematics, and physics (Bert and Malik, 1996a;1996b;Malik and Bert 1996;Tornabene et al., 2015a;Eftekhari 2016a).It is simple to use and also straightforward to implement.However, in spite of its many desirable features, 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 mainly caused by the fact that the Dirac-delta function cannot be directly discretized by the DQM.
A way for overcoming the above-mentioned difficulty is to approximate (or regularize) the Dirac-delta function with simple mathematical 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.Jung (2009), Jung and Don (2009), Jung et al. (2009), Wang et al. (2014), Eftekhari (2015a), and Tornabene et al. (2015b) have successfully applied this technique in conjunction with the point discretization methods to solve various partial differential equations involving the Dirac-delta function.However, this technique involves a regularization parameter that should be carefully adjusted before the problem being solved.This parameter has also been shown to be a problem-dependent parameter (Eftekhari, 2015a).
On the other hand, Eftekhari (2015bEftekhari ( , 2016b) ) proposed a combined application of the DQM and integral quadrature method (IQM) to handle the Dirac-delta function in vibration problem of beams and rectangular plates subjected to a moving (or static) point load.In this combined approach, the particular properties of the Dirac-delta function are taken into account by application of the IQM.The coupled DQM/IQM approach has been shown to be highly accurate and reliable.However, this approach requires the help of another method (i.e., the IQM) to treat the Dirac-delta function.In general, it is more desirable to solve such type of problems by the DQM itself, not by combining the DQM with other methods or techniques.
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(x), is the Dirac-delta function, δ(x), in the distribution sense, namely, dH(x)/dx = δ(x).The method was referred to as the direct projection method.It has been shown that although the direct projection of the Heaviside function gives 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.Recently, Wang et al. (2014) studied at length various ways of approximation of the Dirac-delta function.They considered the sine-Gordon equation with impulsive forcing and concluded that the domain decomposition method is the best among the methods considered in their study.
In this study, a differential quadrature procedure based on the direct projection of the Heaviside function is proposed for the numerical solution of moving load problem, wherein the motion of the point load is described by a time-dependent Dirac-delta function.Since the original direct projection approach shows some oscillations for approximation of the Dirac-delta function, we also introduce a modified direct projection approach that eliminates such difficulty.To demonstrate the applicability and reliability of the proposed approach, it is applied herein to some moving load problems of Latin American Journal of Solids and Structures 13 (2016) 1763-1781 beams and rectangular plates.Numerical results are presented and compared with analytical and numerical results 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.

DQM
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 function 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 , are given by (Quan and Chang, 1989) where The weighting coefficients of the higher-order derivatives (r > 1) can be computed from the following recurrence relationship (Bert and Malik, 1996a) 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 .

DIRECT PROJECTION APPROACH
As we discussed before in introduction, the direct discretization of the Dirac-delta function using point discretization methods like the DQM is not feasible.This is mainly due to particular characteristics of the Dirac-delta function.For instance, the one-dimensional Dirac-delta function has the following characteristics: 0 where f(x) may be any arbitrary function.As we mentioned in introduction, there has been a considerable effort to overcome this difficulty.In this paper, we have tried to solve this difficulty by the help of direct projection of the Heaviside function.The details of this approach will be given in the following sub-sections.
3.1 Approximation of the One-dimensional Dirac-delta Function

Conventional Approach
The direct projection approach is based on the fact that the derivative of the Heaviside function is the Dirac-delta function in the distribution sense, namely, where is the Heaviside function defined as (Jung, 2009) By the help of above definition, the Dirac-delta function can now be directly and simply discretized using the DQM.Let x 1 , x 2 , …, x n be a set of grid points in the x-direction.Application of the DQM to Eq. ( 9) gives Latin American Journal of Solids and Structures 13 (2016) 1763-1781 where [A] (1) is the DQM weighting coefficient matrix of the first-order derivative, and Figure 1: Approximation of the one-dimensional Dirac-delta function using conventional direct projection approach.
Figure 1 shows the results of conventional direct projection approach for approximation of the Dirac-delta function for different number of grid points (n).The coordinates of the grid points are calculated using ).As it can be seen, the results of this approach are high oscillatory.The amplitude of the oscillation becomes even more remarkable near the singular point (x = 0).This trend of solutions has also been reported earlier by Jung (2009) and Jung et al. (2009).

Proposed Approach
One of the main limitations of the conventional direct projection approach is that the case of arbitrary location of the singular point (x 0 ) cannot be accurately handled using this approach.Even if the location of the singular point being very close to one of the DQM grid points, this approach is found to show remarkable oscillation for approximation of the Dirac-delta function (see Figure 1).This limitation can be overcome if the weighted distribution of the source point is considered on the two nearest nodes around it (see Figure 2).By doing so, the following definition can be introduced for the Dirac-delta function where Now, introducing Eq. ( 9) to Eq. ( 14) gives Application of the DQM to Eq. ( 16) gives where } {δ is defined already in Eq. ( 12) and When above definition is used, Figure 3 illustrate the results.Comparing these results with those of Figure 1, one sees that the proposed method gives much better approximation for the Dirac-delta function as compared with the conventional direct projection approach.The amplitude of Latin American Journal of Solids and Structures 13 (2016) 1763-1781 oscillations is damped out considerably.However, there are some small amplitude oscillations near the singular point which are known to be due to the Gibbs phenomenon (for more details about this phenomenon see Myint-U (1980Myint-U ( , 2007))).

Approximation of the Two-Dimensional Dirac-Delta Function
The approaches presented in previous section (section 3.1) can be easily extended to discretize the two-dimensional Dirac-delta function.The details are given in the following sub-sections.

Conventional Approach
Using Eq. ( 9), the two-dimensional Dirac-delta function can be expressed as where is defined in Eq. ( 10) and 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.Application of the DQM to Eq. ( 21) gives where ] [A are the first-order DQM weighting coefficient matrices associated with the x-and y-derivatives, respectively.Furthermore, .Using Eqs. ( 9) and ( 14), the two-dimensional Dirac-delta function can be expressed as are defined in Eqs. ( 15) and ( 17).Furthermore, , 0 Now, application of the DQM to Eq. ( 28) gives where Latin American Journal of Solids and Structures 13 (2016) 1763-1781 In Figures 4 and 5 the results of proposed approach are compared with those of conventional approach for approximation of the two-dimensional Dirac-delta function.These results are obtained using n = m = 51.The results of the conventional approach are observed to be highly oscillatory.The results of the proposed approach, however, do not show such an oscillatory behavior.

FORMULATION FOR FORCED VIBRATION OF BEAMS CARRYING MOVING LOADS
Consider an Euler-Bernoulli beam with length L, mass per unit length ρA, and bending stiffness EI subjected to a point 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 transverse deflection of the beam, and x f (t) is the position coordinate of the moving point load.
Consider n grid points with coordinates x 1 , x 2 , …, x n in the x-direction.In view of Eq. ( 16), Eq.
(38) can be rewritten as ( Satisfying Eq. ( 39) at any grid point x = x i (i = 1, 2, …, n) and substituting the quadrature rule into results gives where [I] is an identity matrix of order n × n, [A] (r) (r = 1, 4) are the DQM weighting coefficient matrix of the rth-order derivative, and the matrices } { L H and } { R H are defined already in Eqs.
(19) and ( 20).It should be noted that since the position of the moving point load changes from time to time, the variables x r and x l are time-dependent variables.For the same reason, the vectors are time-dependent vectors.As a result, the load vector )} ( { t F is also a timedependent vector. After applying the boundary conditions of the beam, Eq. ( 40) can be solved for the unknown gridpoint 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) and also in Eftekhari (2015c).

LOADS
Consider a rectangular thin plate with length a, width b, mass per unit area ρh, and bending rigidity D under a moving point load f.The governing differential equation of motion of the rectangular plate is given by (Rao, 2007) )), ( where w(x, y, t) is the transverse deflection of the plate, x f (t) and y f (t) are position coordinates of the moving point load in the x-and y-direction, respectively, and v x and v y are the corresponding moving speeds.
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.In view of Eq. ( 28), Eq. ( 45) can be rewritten as Satisfying Eq. ( 46) at any grid point x = x i and substituting the quadrature rule1), into results gives where [I] is an identity matrix of order n × n, and [A] (r) (r = 2, 4) are the DQM weighting coefficient matrices of the rth-order x-derivative.Furthermore, where the matrices are defined in Eqs. ( 34) and (35).
Satisfying Eq. ( 47) at any grid point y = y i and substituting the quadrature rule into results gives where ij I are the elements of m × m identity matrix, and ) (r ij A (r = 2, 4) are the DQM weighting coefficients of the rth-order y-derivative.Furthermore,

H
are defined in Eqs. ( 36) and ( 37).It should be noted that since the location of the moving load changes from time to time, the variables x r , x l , y r , and y l also After applying the plate boundary conditions, one can solve the resulting system of ordinary differential equations for unknowns using various time integration schemes.The details of implementation of the plate boundary conditions can be found in Malik and Bert (1996) and also in Eftekhari (2015c).In this study, the Newmark method (Bathe and Wilson, 1976) is employed to solve Eqs. ( 40) and (50).

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 beams and rectangular plates subjected to point loads.These numerical examples are presented herein to better verify the accuracy and convergence of the method.The same problems are considered in the third and fourth examples, but the location of the point load is assumed to vary with time.

Deflection of a Simply Supported Beam Due to a Concentrated Point Load
Consider the bending problem of a simply supported beam subjected to a concentrated load f acting at a point x = x f .Figure 6  using the conventional approach are also shown for comparison purposes.It can be seen that the results of proposed approach converge uniformly to their final values.The results of the conventional approach, however, are somewhat erratic and their behavior is somehow oscillatory.In Figure 7 the results of proposed approach are compared with those of the DQMR (DQM with regularized Dirac-delta function).The DQMR results are shown for three different values of α.It is noted that in DQMR, the Dirac-delta function is approximated as (Eftekhari, 2015a) where α is the regularization parameter.It can be seen from Figure 7 that the results of the DQMR is highly dependent on the values of α.For large values of α, the results of the DQMR are not very accurate, and for small values of this parameter the DQMR solutions show some oscillation for small number of grid points.The results of the proposed approach are found to be better than those of the DQMR in terms of both accuracy and convergence.Consider the bending problem of a simply supported rectangular plate subjected to a concentrated point load f acting at a point (x, y) = (x f , y f ). Figure 8 shows the convergence and accuracy of the dimensionless central deflection of the plate for different locations of the applied load.It can be seen that the results by present method show better convergence trend than the conventional approach.On the other hand, by comparing these results with those of Figure 6, one sees that the behavior of solutions in this case (plate problem) is very similar to that of the beam problem.The only difference is that the order of error in solutions of the plate problem is larger than that of the beam problem.
In Figure 9 the results are compared with those of the DQMR.Again, one sees that the results of proposed approach are better than those of the DQMR in terms of both accuracy and convergence.The results are also compared with those obtained using the conventional approach.It can be seen from Figures 10 and 11 that the results generated by the proposed method converge rapidly to their final values and agree well with analytical solutions.The results of the conventional approach, however, are highly oscillatory and erratic.Note that the error in solutions of the conventional approach never decays even if a large number of grid points are used.This implies that the conventional approach is not suitable for studying the moving load-type problems.In Table 1, the numerical results are given for dynamic magnification factors (DMFs) (maximum value of normalized deflection at the beam center) of a simply supported beam for different values of moving speed.The analytical solution values are also shown for comparison.It can be seen from Table 1 that the results of proposed approach converge uniformly to those of exact solution as the number of grid points increases.Besides, they agree very closely with analytical solution values.In Table 2, the results of proposed approach are compared with those of other DQM procedures.It can be seen from Table 2 that the results generated by the DQMR are somewhat oscillatory.Note that the results of this approach with n = 25 are erroneous.By increasing the number of grid points, however, the DQMR results approach to the exact solution results.Compared with the DQMR approach, the proposed approach provides better convergence and, in most cases, leads to better accuracy.The DQM/IQM approach is found to be more efficient than other two approaches.But, this approach requires the help of another method (i.e., the IQM) to treat the Dirac-delta function, and this is not very desirable.Now consider a clamped beam under a moving point load.The numerical results for this case are given in Table 3 together with the finite element solution results of Rieker et al. (1996).Note that, in this case, the central deflections of the beam are normalized by the static deflection w cs = fL 3 /(192EI).It can be seen that the present results converge uniformly to those of Rieker et al. (1996) as the number of grid points increases.In Table 4, the results of proposed method are compared with those of DQMR and DQM/IQM approaches.As it can be seen, the differences between the solutions of various DQM procedures are negligible.But, when n = 31, the solutions of the Latin American Journal of Solids and Structures 13 (2016) 1763-1781 DQM/IQM approach are found to have closer agreement with the finite element solution results of Rieker et al. (1996).It can be seen from Figure 12 that the computed solutions converge monotonically to the exact solutions as the number of grid points increases.Convergence is also more rapid when the moving speed is high.Figure 13 shows again the better performance of the proposed approach over the DQMR approach.Note that the DQMR approach may produce erroneous results with small number of grid points (n = m = 25) when the value of the regularization parameter (α) is too small (α = 0.035).Although this can be avoided by increasing the number of grid points (see the results correspond to n = m = 41), the regularization parameter should be carefully adjusted before the problem being solved.In the proposed approach, there is no need to adjust any parameter before solving the problem.Therefore, the proposed approach is more straightforward to implement than the DQMR approach.

CONCLUSIONS
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 use an alternative definition of the Dirac-delta function that the derivative of the Heaviside function in the distribution sense is the Dirac-delta function.This approach has been referred in the literature as the direct projection approach.The conventional direct projection approach, however, shows some oscillations for approximation of the Dirac-delta function.To overcome this limitation, this paper introduces a modified direct projection approach that provides better approximation for the Dirac-delta function while eliminates such difficulty.The applicability and reliability of the proposed method are demonstrated herein through solution of some moving load problems of beams and rectangular plates.The proposed approach is shown to be superior over the DQMR approach in terms of both accuracy and convergence.

Figure 2 :
Figure 2: Representation of a point source by two statically equivalent point sources at its adjacent grid points.

Figure 3 :
Figure 3: Approximation of the one-dimensional Dirac-delta function using proposed direct projection approach.

Figure 4 :
Figure 4: Approximation of the two-dimensional Dirac-delta function using conventional direct projection approach (n = m = 51).

Figure 5 :
Figure 5: Approximation of the two-dimensional Dirac-delta function using proposed direct projection approach (n = m = 51).
presents the variations of the percent error in numerical solutions (den for different values of x f .The numerical results obtained Latin American Journal of Solidsand Structures 13 (2016) 1763-1781

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

Figure 7 :
Figure 7: Comparison of the dimensionless deflection (w/β, β = fL 3 /EI) of a simply supported beam subjected to a point load for different locations of the applied load.
of a Simply Supported Rectangular Plate Due to a Concentrated Point Load

Figure 8 :
Figure 8: Convergence and accuracy of the dimensionless central deflection (w/γ, γ = fa 2 /D) of a simply supported square plate subjected to a point load for different locations of the applied load.

Figure 9 :
Figure 9: Comparison of the dimensionless deflection (w/γ, γ = fa 2 /D) of a simply supported square plate subjected to a point load for different locations of the applied load.
of a Beam Due to a Moving Point Load Consider the vibration problem of a simply supported beam subjected to a moving point load.The parameters of the problem are assumed to be (Eftekhari, 2015a): ρ = 10686.9kg/m 3 , b = 0.635 cm, h = 0.635 cm, E = 2.068 × 10 11 Pa, L = 10.16 cm, f = 4.45 N where b and h are the width and thickness of the beam.The dynamic central deflections of the beam, w cd , are evaluated for different values of moving speed (v) and normalized by the static deflection w cs = fL 3 /(48EI).Figures 10 and 11 present the convergence of solutions for different load speeds.

Figure 10 :
Figure 10: Convergence and accuracy of numerical results for normalized central deflection of a simply supported beam subjected to a moving point load (v = 31.2m/s).

Figure 11 :
Figure 11: Convergence and accuracy of numerical results for normalized central deflection of a simply supported beam subjected to a moving point load (v = 62.4 m/s).

6. 4
Figure12illustrates the convergence of solutions and Figure13presents a comparison of the results with the DQMR results.

Figure 12 :
Figure 12: Convergence and accuracy of numerical results for normalized central deflection of a simply supported rectangular plate subjected to a moving point load.

Figure 13 :
Figure 13: Comparison of numerical results for normalized central deflection of a simply supported rectangular plate subjected to a moving point load (v = 0.45 m/s).

Table 1 :
Convergence and accuracy of DMFs of a simply supported beam subjected to a moving point load for different load speeds.

Table 2 :
Comparison of DMFs of a simply supported beam subjected to a moving point load for different load speeds.

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

Table 4 :
Comparison of DMFs of a clamped beam subjected to a moving point load for different load speeds.