COMPUTER METHODS NORTH-HOLLAND
IN APPLIED
MECHANICS
AND ENGINEERING
63 (1987) 133-154
MIXED PETROV-GALERKIN METHODS FOR THE TIMOSHENKO BEAM PROBLEM* Abimael
F.D. LOULA,**
Division of Applied Mechanics,
Thomas J.R. HUGHES, Leopold0 Isidoro MIRANDA* **
P. FRANCA
and
Durand Building, Stanford University, Stanford, CA 94305, U.S.A.
Received 12 February
1987
A new mixed Petrov-Galerkin method is presented for the Timoshenko beam problem. The method has enhanced stability compared to the Galerkin formulation, allowing new combinations of interpolation, in particular, equal-order stress and displacement fields. The methodology is easily generalizable for multi-dimensional Hellinger-Reissner systems.
1. Introduction
Mixed finite element methods have been widely used in constrained problems. The introduction of Lagrange multipliers overcomes the severe limitation of finding interpolation functions satisfying the constraints a priori. But the construction of mixed finite element approximations is a problem of interpolation involving two spaces: the primal space (IV,) and the space of Lagrange multipliers (V,). And it is well known that, in order to satisfy convergence criteria (ellipticity and the Babuka-Brezzi (LBB) condition) we are limited to very few specific combinations of W,, and V,. In particular, a rather convenient choice-qualorder interpolation-is ruled out for large classes of problems. In this paper we address the question of constructing finite element approximations based on the Hellinger-Reissner principle for the Timoshenko beam. We first identify the difficulty in generating convergent approximations for this problem: lack of stability. Here, the limitation on the choice of the finite element spaces is not due to the LBB condition, but due to the lack of ellipticity of the bilinear form on the primal space. To circumvent this difficulty we enhance stability by perturbing the original Galerkin formulation. Our approach is consistent and preserves the classical mixed format of the variational equations. For this reason we denote it a mixed Petrov-Galerkin finite element method. The additional stability obtained in this method makes new choices of W,, and V, possible, in particular, equal-order interpolation which is unstable within the Galerkin approach. *This research was sponsored by the U.S. Office of Naval Research under Contract N00014-84-K-0600 (Scientific Officer: Dr. Alan Kushner). ** Present address: Laboratorio National de Computacgo Cientifica-LNCC/CNPq, Rua Laura Muller 455, 22290 Rio de Janeiro, Brazil. *** Present address: Escuela Superior de Ingenieros Industriales, Urdaneta 7, San Sebastian, Spain. 0045-7825/87/$3.50
0
1987, Elsevier Science Publishers
B.V. (North-Holland)
134
A. F. D. Loula
et al.,
Mixed
PG methods
for the Timoshenko
beam
In a different context, but with the same purpose of improving stability, some PetrovGalerkin methods have been successfully developed in fluid mechanics (see e.g. [4,7-141). In solid mechanics, despite the fact that a large number of problems have unstable Galerkin formulations, Petrov-Galerkin methods have not yet been used. Only recently, Loula, Hughes and Franca [16] presented a Petrov-Galerkin formulation for the Timoshenko beam problem in its displacement-rotation version based on the concept of an optimal weighting function which attains nodal exactness and optimal rate of convergence. Unfortunately that formulation does not seem easy to generalize. On the other hand, the method presented hcrc is rather general. In principle, it can be applied to any mixed formulation. An outline of the paper follows. In Section 2 we present an analysis of our model problem in its mixed variational form. and a summary of the main results on mixed finite element approximations based on Brezzi’s theory. In Section 3 we describe some facts about the standard Galerkin formulation. The new Petrov-Galerkin formulation is presented in Section 4, and numerical results and conclusions are given in Sections 5 and 6. respectively.
2. Preliminaries This section is devoted to a statement and analysis of the Timoshenko beam problem in the mixed form of the Hellinger-Reissner principle. We also present here the main results on the theory of mixed finite element methods. 2.1. Statement of the problem Let us consider the in-plane bending of a clamped uniform beam of length L, cross-section A, moment of inertia I, Young’s modulus E, and shear modulus G, subjected to a distributed load p(X), and a distributed moment m(x), with X E (0, L) representing the independent variable. According to Timoshenko system of differential equations:
-.-=dQ d.?
”
+Q=m;
-- Q KGA
beam theory, this problem is governed
by the following
(1)
(2)
+%0=0, dx
(3) (4)
with boundary
conditions:
w(O)=w(L)=O,
(3
O(O)= 8(L) = 0,
(6)
A. F. D. Loula et al., Mixed PG methods for the Timoshenko
135
beam
where Q(T) is the shear force; M(x) is the bending moment; w(X) is the transverse displacement; 19(x) is the cross-section rotation; and K is the shear correction factor. To explicate the dependence of this problem on a small parameter, EI
E2 =
KGAL~
(7)
’
the following change of variables:
we introduce
W ul=-’
u2=0,
L
(8)
(10) which reduces the original problem {a,(X), a,(x)}‘, x E (0,l) satisfying,
-a; =f, -a; -
to
finding
and
U(X) =
>
(11)
- a, = f, )
(12)
E2cT, + u; -
0
u2 =
(13)
)
-a,+u;=o, with boundary
U(X) = {U,(X), Us}’
(14)
conditions:
u,(O) = $(l)
=0)
(15)
u*(O) = U2(1) = 0 )
(16)
where the prime superscript denotes differentiation with respect to the nondimensional variable x = X/L. We observe that the nondimensional problem above depends explicitly on a parameter E, proportional to the ratio of thickness to length. In many applications E 4 1. In this limit the construction of finite element approximations is delicate. 2.2. Variational formulation Let W and V be the following Hilbert spaces:
w= CT= h >7J
E (L,(O, W2},
v= {u = {q, U2}tE (#(O,
l))‘} ) .
(17) (18)
136
A. F. D. Loula
et al.,
Mixed PG methods
for the Timoshenko
beam
functions on the unit interval with inner
where L,(O, 1) denotes the space of square-integrable product,
and &x0,1)
= {g E L,(O, 1); g’ E L,(O, 1); g(0) = g(1) = O} .
The norms on L,(O, l), H,!)(O, I), IV, and V are denoted
(2(j)
by, respectively.
Ilfll’ = (“6 f>
v.fE
Ildl:=(&s>+wd>
b~f4xOJL
ll~ll’w = II?II2+ lIdI
if7 E w ’
(22) (23)
IMIZ= IMI: + II%llf
VcJ Ev.
(24)
We now present problems
(II)-(
L(O, 1) ,
(21)
16) in an abstract variational
form:
Problem M: find (a, U) E W x V such that
w
a(a, 7) + b(7, u) = g(7)
VT E
b(a, rJ>= f(v)
YUEV,
where the bilinear forms a : W X Wa(a, 7) = - E2( @1r 5) b(7, u) = (7,)
u; -
R
(25)
)
(26)
and b
- (% 4
:W x
V+ IR are
VU,TEW,
u*) + (TX, u;> VT E
and f(u) and g(T) are given linear functionals.
(27)
w vu E v
(28)
)
In our problem,
with f= { fi, f2}’ E V’, the dual of V. We now show that the bilinear forms a and b in (27) and (28) are continuous,
l&T 41d II~IIwIl4v
vu, 7 E w ’
Ib(7,U)I~~~~T~~,)IUII~
VTEW
VuEV
since. (31)
3
(32)
A.F.D.
Loula et al., Mixed PG methods for the Timoshenko
beam
137
and
The abstract form (25)) (26) is shared by a wide class of problems. For instance, in elasticity, the symbols u, u represent the stress tensor and the displacement vector, respectively. The associated variational formulation is termed the Hellinger-Reissner principle [6,18]. From the mathematical point of view, Problem M fits in a general theory of abstract mixed problems studied by Brezzi [3], BabuSka [2], and many others. To make use of known general results we need to characterize the following subspaces of W: K(f) = (7 E W; b(T, u) = f(u)
Vu E V} ,
(35)
K = K(0) = Ker B = {T E W, b(T, u) = 0 Vu E V} ,
where B is the linear operator, (BT, u) = b(T, u)
(36)
B : WG V’,
VT E W Vu E V .
With respect to the subspace K(f)
(37)
we may define:
Problem P: find u E K(f) a(a, T) = g(7)
such that v7 E K )
which under the hypotheses e.g.,
of Theorem
(38)
2.1, presented
next, is equivalent to Problem M (see
PI).
2.3. Analysis of the problem Existence and uniqueness of solution for the class of problems governed by (25), (26) are assured by the following theorem due to Brezzi [3]. THEOREM
2.1. Zf there exist constants (Y > 0 and p > 0 such that
then Problem M has a unique solution (a, u) E W x V.
138
A. F. D. Loulu et al., Mixed PG methods for the Timoshenko
beum
We now prove that our problem satisfies hypotheses Hl and H2. For E > 0 in (27), the ellipticity condition (39) holds for the whole space W by choosing (Y= min(1, Ed). In the limit as F-+ 0, IV-ellipticity is lost. But it is important to note that from Theorem 2.1 we only need to verify (39) in K k W given by, . K={eW;
(~,,u;-u~)+(~~,u;)=O
VuEV),
(41)
which implies that every 7 E K must satisfy (r&)=0
vu, E #,(O, 1) ,
(42)
(72, u;> - (71, u2) = 0
v’v, E H:,(O,
(43)
1) .
The general solution of (42) is T, = my= constant.
Selecting u2 = ~ycp(x) in (43). with
cp(x) = 6x(1 - x) , (&),
1) = 1 1
(44) (45)
114112 = 12,
yields the following estimate, 1171112 = I/~;[l’d 12j/72/12 VT E K.
(46)
By definition, MT
711= ~211d12+ 1172112 ?=lIdI
which when combined
’
(47)
with (46) gives,
la(q+iljj/~11”w
VTEK,
which proves K-ellipticity with (Y= &. Hypothesis H2, the well-known BabuSka-Brezzi each u E V, a ? E W such that,
(48)
condition,
is easily verified by choosing for
_ 7r = u; )
(49)
fZ = u:‘. .
WV
For this choice, b(C v) G=4 Ilr$ ll~ll’, s Il4l;
)
,
(51)
(52)
A. F. D. Loula et al., Mixed PG methods for the Timoshenko
beam
139
and consequently,
(53) Thus conditions Hl and H2 are satisfied with (Y= &, and /? = $, independent of E. From Theorem 2.1 we conclude that Problem M has a unique solution ((T, U) E W x V for any E. 2.4. Finite element approximation results Let W,, and V, be finite element subspaces of W and V, respectively. define:
In these subspaces we
Problem M, : find (uh, u,,) E W,, x V, such that
4Uh7Th)+ Qh> Uh)= &h)
VThE Wh>
(54)
WJhTUh)=f@,>
vu, E v,;
(55)
and: Problem P,, : find ah E K,,(f) 4% ) T/l > = &I
such that 1
bl
(56)
E JGI
where
Kh(f)=bhEwh; Kh=Kh(0)={ThEWh;
4 'h,
'h)
=fkh>
"h
E 'h>
9
b(Th,Uh)=O Vu,EV,}.
(57) (58)
In general K,gK and Kh(f)gK(f), th us Problem Ph is typically a nonconforming approximation of Problem P. Therefore, hypotheses Hl and H2 of Theorem 2.1 must be verified for each particular finite element approximation. The cause of the poor performance of some apparently natural combinations of W, and V, is the inability to fulfill the discrete versions of hypotheses Hl and H2. The following theorems present some general results on convergence for Problems P, and M,. Their proofs can be found in, for example, [5]. THEOREM
2.2.
If there exists a constant (Y > 0 such that,
140
A. F. D. Loula
et ~1.. Mixed
then Problem P,, has a unique solution a,,
THEOREM
2.3.
PG methods
for the Titnoshenko
K,,(f)
and the following
E
helitt~
estimate holds:
If there exist constants (Y> 0 and B > 0 such that,
then Problem M,, has a unique solution (a,, , u,,) E W,7x V,, and the following
estimates hold: (63)
(64) Standard results from finite element interpolation theory and Theorems 2.2 and 2.3 yield error estimates for the discrete solution. Unfortunately, there is still the possibility of finding for particular finite element approximations (Y and/or p depending on the mesh parameter h. This will prevent convergence if the power of h from interpolation theory is not sufficiently large to overcome the division by (Y and/or p in (63) and (64).
3. Galerkin method The standard Galerkin finite element method consists of Problem M,, with particular choices of W, and V,. Let Qt be the space of C’ piecewise polynomial finite element interpolations of degree k and S:, = QL II H:,(O, 1) be the space of C” piecewise polynomials of degree 1. is the mesh parameter, h, is the diameter of the eth element. Here, h = max h,, 1 G e G y1,,, and IZ,, is the number of elements. We take W,, = (Q:)’ and V,,($,)“. 3.1.
Analysis
LEMMA
3.1 (Discrete BabGka-Brezzi
condition). Assume k b I- 1. Then
bh, uh)3 #kIIIUhllL’ vql E v,,, 7::;,, IIT,,IIw with B, a positive constant, independent of
E
and h.
(65)
A.F.D.
PROOF.
Loula et al., Mixed PG methods for the Timoshenko
beam
141
Since ul, E W,,, for each v, E V,, we can find a 7,, E W, such that
For this choice we have, as in the continuous
problem,
(68) which proves Lemma 3.1 with p = 4. LEMMA 3.2 (Characterization of the shear component, T,,~ of r,, = (Tag, T,,~) E K,,). Let pi, denote the mean value of the component 7h1 of r,, E K,, i.e., 1
0 Thl=
Define,
I
Th,
‘hl
-
o
* -
‘hl
dx .
0 ‘hl
(69)
(70)
.
Then, for any &, E QL ,
I where
T:,
PROOF.
h, o T:&dx=O,
e=l,2
is the restriction of
TX, E
,...,
n,,,
Q,” to element e, and p, is the restriction of CL,,to element e.
to the definition of Kh, any component
According
(71)
7h1
of 7h E Kh satisfies
1
I0
Th,U;,
dx = 0
(72)
vu,, E s:, .
It follows from (69) and (70) that, 1
(73) Given any ph E QL , there exists uhl E Si and a constant c such that (74)
u’el =pL:+cY
where u,* is the restriction
of vhl to element e. Substituting
(74) into (73) yields,
142
A. F. D. Loula et al., Mixed PG methods for the Timoshenko beam
By (69) and (70), (76)
and therefore (77) Invoking the independence of p: on each element, vanishes, which completes the proof of the lemma. COROLLARY element, that
3.3. It follows
it follows that each term in the sum
(77), by selecting pi, an arbitrary constant on each
from
he
I
0
rz,dx=O,
e=1,2
,...,
y1,,.
(78)
Remarks
(1) We want to emphasize that for this problem the BabuSka-Brezzi condition is trivially satisfied for k b I- 1. It is a simple consequence of choosing the stress discontinuous at element boundaries and displacements continuous. It is not with the inf-sup condition that one has difficulties in proving convergence, but with K,-ellipticity. (2) From Lemma 3.2 any shear component rl,, of 7h E K,, can be decomposed into two orthogonal components: rE1 and rz,. It is clear from (71) that if k 3 1, T:, is nontrivial. The component rz, is under control for any I, as shown in Section 2.3. The real problem here is the ri,‘s, often referred to as ‘spurious modes’. From (71) we obtain ~z, as shown in Fig. 1 for linear and quadratic interpolations. These plots show clearly how reduced integration filters the spurious modes. (3) One solution is to choose W, = (T,)’ with T,, = {‘hl
t
‘h
I =
$A
+
c
v$h
E
s;}
(79)
,
where SL is as defined before, and c is an arbitrary constant. Note that T, = &’ For this choice K,-ellipticity is immediate by substituting rh, = I,#~+ c in (72) and selecting u/,, = G,,, yielding rh, = c Consequently,
Vr, E K, .
(80)
rh*, = 0, and the coercivity of u( . , -) in K, follows:
with (Yindependent
of h and F. Therefore
for this approximation
we have from Theorem
2.3
A. F. D. Loula et al., Mixed PG methods for the Timoshenko
Fig. 1. ‘Spurious modes’: T:~ for linear (top) and quadratic
the following quasi-optimal
with c constant independent
(bottom)
beam
143
elements.
estimate:
of h and E.
3.2. Error estimates In what follows we consider the approximations of a, and uh as being of degree k and 1, respectively. By the usual interpolation theory (e.g., see [17]) we have the following. Given functions (T E (Hk+‘(O, 1))’ and u E (H’+‘(O, 1))’ there exist interpolants ah E IV, and u”hE V, such that,
IIU - a&/ =sCl(a)hk+l
lb - u”hll” s C*W’ 3
)
(83)
(84)
where c,(a) and C,(U) are independent of h. For the case k = I- 1, by choosing T,, = 0, and uh = L,, in (82) we have from (83) and (84):
A. F. D. Loula et al., Mixed PG methods for the Timoshenko
144
beam
where c~((T, u) and ~,(a, u) are independent of h and z:. Equations (X5) and (86) show that stresses and displacements converge with optimal rate. A similar result was obtained by Arnold [l] for the shear-only/ displacement mixed formulation.
4. Petrov-Galerkin
method
In this section we present a new finite element formulation for Problem M, termed a mixed Petrov-Galerkin finite element method, which preserves the good properties of the Galerkin method, presented in Section 3, but without its disadvantages. With I&‘,,and V,, as presented in Section 3.1, we define the following modification of Problem M,, . Problem PG, : find (a,, u,) E W, x V,, such that
in which 6 is a nondimensional
positive scalar parameter,
and ( - , * ),, is defined by
(91) where pC,and I,!I~are the restrictions of puI,and $J~to element e. To derive this formulation we view
as the residual form of the shear equilibrium equation, for the shear term, as we show in Lemma 4.2.
which provides the necessary srabifity
4.1. Analysis We first note that the new formulation is consistent in the sense that if (a, u) E W X V is the solution of Problem M, then (a, u) also satisfies Problem PG,, i.e.,
Mu, u,,)=f(u,,)
4,
E 4, .
(94)
A. F. D. Loula et al., Mixed PG methods for the Timoshenko
145
beam
We also note that W,,, V, and the bilinear form b : W,, x Vh += [ware the same as in Problem M. Consequently, the continuity of b, the LBB condition, and the characterization of K,, are preserved. Therefore Lemmas 3.1 and 3.2 are still valid. The K,-ellipticity of the new bilinear form ah : W, x W, + [wis established in the following lemmas. LEMMA
PROOF.
4.1
(Shear stability). With pi, as defined in Lemma 3.2.,
From (70) and (78) we can use the Poincare inequality in each element to get,
II
II6
e=1,2
,...,
It,,,
and the lemma follows by summing over elements. LEMMA that
4.2
(I&-ellipticity of a,,). There exists a constant cx > 0, independent of E and h, such
PROOF.
From the definition of a,, we have,
Using Lemma 4.1 yields,
with ri, as presented in Lemma 3.2. So far we have proved K,-ellipticity of a,, up to a global constant 7i1. To show that ril is also under control in K,, we proceed as we did in the continuous case, using the fact that any T,, E Kh must satisfy, (52, 42) - hd~ rJh2)= 0
vu,, E s’h
7
(99)
and selecting 0
‘h2
=
‘hl’ph
(100)
)
with (phE SL, such that
where p(x) is given by (44). For quadratic or higher-order
displacement
interpolations
(ph= cp.
146
A.F.D.
For linear interpolation
Loula
et al.,
Mixed PG methods
for the Timoshenko
beam
(ph# q, but the following estimates hold:
((Phr1)3(1-hZ),
lIdI = 12.
II4’~
Substituting
(102)
(103)
(100) into (99) and taking into account (102) and (103) we conclude that,
with c,, = a/( 1 - h’), for the linear interpolation. For quadratic and higher-order displacement interpolations (104) holds with c,, = v%!. The K,-ellipticity of all, with LYindependent of E and h, follows immediately from (98) and (104). It is interesting to note that in the linear interpolation case, due to (102), the finite element mesh must consist of at least two elements to satisfy inequality (96). The appropriate notion of continuity for the bilinear form a,, involves the mesh-dependent norm,
IIul(h,w = llallw +
sup *,st w,, 7il, # II
lM1r= (6 ’ 6 r
.
(
VUEW,
where
105)
( 106)
We note that, by definition.
l14.w 3 ll~ll1v ll~/,Ilh,W w + ~Mlw with n constant,
(107)
b7 E WI .
(108)
which is a consequence
h21&I;d r7h, II2 Vq,j E The required
V7EW,
continuity
condition
of the following inverse estimate: (109)
Qt
is,
I%(@? 0 ss~,Il4.,ll~/,IIU’
vu E w VT,,E w, ?
(110)
which follows from (109) with M, = max{ 1, Sfi}. Continuity of b : W x V-+ R is preserved with respect to the mesh-dependent norm by virtue of (32) and (107). The properties established in this section are sufficient to prove the following theorem. THEOREM 4.3. Assume W,, X V,, and the following
k 2 I - 1. Then Problem error estimates hold,
PG,, has a unique
solution
(a,, . u,,) E
lb - ~/JIh,W + lb - 4” =s4b - TA,.U~ + IIU- %ll”) VT, E w, vu/lE v,z9 (111) with c independent
of E and h.
A. F. D. Loula et al., Mixed PG methods for the Timoshenko beam
147
The proof of Theorem 4.1 is long, and therefore is omitted. Invoking standard interpolation results, it follows from (111) that,
II~-o;,l/,,~+Il~-~,ll”=~~~‘~~
(112)
is also O(h’). This estimate indicates that for all k a 1 - 1 displacements BY W), Ilo - ~hllW converge with the optimal rate in the V norm, for k = I- 1 stresses converge with optimal rates in the W norm, as in the Galerkin method, and for k = 1, stresses converge with suboptimal rate with gap 1 in the W norm. These error estimates are independent of E. Therefore equal-order interpolation also ‘works’ for thin beams in the new formulation. Remarks
(1) In our problem K C Kh, but this inclusion does not occur in general. The additional term in the bilinear form ah provides stability in the subset K* = K,, - K,, fl K but not in K” = K,, n K. If the bilinear form a : W X W+ [w is K-elliptic, a,, : W, x W,, -+ [w is also K,-elliptic and stability is guaranteed. (2) Since we use discontinuous interpolation for TV, for E > 0 we can always eliminate 7, at element level to generate a displacement method. Assuming f,(x) = q = constant, f*(x) = 0 and equal-order linear interpolation, after eliminating the ~=‘s we obtain
(
I&+&&
>
E2
d=F,
where d is the displacement vector, F is the force vector, K, is the bending stiffness, identical to that obtained in the Galerkin method, and K, is the shear stiffness. The element contributions are,
K’, = l/k
7
K; =
llh l/2 -l/h l/2
l/2 c_xh
-l/2 Ph
-l/h -l/2 llh -l/2
l/2 ph
-l/2 ah
and Fe=
qh
Sqh*
i
2
126 + E2
(114)
i- 01
where (Y=
P= When ~*+a, (Y +
96 + E2 36s + 3~~ ’
(115)
186 + E* 726 + 6~* .
(116)
from (115), (116), a,
P +
a,
A. F. D. Loulu et ul., Mixed PG methods ,for the Timoshenko
148
hecm
and the matrix Kz reduces to the well-known one-point quadrature shear matrix for the standard Galerkin method. We also observe an extra loading term in F”. This term does not appear in the standard Galerkin method and it is responsible for a better approximation when the element distributed loading is nonzero. It is also interesting to note that when F’ e 8, 8”’ reduces to the consistent load for Bernoulli-Euler beam theory with Hermite cubic interpolation. (3) The mixed formulation of the Timoshenko beam problem in terms of shear resultant, displacement and rotation, as studied by Arnold [l], was also studied by the authors. Another Petrov-Galerkin formulation was developed by selecting as primal variables, */, =
and Lagrange
CT,,,u,J E 4, = 3
Q:,x $, *
(117)
multiplier,
In this formulation
we have the following
bilinear
forms,
and the same expressions for ~(zJ,~) and g,?(~,~) as those in Problem PC;,, . Excellent results were obtained for this formulation as well. The analysis is slightly different be presented together with the numerical results in future work.
numerical and it will
5. Numerical results In this section we report on some representative tests we performed with the mixed Petrov-Galerkin method using equal-order interpolation elements with u discontinuous and u continuous. In all cases presented we took E = 1OY’. The contribution to stability of the term 6 is illustrated in Figs. 2 and 3 for different physical situations. In both figures we plot the a,,, (proportional to shear force) and u,,* (proportional to bending moment) obtained for the linear element for different values of 6. tn all cases spurious oscillations in the shear component disappear for 6 sufficiently large. For S > 1 the solutions are the same as the solution for 6 = 1, and therefore are not plotted. In Fig. 2 we have a beam with both ends clamped (boundary conditions of (15), (16)) subjected to three different loads: unit load at the midpoint; unit uniform load; linear load varying from 0 to 1. Note the excellent results for just 4 elements in both cases when S = 1. In Fig. 3 we have a clamped beam subjected to unit tip displacement for discretizations of 4 and 16 elements. The solution for S = 1 is free from spurious oscillations and approaches the exact solution as we refine the mesh. These results confirm the important role of the term 8 (see (91)) in the stabilization of an
A.F. D. Loula et al., Mixed PG methods for the Timoshenko
N
I
Midspan unit load
Midspan unit load
3k ....
d
,.: A:
‘.
,_:
‘. . .
2
..I’
._ ‘.
;
00
“‘A.”
.A..
_:’
q
..
.:’
: ,::’ :.: A
exact
66 ..
t_ ‘...
El '.
‘..
E2 M
d
6=1.
A
6=l.e-03 I 0.25
A
N
I
I
0.50
0.75
6=l.e-03
dl I
0
0 6=1.
I
0
“:
X(,
I
I
1
I
0
I
0.25
I
0.50
/
0.75
I 1
x
X
N
149
beam
I
I
1
I
Uniform unit load
I
Uniform unit load
2
aJ0
2 3
I
N I
0
6=1.
A
d=l.e-3
0
A
I
I
I
0.25
0.50
0.75
1
6=l.e-3
0
0.25
/
I
I
I
0
6=l.e-3
0
exact 6=1.
A
6=l.e-3
I
I
0.25
0.25
Fig. 2. Shear
0.50 x
(ukI) and moment
0.75
I
Linear load
Linear load
A
0.50
X
X
0.75
1
0
u.50
0.75
X (@,,h2)along the axis for a beam with both ends clamped
(4 elements).
1
A. F. D. Loula et al., Mixed PG methods for the Timoshenko
150
A
0 0
I
0
0.25
1
0.75
0.50
0
0.25
2I
1
0.50
0.75
1
6=l.e-03
/ 0
0.25
X Fig. 3. Shear elements.
/
0.75
6=l.e-03
A 0.25
/
0.50
X
0
0
I
6=l.e-3
X
A
beam
(cT~,) and moment
0.50
0.75
1
X (cT~_) along the axis for a beam with unit tip displacement.
Results
for 4 and 16
unstable Galerkin method. Now we show that our approach gives improved results even when applied to an originally stable formulation. To this end we introduce an additional perturbation in (87) related to the moment equilibrium equation. In this Petrov-Galerkin method WC have, a/@,>
Q) = a(ajl? 7,) - 6h2(6,,
3 T;,1I,, -
&A4
= 6h2(f,>
d2
4,>,,
+
Vh2(f2>
+
?71)h
J442 )
+ a/,, > &
+ T/l, ),, 7
(121) ( 122)
in place of Problem PG, . Of course this additional perturbation is also consistent, and since the K,-ellipticity and the discrete LBB condition are verified independently of h, the estimate (112) is still valid. Figure 4 shows significant improvement in the approximation of the moment due to the v-term. We also observe a gain in the shear accuracy, when possible. In Fig. 5 we perform an error study for the beam with both ends clamped subjected to unit
A. F. D. Loula et al., Mixed PG methods for the Timoshenko
I
I
151
beam
1
Unit displacement
at tip
Unit displacement
at tip
exact 0
6=1., v=o.
A
ckl.,
v=l.
0
a
:
I
I
I
0.25
0.50
0.75
1
0
0.50
0.25
X
Midspan
0.75
1
0.75
1
X
unit load
Midspan
0
0.25
unit load
0.50 x
Uniform
Uniform
unit load
unit load
1
exact
0
0
6=1., v=o.
0
d=l.,
v=o.
A
d=l.,
A
6=1.,
v=l. /
v=l. 0.25
0.50
0.75
0
1
0.25
x
Fig. 4. Comparison additional parameter
of shear V.
I
I
0.50
0.75
-I
1
x
(uhI)
and
moment
(c,,*)
along
the
axis for the
Petrov-Galerkin
method
with
A. F. D. Lo&
et al., Mixed PG methods for the Timoshenko
beam
0 = linear 0 = quadratic
“: b
0 = linear 0 = quadratic
0 = linear 0 = quadratic
0 = linear s 1 0 = ,quadr;tic
0
0
1
3
2
0
log N elements
0 = linear 0 = quadratic
0
/
I
1
2
log N elements
1
2
3
7
0
I
2
log N elements
0 = linear 0 = quadratic
0 = linear 0 = quadratic
0 3
I
log N elements
3
0
1
1
2
3
log N elements
Fig. 5. Log of the error of the variables and its derivatives both ends clamped subjected to a uniform load.
0
1
2
3
I
0
I
I
1
2
log N elements
log N elements
versus log of the number
3
0 = linear 0 = quadratic 0
0
2
log N elements
of elements
for a beam with
3
A.F.D.
Loula et al., Mixed PG methods for the Timoshenko
beam
153
uniform load for 6 = 1. We plot the &-norm of the error of each variable and its derivative for both linear and quadratic elements. We obtain optimal rate for each displacement component and suboptimal rate with gap 1 for each stress component as predicted in the error estimates in Section 4.2. Confirmation of the error analysis was also obtained for other boundary conditions and loading cases.
6. Conclusions Stated in Hellinger-Reissner form, classical elasticity problems usually do not have difficulty in fulfilling the BabuSka-Brezzi condition for both continuous and discrete formulations. In this framework the dijjkulty is the satisfaction of the ellipticity condition in parameterdependent problems, such as nearly-incompressible elasticity, thin beams, and thin structures in general. It is well known that for the standard Galerkin method stability of the finite element discretization is often a question of interpolation, i.e., the choice of W,, and V,, is crucial. This is, of course, a serious limitation for the use of mixed finite element methods. In this paper we presented a new approach for overcoming this difficulty. By perturbing the Galerkin method we recover stability for various combinations of W, and V, , even equal-order interpolation. Though presented for a simple model problem we have already extended this methodology successfully to compressible and nearly-incompressible elasticity (including the Stokes problem), arch problems, axisymmetric shells, etc. These will be presented in forthcoming papers.
Acknowledgment
During fellowship Brazilian supported
the course of this work Abimael F.D. Loula was partially supported by a CNPq from the Brazilian Government, Leopold0 P. Franca was partially supported by Government fellowship CAPES Proc. 4093/82-5 and Isidoro Miranda was partially by a Fulbright/MEC fellowship from the Spanish Government.
References [l] D.N. Arnold, Discretization by finite elements of a model parameter dependent problem, Numer. Math. 37 (1981) 405421. [2] I. BabuSka, Error bounds for finite element method, Numer. Math. 16 (1971) 322-333. [3] F. Brezzi, On the existence, uniqueness and approximation of saddle-point problems arising from Lagrange multipliers, RAIRO Anal. NumCr. R-2 (1974) 129-151. . [4] A.N. Brooks and T.J.R. Hughes, Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations, Comput. Meths. Appl. Mech. Engrg. 32 (1982) 199-259. [5] V Girault and P.A. Raviart, Finite Element Approximation of the Navier-Stokes Equations, Lecture Notes in Mathematics 749 (Springer, Berlin, 1979).
154
A. F. D. Loula et al., Mixed PG methods for the Timoshenko
beam
[6] E. Hellinger, Der Allgemeine Ansatz der Mechanik der Kontinua. in: Encyclopadie der Mathematischen Wissenschaften 4 (4) (1914) 602-694. [7] T.J.R. Hughes and A.N. Brooks, A theoretical framework for Petrov-Galerkin methods with discontinous weighting functions: Application to the streamline upwind procedure, in: R.H. Gallagher et al.. eds., Finite Elements in Fluids IV (Wiley, London, 1982) 46-65. [8] T.J.R. Hughes, L.P. Franca and M. Balestra, A new finite element formulation for computational fluid dynamics: V. Circumventing the BabuSka-Brezzi condition: A stable Petrov-Galerkin formulation of the Stokes problem accommodating equal-order interpolations, Comput. Meths. Appl. Mech. Engrg. 59 (1986) 85-99. [9] T.J.R. Hughes and M. Mallet, A new finite element formulation for computational fluid dynamics: III. The generalized streamline operator for multidimensional advection-diffusion systems, Comput. Meths. Appl. Mech. Engrg. 58 (1986) 305-328. [lo] T.J.R. Hughes and M. Mallet, A new finite element formulation for computational fluid dynamics: IV. A discontinuity-capturing operator for multidimensional advective-diffusive systems. Comput. Meths. Appl. Mech. Engrg. 58 (1986) 329-336. [il] T.J.R. Hughes, M. Mallet and L.P. Franca, Entropy-stable finite element methods for compressible fluids: Application to high Mach number flows with shocks, in: P.G. Bergan, K.J. Bathe and W. Wunderlich, eds.. Finite Element Methods for Nonlinear Problem (Springer, Berlin, 1985) 761-773. [12] T.J.R. Hughes, M. Mallet and L.P. Franca. New finite element methods for the compressible Euler and Navier-Stokes equations, in: R. Glowinski and J.L. Lions, eds., Computing Methods in Applied Sciences and Engineering VII (North-Holland. Amsterdam, 1987) 339-360. [13] C. Johnson, Streamline diffusion methods for problems in fluid mechanics, in: R.H. Gallagher et al.. eds., Finite Elements in Fluids VI (Wiley, London, 1986) 251-261. [14] C. Johnson and Saranen, Streamline diffusion methods for the incompressible Euler and Navier-Stokes equations, Math. Comput. 47 (1986) 1-18. [15] C. Johnson, U. NCvert and J. Pitkaranta, Finite element methods for linear hyperbolic problems, Comput. Meths. Appl. Mech. Engrg. 45 (1985) 285-312. [16] A.F.D. Loula, T.J.R. Hughes and L.P. Franca. Petrov-Galerkin formulations of the Timoshenko beam problem, Comput. Meths. Appl. Mech. Engrg. 63 (1987) 115-132 (this issue). [17] J.T. Oden and G.F. Carey, Finite Elements: Mathematical Aspects IV (PrenticeHall, Englewood Cliffs, NJ, 1983). [18] E. Reissner, On a variational theorem in elasticity. J. Math. Phys. 2Y (2) (1950) YO-YS.