- Email: [email protected]
Application of generalized differential quadrature to vibration analysis
Journal of Sound and Vibration (1995) 181(2), 279–293
APPLICATION OF GENERALIZED DIFFERENTIAL QUADRATURE TO VIBRATION ANALYSIS H. D, M. K. L R. M. L School of Mechanical and Production Engineering, Nanyang Technological University, Singapore 2263 (Received 12 February 1993, and in final form 4 January 1994) This paper presents a first endeavour to exploit the generalized differential quadrature method as an accurate, efficient and convenient numerical technique for vibration analysis. The generalized differential quadrature (GDQ) overcomes some major drawbacks existing in the method of differential quadrature (DQ). In GDQ, the weighting coefficients for derivative approximations are given by a simple algebraic expression and a recurrence relationship, together with arbitrary choice of grid points. GDQ is then applied to solve some free vibration problems. The accuracy, efficiency and convenience of GDQ are demonstrated throughout the numerical examples.
1. INTRODUCTION
Often, numerical approximation methods have to be sought to solve vibration problems due to the complexity of the problems. Classical techniques such as finite element and finite difference methods are well developed and well known. These methods can provide very accurate results by using a large number of grid points; thus they are computationally expensive. In a large number of cases, only a limited number of frequencies and mode shapes or dynamic responses at only a limited number of points is needed to be found. However, for acceptable accuracy, conventional numerical methods still require the use of a large number of grid points to obtain accurate solutions. Consequently, the requirement for computer capacity is often unncessarily large in such cases. An alternative method called the method of differential quadrature may be used for such purposes. The method of differential quadrature was first introduced by Bellman et al. [1, 2] for numerically solving partial differential equations using considerably fewer grid points. This method is based on the idea that the derivative of a function with respect to a coordinate direction can be expressed as a weighted linear sum of all the function values at all mesh points along that direction. Applications of this method to various problems showed that it has potential as an attractive approximation technique [3–9]. However, there exist some major drawbacks to be overcome in the originally proposed differential quadrature. As is known, the most important part of differential quadrature is to determine the weighting coefficients for discretization of any order partial derivative. There were two methods in use which were proposed by Bellman et al. [2] to obtain the weighting coefficients. One method is to solve a set of algebraic equations which satisfy exactly the linear constrained relation for all polynomials of degree less than or equal to N − 1. This set of equations has a unique solution because the matrix elements are composed of a Vandermonde matrix. Unfortunately, when N is large the Vandermonde matrix is ill-conditioned and the inversion of this matrix becomes difficult. Moreover, a set of N × N linear algebraic equations has to be solved for each order derivative even if 279 0022–460X/95/120279 + 15 $08.00/0
7 1995 Academic Press Limited
. .
280
the equations are solvable. The other method is to compute the weighting coefficients by a simple algebraic formula, but with the coordinates of grid points chosen as the roots of an Nth order shifted Legendre polynomial. This means that if N is specified, the distribution of grid points is fixed even for different physical problems or different boundary conditions. This also creates a major drawback and restricts the application of differential quadrature, since some practical problems may need more grid points near the boundary, while some others may not. Obviously, both of the methods originally proposed for determining the weighting coefficients in the method of differential quadrature have some major drawbacks. In order to overcome these drawbacks, Shu presented a generalized differential quadrature in connection with solving partial differential equations in fluid mechanics [10]. Preliminary results in solving some fluid mechanics problems have shown that the method is very efficient and convenient [11, 12]. The objective of this paper is to explore the great potential of the generalized differential quadrature method as an accurate, efficient and convenient numerical method for vibration analysis. The drawbacks existing in the original differential quadrature (DQ) method are first summarized and quantified in relation to its application to vibration analysis. Then, the generalized differential quadrature method (GDQ) is introduced. The generalized differential quadrature has overcome the possible singularity problem of the original differential quadrature in obtaining the weighting coefficients. Furthermore, in GDQ, the weighting coefficients for the first derivative and for the higher order derivatives are given by a simple algebraic expression and a recurrence relationship respectively, and without any restriction on arbitrary choice of grid points. GDQ is finally applied to solve for natural frequencies of some structures under various boundary conditions. The ease of use and the accuracy of the GDQ are demonstrated through the numerical examples. 2. DIFFERENTIAL QUADRATURE
The method of differential quadrature is based on the idea that the partial derivative of a function with respect to a space variable at a given discrete point can be expressed as a weighted linear sum of the function values at all discrete points in the domain of that variable. Let us take the first derivative of a one dimensional function u(x, t) as an example: the higher order partial derivatives will have essentially the same formation. A differential quadrature approximation of the first derivative of the function u(x, t) at the ith discrete point on a grid is N
ux (xi , t) = s cij(1) u(xj , t)
for i = 1, 2, . . . , N,
(1)
j=1
where ux (xi , t) is the first derivative of u(x, t) with respect to x at xi , and N is the number of discrete grids. cij(1) are the weighting coefficients for the first derivative approximation, which have to be pre-determined. There are two approaches used in the method of differential quadrature to determine the weighting coefficients cij(1) [2]. The first one is to let equation (1) be exact for all polynomials of degree less than or equal to (N − 1), g(x) = x k, k = 0, 1, . . . , N − 1. Then, one has a set of N × N linear algebraic equations: N
s cij xjk = kxik − 1 j=1
for k = 0, 1, . . . , N − 1
and
i = 1, 2, . . . , N.
(2)
281
This set of equations has a unique solution, since its matrix is of Vandermonde form. Unfortunately, it has been found that this set of equations becomes ill-conditioned and is difficult to solve when N is large. In order to quantify this singularity, weighting coefficients have been calculated for equal space grids based on equation (2) for various numbers of grid points N. From the computation, it is found that the maximum number of grid points is N = 22. Once the grid number if greater than 22, the set of linear algebraic equations become singular and cannot be solved at all. The computed weighting coefficients are also compared with the accurate coefficients obtained from the generalized differential quadrature method to be introduced later; it is also found that the results from equations (2) have some errors when the grid number is more than 17. Therefore, the maximum number of grid points is practically 17 for equal space grids if this method is used for determining the weighting coefficients. In addition, one has to solve a set of N × N linear equations for each order of the derivatives in the governing equations even when the equations are solvable (N E 17). The other approach to determine the weighting coefficients is similar to the first one with the exception that a different set of test functions g(x) is chosen for satisfying equation (1) exactly: g(x) = LN (x)/(x − xj )L'N (xj ) for j = 1, 2, . . . , N. (3) Here N is the number of grid points, LN (x) is the Nth order Legendre polynomial and L'N (x) is the first derivative of LN (x). By choosing xj to be the roots of the shifted Legendre polynomial and substituting equation (3) into equation (1), Bellman et al. [2] obtained a direct simple algebraic expression for the weighting coefficients cij(1) = L'N (xi )/(xi − xj )L'N (xj ) (1) ii
c = (1 − 2xi )/2xi (xi − 1),
for i $ j; i, j = 1, 2, . . . , N.
(4a, b)
In this approach, the weighting coefficients are easy to obtain without solving algebraic equations or having a singularity problem. However, it is obvious that once the number of grid points N is specified, the roots of the shifted Legendre polynomial are given; thus the distribution of the grid points is fixed no matter what problems are considered. This puts a major restriction on this method to problems in structural dynamic analysis, since all sorts of boundary conditions could appear and different mesh grids may be needed for different boundary conditions and different structure geometry. From the above discussion, some deficiencies have been identified in the originally proposed differential quadrature. These deficiences put some restrictions on the application of the DQ method to structural dynamics problems. This is probably one of the reasons that the method of differential quadrature is not widely used in structural analysis. In order to overcome the deficiencies, a generalized differential quadrature, which was recently proposed by Shu [10] for solving partial differential equations in fluid mechanics, will be introduced and applied to solve some problems in vibration analysis. The ease of use and the accuracy of the GDQ will be demonstrated through numerical examples. 3. GENERALIZED DIFFERENTIAL QUADRATURE
As pointed out above, two approaches have been proposed by previous researchers for the differential quadrature method for determining the weighting coefficients. Both of them have some drawbacks. The first restricts one to a small number of grids points and the need to solve sets of linear equations. The second limits the distribution of the grid points, which is critical in structural dynamic analysis. To remedy these deficiencies, what one
. .
282
wants is to find a method to determine the weighting coefficients without any limitation on the choice of grid meshes and which yet provides a simple algebraic expression to determine the weighting coefficients. Such a method was recently proposed by Shu [10] in relation to solving some partial differential equations in fluid mechanics. In order to find a simple algebraic expression for calculating the weighting coefficients without restricting the choice of grid meshes, Shu chose Lagrange interpolated polynomials as the set of tests functions g(x) instead of using power polynomials or Legendre polynomials: gi (x) = M(x)/(x − xi )M (1)(xi ) for i = 1, 2, . . . , N. (5) Here N
M(x) = t (x − xj ),
(6)
j=1
M (1)(x) is the first derivative of M(x), N
M (1)(xi ) = t
(xi − xj ),
(7)
j = 1, j $ i
and N is the number of grid points. When equation (1) is satisfied for all the test functions gi (x), a simple algebraic expression can be obtained to determine the weighting coefficients: cij(1) = M (1)(xi )/(xi − xj ) · M (1)(xj ) (1) ii
(2)
(1)
c = M (xi )/2M (xi ),
for i $ j;
i, j = 1, 2, . . . , N.
(8a, b)
Equations (8) provide simple expressions for computing cij(1) without any restriction in choice of the co-ordinates of the grid points xi . It is obvious that once the grid points xi are given, M (1)(x) is very easily obtained from equation (7). Hence, cij(1) can be easily calculated for i $ j. The calculation of cii(1) is based on the calculation of the second derivative of M(x) which is more difficult to obtain. Instead of using equation (8b), a more convenient relationship can be obtained and used for calculating cii(1) . It can be shown by using a Taylor series expansion that the following relationship exists for cij(1) : N
s cij(1) = 0
for i = 1, 2, . . . , N.
(9)
j=1
Thus, from equation (9), the coefficient cii(1) can be calculated from cij(1) (i $ j): that is, N
cii(1) = − s
cij(1)
for i = 1, 2, . . . , N.
(10)
j = 1,j $ i
For the approximation of higher order derivatives, the generalized differential quadrature gives a form similar to equation (1), N
ux(m) (xi , t) = s cij(m) u(xj , t)
for i = 1, 2, . . . , N,
m = 2, 3, . . . , N − 1,
(11)
j=1
where ux(m) (xi , t) is the mth derivative of u(x, t) with respect to x at xi and N is the number of discrete grid points. cij(m) are the weighting coefficients for the mth derivative approximation.
283
The weighting coefficients for second and higher order derivatives can be similarly obtained. Amazingly, a recurrence relationship has been found for the mth order weighting coefficients cij(m) :
0
cij(m) = m cii(m − 1)cij(1) −
cij(m − 1) xi − xj
1
for i $ j,
m = 2, 3, . . . , N − 1,
i, j = 1, 2, . . . , N.
(12) Thus the cij(m) can be derived from the (m − 1)th order weighting coefficients cij(m − 1) . The cii(m) can be obtained from a relationship similar to equation (10): N
cii(m) = − s
cij(m)
for i = 1, 2, . . . , N.
(13)
j = 1,j $ i
Thus equations (12) and (13) together with equation (8a) and equation (10) give a convenient and general form for determining the weighting coefficients for the derivatives of orders one through N − 1. There are no restrictions on the co-ordinates of the chosen grid points. There is no need to solve for the weighting coefficients from a set of algebraic equations which could be ill-conditioned when the number of grids is large. Furthermore, this set of expressions for the determination of the weighting coefficients is so compact and simple that it is very easy to be implemented in formulating and programming because of the recurrence feature. All these features are very convenient for solving practical problems in structural dynamic analysis. Application of this generalized differential quadrature method to time dependent dynamic problems will result in a set of ordinary differential equations with the time dependent function values at the grid points as unknowns. This leads to an eigenvalue problem for free vibration analysis or to a set of time dependent ordinary differential equations for forced vibration problems. The natural frequencies and mode shapes or time history response can then be obtained. Finally, once the function values at all grid points have been obtained, it is very easy to determine the function values in the overall domain in terms of a polynomial approximation: that is, N
u(x) = s u(xi )gi (x),
(14)
i=1
where gi (x) are the Lagrange interpolated polynomials of equation (5). 4. APPLICATION OF GDQ TO VIBRATION ANALYSIS
The method of generalized differential quadrature is used to analyze the free vibration problems of beams and plates in this section. The first case is calculating the natural frequencies of slender beams. The second and third cases are calculating the natural frequencies of a circular and a rectangular plate, respectively. The formulating and programming procedures are shown to be very straightforward and simple. The boundary conditions are easy to be implemented. 4.1. The governing equation of a Bernoulli–Euler beam in bending is
$
%
12 1 2w 1 2w = 0, 2 EI 2 + rA 1x 1t 2 1x
0 Q x Q L,
(15)
where EI is the beam’s flexural rigidity, rA is the mass per unit length and L is the length of the beam.
. .
284
For free vibration, the transverse displacement w is assumed as: w(x, t) = W(x) eivt.
(16)
Substituting expression (16) into equation (15) and normalizing the equation yields 1 2s(X) 1 2W 1s(X) 1 3W 1 4W = v¯ 2W, 2 2 +2 3 + s(X) 1X 1X 1X 1X 1X 4
(17)
where s(X) = EI/EI0 , v¯ 2 = (rAL 4/EI0 )v 2 and X = x/L. Given the number of grid points N (see Figure 1), and applying the generalized differential quadrature approximation (11) to equation (17) at each discrete point on the grid, one obtains 1 2s(Xi ) 1X 2
0
1
N
1s(Xi ) 1X
s cij(2) Wj + 2 j=1
0
1
N
0
1
N
s cij(3) Wj + s(Xi ) s cij(4) Wj = v¯ 2Wi , j=1
j=1
(18)
for i = 1, 2, . . . , N. Equations (18) give an eigenvalue problem. Solving this eigenvalue problem together with appropriate boundary conditions, one can obtain the natural frequencies of the structure. Consider a beam with left end fixed and right end free. The boundary conditions then are W = 1W/1X = 0 at X = 0, (19a) EI 1 2W/1x 2 = 0,
(1/1x)(EI 1 2W/1x 2) = 0
at x = L,
or s(X) 1 2W/1X 2 = 0,
{1s(X)/1X} 1 2W/1X 2 + s(X) 1 3W/1X 3 = 0
at X = 1.
(19b)
Applying the differential quadrature approximation (11) to the boundary conditions (19) yields N
W1 = 0,
s cij(1) Wj = 0,
(20a)
j=1
0
N
1
(2) s(XN ) s cNj Wj = 0, j=1
0
N
1
0
N
1
(2) (3) {1s(XN )/1X} s cNj Wj + s(XN ) s cNj Wj = 0. j=1
j=1
(20b)
In order to introduce the boundary conditions into equations (18), first equations (20) are used to solve for W1 , W2 , WN − 1 and WN in terms of the variables W3 , W4 , . . . , WN − 2 . The expressions for W1 , W2 , WN − 1 and WN in terms of the variables W3 , W4 , . . . , WN − 2 are then substituted into equations (18) to eliminate the variables W1 , W2 , WN − 1 and WN , and only the discretized equations at the points i = 3, 4, . . . , N − 2 are to be used in equations (18). Finally, by solving the remaining eigenvalue problem, the natural frequencies of the structure can be obtained. For a beam with its left end fixed and its right end simply supported, the boundary conditions are W = 1W/1X = 0
at X = 0,
W = s(X) 1 2W/1X 2 = 0
Figure 1. The grid points of a slender beam.
at X = 1.
(21a, b)
285
T 1 Natural frequencies of a uniform beam (clamped–free)
Exact [15] N = 8 [8] N = 8, uniform grid N = 10 uniform grid N = 12, uniform grid N = 14, uniform grid N = 12, non-uniform grid
v¯ 1
% error
v¯ 2
3·5160 3·5226 3·5224
— 0·19 0·18
22·035 — 20·837
3·5161
0·003
3·5160
% error
v¯ 3
% error
v¯ 4
% error
— — 5·4
61·697 — —
— — —
120·90 — —
— — —
22·213
0·81
—
—
—
—
0·0
22·025
0·045
61·298
0·65
131·93
9·1
3·5160
0·0
22·035
0·0
61·747
0·081
116·75
3·4
3·5160
0·0
22·033
0·009
61·604
0·15
125·63
3·9
Applying the differential quadrature approximation (11) to the boundary conditions (21) gives N
s cij(1)Wj = 0,
W1 = 0,
(22a)
j=1
WN = 0,
0
N
1
(2) s(XN ) s cNj Wj = 0. j=1
(22b)
Again, equations (22) can be used to eliminate the variables W1 , W2 , WN − 1 and WN in equations (18), as discussed above. The natural frequencies are then solved from the remaining eigenvalue problem. Calculations have been done for both a uniform beam s(X) = 1 and a non-uniform beam under two sets of different boundary conditions, namely clamped–free (C–F), and clamped – simply supported (C–S). The natural frequencies for the first four modes of a uniform beam are presented in Table 1 for clamped–free support and in Table 2 for clamped – simply supported support together with exact solutions. The first natural frequency obtained by Bert et al. [8] using the original differential quadrature method is also presented in Table 1. The present result for N = 8 is found to be consistent with it. T 2 Natural frequencies of a uniform beam (clamped – simply supported) v¯ 1 Exact [15] N = 8, uniform grid N = 10, uniform grid N = 12, uniform grid N = 14, uniform grid N = 12, non-uniform grid
% error
v¯ 2
% error
v¯ 3
% error
v¯ 4
% error
15·418 15·637
— 1·4
49·965 54·430
— 8·9
104·25 —
— —
178·27 —
— —
15·410
0·052
49·636
0·66
—
—
—
—
15·419
0·006
49·978
0·026
99·128
4·9
153·68
13·8
15·418
0·0
49·965
0·0
105·17
0·88
187·03
4·9
15·418
0·0
49·966
0·002
103·80
0·43
174·21
2·3
286
. .
Figure 2. The mode shapes of a clamped–free beam.
Results were obtained by using various numbers of grid points. From the tables, the convergence of this method is seen to be very good. Accurate results can be achieved by using very few grid points. Accurate natural frequencies for the first two modes can be obtained by using ten uniformly spaced grid points, while 14 uniformly spaced grid points are sufficient to give accurate frequencies for the first four modes. It is especially interesting to note that the accuracy can be improved by using non-uniform meshes with finer grids closer to the boundaries. The non-uniformly spaced grid points for N = 12 are chosen as X = 0·0, 0·05, 0·1, 0·15, 0·25, 0·4, 0·6, 0·75, 0·85, 0·9, 0·95, 1·0. With only 12 points, one can obtain very good results for the first four modes. The computing time is tiny since one only has to solve an eigenvalue problem of an 8 × 8 matrix after imposing the boundary conditions. To demonstrate better the accuracy of the solutions, the mode shapes of the first four modes were also obtained for both sets of boundary conditions by using 12 non-uniformly spaced grid points. The mode shapes are presented in Figures 2 and 3.
Figure 3. The mode shapes of a clamped – simply supported beam.
287
T 3 Natural frequencies of a non-uniform beam (N = 12 non-uniformly spaced grids) s(X) = 1 + 0·5X (C–F) s(X) = 1 − X (C–F) s(X) = 1 + 0·5X (C–S) s(X) = 1 − X (C–S)
v¯ 1
v¯ 2
v¯ 3
v¯ 4
3·6733 3·2746 16·928 10·739
24·079 16·871 55·409 31·693
67·978 43·700 115·48 65·311
136·17 82·761 191·20 110·92
In Table 3 are presented the natural frequencies of the first four modes for two non-uniform beams with different stiffness distributions, s(X) = 1 + 0·5X and s(X) = 1 − X, respectively. Again, two sets of different boundary conditions (C–F and C–S) were considered for each beam. Twelve non-uniformly spaced grid points were chosen in the calculation. The varying cross section stiffness s(X) of the beams can be very easily implemented from equations (18). The programming is also easy and straightforward for such more complex problems. The computing effort is still small since one has to solve an eigenvalue problems of a matrix of dimension 8 × 8 only. 4.2. The governing equation for a thin circular plate of isotropic material and undergoing axisymmetric motion is 1 4w 2 1 3w 1 1 2w 1 1w rh 1 2w + + − + = 0. 1r 4 D 1t 2 r 1r 3 r 2 1r 2 r 3 1r
(23)
Substituting equation (16) into equation (23) and normalizing it yields 1 4W 2 1 3W 1 1 2W 1 1W + − + − v¯ 2W = 0, 1R 4 R 1R 3 R 2 1R 2 R 3 1R
(24)
where R = r/a, v¯ 2 = rha 4v 2/D, a is the outside radius of the plate, h the thickness and D the flexural rigidity. The regularity condition at the centre of the plate is 1W/1R = 0
at R = 0.
(25)
The regularity condition is necessary to assure that the plate slope is zero at the origin to avoid a singularity at this location. The boundary conditions for a clamped outside edge are W = 0, 1W/1R = 0 at R = 1. (26) Applying the generalized differential quadrature approximation (11) to equation (24) at each discrete point on the grid (see Figure 4) yields N
s cij4 Wj + j=1
2 N (3) 1 N 1 N s cij Wj − 2 s cij(2) Wj + 3 s cij(1)Wj = v¯ 2Wi Ri j = 1 R1 j = 1 R1 j = 1
for i = 1, 2, . . . , N, (27)
and applying the generalized differential quadrature approximation (11) to the boundary conditions (25) and (26) yields N
s c1j(1) Wj = 0 j=1
N
and
WN = 0,
(1) s cNj Wj = 0.
(28a, b)
j=1
Again, using equations (28) to eliminate the variables most close to the boundaries, W1 , WN − 1 and WN , in equations (27) as discussed in the previous example, one can solve the remaining eigenvalue problem to obtain the natural frequencies.
. .
288
Figure 4. The grid points of a circular plate.
The natural frequencies obtained for the first four modes are given in Table 4, together with exact solutions available for comparison. The corresponding mode shapes obtained by using 12 non-uniformly spaced grid points are plotted in Figure 5. The first natural frequency obtained by Bert et al. [8] using the original differential quadrature method is also presented in Table 4. The present result for N = 8 is found to be consistent with it. As observed from Table 4, the first four natural frequencies obtained by using 12 non-uniformly spaced grid points are very accurate. The computing time is very small for solving the eigenvalues of the resulting 8 × 8 matrix. 4.3. The governing equation for a thin rectangular plate of isotropic material is given as: D 1 4w/1x 4 + 2D 1 4w/1x 2 1y 2 + D 1 4w/1y 4 + rh 1 2w/1t 2 = 0,
(29)
where r is the density of the material, h the thickness of the plate and D the flexural rigidity of the plate. The transverse displacement w for free vibration is taken as w(x, y, t) = W(x, y) eivt.
(30)
Substituting equation (30) into equation (29), one obtains the normalized equation 1 4W/1X 4 + 2b 2 1 4W/1X 2 1Y 2 + b 4 1 4W/1Y 4 − v¯ 2W = 0,
(31)
T 4 Natural frequencies of a circular plate
Exact [15] N = 8 [8] N = 8, uniform grid N = 12, non-uniform grid
v¯ 1
% error
v¯ 2
% error
v¯ 3
% error
v¯ 4
% error
10·22 10·20 10·20
— 0·19 0·19
39·77 — 44·68
— — 12
89·10 — 80·24
— — 9·9
158·2 — —
— — —
10·22
0·0
39·77
89·22
0·13
156·6
1·0
0·0
289
Figure 5. The mode shapes of a circular plate with outside edge clamped.
where b = a/b, v¯ 2 = rha 4v 2/D, X = x/a, Y = y/b, a is the length of the plate and b its width. Applying the generalized differential quadrature approximation (11) to equation (31) at each discrete point on the grid (see Figure 6) gives Nx
Ny
Nx
Ny
k=1
m=1
k=1
k=1
(2) s cik(4) Wkj + 2b 2 s cjm s cik(2) Wkm + b 4 s cjk(4) Wik = v¯ 2Wij ,
i = 1, 2, . . . , Nx ,
j = 1, 2, . . . , Ny ,
(32)
where Nx and Ny are the number of grid points along the X and Y directions, respectively (see Figure 6). The boundary conditions for a plate clamped on all four edges are: W(X, 0) = W(X, 1) = W(0, Y) = W(1, Y) = 0,
(33a)
(1W/1Y)(X, 0) = (1W/1Y)(X, 1) = (1W/1X)(0, Y) = (1W/1X)(1, Y) = 0.
(33b)
Applying the GDQ expression (11) to the boundary conditions (33), one has W1j = WNj = Wi1 = WiN = 0, Nx
Nx
Ny
Ny
k=1
k=1
k=1
k=1
(1) (1) (1) (1) s c1k Wkj = s cNk Wkj = s c1k Wik = s cNk Wik = 0,
Figure 6. The grid points of a rectangular plate.
(34a) (34b)
. .
290
T 5 Natural frequencies of a square plate (C–C–C–C)
Exact [15, 16] 8 × 8 [8] 7 × 7, uniform grid 9 × 9, uniform grid 7 × 7, non-uniform grid 9 × 9, non-uniform grid
v¯ 1
% error
v¯ 2
% error
v¯ 3
% error
v¯ 4
% error
35·99 35·78 36·83
— 0·55 2·3
73·41 — 65·03
— — 11
73·41 — 65·03
— — 11
108·3 — 90·12
— — 17
36·03
0·11
79·70
8·6
79·70
8·6
115·9
7·0
36·07
0·22
75·81
3·3
75·81
3·3
110·4
1·9
35·98
0·03
73·34
0·10
73·34
0·10
108·0
0·28
for i = 1, 2, . . . , Nx and j = 2, 3, . . . , Ny − 1. Again, after introducing the boundary conditions (34) into equations (32), the remaining eigenvalue problem can be solved to obtain the natural frequencies. For a plate with all four edges simply supported, the boundary conditions are W(X, 0) = W(X, 1) = W(0, Y) = W(1, Y) = 0,
(35a)
(1 2W/1Y 2)(X, 0) = (1 2W/1Y 2)(X, 1) = (1 2W/1X 2)(0, Y) = (1 2W/1X 2 )(1, Y) = 0.
(35b)
Applying GDQ (11) to the boundary conditions (35) yields W1j = WNj = Wi1 = WiN = 0, Nx
Nx
Ny
Ny
k=1
k=1
k=1
k=1
(36a)
(2) (2) (2) (2) s c1k Wkj = s cNk Wkj = s c1k Wik = s cNk Wik = 0,
(36b)
for i = 1, 2, . . . , Nx and j = 2, 3, . . . , Ny − 1. As before, the natural frequencies of a simply supported rectangular plate can be found by combining the boundary conditions (36) with equations (32). The natural frequencies of the first four modes are presented in Table 5 for a square plate with all the four edges clamped (C–C–C–C) and in Table 6 for a square plate with all four edges simply supported (S–S–S–S). In each case, results were obtained by using various numbers of grid points, both uniformly spaced and non-uniformly spaced. The T 6 Natural frequencies of a square plate (S–S–S–S)
Exact [15, 16] 7 × 7, uniform grid 9 × 9, uniform grid 7 × 7, non-uniform grid 9 × 9, non-uniform grid
v¯ 1
% error
v¯ 2
% error
19·74 19·83
— 0·46
49·35 42·21
— 14
49·35 42·21
— 14
78·96 66·59
— 16
19·74
0·0
51·85
5·1
51·85
5·1
81·85
3·7
19·75
0·05
49·17
0·36
49·17
0·36
79·55
0·75
19·74
0·0
49·68
0·67
49·68
0·67
79·45
0·62
v¯ 3
% error
v¯ 4
% error
291
non-uniformly spaced grid points were chosen as X, Y = 0·0, 0·08, 0·2, 0·5, 0·8, 0·92, 1·0 for Nx = Ny = 7, and X, Y = 0·0, 0·05, 0·1, 0·25, 0·5, 0·75, 0·9, 0·95, 1·0 for Nx = Ny = 9. From the results, the non-uniform grids with finer meshes close to boundaries were found to be more accurate than uniform grids. The natural frequencies for the first four modes are reasonably accurate with only 7 × 7 grids on the plate. The natural frequencies were obtained by solving the eigenvalue problem of a 9 × 9 matrix for the 7 × 7 grids or the eigenvalue problem of a 25 × 25 matrix for the 9 × 9 grids. Obviously, the computing time and the storage required are very small. The corresponding mode shapes for the first four modes, obtained by using 9 × 9 grids, are presented in Figures 7 and 8. In Table 7 are given the natural frequencies for the first four modes of two rectangular plates with different ratios of length to width, b = a/b = 0·4 and 1·5, respectively. Two sets of boundary conditions were considered for each plate. Exact results are also given for comparison. The results were obtained by using only nine grid points along both the X and Y direction (i.e., Nx = Ny = 9). Again, the results are found to be in very good agreement with exact solutions. The programming and formulating procedures are also found to be very straightforward and simple and the computing time required very small. 5. CONCLUSIONS
Deficiencies existing in the method of differential quadrature have been evaluated and discussed. An improved and generalized differential quadrature (GDQ) method which
Figure 7. The mode shapes of a square plate with all edges clamped.
. .
292
Figure 8. As Figure 7, but all edges simply supported.
overcomes the drawback has been introduced and applied to solve some vibration problems. GDQ provides a very simple algebraic formula to determine the weighting coefficients required by the differential quadrature approximation without restricting the choice of mesh grids. The numerical applications of GDQ to some vibration problems show that accurate results can be obtained by using relatively few grid points, and relatively very little storage and computing time. The discretizing and the programming procedures are straightforward and easy. The boundary conditions are easily incorporated in the T 7 Natural frequencies of rectangular plates (Nx = Ny = 9, non-uniformly spaced grids) Boundary
b = a/b
v¯ 1
% error
v¯ 2
% error
v¯ 3
% error
v¯ 4
% error
C–C–C–C
Exact GDQ Exact GDQ
(b = 0·4) (b = 0·4) (b = 1·5) (b = 1·5)
23·65 23·65 60·77 60·79
— 0·0 — 0·03
27·82 27·92 93·86 94·99
— 0·36 — 1·2
35·45 36·22 148·8 147·3
— 2·2 — 1·0
46·70 46·21 149·7 152·2
— 1·1 — 1·7
S–S–S–S
Exact GDQ Exact GDQ
(b = 0·4) (b = 0·4) (b = 1·5) (b = 1·5)
11·45 11·45 32·08 32·08
— 0·0 — 0·0
16·19 16·19 61·69 61·97
— 0·0 — 0·45
24·08 24·33 98·70 99·51
— 1·0 — 0·82
35·14 34·73 111·0 111·3
— 1·2 — 0·27
293
GDQ. The accuracy, efficiency and convenience of this method indicate a great potential for its wider use in vibration analysis. REFERENCES 1. R. E. B and J. C 1977 Journal of Mathematical Analysis and Applications 34, 235–238. Differential quadrature and long-term integration. 2. R. B, B. G. K and J. C 1972 Journal of Computational Physics 10, 40–52. Differential quadrature: a technique for the rapid solution of nonlinear partial differential equations. 3. R. B and B. G. K 1974 Mathematical Bioscience 19, 1–8. Solution of the partial differential equation of the Hodgkins–Huxley model using differential quadrature. 4. J. O. M 1979 Journal of Mathematical Analysis and Applications 71, 403–411. The method of differential quadrature for transient nonlinear diffusion. 5. G. N, R. E. B, K. M. W and E. S. L 1984 Journal of Mathematical Analysis and Applications 98, 220–235. Differential quadrature and partial differential equations: some numerical results. 6. F. C and C. M. S 1983 Journal of Mathematical Analysis and Applications 93, 206–221. Applications of differential quadrature to transport processes. 7. F. C and C. M. S 1983 International Journal of Numerical Methods in Engineering 19, 711–724. Solution of the Poisson equation by differential quadrature. 8. C. W. B, S. K. J and A. G. S 1988 American Institute of Aeronautics and Astronautics Journal 26, 612–618. Two new approximate methods for analyzing free vibration of structural components. 9. S. K. J, C. W. B and A. G. S 1989 International Journal of Numerical Methods in Engineering 28, 561–577. Application of differential quadrature to static analysis of structural components. 10. C. S 1991 Ph.D. Thesis, Department of Aerospace Engineering, University of Glasgow, U.K., Generalized differential-integral quadrature and application to the simulation of incompressible viscous flows including parallel computation. 11. C. S and B. E. R 1990 Proceedings of the Third International Conference on Advances in Numerical Methods in Engineering: Theory and Applications, Swansea, U.K. II, 978–985. High resolution of natural convection in a square cavity by generalized differential quadrature. 12. C. S and B. E. R 1992 International Journal of Numerical Methods in Fluid Dynamics 15, 791–798. Application of generalized differential quadrature to solve two-dimensional incompressible Navier–Stokes equations. 13. D. Y 1950 Journal of Applied Mechanics 17, 448–453. Vibration of rectangular plates by the Ritz method. 14. S. T and D. H. Y 1955 Vibration Problems in Engineering. New York: Van Nostrand Reinhold Company; (Third edition). 15. R. D. B 1984 Formulas for Natural Frequency and Mode Shapes. Malabur, Florida: Robert E. Krieger. 16. A. W. L 1973 Journal of Sound and Vibration 31, 257–293. The free vibration of rectangular plates.
APPLICATION OF GENERALIZED DIFFERENTIAL QUADRATURE TO VIBRATION ANALYSIS H. D, M. K. L R. M. L School of Mechanical and Production Engineering, Nanyang Technological University, Singapore 2263 (Received 12 February 1993, and in final form 4 January 1994) This paper presents a first endeavour to exploit the generalized differential quadrature method as an accurate, efficient and convenient numerical technique for vibration analysis. The generalized differential quadrature (GDQ) overcomes some major drawbacks existing in the method of differential quadrature (DQ). In GDQ, the weighting coefficients for derivative approximations are given by a simple algebraic expression and a recurrence relationship, together with arbitrary choice of grid points. GDQ is then applied to solve some free vibration problems. The accuracy, efficiency and convenience of GDQ are demonstrated throughout the numerical examples.
1. INTRODUCTION
Often, numerical approximation methods have to be sought to solve vibration problems due to the complexity of the problems. Classical techniques such as finite element and finite difference methods are well developed and well known. These methods can provide very accurate results by using a large number of grid points; thus they are computationally expensive. In a large number of cases, only a limited number of frequencies and mode shapes or dynamic responses at only a limited number of points is needed to be found. However, for acceptable accuracy, conventional numerical methods still require the use of a large number of grid points to obtain accurate solutions. Consequently, the requirement for computer capacity is often unncessarily large in such cases. An alternative method called the method of differential quadrature may be used for such purposes. The method of differential quadrature was first introduced by Bellman et al. [1, 2] for numerically solving partial differential equations using considerably fewer grid points. This method is based on the idea that the derivative of a function with respect to a coordinate direction can be expressed as a weighted linear sum of all the function values at all mesh points along that direction. Applications of this method to various problems showed that it has potential as an attractive approximation technique [3–9]. However, there exist some major drawbacks to be overcome in the originally proposed differential quadrature. As is known, the most important part of differential quadrature is to determine the weighting coefficients for discretization of any order partial derivative. There were two methods in use which were proposed by Bellman et al. [2] to obtain the weighting coefficients. One method is to solve a set of algebraic equations which satisfy exactly the linear constrained relation for all polynomials of degree less than or equal to N − 1. This set of equations has a unique solution because the matrix elements are composed of a Vandermonde matrix. Unfortunately, when N is large the Vandermonde matrix is ill-conditioned and the inversion of this matrix becomes difficult. Moreover, a set of N × N linear algebraic equations has to be solved for each order derivative even if 279 0022–460X/95/120279 + 15 $08.00/0
7 1995 Academic Press Limited
. .
280
the equations are solvable. The other method is to compute the weighting coefficients by a simple algebraic formula, but with the coordinates of grid points chosen as the roots of an Nth order shifted Legendre polynomial. This means that if N is specified, the distribution of grid points is fixed even for different physical problems or different boundary conditions. This also creates a major drawback and restricts the application of differential quadrature, since some practical problems may need more grid points near the boundary, while some others may not. Obviously, both of the methods originally proposed for determining the weighting coefficients in the method of differential quadrature have some major drawbacks. In order to overcome these drawbacks, Shu presented a generalized differential quadrature in connection with solving partial differential equations in fluid mechanics [10]. Preliminary results in solving some fluid mechanics problems have shown that the method is very efficient and convenient [11, 12]. The objective of this paper is to explore the great potential of the generalized differential quadrature method as an accurate, efficient and convenient numerical method for vibration analysis. The drawbacks existing in the original differential quadrature (DQ) method are first summarized and quantified in relation to its application to vibration analysis. Then, the generalized differential quadrature method (GDQ) is introduced. The generalized differential quadrature has overcome the possible singularity problem of the original differential quadrature in obtaining the weighting coefficients. Furthermore, in GDQ, the weighting coefficients for the first derivative and for the higher order derivatives are given by a simple algebraic expression and a recurrence relationship respectively, and without any restriction on arbitrary choice of grid points. GDQ is finally applied to solve for natural frequencies of some structures under various boundary conditions. The ease of use and the accuracy of the GDQ are demonstrated through the numerical examples. 2. DIFFERENTIAL QUADRATURE
The method of differential quadrature is based on the idea that the partial derivative of a function with respect to a space variable at a given discrete point can be expressed as a weighted linear sum of the function values at all discrete points in the domain of that variable. Let us take the first derivative of a one dimensional function u(x, t) as an example: the higher order partial derivatives will have essentially the same formation. A differential quadrature approximation of the first derivative of the function u(x, t) at the ith discrete point on a grid is N
ux (xi , t) = s cij(1) u(xj , t)
for i = 1, 2, . . . , N,
(1)
j=1
where ux (xi , t) is the first derivative of u(x, t) with respect to x at xi , and N is the number of discrete grids. cij(1) are the weighting coefficients for the first derivative approximation, which have to be pre-determined. There are two approaches used in the method of differential quadrature to determine the weighting coefficients cij(1) [2]. The first one is to let equation (1) be exact for all polynomials of degree less than or equal to (N − 1), g(x) = x k, k = 0, 1, . . . , N − 1. Then, one has a set of N × N linear algebraic equations: N
s cij xjk = kxik − 1 j=1
for k = 0, 1, . . . , N − 1
and
i = 1, 2, . . . , N.
(2)
281
This set of equations has a unique solution, since its matrix is of Vandermonde form. Unfortunately, it has been found that this set of equations becomes ill-conditioned and is difficult to solve when N is large. In order to quantify this singularity, weighting coefficients have been calculated for equal space grids based on equation (2) for various numbers of grid points N. From the computation, it is found that the maximum number of grid points is N = 22. Once the grid number if greater than 22, the set of linear algebraic equations become singular and cannot be solved at all. The computed weighting coefficients are also compared with the accurate coefficients obtained from the generalized differential quadrature method to be introduced later; it is also found that the results from equations (2) have some errors when the grid number is more than 17. Therefore, the maximum number of grid points is practically 17 for equal space grids if this method is used for determining the weighting coefficients. In addition, one has to solve a set of N × N linear equations for each order of the derivatives in the governing equations even when the equations are solvable (N E 17). The other approach to determine the weighting coefficients is similar to the first one with the exception that a different set of test functions g(x) is chosen for satisfying equation (1) exactly: g(x) = LN (x)/(x − xj )L'N (xj ) for j = 1, 2, . . . , N. (3) Here N is the number of grid points, LN (x) is the Nth order Legendre polynomial and L'N (x) is the first derivative of LN (x). By choosing xj to be the roots of the shifted Legendre polynomial and substituting equation (3) into equation (1), Bellman et al. [2] obtained a direct simple algebraic expression for the weighting coefficients cij(1) = L'N (xi )/(xi − xj )L'N (xj ) (1) ii
c = (1 − 2xi )/2xi (xi − 1),
for i $ j; i, j = 1, 2, . . . , N.
(4a, b)
In this approach, the weighting coefficients are easy to obtain without solving algebraic equations or having a singularity problem. However, it is obvious that once the number of grid points N is specified, the roots of the shifted Legendre polynomial are given; thus the distribution of the grid points is fixed no matter what problems are considered. This puts a major restriction on this method to problems in structural dynamic analysis, since all sorts of boundary conditions could appear and different mesh grids may be needed for different boundary conditions and different structure geometry. From the above discussion, some deficiencies have been identified in the originally proposed differential quadrature. These deficiences put some restrictions on the application of the DQ method to structural dynamics problems. This is probably one of the reasons that the method of differential quadrature is not widely used in structural analysis. In order to overcome the deficiencies, a generalized differential quadrature, which was recently proposed by Shu [10] for solving partial differential equations in fluid mechanics, will be introduced and applied to solve some problems in vibration analysis. The ease of use and the accuracy of the GDQ will be demonstrated through numerical examples. 3. GENERALIZED DIFFERENTIAL QUADRATURE
As pointed out above, two approaches have been proposed by previous researchers for the differential quadrature method for determining the weighting coefficients. Both of them have some drawbacks. The first restricts one to a small number of grids points and the need to solve sets of linear equations. The second limits the distribution of the grid points, which is critical in structural dynamic analysis. To remedy these deficiencies, what one
. .
282
wants is to find a method to determine the weighting coefficients without any limitation on the choice of grid meshes and which yet provides a simple algebraic expression to determine the weighting coefficients. Such a method was recently proposed by Shu [10] in relation to solving some partial differential equations in fluid mechanics. In order to find a simple algebraic expression for calculating the weighting coefficients without restricting the choice of grid meshes, Shu chose Lagrange interpolated polynomials as the set of tests functions g(x) instead of using power polynomials or Legendre polynomials: gi (x) = M(x)/(x − xi )M (1)(xi ) for i = 1, 2, . . . , N. (5) Here N
M(x) = t (x − xj ),
(6)
j=1
M (1)(x) is the first derivative of M(x), N
M (1)(xi ) = t
(xi − xj ),
(7)
j = 1, j $ i
and N is the number of grid points. When equation (1) is satisfied for all the test functions gi (x), a simple algebraic expression can be obtained to determine the weighting coefficients: cij(1) = M (1)(xi )/(xi − xj ) · M (1)(xj ) (1) ii
(2)
(1)
c = M (xi )/2M (xi ),
for i $ j;
i, j = 1, 2, . . . , N.
(8a, b)
Equations (8) provide simple expressions for computing cij(1) without any restriction in choice of the co-ordinates of the grid points xi . It is obvious that once the grid points xi are given, M (1)(x) is very easily obtained from equation (7). Hence, cij(1) can be easily calculated for i $ j. The calculation of cii(1) is based on the calculation of the second derivative of M(x) which is more difficult to obtain. Instead of using equation (8b), a more convenient relationship can be obtained and used for calculating cii(1) . It can be shown by using a Taylor series expansion that the following relationship exists for cij(1) : N
s cij(1) = 0
for i = 1, 2, . . . , N.
(9)
j=1
Thus, from equation (9), the coefficient cii(1) can be calculated from cij(1) (i $ j): that is, N
cii(1) = − s
cij(1)
for i = 1, 2, . . . , N.
(10)
j = 1,j $ i
For the approximation of higher order derivatives, the generalized differential quadrature gives a form similar to equation (1), N
ux(m) (xi , t) = s cij(m) u(xj , t)
for i = 1, 2, . . . , N,
m = 2, 3, . . . , N − 1,
(11)
j=1
where ux(m) (xi , t) is the mth derivative of u(x, t) with respect to x at xi and N is the number of discrete grid points. cij(m) are the weighting coefficients for the mth derivative approximation.
283
The weighting coefficients for second and higher order derivatives can be similarly obtained. Amazingly, a recurrence relationship has been found for the mth order weighting coefficients cij(m) :
0
cij(m) = m cii(m − 1)cij(1) −
cij(m − 1) xi − xj
1
for i $ j,
m = 2, 3, . . . , N − 1,
i, j = 1, 2, . . . , N.
(12) Thus the cij(m) can be derived from the (m − 1)th order weighting coefficients cij(m − 1) . The cii(m) can be obtained from a relationship similar to equation (10): N
cii(m) = − s
cij(m)
for i = 1, 2, . . . , N.
(13)
j = 1,j $ i
Thus equations (12) and (13) together with equation (8a) and equation (10) give a convenient and general form for determining the weighting coefficients for the derivatives of orders one through N − 1. There are no restrictions on the co-ordinates of the chosen grid points. There is no need to solve for the weighting coefficients from a set of algebraic equations which could be ill-conditioned when the number of grids is large. Furthermore, this set of expressions for the determination of the weighting coefficients is so compact and simple that it is very easy to be implemented in formulating and programming because of the recurrence feature. All these features are very convenient for solving practical problems in structural dynamic analysis. Application of this generalized differential quadrature method to time dependent dynamic problems will result in a set of ordinary differential equations with the time dependent function values at the grid points as unknowns. This leads to an eigenvalue problem for free vibration analysis or to a set of time dependent ordinary differential equations for forced vibration problems. The natural frequencies and mode shapes or time history response can then be obtained. Finally, once the function values at all grid points have been obtained, it is very easy to determine the function values in the overall domain in terms of a polynomial approximation: that is, N
u(x) = s u(xi )gi (x),
(14)
i=1
where gi (x) are the Lagrange interpolated polynomials of equation (5). 4. APPLICATION OF GDQ TO VIBRATION ANALYSIS
The method of generalized differential quadrature is used to analyze the free vibration problems of beams and plates in this section. The first case is calculating the natural frequencies of slender beams. The second and third cases are calculating the natural frequencies of a circular and a rectangular plate, respectively. The formulating and programming procedures are shown to be very straightforward and simple. The boundary conditions are easy to be implemented. 4.1. The governing equation of a Bernoulli–Euler beam in bending is
$
%
12 1 2w 1 2w = 0, 2 EI 2 + rA 1x 1t 2 1x
0 Q x Q L,
(15)
where EI is the beam’s flexural rigidity, rA is the mass per unit length and L is the length of the beam.
. .
284
For free vibration, the transverse displacement w is assumed as: w(x, t) = W(x) eivt.
(16)
Substituting expression (16) into equation (15) and normalizing the equation yields 1 2s(X) 1 2W 1s(X) 1 3W 1 4W = v¯ 2W, 2 2 +2 3 + s(X) 1X 1X 1X 1X 1X 4
(17)
where s(X) = EI/EI0 , v¯ 2 = (rAL 4/EI0 )v 2 and X = x/L. Given the number of grid points N (see Figure 1), and applying the generalized differential quadrature approximation (11) to equation (17) at each discrete point on the grid, one obtains 1 2s(Xi ) 1X 2
0
1
N
1s(Xi ) 1X
s cij(2) Wj + 2 j=1
0
1
N
0
1
N
s cij(3) Wj + s(Xi ) s cij(4) Wj = v¯ 2Wi , j=1
j=1
(18)
for i = 1, 2, . . . , N. Equations (18) give an eigenvalue problem. Solving this eigenvalue problem together with appropriate boundary conditions, one can obtain the natural frequencies of the structure. Consider a beam with left end fixed and right end free. The boundary conditions then are W = 1W/1X = 0 at X = 0, (19a) EI 1 2W/1x 2 = 0,
(1/1x)(EI 1 2W/1x 2) = 0
at x = L,
or s(X) 1 2W/1X 2 = 0,
{1s(X)/1X} 1 2W/1X 2 + s(X) 1 3W/1X 3 = 0
at X = 1.
(19b)
Applying the differential quadrature approximation (11) to the boundary conditions (19) yields N
W1 = 0,
s cij(1) Wj = 0,
(20a)
j=1
0
N
1
(2) s(XN ) s cNj Wj = 0, j=1
0
N
1
0
N
1
(2) (3) {1s(XN )/1X} s cNj Wj + s(XN ) s cNj Wj = 0. j=1
j=1
(20b)
In order to introduce the boundary conditions into equations (18), first equations (20) are used to solve for W1 , W2 , WN − 1 and WN in terms of the variables W3 , W4 , . . . , WN − 2 . The expressions for W1 , W2 , WN − 1 and WN in terms of the variables W3 , W4 , . . . , WN − 2 are then substituted into equations (18) to eliminate the variables W1 , W2 , WN − 1 and WN , and only the discretized equations at the points i = 3, 4, . . . , N − 2 are to be used in equations (18). Finally, by solving the remaining eigenvalue problem, the natural frequencies of the structure can be obtained. For a beam with its left end fixed and its right end simply supported, the boundary conditions are W = 1W/1X = 0
at X = 0,
W = s(X) 1 2W/1X 2 = 0
Figure 1. The grid points of a slender beam.
at X = 1.
(21a, b)
285
T 1 Natural frequencies of a uniform beam (clamped–free)
Exact [15] N = 8 [8] N = 8, uniform grid N = 10 uniform grid N = 12, uniform grid N = 14, uniform grid N = 12, non-uniform grid
v¯ 1
% error
v¯ 2
3·5160 3·5226 3·5224
— 0·19 0·18
22·035 — 20·837
3·5161
0·003
3·5160
% error
v¯ 3
% error
v¯ 4
% error
— — 5·4
61·697 — —
— — —
120·90 — —
— — —
22·213
0·81
—
—
—
—
0·0
22·025
0·045
61·298
0·65
131·93
9·1
3·5160
0·0
22·035
0·0
61·747
0·081
116·75
3·4
3·5160
0·0
22·033
0·009
61·604
0·15
125·63
3·9
Applying the differential quadrature approximation (11) to the boundary conditions (21) gives N
s cij(1)Wj = 0,
W1 = 0,
(22a)
j=1
WN = 0,
0
N
1
(2) s(XN ) s cNj Wj = 0. j=1
(22b)
Again, equations (22) can be used to eliminate the variables W1 , W2 , WN − 1 and WN in equations (18), as discussed above. The natural frequencies are then solved from the remaining eigenvalue problem. Calculations have been done for both a uniform beam s(X) = 1 and a non-uniform beam under two sets of different boundary conditions, namely clamped–free (C–F), and clamped – simply supported (C–S). The natural frequencies for the first four modes of a uniform beam are presented in Table 1 for clamped–free support and in Table 2 for clamped – simply supported support together with exact solutions. The first natural frequency obtained by Bert et al. [8] using the original differential quadrature method is also presented in Table 1. The present result for N = 8 is found to be consistent with it. T 2 Natural frequencies of a uniform beam (clamped – simply supported) v¯ 1 Exact [15] N = 8, uniform grid N = 10, uniform grid N = 12, uniform grid N = 14, uniform grid N = 12, non-uniform grid
% error
v¯ 2
% error
v¯ 3
% error
v¯ 4
% error
15·418 15·637
— 1·4
49·965 54·430
— 8·9
104·25 —
— —
178·27 —
— —
15·410
0·052
49·636
0·66
—
—
—
—
15·419
0·006
49·978
0·026
99·128
4·9
153·68
13·8
15·418
0·0
49·965
0·0
105·17
0·88
187·03
4·9
15·418
0·0
49·966
0·002
103·80
0·43
174·21
2·3
286
. .
Figure 2. The mode shapes of a clamped–free beam.
Results were obtained by using various numbers of grid points. From the tables, the convergence of this method is seen to be very good. Accurate results can be achieved by using very few grid points. Accurate natural frequencies for the first two modes can be obtained by using ten uniformly spaced grid points, while 14 uniformly spaced grid points are sufficient to give accurate frequencies for the first four modes. It is especially interesting to note that the accuracy can be improved by using non-uniform meshes with finer grids closer to the boundaries. The non-uniformly spaced grid points for N = 12 are chosen as X = 0·0, 0·05, 0·1, 0·15, 0·25, 0·4, 0·6, 0·75, 0·85, 0·9, 0·95, 1·0. With only 12 points, one can obtain very good results for the first four modes. The computing time is tiny since one only has to solve an eigenvalue problem of an 8 × 8 matrix after imposing the boundary conditions. To demonstrate better the accuracy of the solutions, the mode shapes of the first four modes were also obtained for both sets of boundary conditions by using 12 non-uniformly spaced grid points. The mode shapes are presented in Figures 2 and 3.
Figure 3. The mode shapes of a clamped – simply supported beam.
287
T 3 Natural frequencies of a non-uniform beam (N = 12 non-uniformly spaced grids) s(X) = 1 + 0·5X (C–F) s(X) = 1 − X (C–F) s(X) = 1 + 0·5X (C–S) s(X) = 1 − X (C–S)
v¯ 1
v¯ 2
v¯ 3
v¯ 4
3·6733 3·2746 16·928 10·739
24·079 16·871 55·409 31·693
67·978 43·700 115·48 65·311
136·17 82·761 191·20 110·92
In Table 3 are presented the natural frequencies of the first four modes for two non-uniform beams with different stiffness distributions, s(X) = 1 + 0·5X and s(X) = 1 − X, respectively. Again, two sets of different boundary conditions (C–F and C–S) were considered for each beam. Twelve non-uniformly spaced grid points were chosen in the calculation. The varying cross section stiffness s(X) of the beams can be very easily implemented from equations (18). The programming is also easy and straightforward for such more complex problems. The computing effort is still small since one has to solve an eigenvalue problems of a matrix of dimension 8 × 8 only. 4.2. The governing equation for a thin circular plate of isotropic material and undergoing axisymmetric motion is 1 4w 2 1 3w 1 1 2w 1 1w rh 1 2w + + − + = 0. 1r 4 D 1t 2 r 1r 3 r 2 1r 2 r 3 1r
(23)
Substituting equation (16) into equation (23) and normalizing it yields 1 4W 2 1 3W 1 1 2W 1 1W + − + − v¯ 2W = 0, 1R 4 R 1R 3 R 2 1R 2 R 3 1R
(24)
where R = r/a, v¯ 2 = rha 4v 2/D, a is the outside radius of the plate, h the thickness and D the flexural rigidity. The regularity condition at the centre of the plate is 1W/1R = 0
at R = 0.
(25)
The regularity condition is necessary to assure that the plate slope is zero at the origin to avoid a singularity at this location. The boundary conditions for a clamped outside edge are W = 0, 1W/1R = 0 at R = 1. (26) Applying the generalized differential quadrature approximation (11) to equation (24) at each discrete point on the grid (see Figure 4) yields N
s cij4 Wj + j=1
2 N (3) 1 N 1 N s cij Wj − 2 s cij(2) Wj + 3 s cij(1)Wj = v¯ 2Wi Ri j = 1 R1 j = 1 R1 j = 1
for i = 1, 2, . . . , N, (27)
and applying the generalized differential quadrature approximation (11) to the boundary conditions (25) and (26) yields N
s c1j(1) Wj = 0 j=1
N
and
WN = 0,
(1) s cNj Wj = 0.
(28a, b)
j=1
Again, using equations (28) to eliminate the variables most close to the boundaries, W1 , WN − 1 and WN , in equations (27) as discussed in the previous example, one can solve the remaining eigenvalue problem to obtain the natural frequencies.
. .
288
Figure 4. The grid points of a circular plate.
The natural frequencies obtained for the first four modes are given in Table 4, together with exact solutions available for comparison. The corresponding mode shapes obtained by using 12 non-uniformly spaced grid points are plotted in Figure 5. The first natural frequency obtained by Bert et al. [8] using the original differential quadrature method is also presented in Table 4. The present result for N = 8 is found to be consistent with it. As observed from Table 4, the first four natural frequencies obtained by using 12 non-uniformly spaced grid points are very accurate. The computing time is very small for solving the eigenvalues of the resulting 8 × 8 matrix. 4.3. The governing equation for a thin rectangular plate of isotropic material is given as: D 1 4w/1x 4 + 2D 1 4w/1x 2 1y 2 + D 1 4w/1y 4 + rh 1 2w/1t 2 = 0,
(29)
where r is the density of the material, h the thickness of the plate and D the flexural rigidity of the plate. The transverse displacement w for free vibration is taken as w(x, y, t) = W(x, y) eivt.
(30)
Substituting equation (30) into equation (29), one obtains the normalized equation 1 4W/1X 4 + 2b 2 1 4W/1X 2 1Y 2 + b 4 1 4W/1Y 4 − v¯ 2W = 0,
(31)
T 4 Natural frequencies of a circular plate
Exact [15] N = 8 [8] N = 8, uniform grid N = 12, non-uniform grid
v¯ 1
% error
v¯ 2
% error
v¯ 3
% error
v¯ 4
% error
10·22 10·20 10·20
— 0·19 0·19
39·77 — 44·68
— — 12
89·10 — 80·24
— — 9·9
158·2 — —
— — —
10·22
0·0
39·77
89·22
0·13
156·6
1·0
0·0
289
Figure 5. The mode shapes of a circular plate with outside edge clamped.
where b = a/b, v¯ 2 = rha 4v 2/D, X = x/a, Y = y/b, a is the length of the plate and b its width. Applying the generalized differential quadrature approximation (11) to equation (31) at each discrete point on the grid (see Figure 6) gives Nx
Ny
Nx
Ny
k=1
m=1
k=1
k=1
(2) s cik(4) Wkj + 2b 2 s cjm s cik(2) Wkm + b 4 s cjk(4) Wik = v¯ 2Wij ,
i = 1, 2, . . . , Nx ,
j = 1, 2, . . . , Ny ,
(32)
where Nx and Ny are the number of grid points along the X and Y directions, respectively (see Figure 6). The boundary conditions for a plate clamped on all four edges are: W(X, 0) = W(X, 1) = W(0, Y) = W(1, Y) = 0,
(33a)
(1W/1Y)(X, 0) = (1W/1Y)(X, 1) = (1W/1X)(0, Y) = (1W/1X)(1, Y) = 0.
(33b)
Applying the GDQ expression (11) to the boundary conditions (33), one has W1j = WNj = Wi1 = WiN = 0, Nx
Nx
Ny
Ny
k=1
k=1
k=1
k=1
(1) (1) (1) (1) s c1k Wkj = s cNk Wkj = s c1k Wik = s cNk Wik = 0,
Figure 6. The grid points of a rectangular plate.
(34a) (34b)
. .
290
T 5 Natural frequencies of a square plate (C–C–C–C)
Exact [15, 16] 8 × 8 [8] 7 × 7, uniform grid 9 × 9, uniform grid 7 × 7, non-uniform grid 9 × 9, non-uniform grid
v¯ 1
% error
v¯ 2
% error
v¯ 3
% error
v¯ 4
% error
35·99 35·78 36·83
— 0·55 2·3
73·41 — 65·03
— — 11
73·41 — 65·03
— — 11
108·3 — 90·12
— — 17
36·03
0·11
79·70
8·6
79·70
8·6
115·9
7·0
36·07
0·22
75·81
3·3
75·81
3·3
110·4
1·9
35·98
0·03
73·34
0·10
73·34
0·10
108·0
0·28
for i = 1, 2, . . . , Nx and j = 2, 3, . . . , Ny − 1. Again, after introducing the boundary conditions (34) into equations (32), the remaining eigenvalue problem can be solved to obtain the natural frequencies. For a plate with all four edges simply supported, the boundary conditions are W(X, 0) = W(X, 1) = W(0, Y) = W(1, Y) = 0,
(35a)
(1 2W/1Y 2)(X, 0) = (1 2W/1Y 2)(X, 1) = (1 2W/1X 2)(0, Y) = (1 2W/1X 2 )(1, Y) = 0.
(35b)
Applying GDQ (11) to the boundary conditions (35) yields W1j = WNj = Wi1 = WiN = 0, Nx
Nx
Ny
Ny
k=1
k=1
k=1
k=1
(36a)
(2) (2) (2) (2) s c1k Wkj = s cNk Wkj = s c1k Wik = s cNk Wik = 0,
(36b)
for i = 1, 2, . . . , Nx and j = 2, 3, . . . , Ny − 1. As before, the natural frequencies of a simply supported rectangular plate can be found by combining the boundary conditions (36) with equations (32). The natural frequencies of the first four modes are presented in Table 5 for a square plate with all the four edges clamped (C–C–C–C) and in Table 6 for a square plate with all four edges simply supported (S–S–S–S). In each case, results were obtained by using various numbers of grid points, both uniformly spaced and non-uniformly spaced. The T 6 Natural frequencies of a square plate (S–S–S–S)
Exact [15, 16] 7 × 7, uniform grid 9 × 9, uniform grid 7 × 7, non-uniform grid 9 × 9, non-uniform grid
v¯ 1
% error
v¯ 2
% error
19·74 19·83
— 0·46
49·35 42·21
— 14
49·35 42·21
— 14
78·96 66·59
— 16
19·74
0·0
51·85
5·1
51·85
5·1
81·85
3·7
19·75
0·05
49·17
0·36
49·17
0·36
79·55
0·75
19·74
0·0
49·68
0·67
49·68
0·67
79·45
0·62
v¯ 3
% error
v¯ 4
% error
291
non-uniformly spaced grid points were chosen as X, Y = 0·0, 0·08, 0·2, 0·5, 0·8, 0·92, 1·0 for Nx = Ny = 7, and X, Y = 0·0, 0·05, 0·1, 0·25, 0·5, 0·75, 0·9, 0·95, 1·0 for Nx = Ny = 9. From the results, the non-uniform grids with finer meshes close to boundaries were found to be more accurate than uniform grids. The natural frequencies for the first four modes are reasonably accurate with only 7 × 7 grids on the plate. The natural frequencies were obtained by solving the eigenvalue problem of a 9 × 9 matrix for the 7 × 7 grids or the eigenvalue problem of a 25 × 25 matrix for the 9 × 9 grids. Obviously, the computing time and the storage required are very small. The corresponding mode shapes for the first four modes, obtained by using 9 × 9 grids, are presented in Figures 7 and 8. In Table 7 are given the natural frequencies for the first four modes of two rectangular plates with different ratios of length to width, b = a/b = 0·4 and 1·5, respectively. Two sets of boundary conditions were considered for each plate. Exact results are also given for comparison. The results were obtained by using only nine grid points along both the X and Y direction (i.e., Nx = Ny = 9). Again, the results are found to be in very good agreement with exact solutions. The programming and formulating procedures are also found to be very straightforward and simple and the computing time required very small. 5. CONCLUSIONS
Deficiencies existing in the method of differential quadrature have been evaluated and discussed. An improved and generalized differential quadrature (GDQ) method which
Figure 7. The mode shapes of a square plate with all edges clamped.
. .
292
Figure 8. As Figure 7, but all edges simply supported.
overcomes the drawback has been introduced and applied to solve some vibration problems. GDQ provides a very simple algebraic formula to determine the weighting coefficients required by the differential quadrature approximation without restricting the choice of mesh grids. The numerical applications of GDQ to some vibration problems show that accurate results can be obtained by using relatively few grid points, and relatively very little storage and computing time. The discretizing and the programming procedures are straightforward and easy. The boundary conditions are easily incorporated in the T 7 Natural frequencies of rectangular plates (Nx = Ny = 9, non-uniformly spaced grids) Boundary
b = a/b
v¯ 1
% error
v¯ 2
% error
v¯ 3
% error
v¯ 4
% error
C–C–C–C
Exact GDQ Exact GDQ
(b = 0·4) (b = 0·4) (b = 1·5) (b = 1·5)
23·65 23·65 60·77 60·79
— 0·0 — 0·03
27·82 27·92 93·86 94·99
— 0·36 — 1·2
35·45 36·22 148·8 147·3
— 2·2 — 1·0
46·70 46·21 149·7 152·2
— 1·1 — 1·7
S–S–S–S
Exact GDQ Exact GDQ
(b = 0·4) (b = 0·4) (b = 1·5) (b = 1·5)
11·45 11·45 32·08 32·08
— 0·0 — 0·0
16·19 16·19 61·69 61·97
— 0·0 — 0·45
24·08 24·33 98·70 99·51
— 1·0 — 0·82
35·14 34·73 111·0 111·3
— 1·2 — 0·27
293
GDQ. The accuracy, efficiency and convenience of this method indicate a great potential for its wider use in vibration analysis. REFERENCES 1. R. E. B and J. C 1977 Journal of Mathematical Analysis and Applications 34, 235–238. Differential quadrature and long-term integration. 2. R. B, B. G. K and J. C 1972 Journal of Computational Physics 10, 40–52. Differential quadrature: a technique for the rapid solution of nonlinear partial differential equations. 3. R. B and B. G. K 1974 Mathematical Bioscience 19, 1–8. Solution of the partial differential equation of the Hodgkins–Huxley model using differential quadrature. 4. J. O. M 1979 Journal of Mathematical Analysis and Applications 71, 403–411. The method of differential quadrature for transient nonlinear diffusion. 5. G. N, R. E. B, K. M. W and E. S. L 1984 Journal of Mathematical Analysis and Applications 98, 220–235. Differential quadrature and partial differential equations: some numerical results. 6. F. C and C. M. S 1983 Journal of Mathematical Analysis and Applications 93, 206–221. Applications of differential quadrature to transport processes. 7. F. C and C. M. S 1983 International Journal of Numerical Methods in Engineering 19, 711–724. Solution of the Poisson equation by differential quadrature. 8. C. W. B, S. K. J and A. G. S 1988 American Institute of Aeronautics and Astronautics Journal 26, 612–618. Two new approximate methods for analyzing free vibration of structural components. 9. S. K. J, C. W. B and A. G. S 1989 International Journal of Numerical Methods in Engineering 28, 561–577. Application of differential quadrature to static analysis of structural components. 10. C. S 1991 Ph.D. Thesis, Department of Aerospace Engineering, University of Glasgow, U.K., Generalized differential-integral quadrature and application to the simulation of incompressible viscous flows including parallel computation. 11. C. S and B. E. R 1990 Proceedings of the Third International Conference on Advances in Numerical Methods in Engineering: Theory and Applications, Swansea, U.K. II, 978–985. High resolution of natural convection in a square cavity by generalized differential quadrature. 12. C. S and B. E. R 1992 International Journal of Numerical Methods in Fluid Dynamics 15, 791–798. Application of generalized differential quadrature to solve two-dimensional incompressible Navier–Stokes equations. 13. D. Y 1950 Journal of Applied Mechanics 17, 448–453. Vibration of rectangular plates by the Ritz method. 14. S. T and D. H. Y 1955 Vibration Problems in Engineering. New York: Van Nostrand Reinhold Company; (Third edition). 15. R. D. B 1984 Formulas for Natural Frequency and Mode Shapes. Malabur, Florida: Robert E. Krieger. 16. A. W. L 1973 Journal of Sound and Vibration 31, 257–293. The free vibration of rectangular plates.