Differential quadrature Trefftz method for irregular plate problems

ARTICLE IN PRESS Engineering Analysis with Boundary Elements 33 (2009) 363–367

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

Differential quadrature Trefftz method for irregular plate problems Xu Liu, Xionghua Wu  Department of Applied Mathematics, Tongji University, Shanghai 200092, PR China

a r t i c l e in fo

abstract

Article history: Received 2 January 2008 Accepted 27 June 2008 Available online 6 September 2008

Differential quadrature Trefftz method (DQTM) is developed to deal with plate problems defined in irregular domains. DQTM divides the solution into two parts, a particular solution for inhomogeneous biharmonic equation and the general solution for homogeneous biharmonic equation. For the former, differential quadrature method based on the interpolation of the highest derivative (DQIHD) is involved. For the latter, polynomial basis functions are adopted instead of fundamental solutions. We will show that DQTM not only keeps the advantages of traditional differential quadrature method (DQM), high efficiency and accuracy, but also has no difficulties to deal with geometrically irregular domains. & 2008 Elsevier Ltd. All rights reserved.

Keywords: DQTM DQIHD Polynomial basis functions Particular solution General solution

1. Introduction

2. Methods

Though bending problems of plates have been extensively studied, it is still very difficult to deal with arbitrary shapes and boundary conditions. In this paper, we try to offer a general method based on differential quadrature Trefftz method (DQTM). It combines differential quadrature method (DQM) [1–3] and Trefftz method [4–7], and the essence is to divide the solution into a particular solution for inhomogeneous biharmonic equation and the general solution for homogeneous biharmonic equation. For irregular plates, a large enough rectangular domain containing the original one is set as the computational domain. For a particular solution without boundary conditions, DQM based on the interpolation of the highest derivative (DQIHD) [8] is used, which is different from the traditional DQM. Because it does not involve the numerical differentiation process, the accuracy can be obviously heightened. For the general solution, we will express it by the linear combination of polynomial basis functions instead of fundamental solutions, so that there is no singularity [9,10]. Then the collocation method is used to determine the unknown coefficients. On the one hand, this method can keep the advantages of DQM, such as high accuracy, efficiency and good convergence [3]. On the other hand, it is especially efficient for the plates whose shapes are irregular or whose boundary conditions are complex.

2.1. The linear bending problems of thin plates

 Corresponding author. Tel.: +86 21 65980478; fax: +86 21 65982341.

E-mail address: [email protected] (X. Wu). 0955-7997/$ - see front matter & 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.enganabound.2008.06.008

According to Kirchhoff’s three assumptions [11], the governing equation is a biharmonic equation. Here consider the standard form as follows:

D2 u ¼ f ,

(2.1)

where f ¼ q(x, y)/D, q(x, y) is a known transverse load function, D is the bending stiffness, u is the deflection w to solve. The boundary conditions in this paper involve clamped edges, simple supported edges and free edges, denoted by G1, G2 and G3, respectively. So the boundary of O is G ¼ qO ¼ G1+G2+G3. Here the boundary conditions are expressed as follows: 8 u ¼ g 1 ; on G1 [ G2 ; > > > > < Cu ¼ g 2 ; on G1 ; (2.2) Su ¼ g 3 ; on G2 [ G3 ; > > > > : Fu ¼ g ; on G3 ; 4 where g1, g2, g3 and g4 are known functions and generally equal to 0. C, S and F are three operators defined as follows:

q q þm , qx qy q2 q2 q2 2 2 S ¼ ðl þ nm2 Þ 2 þ ðnl þ m2 Þ 2 þ ð2  nÞlm , qxqy qx qy 3 q q3 2 F ¼  ½ð1  nÞm2 þ 1l 3 þ ð1  nÞð2l  m2  1Þm 2 qx qx qy 3 q q3 2 2  ½ð1  nÞl þ 1m 3 . þ ð1  nÞð2m2  l  1Þl 2 qxqy qy

C¼l

ARTICLE IN PRESS 364

X. Liu, X. Wu / Engineering Analysis with Boundary Elements 33 (2009) 363–367

It is easy to find that they are only related to the Poisson ratio of the plate n and l ¼ cos(n, x), m ¼ sin(n, x), where n is the outward normal direction.

0

2.2. DQIHD For the boundary value problems (2.1) and (2.2), DQTM decomposes the solution as u ¼ up þ uh ,

(2.3)

where up is a particular solution for (2.1) without considering the boundary conditions, so it is not unique. uh is the general solution for homogeneous biharmonic equation:

D2 u ¼ 0.

(2.4)

For a particular solution for (2.1), DQIHD is adopted here. It is introduced briefly in the following part and detailed in Ref. [8]. The essence of the traditional DQM is that the partial derivative of a function with respect to a variable can be approximated by a weighed sum of functional values at all discrete points in that direction. And the weighing coefficients depend only on the grid space. But DQIHD tries to use the values of the partial derivative at discrete points to express the functional values. Take the four-order ordinary differential equation for example. Firstly, u(4) can be approximated by the interpolation of the values at discrete grid points {xi} as uð4Þ ðxÞ ¼

n X

lj ðxÞuð4Þ ðxj Þ,

(2.5)

j¼1

where lj(x) is the basis function of Lagrange interpolation on the discrete points {xi}. In this paper, we choose Gauss–Chebyshev points, expressed as   2i  1 xi ¼ cos p ; 1pipn. (2.6) 2n Next, by integration we have Z x n X uð4Þ ds ¼ aj uð4Þ ðxj Þ þ uð3Þ ð0Þ, uð3Þ ¼ 0

where aj ¼ Denote

(2.7)

j¼1

Rx

0lj ds.

Repeat this process until the order is 0.

A ¼ ðaij Þ; B ¼ ðbij Þ; C ¼ ðcij Þ; D ¼ ðdij Þ, Z xi aj ðsÞds, aij ¼ aj ðxi Þ; bij ¼ bj ðxi Þ ¼ 0 Z xi Z bj ðsÞds; dij ¼ dj ðxi Þ ¼ cij ¼ cj ðxi Þ ¼ 0

0

cj ðsÞds.

U ð4Þ ¼ E U; U ð3Þ ¼ A U; U ð2Þ ¼ B U; U ð1Þ ¼ C U; U ¼ D U ,

(2.8)

where uð1Þ ð0Þ

uð2Þ ð0Þ

uð3Þ ð0Þ

uð0; 0Þ B u1;0 ð0; 0Þ B B B u2;0 ð0; 0Þ B B u ð0; 0Þ 3;0 U¼B B B u4;0 ðx1 ; 0Þ B B .. B . @ u4;0 ðxn ; 0Þ

u0;1 ð0; 0Þ u1;1 ð0; 0Þ

u0;2 ð0; 0Þ u1;2 ð0; 0Þ

u0;3 ð0; 0Þ u1;3 ð0; 0Þ

u0;4 ð0; y1 Þ u1;4 ð0; y1 Þ

 

u2;1 ð0; 0Þ u3;1 ð0; 0Þ u4;1 ðx1 ; 0Þ .. . u4;1 ðxn ; 0Þ

u2;2 ð0; 0Þ u3;2 ð0; 0Þ u4;2 ðx1 ; 0Þ .. . u4;2 ðxn ; 0Þ

u2;3 ð0; 0Þ u3;3 ð0; 0Þ u4;3 ðx1 ; 0Þ .. . u4;3 ðxn ; 0Þ

u2;4 ð0; y1 Þ u3;4 ð0; y1 Þ u4;4 ðx1 ; y1 Þ .. . u4;4 ðxn ; y1 Þ

   .. . 

1 u0;4 ð0; yn Þ u1;4 ð0; yn Þ C C C u2;4 ð0; yn Þ C C u3;4 ð0; yn Þ C C, C u4;4 ðx1 ; yn Þ C C C .. C . A u4;4 ðxn ; yn Þ

(2.10) where us,t(x, y) (0ps, tp4) is short for uxs yt ðx; yÞ. Generally, the boundary valuesuxs yt ðxi ; 0Þ anduxs yt ð0; yj Þ are given via the boundary conditions. But in this paper, they are also unknowns. DQIHD uses the equations as follows for a particular solution: Ex U DTy þ 2Bx U BTy þ Dx U ETy ¼ F .

(2.11)

It can be transformed into ðDy  Ex þ 2By  Bx þ Ey  Dx ÞVecðUÞ ¼ VecðFÞ,

(2.12)

where  is Kronecker-product, Vec(  ) means to reset the matrix as a column vector. (2.12) may have many solutions. The least-square method is used here to get one and according to the relationship between U and U xs yt , all derivatives are solved. Compared with the traditional DQM, DQIHD involves no numerical differentiation, which is very sensitive to even a small level of errors. It is not only for particular solutions. For bending moments and shearing forces which also involve derivatives, the process is almost the same. If the traditional DQM is chosen for particular solutions, numerical differentiation will be unavoidable, because we have to let C, S, F act on the particular solution up. So DQIHD can avoid numerical differentiation in two places.

2.3. General solution for homogeneous biharmonic equation The general solution uh can be expressed by the linear combination of a series of basis functions fuk gk¼1;2::: , i.e. X uh ¼ bk uk . (2.13) k

xi

Then the process can be expressed as

U ¼ ½ uð0Þ

integration instead. Here U is denoted as the value matrix of ðup Þx4 y4 and called the highest derivative of Up. U can be expressed as follows:

uð4Þ ðx1 Þ



uð4Þ ðxn Þ T Theorem 4.1. [10] A biharmonic function u(x, y) in a plane simply connected domain O certainly can be expressed as

E ¼ ½0; 0; 0; 0; E; A ¼ ½0; 0; 0; I; A; B ¼ ½0; 0; I; X; B, C ¼ ½0; I; X; X 2 =2; C; D ¼ ½I; X; X 2 =2; X 3 =6; D,

(2.9) (p)

And general solutions of bending moments and shearing forces can be expressed by linear combinations of relevant derivatives of basis functions, with the same coefficients. There can be different basis functions. Usually fundamental solutions are adopted. But here we choose the polynomial basis functions, because they have no singularity. They can be obtained similarly according to the following theorem.

where E is an identity matrix, the terms like u (0) (0ppp3) are the values or derivatives of u on the boundary points, I, X, X2/2, X3/6 are column vectors making up of the values of functions 1, x, x2/2, x3/6 at {xi}. Notice that matrices with underlines are no longer square. Extend to two-dimensional cases. Since the biharmonic operator contains three terms: q4/qx4, q4/qx2qy2 and q4/qy4 and the highest partial derivatives along x and y are both of four order, all of them can be expressed by ðup Þx4 y4 via numerical

uðx; yÞ ¼ Re½¯zjðzÞ þ cðzÞ,

(2.14)

where j(z) and c(z) are two analytic function, z ¼ x+yi; conversely, for two arbitrary analytic functions j(z) and c(z) in O, the function u(x, y) given by (2.14) is certainly a biharmonic function. Let O be a plane bounded domain. In order to obtain a series of simple biharmonic functions, we can first let j(z) and c(z) as polynomials of z, then substitute them into (2.14) and take the real parts. Distinctly, all of the harmonic functions are also biharmonic

ARTICLE IN PRESS X. Liu, X. Wu / Engineering Analysis with Boundary Elements 33 (2009) 363–367

functions. It is the case j(z) ¼ 0. u1 u2 u3

cðzÞ ¼ 1; cðzÞ ¼ z; cðzÞ ¼ iz; .. .

cðzÞ ¼ zm ; cðzÞ ¼ izm ;

u2m

Table 1 The relative errors e1, e2 (q0/D) for Example 1

¼ RecðzÞ ¼ 1; ¼ RecðzÞ ¼ x;

¼ RecðzÞ ¼ y; .. . ½m=2 P m2l 2l ¼ RecðzÞ ¼ ð1Þl C 2l y ; mx l¼0

u2mþ1 ¼ RecðzÞ ¼

365

m

6

8

10

12

e1 e2 m

3.6018e4 7.3737e3 14

6.1198e4 5.6106e3 16

6.7119e5 6.3014e4 18

7.6528e5 3.4718e4 20

e1 e2

4.6704e6 1.3914e4

1.3698e5 5.8949e5

1.1041e6 4.8924e5

3.5126e6 1.4895e5

½mþ1=2 P l¼1

m2lþ1 2l1 ð1Þlþ1 C 2l1 y : m x

.. .

.. .

Now let c(z) ¼ 0 and j(z) be polynomials of z to obtain another part of biharmonic functions. u2mþ2 ¼ Re½¯zjðzÞ ¼ x2 þ y2 ;

jðzÞ ¼ z; jðzÞ ¼ z2 ; jðzÞ ¼ iz2 ;

u2mþ3 ¼ Re½¯zjðzÞ ¼ ðx2 þ y2 Þx; u2mþ4 ¼ Re½¯zjðzÞ ¼ ðx2 þ y2 Þy; ...

.. .

jðzÞ ¼ zmþ1 ; jðzÞ ¼ izmþ1 ;

u4mþ1 ¼ Re½¯zjðzÞ ¼ z¯zu2m ¼ ðx2 þ y2 Þu2m ; u4mþ2 ¼ Re½¯zjðzÞ ¼ z¯zu2mþ1 ¼ ðx2 þ y2 Þu2mþ1 : .. .

.. .

Here we denote them by Part II. The two parts together are the Trefftz basis functions. It is easy to prove that they are linearly independent. In fact, according to Refs. [11,12], the set is also complete. Here the functions whose degree is no more than m(mX4) are used. There are 4m2 functions totally, where 2m+1 in Part I and 2m3 in Part II. To determine {bk}k ¼ 1,2,y4m2, the collocation method is used. We set {(x1j, y1j)}j ¼ 1,2ya on G1, {(x2k, y2k)}k ¼ 1,2yb on G2, {(x3l, y3l)}l ¼ 1,2yc on G3, and totally N ¼ a+b+c. The algebraic equations are  0 1     ½ u1 u2    u4m2  ðx ;y Þ2 G B 1 C 1j 1j B C0 1 ðx2k ;y2k Þ2G2 B C b1  B C  B C B C½ u u    u C 1 2 4m2  B CB b 2 C ðx1j ;y1j Þ2G1 CB C B C B CB  B CB .. C B S½ u1 u2    u4m2  CB . C @ A B C ðx2k ;y2k Þ2G2 C B ðx3l ;y3l Þ2G3 B C b4m2 B C  @ A     F½ u1 u2    u4m2  0

 g 1  up ðx

1j ;y1j Þ2G1

ðx3l ;y3l Þ2G3

1

B C ðx2k ;y2k Þ2G2 B C  B C B g 2  Cup ðx ;y Þ2G C B 1 C 1j 1j B C. ¼B  C B g 3  Sup ðx2k ;y Þ2G2 C 2k B C ðx3l ;y3l Þ2G3 B C @ A  g  Fup  4

Fig. 1. The convergence rate of log10(e1) for Example 1.

(2.15)

Fig. 2. The distribution of collocation points for Example 2.

ðx3l ;y3l Þ2G3

Generally, the number of equations 2N should be no less than the number of unknowns 4m2, then the least-square method can be used to solve (2.15).

Table 2 The relative errors at (3/8, 3/8) for Example 2 (m ¼ 12) The number of collocation points

3. Numerical examples

Example 1. The linear bending problem of the S–S–S–S rectangular plate ([0, 1]  [0, 1]) with uniform load q0. All boundaries are clamped and w ¼ 0.005q0/D, v ¼ 0.3.

AB, CD, DA

BC

8 7 6 5

4 7 10 13

e3 (q0/D)

e4 (q0)

e5 (q0)

8.8149e6 4.1343e5 1.2942e6 2.6471e6

7.5566e5 2.7744e5 8.4574e7 1.5110e6

3.5486e4 6.3387e4 9.1649e6 2.7677e6

ARTICLE IN PRESS 366

X. Liu, X. Wu / Engineering Analysis with Boundary Elements 33 (2009) 363–367

Table 3 The relative errors (q0/D) of w in the domain for Example 2 r/y

0.10p

0.15p

0.20p

0.25p

0.30p

0.35p

0.40p

0.4 0.5 0.6 0.7 0.8 0.9 1.0

3.2231e6 4.7396e6 6.6307e6 9.0816e6 1.2815e5 2.0012e5 4.0974e5

3.3667e6 4.4196e7 5.7300e6 7.4418e6 9.8616e6 1.3635e5 2.0376e5

3.3778e6 4.2620e6 5.3361e6 6.6937e6 8.4893e6 1.0980e5 1.4620e5

3.4594e6 4.2936e6 5.2970e6 6.5361e6 8.1109e6 1.0169e5 1.2943e5

6.3651e6 4.5385e6 5.6373e6 7.0293e6 8.8582e6 1.1356e5 1.4917e5

3.8896e6 4.9840e6 6.3774e6 8.2338e6 1.0917e5 1.5209e5 2.3089e5

3.9795e6 5.6447e6 7.6671e6 1.0287e5 1.4307e5 2.2072e5 4.4598e5

According to Navier’s double trigonometric method, the exact solution at (x, y) ¼ (0.5, 0.5) is wmaxE0.00406235266067q0/D [13]. The highest degree of basic functions is denoted as m and the number of collocation points on each boundary is denoted as n. After some numerical experiments, we find that it is of great efficiency when m and n have the relationship n ¼ [m/2]+1 and when m is even, the solution is almost the same to that for the next odd number. The relative errors of the deflection wmax are denoted as e1 and the biggest error on boundaries are denoted as e2, both given in Table 1, while log10(e1) is shown in Fig. 1. The dotted and dashed line are the linear regressions of e1 when m take {4i+2} and {4i+4} (i ¼ 1,2,y), which are marked by stars and triangles, respectively. It is clear that the errors are of exponential convergence. Although e2 is bigger, it is proportional to e1 in the main. Example 2. The linear bending problem of the irregular plate as Fig. 2, with the load q ¼ q0 sin px sin py. AB and CD are simple supported, DA is clamped and BC is free. BC is a circinal arc and its center is at the coordinate origin.

Table 4 The deflection (102q0/D) on the free boundary AB for Example 3 x¼

0

1/8

1/4

3/8

1/2

DQTM m ¼ 12 m ¼ 16 m ¼ 20 Ref. [18] (73 Eqs.) Ref. [13] (20 Eqs.) kahtopobn# Finite element

13.1927 13.0312 12.9741 12.933 12.765 12.11 12.708

13.2532 13.0946 13.0382 12.998 12.819 – 12.788

13.3009 13.1456 13.0891 13.056 12.910 – 12.815

13.3317 13.1774 13.1207 13.091 12.971 – 12.892

13.3425 13.1878 13.1313 13.102 12.991 11.92 12.905

should satisfy

q2 w ¼ 0. qxqy

(3.2)

Here we take the relationship n ¼ [m/2] and calculate deflections on AB (y ¼ 0). The numerical results and some solutions of other methods for comparison are given in Table 4 [13].

In order to check the accuracy of DQTM, we use w¼

q0 4p

4D

4. Conclusions sin px sin py

(3.1)

as the exact solution and the deductive function values and derivatives on boundaries as boundary conditions to calculate the numerical solutions, then compare them with (3.1). Because of different curvatures, collocation points on straight and curve boundaries are set different to show their influences on the numerical solutions. Here we use 28 points when m ¼ 12. Besides w, we also calculate the bending moment Mx and shearing force Qx (both along x direction) at (3/8, 3/8). The relative errors are denoted as e3, e4 and e5 respectively and given in Table 2, corresponding to different collocation methods. It is easy to find that generally straight boundaries need fewer points, while curve ones need more. So the collocation points should be chosen wisely to improve accuracy with lower computation cost. And the method is also very accurate for Mx and Qx. The relative errors of w in the domain are shown in Table 3 in polar coordinates. Here 13 points on BC and 5 points on the other three are used. It is shown that the errors are good enough in most of the domain, except near corners. This is because the exact solution is closed to 0. Example 3. The linear bending problem of the rectangular cantilever plate ([0, 1]  [0, 1]) with uniform load q0, where the boundary CD (y ¼ 1) is clamped and the others are free. Cantilever plates are famous difficult problems in classical theories of elastic plates, but they are very important in engineering. Notice that at A (0, 0) and B (a, 0), the solution

In this paper, a kind of highly accurate and general method based on DQTM for linear bending problems of thin plates is discussed. It is especially effective for the following cases: geometrically irregular domains; complex boundary conditions; the calculation of bending moments and shearing forces. Because the Trefftz basis functions are of good property, the numerical solutions appear some disciplinarians, which are helpful for collocation points.

Acknowledgment The support from the National Natural Science Foundation of China (nos. 10671146 and 50678122) is fully acknowledged. References [1] Bellman R, Casti J. Differential quadrature and long term integration. Math Anal Appl 1971;34:235–8. [2] Bellman R, Kashef B, Vasudevan R. The inverse problem of estimating heart parameters from cardiograms. Math Biosci 1974;19. [3] Bert CW, Malik M. Differential quadrature method in computational mechanics: a review. Appl Mech Rev 1996;49:1–28. [4] Kita Eisuke, Kamiya Norio. Trefftz method: an overview. Adv Eng Software 1995;24:3–12. [5] Kita Eisuke, Ikeda Youichi, Kamiya Norio. Indirect Trefftz method for boundary value problem of Poisson equation. Eng Anal Boundary Elem 2003;27:825–33. [6] Jin WG, Cheung YK, Zienkiewicz OC. Application of the Trefftz method in plane elasticity problems. Int J Numer Methods Eng 1990;30:1147–61. [7] Trefftz E. Ein gegenstu¨ck zum ritzschen verfahren. In: Proceedings of the 2nd international congress of applied mechanics, 1926.

ARTICLE IN PRESS X. Liu, X. Wu / Engineering Analysis with Boundary Elements 33 (2009) 363–367

[8] Wu Xionghua, Ren Yu-e. Differential quadrature method based on the highest derivative and its applications. J Comput Appl Math 2007;205:239–50. [9] Wu Xionghua, Du Huijing, Kong Wenbin. Differential quadrature Trefftz method for Poisson-type problems on irregular domains. Eng Anal Boundary Elem 2008;32:413–23. [10] Yu Dehao. Natural boundary integral method and its applications. Beijing: Science Press; 2002.

367

[11] Xia Yongxu. Weighted residual methods for mechanics in plates and shells. Xi’an: Northwestern Polytechnical University Press; 1994. [12] Liu Chein-Shan. A highly accurate collocation Trefftz method for solving the Laplace equation in the doubly connected domains. Numer Methods Partial Differential Equations 2007;27:179–92. [13] Qu Qingzhang, Zhang Quan, Ji Qiuzhi, Liang Xingfu. Theories of elastic plates. Beijing: China Communications Press; 2000.