Mesh Free Methods for Shells

Previous chapters have introduced a number of MFree methods for solids, fluids, beams,
and plates. This chapter formulates MFree methods for shell structures.

Spatial thin shell structures are used very extensively in many engineering structures,

including aircraft, pressure vessels, storage tanks, and so on, due to their outstanding
efficiency in utilizing materials. Many shells are designed using advanced computer-aided
design (CAD) technology to be the main load-carrying structure. Advanced CAD design
of shell structures must consider both static and dynamic loads. Hence, both static and
frequency analysis are very important and are integral parts of the shell design. Because of
the complex nature, both structurally and mechanically, numerical means have to be utilized
for analyses of shells. Advances in computing technology have made this possible and
have led to the widespread use of numerical analysis of complex shells. The finite element
method (FEM) remains the most popular numerical technique as evident from numerous
publications since the 1970s (Yang et al., 1990). However, FEM requires the design of
meshes, which is an extremely tedious, costly, and time-consuming process. In addition,
meshing a shell structure into elements carries a number of risks that are often hidden by
colorful pictures, as the curved surface of the shell is usually modeled by shell elements
that are flat. This kind of geometric simplification in FEM can lead to severely erroneous
results, due to the mechanics coupling effects of bending forces and membrane forces.

MFree methods present a promising alternative to FEM, as they offer opportunities to

relieve the manual meshing process in modeling a structure. This is particularly important
for shells, as shell structures are very complex both in field variable variation and in
geometric configuration. MFree methods can also offer a very important additional capa-
bility in representing the complex curved geometry of shell structures using MFree approx-
imation techniques. These MFree approximations both in field variables and in the
structure itself provide much more accurate results compared with conventional FEM.

Because MFree methods are still very young, very few works have been reported in the

development of MFree methods for shell structures. The first contribution in this regard
was made by Krysl and Belytschko (1996b) based on the thin shell theory using moving
least square (MLS) approximation with Lagrange multipliers for essential boundary con-
ditions. Noguchi et al. (2000) developed a formulation for thick shell using MLS approx-
imation with penalty method for handling essential boundary conditions. Li et al. (2000)
formulated an MFree method based on the reproducing kernel particle method for thin
shells with large deformation. In this work, the essential boundary condition is imposed
by modifying shape functions for nodes near and on the essential boundaries. These three
papers discuss only static deformation problems.

This chapter formulates MFree methods for shell structures based on the element free

Galerkin (EFG) method introduced in

Chapter 6,

and the point interpolation method (PIM)

discussed in

Chapter 8

for two-dimensional (2D) solids. The work has been performed

by L. Liu, G. R. Liu, and V. B. C. Tan during the past 2 years. This chapter covers the
following topics:

1238_Frame_C12.fm Page 501 Wednesday, June 12, 2002 5:24 PM

• MFree formulations for both thin and thick shells based on Galerkin formulation
• Use of both MLS and RPIM shape functions
• Both static and dynamic problems
• Approximation of both field variables and the geometry of the shell using MFree

shape functions

12.1 EFG Method for Spatial Thin Shells

This section formulates an EFG method for thin shells governed by Kirchhoff–Love shell
theory. In the EFG method, the generalized displacement (deflections and rotations) at an
arbitrary point is approximated from nodal displacements using MLS approximation. A
compact support domain is used to determine field nodes to be used for constructing MLS
shape functions. As discussed in

Chapter 6,

use of MLS approximation requires special

treatment for essential boundary conditions for shell structures because EFG interpolations
lack the Kronecker delta property. These techniques include the penalty method (Noguchi
et al., 2000), Lagrange multipliers (Krysl and Belytschko, 1996b), and a method that
modifies shape functions for nodes near and on the essential boundaries (Li et al., 2000).
This chapter discusses both the penalty method and Lagrange multipliers method for
analyzing static problems. For dynamic analysis, the Lagrange multipliers method is used
for transient analyses and the matrix transformation method is used for free-vibration
analyses, as discussed in Chapters 10 and 11 for beams and plates.

The formulation presented in this section is based on the work by L. Liu et al. (2001).

It begins with a brief discussion of MLS approximation, which is basically the same as
that provided in

Chapter 5,

except that a higher order of polynomial basis needs to be

included. The governing equations for the analysis of general shells and membrane struc-
tures are then introduced. Numerical formulation based on a geometrically exact theory
accounting for the Kirchhoff hypothesis is presented. This is followed by the definition of
curved surfaces, kinematics of shells, stress and strain measures, and the constitutive
relations adopted in the formulation. The final discrete equations for the entire shell
structure for static, free vibration, and transient vibration then are obtained. For free
vibration, the essential boundary conditions are imposed using orthogonal transform
techniques to solve the eigenvalue equation (Ouatouati and Johnson, 1999; G. R. Liu and
Chen, 2001). For static problems, essential boundary conditions are imposed through the
Lagrange multipliers method and the penalty method. Finally, the method is applied to
several numerical examples of shells of different geometries to illustrate the efficiency and
accuracy of the present EFG method.

12.1.1

Moving Least Squares Approximation

The derivation of shape functions from the MLS approximation method is the same as
that provided in Chapter 5, except that a higher order of polynomial basis needs to be
included. A 2D field approximation is needed for modeling thin shells. A component of
the displacement vector is approximated by a polynomial function as follows:

(12.1)

( )

( )u

∑

1238_Frame_C12.fm Page 502 Wednesday, June 12, 2002 5:24 PM

where shape function

(

) is the MLS shape function obtained following the procedure

described in Section 5.4 using polynomial basis. The formulations for computing the
derivatives of the shape functions are given in Equations 5.62 through 5.67. Note that for
shell structures the order of the polynomial basis should be higher than that for 2D solids.
In this chapter, the two different orders of polynomial basis are primarily used. The
following six terms of basis functions up to quadratic terms are used for shells where
shear effect is significant.

(12.2)

The following 15 terms of basis functions up to quartic terms are used to ameliorate
membrane locking in bending-dominated cases:

(12.3)

In computing the MLS shape functions, the quartic spline weight function given in

Equation 5.12 is used in this section because of requirements on the continuity of the MLS
shape functions and their derivatives.

12.1.2

Governing Equation for Thin Shell

The shells considered in this section are assumed to be thin with arbitrarily deep and
Gaussian curvature governed by Kirchhoff–Love theory. The governing equations used in
this section are based on geometrically exact theory of shells formulated by Simo and Fox
(1989) with some modifications to account for the Kirchhoff hypothesis. Here, we outline
only the basic concepts of the formulation. Details can be found in the paper by Simo and
Fox (1989).

Kinematic Description of Shell

The reference frame coordinates are illustrated in Figure 12.1. The shell in 3D space is
described in a global Cartesian coordinate system,

. The Gauss intrinsic coordinates

defined locally are used to describe the configuration of the shell.

ϕϕϕϕ

(

) gives the

position of the point on the shell neutral surface, and

(

) is a direction unit vector

FIGURE 12.1

Reference frames of coordinates on the neutral surface of a thin shell.

( )

1, x, y, x

, xy, y

{

}

( )

1, x, y, x

, xy, y

, x

y, xy

, y

, x

y, x

, xy

, y

{

}

1238_Frame_C12.fm Page 503 Wednesday, June 12, 2002 5:24 PM

normal to the shell neutral surface both in the unformed reference and deformed states
according to the Kirchhoff–Love hypothesis. The pair (

ϕϕϕϕ

) defines the position of an

arbitrary point in the shell. The configuration of the shell can be expressed mathematically
as

(12.4)

Here

denotes the parametric space for the shell; (

−

) are the distances of the “lower”

and “upper” surfaces of the shell measured from the shell neutral surface.

The convective basis vectors

are defined as

(12.5)

where denotes the gradient operator,

denotes tensor product, and

is the unit basis

vector in the global Cartesian basis. A contravariant basis

can be obtained from the

standard relation of

(12.6)

The unit area in the neutral surface is defined by the differential two-form:

(12.7)

where “

” denotes vector cross product. The determinants of the tangent maps in the

deformed and reference configuration will be denoted subsequently as

and

, respectively,

with and

denoting the Jacobians on the reference surface

(12.8)

where the superscript “0” is used to denote quantities in the reference configuration. The
surface convected frame, which spans the tangent space to the neutral surface, is defined
as

ϕϕϕϕ

(

= 1, 2). Hence, the first fundamental form on the reference surface is

(12.9)

where

(

1, 2) denotes the dual surface convected basis through the standard relation

⋅

and

⋅

denotes dot product. The second fundamental form is defined as

(12.10)

Strain Measures

The linear membrane and bending strain measures can be derived from the kinematic
variables in Equations 12.9 and 12.10 as

(12.11)

ϕϕϕϕ

(

)

ξ t ξ

(

) with

and

−

, h

〈

〉

∈

{

}

∇x

∂

-------

≡

∇

⋅

ϕϕϕϕ

det

∇x

[

], j

det

∇x

[

], j

ξ =0

, j

ξ =0

αβ

ϕϕϕϕ

⋅

αβ

ϕϕϕϕ

⋅

αβ

1
2

---

ϕϕϕϕ

⋅

(

)

1238_Frame_C12.fm Page 504 Wednesday, June 12, 2002 5:24 PM

and

(12.12)

where only the symmetric part of the bending strain measure is considered.

The Kirchhoff–Love hypothesis needs to be introduced explicitly to obtain the definite

forms for the strain measures. The mathematical form of this hypothesis is expressed as

(12.13)

where

denotes the norm of a vector. Hence, the derivatives of the normal vector in

the reference configuration t

and partial derivatives of the increment

∆t can be

(12.14)

(12.15)

Note that the membrane strain measures of Equation 12.11 are not affected by the

introduction of the Kirchhoff–Love hypothesis. Considering the symmetry with respect
to partial differentiation

and u

,12

= u

,21

, the bending strain measures can be

rewritten as

(12.16)

(12.17)

(12.18)

Stress Resultants and Stress Couples
A section in the current configuration is described by

(12.19)

The stress resultants and resultant couples are defined by normalizing the force and

torque with the surface Jacobian

as follows

(12.20)

(12.21)

αβ

1
2

---

ϕϕϕϕ

⋅

∆

⋅

∆

⋅

(

)

( )

–

ϕϕϕϕ

(

( )

–

ϕϕϕϕ

(

)

∆t

( )

–

ϕϕϕϕ

(

)

ϕϕϕϕ

,12

ϕϕϕϕ

,21

,11

( )

–

ϕϕϕϕ

,11

ϕϕϕϕ

(

) u

ϕϕϕϕ

,11

(

)

⋅

[

]

⋅

–

,22

( )

–

ϕϕϕϕ

,22

ϕϕϕϕ

(

) u

ϕϕϕϕ

,22

(

)

⋅

[

]

⋅

–

1
2

---

,12

,21

(

) t

1
2

--- j

( )

–

ϕϕϕϕ

,12

ϕϕϕϕ

,21

(

) ϕϕϕϕ

(

) u

ϕϕϕϕ

,12

ϕϕϕϕ

,21

(

)

(

)

⋅

[

]

⋅

–

,12

( )

–

ϕϕϕϕ

,12

ϕϕϕϕ

(

) u

ϕϕϕϕ

,12

(

)

⋅

[

]

⋅

–

=const

∈













1, 2

ϕϕϕϕ

( )

–

σg

ξ, α

−

∫

1, 2

( )

–

ϕϕϕϕ

–

(

)

σg

ξ, α

−

∫

1, 2

1238_Frame_C12.fm Page 505 Wednesday, June 12, 2002 5:24 PM

The director stress couple

can also be defined through the expression

(12.22)

The through-thickness stress resultant has been omitted because it does not play a role

in Kirchhoff–Love theory.

Constitutive Equations
For the isotropic elastic shell structures, the effective membrane and stress couple resultant
for the isotropic hyperelastic material can be written as

(12.23)

(12.24)

where the matrix C is given by

(12.25)

Here a

αβ

are the components of the first fundamental form in the dual basis.

12.1.3

Strain–Displacement Relations

The displacement vector can be expressed in the global Cartesian basis E

(12.26)

where n is the number of nodes in the neighborhood (support domain) of

ζ; u

, v

, and w

are the components of the displacement vector at the Ith node in E

, E

, and E

directions,

respectively. The membrane strain–displacement relation for the Ith node is obtained by
substituting the displacement approximation Equation 12.26 into Equation 12.11 to give

(12.27)

⇒

( )

–

ξσg

−

∫

n˜

















–

--------------

































12 1

–

(

)

--------------------------

















011

(

)

νa

011

022

–

(

) a

012

(

)

(

)

011

012

022

(

)

022

012

symm

1
2

(

) a

012

(

)

–

(

011

022

(

)

( )

( ) u

[

]

∑

















I,1

ϕϕϕϕ

⋅

I,1

ϕϕϕϕ

⋅

I,1

ϕϕϕϕ

⋅

I,2

ϕϕϕϕ

⋅

I,2

ϕϕϕϕ

⋅

I,2

ϕϕϕϕ

⋅

I,1

ϕϕϕϕ

I,2

ϕϕϕϕ

(

) E

⋅

I,1

ϕϕϕϕ

I,2

ϕϕϕϕ

(

) E

⋅

I,1

ϕϕϕϕ

I,2

ϕϕϕϕ

(

) E

⋅

 

 

 

 

 

1238_Frame_C12.fm Page 506 Wednesday, June 12, 2002 5:24 PM

The bending strain–displacement matrix [B

(b)I

] for the Ith node is obtained by substitut-

ing the displacement approximation Equation 12.26 into Equation 12.12

(12.28)

where the elements of the strain–displacement matrix [B

(b)I

] are given by

(12.29)

12.1.4

Principle of Virtual Work

The effective membrane and bending forces are defined to describe the weak formulation
of shells

(12.30)

(12.31)

By making use of the basic kinematic assumption (Equation 12.4), the weak form of the

governing equation for thin shells under static load can be written as

(12.32)

Here W

ext

is the virtual work of the external loading given by

(12.33)

where is the applied resultant force per unit length and

is the applied direct couple per

unit length. and

are the prescribed resultant force and the prescribed director couple

on the boundaries

and

, respectively. For static analysis of thin shells, the penalty

method is used to enforce essential boundary conditions by adding an additional bound-
ary condition term to Equation 12.33 to obtain

(12.34)

















( )I

[

]

 

 

 

 

 

( )I

[

]

( )

–

I,11

I,1

ϕϕϕϕ

,11

ϕϕϕϕ

(

)

I,2

ϕϕϕϕ

,11

(

)

[

] E

⋅

( )I

[

]

( )

–

I,22

I,1

ϕϕϕϕ

,22

ϕϕϕϕ

(

)

I,2

ϕϕϕϕ

,22

(

)

[

] E

⋅

( )I

[

]

2 j

( )

–

I,12

I,1

ϕϕϕϕ

,12

ϕϕϕϕ

(

)

I,2

ϕϕϕϕ

,12

(

)

[

] E

⋅

n˜

βα

Sta

δ x

(

)

n˜

βα

δε

βα

δρ

βα

⋅

[

]d W

ext

δ x

(

)

–

∫

ext

n δϕ

ϕϕϕ m˜

δ t

⋅

[

δ ϕϕϕϕ j dΓ

δ t j dΓ

∫

⋅

∫

Sta

δx

( )

ααα

α u u

–

(

)

δudΓ

⋅

∫

–

1238_Frame_C12.fm Page 507 Wednesday, June 12, 2002 5:24 PM

Here

are the nodal displacement vector and the prescribed displacement vector on

all the essential boundary

, and

ααα

α is a diagonal matrix of penalty factors, which are

usually very large numbers (see

Chapter 4

For free-vibration analyses of the thin shells, the discrete system equation and the

boundary conditions are formulated separately. The variational form of the elastic dynamic
undamped equilibrium equation can be written as follows:

(12.35)

The weak form of the essential boundary conditions with Lagrange multipliers is employed

to obtain the discretized essential boundary conditions

(12.36)

where

λλλλ is a matrix of Lagrange multipliers, each of them can be interpolated on the

essential boundary using its nodal values.

(12.37)

and

(12.38)

where s is the arc-length along the essential boundaries, and N

(s) are the Lagrange

interpolates discussed in Equation 6.10.

12.1.5

Surface Approximation

The (neutral) surface of the shell is also approximated using MLS shape functions. The
procedure is exactly the same as for approximating the displacement field variables. The
approximated surface can be described by

(12.39)

where

ζζζζ is the coordinate in the parameter space for the neutral surface of the shell. A

deficiency of this approximation is that the constructed surface does not pass through the
prescribed points, unlike that in finite element meshes. This is due to the use of MLS shape
functions.

12.1.6

Discretized Equations

For static analyses of shells, the penalty term to impose essential boundary conditions in
Equation 12.34 can be discretized as follows:

(12.40)

, u

Dya

δx

( )

n˜

βα

δε

βα

δρ

βα

⋅

[

δu ρu˙˙d W

ext

δx

( )

–

⋅

∫

δ λλλλ

–

(

)dΓ

⋅

∫

λ x

( )

∈

∑

δλ x

( )

δλ

∈

∑

ϕϕϕϕ ζζζζ

( )

ζζζζ

( )x

ααα

δudΓ

⋅

∫

δu˜

ααα

αΦ

u˜

∑

1238_Frame_C12.fm Page 508 Wednesday, June 12, 2002 5:24 PM

and

(12.41)

where n is the number of sampling points for integration on surface

; is the nodal

vector corresponding to a degree of freedom in the translation and rotation fields. The
penalty matrix in Equation 12.40 and the penalty vector in Equation 12.41 are assembled
into the global stiffness matrix and the global external force vector, respectively.

For dynamic analysis of shells, Equations 12.1, 12.23, 12.24, 12.27, and 12.28 are substi-

tuted into the variational weak form (Equation 12.35). This gives the dynamic discrete
equation of

(12.42)

where K, M, and F are, respectively, the global stiffness, global mass matrices, and global
force vector, which are assembled using the corresponding nodal matrices and vectors
formed in the similar manner as those for plates (see

Chapter 10).

12.1.7

Static Analysis

For static analysis of shells, all the terms in Equation 12.42 related to dynamic effects
should vanish, and the equation system is simplified as

(12.43)

which is a set of linear algebraic equations that can be solved for the deflection using
standard routines of equation solvers.

12.1.8

Free Vibration

For free-vibration analysis, the external force vector F should vanish, and the damping
should not be considered; we then have

(12.44)

Considering harmonic vibration, the eigenvalue equations for shells derived from Equa-
tion 12.44 is of the form

(12.45)

where

ω is the natural frequency and Q is a vector that collects the nodal values corre-

sponding the amplitudes of the displacements given by

(12.46)

ααα

δudΓ

⋅

∫

δu˜

ααα

αΦ

∑

u˜

˙˙

+ U˙ KU

˙˙

–

(

…,Q

{

}

1238_Frame_C12.fm Page 509 Wednesday, June 12, 2002 5:24 PM

We now derive the discretized form of essential boundary conditions. Substituting the

displacement field u of Equation 12.1 into the weak form (Equation 12.36) yields a set of
linear algebraic constraint equations:

(12.47)

where

(12.48)

In general, B is a very sparse and singular matrix. For eigenvalue analysis, the essential
boundary conditions are homogeneous, i.e.,

In treating the essential boundary conditions, we follow the procedure described in

Section 11.1. An orthogonal transformation is performed to produce a positive-definite
stiffness matrix for the eigenvalue equation in computing natural frequencies. Using
singular-value decomposition, B can be decomposed as

(12.49)

where U and V are orthogonal matrices. S has a diagonal form whose diagonal elements
are equal to singular values of B. The matrix V can be written as

(12.50)

where the rank r of B is equal to the number of independent constraints, and the others
are redundant. The following change of coordinates satisfied the constraint Equation 12.47:

(12.51)

By substituting Equation 12.51 into Equation 12.45 and left-multiplying the result by

the transpose of

, the reduced order eigenvalue problem for the structure is

obtained

(12.52)

where

and

are the dimension reduced stiff-

ness and mass matrices. Equation 12.52 is now non-negative definite and can be solved
using standard eigenvalue solvers for eigenvalues that relate to the natural frequencies
and eigenvectors that relate to the vibration modes of the thin shells.

12.1.9

Forced (Transient) Vibration

For transient response analysis, Equation 12.42 can be solved using conventional direct
integration techniques based on finite difference approaches, where the time space is
divided into small time steps. The Newmark method is applied in the section for time
integration to obtain the time history of the displacement response of the shell. The
procedure is the same as that described in

Chapter 7

for 2D solids.

( ),

u x

( ), I Γ

, J

∈

USV

×r

, V

× n−r

(

)

(

)

× n−r

(

)

× n−r

(

)

–

(

)Q˜

−r

(

)×n

× n−r

(

)

−r

(

)×n

× n−r

(

)

1238_Frame_C12.fm Page 510 Wednesday, June 12, 2002 5:24 PM

12.1.10

Numerical Example for Static Problems

Example 12.1

Static Deflection of a Barrel Vault Roof under Gravity Force

The present EFG method is used to investigate the response of a barrel vault roof under
self-weight that is calculated simply using the density of gravity. The problem has been
analyzed by several researchers using the EFG method (Krysl and Belytschko, 1996b;
Noguchi et al., 2000) and FEM (Simo et al., 1989). The barrel vault roof is a standard
benchmarking test because it undergoes complex membrane and inextensional bending
states of stress. The example is used to evaluate the convergence and accuracy of the
present EFG method for the static analysis of shells. Figure 12.2 shows the barrel roof and
defines the parameters used in the description of its geometry. The two curved edges of
the roof are diaphragm supported, which allows displacement in the axial direction and
rotation about the tangent to shell boundary. The following parameters are used in the
computation.

Due to symmetry, only a quarter of the vault roof is modeled, and symmetric boundary

conditions are introduced along the planes of symmetry. To evaluate the effectiveness of
the present method, both regular and irregular nodal arrangements in the parametric
space shown in

Figure 12.3

are used in the analyses. Both quadratic (m

= 6) and quartic

FIGURE 12.2
Barrel vault roof and the coordinate system.

Length

Radius

Thickness

Semi-Span

Angle

Young’s

Modulus

Poisson’s

Ratio

Mass

Density

= 600

= 300

= 3.0

θ = 40°

= 3.0 × 10

ν = 0

ρ = 0.20833

1238_Frame_C12.fm Page 511 Wednesday, June 12, 2002 5:24 PM

= 15) basis functions are used in the analysis.

Figures 12.4

and

12.5

show, respectively,

the distributions of the vertical deflection along edges AB and AC, where the FEM solution
is obtained using a general-purpose program, i.e., ANSYS



. The EFG results using 81

nodes agree well with the finite element results of ANSYS using 2520 nodes. The reason
for the demand by FEM for a very fine element mesh may be the simplification of the
geometry, because the field variable (deflection) do not vary drastically and, as shown in
Figures 12.4 and 12.5, a coarse mesh should produce a good approximation. However, in
ANSYS, the shell is modeled using flat shell elements. This simplification cannot model
the mechanics coupling effects of bending forces and membrane forces, which is very
significant for shell structures, unless a very fine mesh is used. In our EFG code, however,
the geometry of the shell is modeled using MLS approximation, which very accurately
represents the curvature of the geometry of the shell, and hence produces very accurate
results using very coarse nodes. This example clearly demonstrates the advantage of the
MFree method in modeling shell structures. Note also that under the action of gravity
force, which is downward, the central point on the roof moves upward against the direc-
tion of body force. This fact provides very clear evidence of the importance of the coupling
of bending force and membrane force. The existence of the membrane force causes the
rise of the rooftop.

The convergence of the EFG results is also investigated in detail, and the results are

summarized in

Figure 12.6,

where the curves are normalized using the converged value

of the vertical deflection at B, that is,

−3.618. It can be seen that the performance of the EFG

method compares well with the finite element solution of Belytschko et al. (1985) and
Simo et al. (1989). The results given by quartic polynomial basis are usually better than
that by quadratic basis when nodal points are sufficiently dense. However, worse results
are obtained when low nodal densities are employed. The reason could be again the
approximation of the geometry of the shell. When the nodes are too coarse, even EFG will
fail to approximate the geometry well.

Lagrange Multiplier Method vs. Penalty Method
The penalty method and Lagrange multipliers can both be used to enforce the essential
boundary conditions. Here we discuss the advantages and disadvantages of these two
methods using the example of barrel vault roof.

Figure 12.5

shows the distribution of the

vertical deflection along edge DF, where deflection is supposed to be zero because a zero
essential boundary condition was imposed in the EFG analysis. It can be seen from Figure 12.5

FIGURE 12.3
Nodal arrangement in parameter space.

(a) (b)

1238_Frame_C12.fm Page 512 Wednesday, June 12, 2002 5:24 PM

that both methods could not give exact results for the zero displacement. Very small errors
are observed especially at the two corners. The Lagrange multipliers method is more
accurate than that of the penalty method. In the penalty method, a penalty factor of 5000
times the value of the elastic modulus is used. This confirms the discussion given in

Chapter 4

on the comparison of the Lagrange multipliers method and penalty method.

Note again that the penalty method is much more computationally efficient.

FIGURE 12.4
(a) Vertical deflection of the barrel vault roof along section AB (see Figure 12.2); (b) vertical deflection along
section AC.

-4.0

-3.0

-2.0

-1.0

0.0

1.0

FEM-ANSYS (1353 nodes)

FEM-ANSYS (2520 nodes)

EFG(81nodes,quartic,regular)

EFG(81nodes,quartic,irregular)

Angle (degree)

Deflection

(a)

0.0

0.1

0.2

0.3

0.4

0.5

0.6

100

150

200

250

300

FEM-ANSYS (1353 nodes)

FEM-ANSYS (2520 nodes)

EFG(81nodes,quartic,regular)

EFG(81nodes,quartic,irregular)

Length

Deflection

(b)

1238_Frame_C12.fm Page 513 Wednesday, June 12, 2002 5:24 PM

12.1.11

Numerical Examples for Free Vibration of Thin Shells

The performance of the EFG method for free-vibration analysis is also evaluated. The
results of several examples are presented in comparison with analytical solutions and
other results.

Example 12.2

Free Vibration of a Clamped Cylindrical Shell Panel

A panel of cylindrical thin shallow shell, which has been investigated by Petyt (1971), is
also examined here to benchmark our EFG code for free-vibration analysis. The thin shell
panel is clamped at all edges, and the natural frequencies are computed using the present
EFG. The geometry and boundary conditions of the shell panel are shown in

Figure 12.7

and the following parameters are used in the analysis.

FIGURE 12.5
EFG solution of vertical deflection of the barrel vault roof along section DF. The vertical deflections are supposed
to be zero. The numerical solution gives very small values. The largest errors were observed at corners.
The Lagrange multipliers method is more accurate than the penalty method. (a) Lagrange multipliers method;
(b) penalty method.

1238_Frame_C12.fm Page 514 Wednesday, June 12, 2002 5:24 PM

Length
L (mm)

Radius

R (mm)

Thickness

h (mm)

Span

Angle

θθθθ

Young’s

Modulus

E (N/m

)

Poisson’s

Ratio

νννν

Mass

Density

ρρρρ (kg/m

)

76.2

762

0.33

7.64

6.8948

× 10

0.33

2657.3

FIGURE 12.6
Convergence of results of vertical deflection at B.

FIGURE 12.7
Clamped panel of cylindrical thin shell.

0.95

1.00

1.05

1.10

1.15

1.20

100

150

200

250

300

EFG(quartic)

EFG(quadratic)

FEM(Belytschko et al. 1985)

FEM(Simo et al. 1989)

Total number of nodes

Deflection

1238_Frame_C12.fm Page 515 Wednesday, June 12, 2002 5:24 PM

Regular nodes of different densities are used to investigate the convergence character-

istics of natural frequencies. Results are given in terms of a frequency parameter defined
as

λ =

= (

ρhR

/D)

1/8

, where D

= (Eh

/(12(1 −

)), and are shown in Figure 12.8. As

can be seen from Figure 12.8, very high convergence rates for the various vibration modes
are obtained. The values of

λ are also tabulated in Table 12.1 together with experimental

results and results obtained by other, different numerical methods (Petyt, 1971). In this
table, ERR denotes the extended Raleigh–Ritz method, FET stands for the triangular finite
element method, and FER denotes rectangular finite element method. It is found that the
EFG results show good convergence and good agreement with other methods.

FIGURE 12.8
Convergence of results for natural frequency parameter

λ of a clamped panel of a cylindrical shallow thin shell.

TABLE 12.1

Natural Frequency Parameter

λ of a Clamped Panel of Cylindrical Thin Shell

Mode

Experimental

Results

ERR

FET

FER

Present Method (EFG)

×××× 6

×××× 9

×××× 12

×××× 16

13.06

13.28

13.35

13.32

13.28

13.54

13.60

13.65

13.67

13.62

13.60

14.57

14.65

14.71

14.79

14.68

14.65

14.70

14.86

14.85

14.88

14.87

14.86

14.85

15.09

15.06

17.21

15.11

15.08

15.06

15.93

15.82

15.83

15.87

15.95

15.85

15.83

15.82

15.78/15.86

15.91

15.88

15.96

16.67

16.03

15.89

15.88

16.55

16.46

16.49

16.95

16.61

16.47

16.46

16.79

16.78

16.81

17.25

16.91

16.80

16.78

16.89

16.93

16.92

16.96

17.77

17.14

16.95

16.92

ERR: extended Raleigh–Ritz method; FET: triangular finite element method;
FER: rectangular finite element method.

13.0

13.5

14.0

14.5

15.0

15.5

16.0

16.5

17.0

17.5

18.0

100

150

200

250

300

No. 10

No. 9

No. 8

No. 7

No. 6

No. 5

No. 4

No. 3

No. 2

No. 1

Number of nodes

Nondimensional frequency

Ω

1238_Frame_C12.fm Page 516 Wednesday, June 12, 2002 5:24 PM

Example 12.3

Free Vibration of a Hyperbolical Shell

Free-vibration analysis of the clamped-free hyperbolical shell shown in Figure 12.9 is
performed. The shell geometry is defined by the following equation of its meridian:

(12.53)

where b is a characteristic dimension of the shell, which can be calculated as b

= ad/

The following geometric and material properties are used in the analysis:

A regular nodal arrangement in the axial and circumferential direction is used in the
analysis. The results of the present study using different numbers of nodes are given in

Table 12.2.

Results obtained by Carter et al. (1969) using a numerical integration technique

and Ozakca and Hinton (1994) using an FEM are also listed in the same table for an easy
comparison. Close agreements are again observed.

Example 12.4

Free Vibration of a Cylindrical Shell

The natural frequencies of a cylindrical shell are examined using the present EFG code.
The geometry of the shell is schematically drawn in

Figure 12.10.

The shell is clamped at

one edge and free at the other. The following geometry and material properties are used:

FIGURE 12.9
Geometry of a hyperbolical shell that is clamped at the bottom circular edge and free at the top circular edge.

Length
L (m)

Height

h (m)

Height

d (m)

Radius

a (m)

Radius

b (m)

Young’s

Modulus

E (N/m

)

Poisson’s

Ratio

νννν

Mass

Density

ρρρρ (kg/m

)

100.787

0.127

82.194

25.603

63.906

2.069

× 10

0.15

2405

Length
L (mm)

Thickness

h (mm)

Radius

R (mm)

Young’s Modulus

E (N/m

)

Poisson’s

Ratio

νννν

Mass Density

ρρρρ (kg/m

)

226.8

1.021

106.1

2.069

× 10

0.3

7868

----

 

 

L d

–
b

------------









–

1238_Frame_C12.fm Page 517 Wednesday, June 12, 2002 5:24 PM

The nodes used in our EFG are regularly arranged in the axial and circumferential direc-
tions. The EFG results are tabulated in Table 12.3. The results obtained by Tou and Wong
(1987) and Subir and Gould (1974) using high-precision finite element method (HPFEM)
and standard FEM, respectively, are also listed in the same table for comparison. Again,
good agreement is evident.

TABLE 12.2

Natural Frequencies (Hz) of a Clamped-Free Hyperbolical Shell

Mode

Ozakca and

Hinton (1994)

Carter et al.

(1969)

Present Method (EFG)

×××× 16

×××× 24

1.0354

1.0348

1.0351

1.0325

1.1508

1.1467

1.1486

1.1450

1.1826

1.1808

1.1809

1.1780

1.3061

1.3015

1.3043

1.2998

1.3293

1.3231

1.3254

1.3223

1.3758

1.3749

1.3799

1.3753

1.4329

1.4293

1.4284

1.4259

1.4488

1.4475

1.4497

1.4470

FIGURE 12.10
Geometry of a cylindrical shell.

TABLE 12.3

Natural Frequencies (Hz) of the Clamped-Free Cylinder

Mode

FEM (Subir and

Gould, 1974)

HPFEM (Tou and

Wong, 1987)

Present Method (FEG)

×××× 16

×××× 24

487

482

490

483

565

561

564

562

621

616

629

624

875

869

982

981

979

980

= no results were given.

1238_Frame_C12.fm Page 518 Wednesday, June 12, 2002 5:24 PM

12.1.12

Numerical Examples for Forced Vibration of Thin Shells

To investigate the accuracy as well as the capability of the EFG method for forced vibration
problems for thin shells, numerical examples for thin shells subjected to different transient
excitations are presented, and the results obtained are compared with those of ordinary
FEM and other numerical methods.

Example 12.5

Clamped Circular Plate Subject to an Impulsive Load

A code developed for shells should also be able to work for plates. The first test case is a
circular plate subjected to rectangular impulsive force. The parameters used in this calcu-
lation are given as follows:

The natural frequencies are first computed by the EFG method, and the results are

compared with those obtained using the boundary element method (BEM) and FEM. We
state without showing the results that they are in very good agreement.

The circular plate is loaded by a concentrated vertical impulsive force F(t) at its center.

The magnitude of the force is 1 N, and the duration of pulse is t

= 0.121 s as shown in

Figure 12.11.

In time stepping, the time increment used is

∆t = 0.001 s.

Figure 12.12

shows

the history of the vertical deflection at the center of the plate obtained by the presented
method, the analytical solution of Sneddon (1945), SAP IV (Bathe et al., 1973), and the
domain-boundary element method by Beskos (1990). It is observed that the results
obtained using the present EFG code agree well with all these methods including the
analytical results by Sneddon (1945).

Example 12.6

Clamped Cylindrical Shell Subject to a Sine Load

The transient response of a clamped cylindrical shell as shown in

Figure 12.10

is investi-

gated. The following geometry and material properties are used:

Figure 12.13

shows the history of the transient excitation of half a sine function of time.

The excitation is at the center of the meridian (y

= L/2, z = R). The force is expressed as

(12.54)

FIGURE 12.11
Transient force of a rectangular pulse applied to the circular plate.

Thickness

h (mm)

Radius

R (mm)

Young’s Modulus

E (N/m

)

Poisson’s

Ratio

νννν

Mass Density

ρρρρ (kg/m

)

0.05

1.1

1000.0

0.3

0.229

Length
L (mm)

Thickness

h (mm)

Radius

R (mm)

Young’s Modulus

E (N/m

)

Poisson’s

Ratio

νννν

Mass Density

ρρρρ (kg/m

)

600

3.0

300.0

2.1

× 10

0.3

7868

0.0

1.0

2.0

F (t)

1000t

(

)

sin

1238_Frame_C12.fm Page 519 Wednesday, June 12, 2002 5:24 PM

where F

= 1000.0 N. Good agreement is observed. In performing the time marching, the

time step

∆t = 2.5e − 5s is used, which is almost the 1/35 of the fundamental period of the

cylinder. Regular nodes (12

× 16) are arranged in the axial and circumferential directions.

Figure 12.14

shows the transient response of a clamped cylinder subjected to a transient

excitation of half a sine function of time. The results are obtained using both the present
EFG code and FEM.

FIGURE 12.12
The time history of the deflection at the center of a clamped circular plate subjected to a rectangular impulsive
excitation.

FIGURE 12.13
Time history of the external force excitation of half a sine function of time.

-5.0

0.0

5.0

10.0

15.0

20.0

0.0

2.0

4.0

6.0

8.0

10.0

12.0

Beskos, 1990

SAP IV

Sneddon, 1945

EFG

)

(1000 × wD)

25x

(

)

200

400

600

800

1000

1200

0.0005

0.001

(N)

t (s)

1238_Frame_C12.fm Page 520 Wednesday, June 12, 2002 5:24 PM

Figure 12.15 shows the transient response for the cylindrical shell, but the external

excitation force becomes F

= F

sin(2000t). Again, very good agreement between the results

of EFG and FEM is obtained.

Example 12.7

Clamped Spherical Shell Subject to a Sine Curve Load

A spherical cap of thin shell with a central angle of

θ = 120° shown in

Figure 12.16

are

investigated. The shell is subjected to an excitation of a force of half a circle sine function of
time at the apex. The time history of the loading is given by Equation 12.54, but F

= 200.0 N.

FIGURE 12.14
Transient response of the vertical displacement at the central point in the meridian of the cylindrical shell (y

/2, z = R).

FIGURE 12.15
Same as Figure 12.14, but F

= F

sin(2000t).

-2.5E-05

-2.0E-05

-1.5E-05

-1.0E-05

-5.0E-06

0.0E+00

5.0E-06

0.0005

0.001

0.0015

0.002

0.0025

0.003

0.0035

0.004

EFG

FEM

w(m)

t(s)

-3.0E-05

-2.5E-05

-2.0E-05

-1.5E-05

-1.0E-05

-5.0E-06

0.0E+00

5.0E-06

1.0E-05

1.5E-05

0.0005

0.001

0.0015

0.002

EFG

FEM

w(m)

t(s)

1238_Frame_C12.fm Page 521 Wednesday, June 12, 2002 5:24 PM

This spherical cap is clamped on the circular boundary at the bottom. The following
geometry and material properties are used in the computation:

The time step used is

∆t = 5.0e − 5s, which is almost the 1/25 of the fundamental period

of the spherical cap; 185 nodes are used in the calculation. Figure 12.17 shows the history
of the vertical displacement response at the apex of the spherical cap obtained by both
the present EFG method and FEM. Very good agreement is obtained.

FIGURE 12.16
Cross section of the spherical cap of a thin shell.

FIGURE 12.17
Transient response of the vertical displacement at the apex of the spherical cap.

Thickness

h (mm)

Radius

R (mm)

Young’s Modulus

E (N/m

)

Poisson’s

Ratio

νννν

Mass Density

ρρρρ (kg/m

)

900.0

2.0e

0.3

7800

-8.0E-06

-6.0E-06

-4.0E-06

-2.0E-06

0.0E+00

2.0E-06

4.0E-06

0.0005

0.001

0.0015

0.002

0.0025

0.003

0.0035

0.004

EFG

FEM

w(m)

t(s)

1238_Frame_C12.fm Page 522 Wednesday, June 12, 2002 5:24 PM

12.1.13

Remarks

The EFG method has been formulated for static, free-vibration, and forced vibration analysis
of spatial thin shell structures. In the EFG method, the MLS technique is used for approx-
imations of both the surface geometry of the shell and the field variables.

The present EFG results are benchmarked with a number of examples by comparison

with the results obtained by other methods. We have more than enough evidence to show
that EFG works for shells.

The present EFG method offers distinct computational advantages over classical FEM:

(1) mesh free, (2) accurate geometry representation, (3) fast convergence.

12.2 EFG Method for Thick Shells

The Kirchhoff–Love shell theory works only for very thin shells. The reason is very much
similar to the reason we argued for the theories for plates (see

Chapter 11

For thick shells,

formulations need to be based on thick shell theories. The formulation for thick shells in
this section is based on the geometrically exact theory of flexible shells proposed by Simo
et al. (1989). We collect the necessary equations in the following, but refer the reader to
the articles by Simo et al. (1989) for details.

12.2.1

Fundamental Relations

The Gauss intrinsic coordinates are used to describe the configuration of the shell as
described in previous sections. Making use of the definition of spatial tensors, the corre-
sponding linearized strain measures are defined relative to the dual spatial surface basis as

(12.55)

(12.56)

(12.57)

Here,

∆t is the incremental spatial rotation defined by

(12.58)

where

(12.59)

and is the (3

× 2) matrix given by

(12.60)

αβ

1
2

---

ϕϕϕϕ

⋅

(

)

ϕϕϕϕ

∆t u

⋅

(

)

αβ

1
2

---

ϕϕϕϕ

∆t

ϕϕϕϕ

∆t

⋅

(

)

∆t

Λ ∆T

⋅

∆T

[

]

1238_Frame_C12.fm Page 523 Wednesday, June 12, 2002 5:24 PM

which can be obtained by deleting the third column of , where

is the

orthogonal transformation, and can be expressed as

(12.61)

The constitutive relations for the effective membrane and for the stress couple result-

ant have the form of Equations 12.23 and 12.24. The constitutive relations for the
effective shear stress resultants can be written as

(12.62)

Here is the shear reduction coefficient and G is the shear modulus.

In calculating the effective stress couple resultant

, Equations 12.24 and 12.31 should

be used. In calculating the effective membrane , Equation 12.30 and the following
equation should be used.

(12.63)

In calculating the effective shear resultant forces, we use

(12.64)

where

is defined as

(12.65)

In calculating t

, we use

(12.66)

The above equation shows that the difference between thick shell theory and thin shell

theory is that the shear effects have been taken into account in the thick shell theory.

12.2.2

Principle of Virtual Work

The weak form of the governing equation for the thick shells under static load can be
written as

(12.67)

Here

is the current surface measure.

is the virtual work of the external

loading given by

(12.68)

where is the applied resultant force per unit length and

is the applied direct couple

per unit length. and

are the prescribed resultant force and the prescribed director

couple on their corresponding boundaries

and

, respectively.

Λ:

⊂

→

t E

⋅

(

[

]

t E

⋅

-------------------

(

)

(

)

⊗

n˜

q˜

 

 

 

κ Gh

 

 

 

n˜

βα

αµ

–

q˜

αµ

–

Sta

δx

( )

n˜

βα

δε

βα

⋅

βα

δκ

βα

q˜

δγ

⋅

[

]d W

ext

δx

( )

–

∫

j d

ext

δϕϕϕϕ

⋅

δ t

⋅

[

δϕϕϕϕ j dΓ

⋅

∫

δt j dΓ

∫

1238_Frame_C12.fm Page 524 Wednesday, June 12, 2002 5:24 PM

The MLS shape functions are then used to approximate the displacement fields. Substi-

tuting Equation 12.1 into the weak forms leads to a set of global system equations. The
procedure is very much the same as that in the previous section.

The treatment of essential boundaries is also very much the same as that in the previous

section. We need to use different ways to handle the essential boundary condition for
different problems, whenever MLS approximation is employed.

12.2.3

Numerical Examples

Example 12.8

Static Deflection of a Barrel Vault Roof under Gravity Force

Example 12.1 is reexamined using the EFG formulation for thick shells. The barrel vault
roof is shown in

Figure 12.2.

All the parameters and conditions are exactly the same as

those in Example 12.1. In constructing the MLS shape function, both quadratic and quartic
basis functions are used.

Figure 12.18

shows the results of convergence of vertical dis-

placement at point B in the barrel vault roof. Figure 12.18 shows that the convergence by
the present analysis is excellent in comparison with the results of using high-performance
FEM by Simo et al. (1989) and EFG for thin shells by Krysl and Belytschko (1996b). It can
also be seen that the convergent rate for quartic basis always exceeds that for quadratic
basis. This was not very clear from the results using thin shell theory (see Figure 12.6).
The result for quartic basis functions approaches the exact result from below with the
increase of EFG nodes, while that for quadratic basis functions oscillates along the exact
result with the increase of EFG nodes.

Figure 12.19

shows the variation of fractions of shear, membrane, and bend energies in

total energy, with respect to the size of the support domain. The computation is performed
using quartic basis function. It can be seen that the membrane and bending energies are
stable with the increase of the support domain; the difference between the membrane and
bending energies are no more than 4.0%. The shear energy is very small and approaches
zero with the increase of EFG nodes, which means that the shear stress plays a very minor
role in deformation of the vault roof.

FIGURE 12.18
Convergence of vertical displacement at point B in the barrel vault roof.

Number of elements/side

0.92

0.94

0.96

0.98

1.00

1.02

1.04

1.06

1.08

1.10

Normalized displacement

Simo et al., 1989

EFG-Krysl and Belytschko, 1996b

EFG(Quadratic)

EFG(Quartic)

1238_Frame_C12.fm Page 525 Wednesday, June 12, 2002 5:24 PM

Figures 12.20 and 12.21 show the distribution of the vertical deflection in the roof along

the edges AB and AC, respectively. The central point A moves upward, that is, in the
reverse direction of the body force, for the effects of the membrane forces. It can be seen
that thousands of nodes are needed to obtain the converged solution using a general-
purpose program, e.g., ANSYS, for the reason mentioned in Example 12.1. Many fewer
nodes are needed for the present EFG method.

Example 12.9

Pinched Cylindrical Shell

The deformation of the cylindrical shell shown in

Figure 12.22

is analyzed using the

present EFG method. The shell is loaded by a pair of centrally located and diametrically
opposed concentrated forces. The cylindrical boundaries are supported by a rigid diaphragm

FIGURE 12.19
Variation of membrane, shear, and bending energy with respect to size of the support domain.

FIGURE 12.20
Vertical displacement in the barrel vault roof along AB.

0.0%

10.0%

20.0%

30.0%

40.0%

50.0%

60.0%

Shear Energy

Membrane Energy

Bending Energy

Normalized size of support domain

Fraction of total energy

Angle (degree)

-4.0

-3.0

-2.0

-1.0

0.0

1.0

FEM-ANSYS (1353 nodes)

FEM-ANSYS (2520 nodes)

EFG-Present (m=15)

Vertical displacement

1238_Frame_C12.fm Page 526 Wednesday, June 12, 2002 5:24 PM

that allows displacement in the axial direction and rotation about the tangent to the shell
boundary. The parameters for this problem are given in the following table:

This problem is one of the most critical tests for numerical methods for both inextensional
bending and complex membrane states of stress. Because of its double symmetry, only 1

of the cylinder has been modeled. The results are shown in

Figure 12.23

which are normal-

ized by the value of 1.8248e-5, which is the convergent numerical solution of magnitude
of the radial displacement at loaded points. The convergence by the present method is
excellent in comparison with the results by FEM (Simo et al., 1989) and by EFG for thin

FIGURE 12.21
Vertical displacement in the barrel vault roof along AC.

FIGURE 12.22
Pinched cylindrical shell.

Length
L (mm)

Thickness

h (mm)

Radius

R (mm)

Young’s Modulus

E (N/m

)

Poisson’s

Ratio

νννν

600

3.0

300.0

3.0

× 10

0.3

Length

0.0

0.1

0.2

0.3

0.4

0.5

0.6

100

150

200

250

300

Vertical displacement

EFG-Present (m=15)

FEM-ANSYS (1353 nodes)

FEM-ANSYS (2520 nodes)

sym

F =1.0

L/2

1238_Frame_C12.fm Page 527 Wednesday, June 12, 2002 5:24 PM

shells (Krysl and Belytschko 1996b). It is seen that the EFG results reported by Krysl and
Belytschko (1996b) converge faster than the present EFG code.

The effects of the dimension of the support domain on the radial displacement at the

loading point are also investigated, and the results are shown in Figure 12.24. The shortest
nodal distance between nodes is used to normalize the size of the support domain. The
fluctuation of the displacement is visible but it is less than 3% when the support domain
is larger than 3.8. Effects of the dimension of the support domain on the energy compo-
nents (bending, membrane, and shear energies) are also investigated, and the results are
summarized in

Figure 12.25.

The energies show a trend of fluctuation similar to that of

the displacement.

Figure 12.26

shows the convergence of membrane (M), shear (S), and

FIGURE 12.23
Convergence of the results of the radial displacement at the center of the cylindrical shell.

FIGURE 12.24
Effects of the dimension of the support domain on the radial displacement at the loading point.

40.0%

50.0%

60.0%

70.0%

80.0%

90.0%

100.0%

110.0%

Number of elements/side

Normalized displacement

Simo et al., 1989

Krysl and Belytschko, 1996b

EFG-Present (m=15)

60.0%

70.0%

80.0%

90.0%

100.0%

110.0%

3.5

4.5

5.5

6.5

7.5

8.5

Normalized size of support domain

Normalized displacement

1238_Frame_C12.fm Page 528 Wednesday, June 12, 2002 5:24 PM

bending (B) energy fractions in the total energy. The results are calculated using both
quartic and quadratic polynomial basis functions. As can be seen, for quartic polynomial
basis, the shear energy converges to a value of less than 4.0% as the field nodes increase,
which means shear stress plays but a small role in the deformation. The membrane energy
fraction approaches approximately 40.0% which gives very satisfactory results in compar-
ison with the results using nine-noded gamma FEM by Belytschko et al. (1985). For the
quadratic polynomial basis, the membrane fraction approaches 40.0%, same as that for
the quartic polynomial. However, the shear energy increases with the increase of field
nodes. It can also be observed that for coarse field nodes, membrane locking is dominant
while for the fine field nodes, shear locking becomes more important.

FIGURE 12.25
Effects of the dimension of the support domain on the energy components (from top to bottom: bending,
membrane, and shear energy).

FIGURE 12.26
Convergence of membrane (M), shear (S), and bending (B) energies.

0.0%

10.0%

20.0%

30.0%

40.0%

50.0%

60.0%

70.0%

3.5

4.5

5.5

6.5

7.5

8.5

Normalized size of support domain

Fraction of total energy

Number of elements/side

Fraction of total energy

0.0%

10.0%

20.0%

30.0%

40.0%

50.0%

60.0%

70.0%

80.0%

M (m=6)

M (m=15)

S (m=6)

S (m=15)

B (m=6)

B (m=15)

1238_Frame_C12.fm Page 529 Wednesday, June 12, 2002 5:24 PM

Example 12.10

Pinched Hemispherical Shell

In this example, a pinched hemispherical shell shown in Figure 12.27 is analyzed using
the EFG method. The shell is pinched by two pairs of opposed radial point loads of
magnitude P

= 2.0. The bottom circumferential edge of the hemisphere is free. The material

parameters are E

= 6.825 × 10

N/m

and

ν = 0.3. The sphere radius is R = 10.0 and the

thickness t

= 0.04.

This problem is also a challenging benchmarking problem to check whether the formu-

lation of the shell can represent inextensional modes of deformation, as it exhibits almost
no membrane strains. This problem is a less critical test with regard to inextensional bending,
compared to the pinched cylinder problem. However, it is a very useful problem for checking
the ability of the present EFG method to handle rigid body rotations normal to the shell
surface. Large sections of this shell rotate almost as rigid bodies in response to this load,
and the ability to accurately model rigid body motion is essential for any numerical method.

Due to symmetry, only a quarter of the hemisphere is modeled. The results are shown

Figure 12.28.

The radial displacements in the direction of loads at loaded points are the

same and are found analytically to be 0.0924. This value is used to normalize the numerical
results presented in Figure 12.28. The convergence is even better than that of the high-
performance results by Simo et al. (1989) using mixed formulation. Among the EFG
solutions, the present method also gives better accuracy than that of the thin shell formu-
lation by Krysl and Belytschko (1996b). The results obtained using the quartic polynomial
basis functions are better than those using quadratic polynomial basis functions. However,
poor accuracy is obtained when very coarse nodes are used in EFG compared with FEM
using the same density of nodes.

Figure 12.29

shows the effects of the dimension of the

support domain on the radial displacement at the loading point on the hemispherical
shell. The present EFG results using quartic basis functions do not depend on the size of
the support domain and gives an almost exact solution, in the range of our investigation.
The results obtained using quadratic basis functions are stable and accurate only in a
narrow range (

to 6.5).

Figure 12.30

shows the effects of the dimension of the

support domain on membrane (M) and shear (S) energies in the hemisphere shell.

Figure 12.31

shows the same effects but for the bending (B) energy in the hemisphere shell.

FIGURE 12.27
Pinched hemisphere shell. One quarter of the shell is shown.

F =1.0

sym

free

4.5

1238_Frame_C12.fm Page 530 Wednesday, June 12, 2002 5:24 PM

It can be seen that both the shear and the membrane energy are very small when the
quartic basis functions are used. Although the shear energy obtained using the quadratic
basis functions is very small, the membrane energy increases with the increase of the
support domain. This means that the membrane energy is overpredicted with the increase
of field nodes, which is an indication of membrane locking.

Figures 12.32

and

12.33

show the fractions of shear, membrane, and bending energies in

the total strain energy calculated using both quartic and quadratic polynomial basis functions.

FIGURE 12.28
Convergence the radial displacement at the loading point on the hemisphere shell.

FIGURE 12.29
Effects of the dimension of the support domain on the radial displacement at the loading point on the hemisphere
shell.

Number of elements/side

60.0%

65.0%

70.0%

75.0%

80.0%

85.0%

90.0%

95.0%

100.0%

105.0%

Normalized displacement

Simo et al., 1989

Krysl and Belytschko, 1996b

EFG-Present (m=15)

EFG-Present (m=6)

90.0%

95.0%

100.0%

105.0%

Normalized size of support domain

EFG-Present (m=6)

EFG-Present (m=15)

Normalized displacement

1238_Frame_C12.fm Page 531 Wednesday, June 12, 2002 5:24 PM

It is found that when quartic polynomial basis functions are used, both membrane and
shear energy tend to zero as the field nodes increase. This implies that the shear stress
and membrane stress play trivial roles in the deformation of the shell, and the bending
energy is dominant. It is also found that membrane locking is dominant for coarse field
nodes, but is quickly eliminated with the increase of field nodes.

12.2.4

Remarks

The EFG method has been extended to analyze thick spatial shells based on the stress-
resultant shell theory, which has five degrees of freedom assigned to every point of the
shell.

FIGURE 12.30
Effects of the dimension of the support domain on membrane (M) and shear (S) energies in the hemisphere shell.

FIGURE 12.31
Effects of the dimension of the support domain on bending (B) energy in the hemisphere shell.

Normalized size of support domain

0.0%

2.0%

4.0%

6.0%

8.0%

10.0%

12.0%

14.0%

Fraction of total energy

M (m=15)

M (m=6)

S (m=15)

S (m=6)

Normalized size of support domain

Fraction of total energy

80.0%

85.0%

90.0%

95.0%

100.0%

B (m=15)

B (m=6)

1238_Frame_C12.fm Page 532 Wednesday, June 12, 2002 5:24 PM

For thick shells, the formulation allows transverse shear strain and results in the Mindlin

plate when the surface becomes flat and the membrane state is negligible. To avoid shear
locking and membrane locking, quartic basis function is recommended for MLS approx-
imation. Using high-order basis functions is an effective way to avoid shear locking. A
good alternative is to use different orders of shape function for deflection and rotation,
as discussed in

Chapter 11

for plates (see Example 11.16).

FIGURE 12.32
Convergence of bending energy.

FIGURE 12.33
Convergence of membrane and shear energies.

Number of elements/side

Fraction of total energy

60.0%

65.0%

70.0%

75.0%

80.0%

85.0%

90.0%

95.0%

100.0%

B (m=6)

B (m=15)

Fraction of total energy

Number of elements/side

0.0%

5.0%

10.0%

15.0%

20.0%

25.0%

30.0%

M (m=6)

M (m=15)

S (m=6)

S (m=15)

1238_Frame_C12.fm Page 533 Wednesday, June 12, 2002 5:24 PM

12.3 RPIM for Thick Shells

12.3.1 Formulation

Procedure

In the previous two sections, we presented EFG formulations for both thin and thick shells,
using MLS approximations. We have shown that special treatment is needed in dealing
with the essential boundary conditions, because of the lack of the Kronecker delta function
property in the MLS shape functions. The natural progress of our research is then to use
PIM shape functions, as we have already seen the advantages of PIM shape functions in
Chapters 8 to 11 for solids, fluids, beams, and plates. This section therefore presents RPIM
methods that use Galerkin formulation and RPIM shape functions for thick shells. RPIM
is then used for static analysis of several benchmark problems to illustrate the validity
and efficiency.

Formulation of RPIM was presented in detail in

Chapter 5,

and the formulation of shells

has been detailed in the previous two sections of this chapter. In formulating RPIM for
shells, all we need do is replace the MLS shape functions with RPIM shape functions, and
ignore all the terms related to essential boundary conditions. Therefore, we need not repeat
the process, and go directly to examination of the performance of the RPIM code we have
developed. Note that in the following examples nonconforming RPIM or NRPIM is used
in the computation. For results of CPIM, readers are referred to the recent work by G.R.
Liu et al. (2002c) where excellent convergence of CPIM is reported.

12.3.2 Numerical

Examples

Example 12.11

Barrel Vault Roof

Example 12.1 is reexamined using the RPIM formulation for thick shells. The barrel vault
roof is shown in

Figure 12.2.

All the parameters and conditions are exactly the same as

those in Example 12.1. The background of this benchmark problem is also mentioned in
Example 12.1. This example is used here for investigating the effects of shape parameters
in the radial basis functions. The shell roof is loaded by its own weight of uniform vertical
gravity load and is supported by diaphragms at the ends, but is free along the sides.
Because of symmetry, only one quarter of panel is modeled by the present RPIM method.
In presenting the results, we use the reference solution of

−3.618 for the vertical deflection

at point A to normalize the results calculated using different methods.

Both multiquartic (MQ) radial functions and Gaussian (EXP) radial basis functions are

used to construct the RPIM shape function (see

Chapter 5).

In using radial basis functions,

one needs to examine and fine-tune the shape parameters. We therefore use regularly
distributed nodes to investigate the shape parameters of these radial basis functions, as
defined in

Table 5.3.

Effect of Shape Parameters of MQ Radial Basis Functions
When q

= 0.5 and q = −0.5 (see Table 5.3), the MQ radial function becomes the original

MQ and the reciprocal MQ (RMQ), which have been traditionally used for many appli-
cations. The characteristic length d

is taken to be the “average” nodal spacing for nodes

in the support domain of the quadrature point (see Section 2.10.3).

We now start to investigate the effects of

. In our first study, we let

, where

is the normalized dimension of the support domain.

Figure 12.34

shows the effects of the

size of the support domain on the results of vertical displacement at the central top of the roof.

1238_Frame_C12.fm Page 534 Wednesday, June 12, 2002 5:24 PM

Figure 12.35

shows the effects of the size of the support domain on the fractions of shear,

membrane, and bending energies in the roof. The membrane-locking phenomenon exists
when too small a domain is utilized for approximation. It can be seen that the membrane
and bending energies become stable with the increase of the support domain; the difference
between the membrane and bending energies is no more than 4.0%. The shear energy is
very small and approaches zero with the increase of the support domain, which means
that the shear stress plays an almost negligible role in the deformation of the roof. The
fractions of membrane and bending energies become oscillating when the normalized

FIGURE 12.34
Effects of the size of the support domain on the results of vertical displacement at the central top of the roof.
MQ radial basis function is used with the parameter C the dimension of the support domain (

80%

90%

100%

110%

120%

130%

60%

80%

100%

120%

140%

160%

Normalized size of support domain

q =

- 0.5

Normalized displacement

= 0.5

Normalized size of support domain

Normalized displacement

1238_Frame_C12.fm Page 535 Wednesday, June 12, 2002 5:24 PM

support domain is greater than 12.0. Correspondingly, the error displacement of the center
point also becomes large and the results are erroneous. Based on the results shown in

Figures 12.34

and 12.35, we state that when the shape parameter

is chosen the same

as the dimension of the support domain, i.e.,

, the recommended dimension of the

support domain for the thin shell is approximately 8.3 to 11.4 and 6.2 to 11.2 times the
smallest nodal spacing for q

= −0.5 and q = 0.5, respectively.

FIGURE 12.35
Effects of the size of the support domain on the fractions of energies in the roof. MQ radial basis function is
used with the parameter C the dimension of the support domain (

10%

20%

30%

40%

50%

60%

70%

80%

90%

100%

Membrane energy

Shear Energy

Bending Energy

10%

20%

30%

40%

50%

60%

70%

80%

90%

Membrane energy

Shear Energy

Bending Energy

= 0.5

Fraction of total energy

Normalized size of support domain

q =

- 0.5

Fraction of total energy

Normalized size of support domain

1238_Frame_C12.fm Page 536 Wednesday, June 12, 2002 5:24 PM

Next, we fix the support domain at

= 10.0, and vary the dimensionless shape parameter

. Figure 12.36 shows the effects of the dimensionless shape parameter

on the dis-

placement results. This figure shows that to obtain good results, the range of parameter

should be 0.65 to 0.85 and 0.45 to 0.78 for q

= −0.5 and q = 0.5, respectively. For

convenience,

= 0.7 is used in the following studies.

Then, we fix

= 0.4, 0.6, 0.8, and 1.0, and vary

, the dimension of the support domain.

Figure 12.37

shows the effects of the shape parameter

on the displacement results for

FIGURE 12.36
Effects of the shape parameter

of MQ on the displacement results.

50%

100%

150%

200%

250%

300%

350%

0.2

0.4

0.6

0.8

20%

40%

60%

80%

100%

120%

140%

160%

180%

200%

0.2

0.4

0.6

0.8

Normalized displacement

q =

- 0.5

Normalized displacement

= 0.5

1238_Frame_C12.fm Page 537 Wednesday, June 12, 2002 5:24 PM

different dimensions of support domains used in RPIM. We found that the best value of

is between 0.6 and 0.8, which supports the finding of 0.7 from Figure 12.36. For an

between 0.6 and 0.8, the support domain

should be approximately 8.4 to 12.0 and 6.3

to 12.0 for q

= −0.5 and q = 0.5, respectively.

Then the variation of displacement with respect to parameter q in the MQ radial basis

function is shown in

Figure 12.38.

When q

= 0 and q = 1.0, the moment matrix may be

singular and the construction of the RPIM shape function may fail, and therefore they
should not be considered.

FIGURE 12.37
Effects of the shape parameter C of MQ on the displacement results for different dimensions of support domains
used in RPIM.

= −0.5

= 0.5

Normalized size of support domain

Normalized displacement

60%

100%

140%

180%

220%

C =0.4

C =0.6

C =0.8

C =0.9

Normalized displacement

Normalized size of support domain

20%

40%

60%

80%

100%

120%

140%

160%

180%

200%

C =0.4

C =0.6

C =0.8

C =1.0

1238_Frame_C12.fm Page 538 Wednesday, June 12, 2002 5:24 PM

Based on the above analyses, the following remarks can be made.

• When q

< −0.7, the results increase smoothly with the increase of q and approach

the accurate solution. The error in this range is too large, and should not be
considered.

• When

−0.7 ≤ q < 0, the results decrease slowly and smoothly, and the error in

the results is no more than 0.5%. A q in this range can be used. The best results
are observed at q

≅ −0.38.

• When 0 < q

≤ 0.7, the results slightly oscillate, but the results are also very good.

The errors are less than 0.4%. A q in this range can also be used.

• When 0.7

≤ q <

1.0 and 1.0 < q < 1.5, the results are unstable, and oscillate violently.

A q in this range should not be considered.

In the previous chapters, a number of results from different studies for different prob-

lems have suggested using q

= 0.98 and 1.03. These two values do not fail in the above-

recommended stable range, and we find they perform well for shells. It would be inter-
esting to try out q

= −0.38 for other problems.

Effect Shape Parameter of EXP Radial Function
There is only one parameter c in the EXP radial basis function, as shown in

Table 5.2.

The

variation of the displacement at the center point with respect to the size of the support
domain is shown in

Figure 12.39.

The results are computed using RPIM-EXP with different

shape parameters, c. It can be seen that the results decrease and approach the accurate
solution with the increase of c in the range 0.003 to 0.03 and the result is almost equal to
the convergent solution when c

= 0.003. The results are very stable and accurate when

parameter c is 0.03 to 0.07. The error increases with the increase of the parameter c. When
c

= 0.12, the error in the result is more than 12%, which is unacceptable.

It is also observed that when a smaller support domain is used, say,

= 4.0 to 7.0, a

small c value between 0.02 and 0.009 performs better. This finding supports the use of

FIGURE 12.38
Variation of displacement with respect to parameter q in the MQ radial basis function.

96%

97%

98%

99%

100%

101%

102%

103%

104%

-1.2

-0.8

-0.4

0.4

0.8

1.2

1.6

Normalized displacement

1238_Frame_C12.fm Page 539 Wednesday, June 12, 2002 5:24 PM

0.003 in the previous chapters where smaller support domains were used. Note also that
using too small a b value incurs a risk of instability.

Figure 12.40 shows the variations in displacement with respect to parameter c of an EXP

radial for a fixed support domain. This figure also suggests using a c between 0.03 and 0.07.

Based on the results from the above discussion, the recommended parameter c is approx-

imately 0.03 to 0.07 and the results have acceptable accuracy, if the support domain is
larger than 5.5 times the nodal spacing. A smaller c, say, 0.003, should be used if a small
support domain is used. In the following examples, c

= 0.03 is used.

FIGURE 12.39
Variation displacement with respect to size of the support domain. Computed using RPIM-EXP with different
shape parameters.

FIGURE 12.40
Variation displacement with respect to parameter c of the EXP radial basis function.

Normalized displacement

20%

40%

60%

80%

100%

120%

140%

160%

c=0.009

c=0.021

c=0.030

c=0.06

c=0.12

c=0.30

50%

100%

150%

200%

250%

300%

350%

0.05

0.1

0.15

0.2

0.25

0.3

0.35

Normalized displacement

1238_Frame_C12.fm Page 540 Wednesday, June 12, 2002 5:24 PM

Example 12.12

Pinched Cylindrical Shell

The pinched cylindrical shell studied in Example 12.9 is used again to test our RPIM.
Configuration of the shell is schematically drawn in

Figure 12.22.

This problem is one of

the critical benchmark tests to evaluate numerical methods in representing inextensional
bending modes and membrane states of stress. All the parameters for the shell are exactly
the same as those in Example 12.9. Because of the double symmetry, only 1

/8 of the cylinder

is modeled by the present method. Figure 12.41 shows the radial displacement at loaded
points, normalized using 1.8248e-5, which is regarded as the “true” value. The convergence
by the present method is excellent in comparison with the results by FEM (Simo et al.,
1989) and by another EFG formulation for shells (Krysl and Belytschko 1996b). It can be
seen for this particular case that the convergence rates of the results by MQ are relatively
slow in comparison with those by EXP.

Figure 12.42

shows the distribution of the vertical displacement along the edge AB

computed using both RPIM and FEM. The FEM solution is obtained by commercial
software (ANSYS) using a rather fine mesh with 441 nodes. It can be seen that the results
of the present method agree very well with the FEM solution.

Example 12.13

Pinched Hemispherical Shell

The pinched hemispherical shell studied in Example 12.10 is reexamined here using RPIM.
The configuration of the shell is shown in

Figure 12.27.

All the parameters used here are

exactly the same as those used in Example 12.10. The only difference is we use RPIM to
solve the problem.

Due to symmetry, only a quarter of the hemisphere is taken into consideration for the

analysis.

Figure 12.43

plots the convergence of displacement at load point B computed

using RPIM-MQ and RPIM-EXP. The results are normalized using the analytical value of
0.0924, and are compared with the results by FEM (Simo et al., 1989) and by another EFG

FIGURE 12.41
Convergence of vertical displacement at point A computed using RPIMs and compared with other available
results.

4.0E-01

5.0E-01

6.0E-01

7.0E-01

8.0E-01

9.0E-01

1.0E+00

1.1E+00

Simo et al., 1989

Krysl and Belytschko, 1996b

EXP

Number of elements/side

Normalized displacement

1238_Frame_C12.fm Page 541 Wednesday, June 12, 2002 5:24 PM

formulation for shells (Krysl and Belytschko, 1996b). The results obtained by the RPIM-
MQ and RPIM-EXP methods are in very good agreement with those obtained using FEM
and the EFG methods for thin shell. It also can be seen that the convergence rate by MQ
and EXP methods is better than that by the EFG method given by Krysl and Belytschko
(1996b).

The distribution of the displacement along the edge AB in the x direction is shown in

Figure 12.44,

where the FEM solution is obtained by FEM software (ANSYS) using a fine

FIGURE 12.42
Distribution of the vertical displacement along AB.

FIGURE 12.43
Convergence of displacement at load point B computed using RPIM-MQ and RPIM-EXP. Other available results
are also plotted for comparison.

-2.0E-05

-1.5E-05

-1.0E-05

-5.0E-06

0.0E+00

5.0E-06

1.0E-05

FEM-ANSYS

EXP

Angle (degree)

Vertical displacement

20%

40%

60%

80%

100%

120%

Simo et al., 1989

Krysl and Belytschko, 1996b

EXP

Normalized displacement

1238_Frame_C12.fm Page 542 Wednesday, June 12, 2002 5:24 PM

mesh with 503 nodes. It can be seen that the results of the MQ and EXP methods agree
with the FEM solution very well.

12.3.3

Remarks

In this section, RPIM methods using both MQ and EXP radial basis functions are examined
for analyzing 3D shell structures based on the stress-resultant shell theory, which has only
five degrees of freedom assigned to every point of the shell.

Both the geometry of the neutral surface of the shell and the field variables are interpo-

lated using RPIM shape functions, which are constructed using nodes in the support
domains. The RPIM shape functions possess the Kronecker delta function property. There-
fore, the essential boundary conditions can be easily implemented with ease as in FEM.
Some important parameters on the performance of the present method are investigated
in detail. The findings are summarized as follows:

• The dimensionless shape parameter

should be around 0.7.

• The shape parameter q should be in the range of

−0.7 ≤ q <

0 and 0 < q ≤ 0.7.

• The shape parameter b for EXP radial basis function should be approximately

0.03 to 0.07, if the support domain is larger than 5.5 times the nodal spacing. A
smaller c, say, 0.003, should be used if a small support domain is used. However,
there is a risk of instability in using too small a value for c.

• RPIM-EXP outperforms RPIM-MQ for our benchmark examples.
• Our study is still ongoing, and currently is investigating the effects of adding

polynomial terms. We hope adding polynomials can improve the accuracy and
reduce the sensitivities of the results to the shape parameters.

FIGURE 12.44
Distribution of the displacement in the x direction in the hemispherical shell along AB.

0.0E+00

1.0E-02

2.0E-02

3.0E-02

4.0E-02

5.0E-02

6.0E-02

7.0E-02

8.0E-02

9.0E-02

1.0E-01

3 0

FEM

EXP

Angle (degree)

Displacement in the

x direction

1238_Frame_C12.fm Page 543 Wednesday, June 12, 2002 5:24 PM

12.4 Summary

In this chapter, we formulated the EFG method for both thin and thick shells and RPIM
for thick shells. We have shown through a number of benchmarking problems that they
work very well. In all cases, the methods performed much better than standard commer-
cially available FEM packages. There are cases where some fine-tuned FEM results
reported in the literature performed better than some of the MFree results. In any case,
the author is confident that the MFree method is a better choice for simulating shell
structures.

We formulated MFree methods here based only on Galerkin formulation; however we

can also use Petrov–Galerkin formulation. This is an ongoing project of the author’s
research group.

1238_Frame_C12.fm Page 544 Wednesday, June 12, 2002 5:24 PM

Document Outline

MESH FREE METHODS: Moving beyond the Finite Element Method
- Contents
- Chapter 12: Mesh Free Methods for Shells

Wyszukiwarka

Podobne podstrony:
C12 4
1238 C09
C12
1238 C04
C12 6
1238 C01
(budzet zadaniowy)id 1238 Nieznany (2)
1238 C07
C12 5
1238 C06
C12 2
1238 C08
1238 C03
1238 C13
1238 walc chemiczny golec uorkiestra 377UBYFFMIZM5RAELF3U2F274QSAIVB6WEYLEYY
c12, SGSP, SGSP, cz.1, hydromechanika, Hydromechanika, instrukcje stare
Dz U 08 201 1238 Warunki Techniczne zm
C12 1

więcej podobnych podstron

1238 C12

Document Outline