- Email: [email protected]
Use of Fourier shape functions in the scaled boundary method
Engineering Analysis with Boundary Elements 41 (2014) 152–159
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
Use of Fourier shape functions in the scaled boundary method Yiqian He a, Haitian Yang a, Andrew J. Deeks b,n a State Key Lab of Structural Analysis for Industrial Equipment, Department of Engineering Mechanics, Dalian University of Technology, Dalian 116024, PR China b President's Office, University College Dublin, Belfield Campus, Dublin 4, Ireland
art ic l e i nf o
a b s t r a c t
Article history: Received 21 November 2013 Accepted 20 January 2014 Available online 14 February 2014
The scaled boundary finite element method (SBFEM) is a semi-analytical method, whose versatility, accuracy and efficiency are not only equal to, but potentially better than the finite element method and the boundary element method for certain problems. This paper investigates the possibility of using Fourier shape functions in the SBFEM to form the approximation in the circumferential direction. The shape functions effectively form a Fourier series expansion in the circumferential direction, and are augmented by additional linear shape functions. The proposed method is evaluated by solving three elastostatic and steady-state heat transfer problems. The accuracy and convergence of the proposed method is demonstrated, and the performance is found to be better than using polynomial elements or using an element-free Galerkin approximation for the circumferential approximation in the scaled boundary method. & 2014 Elsevier Ltd. All rights reserved.
Keywords: Scaled boundary method Fourier shape functions Computational accuracy Stress singularities Unbounded domains
1. Introduction The scaled boundary method (SBM) is a semi-analytical method developed relatively recently by Wolf and Song [1]. The method introduces a normalised radial coordinate system based on a scaling centre and a defining curve (usually taken as the boundary). The governing deferential equations are weakened in the circumferential direction and then solved analytically in the normalised radial direction. The SBM has combined the advantages of the Finite Element Method (FEM) and the Boundary Element Method (BEM), and no fundamental solution is required like BEM. In addition, the SBM has been shown to be more efficient than the FEM for problems involving unbounded domains and for problems involving stress singularities or discontinuities [2]. Effective applications of this method have been demonstrated in various problem domains, including fracture problems [3–6] and foundation problems [7–10]. In the scaled boundary method, the discretisation approach used in the circumferential direction has significant influence on the accuracy of the resulting solutions [11]. The most commonly used method for performing this circumferential discretisation is the finite element approach, leading to the method called the scaled boundary finite element method (SBFEM). Vu and Deeks [12–14] investigated the use of higher-order polynomial shape functions in the SBFEM, and demonstrated the SBFEM converged significantly faster under p-refinement than under h-refinement.
n
Corresponding author. E-mail address: [email protected] (A.J. Deeks).
0955-7997/$ - see front matter & 2014 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.enganabound.2014.01.012
The development of meshless methods provided another approach to building circumferential approximations for the scaled boundary method. Deeks and Augarde [11] developed a Meshless Local Petrov–Galerkin method scaled boundary method (MLPG-SBM) and He et al. [15] developed an Element-free Galerkin scaled boundary method (EFG-SBM). This work showed that these two meshless scaled boundary methods gave a higher level of accuracy and rate of convergence than the conventional SBFEM using linear or quadratic elements, with the EFG-SBM performing slightly better than the MLPG-SBM. In this paper, the possibility of using shape functions based on the terms of a Fourier series for the circumferential approximation of the SBFEM is investigated. Fourier interpolations containing trigonometric functions have been applied to both the finite element method (FEM) and the boundary element method (BEM). For example, Guan et al. [16] developed a Fourier series based FEM into for the analysis of tube hydroforming, and showed that this Fourier shape function reduced the number of degrees of freedom required. Khaji et al. [17,18] applied Fourier radial basis functions into the BEM, and showed that of the resulting BEM is much more accurate than the BEM using classic Lagrange shape functions. Although the advantages of Fourier based FEM and BEM have been illustrated in previous work, to date there has been no work reported on the use of Fourier shape functions in the SBFEM. A new Fourier-based scaled boundary method (F-SBM) is presented in this paper. A set of shape functions based on Fourier series expansion is derived, and augmented with linear shape functions. The new shape functions provide good approximation to both trigonometric and polynomial functions in the circumferential direction of the
Y. He et al. / Engineering Analysis with Boundary Elements 41 (2014) 152–159
scaled boundary system. In three numerical examples, the F-SBM is used to solve two-dimensional elastostatic and steady-state heat transfer problems. The accuracy and convergence of F-SBM is compared with the conventional SBFEM using both linear and quadratic elements and with the EFG-SBM. Superior performance in terms of both accuracy and convergence is demonstrated. This paper is organised as follows: The basic equations of scaled boundary method are given in Section 2. Section 3 introduces the Fourier shape functions for use in the scaled boundary method. Three example problems are presented in Section 4 to verify the effectiveness of proposed method, and the paper draws conclusions at the end.
The scaled boundary method introduces a normalised radial coordinate system by scaling a defining curve (usually the domain boundary or a part of the domain boundary) relative to a scaling centre ðx0 ; y0 Þ selected within the domain or at the intersection of two straight sections of the boundary (Fig. 1). The normalised radial coordinate ξ runs from the scaling centre towards the defining curve, and has values of zero at the scaling centre and unity at the defining curve. The other circumferential coordinate s specifies a distance around the defining curve from an origin on the curve. The scaled boundary and Cartesian coordinate systems are related by the scaling equations x ¼ x0 þ ξxs ðsÞ
ð1Þ
y ¼ y0 þ ξys ðsÞ
ð2Þ
Displacement and stress components are retained in the original Cartesian coordinate directions, while position is specified in terms of the scaled boundary coordinates. An approximate solution is sought in the form n
fuh ðξ; sÞg ¼ ∑ ½Ni ðsÞuhi ðξÞ ¼ ½NðsÞfuh ðξÞg
ð3Þ
i¼1
This represents a discretisation of the part of the boundary located at ξ ¼ 1 with the shape function [N(s)]. The unknown vector fuh ðξÞg is a set of n functions analytical in ξ. The shape functions apply for all lines with a constant ξ. (If the scaling centre lies on the boundary, as in Fig. 1, the straight portions of the boundary adjacent to the scaling centre and representing radial lines are not discretised, and in the solution process an analytical solution is found along these lines.) Mapping the linear operator to the scaled boundary coordinate system using standard methods ½L ¼ ½L1
∂ ∂ ∂ 1 2 ∂ 1 þ ½L2 ¼ ½b ðsÞ þ ½b ðsÞ ∂x ∂y ∂ξ ξ ∂s 1
The stresses are obtained by multiplying the strains (obtained form the displacement field using the linear operator) by the elasticity matrix ½D in the form
ð4Þ
2
where ½b ðsÞ and ½b ðsÞ are dependent only on the boundary definition.
1 ξ
sðξ; sÞ ¼ ½D εðξ; sÞ ¼ ½D½B1 ðsÞfuðξÞg;ξ þ ½D½B2 ðsÞ uh ðξÞ
ð5Þ
where 1
ð6Þ
2
ð7Þ
½B1 ðsÞ ¼ ½b ðsÞ½NðsÞ ½B2 ðsÞ ¼ ½b ðsÞ½NðsÞ;s Z V
2. The scaled boundary method
153
In this case the virtual work statement becomes Z fδεðξ; sÞgT fsh ðξ; sÞgdV fδuðsÞgT ftðsÞgds ¼ 0
ð8Þ
S
where the first term represents the internal work and the second term the external work, and ftðsÞg is the external force vector. The virtual strain field is of the form (analogous to Eq. (5)) 1 δεðξ; sÞ ¼ ½B1 ðsÞfδuðξÞg;ξ þ ½B2 ðsÞ δuðξÞ ξ
ð9Þ
where fδuðξÞg is virtual displacement and dV ¼ jJjξdξds ð10Þ where J is the Jacobian at the boundary ðξ ¼ 1Þ. Substituting Eqs. (5), (9) and (10), integrating the area integrals containing fδuðξÞg;ξ with respect to ξ using Green's Theorem, and introducing the coefficient matrices Z ð11Þ ½E0 ¼ ½B1 ðsÞT ½D½B1 ðsÞjJjds S
Z ½E1 ¼
½B2 ðsÞT ½D½B1 ðsÞjJjds
ð12Þ
½B2 ðsÞT ½D½B2 ðsÞjJjds
ð13Þ
S
Z ½E2 ¼ S
the virtual work equation may be expressed as Z n o fδεðξ; sÞgT fsh ðξ; sÞgdV ¼ fδugT ½E0 fuh g;ξ þ ½E1 fuh g V
Z 0
1
1 fδuðξÞgT ½E0 ξfuh ðξÞg;ξξ þ ½½E0 þ ½E1 T ½E1 fuh ðξÞg;ξ ½E2 uh ðξÞ dξ ξ
ð14Þ On substitution of Eq. (3), the external virtual work term becomes Z Z fδuðsÞgT ftðsÞgds ¼ fδugT fNðsÞgftðsÞgds ¼ fδugT fPg ð15Þ S
S
Thus the complete virtual work equation becomes fδugT f½E0 fuh g;ξ þ½E1 fuh gg fδugT fPg Z 1 T n 0 δuðξÞ ½E ξfuh ðξÞg;ξξ þ ½½E0 þ ½E1 T 0 1 ½E1 fuh ðξÞg;ξ ½E2 fuh ðξÞg dξ ¼ 0g ξ
ð16Þ
In order for Eq. (16) to be satisfied for all fδuðξÞg (implying that equilibrium it closely satisfied in the radial direction and in the approximate sense in the circumferential direction), both of the following conditions must be satisfied. fPg ¼ ½E0 fuh g;ξ þ ½E1 T fuh g ½E0 ξ
2
ð17Þ
uh ðξÞ ;ξξ þ ½½E0 þ ½E1 T ½E1 ξ uh ðξÞ ;ξ ½E2 uh ðξÞ ¼ f0g ð18Þ
Fig. 1. Bounded domain with side faces showing scaled boundary coordinate system.
By inspection, solutions to the homogeneous set of Euler– Cauchy differential equations represented by Eq. (18) must be of
154
Y. He et al. / Engineering Analysis with Boundary Elements 41 (2014) 152–159
the form fϕ1 g þ c2 ξ
λ2
fϕ2 g þ :::
ð19Þ
where the exponents λi and corresponding vectors fϕi g may be interpreted as independent modes of deformation which closely satisfy internal equilibrium in the ξ direction. The integration constants ci represent the contribution of each mode to the solution, and are dependent on the boundary conditions. The displacements for each mode take the form (omitting the subscript) fuðξÞg ¼ ξ λ fϕg
trigonometric functions as [17] 1 nπ nπ f ðrÞ ¼ a0 þ ∑ an cos r þ bn sin r Lmax Lmax n¼1
ð30Þ
where a0, an, bn and Lmax represent the Fourier series parameters. Thus on the boundary at ξ ¼ 1, the displacement could be approximated as m jπ jπ s þ bj sin s ð31Þ uh ðsÞ ¼ a0 þ ∑ aj cos L L j¼1
ð20Þ
The vector fϕg can be identified as the modal displacements at the boundary nodes, while λ can be identified as a modal scaling factor for the ‘radial’ direction. Substituting this solution into Eq. (18) yields the quadratic eigenproblem ½λ2 ½E0 λ½½E1 T ½E1 ½E2 fϕg ¼ f0g
ð21Þ
The equivalent nodal forces required at the boundary to equilibrate each displacement mode are obtained by substituting Eq. (20) into Eq. (17) (which is evaluated at ξ ¼ 1) as fqg ¼ ½½E1 T λ½E0 fϕg
ð22Þ
The Eqs. (21)–(22) can be expressed in linear form " #( ) ( ) ϕ ϕ ½E0 1 ½E0 1 ½E1 T ¼ λ q q ½E1 ½E0 1 ½E1 T ½E2 ½E1 ½E0 1
ð23Þ
Solution of this standard eigenproblem yields 2n modes. The eigenvectors contain the modal displacements and the equivalent modal node forces. For a bounded domain only the modes with non-positive real components of λ lead to finite displacements or the scaling centre. This subset of n modal displacements is designated by ½Φ1 , where the vectors in the set from the columns of the matrix. The subset of modal force vectors corresponding to the n modes in ½Φ1 is denoted as ½Q 1 . For any set of boundary node displacements ½uh , the integration constants required to satisfy Eq. (18) on the boundary ðξ ¼ 1Þ are fcg ¼ ½Φ1 1 fuh g
ð24Þ
The equivalent nodal forces required to cause these displacements are fPg ¼ ½Q 1 fcg ¼ ½Q 1 ½Φ1
1
fuh g
where s is the circumferential coordinate in scaled boundary element, L is the length of the boundary and m represents the order of Fourier series. However, for many practical problems, the functions describing the defining curve in the scaled boundary method, xs(s) and ys(s), may be only C0 continuous at certain points along the boundary, rather than fully continuous. These points are termed vertices, in line with conventional usage of the term. Although the derivatives of a solution field with respect to the Cartesian coordinate system may be continuous, derivatives of the Cartesian components of this displacement field with respect to the circumferential coordinate s will not be continuous at the vertices. This is illustrated in Fig. 2. Consequently it is not advantageous to use a continuous Fourier series representation across the vertices. Between the vertices, the derivatives of the solution field are expected to be continuous. The defining curve can thus be defined to consist of edges, over which the solution is expected to be smooth, connected together at vertices, at which the components of the solution field are expected to be C0 continuous with respect to s, but at which the derivatives of these components with respect to s are discontinuous. To preserve the C0 continuity, nodes are introduced at the vertices. Along each
1 0.8 u(x,y) = x.y
fuh ðξÞg ¼ c1 ξ
λ1
0.6 0.4 0.2 0 1
ð25Þ
1
The stiffness matrix of the domain is therefore ½K ¼ ½Q 1 ½Φ1 1
0.5
ð26Þ
0.5
y
0
and the equilibrium requirement is reduced to ½Kfuh g fPg ¼ f0g
ð27Þ
1
The integration constants are then obtained using Eq. (24) and the displacement field is recovered by as n
fuh ðξ; sÞg ¼ ½NðsÞ ∑ ci ξ i¼1
λi
fϕi g
n
0.8
ð28Þ 0.6 u(s)
The stress field is obtained as fsh ðξ; sÞg ¼ ½D ∑ ci ξ λi 1 ½ λi ½B1 ðsÞ þ ½B2 ðsÞfϕi g
s
x
0
ð29Þ
0.4
i¼1
0.2
3. A Fourier shape function
0 0
This paper employs shape functions obtained from the wellknown Fourier series. Based on the theory of the Fourier series, any continuous function f ðrÞ maybe represented by a series of
0.5
1
1.5
2 s
2.5
3
3.5
4
Fig. 2. (a) Smooth function in Cartesian coordinates. (b) Same function plotted against boundary coordinate s.
Y. He et al. / Engineering Analysis with Boundary Elements 41 (2014) 152–159
edge a single element with the shape functions defined using the Fourier approach is introduced. To preserve C0 continuity between the edges or elements, linear polynomial functions terms are added into the standard Fourier approximation as m Ls s jπ jπ þ β1 þ ∑ aj cos s þbj sin s ð32Þ uh ðsÞ ¼ α1 L L j¼1 L L where α1 and β1 represent the values of the function at the end nodes of the element. While it is possible to use the Fourier parameters as the unknown boundary parameters when solving the scaled boundary finite element equations, here the parameters in the Fourier expansion above are transferred to nodal values at equally spaced nodes along each element for ease of applying essential boundary conditions and enforcing C0 continuity between elements. If ð2m þ 2Þ nodes are used, the nodal values vector fug can be related to the parameters in Eq. (32) by ^ fug ¼ ½Tfug
ð33Þ
^ ¼ where fug
n
α1
a1
⋯
am
⋯
b1
bm
β1
oT
, and [T] is a
transfer matrix assembled as T ij ¼ ψ j ðSi Þ
i; j ¼ 1; 2 … 2m þ 2
ð34Þ
where Si is the circumferential coordinate of the ith node, and the component functions of the Fourier expansion are 8 Ls i¼1 > > > L
> > ði 1Þπ > 2rirmþ1 < cos L s
ð35Þ ψ i ðsÞ ¼ > > sin ði 1L mÞπ s m þ 2 ri r 2m þ 1 > > > > :s i ¼ 2m þ 2
4. Performance of the method Here three examples are provided to demonstrate the effectiveness of the proposed method, two elastostatic examples and one steady-state heat transfer example. (The solution procedure of scaled boundary method for steady-state heat transfer problems is similar to that for elastostatic problems provided in the previous sections, and the detailed derivations can be found in references [19,20].)
4.1. Example 1 – a plate of infinite extent subjected to a uniaxial stress The first example is a plate of infinite extent containing a hole of unit radius and subjected to a uniaxial stress of unity. This example demonstrates the ability of the method to model unbounded problems efficiently and accurately. Plane strain conditions are assumed. Considering the symmetry of the problem, one quarter of the plate is represented (Fig. 4(a)). The Young's modulus is take as E ¼ 1000 and Poisson's ratio v ¼ 0:3. The boundary of the hole is considered as a single edge with the length l (Fig. 4(b)). Uniform nodes are introduced along this edge, the nodes spacing is specified as ds. The centre of the hole is used as the scaling centre. The exact solution for the stresses for this problem [21] is 1 3 31 cos 2θ þ cos 4θ þ 4 cos 4θ sxx ¼ 1 2 ð40Þ 2r r 2
^ ¼ ½T 1 fug fug
ð36Þ
Thus the approximation for displacement can be rewritten as uh ðsÞ ¼ fψg½T 1 fug
ð37Þ
The shape functions relating to the nodal displacements are hence fφg ¼ fψg½T 1
ð38Þ
and the shape function matrix for the scaled boundary method then becomes " # φ1 ðsÞ 0 ::: φ2m þ 2 ðsÞ 0 ½NðsÞ ¼ ð39Þ 0 φ2m þ 2 ðsÞ 0 φ1 ðsÞ ::: Fig. 3 plots these Fourier shape functions for m ¼2, where 6 nodes are required. 1.5
1
ϕ2
ϕ1
ϕ3
ϕ4
ϕ5
ϕ6
τxy ¼
ð41Þ
1 1 31 sin 2θ þ sin 4θ þ 4 sin 4θ 2 2r r 2
ð42Þ
where ðr; θÞ are the standard polar coordinates. Table 1 shows the F-SBM solutions for stress sx along the vertical line at x ¼ 0, for approximations of order m ¼2 (6 nodes) and m ¼5 (12 nodes). When the 6 nodes are used, the average error is 1.07%; when 12 nodes are used, the average error is reduced to 0.0044%. It can be seen that the F-SBM obtains high accuracy for the stress solution, even though very few nodes are used, and the accuracy is increased rapidly by increasing the number of nodes. In order to investigate the computing cost, the F-SBM's computing time is compared with the conventional linear SBFEM. 12 nodes are used, and the standard Gauss integration scheme is employed for both F-SBM and linear SBFEM solutions. The calculation is conducted by a PC with I5-2410M 2.3 GHz CPU and 4 GB RAM. The CPU times for the F-SBM and the SBFEM to calculate the coefficient matrix are 0.37 s and 0.34 s, respectively, and the total computing times are 1.04 s and 0.96 s, respectively. The results show computing time of the F-SBM is very close to the conventional SBFEM. In order to evaluate the computed error in the domain, the relative error in energy norm is defined as [11]
0.5
ϕ (s )
1 1 31 cos 2θ cos 4θ 4 cos 4θ 2 2r r 2
syy ¼
L
^ in Eq. (32) can be related to Inverting Eq. (33), the parameters fug the nodal values {uh} by
155
jjujje ¼ 0
jjfsðξ; sÞg fsexact ðξ; sÞgjjΩ 100% jjfsðξ; sÞgjjΩ
ð43Þ
where -0.5
-1 -1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
s Fig. 3. The Fourier shape functions for order m ¼ 2.
0.8
1
jjfsðξ; sÞg fsexact ðξ; sÞgjjΩ Z 1=2 ¼ ffsðξ; sÞg fsexact ðξ; sÞggT ½D 1 ffsðξ; sÞg fsexact ðξ; sÞggdΩ Ω
ð44Þ
156
Y. He et al. / Engineering Analysis with Boundary Elements 41 (2014) 152–159
Fig. 4. A plate of infinite extent subjected to a uniaxial stress: (a) geometry and (b) scaled boundary model.
Table 1 F-SBM solution for sx at x ¼ 0. r
F-SBM (6 nodes)
Error%
F-SBM (12 nodes)
Error%
Exact
1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0
2.981820 2.046881 1.622595 1.403716 1.279443 1.203442 1.154132 1.120568 1.096804 1.079414 1.066330
0.61 1.14 1.39 1.43 1.37 1.26 1.13 1.01 0.90 0.80 0.72
2.999866 2.070494 1.645466 1.424111 1.297140 1.218690 1.167287 1.131973 1.106751 1.088146 1.074045
0.0044 0.0052 0.0059 0.0058 0.0054 0.0049 0.0044 0.0038 0.0033 0.0030 0.0027
3.000000 2.070601 1.645564 1.424194 1.297210 1.218750 1.167338 1.132016 1.106788 1.088179 1.074074
Average error
10
1.07
0.0044
0
Fig. 6. Heat conduction in an L-shape domain.
Error in energy norm
10
10
10
10
10
10
-1
-2
-3
-4
Fourier SBM SBFEM linear element SBFEM quadratic element EFG-SBM quadratic basis
-5
-6
10
0
10
1
10
2
DOF Fig. 5. Convergence of F-SBM, SBFEM and EFG-SBM for infinite plate with circular hole.
is the energy norm of the error and jjfsðξ; sÞgjjΩ ¼
Z Ω
1=2 fsðξ; sÞgT ½D 1 fsðξ; sÞgdΩ
Fig. 7. The scaled boundary model of L-shape plate.
ð45Þ
is the energy norm of the numerical solution. Fig. 5 shows the variation of relative error in energy norm with degrees of freedom (DOF) in the domain ξ ¼ ½1; 3. The solutions of the new F-SBM are compared with the SBFEM using linear and
quadratic elements and the EFG-SBM with a quadratic basis, recently developed by authors [15]. As the number of DOF is the same at each stage of refinement, it is clear that better accuracy and faster convergence can be obtained by using the F-SBM method rather than the linear or quadratic SBFEM approach.
Y. He et al. / Engineering Analysis with Boundary Elements 41 (2014) 152–159
157
Table 2 The results of qr using F-SBM. r
F-SBM (31 nodes)
Error%
F-SBM (55 nodes)
Error%
Exact solutions (W/m oC)
1e 2 1e 4 1e 6 1e 8 1e 10 1e 12 1e 14 1e 16 1e 18 1e 20
2.679823455 12.43870544 57.73566644 267.9866642 1243.890590 5773.659689 26799.09830 124391.0636 577375.2736 2679953.019
3.32e 5 5.70e 4 1.10e 3 1.64e 3 2.18e 3 2.71e 3 3.25e 3 3.79e 3 4.33e 3 4.86e 3
2.679822562 12.43863451 57.73502719 267.9822581 1243.863460 5773.502760 26798.22601 124386.3469 577350.2801 2679822.620
3.73e 8 2.41e 7 4.67e 7 6.72e 7 9.65e 7 1.18e 6 1.42e 6 1.69e 6 1.89e 6 2.13e 6
2.679822563 12.43863448 57.73502692 267.9822563 1243.863448 5773.502692 26798.22563 124386.3448 577350.2692 2679822.563
2.40e 3
Average error
1.07e 6
Relative error of heat flux %
100
10-2
10-4
10-6
10-8
Fourier SBM SBFEM linear EFG-SBM linear
10-10 101
102
DOF
Relative error of heat flux %
100
10-2
Fig. 9. Infinite plate with a through crack: geometry and loads.
10-4
10-6
10-8
Fourier SBM SBFEM linear EFG-SBM linear
10-10 101
102
DOF Fig. 8. Convergence of relative error for qr comparing F-SBM with SBFEM and EFG-SBM. (a) r¼ 1/10, θ¼ 5/4π and (b) r¼ 1/100, θ¼ 3/4π.
Comparing the F-SBM and the EFG-SBM, it can be seen that the accuracy is similar when the number of DOF are small, but the convergence of F-SBM is much faster than EFG-SBM, and thus higher accuracy can be obtained as the number of DOF is increased. The reason for F-SBM's superior performance in this example is that the trigonometric functions in its approximation can capture the variation of the exact solution more precisely in the circumferential direction.
Fig. 10. Scaled boundary model of an infinite plate with a through crack.
4.2. Example 2 – steady-state heat transfer in an L-shape domain The second example deals with the steady-state heat transfer problem in an L-shaped domain illustrated in Fig. 6. The exact
158
Y. He et al. / Engineering Analysis with Boundary Elements 41 (2014) 152–159
Table 3 The results of SIF using F-SBM. Number of nodes
F-SBM
Error%
SBFEM (Linear)
Error%
EFG-SBM (Linear)
Error%
Exact solutions
13 21 28 37 45 53
1.773047828 1.771909603 1.772534634 1.772441208 1.772455534 1.772452866
3.35e 2 3.07e 2 4.56e 3 7.13e 4 9.49e 5 5.55e 5
1.765973947 1.770015354 1.771193770 1.771691632 1.771443822 1.771456952
3.65e 1 1.37e 1 7.11e 2 4.30e 2 5.69e 2 5.62e 2
1.770648223 1.772806367 1.772432519 1.772372556 1.772495119 1.772400534
1.02e 1 1.99e 2 1.20e 3 4.58e 3 2.33e 3 3.01e 3
1.772453851
solution for this example contains a thermal singularity at the re-entrant corner, and is given in the polar coordinate system with the origin O by [22] Tðr; θÞ ¼ r 2=3
2θ π sin 3
ð46Þ
The specified boundary conditions are also indicated in Fig. 5, with zero Dirichlet conditions at part of the boundary that lies on axis and Neumann conditions determined by the exact solution Tðr; θÞ elsewhere. The exact solution of heat flux on the radial direction qr is 2 2θ π qr ¼ k r 1=3 sin 3 3
ð47Þ
where k is the thermal conductivity, taken as unit value in this example. Eq. (47) implies that the radial heat flux qr becomes infinitely large at the re-entrant corner O where r ¼ 0. The SBM model is illustrated in Fig. 7. Nodes are spaced uniformly, and ds is the distance between nodes. Table 2 shows the F-SBM solutions of values of qr in the corner region, where qr is computed along the line with angle θ ¼ π. Two meshes are used. The average relative error is 0.0024% and 0.00000107% for the meshes with 31 nodes and 55 nodes, respectively. Even when the computed position is extremely close to the singular point where qr changes sharply, the computational accuracy is still maintained. Fig. 8(a) and (b) shows the convergence of relative error for qr computed on the points ðr ¼ 1=10; θ ¼ 5π=4Þ and ðr ¼ 1=100; θ ¼ 3π=4Þ, comparing the F-SBM with the SBFEM and the EFG-SBM. The excellent performance of the F-SBM in terms of accuracy and convergence is clearly shown.
4.3. Example 3 – an infinite plate with a through crack The third example refers to the problem of determining the mode I stress intensity factor (SIF) K I for a through crack in an infinite plate, as illustrated in Fig. 9. The applied stress s0 ¼ 1. Due to the symmetry, one quarter of the problem is modelled, as shown in Fig. 10, with the model consisting of a bounded domain, with the scaling centre at the crack tip (point E), and an unbounded domain, with the scaling centre at the middle of crack (point A). The nodes are introduced on the edges AB, BC and CD with uniform pffiffiffiffiffiffiffiffi spacing, ds. The problem has an exact solution, K I ¼ s0 2πa. In Table 3 the F-SBM solutions are compared with the SBFEM with linear elements and the EFG-SBM with linear basis. The results show that the F-SBM achieves high accuracy for SIF, for example, a relative error as low as 0.0000555% can be obtained using 53 nodes. In comparison with SBFEM and EFG-SBM, it can be seen that F-SBM has higher accuracy when the same number of nodes are used.
5. Conclusions A new SBFEM using Fourier shape functions is presented in this paper. The shape functions are based on the Fourier series expansion and augmented with additional linear shape functions terms. By using a transfer matrix, the nodal values are related with Fourier parameters, and in this way the essential boundary conditions can be conveniently handled. In the numerical examples, the new approach has been shown to yield higher accuracy and faster convergence in comparison with the SBFEM using linear or quadratic elements and the EFG-SBM using linear or quadratic basis. In summary a new Fourier scaled boundary method has been developed, and has been shown to offer a convenient and efficient alternative technique for computational mechanics, especially in accurately solving problems involving singularity or unbounded domains. Extension of the method to 3D problems is feasible, and our intention is to examine this in future work.
Acknowledgements The research leading to this paper is supported by National Natural Science Foundation of China [11202046, 10772035, and 11072043], NKBRSF [2010CB832703] and the China Postdoctoral Science Foundation [2012M520617]. References [1] Wolf JP, Song C. Finite-element modelling of unbounded media. Chichester: John Wiley and Sons; 1996. [2] Deeks AJ, Wolf JP. An h-hierarchical adaptive procedure for the scaled boundary finite-element method. Int J Numer Methods Eng 2002;54:585–605. [3] Ooi ET, Yang ZJ. Modelling crack propagation in reinforced concrete using a hybrid finite element-scaled boundary finite element method. Eng Fract Mech 2010;78(2):252–73. [4] Ooi ET, Yang ZJ. A hybrid finite element-scaled boundary finite element method for crack propagation modelling. Comput Methods Appl Mech Eng 2010;199(17–20):1178–92. [5] Ooi ET, Yang ZJ. Modelling dynamic crack propagation using the scaled boundary finite element method. Int J Numer Methods Eng 2011;88 (4):329–49. [6] Yang ZJ, Deeks AJ, Hao H. Transient dynamic fracture analysis using scaled boundary finite element method: a frequency domain approach. Eng Fract Mech 2007;74(5):669–87. [7] Deeks AJ, Wolf JP. Semi-analytical elastostatic analysis of two-dimensional unbounded domains. Int J Numer Anal Methods Geomech Mech 2002;26:1031–57. [8] Doherty JP, Deeks AJ. Scaled boundary finite-element analysis of a nonhomogeneous half-space. Int J Numer Methods Eng 2003;57(7):955–73. [9] Doherty JP, Deeks AJ. Scaled boundary finite-element analysis of a nonhomogeneous axisymmetric domain subjected to general loading. Int J Numer Anal Methods Geomech Mech 2003;27:813–35. [10] Doherty JP, Deeks AJ. Elastic response of circular footings embedded in a nonhomogeneous half-space. Geotechnique 2003;53(8):703–14. [11] Deeks AJ, Augarde CE. A meshless local Petrov–Galerkin scaled boundary method. Comput Mech 2005;36:159–70. [12] Vu TH, Deeks AJ. Use of higher order shape functions in the scaled boundary finite element method. Int J Numer Methods Eng 2006;65:1714–33. [13] Vu TH, Deeks AJ. A p-hierarchical adaptive procedure for the scaled boundary finite element method. Int J Numer Methods Eng 2008;73:47–70. [14] Vu TH, Deeks AJ. A p-adaptive scaled boundary finite element method based on maximization of the error decrease rate. Comput Mech 2008;41:441–55. [15] He YQ, Yang HT, Deeks AJ. An Element-free Galerkin (EFG) scaled boundary method. Finite Elements Anal Des 2012;62:28–36.
Y. He et al. / Engineering Analysis with Boundary Elements 41 (2014) 152–159 [16] Guan Y, Pourboghrat F, Yu WR. Fourier series based finite element analysis of tube hydroforming: a axisymmetric model. Eng Comput 2006;23(7):697–728. [17] Hamzeh Javaran S, Khaji N, Moharrami H. A dual reciprocity BEM approach using new Fourier radial basis functions applied to 2D elastodynamic transient analysis. Eng Anal Bound Elements 2011;35:85–95. [18] Khaji N, Hamzehei Javaran S. New complex Fourier shape functions for the analysis of two-dimensional potential problems using boundary element method. Eng Anal Bound Elements 2013;37:260–72. [19] Song C, Wolf JP. The scaled boundary finite element method – alias consistent infinitesimal finite element cell method for diffusion. Int J Numer Methods Eng 1999;45:1403–31.
159
[20] Lin Z, Liao S. The scaled boundary FEM for nonlinear problems. Commun Nonlinear Sci Numer Simul 2011;16:63–75. [21] Timoshenko S, Goodier JN. Theory of elasticity. 2nd edition. New York: McGraw Hill; 1951. [22] Dörfel MR, Jüttler B, Simeon B. Adaptive isogeometric analysis by local h-refinement with T-splines. Comput Methods Appl Mech Eng 2010; 199(5–8):264–75.
Contents lists available at ScienceDirect
Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound
Use of Fourier shape functions in the scaled boundary method Yiqian He a, Haitian Yang a, Andrew J. Deeks b,n a State Key Lab of Structural Analysis for Industrial Equipment, Department of Engineering Mechanics, Dalian University of Technology, Dalian 116024, PR China b President's Office, University College Dublin, Belfield Campus, Dublin 4, Ireland
art ic l e i nf o
a b s t r a c t
Article history: Received 21 November 2013 Accepted 20 January 2014 Available online 14 February 2014
The scaled boundary finite element method (SBFEM) is a semi-analytical method, whose versatility, accuracy and efficiency are not only equal to, but potentially better than the finite element method and the boundary element method for certain problems. This paper investigates the possibility of using Fourier shape functions in the SBFEM to form the approximation in the circumferential direction. The shape functions effectively form a Fourier series expansion in the circumferential direction, and are augmented by additional linear shape functions. The proposed method is evaluated by solving three elastostatic and steady-state heat transfer problems. The accuracy and convergence of the proposed method is demonstrated, and the performance is found to be better than using polynomial elements or using an element-free Galerkin approximation for the circumferential approximation in the scaled boundary method. & 2014 Elsevier Ltd. All rights reserved.
Keywords: Scaled boundary method Fourier shape functions Computational accuracy Stress singularities Unbounded domains
1. Introduction The scaled boundary method (SBM) is a semi-analytical method developed relatively recently by Wolf and Song [1]. The method introduces a normalised radial coordinate system based on a scaling centre and a defining curve (usually taken as the boundary). The governing deferential equations are weakened in the circumferential direction and then solved analytically in the normalised radial direction. The SBM has combined the advantages of the Finite Element Method (FEM) and the Boundary Element Method (BEM), and no fundamental solution is required like BEM. In addition, the SBM has been shown to be more efficient than the FEM for problems involving unbounded domains and for problems involving stress singularities or discontinuities [2]. Effective applications of this method have been demonstrated in various problem domains, including fracture problems [3–6] and foundation problems [7–10]. In the scaled boundary method, the discretisation approach used in the circumferential direction has significant influence on the accuracy of the resulting solutions [11]. The most commonly used method for performing this circumferential discretisation is the finite element approach, leading to the method called the scaled boundary finite element method (SBFEM). Vu and Deeks [12–14] investigated the use of higher-order polynomial shape functions in the SBFEM, and demonstrated the SBFEM converged significantly faster under p-refinement than under h-refinement.
n
Corresponding author. E-mail address: [email protected] (A.J. Deeks).
0955-7997/$ - see front matter & 2014 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.enganabound.2014.01.012
The development of meshless methods provided another approach to building circumferential approximations for the scaled boundary method. Deeks and Augarde [11] developed a Meshless Local Petrov–Galerkin method scaled boundary method (MLPG-SBM) and He et al. [15] developed an Element-free Galerkin scaled boundary method (EFG-SBM). This work showed that these two meshless scaled boundary methods gave a higher level of accuracy and rate of convergence than the conventional SBFEM using linear or quadratic elements, with the EFG-SBM performing slightly better than the MLPG-SBM. In this paper, the possibility of using shape functions based on the terms of a Fourier series for the circumferential approximation of the SBFEM is investigated. Fourier interpolations containing trigonometric functions have been applied to both the finite element method (FEM) and the boundary element method (BEM). For example, Guan et al. [16] developed a Fourier series based FEM into for the analysis of tube hydroforming, and showed that this Fourier shape function reduced the number of degrees of freedom required. Khaji et al. [17,18] applied Fourier radial basis functions into the BEM, and showed that of the resulting BEM is much more accurate than the BEM using classic Lagrange shape functions. Although the advantages of Fourier based FEM and BEM have been illustrated in previous work, to date there has been no work reported on the use of Fourier shape functions in the SBFEM. A new Fourier-based scaled boundary method (F-SBM) is presented in this paper. A set of shape functions based on Fourier series expansion is derived, and augmented with linear shape functions. The new shape functions provide good approximation to both trigonometric and polynomial functions in the circumferential direction of the
Y. He et al. / Engineering Analysis with Boundary Elements 41 (2014) 152–159
scaled boundary system. In three numerical examples, the F-SBM is used to solve two-dimensional elastostatic and steady-state heat transfer problems. The accuracy and convergence of F-SBM is compared with the conventional SBFEM using both linear and quadratic elements and with the EFG-SBM. Superior performance in terms of both accuracy and convergence is demonstrated. This paper is organised as follows: The basic equations of scaled boundary method are given in Section 2. Section 3 introduces the Fourier shape functions for use in the scaled boundary method. Three example problems are presented in Section 4 to verify the effectiveness of proposed method, and the paper draws conclusions at the end.
The scaled boundary method introduces a normalised radial coordinate system by scaling a defining curve (usually the domain boundary or a part of the domain boundary) relative to a scaling centre ðx0 ; y0 Þ selected within the domain or at the intersection of two straight sections of the boundary (Fig. 1). The normalised radial coordinate ξ runs from the scaling centre towards the defining curve, and has values of zero at the scaling centre and unity at the defining curve. The other circumferential coordinate s specifies a distance around the defining curve from an origin on the curve. The scaled boundary and Cartesian coordinate systems are related by the scaling equations x ¼ x0 þ ξxs ðsÞ
ð1Þ
y ¼ y0 þ ξys ðsÞ
ð2Þ
Displacement and stress components are retained in the original Cartesian coordinate directions, while position is specified in terms of the scaled boundary coordinates. An approximate solution is sought in the form n
fuh ðξ; sÞg ¼ ∑ ½Ni ðsÞuhi ðξÞ ¼ ½NðsÞfuh ðξÞg
ð3Þ
i¼1
This represents a discretisation of the part of the boundary located at ξ ¼ 1 with the shape function [N(s)]. The unknown vector fuh ðξÞg is a set of n functions analytical in ξ. The shape functions apply for all lines with a constant ξ. (If the scaling centre lies on the boundary, as in Fig. 1, the straight portions of the boundary adjacent to the scaling centre and representing radial lines are not discretised, and in the solution process an analytical solution is found along these lines.) Mapping the linear operator to the scaled boundary coordinate system using standard methods ½L ¼ ½L1
∂ ∂ ∂ 1 2 ∂ 1 þ ½L2 ¼ ½b ðsÞ þ ½b ðsÞ ∂x ∂y ∂ξ ξ ∂s 1
The stresses are obtained by multiplying the strains (obtained form the displacement field using the linear operator) by the elasticity matrix ½D in the form
ð4Þ
2
where ½b ðsÞ and ½b ðsÞ are dependent only on the boundary definition.
1 ξ
sðξ; sÞ ¼ ½D εðξ; sÞ ¼ ½D½B1 ðsÞfuðξÞg;ξ þ ½D½B2 ðsÞ uh ðξÞ
ð5Þ
where 1
ð6Þ
2
ð7Þ
½B1 ðsÞ ¼ ½b ðsÞ½NðsÞ ½B2 ðsÞ ¼ ½b ðsÞ½NðsÞ;s Z V
2. The scaled boundary method
153
In this case the virtual work statement becomes Z fδεðξ; sÞgT fsh ðξ; sÞgdV fδuðsÞgT ftðsÞgds ¼ 0
ð8Þ
S
where the first term represents the internal work and the second term the external work, and ftðsÞg is the external force vector. The virtual strain field is of the form (analogous to Eq. (5)) 1 δεðξ; sÞ ¼ ½B1 ðsÞfδuðξÞg;ξ þ ½B2 ðsÞ δuðξÞ ξ
ð9Þ
where fδuðξÞg is virtual displacement and dV ¼ jJjξdξds ð10Þ where J is the Jacobian at the boundary ðξ ¼ 1Þ. Substituting Eqs. (5), (9) and (10), integrating the area integrals containing fδuðξÞg;ξ with respect to ξ using Green's Theorem, and introducing the coefficient matrices Z ð11Þ ½E0 ¼ ½B1 ðsÞT ½D½B1 ðsÞjJjds S
Z ½E1 ¼
½B2 ðsÞT ½D½B1 ðsÞjJjds
ð12Þ
½B2 ðsÞT ½D½B2 ðsÞjJjds
ð13Þ
S
Z ½E2 ¼ S
the virtual work equation may be expressed as Z n o fδεðξ; sÞgT fsh ðξ; sÞgdV ¼ fδugT ½E0 fuh g;ξ þ ½E1 fuh g V
Z 0
1
1 fδuðξÞgT ½E0 ξfuh ðξÞg;ξξ þ ½½E0 þ ½E1 T ½E1 fuh ðξÞg;ξ ½E2 uh ðξÞ dξ ξ
ð14Þ On substitution of Eq. (3), the external virtual work term becomes Z Z fδuðsÞgT ftðsÞgds ¼ fδugT fNðsÞgftðsÞgds ¼ fδugT fPg ð15Þ S
S
Thus the complete virtual work equation becomes fδugT f½E0 fuh g;ξ þ½E1 fuh gg fδugT fPg Z 1 T n 0 δuðξÞ ½E ξfuh ðξÞg;ξξ þ ½½E0 þ ½E1 T 0 1 ½E1 fuh ðξÞg;ξ ½E2 fuh ðξÞg dξ ¼ 0g ξ
ð16Þ
In order for Eq. (16) to be satisfied for all fδuðξÞg (implying that equilibrium it closely satisfied in the radial direction and in the approximate sense in the circumferential direction), both of the following conditions must be satisfied. fPg ¼ ½E0 fuh g;ξ þ ½E1 T fuh g ½E0 ξ
2
ð17Þ
uh ðξÞ ;ξξ þ ½½E0 þ ½E1 T ½E1 ξ uh ðξÞ ;ξ ½E2 uh ðξÞ ¼ f0g ð18Þ
Fig. 1. Bounded domain with side faces showing scaled boundary coordinate system.
By inspection, solutions to the homogeneous set of Euler– Cauchy differential equations represented by Eq. (18) must be of
154
Y. He et al. / Engineering Analysis with Boundary Elements 41 (2014) 152–159
the form fϕ1 g þ c2 ξ
λ2
fϕ2 g þ :::
ð19Þ
where the exponents λi and corresponding vectors fϕi g may be interpreted as independent modes of deformation which closely satisfy internal equilibrium in the ξ direction. The integration constants ci represent the contribution of each mode to the solution, and are dependent on the boundary conditions. The displacements for each mode take the form (omitting the subscript) fuðξÞg ¼ ξ λ fϕg
trigonometric functions as [17] 1 nπ nπ f ðrÞ ¼ a0 þ ∑ an cos r þ bn sin r Lmax Lmax n¼1
ð30Þ
where a0, an, bn and Lmax represent the Fourier series parameters. Thus on the boundary at ξ ¼ 1, the displacement could be approximated as m jπ jπ s þ bj sin s ð31Þ uh ðsÞ ¼ a0 þ ∑ aj cos L L j¼1
ð20Þ
The vector fϕg can be identified as the modal displacements at the boundary nodes, while λ can be identified as a modal scaling factor for the ‘radial’ direction. Substituting this solution into Eq. (18) yields the quadratic eigenproblem ½λ2 ½E0 λ½½E1 T ½E1 ½E2 fϕg ¼ f0g
ð21Þ
The equivalent nodal forces required at the boundary to equilibrate each displacement mode are obtained by substituting Eq. (20) into Eq. (17) (which is evaluated at ξ ¼ 1) as fqg ¼ ½½E1 T λ½E0 fϕg
ð22Þ
The Eqs. (21)–(22) can be expressed in linear form " #( ) ( ) ϕ ϕ ½E0 1 ½E0 1 ½E1 T ¼ λ q q ½E1 ½E0 1 ½E1 T ½E2 ½E1 ½E0 1
ð23Þ
Solution of this standard eigenproblem yields 2n modes. The eigenvectors contain the modal displacements and the equivalent modal node forces. For a bounded domain only the modes with non-positive real components of λ lead to finite displacements or the scaling centre. This subset of n modal displacements is designated by ½Φ1 , where the vectors in the set from the columns of the matrix. The subset of modal force vectors corresponding to the n modes in ½Φ1 is denoted as ½Q 1 . For any set of boundary node displacements ½uh , the integration constants required to satisfy Eq. (18) on the boundary ðξ ¼ 1Þ are fcg ¼ ½Φ1 1 fuh g
ð24Þ
The equivalent nodal forces required to cause these displacements are fPg ¼ ½Q 1 fcg ¼ ½Q 1 ½Φ1
1
fuh g
where s is the circumferential coordinate in scaled boundary element, L is the length of the boundary and m represents the order of Fourier series. However, for many practical problems, the functions describing the defining curve in the scaled boundary method, xs(s) and ys(s), may be only C0 continuous at certain points along the boundary, rather than fully continuous. These points are termed vertices, in line with conventional usage of the term. Although the derivatives of a solution field with respect to the Cartesian coordinate system may be continuous, derivatives of the Cartesian components of this displacement field with respect to the circumferential coordinate s will not be continuous at the vertices. This is illustrated in Fig. 2. Consequently it is not advantageous to use a continuous Fourier series representation across the vertices. Between the vertices, the derivatives of the solution field are expected to be continuous. The defining curve can thus be defined to consist of edges, over which the solution is expected to be smooth, connected together at vertices, at which the components of the solution field are expected to be C0 continuous with respect to s, but at which the derivatives of these components with respect to s are discontinuous. To preserve the C0 continuity, nodes are introduced at the vertices. Along each
1 0.8 u(x,y) = x.y
fuh ðξÞg ¼ c1 ξ
λ1
0.6 0.4 0.2 0 1
ð25Þ
1
The stiffness matrix of the domain is therefore ½K ¼ ½Q 1 ½Φ1 1
0.5
ð26Þ
0.5
y
0
and the equilibrium requirement is reduced to ½Kfuh g fPg ¼ f0g
ð27Þ
1
The integration constants are then obtained using Eq. (24) and the displacement field is recovered by as n
fuh ðξ; sÞg ¼ ½NðsÞ ∑ ci ξ i¼1
λi
fϕi g
n
0.8
ð28Þ 0.6 u(s)
The stress field is obtained as fsh ðξ; sÞg ¼ ½D ∑ ci ξ λi 1 ½ λi ½B1 ðsÞ þ ½B2 ðsÞfϕi g
s
x
0
ð29Þ
0.4
i¼1
0.2
3. A Fourier shape function
0 0
This paper employs shape functions obtained from the wellknown Fourier series. Based on the theory of the Fourier series, any continuous function f ðrÞ maybe represented by a series of
0.5
1
1.5
2 s
2.5
3
3.5
4
Fig. 2. (a) Smooth function in Cartesian coordinates. (b) Same function plotted against boundary coordinate s.
Y. He et al. / Engineering Analysis with Boundary Elements 41 (2014) 152–159
edge a single element with the shape functions defined using the Fourier approach is introduced. To preserve C0 continuity between the edges or elements, linear polynomial functions terms are added into the standard Fourier approximation as m Ls s jπ jπ þ β1 þ ∑ aj cos s þbj sin s ð32Þ uh ðsÞ ¼ α1 L L j¼1 L L where α1 and β1 represent the values of the function at the end nodes of the element. While it is possible to use the Fourier parameters as the unknown boundary parameters when solving the scaled boundary finite element equations, here the parameters in the Fourier expansion above are transferred to nodal values at equally spaced nodes along each element for ease of applying essential boundary conditions and enforcing C0 continuity between elements. If ð2m þ 2Þ nodes are used, the nodal values vector fug can be related to the parameters in Eq. (32) by ^ fug ¼ ½Tfug
ð33Þ
^ ¼ where fug
n
α1
a1
⋯
am
⋯
b1
bm
β1
oT
, and [T] is a
transfer matrix assembled as T ij ¼ ψ j ðSi Þ
i; j ¼ 1; 2 … 2m þ 2
ð34Þ
where Si is the circumferential coordinate of the ith node, and the component functions of the Fourier expansion are 8 Ls i¼1 > > > L
> > ði 1Þπ > 2rirmþ1 < cos L s
ð35Þ ψ i ðsÞ ¼ > > sin ði 1L mÞπ s m þ 2 ri r 2m þ 1 > > > > :s i ¼ 2m þ 2
4. Performance of the method Here three examples are provided to demonstrate the effectiveness of the proposed method, two elastostatic examples and one steady-state heat transfer example. (The solution procedure of scaled boundary method for steady-state heat transfer problems is similar to that for elastostatic problems provided in the previous sections, and the detailed derivations can be found in references [19,20].)
4.1. Example 1 – a plate of infinite extent subjected to a uniaxial stress The first example is a plate of infinite extent containing a hole of unit radius and subjected to a uniaxial stress of unity. This example demonstrates the ability of the method to model unbounded problems efficiently and accurately. Plane strain conditions are assumed. Considering the symmetry of the problem, one quarter of the plate is represented (Fig. 4(a)). The Young's modulus is take as E ¼ 1000 and Poisson's ratio v ¼ 0:3. The boundary of the hole is considered as a single edge with the length l (Fig. 4(b)). Uniform nodes are introduced along this edge, the nodes spacing is specified as ds. The centre of the hole is used as the scaling centre. The exact solution for the stresses for this problem [21] is 1 3 31 cos 2θ þ cos 4θ þ 4 cos 4θ sxx ¼ 1 2 ð40Þ 2r r 2
^ ¼ ½T 1 fug fug
ð36Þ
Thus the approximation for displacement can be rewritten as uh ðsÞ ¼ fψg½T 1 fug
ð37Þ
The shape functions relating to the nodal displacements are hence fφg ¼ fψg½T 1
ð38Þ
and the shape function matrix for the scaled boundary method then becomes " # φ1 ðsÞ 0 ::: φ2m þ 2 ðsÞ 0 ½NðsÞ ¼ ð39Þ 0 φ2m þ 2 ðsÞ 0 φ1 ðsÞ ::: Fig. 3 plots these Fourier shape functions for m ¼2, where 6 nodes are required. 1.5
1
ϕ2
ϕ1
ϕ3
ϕ4
ϕ5
ϕ6
τxy ¼
ð41Þ
1 1 31 sin 2θ þ sin 4θ þ 4 sin 4θ 2 2r r 2
ð42Þ
where ðr; θÞ are the standard polar coordinates. Table 1 shows the F-SBM solutions for stress sx along the vertical line at x ¼ 0, for approximations of order m ¼2 (6 nodes) and m ¼5 (12 nodes). When the 6 nodes are used, the average error is 1.07%; when 12 nodes are used, the average error is reduced to 0.0044%. It can be seen that the F-SBM obtains high accuracy for the stress solution, even though very few nodes are used, and the accuracy is increased rapidly by increasing the number of nodes. In order to investigate the computing cost, the F-SBM's computing time is compared with the conventional linear SBFEM. 12 nodes are used, and the standard Gauss integration scheme is employed for both F-SBM and linear SBFEM solutions. The calculation is conducted by a PC with I5-2410M 2.3 GHz CPU and 4 GB RAM. The CPU times for the F-SBM and the SBFEM to calculate the coefficient matrix are 0.37 s and 0.34 s, respectively, and the total computing times are 1.04 s and 0.96 s, respectively. The results show computing time of the F-SBM is very close to the conventional SBFEM. In order to evaluate the computed error in the domain, the relative error in energy norm is defined as [11]
0.5
ϕ (s )
1 1 31 cos 2θ cos 4θ 4 cos 4θ 2 2r r 2
syy ¼
L
^ in Eq. (32) can be related to Inverting Eq. (33), the parameters fug the nodal values {uh} by
155
jjujje ¼ 0
jjfsðξ; sÞg fsexact ðξ; sÞgjjΩ 100% jjfsðξ; sÞgjjΩ
ð43Þ
where -0.5
-1 -1
-0.8
-0.6
-0.4
-0.2
0
0.2
0.4
0.6
s Fig. 3. The Fourier shape functions for order m ¼ 2.
0.8
1
jjfsðξ; sÞg fsexact ðξ; sÞgjjΩ Z 1=2 ¼ ffsðξ; sÞg fsexact ðξ; sÞggT ½D 1 ffsðξ; sÞg fsexact ðξ; sÞggdΩ Ω
ð44Þ
156
Y. He et al. / Engineering Analysis with Boundary Elements 41 (2014) 152–159
Fig. 4. A plate of infinite extent subjected to a uniaxial stress: (a) geometry and (b) scaled boundary model.
Table 1 F-SBM solution for sx at x ¼ 0. r
F-SBM (6 nodes)
Error%
F-SBM (12 nodes)
Error%
Exact
1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0
2.981820 2.046881 1.622595 1.403716 1.279443 1.203442 1.154132 1.120568 1.096804 1.079414 1.066330
0.61 1.14 1.39 1.43 1.37 1.26 1.13 1.01 0.90 0.80 0.72
2.999866 2.070494 1.645466 1.424111 1.297140 1.218690 1.167287 1.131973 1.106751 1.088146 1.074045
0.0044 0.0052 0.0059 0.0058 0.0054 0.0049 0.0044 0.0038 0.0033 0.0030 0.0027
3.000000 2.070601 1.645564 1.424194 1.297210 1.218750 1.167338 1.132016 1.106788 1.088179 1.074074
Average error
10
1.07
0.0044
0
Fig. 6. Heat conduction in an L-shape domain.
Error in energy norm
10
10
10
10
10
10
-1
-2
-3
-4
Fourier SBM SBFEM linear element SBFEM quadratic element EFG-SBM quadratic basis
-5
-6
10
0
10
1
10
2
DOF Fig. 5. Convergence of F-SBM, SBFEM and EFG-SBM for infinite plate with circular hole.
is the energy norm of the error and jjfsðξ; sÞgjjΩ ¼
Z Ω
1=2 fsðξ; sÞgT ½D 1 fsðξ; sÞgdΩ
Fig. 7. The scaled boundary model of L-shape plate.
ð45Þ
is the energy norm of the numerical solution. Fig. 5 shows the variation of relative error in energy norm with degrees of freedom (DOF) in the domain ξ ¼ ½1; 3. The solutions of the new F-SBM are compared with the SBFEM using linear and
quadratic elements and the EFG-SBM with a quadratic basis, recently developed by authors [15]. As the number of DOF is the same at each stage of refinement, it is clear that better accuracy and faster convergence can be obtained by using the F-SBM method rather than the linear or quadratic SBFEM approach.
Y. He et al. / Engineering Analysis with Boundary Elements 41 (2014) 152–159
157
Table 2 The results of qr using F-SBM. r
F-SBM (31 nodes)
Error%
F-SBM (55 nodes)
Error%
Exact solutions (W/m oC)
1e 2 1e 4 1e 6 1e 8 1e 10 1e 12 1e 14 1e 16 1e 18 1e 20
2.679823455 12.43870544 57.73566644 267.9866642 1243.890590 5773.659689 26799.09830 124391.0636 577375.2736 2679953.019
3.32e 5 5.70e 4 1.10e 3 1.64e 3 2.18e 3 2.71e 3 3.25e 3 3.79e 3 4.33e 3 4.86e 3
2.679822562 12.43863451 57.73502719 267.9822581 1243.863460 5773.502760 26798.22601 124386.3469 577350.2801 2679822.620
3.73e 8 2.41e 7 4.67e 7 6.72e 7 9.65e 7 1.18e 6 1.42e 6 1.69e 6 1.89e 6 2.13e 6
2.679822563 12.43863448 57.73502692 267.9822563 1243.863448 5773.502692 26798.22563 124386.3448 577350.2692 2679822.563
2.40e 3
Average error
1.07e 6
Relative error of heat flux %
100
10-2
10-4
10-6
10-8
Fourier SBM SBFEM linear EFG-SBM linear
10-10 101
102
DOF
Relative error of heat flux %
100
10-2
Fig. 9. Infinite plate with a through crack: geometry and loads.
10-4
10-6
10-8
Fourier SBM SBFEM linear EFG-SBM linear
10-10 101
102
DOF Fig. 8. Convergence of relative error for qr comparing F-SBM with SBFEM and EFG-SBM. (a) r¼ 1/10, θ¼ 5/4π and (b) r¼ 1/100, θ¼ 3/4π.
Comparing the F-SBM and the EFG-SBM, it can be seen that the accuracy is similar when the number of DOF are small, but the convergence of F-SBM is much faster than EFG-SBM, and thus higher accuracy can be obtained as the number of DOF is increased. The reason for F-SBM's superior performance in this example is that the trigonometric functions in its approximation can capture the variation of the exact solution more precisely in the circumferential direction.
Fig. 10. Scaled boundary model of an infinite plate with a through crack.
4.2. Example 2 – steady-state heat transfer in an L-shape domain The second example deals with the steady-state heat transfer problem in an L-shaped domain illustrated in Fig. 6. The exact
158
Y. He et al. / Engineering Analysis with Boundary Elements 41 (2014) 152–159
Table 3 The results of SIF using F-SBM. Number of nodes
F-SBM
Error%
SBFEM (Linear)
Error%
EFG-SBM (Linear)
Error%
Exact solutions
13 21 28 37 45 53
1.773047828 1.771909603 1.772534634 1.772441208 1.772455534 1.772452866
3.35e 2 3.07e 2 4.56e 3 7.13e 4 9.49e 5 5.55e 5
1.765973947 1.770015354 1.771193770 1.771691632 1.771443822 1.771456952
3.65e 1 1.37e 1 7.11e 2 4.30e 2 5.69e 2 5.62e 2
1.770648223 1.772806367 1.772432519 1.772372556 1.772495119 1.772400534
1.02e 1 1.99e 2 1.20e 3 4.58e 3 2.33e 3 3.01e 3
1.772453851
solution for this example contains a thermal singularity at the re-entrant corner, and is given in the polar coordinate system with the origin O by [22] Tðr; θÞ ¼ r 2=3
2θ π sin 3
ð46Þ
The specified boundary conditions are also indicated in Fig. 5, with zero Dirichlet conditions at part of the boundary that lies on axis and Neumann conditions determined by the exact solution Tðr; θÞ elsewhere. The exact solution of heat flux on the radial direction qr is 2 2θ π qr ¼ k r 1=3 sin 3 3
ð47Þ
where k is the thermal conductivity, taken as unit value in this example. Eq. (47) implies that the radial heat flux qr becomes infinitely large at the re-entrant corner O where r ¼ 0. The SBM model is illustrated in Fig. 7. Nodes are spaced uniformly, and ds is the distance between nodes. Table 2 shows the F-SBM solutions of values of qr in the corner region, where qr is computed along the line with angle θ ¼ π. Two meshes are used. The average relative error is 0.0024% and 0.00000107% for the meshes with 31 nodes and 55 nodes, respectively. Even when the computed position is extremely close to the singular point where qr changes sharply, the computational accuracy is still maintained. Fig. 8(a) and (b) shows the convergence of relative error for qr computed on the points ðr ¼ 1=10; θ ¼ 5π=4Þ and ðr ¼ 1=100; θ ¼ 3π=4Þ, comparing the F-SBM with the SBFEM and the EFG-SBM. The excellent performance of the F-SBM in terms of accuracy and convergence is clearly shown.
4.3. Example 3 – an infinite plate with a through crack The third example refers to the problem of determining the mode I stress intensity factor (SIF) K I for a through crack in an infinite plate, as illustrated in Fig. 9. The applied stress s0 ¼ 1. Due to the symmetry, one quarter of the problem is modelled, as shown in Fig. 10, with the model consisting of a bounded domain, with the scaling centre at the crack tip (point E), and an unbounded domain, with the scaling centre at the middle of crack (point A). The nodes are introduced on the edges AB, BC and CD with uniform pffiffiffiffiffiffiffiffi spacing, ds. The problem has an exact solution, K I ¼ s0 2πa. In Table 3 the F-SBM solutions are compared with the SBFEM with linear elements and the EFG-SBM with linear basis. The results show that the F-SBM achieves high accuracy for SIF, for example, a relative error as low as 0.0000555% can be obtained using 53 nodes. In comparison with SBFEM and EFG-SBM, it can be seen that F-SBM has higher accuracy when the same number of nodes are used.
5. Conclusions A new SBFEM using Fourier shape functions is presented in this paper. The shape functions are based on the Fourier series expansion and augmented with additional linear shape functions terms. By using a transfer matrix, the nodal values are related with Fourier parameters, and in this way the essential boundary conditions can be conveniently handled. In the numerical examples, the new approach has been shown to yield higher accuracy and faster convergence in comparison with the SBFEM using linear or quadratic elements and the EFG-SBM using linear or quadratic basis. In summary a new Fourier scaled boundary method has been developed, and has been shown to offer a convenient and efficient alternative technique for computational mechanics, especially in accurately solving problems involving singularity or unbounded domains. Extension of the method to 3D problems is feasible, and our intention is to examine this in future work.
Acknowledgements The research leading to this paper is supported by National Natural Science Foundation of China [11202046, 10772035, and 11072043], NKBRSF [2010CB832703] and the China Postdoctoral Science Foundation [2012M520617]. References [1] Wolf JP, Song C. Finite-element modelling of unbounded media. Chichester: John Wiley and Sons; 1996. [2] Deeks AJ, Wolf JP. An h-hierarchical adaptive procedure for the scaled boundary finite-element method. Int J Numer Methods Eng 2002;54:585–605. [3] Ooi ET, Yang ZJ. Modelling crack propagation in reinforced concrete using a hybrid finite element-scaled boundary finite element method. Eng Fract Mech 2010;78(2):252–73. [4] Ooi ET, Yang ZJ. A hybrid finite element-scaled boundary finite element method for crack propagation modelling. Comput Methods Appl Mech Eng 2010;199(17–20):1178–92. [5] Ooi ET, Yang ZJ. Modelling dynamic crack propagation using the scaled boundary finite element method. Int J Numer Methods Eng 2011;88 (4):329–49. [6] Yang ZJ, Deeks AJ, Hao H. Transient dynamic fracture analysis using scaled boundary finite element method: a frequency domain approach. Eng Fract Mech 2007;74(5):669–87. [7] Deeks AJ, Wolf JP. Semi-analytical elastostatic analysis of two-dimensional unbounded domains. Int J Numer Anal Methods Geomech Mech 2002;26:1031–57. [8] Doherty JP, Deeks AJ. Scaled boundary finite-element analysis of a nonhomogeneous half-space. Int J Numer Methods Eng 2003;57(7):955–73. [9] Doherty JP, Deeks AJ. Scaled boundary finite-element analysis of a nonhomogeneous axisymmetric domain subjected to general loading. Int J Numer Anal Methods Geomech Mech 2003;27:813–35. [10] Doherty JP, Deeks AJ. Elastic response of circular footings embedded in a nonhomogeneous half-space. Geotechnique 2003;53(8):703–14. [11] Deeks AJ, Augarde CE. A meshless local Petrov–Galerkin scaled boundary method. Comput Mech 2005;36:159–70. [12] Vu TH, Deeks AJ. Use of higher order shape functions in the scaled boundary finite element method. Int J Numer Methods Eng 2006;65:1714–33. [13] Vu TH, Deeks AJ. A p-hierarchical adaptive procedure for the scaled boundary finite element method. Int J Numer Methods Eng 2008;73:47–70. [14] Vu TH, Deeks AJ. A p-adaptive scaled boundary finite element method based on maximization of the error decrease rate. Comput Mech 2008;41:441–55. [15] He YQ, Yang HT, Deeks AJ. An Element-free Galerkin (EFG) scaled boundary method. Finite Elements Anal Des 2012;62:28–36.
Y. He et al. / Engineering Analysis with Boundary Elements 41 (2014) 152–159 [16] Guan Y, Pourboghrat F, Yu WR. Fourier series based finite element analysis of tube hydroforming: a axisymmetric model. Eng Comput 2006;23(7):697–728. [17] Hamzeh Javaran S, Khaji N, Moharrami H. A dual reciprocity BEM approach using new Fourier radial basis functions applied to 2D elastodynamic transient analysis. Eng Anal Bound Elements 2011;35:85–95. [18] Khaji N, Hamzehei Javaran S. New complex Fourier shape functions for the analysis of two-dimensional potential problems using boundary element method. Eng Anal Bound Elements 2013;37:260–72. [19] Song C, Wolf JP. The scaled boundary finite element method – alias consistent infinitesimal finite element cell method for diffusion. Int J Numer Methods Eng 1999;45:1403–31.
159
[20] Lin Z, Liao S. The scaled boundary FEM for nonlinear problems. Commun Nonlinear Sci Numer Simul 2011;16:63–75. [21] Timoshenko S, Goodier JN. Theory of elasticity. 2nd edition. New York: McGraw Hill; 1951. [22] Dörfel MR, Jüttler B, Simeon B. Adaptive isogeometric analysis by local h-refinement with T-splines. Comput Methods Appl Mech Eng 2010; 199(5–8):264–75.