Finite Elements in Analysis and Design 41 (2005) 821 – 833 www.elsevier.com/locate/finel
A multivariable wavelet-based finite element method and its application to thick plates Jian-Gang Hana , Wei-Xin Rena,∗ , Yih Huangb a Department of Civil Engineering, Fuzhou University, 523 Gongye Road, Fuzhou, Fujian 350002, PR China b School of Science, Xi’an University of Architecture & Technology, Xi’an, 710055, PR China
Received 2 June 2004; received in revised form 19 October 2004; accepted 7 November 2004 Available online 9 December 2004
Abstract A multivariable wavelet-based finite element method (FEM) is presented to resolve the bending problems of thick plates. The interpolating wavelet functions based on boundary conditions are constructed to represent the generalized field functions of thick plates. The formulation of multivariable wavelet-based FEM is derived by the Hellinger–Reissner generalized variational principle with two kinds of independent variables. The proposed formulation can be solved directly when the stress–strain relations and the differential calculations are not utilized in determining the variables. The applicability of the multivariable wavelet-based FEM is demonstrated by determining the bending solutions of a single thick plate and of an elastic foundation plate. Comparisons with corresponding analytical solutions are also presented. The wavelet-based approach is highly accurate and the wavelet-based finite element has potential to be used as a numerical method in analysis and design. 䉷 2004 Elsevier B.V. All rights reserved. Keywords: Wavelet; Finite element method; Multivariable; Generalized variational principle; Thick plate; Bending
0. Introduction The finite element method (FEM) was developed in the middle of the 20th century. It successfully combined advantages from the conventional energy method and the finite difference method. FEM was based on the energy variational principle and discrete interpolation. It remains an important numerical tool in ∗ Corresponding author. Tel: +86 591 7892454; fax: +86 591 3737442.
E-mail address:
[email protected] (W.-X. Ren) URL: http://bridge.fzu.edu.cn. 0168-874X/$ - see front matter 䉷 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.finel.2004.11.001
822
J.-G. Han et al. / Finite Elements in Analysis and Design 41 (2005) 821 – 833
analysis and design. In the development phase of FEM, the energy variational principle provided a theoretical foundation to develop varieties of finite element formulations [1,2]. Melosh [3] established a finite element displacement formulation using the potential energy principle. Pian [4] created an intercrossed stress finite element formulation utilizing the complementary energy principle. Elias [5] introduced a finite element stress formulation by applying the complementary energy principle. Herrmann [6] created a kind of mixed finite element formulation. Multivariable finite element formulation requires the generalized variational principle with multiple variables. Hellinger [2] first proposed the concept of a generalized variational principle with two kinds of independent variables in 1914. Reissner [10] further clarified this principle in 1950. The generalized variational principle with two kinds of independent variables was therefore called the Hellinger–Reissner variational principle. Shen and He [7,8] as well as Prenter [9] proposed a spline-based FEM by using the Hellinger–Reissner variational principle and applying the formulation to solve engineering problems. Another key issue in developing FEM was the selection of the discrete interpolation function. Wavelet theory is a mathematical tool that was developed recently. As a new developing branch in mathematics, wavelet theory has gained attention in the engineering fields, especially in signal analysis, image handling, pattern recognition, phonetic recognition, quantum physics, earthquake reconnaissance, fluid mechanics, electromagnetic fields, CT imagery, diagnosis and monitoring of machinery defects, fractals, etc. [12–15]. Wavelet analysis based on wavelet transform has extraordinary characteristics that combine the advantages of functional analysis, Fourier transform, spline analysis, harmonic analysis and numerical analysis. Some wavelet transform functions have been used as interpolation functions to develop wavelet-based finite element formulations [16–18]. The objective of this paper is to present a multivariable wavelet-based finite element formulation to resolve the bending problems of thick plates. The scaling functions of interpolating wavelets are selected because the scaling functions in multi-scale wavelet analysis have excellent analysis and computing capabilities. Interpolating wavelets are symmetric and their calculating accuracy is high. The scaling functions of interpolating wavelets are first reconstructed in order to satisfy the boundary conditions of the thick plate. The interpolating wavelet function are then used to represent the generalized field functions of thick plates. The multivariable wavelet-based finite element formulation is derived by using the Hellinger–Reissner generalized variational principle with two kinds of plate variables. One variable is the plate deflection and the other is the generalized stresses of the plate. The proposed multivariable formulation can be solved directly when the stress–strain relations and the differential calculations are not utilized in determining the variables. The applicability of the multivariable wavelet-based FEM is demonstrated by determining the bending solutions of a single thick plate and of an elastic foundation plate. Accuracy and convergence are observed. 1. Hellinger–Reissner generalized variational principle The generalized potential energy principle of a thick plate with two categories of variables is defined under homogeneous boundary conditions [2] as follows: T T T T 1 T −1 2p = − M D W d + Q E W d − WT P d 2 M Db M d − T −1 1 −2 Q C Q d, (1)
J.-G. Han et al. / Finite Elements in Analysis and Design 41 (2005) 821 – 833
823
where the displacement function, W, and the general force functions, M and Q, are two categories of independent variable functions that are unknown. For a thick plate, the eight unknown variables are formulated as follows: W = [w
x
y ]T ,
M = [ Mx
My
Q = [ Qx
Qy ]T ,
where P = [ q D=
0 0
j jx
0 Db =
Mxy ]T ,
0 ]T is a load vector under uniform loading, q, and j j 0 0 1 jx jy j 0 jy , E = −1 0 , Db = Db 1 j j 0 0 0 −1 jy jx
Et 3 , 12(1 − 2 )
C=
0 0
1− 2
,
5Et 12(1 + )
in which Db is the plate bending stiffness, C is the shear stiffness of plate, E is the Young’s modulus and
is the Poisson ratio.
The generalized potential energy variational principle requires that the generalized potential energy function (1) has stationary values for the displacement and general force variations. That is, −1 T T T T −1 2p = M (D W + Db M) d + Q (−E W + C Q) d − QT ET W d T T + M D W d + WPT d = 0. (2)
For a thick plate, the third and fourth terms of Eq. (2) are integrated using the boundary conditions. Taking into consideration that W, M and Q are independent and not equal to zero, the following equations can be obtained −DTW = D−1 b M,
(3)
ETW = C −1 Q,
(4)
DDb DTW − EC ETW + P = 0.
(5)
Eq. (5) is basically the equilibrium equation of a thick plate that is expressed by the general 1displacement components. When a thick plate lies on the Winkler foundation, an additional term, 2 ke w 2 dx dy, is needed in Eq. (1), where ke is the elastic foundation coefficient. 2. Wavelet interpolating functions 2.1. Interpolating wavelets This section describes the wavelet interpolating functions and their principal characteristics. The scaling functions, (x), and wavelet functions, (x), of interpolating wavelets satisfy two scale relations
J.-G. Han et al. / Finite Elements in Analysis and Design 41 (2005) 821 – 833
824
Table 1 Filtering coefficients, p(k) and q(k) k
−4
−3
−2
−1
0
1
2
3
4
p(k) q(k)
0 0
−1/16 0
0 0
9/16 0
1 0
9/16 −1
0 0
−1/16 0
0 0
as follows: M
(x) =
p(k)(2x − k),
(6)
q(k)(2x − k),
(7)
k=−M M
(x) =
k=−M
where p(k) and q(k) are the filtering coefficients and M is a minimum positive integer. The compact domains of scaling and wavelet functions are supp() = [−(M − 1), (M − 1)], supp() = [−(M/2 − 1), (M/2)]. Interpolating wavelets have a number of unique properties that are not possessed by Daubechies-type wavelets that are commonly used [13], namely (1) (2) (3) (4) (5)
Interpolating wavelets have a symmetrical axis. The scaling functions are identical to the physical space representations. Moments and expansion coefficients are easily computed. Interpolating wavelets are more accurate in certain aspects. The basic functions (but not their derivatives) can be computed without solving an eigenproblem.
Specifically, the scaling functions have the following properties. (l) = l
l ∈ Z,
(x) = (−x).
(8) (9)
Due to these properties, determining the integer point values of interpolating wavelets is very convenient. The interpolating wavelets can be obtained by Eq. (6)–(8) without solving an eigenproblem, which simplifies the calculation steps. When M = 4, for instance, the values of the filtering coefficients, p(k) and q(k), are given in Table 1. The scaling and wavelet functions of the interpolating wavelets are plotted in Fig. 1.
J.-G. Han et al. / Finite Elements in Analysis and Design 41 (2005) 821 – 833
825
Fig. 1. Interpolating wavelets: (a) Scaling function (x); (b) Wavelet function (x).
For an arbitrary integer, j ∈ Z, the scaling functions, j,l (x) = 2j/2 (2j x − l), where −∞ < l < ∞ and l ∈ Z, are the orthogonal basis in the subspace, L2 (R). When Vj is f ∈ L2 (R) | f (x) =
∞ l=−∞
cj,l j,l (x), −∞ < x < ∞, {cj,l }∞ l=−∞ ⊂ R
(10)
for any f (x) ∈ Vj , f (x) can be denoted as follows due to the compact support of the scaling functions: f (x) =
∞
cj,l j,l (x)
− ∞ < x < ∞.
(11)
l=−∞
2.2. Higher derivatives of interpolating wavelets Since the wavelets are not explicit expressions, the higher derivative values of interpolating wavelets can be obtained only by the numerical method. The dth-order derivative values of interpolating wavelets (6) can be obtained from M
(d) (x) = 2d
p(k)(d) (2x − k).
(12)
k=−M
For supp(x) = [−(M − 1), (M − 1)], it is found that supp(d) (x) ⊆ [−(M − 1), (M − 1)]. By extending Eq. (12) into the integer points of interval [−(M − 1), (M − 1)], the (2M − 1)th-order system equations can be determined. The equations are rearranged in the matrix form A(d) (x) = 0,
(13)
where (d) (x) = [(d) (−M + 1), (d) (−M + 2), . . . , (d) (M − 1)]T and A is an matrix of the base coefficients. In order for Eq. (13) to have sole solutions, additional restraints have to be applied. The following equation can be obtained from the properties of the interpolating wavelets [13]. xd =
∞ l=−∞
l d (x − l),
d = 0, 1, . . . , M − 1.
(14)
J.-G. Han et al. / Finite Elements in Analysis and Design 41 (2005) 821 – 833
826
The dth-order derivative of Eq. (14) gives ∞
d! =
l d (d) (x − l)
− ∞ < x < ∞.
l=−∞
By setting x =0 and using the symmetrical relationship (d) (−l)=(−1)d (d) (l), the additional restraint condition can be obtained ∞
l d (d) (l) = (−1)d d!.
(15)
l=−∞
2.3. Improved interpolating functions In the current multivariable wavelet-based finite element method, the scaling functions of dual interpolating wavelets are used to represent the field functions of thick plates. Those functions are smooth and can satisfy continuity requirements. In practical applications, uniform meshes are used in the wavelet-based FEM for the interval [x0 , xN ]. If a total of N nodes is selected, the coordinate xi of specific node i is xi = x0 + ih,
h = 2−j =
xN − x 0 , N
where h is the mesh distance. When M = 4, any function S(x) with an interval [x0 , xN ] can be expressed as follows according to Eq. (11): S(x) =
N +2 i=−2
x − x0 −i , Ci h
(16)
where Ci is the wavelet coefficient. By setting x0 = 0 for convenience, the boundary conditions of a thick plate can be presented as • Fixed side: w = 0, n = s = 0, • Simply supported side: w = 0, Mn = 0, s = 0, • Free side: Mns = 0, Qn = 0. The boundary conditions of a thick plate may have elastic supported sides, side forces and bending moments. Specific boundary conditions can be included by setting the specific constraints on the sides of the plate. The boundary conditions w = 0, n = s = 0 at the fixed sides and w = 0, s = 0 at the simply supported sides are all forced boundary conditions. The other boundary conditions of a thick plate are the natural boundary conditions. Once the generalized potential energy (1) is minimized, i.e. if the continuum retains its equilibrium status, the natural boundary conditions will be automatically satisfied. In other words, the basis functions are not required to adapt to the natural boundary conditions. However, the forced boundary conditions have to be satisfied after formulating the system equations. To deal with the plate boundary conditions easily, the five basis functions near the supports are linearly combined by using a trial and error process to determine the scaling basis functions of the interpolating wavelets. They are
J.-G. Han et al. / Finite Elements in Analysis and Design 41 (2005) 821 – 833
827
denoted as follows: S(x) =
N +2
ri i (x) = (x)r.
(17)
i=−2
Eq. (17) is called the scaling basis functions of modified interpolating wavelets where r is a column matrix that is composed of the scaling function coefficients at each node: r = [r−2 , r−1 , r0 , . . . , ri , . . . , rN , rN +1 , rN +2 ]T
(18)
and the scaling basis functions are (x) = [−2 (x), −1 (x), 0 (x), . . . , i (x), . . . , N (x), N +1 (x), N +2 (x)].
(19)
When N 6, the first five functions of the modified scaling basis functions, namely, −2 (x), −1 (x), 0 (x), 1 (x) and 2 (x) are as follows:
3 x +2 + +1 , h h h h 24 h
x x x 3 x −1 + +1 , −2 + −1 . h h h 24 h The middle functions of the modified scaling basis functions can be expressed as x
i (x) = − i , 2 < i < N − 2. h
x
+
x
+1 ,
x
+1 ,
x
The last five functions of the modified scaling basis functions, namely, N −2 (x), N −1 (x), N (x),
N +1 (x) and N +2 (x) are as follows:
x x 3 x −N +2 + −N +1 , −N +1 + −N −1 , h 24 h h h x x
x x x 3 −N −2 + −N −1 , −N −1 , −N + −N −1 . h 24 h h h h The modified scaling functions have the following characteristics: at the left end, x = 0, of the interval
x
k (0) = 0, k (0) = 0,
(k = −2); (k = −2, −1);
at the right end, x = xN , of the interval k (xN ) = 0, k (xN ) = 0,
(k = N + 2); (k = N + 2, N + 1).
These characteristics of the basis functions are useful in dealing with boundary conditions. The boundary conditions in the form of plate displacements and general forces that appear in the multivariable wavelet-based finite element formulations can be treated as easily as in the conventional finite element formulations.
J.-G. Han et al. / Finite Elements in Analysis and Design 41 (2005) 821 – 833
828
3. Multivariable wavelet-based finite element formulations The multivariable field functions of a thick plate are constructed by using the dual scaling basis functions of interpolating wavelets in product form M = 1W1 ,
(20)
Q = 2W2 ,
(21)
W = 3W3 ,
(22)
where W1 = [A1 B1 C1 ]T ,
W2 = [A2 B2 ]T ,
W3 = [A3 B3 C3 ]T
(23)
in which A1 , B1 , C1 , A2 , B2 , A3 , B3 and C3 are all the constant coefficient row matrices and W1 , W2 and W3 are the column matrices that are all unknown quantities to be solved. For example A1 = [A1−2 , A1−1 , A10 , . . . , A1N , A1N +1 , A1N +2 ], 1 1 1 1 1 1 A1i = [a−2i , a−1i , a0i , . . . , aN i , aN +1i , aN +2i ],
The matrices 1 , 2 and follows: (x) ⊗ (y) 1 = 0 0 (x) ⊗ (y) 2 = 0 (x) ⊗ (y) 3 = 0 0
(24) (i = −2, −1, 0, . . . , N, N + 1, N + 2).
(25)
3 are related to the scaling basis functions and their expressions are as
0 (x) ⊗ (y) 0
0 , 0 (x) ⊗ (y)
0 , (x) ⊗ (y) 0 (x) ⊗ (y) 0
0 , 0 (x) ⊗ (y)
(26)
(27)
(28)
where ⊗ denotes the Kronecker product of two matrices, i.e., A ⊗ B = aim bj l .
Eqs. (20)–(22) are substituted into the generalized potential energy Eq. (1) and then they are minimized by the Hellinger–Reissner generalized variational principle with two kinds of plate variables. j2p jWi
= 0,
(i = 1, 2, 3).
The multivariable wavelet-based finite element can be formulated as follows: F H K 0 = , T H
Ke
W3
R
(29)
(30)
J.-G. Han et al. / Finite Elements in Analysis and Design 41 (2005) 821 – 833
829
where F
1 A00 ⊗ A00 y Db (1−2 ) x 00 00 Db (1−2 ) Ax ⊗ Ay
=
A00 ⊗ A00 y Db (1−2 ) x −1 00 ⊗ A00 A y Db (1−2 ) x
0
0
0
0
0
0
0
−2 A00 ⊗ A00 y Db (1−2 ) x
0
0
0
0
0
0
0
0 0 H= 0 A01 ⊗ A00 x y 01 A00 ⊗ A x y
Ke = ke
0
−1 A00 ⊗ A00 x y C
0
0
00 −A01 x ⊗ Ay 0 00 −Ax ⊗ A01 y 00 −A00 ⊗ A x y 0
00 A00 x ⊗ Ay
0 0
0 0 0
0 01 −A00 x ⊗ Ay 00 −A01 x ⊗ Ay , 0 00 −A00 x ⊗ Ay
,
0 −1 A00 ⊗ A00 x y C
0 0 , 0
K = [A1 B1 C1 A2 B2 ]T , R = [PT1 0 0]T
in which Ke is the stiffness matrix of the elastic foundation and L L 00 T 01 Ax = [(x)] [(x)] dx Ax = [(x)]T [ (x)] dx 0
(x = x, y).
(31)
0
Eq. (31) is a square matrix where x can be alternatively substituted by y. The elements of matrices A00 x and A01 can be calculated as follows: x 00 A00 x = Ay = −h/7456566
×
−2582017.11
192565.4 16349.7 −19457.6 −1827.2 −189.6 sym
01 A01 x = Ay = −1/671096
332996.8 −31784.5 −3110.7 3658321 × −16784.4 4078.3 .. .
−890819.6 154050.8 130474.4 25.4 14832.8 24.0 −5673347.1 −1766092.8 −6321676.5
37999 2686.4 −354589 312.9 −145.0 35264.8 149.2 −0.4 3530.1 −35314.3 −3507.6 −48.4 −1106.8 −138.1 464382.4 −256.9 −33.2 −70842.2 .. .. .. . . .
−20235.7 563 562 1 71.3 0.1 300426 −20796.7 −983811 297264.3 −5973211.1 −1021294 .. .
16570.9 −4106.9 259.1 1107.7 258.1 1 136.2 33.3 0.1 −464093.6 71041.6 −4332.7 18.3 −456208.3 70237.8 456517.2 −16.7 −465056.2 .. .. .. . . .
1 562 1 −20796.7 562 1 299864 −20796.7 562 .. .. .. . . .
1 258.1 1 −4332.7 258.1 1 70783.5 −4332.7 258.1 .. .. .. . . .
. 1 .. .
; 1 .. .
J.-G. Han et al. / Finite Elements in Analysis and Design 41 (2005) 821 – 833
830
The load column matrix is L L P1 = q(x, y)(x) ⊗ (y) dx dy 0
(32)
0
when the load is uniformly distributed: L q (x) dx 0 885 357 17 = qh 1939 − 8192 − 4081 1 when the load is linearly distributed: L 291 x (x) dx = h 896 929 − 16384 128 65793
17N − 4081
··· 1
1
128 − 65793
0
− 896 929
271 240
896 929
291 16384
271 240
1501 707
− 357N 8192
17 − 4081
··· N − 3 885N 896 − 929 , 1939
3
when the load is concentrated: L L (x − )(y − )(x) ⊗ (y) dx dy = () ⊗ ( ), 0
1
357 − 8192
271N 240
885 1939
− 1501 707
,
N
(33)
0
where and are the coordinates at the concentrated load point. 4. Numerical examples In this section, the applicability of the multivariable wavelet-based finite element formulations is verified by solving the bending problems of a rectangular thick plate and a thick plate on an elastic foundation. The numerical calculations are implemented in Matlab programming [19] in terms of the above multivariable wavelet finite element formulations. 4.1. Bending of rectangular thick plates Rectangular thick plates simply supported and fixed on four sides under uniform loading, q, are studied by the proposed multivariable wavelet-based finite element method. The thickness and edge length of the plate are represented by t and l, respectively. The Poisson ratio of the plate is = 0.3. Table 2 shows the comparison of the mid point deflections between the current solutions and the exact solutions for the different thickness-span ratios where the 8 × 8 mesh division is used. The results obtained from the multivariable wavelet-based finite element formulations agree well with the exact solutions for both the simply supported and fixed supported boundary conditions of the rectangular thick plates. Table 3 demonstrates the comparison of the deflections and bending moments at the mid point of the rectangular thick plates under uniform loading with different mesh divisions, (t/ l = 0.3). It can be observed that both deflections and bending moments rapidly converge to the exact solutions with an increase in the mesh division. The proposed multivariable wavelet-based FEM is therefore further verified.
J.-G. Han et al. / Finite Elements in Analysis and Design 41 (2005) 821 – 833
831
Table 2 Mid point deflection of rectangular thick plates under uniform loading, (ql 4 /100Db ) Thickness span ratio t/ l
0.1 0.2 0.3 0.35
Simply supported on four sides
Fixed on four sides
Wavelet FEM solution
Exact solution
Wavelet FEM solution
Exact solution
0.4269 0.4887 0.5952 0.6631
0.4273 0.4906 0.5956 0.6641
0.1493 0.2166 0.3224 0.3948
0.1499 0.2167 0.3227 0.3951
Table 3 The deflection and bending moment of the mid point of rectangular thick plates under uniform loading with different mesh divisions (t/ l = 0.3) Mesh of plate
6×6 8×8 10 × 10 Exact solution
Deflection (ql 4 /100Db )
Bending moment (ql 2 /10)
Simply supported square plate
Fixed square plate
Simply supported square plate (mid point)
Fixed square (midpoint)
0.5947 0.5952 0.5953 0.5956
0.3218 0.3224 0.3226 0.3227
0.4932 0.4875 0.4810 0.4789
−0.4087 −0.4126 −0.4182 −0.426
Table 4 The deflection and stresses of a thick plate on an elastic foundation under uniform loading Coordinate (m)
(10−8 m)
(x, y)
Series solution [11]
Wavelet FEM solution
Series solution [11]
Wavelet FEM solution
(0,0) (0.1,0) (0.2,0) (0.3,0) (0.4,0) (0.5,0)
0.576 0.568 0.546 0.530 0.501 0.477
0.566 0.560 0.542 0.524 0.500 0.481
0.928 0.866 0.668 0.352 0.154 0.091
0.825 0.763 0.562 0.260 0.054 0.0
x (N/m2 )
4.2. Bending of a square thick plate on an elastic foundation A square thick plate on an elastic foundation is shown Fig. 2. The length of the plate is a = 1 m. The range of the applied uniform loading, q = 1 N/m2 , is a central half square with the length of b = 0.5 m. The material properties of the plate are = 0.35 and E = 300 MN/m2 . The elastic foundation coefficient of the Winkler foundation is ke =50 MN/m2 . The thickness of the plate is t =0.4. The bending problem of the square thick plate on an elastic foundation has been solved by the current multivariable wavelet-based FEM. Table 4 gives the comparisons of the deflections and stresses obtained from the current solutions
832
J.-G. Han et al. / Finite Elements in Analysis and Design 41 (2005) 821 – 833
Fig. 2. A thick plate on an elastic foundation.
with those obtained from the Fourier series solutions. An 8 × 8 mesh division is used in the multivariable wavelet-based FEM. The boundary conditions of the elastic foundation plate are free on four sides. The bending stresses are calculated using the relationship between stresses and bending moments of the thick plate. The proposed multivariable wavelet-based FEM achieves excellent. 5. Conclusions A multivariable wavelet-based finite element formulation has been developed to solve the bending problems of thick plates based on the Hellinger–Reissner generalized variational principle with two kinds of independent variables. The dual wavelet bases are constructed to represent the general field functions of a thick plate. The proposed multivariable formulation can be solved directly when the stress–strain relations and the differential calculations are not implemented in determining the variables. As the proposed formulation is mixed where the displacements and general forces are independent field functions, the method has high precision and good convergent characteristics over other displacement based approaches. The applicability of the multivariable wavelet-based FEM is demonstrated from the bending solutions of a thick plate and an elastic foundation thick plate. Favorable accuracy and convergence of the deflections, bending moments and stresses of the thick plates with various supports conditions have been observed. The wavelet-based finite element has potential as a numerical method in solving real-world problems. Acknowledgements Support from the Natural Science Foundation of China (NSFC), under grant number 59978038, is greatly acknowledged.
J.-G. Han et al. / Finite Elements in Analysis and Design 41 (2005) 821 – 833
833
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19]
K. Washizu, Variational Method in Elasticity and Plasticity, second ed, Pergamon Press, Oxford, 1975. Hu Haichang, Variational Principle and Application of Elastic Mechanics, Science Press, Beijing, China, 1981. R.J. Melosh, Basis for derivation of matrices for the direct stiffness method, AIAA J. 1 (1963) 1631–1637. T.H.H. Pian, Derivation of element stiffness matrices by assumed stress distributions, AIAA J. 2 (1964) 1333–1336. Z.M. Elias, Theory and Methods of Structural Analysis, Wiley, New York, 1986. Herrmann, Finite element bending analysis for plats, J. Eng. Mech. Div. ASCE 93 (1967) 13–26. Shen Pengcheng, He Peixiang, Bending analysis and spherical shells by multivariable spline element method based generalized variational principle, Comput. Struct. 55 (1995) 151–157. Shen Pengcheng, He Peixiang, Multivariable spline FEM method, Acta Mech. Solida Sinica 15 (1994) 234–243. P.M. Prenter, Splines and Variational Methods, Wiley, New York, 1975. E. Reissner, On a variational theorem in elasticity, J. Math. Phys. (1950) 90–95. D.J. Henwood, Fourier series for a rectangular thick plate with free edges on an elastic foundation, Int. J. Numer. Methods Eng. 18 (1982) 1801–1820. I. Daubechies, Ten lectures on wavelets. CBMS–NSF Regional Conference Series in Applied Mathematics, Department of Mathematics, University of Lowell, MA, Society for Industrial and Applied Mathematics, Philadelphia, 1992. V.A. Barker, Some computational aspects of wavelets, Informatics and Mathematical Modeling, Technical University of Denmark, Denmark, 2001. Du Shoujun, Zhao Guojing, Multiresolution spline plane element, Eng. Mech. 16 (1999) 33–40 (in Chinese). C.S. Burrus, R.A. Gopinath, H. Guo, Introduction to Wavelets and Wavelet Transforms-A primer, Prentice-Hall Inc., Upper Saddle River, NJ, 1998. Ma Junxing, Xue Jijun, Yang Shengjun, He zhengjia, A study of the construction and application of a Daubechies waveletbased beam element, Finite Elements Anal. Design 39 (2003) 956–975. J. Ko, A.J. Kurdila, M.S. Pilant, A class of finite element methods based on orthonormal compactly supported wavelets, Comput. Mech. 16 (1995) 235–244. J. Ko, A.J. Kurdila, M.S. Pilant, Triangular wavelet-based finite elements via multivalued scaling equations, Comput. Methods Appl. Mech. Eng. 146 (1997) 1–17. The Mathworks, Inc. Matlab Release 5, Natick, MA, USA, 1996.