- Email: [email protected]
Meshless methods for shear-deformable beams and plates
Computer methods in applied mechanics and engineerlng Comput.
ELSEVIER
Meshless methods
Methods
Appl.
Mech. Engrg.
152 (1998) 47-71
for shear-deformable
beams and plates
Brian M. Donning ‘, Wing Kam Liu *j2 Dept. of Mechanical Engineering,
Northwestern University, 2145 Sheridan Rd., Evanston, IL 60208, USA
Presented at Professor J. 7: Odenk 60th Birthday Symposium Received 13 February 1997
Abstract A meshless method IS developed to analyze moderately thick and thin structures using Mindlin-Reissner theory. A uniform discretization is used to allow for efficient integration and for the shape functions to be written explicitly. Irregular boundaries are modeled in a straightforward manner. An unmodified displacement-based Galerkin method is used; full integration is used to evaluate all energy terms, and convergence is independent of the thickness. Shear and membrane locking are completely eliminated pointwise at the interpolant level using cardinal splines. The continuity of the splines chosen results in continuous stresses. An extension to general meshless methods is given. Beam and plate examples show the accuracy of this method for coarse discretizations.
Introduction
To make nmneri~cal analysis of moderately thick/thin structures which exhibit shear deformation more practical, a method. needs to have the following properties: complete absence of numerical locking, accuracy for coarse discretizations, continuous stresses, automatic discretization, accurate geometry modeling, and numerical efficiency. Successful approaches to locking in finite element methods (FEM) have focused on modified variational principles. Assumed strain methods [5,7], hybrid stress methods [4], and reduced and selective integration methods [11,14] have all been used to ameliorate locking and are related to a general class of mixed methods [3]. However, these methods often lack robustness and accuracy on coarse meshes. Meshless methods inherently satisfy many of the problems associated with FEM. Different forms of meshless methods have emerged which differ in the construction of the interpolant, but all maintain several common properties. Reproducing kernel particle methods (RKPM) modify interpolants to satisfy consistency requirements via a correction function [9,10], element free Galerkin methods use a moving least squares procedure to construct interpolants [6], and H-p Clouds use a partition of unity and a set of basis functions to construct interpolants [8]. In each method, interpolation fields yielding continuous stresses are easily constructed, a FEM mesh need not be created to discretize the problem, and the geometry is arbitrary since there are no element edges, but only a set of points. However, efficiency and locking have limited the use of meshless methods to model thin structures.
*Corresponding author. ‘Graduate student in Theoretical and Applied Mechanics, Northwestern University. *Professor of Mechanical and Civil Engineering, Northwestern University. Email: [email protected] 0045-7825/98/$19.00 c:, 1998 Elsevier PZZ SOO45-7825(97)00181-3
Science
S.A. All rights reserved
48
B.M. Donning,
WK.
LidComput.
Methods
Appl.
Mech. Engrg. 152 (1998) 47-71
A new method is proposed to maintain the desirable properties of meshless methods, but also eliminate locking and improve efficiency to the level of FEM. The displacement-based Galerkin method for the principle of virtual work is used as the variational framework for this method (linear, static, elastic). Cardinal splines are used as the shape functions, and their interpolation properties greatly enhance the accuracy of the method for coarse meshes. The shape functions are given analytically, and integration is exact in the interior of the problem domains; background integration cells do not need to be created from a FEM mesh. Various examples in linear elasticity, including straight beams, curved beams, and plates of various geometries, show the effectiveness of this method. Also, in the context of general meshless methods, locking free interpolation fields are identified for beams and plates. RKPM is used to demonstrate this. An extension to FEM is also ma’de, using the interpolant for the proposed method as the local interpolant-the resulting elements are free of locking and have continuous stresses.
1. Mindlin-Reissner
beams via mesbless
methods
1.1. Beam model
The beam is isotropic, linear, elastic; the kinematic assumptions and strain-displacement the bending of straight beams are as follows: axial displacement transverse displacement
U(X,Z) = -ze w(x,z) = w Ex
=
u,x =
relations for
-ze,x axial strain shear strain
YZX= w,, -0
The Galerkin form for a rectangular straight beam is, neglecting body forces: SWex’ = gwbe”d
awbend
Ebt3
= -
12 6W shear = Gbkt
+
swshear
@,x dx J8,x -t3)(6w+ -60)d.x J w,, L
0 0
where t = thickness, b = width, k = shear correction factor, E = Young’s modulus, G = shear modulus and L = length of beam. 1.2. Shear locking Shear locking develops in formulations which cannot represent a state of zero shear in the thin beam limit independent of the mesh. This is exhibited as a violation of the Kirchhoff assumption in thin beam theory, stated as w,, - 8 = 0. When a beam is thin, most of the energy of deformation is due to bending. If a formulation does not admit non-trivial zero shear solutions, then most of the energy of deformation goes into shear deformation, resulting in very little bending-hence locking. This can be seen from the variational form; since the coefficient of the bending term is of order t3 and the coefficient for the shear term is of order t (fix L), the latter becomes much larger as t -+ 0. Since the variational form gives minimum energy solutions, the large shear term coefficient forces the minimization of the shear deformation. If the interpolant cannot physically represent such a state, spurious shear develops, giving a grossly inaccurate representation of the deformation state of the beam. To rid the formulation of locking, the interpolation of 8 must be able to approximate the derivative of the lateral displacement w arbitrarily well. The degree of satisfaction of this requirement determines the relative amount of spurious shear energy present and thus the degree of shear locking. Since meshless interpolants can be created from a wide variety of functions, many different combinations potentially
13.M. Donning,
WK. Liu/Comput.
Methods Appl. Mech. Engrg. 152 (1998) 47-71
49
could be used to alleviate locking. However, once an interpolation field that alleviates locking is found, it may not exactly satisfy the Kirchhoff constraint pointwise. This deficiency, while possibly giving accurate displacement results, causes oscillations in the stress fields to develop. The next four secl:ions outline the construction of locking-free interpolation fields for straight MindlinReissner beams. 1.3 Approximation
of spline derivatives
Spline functions are chosen here because of their many interpolation properties [2] (see Appendix A). A spline of order p is denoted by S’. The spline family has an additional property that S$ can be exactly represented by SP-’ for uniform spacing when the splines are spaced such that their piecewise components are aligned. S:(x) = Sp-‘(x + m) - Y’(x
- m)
For convenience, S::(X) is centered at x = 0, and 2m is the spacing between consecutive splines for both fields. Fig. 1 shows the alignment of quadratic and cubic splines fields, and the ability of the quadratic spline to exactly represent the derivative of the cubic spline. The vertical lines denote the division of the piecewise parts of both spline fields. Of many functions which could alleviate locking by approximating the Kirchhoff constraint, splines are very practical since they need relatively low integration orders and can easily give continuous stress fields. Fig. 2 shows the ability of the present formulation to model thin beams. A cantilevered beam with an end point load is modeled using 7 degrees of freedom. Conditioning problems start to degrade the solution at L/t = 106. Whereas traditional FEM solutions to locking in a displacement-based formulation have relied on reduced integration or assumed strain techniques, the present technique uses full integration, and all displacements and strains are related by a derivative at every point in the problem domain. FEM could use a similar technique, but the Co continuity inherent to standard FEM interpolants would cause the continuity for 8 to drop to C-‘. This would create Dirac delta functions in the stress interpolation. In addition, in the interpolation schemes given here, all field variables are at least C’, meaning stresses will be continuous.
0.6
1.6
\ \ \ ! \ \
-0.2-
-0.4--
2 Quadratic Splines
,
-0.6- - -Quad. SqlineDiWlerem
I.'
,
. .
..I
.
. . .
.
. . .
.
. . .
.
..C
P gi.1-
I'
\
,.-..I
x +1.3E e Sl.2‘
,'
.ji I
'i \ \ \
- - 1 CubicSpline
. ...I
/ .I
.P l-
-J
.- Deriv.of CubicSpline -0.B' -2
-1.5
-1
-0.5
0 x
0.5
1
1.5
Fig. 1. Spline alignment Fig. 2. Absence
os-
2
IO0
and derivative
of locking
...._I 10'
".
c ___ IO'
approximation.
in the thin beam
limit.
IO3
u
.,_c ..._.f 10' loa
. ..._a 10'
.....10'
50
B.M. Donning,
WK. Liu/Comput.
Methods
Appl. Mech. Engrg. 152 (1998) 47-71
1.4. Interpolation via cardinal splines
To maintain consistency on a finite domain, the finite domain can be ‘cut out’ of an infinite uniform interpolation field (Fig. 3). Only the parts of the interpolants which overlap with the problem domain are used. The following method is for uniform spacing, yet it still satisfies the interior consistency at the boundary for irregular geometries by extending the mesh outside the physical domain. W(X) =
C
8 (X)
C Sp-’(X -
=
S”(X - xi)di
Xi)hi
Here, di and hi are the unknowns, and Xi is the position of the center of each spline. A cardinal cubic spline for the w interpolant and a cardinal quadratic spline for the 8 interpolant were used in all the beam examples. Three point Gauss quadrature in each nodal interval is needed for exact integration (based on the w nodal spacing). If higher continuity is desired, higher-order splines can be used in this same setting, as long as the 8 field is one order lower than that of w; higher-order quadrature and a larger bandwidth in the system equations would be the tradeoff for the higher continuity in the spline family. 1.5. Interpolation for any C’ meshless method
Since these meshless methods have high continuity, the interpolation field for the rotation of the midline is chosen by differentiation of the field for the lateral displacement. In meshless methods suitable for irregular discretizations, the finite domain includes only those shape functions whose centers are included in the domain; subsequently, the interpolants must be modified to satisfy consistency near the boundary and in the interior for irregular spacing. Using a cubic spline as the interpolant and a uniform discretization, RKPM only modifies the splines on the boundary to maintain consistency (Fig. 3). The approximation level of the modified shape functions near the boundary is decreased due to the elimination of shape functions whose centers are outside the domain, reducing the consistency for the global interpolant. The interior shape functions (chosen to be the cubic spline) have third-order approximation character, whereas on the boundary, the modified splines only have first-order
lnltnite 0.8
spline interpo!mion 1
;i2; -1.5
-1
-0.5
0
0.5
1
1.5
2
Modifiedboundaryinterpoiation
Fig. 3. Two different
boundary
interpolations
for a cubic spline S3.
E1.M.Donning,WK. Liu/Comput. MethodsAppl. Mech. Engrg. 152 (1998) 47-71
character. Also, the modified functions in general will no longer be polynomials for non-uniform so integration of these modified functions must be done carefully. W(X) =
C
51
spacing,
Nil:X)di
i
8
(X) =
C
Ni,x(x)hi
i
Here, N, is an RKI’M or EFG, etc. shape function. The interpolant Ni,X is continuous for C’ or higher fields for W(X). However, this shape function is no longer non-negative, and cancellation can cause numerical problems [l]. 1.6. Interpolation via spline based FEM Only a suitable mapping field must be chosen to implement a locking-free spline based element. Here, the mapping is assumed to be sufficiently continuous. x = f(s)
($)-lhi
f3(X)=CN,(t)
i
Here, the properties
of cardinal splines are used:
N&L!) = SP(S - &)
These are locking-free
interpolants,
since they exactly satisfy the Kirchhoff constraint.
1.7. Bending of curved Mindlin-Reissner
beams
Curved beams exhibit both shear and membrane locking. The following method follows the same ideas as the straight beam. The beam considered is a circular quarter arc (rectangular cross-section), loaded in its plane of curvature (see Fig. 4). The midline of the beam is parameterized with respect to the global Cartesian coordinates xi. The curvilinear coordinates c& are orthogonal. x1 =x=Rcos-Sx2 =y
=Rs_inL
a1 = t
R
a2 = s
R
a3 = z
x3 = z
A point on the midline is given as
The outward normal to the midline is n(“)=cosSe1+sinSe2 R
R
B.M. Donning, WK. Liu/Comput.
Methods Appl. Mech. Engrg. 152 (1998) 47-71
-RR Fig. 4. Geometry of circular arc.
A point
the beam is given by
in
r(O) = r!)(s)
+ tn(O)(s)
where (s, t) are the coordinate
parameterization.
1.8. Continuum mechanics The strains for the curved beam are derived from displacement components defined in the curvilinear system which are analogous to the Mindlin-Reissner kinematic assumptions used in the straight beam
PI. The base vector along (Y*is
i3rp -=thf2
(&O) =
a,(y)
II-I
-sin Se, + cos Se2 R R
da2
The base vectors g, in the curvilinear coordinates g1
are
dr co) = -da’ = J&O) dr co)
t3-r’
dd”)
g,=g=da2+tda2=
t +1
( 1 u(oj
E
g3 = e3.
The volume transformation
between Cartesian and curvilinear is
dV = &dsdtdz where ,,& = g, ’ (g3 x gI) = (t/R + 1). The position of a point after deformation r =
r(O)
+ u
where (u, 0 and w are the displacement u =
(u+te)a(“)+wn(o)
=
variables) (u+t~)g*+w&+a3.
is
B.M. Donning,
WK. Liu/Comput.
Metho&
Appl.
Mech. Engrg. 152 (1998) 47-71
53
The displacement ;Z is not defined because it is left to satisfy the plane stress conditions. The metric tensor gii = gi . gj is
t
1
0
[gijl= 0
0 2
() +1
0
0
1
3
0
1
I
The strains are derived as follows. Let u = uigi and the deformed Gi . Gj. Then define the strain as
base vectors Gi = r,i with Gij =
eij = k (G,j - gij) = i (ui(j + Ujii + UkliUkIj) and linearize in the displacements Eij
=
k
(gikUk
j
+
gkjL&Ji)
where ui lj = u:j + { Gj}r_.8. The Christoffel symbol of the second kind is denoted by {ii}. The following two expressions for the linearized strain in curvilinear coordinates do not follow the summation convention.
2Eij = yij = g,i g
+ gjj g
i # j
The nonzero in-plane strain components Es s! E22 =
( I( J jr
+l
are
u,, +to,, + w
R>
‘y,$E 2E,2 = tv,, + 8 - ;. The following substitution is made to separate the rigid body rotation (u/R) from the total rotation (0) of the midline normal:
The strains are now rewritten in terms of U, w and 4. Yts = w,s + 4
The membrane Em=
strain is
(;+1) (u,,+i)
The bending strain. is
For the constitutive relations, a local Cartesian system coincident with the orthogonal curvilinear system ci is used, since the plane stress assumptions through the thickness of the beam are aligned
54
B. M. Donning, WK. Liu/Comput.
Methods Appl. Mech. Engrg. 152 (1998) 47-71
with the local Cartesian system. The material is linear and isotropic. The stress, strain and displacement variables in this local system are denoted with an overbar (no sum). ii, = 6
ui 1
Since ~~j~Ejj = &jSeii, the constitutive a22 =
relation will only be developed
in the local Cartesian system:
E~22
EEZ2= -
g22
1.9. Variational form Neglecting body forces and using the curvilinear coordinate
J Js( Jl( &j
SEij dV
zz
system:
swext
V
EQ
8~
+ 2Gkr12 8~
k22j2
E~22
8~22
+ 2Gke2, Se2, g22
g22 + Gb2
h12
g22
k2212
JJ(E(&)’
&stGk(-&j
&
&ds
dt = Z5Wex’
)
ds dt = 8WeX’
1
y&rs)
dsdt=sW”“’
where the width of the beam is 1 unit. Using the interpolations for u, 8 and w, and integrating analytically in t, the final matrix form is found in terms of the arc-length parameterization s.
1.10. Shear and membrane
locking
Shear deformation is important to model moderately-thick to thick beams. Most of the strain energy in a thick beam goes into shear deformation, whereas in a thin beam, most energy goes into bending deformation. Shear locking appears in thin beams if the interpolation does not allow bending without the generation of spurious shear strains. Membrane action comes from the coupling of bending deformation and membrane deformation in curved beams. Membrane locking [18] appears when the coupling is too strong-the interpolation does not allow for a state of bending without corresponding membrane deformation. This is not necessarily non-physical, since curved beams are supposed to show this coupling; however, since they are dependent in a specific way determined by the interpolants, there is no flexibility to allow different amounts of coupling which appear in different beam configurations (esp. thin beams). Let the beam become very thin so membrane and shear locking become issues: t/R -+ 0 : fix R, t + 0. Then the variational form for the virtual work becomes
B.M. Donning,
WK. LidComput.
Methods Appl. Mech. Engrg. 152 (1998) 47-71
55
shear membrane
+Et’($+&,)
dsdt
(%+a+,,
. bending
The coefficients of the terms in the minimum energy functional act as penalties [14] on the membrane and shear energy terms: E>Gk>>Et=
The appearance of locking comes from the inability of the interpolation to represent certain states of strain. When a curved beam becomes thin (t/R + 0) , almost all of the deformation comes from bending, and two constraints appear. (1) ‘yIs= w,, + $ = 0
shear constraint
(2) 6, = u,s .t w = 0
membrane
R
constraint
Traditional finite element interpolation fields only satisfy these constraints on coarse meshes with a zero solution for all fields, hence the term locking. Using higher continuity fields, both membrane and shear locking can be avoided completely (using exact integration) for all meshes with an appropriate choice of interpolants. 1.11. Interpolation via cardinal splines
The two constraints from the previous section define the relationship between the interpolants for u, w and 4. For comtraint (1) to hold, 4 must be able to approximate the derivative of w. For constraint (2) to hold, w must be able to approximate the derivative of u. u=
cS”(s -
w =
c
#I =
c sqs - Sj)rni
sj)d,
sp-’ (s - s;)h;
The properties of the spline family on a uniform discretization are exploited for this choice of interpolants. Since one member of the spline family can exactly interpolate the derivative of the next higher spline if the piece-wise components of the splines are aligned, one only needs to choose three ‘consecutive’ members of the spline family for u, w and 4 in descending order. The choice of the order for u automatically defines the order of w and 4. A quartic spline is the lowest-order spline for L(which results in continuous stresses (p 3 4). 1.12. Interpolation for any C2 meshless method
In the context of a general meshless method, a C2 interpolation is needed for a set of nodes for u (a C3 or higher field is needed for all stresses to be continuous). Then the shape functions for w and 4 are derived simply through differentiation of the shape functions for u.
56
B.M. Donning,
u = C
WK. Liu/Comput.
Methods
Appl. Mech. Engrg. 152 (1998) 47-71
Ni(S)d;
1.13. Interpolation via spline based FEM Similar to the straight beam, a spline based locking-free element can be constructed. s = f(t) ~(3) = C Sp(t - &)di
W(S)
C Sp-‘(5 - &)
=
($) -’hi
i
These are locking-free interpolants, since they exactly satisfy the shear and membrane constraints. Because the coordinate mapping creates an extra term in the derivative, twice as many degrees of freedom are needed to interpolate 4; this can be avoided by using the true derivative as the shape function, but conditioning problems may appear.
2. Mindlin-Reissner
plates via meshless
The classical Mindlin-Reissner setting and removed.
methods
plate bending model is used. Shear locking is identified in this 2D
2.1. Plate formulation A static, linear, isotropic plate is formulated as follows [ll]. The 3D plate domain V is reduced to a 2D domain by defining all variables with respect to the plate’s midsurface A. V =
{
(x,y,z)
The displacements u1 =
-Z&(&Y)
u2
=
-ZW,Y)
u3
=
W(X,Y)
E R31z E
-;, 1
;
are defined as
1
, (x,y) E A c R2
For an isotropic material (sum on k) uii = 2peij + A6ijekk (i, j, k = 1,2,3)
>
B.M. Donning,
WK. Liu/Comput.
Methods
Appl.
Mech. Engrg. 152 (1998) 47-71
57
The plane stress condition imposes CT33=
0 +E33=-hf2P
h(
(1)
El’ + E22)
The stresses are then simplified by eliminating ~~~~The strain-displacement Eij
=
k
(Ui,j
--
relations are
Uj,i)
This implies l33 = w,~ = 0, but this contradicts the plane stress assumption (1). However, this discrepancy does not enter into the variational form since the virtual work done by the true deformation across the thickness is zero (8~~~r33 = 0). The variational form is (a, p = 1,2):
J
(SE,~ aUp + 2k Sea3 uG3) dV =
V
SwT,dS+C6wF, ss
i
Here, Fi are the point loads normal to the midsurface, T, are the applied tractions midsurface and k is the shear correction factor. Body forces are neglected.
normal to the
2.2. Shear locking In this plate formulation,
two shear locking constraints appear as t + 0 (c~ = 1,2):
w>, -l&=0 Since these constraints are each associated with an independent direction, they can be treated separately in succession. Once the interpolation for w is chosen, applying each constraint in turn generates 13’and 02.
2.3. Interpolation via products of cardinal splines Splines for highker dimensional approximation are derived via a product rule from 1D splines. As a consequence, the higher dimensional splines behave exactly like the 1D splines in each direction. In 2D: P(x,y)
= SP(x)S4(y)
The resulting functions have square support and consist of piecewise polynomial patches. The consistency properties and ability to approximate derivatives of higher-order splines remains intact in each direction. Combinations of different order splines in each direction can approximate mixed and other derivatives of higher-order spline products.
sjJ$, y) = SP-‘(X + m)S4(y) - Y’(X
- m)Sq(y)
Sp::(x,y) = SP-‘( X + m)Y-’ 0, + m) - s-’
(X + m)W
- 9-l (X - m)S@ (y + m) + 9-l (x - m)W
(y - m) 0, - m)
These propertie:s are all based on rectangular uniform collections of spline functions with nodal spacing m in each direction. The ability of such a collection to serve as an interpolation field for irregular geometries relies on the arbitrariness of boundaries in meshless methods. Meshless methods usually choose some polygonalization of the boundary nodes to approximate the problem geometry. This is purely arbitrary, since the interpolant is not defined by the polygonalization, but merely truncated by it. For the uniform splines, the exact geometry can be used (likewise in any meshless method) in an analogous manner. Qualitatively, a uniform spline interpolation field is created which extends outside the boundary; the exact geometry is superimposed on the interpolation field and the interpolant which intersects with the domain is used in the interpolation field (Fig. 5). This maintains all the properties of the splines on a finite domain without modification near the boundary. This results in the inclusion of some functions whose nodes (center of support) are exterior to the problem domain. However, only
B.M. Donning,
WK. Liu/Comput.
Methods Appl. Mech. Engrg. 152 (1998) 47-71
2D Problem Boundary
l
0
0
Fig. 5. 2D uniform
discretization
Interior nodes
Non-zero shape function region for center node
of irregular
domains.
the nonzero region extending into the domain is used in the Galerkin procedure. The accuracy of the boundary is determined by the accuracy of the integration scheme around the boundary. Interpolation via cardinal splines: 8, =
c
02 =
C Sq-’(J
W
C Sp(X - Xi)S’O, - yi)mi
=
sp-’ (x - x;)Y
0, - Yi)hi
- yj)S’(X - Xi)gi
Integration of these uniformly spaced splines is exact in the interior for an appropriate Gauss quadrature rule. Since the interpolation field consists of piecewise regions of polynomials, each such region can define a cell for the integration. Product rules of 1D Gauss quadrature formulas integrate these functions exactly in the interior. Near the boundary, the problem of integrating a general multidimensional domain appears. The interpolant is still piecewise polynomial in character, but the area or volume of the cell is no longer rectangular. Subdivision or polygonalization of boundary cells reduces quadrature error, and exact Gaussian product rules for right triangles are easily constructed, but many more quadrature points are required for the right triangle than for a rectangle. 2.4. Interpolation for any C’ meshless method Interpolation
6
=
via any sufficiently smooth method:
CNi,y(X,Y)gi i
B.M. Donning,
2.5. interpolation
WK. LidComput.
Methods
Appl. Mech. Engrg. 152 (1998) 47-71
59
via spline based FEM
A sufficiently smooth coordinate x = f(5777)
mapping in 2D is chosen:
Y = g((, 77)
J=axay__axdy 86 877 877 a,$
e,(x) =
3 :~Sq(7] [
77i)Sp-‘(5 - ~~)~aY
(1)
hi
-
C
Sp(~ -
E)Sq-‘(rl -
77i) ~
h12)
i
i
Because the coordinate mapping creates an extra term in the derivatives, twice as many degrees of freedom are needed to interpolate 61 and 02; this can be avoided by using the true derivative as the shape function.
3. Essential bounldary conditions The displacement boundary conditions were enforced via Lagrange multipliers and a modified Lagrange multiplier [15] method in the examples given. The summation convention is used in this section. 3.1. Collocation
Collocation was tested in the two straight beam problems considered in the numerical examples section. The performance of collocation on coarse meshes was very poor, since a collocation equation replaces an equation from the variational form. The lost information has a profound effect on the solution near the boundary. The two following methods were implemented to impose boundary conditions while retaining all of the stiffness equations. 3.2. Lagrange snultiplier method
This approach adds constraint equations to the system of equations, coupling degrees of freedom via Lagrange multipl:iers. A term is added to the internal work of the variational equation as follows: @int =
wint
_
hi (ui - Ui) dS .%‘ The Lagrange multiplier fields hi are arbitrary fields on the plate boundary where the displacement boundary conditilons Ui are to be imposed. This approach has a few disadvantages. The constraints augment the size of the system matrix; the form of the Lagrange multipliers renders the system of equations indefinite, requiring the use of a special solver for large systems;. the conditioning of the equations (in the sense of the LINPACK condition estimator) is large for every value of t. Also, the choice of functions for the hi fields, while still providing good satisfaction of the boundary conditions, affects the interior solution markedly in coarse meshes. However, the arbitrariness of the Lagrange multiplier fields potentially allows for more accurate satisfaction of BC’s on non-smooth boundaries. I
B.M. Donning,
60
WK. Liu/Comput.
Methods
Appl. Me&.
Engrg. I52 (1998) 47-71
3.3. Modified Lagrange multiplier method
The modified Lagrange multiplier method [15] solves some of the problems associated with the classical Lagrange multiplier approach. Here, h; is replaced with its physical counterpart (the equivalent reaction tractions of the properly constrained system), called T,R, arising from the Euler-Lagrange equations of the variational principle. In this form, the constraints do not add any additional degrees of freedom and the arbitrariness of the choice of Ai is no longer an issue. Also, the positive-definiteness of the original system is retained [lS]. @in, = jyint
_
J
Tr(ui - U;) dS
s,,
Here, nj is the outward normal to the boundary and cij is related to the displacement degrees of freedom via the constitutive and strain-displacement relations. Since the procedure modifies stiffness terms of the boundary associated shape functions, the accuracy of this method depends on the relative accuracy of integration of both the stiffness terms of the boundary interpolants and the boundary integration for the essential BC term. The numerical behavior of this method is very good for coarse meshes for the plate problems solved. The conditioning of the system matrix depends on t as is expected for a pure displacement formulation. For Co geometries, the accuracy of the BC satisfaction of this method greatly deteriorates near the non-smooth points relative to the accuracy on the smooth regions; this inaccuracy shows up as stress oscillations.
4. Numerical results 4.1. Straight beam Two thin beam problems are solved: problem 1 is a cantilever beam with a point end load (Fig. 6), and problem 2 is a clamped-clamped beam with a center point load (Fig. 7). Lagrange multipliers are
=:[[I 0
2
1
3
4
2 X 1O-6 *r-
'---'
Fig. 6. Cantilever Fig. 7. Clamped/clamped
7 ii_
9
10 .
- - w,x w,xexacl
. . '\ '..
0'
beam with end point load (7 DOFs). beam with center
6
._
/' / E
5 x
point load (9 DOFs).
\
,
B.M. Donning, WK. LidComput.
Methods Appl. Mech. Engrg. 152 (1998) 47-71
used to enforce essential boundary conditions. Both problems use the following set of parameters, w is interpolated via cubic splines, and ~9is interpolated via quadratic splines: V = 0.3
L = lo4
E’ = lo5
61
and
t=b=l
The exact thin beam (Euler-Bernoulli)
solution for problem 1 is
w(x) = mQ (c-X3 + 3Lx2 \ > and for problem 2 H(x - L/2)
x--)I L
3
(
3
2
where Q is the magnitude of the point force and the Heaviside unit step function H(t) is defined as zero for t < 0, 1 for t 3 0. The exact solution for problem 1 is reproduced using the minimum number of degrees of freedom to discretize the beam since the approximation order of the interpolant is third order, and the the straindisplacement equations are satisfied everywhere. Problem 2 is piecewise cubic and is reproduced exactly by a piecewise cubic field. 4.2. Curved beam The 2D plane stress elasticity solution for the quarter circular arc with an end point load (Fig. 8 ) is taken from [13]. 20 u = -Ft?
D(1-u)lnr+A(1-3u)r2+B(~~V)
case +
+Ksin8+LcosB I
u=z/jsin(j - Cos[A(5+v)r2-D(l-v)lnr+ E
-L sin 8 -t D’y
E
”
B(l+ v) r2
1
+Kcos8-LsinO+...
cos 8
X
Fig. 8. Geometry
for exact solution.
62
B.M. Donning,
A=2N
WK. Liu/Comput.
P
Methods
Appl. Mech. Engrg. 152 (1998) 47-71
N = a2 - b2 + (a2 + b2) In i
D = -$(a2+b2) B(l + v)
r2
D(l-v)lnr+A(l-3v)r2+
I The solution in Fig. 9 is in the coordinates of Fig. 4. The solution uses the same material parameters as for the straight beam with R/t = 2/7r x 104. The solution breaks down for R/t 3 105, possibly because the conditioning becomes very poor in the thin limit. The accuracy for this problem is uniform in t. Full integration is used (4 pt. Gauss quadrature) and a quartic spline is used for u, a cubic spline is used for w and a quadratic spline is used for 4. A uniform spline mesh is used and essential BC’s are imposed via Lagrange multipliers. Here the solution is no longer a low-order polynomial and the splines used give good agreement with the exact solution for the prime displacement variables and their derivatives with a minimal number of nodes. 4.3. Plates In each example, the whole plate is modeled and, where possible, symmetric discretizations and quadrature meshes are used (for symmetric problems). A product of cubic splines is used to interpolate w and corresponding products of cubic and quadratic splines are used to interpolate 8i and 0,. The boundary conditions employed are as follows: Displacement CL SS
BC’s
Clamped Simply Supported
W = 0,8, = 0, 0, = 0 w= 0
Traction BC’s C Concentrated center load U Uniform lateral load
El El --._ F-J 4X1O-6
0.25
I
El
0.2
, /’
- wcmlp.
w exact
0.15
i 0.1
,r
,
I’
,
22
I’
0
*’
I’ 1.’
/’
0
I
5000
00
10000
4.1
-0.15
‘\
0
‘\ \
‘\
I -
0
-.
mmp.
w.0
oxacl
loo00
-1
‘\
-
\\
‘\
“3
‘\\
-2
\
El ‘W
-.
5
\
u wmp.
‘U
-0.2
W.8
,XI04
-0.05 =
-El
5ow s
s
0
‘1
/.’
,//
,’
0.05
- -.._
_I
/’
-3
axact
exed
amp. ‘\
‘\
‘\ 1oMx)
-4 : 0
IWO0
s
s
Fig. 9. Numerical
solution
for a curved
beam
(15 DOFs).
B.M. Donning,
WK. LitdComput.
Methods Appl. Mech. Engrg. 152 (1998) 47-71
63
For a thin square plate, only the uniform spline method was used, since RKPM is closely related for rectangular geometries. The convergence of the center displacement is shown in Fig. 11. A symmetric uniform discretization is used (Fig. 10). The following parameters were used in the problem: L = 10
t ==1 x 1o-4
E = 10’
v = 0.3
k = 1.19
To compare with finite elements, the convergence plots are reported in equivalent linear 4-node elements for a quarter of a uniform square mesh using the following relation:
(2) This relation relates the degrees of freedom (DOF) of a given discretization of splines to the number of elements for a quarter symmetric mesh of linear elements containing the same number of degrees of freedom. The uniform spline method is used to solve a circular plate (Fig. 14). Oscillations in the convergence plots are due to the inaccurate integration near the boundary (rectangular Gauss product rules were used, and the integration weights of points outside the domain were set to zero). A symmetric discretization (with respect to a quarter of a circle, Fig. 12) and integration mesh are used. The following parameters are used: R=5
t = 1 x lop4
v = 0.3
E = lo5
k = 1.19
To compare with finite elements, the convergence plots are reported in equivalent linear 4-node elements for a quarter of a uniform circular mesh using the following relation: equiv.elementsfor$mesh=
(DoL-l)
- (y)
(3)
The moments for a thin circular CL-U plate (676 DOF) are plotted with the exact solution [16] in Fig. 13. The thin circular plate is also solved using RKPM (Fig. 16). A non-uniform discretization is used (Fig. 15). A circular cubic spline window is employed as follows:
Thin Square Plate 0
I.&-
~_________-+_________+-___-____~ 8
Q
I
I
I
I
I
I
,
1.01 -
zB
_______ 1
r
D
go.99 b
6 o
0
I
,
I
I
0.96
= ; g 0.97 -
b_________~_--______~___-______6
Q 0.96 -
0
Q
Q
0.95 IO"
Q
0
Fig. 10. Uniform Fig. 11. Convergence
10’
Q
equiv. 4-node elements
spline mesh for a square
plot uniform
spline method-thin
plate. square
plate.
10'
B.M. Donning, WK. LidComput.
64
Methods Appl. Mech. Engrg. 152 (1998) 47-71
3
2
1 t i 0 5
0
5
-1
i -2
-3
-4 0
0.2
0.4
0.6
0.8
r/R
Fig. 13. Moments
along the x-axis for a thin circular
plate
(CL-U).
Thin Circular Plate 1.02
I
,
1.01
‘\
I’
0.94
‘\
’ 10’
102 equn. 4-m&
Fig.
14. Convergence
for uniform
elements
sphnedircular
plate.
1
KM. Donning,
WK. Liu/Comput.
Fig. 15. RKPM Fig. 16. Convergence
65
Methods Appl. Mech. Engrg. 152 (1998) 47-71
mesh for a circular of RKPM
plate.
for thin circular
plate.
o.ooc25
computed analytical
----.
5e-a5
0
0.2
0
0.4
0.6
0.8
1
r/R
Fig. 17. Thick circular
S3(X
-Xj,Y
-,Yj)
=
plate profile.
S3 (f)
where r =
d(X -
Xj)”
+
(Y -
Yj)’
The quadrature mesh is based on a uniform pattern of cells. A FEM background create the cells. The following parameters were used: R=5
t =
1 x 1o-2
v = 0.3
E = lo5
a =
1.5
mesh is not used to
66
B.M. Donning,
WK. Liu/Comput.
Methods Appl.
Mech. Engrg. 152 (1998) 47-71
This plate is thicker than for the uniform spline problem because conditioning problems were encountered at a lower L/t ratio. A clamped thick circular plate with a uniform load is modeled with a uniform mesh. The computed solution for the displacements are plotted against the exact solution [16] in Fig. 17. The following parameters were used: R=5
t=2
DOF = 1149
Using Eq. (3), this solution corresponds to 90 4-node linear elements on a quarter symmetric mesh. A simply supported thin rhombic plate with a uniform load is modeled using a uniform symmetric spline discretization (Fig. 18). The integration for this problem is exact. A combination of Gauss product rules on a rectangle and a right triangle is used. The following parameters are used: L=loo
t = 0.01
DOF = 1275
v = 0.3
i
\
I’ / I
Q Fig. 18. Uniform
Q
Q
\ \ \ \
0
0
spline mesh for a rhombic
plate.
B.M. Donning, U!K. Liu/Comput.
61
Methods Appl. Mech. Engrg. 152 (1998) 47-71
The angle of the obtuse vertex is 150 degrees. Fig. 19 shows the moments near the singular obtuse vertex of the rhombic p1at.e. The solution is compared to a mesh of 384 Hughes-Tezduyar 4-node Tl elements Ill]. The computed solution agrees well qualitatively with the analytical solution [17] and captures the character of the singularity better than the FEM results. The inaccuracy near the boundary is attributed to the presence of a stress singularity there and the inaccuracy of the satisfaction of the essential BC’s there due to MLM (Fig. 20).
mom11 camp. mom22 camp. mom11 FEM mom22 FEM mom11 analytical mom2 2 analytical
-0.38
-
_-__. Q + --~--.-.--.--.
’
I
I
I
I
I
I
0
0.05
0.1
0.15 x/L
0.2
0.25
0.3
Fig. 19. Moments along the x-axis near the singular obtuse vertex of a rhombic plate.
o.c’02 I O.C’Ol
0
-O.C~Ol
-0.002
-0.c103
_O.C'04
0
I I I I I I I 50 100 150 200 250 300 350 length along boundary clockwise from the left obtuse vertex
Fig. 20. SS BC error along the rhombic plate boundary.
I 400
68
B.M. Donning,
WK. Liu/Comput.
Methods
Appl.
Mech. Engrg. 152 (1998) 47-71
5. Conclusion
The method proposed completely eliminates shear and membrane locking in beams and plates using an unmodified displacement-based variational principle. The efficiency and accuracy of the method rivals that of finite elements; however, mesh creation and post-processing to smooth stresses are eliminated. Error estimation is greatly simplified in the framework of the uniform discretization and adaptivity is accommodated through the insertion of nodes over the background uniform mesh. An extension to general shells (including composites) is being studied. In the framework of meshless methods, the satisfaction of essential boundary conditions remains a substantial problem and the effect on boundary stress results are apparent. Acknowledgements
We would like to thank the Office of Naval Research for partial support of this work. This paper is based on the Master of Science thesis of Brian Donning at Northwestern University. We would also like to thank Tom Black for his insight into numerical methods. Appendix
A. 1D uniform cardinal spline interpolation
A.l. Derivation of splines
Splines can be derived by solving a system of equations representing systematically through a convolution of lower-order ‘splines’. Let S’ be a spline of order p with
so=
1
-;
0
elsewhere
= (S” * S”)(x) = S n+m+l
IN S”Yy)S”(x --oo
their desired properties
or more
- y) dy
The property of the convolution (denoted by *) in Fourier space ( F(f * g) = FCf)F(g)) order splines to be computed easily using computer algebra software.
allows higher-
A.2. Continuity, consistency and general spline properties
An S’ spline has CP-’ continuity and pth order consistency. The ability of the interpolation field to exactly interpolate certain order polynomials will be termed ‘consistency’. This definition means that the interpolating functions, with the appropriate choice of coefficients (the coefficients themselves do not need to interpolate that polynomial), exactly interpolate every polynomial up to its order of consistency. The consistency and continuity of a spline are interrelated-the consistency would suggest higher continuity, but this is only true when the splines are interpolating an exact polynomial of that order. The consistency is based on uniform spacing on an infinite domain-a uniform mesh aligns the regions of piecewise functions which make up the spline. The basic spline can be scaled with respect to the mesh to create a more general spline. Spline ;‘: S2 (quadratic) S3 (cubic) . ..
Consistency 1 1,x 1,x,x2 1, x, x2, x3 .. .
Continuity c-’ co C’ C2 .. .
r
13M. Donning,
WK. Liu/Conzput.
Methods
Appl. Mech. Engrg. 152 (1998) 47-71
69
A.3. Explicit form of selected cardinal splines
Quadratic Spline, defined on
P(x)=
(symmetric about x = 0)
03 -. X2
X<--
;+;,+;,2
_2,,,_; 2
4
-;
2 3
Cubic Spline, defined on [-2,0] (symmetric about x = 0) x < -2
10 P(x)
=
4
+2x+x2+ -x3 1 3 6 2 1 - -.x2 - -X3 3 2
Quartic Spline, defined on 0, s [
-2
1(symmetric
about x = 0)
f 0
625 S4(x)
=
384:
X3, 5
__x3+_x~_-x+-x4
12
-+-x--x~+_X3__X4 55 5
5
25 16
96 24 4 115 5 1, 192;-8x2+4x \
125 5
6
48
5
1 1
24
6
Fig. 21 shows the splines used in the numerical examples. They are based on unit spacing. These splines are scaled with respect to the mesh spacing A_xby replacing x with x/Ax in the above definitions.
Fig. 21. Cardinal splines with unit spacing, p = 2,3,4.
70
B.M. Donning,
WK.
LidComput.
Methods
Appl.
Mech. Engrg. 152 (1998) 47-71
Appendix B. Boundary splines The spline interpolation can be modified on the boundary to eliminate splines entering from ‘outside’ the domain. Such functions are expected to improve conditioning since most of the function lies in the domain-there are no cases where a small part of the function contributes to the domain and creates numerically small stiffness entries. There is a unique solution for such a modification for each spline member which retains the following properties: exact derivative approximation of the next higher-order modified spline, the same polynomial character, the same continuity. However, the consistency drops in the boundary region due to the reduced degrees of freedom present there. The modified quadratic and cubic splines also have the interesting property of a FEM-like boundary treatment for rectangular domains and spline-based FE: only the boundary spline is nonzero on the boundary. However, due to the lower approximation order (consistency) on the boundary, the character of the derivatives loses generality. More specifically, for the modified quadratic spline, the derivative at the boundary is always zero-this results in very poor stress approximations near boundaries in structural formulations, areas where stresses are often most complicated. The character of these boundary splines also gives insight into other meshless methods such as RKPM which attempt to modify interpolation theories for infinite domains for use on finite domains. In lD, RKPM is identical to the above modified spline for uniform spacing using the cubic spline. Fig. 3 shows both the spline and RKPM interpolations. The modified quadratic spline (defined with the ‘center’ of support xj = 0, modified for the left boundary) is given by (Fig. 22) 7 1 1, S&(x) = jj - 2x - 5”
;<*
The modified cubic splines (defined with the ‘center’ of support xi = 0, modified for the left boundary) are given by (Fig. 23)
The modified quartic splines (defined with the ‘center’ of support xi = 0, modified for the left boundary) are given by (Fig. 24)
1 0.9 -
.
.\
’ri;
\ \ \ \ \
0.6 -
\
\
i”;. \
0.7
0.7-
Fi
\
:
1
0.6 0.5 0.4 -
.
.
.
.
.
.
-
0.3 0.2 0.1
-!.5
-1
-0.5
0
0.5
1
1.5
Fig. 22. Modified Fig. 23. Modified
2
quadratic
-2
-1.5
spline boundary
cubic spline boundary
-1
function. functions.
-0.5
0
0.5
1
1.5
2
B. M. Donning,
WK. Liu/Comput.
Fig. 24. Modified
463 71 7 s~,)(x)=~---_+--_2+-x3--x4 384 48 :: s~2~(x)=~+~x-_162+~x3+~x4 st,(,,=~-~.-~.2-a~3-~X4
Methods Appl. Mech. Engrg. 152 (1998) 47-71
quartic
1
1
12
24
spline boundary
71
functions.
-;
These modified boundary splines are defined using a unit mesh spacing, and they are scaled with respect to a non-unit mesh spacing Ax by replacing x with x/Ax. References [1] [2] [3] [4]
D. Kahaner, C. Moler and S. Nash, Numerical Methods and Software (Prentice-Hall, Englewood Cliffs, NJ, 1977). C. de Boor, A Practical Guide to Splines, (Springer-Verlag, New York, 1978). F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods (Springer-Verlag, New York, 1991). T.H.H. Pian and K. Sumihara, Rational approach for assumed stress elements, Int. J. Numer. Methods Engrg. 20 (1984) 168551695. [S] J.C. Simo and T.J.F!. Hughes, On the variational foundations of assumed strain methods, J. Appl. Mech. ASME 53 (1986) 51-54. [6] T. Belytschko, Y.Y. Lu and L. Gu, Element free Galerkin methods, Int. J. Numer. Methods Engrg. 37 (1994) 229-256. [7] J.C. Simo and M.S. Rifai, A class of mixed assumed strain methods and the method of incompatible modes, Int. J. Numer. Methods Engrg. 29 (1990) 1595-1638. [S] C.A. Duarte and J.T. Oden, An h-p adaptive method using clouds, Comput. Methods Appl. Mech. Engrg. 139 (1996) 237-262. [9] W.K. Liu and Y. Chen, Wavelet and multiple scale reproducing kernel methods, Int. J. Numer. Methods Fluids 21 (1995) 901-932. [lo] W.K. Liu, S. Jun and Y.F. Zhang, Reproducing kernel particle methods, Int. J. Numer. Methods Fluids 20 (1995) 1081-1106. [ll] T. Hughes, The Finlte Element Method (Prentice-Hall, Englewood Cliffs, NJ, 1987). [12] K. Washizu, Variational Methods in Elasticity and Plasticity, 1st ed., Chs. 4 & 9 (Pergamon Press, 1968). [13] S. Timoshenko and J. Goodier, Theory of Elasticity, 3rd ed. (McGraw-Hill, 1970). [14] R. Cook, D. Malkus and M. Plesha, Concepts and Applications of Finite Element Analysis, 3rd ed. (John Wiley & Sons, 1989). [15] Y.Y. Lu, T. Belytschko and L. Gu, A new implementation of the EFG method, Comput. Methods Appl. Mech. Engrg. 113 (1994) 397414. [16] S. Timoshenko and S. Woinowsky-Krieger, Theory of Plates and Shells, 2nd ed. (McGraw-Hill, 1959). [17] L. Morley, Skew Plates and Structures (MacMillan, NY, 1963) 89-99. [18] H. Stolarski and T. Belytschko, Membrane locking and reduced integration for curved elements, J. Appl. Mech. 49 (1982) 172-177.
ELSEVIER
Meshless methods
Methods
Appl.
Mech. Engrg.
152 (1998) 47-71
for shear-deformable
beams and plates
Brian M. Donning ‘, Wing Kam Liu *j2 Dept. of Mechanical Engineering,
Northwestern University, 2145 Sheridan Rd., Evanston, IL 60208, USA
Presented at Professor J. 7: Odenk 60th Birthday Symposium Received 13 February 1997
Abstract A meshless method IS developed to analyze moderately thick and thin structures using Mindlin-Reissner theory. A uniform discretization is used to allow for efficient integration and for the shape functions to be written explicitly. Irregular boundaries are modeled in a straightforward manner. An unmodified displacement-based Galerkin method is used; full integration is used to evaluate all energy terms, and convergence is independent of the thickness. Shear and membrane locking are completely eliminated pointwise at the interpolant level using cardinal splines. The continuity of the splines chosen results in continuous stresses. An extension to general meshless methods is given. Beam and plate examples show the accuracy of this method for coarse discretizations.
Introduction
To make nmneri~cal analysis of moderately thick/thin structures which exhibit shear deformation more practical, a method. needs to have the following properties: complete absence of numerical locking, accuracy for coarse discretizations, continuous stresses, automatic discretization, accurate geometry modeling, and numerical efficiency. Successful approaches to locking in finite element methods (FEM) have focused on modified variational principles. Assumed strain methods [5,7], hybrid stress methods [4], and reduced and selective integration methods [11,14] have all been used to ameliorate locking and are related to a general class of mixed methods [3]. However, these methods often lack robustness and accuracy on coarse meshes. Meshless methods inherently satisfy many of the problems associated with FEM. Different forms of meshless methods have emerged which differ in the construction of the interpolant, but all maintain several common properties. Reproducing kernel particle methods (RKPM) modify interpolants to satisfy consistency requirements via a correction function [9,10], element free Galerkin methods use a moving least squares procedure to construct interpolants [6], and H-p Clouds use a partition of unity and a set of basis functions to construct interpolants [8]. In each method, interpolation fields yielding continuous stresses are easily constructed, a FEM mesh need not be created to discretize the problem, and the geometry is arbitrary since there are no element edges, but only a set of points. However, efficiency and locking have limited the use of meshless methods to model thin structures.
*Corresponding author. ‘Graduate student in Theoretical and Applied Mechanics, Northwestern University. *Professor of Mechanical and Civil Engineering, Northwestern University. Email: [email protected] 0045-7825/98/$19.00 c:, 1998 Elsevier PZZ SOO45-7825(97)00181-3
Science
S.A. All rights reserved
48
B.M. Donning,
WK.
LidComput.
Methods
Appl.
Mech. Engrg. 152 (1998) 47-71
A new method is proposed to maintain the desirable properties of meshless methods, but also eliminate locking and improve efficiency to the level of FEM. The displacement-based Galerkin method for the principle of virtual work is used as the variational framework for this method (linear, static, elastic). Cardinal splines are used as the shape functions, and their interpolation properties greatly enhance the accuracy of the method for coarse meshes. The shape functions are given analytically, and integration is exact in the interior of the problem domains; background integration cells do not need to be created from a FEM mesh. Various examples in linear elasticity, including straight beams, curved beams, and plates of various geometries, show the effectiveness of this method. Also, in the context of general meshless methods, locking free interpolation fields are identified for beams and plates. RKPM is used to demonstrate this. An extension to FEM is also ma’de, using the interpolant for the proposed method as the local interpolant-the resulting elements are free of locking and have continuous stresses.
1. Mindlin-Reissner
beams via mesbless
methods
1.1. Beam model
The beam is isotropic, linear, elastic; the kinematic assumptions and strain-displacement the bending of straight beams are as follows: axial displacement transverse displacement
U(X,Z) = -ze w(x,z) = w Ex
=
u,x =
relations for
-ze,x axial strain shear strain
YZX= w,, -0
The Galerkin form for a rectangular straight beam is, neglecting body forces: SWex’ = gwbe”d
awbend
Ebt3
= -
12 6W shear = Gbkt
+
swshear
@,x dx J8,x -t3)(6w+ -60)d.x J w,, L
0 0
where t = thickness, b = width, k = shear correction factor, E = Young’s modulus, G = shear modulus and L = length of beam. 1.2. Shear locking Shear locking develops in formulations which cannot represent a state of zero shear in the thin beam limit independent of the mesh. This is exhibited as a violation of the Kirchhoff assumption in thin beam theory, stated as w,, - 8 = 0. When a beam is thin, most of the energy of deformation is due to bending. If a formulation does not admit non-trivial zero shear solutions, then most of the energy of deformation goes into shear deformation, resulting in very little bending-hence locking. This can be seen from the variational form; since the coefficient of the bending term is of order t3 and the coefficient for the shear term is of order t (fix L), the latter becomes much larger as t -+ 0. Since the variational form gives minimum energy solutions, the large shear term coefficient forces the minimization of the shear deformation. If the interpolant cannot physically represent such a state, spurious shear develops, giving a grossly inaccurate representation of the deformation state of the beam. To rid the formulation of locking, the interpolation of 8 must be able to approximate the derivative of the lateral displacement w arbitrarily well. The degree of satisfaction of this requirement determines the relative amount of spurious shear energy present and thus the degree of shear locking. Since meshless interpolants can be created from a wide variety of functions, many different combinations potentially
13.M. Donning,
WK. Liu/Comput.
Methods Appl. Mech. Engrg. 152 (1998) 47-71
49
could be used to alleviate locking. However, once an interpolation field that alleviates locking is found, it may not exactly satisfy the Kirchhoff constraint pointwise. This deficiency, while possibly giving accurate displacement results, causes oscillations in the stress fields to develop. The next four secl:ions outline the construction of locking-free interpolation fields for straight MindlinReissner beams. 1.3 Approximation
of spline derivatives
Spline functions are chosen here because of their many interpolation properties [2] (see Appendix A). A spline of order p is denoted by S’. The spline family has an additional property that S$ can be exactly represented by SP-’ for uniform spacing when the splines are spaced such that their piecewise components are aligned. S:(x) = Sp-‘(x + m) - Y’(x
- m)
For convenience, S::(X) is centered at x = 0, and 2m is the spacing between consecutive splines for both fields. Fig. 1 shows the alignment of quadratic and cubic splines fields, and the ability of the quadratic spline to exactly represent the derivative of the cubic spline. The vertical lines denote the division of the piecewise parts of both spline fields. Of many functions which could alleviate locking by approximating the Kirchhoff constraint, splines are very practical since they need relatively low integration orders and can easily give continuous stress fields. Fig. 2 shows the ability of the present formulation to model thin beams. A cantilevered beam with an end point load is modeled using 7 degrees of freedom. Conditioning problems start to degrade the solution at L/t = 106. Whereas traditional FEM solutions to locking in a displacement-based formulation have relied on reduced integration or assumed strain techniques, the present technique uses full integration, and all displacements and strains are related by a derivative at every point in the problem domain. FEM could use a similar technique, but the Co continuity inherent to standard FEM interpolants would cause the continuity for 8 to drop to C-‘. This would create Dirac delta functions in the stress interpolation. In addition, in the interpolation schemes given here, all field variables are at least C’, meaning stresses will be continuous.
0.6
1.6
\ \ \ ! \ \
-0.2-
-0.4--
2 Quadratic Splines
,
-0.6- - -Quad. SqlineDiWlerem
I.'
,
. .
..I
.
. . .
.
. . .
.
. . .
.
..C
P gi.1-
I'
\
,.-..I
x +1.3E e Sl.2‘
,'
.ji I
'i \ \ \
- - 1 CubicSpline
. ...I
/ .I
.P l-
-J
.- Deriv.of CubicSpline -0.B' -2
-1.5
-1
-0.5
0 x
0.5
1
1.5
Fig. 1. Spline alignment Fig. 2. Absence
os-
2
IO0
and derivative
of locking
...._I 10'
".
c ___ IO'
approximation.
in the thin beam
limit.
IO3
u
.,_c ..._.f 10' loa
. ..._a 10'
.....10'
50
B.M. Donning,
WK. Liu/Comput.
Methods
Appl. Mech. Engrg. 152 (1998) 47-71
1.4. Interpolation via cardinal splines
To maintain consistency on a finite domain, the finite domain can be ‘cut out’ of an infinite uniform interpolation field (Fig. 3). Only the parts of the interpolants which overlap with the problem domain are used. The following method is for uniform spacing, yet it still satisfies the interior consistency at the boundary for irregular geometries by extending the mesh outside the physical domain. W(X) =
C
8 (X)
C Sp-’(X -
=
S”(X - xi)di
Xi)hi
Here, di and hi are the unknowns, and Xi is the position of the center of each spline. A cardinal cubic spline for the w interpolant and a cardinal quadratic spline for the 8 interpolant were used in all the beam examples. Three point Gauss quadrature in each nodal interval is needed for exact integration (based on the w nodal spacing). If higher continuity is desired, higher-order splines can be used in this same setting, as long as the 8 field is one order lower than that of w; higher-order quadrature and a larger bandwidth in the system equations would be the tradeoff for the higher continuity in the spline family. 1.5. Interpolation for any C’ meshless method
Since these meshless methods have high continuity, the interpolation field for the rotation of the midline is chosen by differentiation of the field for the lateral displacement. In meshless methods suitable for irregular discretizations, the finite domain includes only those shape functions whose centers are included in the domain; subsequently, the interpolants must be modified to satisfy consistency near the boundary and in the interior for irregular spacing. Using a cubic spline as the interpolant and a uniform discretization, RKPM only modifies the splines on the boundary to maintain consistency (Fig. 3). The approximation level of the modified shape functions near the boundary is decreased due to the elimination of shape functions whose centers are outside the domain, reducing the consistency for the global interpolant. The interior shape functions (chosen to be the cubic spline) have third-order approximation character, whereas on the boundary, the modified splines only have first-order
lnltnite 0.8
spline interpo!mion 1
;i2; -1.5
-1
-0.5
0
0.5
1
1.5
2
Modifiedboundaryinterpoiation
Fig. 3. Two different
boundary
interpolations
for a cubic spline S3.
E1.M.Donning,WK. Liu/Comput. MethodsAppl. Mech. Engrg. 152 (1998) 47-71
character. Also, the modified functions in general will no longer be polynomials for non-uniform so integration of these modified functions must be done carefully. W(X) =
C
51
spacing,
Nil:X)di
i
8
(X) =
C
Ni,x(x)hi
i
Here, N, is an RKI’M or EFG, etc. shape function. The interpolant Ni,X is continuous for C’ or higher fields for W(X). However, this shape function is no longer non-negative, and cancellation can cause numerical problems [l]. 1.6. Interpolation via spline based FEM Only a suitable mapping field must be chosen to implement a locking-free spline based element. Here, the mapping is assumed to be sufficiently continuous. x = f(s)
($)-lhi
f3(X)=CN,(t)
i
Here, the properties
of cardinal splines are used:
N&L!) = SP(S - &)
These are locking-free
interpolants,
since they exactly satisfy the Kirchhoff constraint.
1.7. Bending of curved Mindlin-Reissner
beams
Curved beams exhibit both shear and membrane locking. The following method follows the same ideas as the straight beam. The beam considered is a circular quarter arc (rectangular cross-section), loaded in its plane of curvature (see Fig. 4). The midline of the beam is parameterized with respect to the global Cartesian coordinates xi. The curvilinear coordinates c& are orthogonal. x1 =x=Rcos-Sx2 =y
=Rs_inL
a1 = t
R
a2 = s
R
a3 = z
x3 = z
A point on the midline is given as
The outward normal to the midline is n(“)=cosSe1+sinSe2 R
R
B.M. Donning, WK. Liu/Comput.
Methods Appl. Mech. Engrg. 152 (1998) 47-71
-RR Fig. 4. Geometry of circular arc.
A point
the beam is given by
in
r(O) = r!)(s)
+ tn(O)(s)
where (s, t) are the coordinate
parameterization.
1.8. Continuum mechanics The strains for the curved beam are derived from displacement components defined in the curvilinear system which are analogous to the Mindlin-Reissner kinematic assumptions used in the straight beam
PI. The base vector along (Y*is
i3rp -=thf2
(&O) =
a,(y)
II-I
-sin Se, + cos Se2 R R
da2
The base vectors g, in the curvilinear coordinates g1
are
dr co) = -da’ = J&O) dr co)
t3-r’
dd”)
g,=g=da2+tda2=
t +1
( 1 u(oj
E
g3 = e3.
The volume transformation
between Cartesian and curvilinear is
dV = &dsdtdz where ,,& = g, ’ (g3 x gI) = (t/R + 1). The position of a point after deformation r =
r(O)
+ u
where (u, 0 and w are the displacement u =
(u+te)a(“)+wn(o)
=
variables) (u+t~)g*+w&+a3.
is
B.M. Donning,
WK. Liu/Comput.
Metho&
Appl.
Mech. Engrg. 152 (1998) 47-71
53
The displacement ;Z is not defined because it is left to satisfy the plane stress conditions. The metric tensor gii = gi . gj is
t
1
0
[gijl= 0
0 2
() +1
0
0
1
3
0
1
I
The strains are derived as follows. Let u = uigi and the deformed Gi . Gj. Then define the strain as
base vectors Gi = r,i with Gij =
eij = k (G,j - gij) = i (ui(j + Ujii + UkliUkIj) and linearize in the displacements Eij
=
k
(gikUk
j
+
gkjL&Ji)
where ui lj = u:j + { Gj}r_.8. The Christoffel symbol of the second kind is denoted by {ii}. The following two expressions for the linearized strain in curvilinear coordinates do not follow the summation convention.
2Eij = yij = g,i g
+ gjj g
i # j
The nonzero in-plane strain components Es s! E22 =
( I( J jr
+l
are
u,, +to,, + w
R>
‘y,$E 2E,2 = tv,, + 8 - ;. The following substitution is made to separate the rigid body rotation (u/R) from the total rotation (0) of the midline normal:
The strains are now rewritten in terms of U, w and 4. Yts = w,s + 4
The membrane Em=
strain is
(;+1) (u,,+i)
The bending strain. is
For the constitutive relations, a local Cartesian system coincident with the orthogonal curvilinear system ci is used, since the plane stress assumptions through the thickness of the beam are aligned
54
B. M. Donning, WK. Liu/Comput.
Methods Appl. Mech. Engrg. 152 (1998) 47-71
with the local Cartesian system. The material is linear and isotropic. The stress, strain and displacement variables in this local system are denoted with an overbar (no sum). ii, = 6
ui 1
Since ~~j~Ejj = &jSeii, the constitutive a22 =
relation will only be developed
in the local Cartesian system:
E~22
EEZ2= -
g22
1.9. Variational form Neglecting body forces and using the curvilinear coordinate
J Js( Jl( &j
SEij dV
zz
system:
swext
V
EQ
8~
+ 2Gkr12 8~
k22j2
E~22
8~22
+ 2Gke2, Se2, g22
g22 + Gb2
h12
g22
k2212
JJ(E(&)’
&stGk(-&j
&
&ds
dt = Z5Wex’
)
ds dt = 8WeX’
1
y&rs)
dsdt=sW”“’
where the width of the beam is 1 unit. Using the interpolations for u, 8 and w, and integrating analytically in t, the final matrix form is found in terms of the arc-length parameterization s.
1.10. Shear and membrane
locking
Shear deformation is important to model moderately-thick to thick beams. Most of the strain energy in a thick beam goes into shear deformation, whereas in a thin beam, most energy goes into bending deformation. Shear locking appears in thin beams if the interpolation does not allow bending without the generation of spurious shear strains. Membrane action comes from the coupling of bending deformation and membrane deformation in curved beams. Membrane locking [18] appears when the coupling is too strong-the interpolation does not allow for a state of bending without corresponding membrane deformation. This is not necessarily non-physical, since curved beams are supposed to show this coupling; however, since they are dependent in a specific way determined by the interpolants, there is no flexibility to allow different amounts of coupling which appear in different beam configurations (esp. thin beams). Let the beam become very thin so membrane and shear locking become issues: t/R -+ 0 : fix R, t + 0. Then the variational form for the virtual work becomes
B.M. Donning,
WK. LidComput.
Methods Appl. Mech. Engrg. 152 (1998) 47-71
55
shear membrane
+Et’($+&,)
dsdt
(%+a+,,
. bending
The coefficients of the terms in the minimum energy functional act as penalties [14] on the membrane and shear energy terms: E>Gk>>Et=
The appearance of locking comes from the inability of the interpolation to represent certain states of strain. When a curved beam becomes thin (t/R + 0) , almost all of the deformation comes from bending, and two constraints appear. (1) ‘yIs= w,, + $ = 0
shear constraint
(2) 6, = u,s .t w = 0
membrane
R
constraint
Traditional finite element interpolation fields only satisfy these constraints on coarse meshes with a zero solution for all fields, hence the term locking. Using higher continuity fields, both membrane and shear locking can be avoided completely (using exact integration) for all meshes with an appropriate choice of interpolants. 1.11. Interpolation via cardinal splines
The two constraints from the previous section define the relationship between the interpolants for u, w and 4. For comtraint (1) to hold, 4 must be able to approximate the derivative of w. For constraint (2) to hold, w must be able to approximate the derivative of u. u=
cS”(s -
w =
c
#I =
c sqs - Sj)rni
sj)d,
sp-’ (s - s;)h;
The properties of the spline family on a uniform discretization are exploited for this choice of interpolants. Since one member of the spline family can exactly interpolate the derivative of the next higher spline if the piece-wise components of the splines are aligned, one only needs to choose three ‘consecutive’ members of the spline family for u, w and 4 in descending order. The choice of the order for u automatically defines the order of w and 4. A quartic spline is the lowest-order spline for L(which results in continuous stresses (p 3 4). 1.12. Interpolation for any C2 meshless method
In the context of a general meshless method, a C2 interpolation is needed for a set of nodes for u (a C3 or higher field is needed for all stresses to be continuous). Then the shape functions for w and 4 are derived simply through differentiation of the shape functions for u.
56
B.M. Donning,
u = C
WK. Liu/Comput.
Methods
Appl. Mech. Engrg. 152 (1998) 47-71
Ni(S)d;
1.13. Interpolation via spline based FEM Similar to the straight beam, a spline based locking-free element can be constructed. s = f(t) ~(3) = C Sp(t - &)di
W(S)
C Sp-‘(5 - &)
=
($) -’hi
i
These are locking-free interpolants, since they exactly satisfy the shear and membrane constraints. Because the coordinate mapping creates an extra term in the derivative, twice as many degrees of freedom are needed to interpolate 4; this can be avoided by using the true derivative as the shape function, but conditioning problems may appear.
2. Mindlin-Reissner
plates via meshless
The classical Mindlin-Reissner setting and removed.
methods
plate bending model is used. Shear locking is identified in this 2D
2.1. Plate formulation A static, linear, isotropic plate is formulated as follows [ll]. The 3D plate domain V is reduced to a 2D domain by defining all variables with respect to the plate’s midsurface A. V =
{
(x,y,z)
The displacements u1 =
-Z&(&Y)
u2
=
-ZW,Y)
u3
=
W(X,Y)
E R31z E
-;, 1
;
are defined as
1
, (x,y) E A c R2
For an isotropic material (sum on k) uii = 2peij + A6ijekk (i, j, k = 1,2,3)
>
B.M. Donning,
WK. Liu/Comput.
Methods
Appl.
Mech. Engrg. 152 (1998) 47-71
57
The plane stress condition imposes CT33=
0 +E33=-hf2P
h(
(1)
El’ + E22)
The stresses are then simplified by eliminating ~~~~The strain-displacement Eij
=
k
(Ui,j
--
relations are
Uj,i)
This implies l33 = w,~ = 0, but this contradicts the plane stress assumption (1). However, this discrepancy does not enter into the variational form since the virtual work done by the true deformation across the thickness is zero (8~~~r33 = 0). The variational form is (a, p = 1,2):
J
(SE,~ aUp + 2k Sea3 uG3) dV =
V
SwT,dS+C6wF, ss
i
Here, Fi are the point loads normal to the midsurface, T, are the applied tractions midsurface and k is the shear correction factor. Body forces are neglected.
normal to the
2.2. Shear locking In this plate formulation,
two shear locking constraints appear as t + 0 (c~ = 1,2):
w>, -l&=0 Since these constraints are each associated with an independent direction, they can be treated separately in succession. Once the interpolation for w is chosen, applying each constraint in turn generates 13’and 02.
2.3. Interpolation via products of cardinal splines Splines for highker dimensional approximation are derived via a product rule from 1D splines. As a consequence, the higher dimensional splines behave exactly like the 1D splines in each direction. In 2D: P(x,y)
= SP(x)S4(y)
The resulting functions have square support and consist of piecewise polynomial patches. The consistency properties and ability to approximate derivatives of higher-order splines remains intact in each direction. Combinations of different order splines in each direction can approximate mixed and other derivatives of higher-order spline products.
sjJ$, y) = SP-‘(X + m)S4(y) - Y’(X
- m)Sq(y)
Sp::(x,y) = SP-‘( X + m)Y-’ 0, + m) - s-’
(X + m)W
- 9-l (X - m)S@ (y + m) + 9-l (x - m)W
(y - m) 0, - m)
These propertie:s are all based on rectangular uniform collections of spline functions with nodal spacing m in each direction. The ability of such a collection to serve as an interpolation field for irregular geometries relies on the arbitrariness of boundaries in meshless methods. Meshless methods usually choose some polygonalization of the boundary nodes to approximate the problem geometry. This is purely arbitrary, since the interpolant is not defined by the polygonalization, but merely truncated by it. For the uniform splines, the exact geometry can be used (likewise in any meshless method) in an analogous manner. Qualitatively, a uniform spline interpolation field is created which extends outside the boundary; the exact geometry is superimposed on the interpolation field and the interpolant which intersects with the domain is used in the interpolation field (Fig. 5). This maintains all the properties of the splines on a finite domain without modification near the boundary. This results in the inclusion of some functions whose nodes (center of support) are exterior to the problem domain. However, only
B.M. Donning,
WK. Liu/Comput.
Methods Appl. Mech. Engrg. 152 (1998) 47-71
2D Problem Boundary
l
0
0
Fig. 5. 2D uniform
discretization
Interior nodes
Non-zero shape function region for center node
of irregular
domains.
the nonzero region extending into the domain is used in the Galerkin procedure. The accuracy of the boundary is determined by the accuracy of the integration scheme around the boundary. Interpolation via cardinal splines: 8, =
c
02 =
C Sq-’(J
W
C Sp(X - Xi)S’O, - yi)mi
=
sp-’ (x - x;)Y
0, - Yi)hi
- yj)S’(X - Xi)gi
Integration of these uniformly spaced splines is exact in the interior for an appropriate Gauss quadrature rule. Since the interpolation field consists of piecewise regions of polynomials, each such region can define a cell for the integration. Product rules of 1D Gauss quadrature formulas integrate these functions exactly in the interior. Near the boundary, the problem of integrating a general multidimensional domain appears. The interpolant is still piecewise polynomial in character, but the area or volume of the cell is no longer rectangular. Subdivision or polygonalization of boundary cells reduces quadrature error, and exact Gaussian product rules for right triangles are easily constructed, but many more quadrature points are required for the right triangle than for a rectangle. 2.4. Interpolation for any C’ meshless method Interpolation
6
=
via any sufficiently smooth method:
CNi,y(X,Y)gi i
B.M. Donning,
2.5. interpolation
WK. LidComput.
Methods
Appl. Mech. Engrg. 152 (1998) 47-71
59
via spline based FEM
A sufficiently smooth coordinate x = f(5777)
mapping in 2D is chosen:
Y = g((, 77)
J=axay__axdy 86 877 877 a,$
e,(x) =
3 :~Sq(7] [
77i)Sp-‘(5 - ~~)~aY
(1)
hi
-
C
Sp(~ -
E)Sq-‘(rl -
77i) ~
h12)
i
i
Because the coordinate mapping creates an extra term in the derivatives, twice as many degrees of freedom are needed to interpolate 61 and 02; this can be avoided by using the true derivative as the shape function.
3. Essential bounldary conditions The displacement boundary conditions were enforced via Lagrange multipliers and a modified Lagrange multiplier [15] method in the examples given. The summation convention is used in this section. 3.1. Collocation
Collocation was tested in the two straight beam problems considered in the numerical examples section. The performance of collocation on coarse meshes was very poor, since a collocation equation replaces an equation from the variational form. The lost information has a profound effect on the solution near the boundary. The two following methods were implemented to impose boundary conditions while retaining all of the stiffness equations. 3.2. Lagrange snultiplier method
This approach adds constraint equations to the system of equations, coupling degrees of freedom via Lagrange multipl:iers. A term is added to the internal work of the variational equation as follows: @int =
wint
_
hi (ui - Ui) dS .%‘ The Lagrange multiplier fields hi are arbitrary fields on the plate boundary where the displacement boundary conditilons Ui are to be imposed. This approach has a few disadvantages. The constraints augment the size of the system matrix; the form of the Lagrange multipliers renders the system of equations indefinite, requiring the use of a special solver for large systems;. the conditioning of the equations (in the sense of the LINPACK condition estimator) is large for every value of t. Also, the choice of functions for the hi fields, while still providing good satisfaction of the boundary conditions, affects the interior solution markedly in coarse meshes. However, the arbitrariness of the Lagrange multiplier fields potentially allows for more accurate satisfaction of BC’s on non-smooth boundaries. I
B.M. Donning,
60
WK. Liu/Comput.
Methods
Appl. Me&.
Engrg. I52 (1998) 47-71
3.3. Modified Lagrange multiplier method
The modified Lagrange multiplier method [15] solves some of the problems associated with the classical Lagrange multiplier approach. Here, h; is replaced with its physical counterpart (the equivalent reaction tractions of the properly constrained system), called T,R, arising from the Euler-Lagrange equations of the variational principle. In this form, the constraints do not add any additional degrees of freedom and the arbitrariness of the choice of Ai is no longer an issue. Also, the positive-definiteness of the original system is retained [lS]. @in, = jyint
_
J
Tr(ui - U;) dS
s,,
Here, nj is the outward normal to the boundary and cij is related to the displacement degrees of freedom via the constitutive and strain-displacement relations. Since the procedure modifies stiffness terms of the boundary associated shape functions, the accuracy of this method depends on the relative accuracy of integration of both the stiffness terms of the boundary interpolants and the boundary integration for the essential BC term. The numerical behavior of this method is very good for coarse meshes for the plate problems solved. The conditioning of the system matrix depends on t as is expected for a pure displacement formulation. For Co geometries, the accuracy of the BC satisfaction of this method greatly deteriorates near the non-smooth points relative to the accuracy on the smooth regions; this inaccuracy shows up as stress oscillations.
4. Numerical results 4.1. Straight beam Two thin beam problems are solved: problem 1 is a cantilever beam with a point end load (Fig. 6), and problem 2 is a clamped-clamped beam with a center point load (Fig. 7). Lagrange multipliers are
=:[[I 0
2
1
3
4
2 X 1O-6 *r-
'---'
Fig. 6. Cantilever Fig. 7. Clamped/clamped
7 ii_
9
10 .
- - w,x w,xexacl
. . '\ '..
0'
beam with end point load (7 DOFs). beam with center
6
._
/' / E
5 x
point load (9 DOFs).
\
,
B.M. Donning, WK. LidComput.
Methods Appl. Mech. Engrg. 152 (1998) 47-71
used to enforce essential boundary conditions. Both problems use the following set of parameters, w is interpolated via cubic splines, and ~9is interpolated via quadratic splines: V = 0.3
L = lo4
E’ = lo5
61
and
t=b=l
The exact thin beam (Euler-Bernoulli)
solution for problem 1 is
w(x) = mQ (c-X3 + 3Lx2 \ > and for problem 2 H(x - L/2)
x--)I L
3
(
3
2
where Q is the magnitude of the point force and the Heaviside unit step function H(t) is defined as zero for t < 0, 1 for t 3 0. The exact solution for problem 1 is reproduced using the minimum number of degrees of freedom to discretize the beam since the approximation order of the interpolant is third order, and the the straindisplacement equations are satisfied everywhere. Problem 2 is piecewise cubic and is reproduced exactly by a piecewise cubic field. 4.2. Curved beam The 2D plane stress elasticity solution for the quarter circular arc with an end point load (Fig. 8 ) is taken from [13]. 20 u = -Ft?
D(1-u)lnr+A(1-3u)r2+B(~~V)
case +
+Ksin8+LcosB I
u=z/jsin(j - Cos[A(5+v)r2-D(l-v)lnr+ E
-L sin 8 -t D’y
E
”
B(l+ v) r2
1
+Kcos8-LsinO+...
cos 8
X
Fig. 8. Geometry
for exact solution.
62
B.M. Donning,
A=2N
WK. Liu/Comput.
P
Methods
Appl. Mech. Engrg. 152 (1998) 47-71
N = a2 - b2 + (a2 + b2) In i
D = -$(a2+b2) B(l + v)
r2
D(l-v)lnr+A(l-3v)r2+
I The solution in Fig. 9 is in the coordinates of Fig. 4. The solution uses the same material parameters as for the straight beam with R/t = 2/7r x 104. The solution breaks down for R/t 3 105, possibly because the conditioning becomes very poor in the thin limit. The accuracy for this problem is uniform in t. Full integration is used (4 pt. Gauss quadrature) and a quartic spline is used for u, a cubic spline is used for w and a quadratic spline is used for 4. A uniform spline mesh is used and essential BC’s are imposed via Lagrange multipliers. Here the solution is no longer a low-order polynomial and the splines used give good agreement with the exact solution for the prime displacement variables and their derivatives with a minimal number of nodes. 4.3. Plates In each example, the whole plate is modeled and, where possible, symmetric discretizations and quadrature meshes are used (for symmetric problems). A product of cubic splines is used to interpolate w and corresponding products of cubic and quadratic splines are used to interpolate 8i and 0,. The boundary conditions employed are as follows: Displacement CL SS
BC’s
Clamped Simply Supported
W = 0,8, = 0, 0, = 0 w= 0
Traction BC’s C Concentrated center load U Uniform lateral load
El El --._ F-J 4X1O-6
0.25
I
El
0.2
, /’
- wcmlp.
w exact
0.15
i 0.1
,r
,
I’
,
22
I’
0
*’
I’ 1.’
/’
0
I
5000
00
10000
4.1
-0.15
‘\
0
‘\ \
‘\
I -
0
-.
mmp.
w.0
oxacl
loo00
-1
‘\
-
\\
‘\
“3
‘\\
-2
\
El ‘W
-.
5
\
u wmp.
‘U
-0.2
W.8
,XI04
-0.05 =
-El
5ow s
s
0
‘1
/.’
,//
,’
0.05
- -.._
_I
/’
-3
axact
exed
amp. ‘\
‘\
‘\ 1oMx)
-4 : 0
IWO0
s
s
Fig. 9. Numerical
solution
for a curved
beam
(15 DOFs).
B.M. Donning,
WK. LitdComput.
Methods Appl. Mech. Engrg. 152 (1998) 47-71
63
For a thin square plate, only the uniform spline method was used, since RKPM is closely related for rectangular geometries. The convergence of the center displacement is shown in Fig. 11. A symmetric uniform discretization is used (Fig. 10). The following parameters were used in the problem: L = 10
t ==1 x 1o-4
E = 10’
v = 0.3
k = 1.19
To compare with finite elements, the convergence plots are reported in equivalent linear 4-node elements for a quarter of a uniform square mesh using the following relation:
(2) This relation relates the degrees of freedom (DOF) of a given discretization of splines to the number of elements for a quarter symmetric mesh of linear elements containing the same number of degrees of freedom. The uniform spline method is used to solve a circular plate (Fig. 14). Oscillations in the convergence plots are due to the inaccurate integration near the boundary (rectangular Gauss product rules were used, and the integration weights of points outside the domain were set to zero). A symmetric discretization (with respect to a quarter of a circle, Fig. 12) and integration mesh are used. The following parameters are used: R=5
t = 1 x lop4
v = 0.3
E = lo5
k = 1.19
To compare with finite elements, the convergence plots are reported in equivalent linear 4-node elements for a quarter of a uniform circular mesh using the following relation: equiv.elementsfor$mesh=
(DoL-l)
- (y)
(3)
The moments for a thin circular CL-U plate (676 DOF) are plotted with the exact solution [16] in Fig. 13. The thin circular plate is also solved using RKPM (Fig. 16). A non-uniform discretization is used (Fig. 15). A circular cubic spline window is employed as follows:
Thin Square Plate 0
I.&-
~_________-+_________+-___-____~ 8
Q
I
I
I
I
I
I
,
1.01 -
zB
_______ 1
r
D
go.99 b
6 o
0
I
,
I
I
0.96
= ; g 0.97 -
b_________~_--______~___-______6
Q 0.96 -
0
Q
Q
0.95 IO"
Q
0
Fig. 10. Uniform Fig. 11. Convergence
10’
Q
equiv. 4-node elements
spline mesh for a square
plot uniform
spline method-thin
plate. square
plate.
10'
B.M. Donning, WK. LidComput.
64
Methods Appl. Mech. Engrg. 152 (1998) 47-71
3
2
1 t i 0 5
0
5
-1
i -2
-3
-4 0
0.2
0.4
0.6
0.8
r/R
Fig. 13. Moments
along the x-axis for a thin circular
plate
(CL-U).
Thin Circular Plate 1.02
I
,
1.01
‘\
I’
0.94
‘\
’ 10’
102 equn. 4-m&
Fig.
14. Convergence
for uniform
elements
sphnedircular
plate.
1
KM. Donning,
WK. Liu/Comput.
Fig. 15. RKPM Fig. 16. Convergence
65
Methods Appl. Mech. Engrg. 152 (1998) 47-71
mesh for a circular of RKPM
plate.
for thin circular
plate.
o.ooc25
computed analytical
----.
5e-a5
0
0.2
0
0.4
0.6
0.8
1
r/R
Fig. 17. Thick circular
S3(X
-Xj,Y
-,Yj)
=
plate profile.
S3 (f)
where r =
d(X -
Xj)”
+
(Y -
Yj)’
The quadrature mesh is based on a uniform pattern of cells. A FEM background create the cells. The following parameters were used: R=5
t =
1 x 1o-2
v = 0.3
E = lo5
a =
1.5
mesh is not used to
66
B.M. Donning,
WK. Liu/Comput.
Methods Appl.
Mech. Engrg. 152 (1998) 47-71
This plate is thicker than for the uniform spline problem because conditioning problems were encountered at a lower L/t ratio. A clamped thick circular plate with a uniform load is modeled with a uniform mesh. The computed solution for the displacements are plotted against the exact solution [16] in Fig. 17. The following parameters were used: R=5
t=2
DOF = 1149
Using Eq. (3), this solution corresponds to 90 4-node linear elements on a quarter symmetric mesh. A simply supported thin rhombic plate with a uniform load is modeled using a uniform symmetric spline discretization (Fig. 18). The integration for this problem is exact. A combination of Gauss product rules on a rectangle and a right triangle is used. The following parameters are used: L=loo
t = 0.01
DOF = 1275
v = 0.3
i
\
I’ / I
Q Fig. 18. Uniform
Q
Q
\ \ \ \
0
0
spline mesh for a rhombic
plate.
B.M. Donning, U!K. Liu/Comput.
61
Methods Appl. Mech. Engrg. 152 (1998) 47-71
The angle of the obtuse vertex is 150 degrees. Fig. 19 shows the moments near the singular obtuse vertex of the rhombic p1at.e. The solution is compared to a mesh of 384 Hughes-Tezduyar 4-node Tl elements Ill]. The computed solution agrees well qualitatively with the analytical solution [17] and captures the character of the singularity better than the FEM results. The inaccuracy near the boundary is attributed to the presence of a stress singularity there and the inaccuracy of the satisfaction of the essential BC’s there due to MLM (Fig. 20).
mom11 camp. mom22 camp. mom11 FEM mom22 FEM mom11 analytical mom2 2 analytical
-0.38
-
_-__. Q + --~--.-.--.--.
’
I
I
I
I
I
I
0
0.05
0.1
0.15 x/L
0.2
0.25
0.3
Fig. 19. Moments along the x-axis near the singular obtuse vertex of a rhombic plate.
o.c’02 I O.C’Ol
0
-O.C~Ol
-0.002
-0.c103
_O.C'04
0
I I I I I I I 50 100 150 200 250 300 350 length along boundary clockwise from the left obtuse vertex
Fig. 20. SS BC error along the rhombic plate boundary.
I 400
68
B.M. Donning,
WK. Liu/Comput.
Methods
Appl.
Mech. Engrg. 152 (1998) 47-71
5. Conclusion
The method proposed completely eliminates shear and membrane locking in beams and plates using an unmodified displacement-based variational principle. The efficiency and accuracy of the method rivals that of finite elements; however, mesh creation and post-processing to smooth stresses are eliminated. Error estimation is greatly simplified in the framework of the uniform discretization and adaptivity is accommodated through the insertion of nodes over the background uniform mesh. An extension to general shells (including composites) is being studied. In the framework of meshless methods, the satisfaction of essential boundary conditions remains a substantial problem and the effect on boundary stress results are apparent. Acknowledgements
We would like to thank the Office of Naval Research for partial support of this work. This paper is based on the Master of Science thesis of Brian Donning at Northwestern University. We would also like to thank Tom Black for his insight into numerical methods. Appendix
A. 1D uniform cardinal spline interpolation
A.l. Derivation of splines
Splines can be derived by solving a system of equations representing systematically through a convolution of lower-order ‘splines’. Let S’ be a spline of order p with
so=
1
-;
0
elsewhere
= (S” * S”)(x) = S n+m+l
IN S”Yy)S”(x --oo
their desired properties
or more
- y) dy
The property of the convolution (denoted by *) in Fourier space ( F(f * g) = FCf)F(g)) order splines to be computed easily using computer algebra software.
allows higher-
A.2. Continuity, consistency and general spline properties
An S’ spline has CP-’ continuity and pth order consistency. The ability of the interpolation field to exactly interpolate certain order polynomials will be termed ‘consistency’. This definition means that the interpolating functions, with the appropriate choice of coefficients (the coefficients themselves do not need to interpolate that polynomial), exactly interpolate every polynomial up to its order of consistency. The consistency and continuity of a spline are interrelated-the consistency would suggest higher continuity, but this is only true when the splines are interpolating an exact polynomial of that order. The consistency is based on uniform spacing on an infinite domain-a uniform mesh aligns the regions of piecewise functions which make up the spline. The basic spline can be scaled with respect to the mesh to create a more general spline. Spline ;‘: S2 (quadratic) S3 (cubic) . ..
Consistency 1 1,x 1,x,x2 1, x, x2, x3 .. .
Continuity c-’ co C’ C2 .. .
r
13M. Donning,
WK. Liu/Conzput.
Methods
Appl. Mech. Engrg. 152 (1998) 47-71
69
A.3. Explicit form of selected cardinal splines
Quadratic Spline, defined on
P(x)=
(symmetric about x = 0)
03 -. X2
X<--
;+;,+;,2
_2,,,_; 2
4
-;
2 3
Cubic Spline, defined on [-2,0] (symmetric about x = 0) x < -2
10 P(x)
=
4
+2x+x2+ -x3 1 3 6 2 1 - -.x2 - -X3 3 2
Quartic Spline, defined on 0, s [
-2
1(symmetric
about x = 0)
f 0
625 S4(x)
=
384:
X3, 5
__x3+_x~_-x+-x4
12
-+-x--x~+_X3__X4 55 5
5
25 16
96 24 4 115 5 1, 192;-8x2+4x \
125 5
6
48
5
1 1
24
6
Fig. 21 shows the splines used in the numerical examples. They are based on unit spacing. These splines are scaled with respect to the mesh spacing A_xby replacing x with x/Ax in the above definitions.
Fig. 21. Cardinal splines with unit spacing, p = 2,3,4.
70
B.M. Donning,
WK.
LidComput.
Methods
Appl.
Mech. Engrg. 152 (1998) 47-71
Appendix B. Boundary splines The spline interpolation can be modified on the boundary to eliminate splines entering from ‘outside’ the domain. Such functions are expected to improve conditioning since most of the function lies in the domain-there are no cases where a small part of the function contributes to the domain and creates numerically small stiffness entries. There is a unique solution for such a modification for each spline member which retains the following properties: exact derivative approximation of the next higher-order modified spline, the same polynomial character, the same continuity. However, the consistency drops in the boundary region due to the reduced degrees of freedom present there. The modified quadratic and cubic splines also have the interesting property of a FEM-like boundary treatment for rectangular domains and spline-based FE: only the boundary spline is nonzero on the boundary. However, due to the lower approximation order (consistency) on the boundary, the character of the derivatives loses generality. More specifically, for the modified quadratic spline, the derivative at the boundary is always zero-this results in very poor stress approximations near boundaries in structural formulations, areas where stresses are often most complicated. The character of these boundary splines also gives insight into other meshless methods such as RKPM which attempt to modify interpolation theories for infinite domains for use on finite domains. In lD, RKPM is identical to the above modified spline for uniform spacing using the cubic spline. Fig. 3 shows both the spline and RKPM interpolations. The modified quadratic spline (defined with the ‘center’ of support xj = 0, modified for the left boundary) is given by (Fig. 22) 7 1 1, S&(x) = jj - 2x - 5”
;<*
The modified cubic splines (defined with the ‘center’ of support xi = 0, modified for the left boundary) are given by (Fig. 23)
The modified quartic splines (defined with the ‘center’ of support xi = 0, modified for the left boundary) are given by (Fig. 24)
1 0.9 -
.
.\
’ri;
\ \ \ \ \
0.6 -
\
\
i”;. \
0.7
0.7-
Fi
\
:
1
0.6 0.5 0.4 -
.
.
.
.
.
.
-
0.3 0.2 0.1
-!.5
-1
-0.5
0
0.5
1
1.5
Fig. 22. Modified Fig. 23. Modified
2
quadratic
-2
-1.5
spline boundary
cubic spline boundary
-1
function. functions.
-0.5
0
0.5
1
1.5
2
B. M. Donning,
WK. Liu/Comput.
Fig. 24. Modified
463 71 7 s~,)(x)=~---_+--_2+-x3--x4 384 48 :: s~2~(x)=~+~x-_162+~x3+~x4 st,(,,=~-~.-~.2-a~3-~X4
Methods Appl. Mech. Engrg. 152 (1998) 47-71
quartic
1
1
12
24
spline boundary
71
functions.
-;
These modified boundary splines are defined using a unit mesh spacing, and they are scaled with respect to a non-unit mesh spacing Ax by replacing x with x/Ax. References [1] [2] [3] [4]
D. Kahaner, C. Moler and S. Nash, Numerical Methods and Software (Prentice-Hall, Englewood Cliffs, NJ, 1977). C. de Boor, A Practical Guide to Splines, (Springer-Verlag, New York, 1978). F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods (Springer-Verlag, New York, 1991). T.H.H. Pian and K. Sumihara, Rational approach for assumed stress elements, Int. J. Numer. Methods Engrg. 20 (1984) 168551695. [S] J.C. Simo and T.J.F!. Hughes, On the variational foundations of assumed strain methods, J. Appl. Mech. ASME 53 (1986) 51-54. [6] T. Belytschko, Y.Y. Lu and L. Gu, Element free Galerkin methods, Int. J. Numer. Methods Engrg. 37 (1994) 229-256. [7] J.C. Simo and M.S. Rifai, A class of mixed assumed strain methods and the method of incompatible modes, Int. J. Numer. Methods Engrg. 29 (1990) 1595-1638. [S] C.A. Duarte and J.T. Oden, An h-p adaptive method using clouds, Comput. Methods Appl. Mech. Engrg. 139 (1996) 237-262. [9] W.K. Liu and Y. Chen, Wavelet and multiple scale reproducing kernel methods, Int. J. Numer. Methods Fluids 21 (1995) 901-932. [lo] W.K. Liu, S. Jun and Y.F. Zhang, Reproducing kernel particle methods, Int. J. Numer. Methods Fluids 20 (1995) 1081-1106. [ll] T. Hughes, The Finlte Element Method (Prentice-Hall, Englewood Cliffs, NJ, 1987). [12] K. Washizu, Variational Methods in Elasticity and Plasticity, 1st ed., Chs. 4 & 9 (Pergamon Press, 1968). [13] S. Timoshenko and J. Goodier, Theory of Elasticity, 3rd ed. (McGraw-Hill, 1970). [14] R. Cook, D. Malkus and M. Plesha, Concepts and Applications of Finite Element Analysis, 3rd ed. (John Wiley & Sons, 1989). [15] Y.Y. Lu, T. Belytschko and L. Gu, A new implementation of the EFG method, Comput. Methods Appl. Mech. Engrg. 113 (1994) 397414. [16] S. Timoshenko and S. Woinowsky-Krieger, Theory of Plates and Shells, 2nd ed. (McGraw-Hill, 1959). [17] L. Morley, Skew Plates and Structures (MacMillan, NY, 1963) 89-99. [18] H. Stolarski and T. Belytschko, Membrane locking and reduced integration for curved elements, J. Appl. Mech. 49 (1982) 172-177.