Improved hybrid boundary solution for shell elements

Engineering Analysis with Boundary Elements 71 (2016) 70–78

Contents lists available at ScienceDirect

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

Improved hybrid boundary solution for shell elements Taha H.A. Naga a, Youssef F. Rashed b,n a b

Department of engineering mathematics and physics, Faculty of engineering – Shoubra, Benha University, Egypt Supereme Councile of Universities, Egypt, and Department of Structural Engineering, Cairo University, Giza, Egypt

art ic l e i nf o

a b s t r a c t

Article history: Received 1 June 2016 Received in revised form 12 July 2016 Accepted 19 July 2016

In this paper, a new finite element for shell structures is developed using an improved hybrid boundary solution. First the variational boundary integral equation for shear-deformable plate bending problems is developed based on quadratic boundary elements. Hence such a formulation is coupled to a similar formulation for 2D plane stress problems to produce the developed shell elements. Numerical examples are presented to demonstrate the accuracy and validity of the proposed formulation. & 2016 Published by Elsevier Ltd.

Keywords: Finite element method Boundary element method Variational formulation Stiffness matrix Shear-deformable plates Shell structures

1. Introduction Shell structures gain huge attention especially in structural engineering [1]. Engineers, nowadays, use several commercially based finite element packages (such as, SAP [2], ANSYS [3], etc.) to analyze such problems. Most of shell elements used in these packages are based on the hybrid integral equation formulation of Pian and Tong [4]. Alternatively, few packages were developed based on alternative technology, the boundary element method, such as the PLPAK [5]. The mathematical base behind this package [5] is the direct boundary integral equation formulation based on the well-known virtual work statement, in which the virtual state is the fundamental problem, or the solution due to unit generalized load in an infinite domain. This solution is called the fundamental solution [6]. It has to be noted that Zienkiewicz [7] (page 646) in 1977, was the first textbook referred to use fundamental solutions as trial functions within finite elements to form the “Boundary Solution procedures”. However, the solution presented in Ref. [7] was based on the traditional energy based functional (without the term that forces the compatibility between the domain and boundary displacements along the boundary). This term was introduced in Tong [4]. Zienkiewicz in Ref. [7] solved 2D problems and since that time, most of finite element formulations employed polynomial based n

Corresponding author. E-mail addresses: [email protected] (T.H.A. Naga), [email protected], [email protected] (Y.F. Rashed). http://dx.doi.org/10.1016/j.enganabound.2016.07.011 0955-7997/& 2016 Published by Elsevier Ltd.

trial functions. In 1989, Dumont [8] and Tania [9] developed boundary element formulation based on the variational functional of Tong [4]. Their formulation has several advantages, such as it leads to similar formulation to finite elements in terms of sparse symmetric stiffness matrices. They considered this formulation for potential problems, and 2D elasticity problems. One of the disadvantages of their formulation is the appearance of higher order singularities; which led to use equally divided boundary elements (constant or quadratic) to cancel such singularities naturally. Several papers have then been appeared in this context, such as the work of Leung et al. [10] to couple boundary element method to finite element method, Dumont [11] for hybrid stress formulation, Naga and Rashed [12] for shear-deformable plate bending problems with constant elements. It has to be noted that all these publications considered this variational formulation as new formulation inside the boundary element method. Qin [13–17] developed new fundamental solution-based finite elements; which employs the formulation of Tong [4]. For example, in Ref. [13], Qin developed a new type of hybrid finite element formulation for solving 2D heat conduction problems. In Ref. [14], Qin considered the solution of two-dimensional orthotropic elasticity. Qin [15] developed a new type of finite elements with special fundamental solutions for analyzing plane elastic problems containing holes. In Ref. [16], Qin developed the same elements for three-dimensional elastic problem. In Ref. [17], Qin classified his developed fundamental solution-based finite elements as a fourth type of the finite element method. He

T.H.A. Naga, Y.F. Rashed / Engineering Analysis with Boundary Elements 71 (2016) 70–78

71

demonstrated the main advantage of his new elements compared to the traditional finite element methods. It has to be noted that in all of the previous work [13–17], the so-called “regular collocation” or placing to collocation point outside the element domain is considered to avoid computing singular integrals. Such a technique has proven to be efficient. The purpose of this paper is to derive “Hybrid Boundary solution” procedures for shell elements (or fundamental solutionbased shell elements). The formulation for the uncoupled 2D membrane and plate bending will be based on the extended Tong hybrid variational functional. Fundamental solutions are employed for both 2D membrane and plate in bending. Quadratic boundary elements are used to represent the boundaries of the new developed shell elements. In this work, regular collocation is also used to avoid singularities; especially in the 3D spatial space for generally inclined shell element. The derived element is going to be tested via numerical examples. Fig. 2. The proposed shell finite element.

2. The proposed new finite element for plate bending Without losing the generality, the proposed element could be formed from any multi-line/curve closed polygon; however, in this paper a traditional eight-node quadrilateral shape is considered. Consider the new plate finite element of domain Ω(y) with eight nodes and 24 degrees of freedom bounded by a closed boundary surface Γ as shown in Fig. 1. A generalized form of the hybrid energy functional ( Π ) of the type proposed by Tong [4] for Reissner's plate bending could be written as follows:

(

Π ui ( y), u˜ i ( x), p˜ i ( x)

Simply Supported Thick. 250mm

A

A

Thick. 500mm

)

1 Mαβ ( y) χαβ ( y) + Q 3α ( y) ψ3α ( y) dΩ ( y) Ω ( y) 2

(

)

=

∫

−

∫Ω( y) bi ( y) ui ( y) dΩ ( y) − ∫Γ ( x) p¯ i ( x) u˜ i ( x) d Γ ( x) P

−

2.50m

∫Γ ( x,y) p˜ i ( x) ( u˜ i ( x) − ui ( y) ) d Γ ( x, y)

(1)

2.50m

Fig. 3. The Multi-thickness circular plate in Example (1).

Where Mαβ ( y) , Q 3α ( y) , bi ( y) , ui ( y) are the moment, shear stress resultants, body loads and generalized displacements respectively. The symbols χαβ ( y) andψ3α ( y) denote bending and shear strain. ΓP is the portion of the boundary where generalized tractions are prescribed. The over bars ( p̅ i ( x)) indicate prescribed boundary values for tractions and the over tilde ( ũ i ( x) andp̃ i ( x)) represents a displacement value or traction value which is defined only on the

ξ6

ξ7 x

ξ5

6

7

5

Γ(x,y) ξ8

8

Ω(y)

3

1

ξ1

ξ4

4

2

ξ2 Fig. 1. The proposed finite element for plate bending.

ξ3 Fig. 4. Mesh (A) in Example (1), 80 elements with 8 nodes each.

72

T.H.A. Naga, Y.F. Rashed / Engineering Analysis with Boundary Elements 71 (2016) 70–78

2.1. Approximation of the domain variables The first four integrals in Eq. (2) contain the unknown domain variables ( ui ( y) , pi ( y)). which are approximated as the summation of product of fundamental solution ( U* ( y , ξ )) [18] and an unki

n

known set of fictitious concentrated generalized tractions ( γk ( ξn )) located at a set of arbitrary source points ( ξn ); where n¼ 1 to 8 (recall Fig. (1)) [19], as follows: n= 8

ui ( y)=

∑ U*ki ( y, ξ n ) γk ( ξ n )

(3)

n= 1

Where γk ( ξn )is the set of fictitious generalized tractions applied along the direction (xk). In this paper, (n) is taken equal to 8 to represent the same number of the boundary points and placed along offset equal to the element length as shown in Fig. (1). In a similar way, the generalized traction components can be written in matrix form as follows: n= 8

pi ( y)=

∑ P*ki ( y, ξ n ) γk ( ξ n )

(4)

n= 1

(

)

Where P*ki ( y,ξn ) is the traction fundamental solution kernel [18]. Fig. 5. Mesh (B) in Example (1), 192 elements with 8 nodes eac.

2.2. Approximation of the boundary variables

problem boundary. Integrating by parts the first domain integral on the right hand side of Eq. (1), it gives [18]:

(

1 ui ( y) pi ( y) d Γ ( y)− 2

) ∫

Π ui ( y), ũ i ( x), p̃ i ( x) =

Γ ( y)

uα ( y) dΩ ( y)−

∫ Ω ( y)

∫

Ω ( y)

1 M ( y) 2 αβ, β

m=3

ũ i ( x)=

1 Q ( y) u3 ( y) dΩ ( y)− 2 3α, α

∫

bi ( y) ui ( y) dΩ ( y)−

Ω ( y)

Γ+

∫

The unknown generalized boundary displacement and traction vectors denoted by ũ i ( x) andp̃ i ( x) are approximated using quadratic shape function ϕm ( η), as follows:

m=3

p̃ i ( x)=

Γ P ( x)

∫ Γ ( x)

m

ϕ m ( η) ui ( x e) ∀ xinΓe

(5)

m=1

ũ i ( x) p̅ i ( x) d

∑

m

ϕ m ( η) pi ( x e) ∀ xinΓe

(6)

m=1

∫

p̃ i ( x) ũ i ( x) d Γ ( x)−

∑

p̃ i ( x)

m m ui ( x e) andpi (x e)

Where, are vectors of the nodal values of the element ( x e) for boundary displacements and boundary tractions.

Γ ( x, y)

ui ( y) d Γ ( x, y)

(2) 2.3. The proposed stiffness equation

The main idea of the present formulation is to approximate the field variables ( ui ( y) , pi ( y) , ũ i ( x) , p̃ i ( x)) before applying the variational operator to the functional in Eq. (2).

Using the representation given in Eqs. (3)–(6), Eq. (2) could be re-written as follows:

Distance along the strip (m)

0 -5.5

-3.5

-1.5

0.5

2.5

-0.001 -0.0015 -0.002

Deflection (m)

-0.0005

-0.0025 -0.003 -0.0035

Analytical [21] FEM Mesh (A)

-0.004

FEM Mesh (B) Present Formulation Mesh (A) Ref. [21] (441 division)

-0.0045

-0.005

Fig. 6. Deflections along section (A-A) in Example (1).

4.5

T.H.A. Naga, Y.F. Rashed / Engineering Analysis with Boundary Elements 71 (2016) 70–78

73

0 -5.5

-4.5

-3.5

-2.5

-1.5

-0.5 -2 -4

Bending moment(t.m)

Distance along the strip (m) 2

0.5

1.5

2.5

3.5

4.5

5.5

-6 -8 Analytical [21] -10

FEM Mesh (A) FEM Mesh (B)

-12 Present Formulation Mesh (A) Ref. [21] (441 division)

-14 -16

Fig. 7. Bending moment along section (A-A) in Example (1).

Π3 =

T T 1 { γ}24×1 ⎡⎣ F ⎤⎦24×24 { γ}24×1– { u}24×1 P¯ 2

{ }24×1

T

T

T

+ { p}24 × 1 ⎡⎣ L ⎤⎦24 × 24 { u}24 × 1 − { p}24 × 1 ⎡⎣ G⎤⎦ { γ}24×1 24 × 24 T

− { γ}24 × 1 { B}24 × 1

(7)

Where

⎡⎣ F ⎤⎦ = 24 × 24 T

⎡⎣ G⎤⎦ = 24 × 24

∫Γ( y) U*ki ( y, ξn ) P*mi ( y, ξn ) d Γ( y)

∫Γ ( x ) ϕ j( η) U*ki ( x e, ξn ) d Γ( x e) e

(9)

e

{ P¯ }24×1 = ∫

ϕ j ( η) p¯ i (x e) d Γ(x e) (10)

Γe ( x e)

⎡⎣ L ⎤⎦ = 24 × 24 Fig. 8. The multi–thickness cantilever plate in Example (2).

(8)

∫Γ ( x ) ϕ j ( η) ϕ j ( η)T d Γ( x e)

(11)

{ B}24×1 = ∫Ω( y) U*ki ( y, ξ n ) bi ( y) dΩ( y)

(12)

e

e

At equilibrium conditions, the functional Π3 is stationary i.e. its Distance along the strip (m) 0 9

8

7

6

5

4

3

2

1

0

-0.5

Analytical [21]

-1 FEM Mesh (A) FEM Mesh (B) Ref. [21] (5x20 Divisions) Ref. [21] (5x25 Divisions)

-1.5

Ref. [22] Proposed Formulation Mesh (A)

-2

-2.5

Fig. 9. Deflection of the multi-thickness plate along the center line, in Example (2).

Deflection (m)

10

74

T.H.A. Naga, Y.F. Rashed / Engineering Analysis with Boundary Elements 71 (2016) 70–78

120

Analytical [21]

100

FEM Mesh (A) FEM Mesh (B)

Ref. [21] (5x25 Divisions)

60 Proposed Formulation Mesh (A)

Bending moment(t.m)

80

Ref. [21] (5x20 Divisions)

40

20

0 10

9

8

7

6

5

4

3

2

1

0

Distance along the strip (m)

Fig. 10. Bending moment of the multi-thickness slab along the center line, in Example (2).

A q A

Fixed support Fig. 11. The L-shaped plate structure in Example (3).

1.20 1.10 1.00

Fig. 13. Bending moment of L-shaped plate structure along section (A-A), in Example (3).

0.90 0.80

Fixed support

Analytical [22]

0.70

FEM Mesh (A) FEM Mesh (B)

0.60

Proposed Formulation Mesh (A) Undeformed shape

0.50

Ref. [22]

P1

0.40 0.30

P2

0.20 0.10 0.00 -0.10

0.00

0.10

0.20

0.30

0.40

0.50

0.60

0.70

0.80

0.90

1.00

1.10

Distance along the strip (m) Fig. 12. Deformed shape of L-shaped plate structure along section (A-A), in Example (3).

first variation δΠ3 vanishes for any arbitrary values of(δγ, δu and δp). Therefore the corresponding generalized Euler-Lagrange equations are:

Fig. 14. The hollow cantilever box considered in Example (4).

⎡⎣ F ⎤⎦ { γ}24 × 1 − ⎡⎣ G⎤⎦24 × 24 { p} − { B}24 × 1 = 0 24 × 24 T

⎡⎣ L⎤⎦ { p}24 × 1 − P¯ 24 × 24

{ }24×1 = 0

(13)

(14)

T.H.A. Naga, Y.F. Rashed / Engineering Analysis with Boundary Elements 71 (2016) 70–78

75

Fig. 15. Vertical displacements of the box in Example (4) due to bending loading.

0.09

0.07

0.05

0.03

z/L2

0.01 -0.25

-0.2

-0.15

-0.1

-0.05

-0.01

0

0.05

0.1

0.15

0.2

0.25

-0.03 FEM mesh (A) FEM mesh (B)

-0.05

FEM mesh (C) Proposed Formulaon mesh (A) Proposed Formulaon mesh (B)

-0.07

Ref. [22] mesh (D) -0.09

x/L2

Fig. 16. Deformed shape of cross section at y ¼2 in the box of Example (4) due to torsional loading.

T

⎡⎣ L ⎤⎦ { u}24 × 1 − ⎡⎣ G⎤⎦ { γ}24 × 1 = 0 24 × 24 24 × 24

(15)

Eqs. (13)–(15) can be written in terms of the boundary displacements only as follows:

⎡⎣ K ⎤⎦ { u}24 × 1 = { Q}24 × 1 24 × 24

(16)

Where:

⎡⎣ K ⎤⎦ ⎡⎣ F⎤⎦ ⎡⎣ ⎤⎦ = ⎡⎣ R⎤⎦ 24 × 24 24 × 24 R 24 × 24 24 × 24 T ⎡ ⎤−1 ⎡⎣ R ⎤⎦ = ⎢ ⎡⎣ G⎤⎦ ⎥ ⎡⎣ L ⎤⎦24 × 24 24 × 24 24 × 24 ⎦ ⎣

(17)

(18)

In which, { u}24 × 1 is the generalized displacement vector, and given by: T

And T

{ Q}24×1

⎡⎣ Q⎤⎦ = 24 × 1

(20)

is the force vector, and could be written as:

{ { Q} { Q} { Q} i1

i2

i3

⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅⋅{ Q}i8

}

(21)

Where { Q}ij j = 1 to 8 is the component of the force, postulated at each degree of freedom, for example:

T

⎡⎣ u⎤⎦ = 24 × 1

{ u}i4 = { θ14 θ24 u34 }

{

{ u}i1 { u}i2 { u}i3 ⋯⋯⋯⋯⋯{ u}i8 }1 × 8

(19)

Where { u}ij ,j = 1to8is the component of the generalized displacement, postulated at each degree of freedom, for example:

{ Q}i4 = { M14 M24 Q 34}

(22)

It has been noted that the obtained [K] or the stiffness matrix is symmetric, positive definite and similar to the one obtained from the finite element method [7]. As previously mentioned the source points ( ξn ) are located outside the proposed element domain; therefore the integrals involved in computing values of matrices [F], [G], [ P̅ ], [L] and [B] are regular, and could be computed numerically using the standard Gauss quadrature scheme. It has to be noted that the domain integral { B} is evaluated using equivalent boundary integrals as given in Ref. [18].

76

T.H.A. Naga, Y.F. Rashed / Engineering Analysis with Boundary Elements 71 (2016) 70–78

⎡⎣ K̅ b⎤⎦ { u̅ b }24×1= Q̅ b 24 × 24

{ }24×1

(24)

Where, the subscripts “p” and “b” denote quantities for 2D membrane and for plate bending, respectively. Combining { u̅ p} and { u̅ b} to give:

{ u̅ }={ u1u2 u3θ1θ2θ3} And combining

(25)

{ Q̅ } and { Q̅ } to give: p

b

{ Q̅ }={ Q1Q 2Q 3M1M2 M3}

(26)

The relationship between the nodal local boundary generalized displacement vector to the nodal boundary generalized force vector of the shell can be written as follows:

⎡⎣ K̅ ⎤⎦ { u̅ }48×1= Q̅ 48 × 48

{ }48×1

(27)

Where ⎡⎣ K̅ ⎤⎦ is the shell stiffness in the local coordinate sys48 × 48 tem; which can be constructed as follows:

⎡ ⎡ K̅ ⎤ ⎢ ⎣ p⎦ ⎡⎣ K̅ ⎤⎦ =⎢ ⎡⎣ 0⎤⎦ 48 × 48 ⎢ ⎢ ⎡⎣ 0⎤⎦ ⎣ Fig. 17. The cylindrical barrel roof considered in Example (5).

3. The proposed shell finite element The previously derived stiffness matrix for the plate bending in Section (2) is coupled to the previous work for 2D membrane formulation in Ref. [9] to produce shell element. Formulation of generally inclined element is presented in this section. Fig. (2) demonstrates general configuration of the proposed shell finite element. 3.1. Local stiffness matrix of the proposed shell finite element The local stiffness matrix ⎡⎣ K̅ b⎤⎦ and ⎡⎣ K̅ p⎤⎦ of the plate bending and 2D membrane respectively are initially constructed from the previously stated formulation and from Ref. [9], relating the nodal boundary generalized displacement vector to the nodal boundary generalized force vector by the following relationships:

⎡⎣ K̅ p ⎤⎦ { u̅ p }16×1= Q̅ p 16 × 16

{ }

(23)

16 × 1

⎡⎣ 0⎤⎦ ⎡⎣ 0⎤⎦⎤ ⎥ ⎡⎣ K̅ b⎤⎦ ⎡⎣ 0⎤⎦⎥ ⎥ ⎡⎣ 0⎤⎦ K̅ θ3 ⎥⎦ 48 × 48

(28)

It can be seen that ⎡⎣ K̅ ⎤⎦ has zero in-plane rotational stiff48 × 48 ness. The value of K̅ θ3is set equal to 10  7 times the smallest bending stiffness [20]. This value is chosen to remove the in-plane rotational singularity from the element stiffness matrix when the local x1, x2;x3 axes coincide with the global X-Y-Z axes [20]. 3.2. Transformation of the stiffness matrix The stiffness matrix ⎡⎣ K̅ ⎤⎦ obtained in Eq. (28) represents the stiffness in the local coordinate system x1, x2, x3. Before the assembly of the over-all stiffness matrix of the structure, the shell stiffness matrix should be transformed to the global coordinate system X, Y, Z. The relationship between generalized displacement and generalized force vector in the global and local coordinates system could be written as follows

{ u̅ }48×1=⎡⎣ T⎤⎦48×48 { u}48×1

(29)

{ Q̅ }48×1=⎡⎣ T⎤⎦48×48 { Q}48×1

(30)

0.1

0.05

Angle along line (A-B) 0 0

5

10

15

20

25

Deflection (ft)

-0.05

-0.1

-0.15

Analytical [1] OFFSET =0.40

-0.2

OFFSET =0.50 to 8.50 OFFSET =8.75

-0.25

-0.3

-0.35

Fig. 18. Effect of source point location on results.

30

35

40

T.H.A. Naga, Y.F. Rashed / Engineering Analysis with Boundary Elements 71 (2016) 70–78

Where ⎡⎣ T⎤⎦48 × 48 is the transformation matrix (also known as the direction cosines matrix) [20]. Finally, the stiffness matrix ⎣⎡ K⎦⎤

(31)

distributed load (q) of intensity 1t/m as shown in Fig. (11). (Figs. (12) and 13) demonstrate the deformed shape and the bending moment along section (A-A). Results are plotted from the proposed new finite element together with those obtained from the traditional finite element method [2] and the conventional direct boundary element method [22]. The results are compared also to the analytical values. Three meshes are considered:

(32)

 Mesh (A): consist of one element with 8 nodes for each surface.  Mesh (B): consist of 16 elements with 4 nodes.  Mesh (C): consist of one region with 12 boundary elements for

48 × 48

of the shell element with reference to the global coordinate system can be obtained by substituting Eqs. (29) and (30) into Eq. (27), to give:

⎡⎣ K⎤⎦ { u}48 × 1={ Q}48 × 1 48 × 48 where T ⎡⎣ K̅ ⎤⎦ ⎡⎣ T⎤⎦ ⎡⎣ K⎤⎦ =⎡⎣ T⎤⎦ 48 × 48 48 × 48 48 × 48 48 × 48

77

each surface. 4. Numerical examples In this section, several numerical examples are presented to demonstrate the accuracy and the validity of the proposed finite element for the plate bending problems and for shell problems. 4.1. Example 1: multi- thickness simply supported circular plate A multi-thickness circular plate are simply supported along its external boundary, as shown in Fig. (3), is considered in this example. The plate is subjected to a uniformly distributed load with intensity 1 t/m2. The Young's modulus for the plate material is E¼ 3  106 t/m2 and Poisson's ratio ν is 0.20. Two meshes (A and B) are used, as shown in (Figs. (4) and 5). (Figs. (6) and 7) demonstrate the deflection and bending moment along section (A-A) of plate using meshes (A) and (B). Results are plotted from the proposed new finite element together with analytical values; results obtained from the traditional finite element method [2] and the stiffness boundary element formulation in Ref. [21]. It can be seen that, the present formulation is accurate compared to analytical solution even when using the coarse mesh (A). 4.2. Example 2: multi-thickness cantilever plate The cantilever multi-thickness plate, as shown in Fig. (8), is considered in this example. The plate is subjected to a uniform distributed load with intensity  1 t/m2. The Young's modulus for the plate material is E ¼1  105 t/m2 and Poisson's ratio ν is set to zero. (Figs. (9) and 10) demonstrate the deflection and bending moment along the plate center line. Results are plotted from the proposed new finite element together with the traditional finite element method [2] and those of two boundary element strategies: the first is the traditional boundary element method with sub regions [22] and the second is based on stiffness formulation of in Ref. [21]. The results are also compared with the analytical values. Three meshes are considered:

 Mesh (A): consist of two elements with 8 nodes.  Mesh (B): consist of 64, four-node quadrilateral elements.  Mesh (C): consist of two sub regions with 12 boundary element. (Figs. (9) and 10) demonstrate the superiority of the present element in term of accuracy; even with two elements. 4.3. Example 3: L-shaped shell structure The previous two examples verified the present formulation for plate bending problems. This example and the next will test the present element for solution of shell structures. The L-shaped plate shown in Fig. (11) with thickness 0.1 m, Young's modulus is E¼1  105 t/m2 and Poisson's ratio ν is 0.0 is considered in this example. The structure is subjected to

(Figs. (12) and 13) demonstrate the accuracy of the present formulation. 4.4. Example 4: hollow cantilever box The hollow cantilever box of span (L2 ¼2000 mm), shown in Fig. (14), and of thickness of 2mm is considered in this example. The box is subjected to different load cases: 1. Case (1): bending loading, in which P1 ¼P2 ¼ 5000N 2. Case (2): torsion loading, in which P1 ¼ – P2 ¼ 5000N The Young's modulus for the plate material is E¼ 7.0  104 MPa and Poisson's ratio ν is 0.30. (Figs. (15) and 16) demonstrate the vertical displacement and the deformed shape at y ¼2 m for the bending and torsion cases of loading, respectively. Results are plotted from the proposed new finite element together with results of the traditional finite element method [2] and the conventional direct boundary element method with sub regions [22]. The results are also compared to the analytical values. Five meshes are used:

    

Mesh (A): consist of 80 elements with 8 nodes. Mesh (B): consist of 160 elements with 8 nodes. Mesh (C): consist of 320, four-node quadrilateral elements. Mesh (D): consist of one sub-region with 16 boundary element for each region. Mesh (E): consist of 4608, four-node quadrilateral elements.

Fig. (15) demonstrates that the results of the proposed formulation for mesh (A) and mesh (B) are close to the analytical solution. From Fig. (16), it can be seen that, the result of the proposed formulation for both meshes (A) and (B) are close to results of the finest finite element mesh (C). 4.5. Example 5: Effect of source point location The previous examples verified the present formulation for plate bending and shell problems. This example demonstrates the effect of source point location on the accuracy of the results. The cylindrical barrel roof shown in Fig. (17) is considered. It has a thickness of 0.25 ft, Young's modulus of E ¼4.3  108 Ib/ft2 and Poisson's ratio ν is set to zero. The roof is loaded by a distributed load of intensity 90 Ib/ft2. Only one quadrant is modeled with a 6  12 shell elements. Fig. (18) demonstrates the variation of the deflection along line A-B for different source point offset from (0.40 to 8.75) times the used element side length. It can be seen that the result are accurate when the source point offset is within the range (0. 5 to 8.5) times the used element side length. 5. Conclusions In this paper, a new finite element for shell structures was

78

T.H.A. Naga, Y.F. Rashed / Engineering Analysis with Boundary Elements 71 (2016) 70–78

developed using an improved hybrid boundary solution. The formulation uses the fundamental solutions for 2D plane and plate bending problems as trial functions. No singular integrals was involved as the place of source points were chosen outside the element domain. It was demonstrated that the choice of source point locations has almost no effect on the results within a recommended range (0. 5 to 8.5) times the used element side length. Several examples with different boundary conditions and geometries were tested. It was demonstrated that the present formulation results were very accurate even with small number of elements. The present formulation combines both advantages of boundary and finite element methods.

References [1] Timoshenko s, Woinowsky-Krieger S. Theory of plates and shells. 2nd ed.. New York: McGraw-Hill; 1959 1959. [2] SAP2000 V16. Structural analysis program software, integrated finite element analysis and design of structures, computer and structures, Inc. Berkeley, California: University Avenue; 2013. [3] ANSYSs Academic Research, Release 16.0, Help System, Coupled Field Analysis Guide, ANSYS, Inc. [4] Pian THH, Tong P. Basis of finite element method for solid continua. Int J Numer Methods Eng 1969;1:3–28. [5] Rashed YF, Mobasher ME. Products and practice spotlight: a new tool for structural designers. ACI Concr Int 2012;34:54–5. [6] Rashed YF, 2000, Boundary Element Formulations for Thick Plates, Topics in

Engineering, Vol. 35., WIT Press, Southampton, UK. [7] Zienkiewicz OC. The finite element method.NewYork: McGraw-Hill; 1977. [8] Dumont NA. The hybrid boundary element method: an alliance between mechanical consistency and simplicity. ASME Appl Mech Rev 1989;42:S54–63. [9] DeFigueiredo, TGB, 1991. A New Boundary Element Formulation in Engineering (Lecture Notes in Engineering, Vol. 68), Springer, New York. [10] Leung KL, Zavareh PB, Beskos DE. 2-D elastostatic analysis by a symmetric BEM/FEM scheme. Eng Anal Bound Elem 1995;15(1):67–78. [11] Dumont NA. Variationally-based, hybrid boundary element methods. Comp Assist Mech Eng Sci 2003;10:407–30. [12] Naga THA, Rashed YF. A variational boundary element formulation for sheardeformable plate bending problems. J Appl Mech 2013:80. [13] Wang H, Qin QH. Hybrid FEM with fundamental solutions as trial functions for heat conduction simulation. Acta Mech Solida Sin 2009;22(5):487–98. [14] Wang H, Qin QH. Fundamental-solution-based finite element model for plane orthotropic elastic bodies. Eur J Mech-A/Solids 2010;29(5):801–9. [15] Wang H, Qin QH. Fundamental-solution-based hybrid FEM for plane elasticity with special elements. Comput Mech 2011;48(5):515–28. [16] Cao C, Qin QH, Yu A. A new hybrid finite element approach for three-dimensional elastic problems. Arch Mech 2012;64:261–92. [17] Qin QH. Fundamental solution based finite element method. J Appl Mech Eng 2013;2:e118. http://dx.doi.org/10.4172/2168-9873.1000e118. [18] Transformation of Domain Effects to the Boundary. In: Rashed YF, Brebbia CA, editors. Southampton, UK: WIT Press; 2003. [19] Mohareb SW, Rashed YF. A dipole method of fundamental solutions applied to Reissner's plate bending theory. Mech Res Commun 2009;36:939–48. [20] Bathe KJ, Ho LW. A simple and effective element for analysis of general shell structures. Comput Struct 1981;13:673–81. [21] Wagdy M, Rashed YF. Boundary element analysis of multi-thickness sheardeformable slabs without sub-regions. Eng Analy Bound Elem 2014;43:86–94. [22] Dirgantara T, Aliabadi MH. Boundary element analysis of assembled plate structures. Commun Numer Methods Eng 2001;17(10):749–60.