11

Iso-P Quadrilateral
Ring Elements

11–1
Chapter 11: ISO-P QUADRILATERAL RING ELEMENTS 11–2

TABLE OF CONTENTS

Page
§11.1. INTRODUCTION 11–3
§11.2. ORGANIZATION OF ELEMENT MODULES 11–3
§11.3. TWO-DIMENSIONAL 2D GAUSS RULES 11–4
§11.4. THE 4-NODE QUADRILATERAL RING ELEMENT 11–5
§11.4.1. Shape Function Evaluation . . . . . . . . . . . . . 11–5
§11.4.2. Element Stiffness . . . . . . . . . . . . . . . 11–6
§11.4.3. Stiffness Test . . . . . . . . . . . . . . . . . 11–7
§11.4.4. Test of Rectangular Geometry . . . . . . . . . . . 11–7
§11.4.5. Consistent Node Forces for Body Force Field . . . . . . . 11–9
§11.4.6. Recovery of Element Corner Stresses . . . . . . . . . 11–11
§11.5. THE 9-NODE QUADRILATERAL RING ELEMENT 11–14
§11.6. THE 16-NODE QUADRILATERAL RING ELEMENT 11–14
EXERCISES . . . . . . . . . . . . . . . . . . . . . . 11–15

11–2
11–3 §11.2 ORGANIZATION OF ELEMENT MODULES

§11.1. INTRODUCTION
This Chapter illustrates the computer implementation of isoparametric quadrilateral elements for
the axisymmetric problem. These are called ring elements. Triangles, which present some pro-
gramming quirks, are covered in the next Chapter.
The Chapter covers 4-node, 9-node and 16-node quadrilaterals, with most detail devoted to the
4-node model. It describes the computation of the element stiffness matrix, consistent node force
vector for a body force field and recovery of element stresses from displacements.

§11.2. ORGANIZATION OF ELEMENT MODULES

We consider the implementation of the 4-node blinear quadrilateral for axisymmetric analysis. The
element cross section is depicted in Figure 11.1.

η
z 3 (r 3 , z 3)
η=1
4 (r 4 , z 4)
ξ=1

ξ=−1 ξ

1 (r 1 , z 1) η=−1
2 (r 2 , z 2)

Figure 11.1. The 4-node bilinear quadrilateral ring element.

For 2D and 3D elements iso-P elements it is convenient to break up the implementation into
application dependent and application independent modules, as sketched in Figure 11.2. The
application independent modules can be “reused” in other FEM applications, for example to form
thermal, fluid or electromagnetic elements.
For the 4-node quadrilateral studied here, the subdivision of Figure 11.2 is done through the fol-
lowing modules:

Stiffness4NodePlaneStressQuad - forms Ke of 4-node bilinear quad in plane stress

QuadGaussRuleInfo - returns Gauss quadrature product rules of order 1-4
IsoQuad4ShapeFunctions - evaluates shape functions and their x/y derivatives
These modules are described in further detail in the following subsections, in a “bottom up” fashion.

11–3
Chapter 11: ISO-P QUADRILATERAL RING ELEMENTS 11–4

APPLICATION DEPENDENT

Element
Stiffness Module

Shape Gauss Quadrature

Function Module Information Module

APPLICATION INDEPENDENT

Figure 11.2. Organization of element stiffness modules.

QuadGaussRuleInfo[{rule_,numer_},point_]:= Module[
{ξ ,η,p1,p2,i1,i2,w1,w2,k,info=Null},
If [Length[rule]==2, {p1,p2}=rule,p1=p2=rule];
If [Length[point]==2,{i1,i2}=point,
k=point; i2=Floor[(k-1)/p1]+1; i1=k-p1*(i2-1) ];
{ξ ,w1}= LineGaussRuleInfo[{p1,numer},i1];
{η,w2}= LineGaussRuleInfo[{p2,numer},i2];
info={{ξ ,η},w1*w2};
If [numer,Return[N[info]],Return[Simplify[info]]];
];
LineGaussRuleInfo[{rule_,numer_},point_]:= Module[
{g2={-1,1}/Sqrt[3],w3={5/9,8/9,5/9},
g3={-Sqrt[3/5],0,Sqrt[3/5]},
w4={(1/2)-Sqrt[5/6]/6,(1/2)+Sqrt[5/6]/6,
(1/2)+Sqrt[5/6]/6,(1/2)-Sqrt[5/6]/6},
g4={-Sqrt[(3+2*Sqrt[6/5])/7],-Sqrt[(3-2*Sqrt[6/5])/7],
Sqrt[(3-2*Sqrt[6/5])/7],Sqrt[(3+2*Sqrt[6/5])/7]},
i=point,info=Null},
If [rule==1,info={0,2}];
If [rule==2,info={g2[[i]],1}];
If [rule==3,info={g3[[i]],w3[[i]]}];
If [rule==4,info={g4[[i]],w4[[i]]}];
If [numer,Return[N[info]],Return[Simplify[info]]];
];

Figure 11.3. Module to get Gauss-product quadrature information for a quadrilateral.

§11.3. TWO-DIMENSIONAL 2D GAUSS RULES

Module QuadGaussRuleInfo, listed in Figure 11.3, implements the two-dimensional product
Gauss rules with 1 through 4 points in each direction. The number of points in each direction may
be the same or different.
If the rule has the same number of points p in both directions the module is called in either of two
ways:

{ { xii,etaj },wij }=QuadGaussRuleInfo[{ p,numer }, { i,j }]

(11.1)
{ { xii,etaj },wij }=QuadGaussRuleInfo[{ p,numer }, k ]

11–4
11–5 §11.4 THE 4-NODE QUADRILATERAL RING ELEMENT

The first form is used to get information for point {i, j} of the p × p rule, in which 1 ≤ i ≤ p and
1 ≤ j ≤ p. The second form specifies that point by a “visiting counter” k that runs from 1 through
p 2 ; if so {i, j} are internally extracted1 as j=Floor[(k-1)/p]+1; i=k-p*(j-1).
If the integration rule has p1 points in the ξ direction and p2 points in the η direction, the module
may be called also in two ways:

{ { xii,etaj },wij }=QuadGaussRuleInfo[{ { p1,p2 },numer }, { i,j }]

(11.2)
{ { xii,etaj },wij }=QuadGaussRuleInfo[{ { p1,p2 },numer }, k ]

The meaning of the second argument is as follows. In the first form i runs from 1 to p1 and j from 1 to
p2 . In the second form k runs from 1 to p1 p2 ; if so i and j are extracted by j=Floor[(k-1)/p1]+1;
i=k-p1*(i-1).
In all four forms flag numer is set to True if numerical information is desired and to False if exact
information is desired. The module returns ξi and η j in xii and etaj, respectively, and the weight
product wi w j in wij.

§11.4. THE 4-NODE QUADRILATERAL RING ELEMENT

Quad4IsoPShapeFunDer[ncoor_,qcoor_]:= Module[
{Nf,dNr,dNz,dNξ ,dNη,i,J11,J12,J21,J22,Jdet,ξ ,η,
r1,r2,r3,r4,z1,z2,z3,z4,r,z},
{ξ ,η}=qcoor; {{r1,z1},{r2,z2},{r3,z3},{r4,z4}}=ncoor;
Nf={(1-ξ )*(1-η),(1+ξ )*(1-η),(1+ξ )*(1+η),(1-ξ )*(1+η)}/4;
dNξ ={-(1-η), (1-η),(1+η),-(1+η)}/4;
dNη= {-(1-ξ ),-(1+ξ ),(1+ξ ), (1-ξ )}/4;
r={r1,r2,r3,r4}; z={z1,z2,z3,z4};
J11=dNξ .r; J12=dNξ .z; J21=dNη.r; J22=dNη.z;
Jdet=Simplify[J11*J22-J12*J21];
dNr= ( J22*dNξ -J12*dNη)/Jdet; dNr=Simplify[dNr];
dNz= (-J21*dNξ +J11*dNη)/Jdet; dNz=Simplify[dNz];
Return[{Nf,dNr,dNz,Jdet}]
];

Figure 11.4. Shape function module for 4-node bilinear quadrilateral.

§11.4.1. Shape Function Evaluation

Quad4IsoPShapeFunDer is an application independent module that computes the shape functions
Ni(e) , i = 1, 2, 3, 4 and its global partial derivatives at a point of coordinates {ξ, η}. The module
is listed in Figure 11.4. It is identical to that used in the IFEM course except for relabeling of
the coordinates {x, y} to {r, z}. The arguments of the module are the {x, y} quadrilateral corner
coordinates, which are passed in ncoor, and the two quadrilateral coordinates {ξ, η}, which are
passed in qcoor.

1 Indices i and j are denoted by i1 and i2, respectively, inside the module.

11–5
Chapter 11: ISO-P QUADRILATERAL RING ELEMENTS 11–6

Quad4IsoPRingStiffness[ncoor_,Emat_,Kfac_,options_]:=
Module[{k,p=2,numer=False,qcoor,w,c,N1,N2,N3,N4,
dNr1,dNr2,dNr3,dNr4,dNz1,dNz2,dNz3,dNz4,Jdet,
r1,r2,r3,r4,z1,z2,z3,z4,rk,Be,Ke=Table[0,{8},{8}],
modname="Quad4IsoPRingStiffness: "},
If [Length[options]==2,{numer,p}=options,{numer}=options];
If [p<1||p>4,Print[modname,"p out of range"];Return[Null]];
{{r1,z1},{r2,z2},{r3,z3},{r4,z4}}=ncoor;
For [k=1,k<=p*p,k++,
{qcoor,w}= QuadGaussRuleInfo[{p,numer},k];
{{N1,N2,N3,N4},{dNr1,dNr2,dNr3,dNr4},{dNz1,dNz2,
dNz3,dNz4},Jdet}=Quad4IsoPShapeFunDer[ncoor,qcoor];
rk=r1*N1+r2*N2+r3*N3+r4*N4; c=Kfac*w*Jdet*rk;
Be={{dNr1, 0,dNr2, 0,dNr3, 0,dNr4, 0},
{ 0,dNz1, 0,dNz2, 0,dNz3, 0,dNz4},
{N1/rk, 0,N2/rk, 0,N3/rk, 0,N4/rk, 0},
{dNz1,dNr1,dNz2,dNr2,dNz3,dNr3,dNz4,dNr4}};
If [numer,Be=N[Be]];
Ke+=Simplify[c*Transpose[Be].(Emat.Be)];
]; Return[Ke]
];

Figure 11.5. Element stiffness formation module for 4-node iso-P quadrilateral ring.

The quadrilateral coordinates define the element location at which the shape functions and their
derivatives are to be evaluated. For the stiffness formation these are Gauss points, but for strain and
stress computations these may be other points, such as corner nodes.
Quad4IsoPShapeFunDer returns as function value the two-level list { Nf,dNr,dNz,Jdet }, in
which the first three are 4-entry arrays. List Nf collects the shape function values, dNr the shape
function r -derivatives, dNz the shape function z-derivatives, and Jdet is the Jacobian determinant.

§11.4.2. Element Stiffness

Module Quad4IsoPRingStiffness computes the stiffness matrix of a four-noded isoparametric

quadrilateral ring element to model an axisymmetric solid. The module logic is listed in Figure
11.5. The logic follows the procudure described in the previous chapter.
The arguments of the module are:
ncoor Quadrilateral node coordinates arranged in two-dimensional list form:
{ { r1,z1 },{ r2,z2 },{ r3,z3 },{ r4,z4 } }.
Emat The 4 × 4 matrix of elastic moduli:
 
E 11 E 12 E 13 E 14
E E 22 E 23 E 24 
E =  12 , (11.3)
E 13 E 23 E 33 0
E 13 E 23 0 E 44

Note that E 34 = 0 to satisfz axisymmetry assumptions. If the material is

11–6
11–7 §11.4 THE 4-NODE QUADRILATERAL RING ELEMENT

isotropic:
 
1+ν ν ν 0
E  ν 1 ν 0 
E=  . (11.4)
(1 + ν)(1 − 2ν) ν ν 1+ν 0
0 0 0 1
2
−ν

Kfac A ring-span factor by which the stiffness matrix will be scaled. Typically
Kfac=1 to take the ring element to span one radian, Kfac=2π to take a span of
a complete circle.
options Processing options. This list may contain two items: { numer,p } or one:
{ numer }.
numer is a logical flag with value True or False. If True, the computa-
tions are forced to proceed in floating point arithmetic. For symbolic or exact
arithmetic work set numer to False.
p specifies the Gauss product rule to have p points in each direction. It may
be 1 through 4. For rank sufficiency, p must be 2 or higher. If p is 1 the element
will be rank deficient by two. If omitted p = 2 is assumed.
The module returns Ke as an 8 × 8 symmetric matrix pertaining to the following arrangement of
nodal displacements:

u(e) = [ u r 1 u z1 ur 2 u z2 ur 3 u z3 ur 4 u z4 ]T . (11.5)

§11.4.3. Stiffness Test

The stiffness module is tested on the rectangular geometry identified in Figure 11.6. It is a rectangle
dimensioned a × b with sides parallel to the {r, z} axes. The distance of the leftmost side to the z
axis is d. The material is isotropic with modulus E and Poisson’s ratio ν.

Isotropic material
modulus E, z
Poisson's ratio ν
z 4 3 Ring element

1 2
d a

Figure 11.6. Test quadrilateral ring element geometry.

11–7
Chapter 11: ISO-P QUADRILATERAL RING ELEMENTS 11–8

Quad4IsoPRingBodyForces[ncoor_,Emat_,Kfac_,options_,bfor_]:=
Module[{k,m,p=2,numer=False,qcoor,w,c,N1,N2,N3,N4,dNr,dNz,
Jdet,r1,r2,r3,r4,z1,z2,z3,z4,rk,br,bz,br1,bz1,br2,bz2,
br3,bz3,br4,bz4,brc,bzc,bk,fe=Table[0,{8}],
modname="Quad4IsoPRingBodyForces: "}, m=Length[bfor];
If [Length[options]==2,{numer,p}=options,{numer}=options];
If [p<1||p>4,Print[modname,"p out of range"];Return[Null]];
If [m==2,{br,bz}=bfor;br1=br2=br3=br4=br;bz1=bz2=bz3=bz4=bz];
If [m==4,{{br1,bz1},{br2,bz2},{br3,bz3},{br4,bz4}}=bfor];
{{r1,z1},{r2,z2},{r3,z3},{r4,z4}}=ncoor;
For [k=1,k<=p*p,k++,
{qcoor,w}= QuadGaussRuleInfo[{p,numer},k];
{{N1,N2,N3,N4},dNr,dNz,Jdet}=Quad4IsoPShapeFunDer[ncoor,qcoor];
rk=r1*N1+r2*N2+r3*N3+r4*N4; c=Kfac*w*Jdet*rk;
brk=br1*N1+br2*N2+br3*N3+br4*N4;
bzk=bz1*N1+bz2*N2+bz3*N3+bz4*N4;
bk={N1*brk,N1*bzk,N2*brk,N2*bzk,N3*brk,N3*bzk,N4*brk,N4*bzk};
If [numer,bk=N[bk]]; fe+=c*bk;
]; Return[Simplify[fe]]
];

Figure 11.8. Driver for stiffness calculation of trapezoidal element

of Figure 11.6(b) for four Gauss integration rules.

§11.4.4. Test of Rectangular Geometry

The script of Figure 11.8 compute and print the stiffness of the rectangular element shown in Figure
11.6, for E = 96, ν = 1/3, a = 4, b = 2 d = 0 and Kfac=1. Nodes 1 and 2 sit on the z axes. The
value of p is changed in a loop. The flag numer is set to True to use floating-point computation for
speed (see Remark 11.1). The computed entries of K(e) are exact integers for all values of p:

 
72 18 36 −18 −36 −18 0 18
 18 153 −54 135 −90 −153 −18 −135 
 36 −54 144 −90 54 −36 90 
 72 
 −18 −90 −54 −135 −153 
= 
135 153 18
K(e)  −36 −90
1×1
 72 −54 144 90 36 54 

 −18 −153 54 −135 90 153 18 135 
 0 −18 −36 18 36 18 72 −18

18 −135 90 −153 54 135 −18 153
 
168 −12 24 12 −24 −36 48 36
 −12 108 −24 84 −72 −102 −36 −90 
 24 −24 −120 72 −24 72 
 216 0 
 12 84 −120 −72 −282 36 −102 
= 
300
K(e)  −24 −72
2×2
 0 −72 216 120 24 24 

 −36 −102 72 −282 120 300 −12 84 
 48 −36 −24 36 24 −12 168 12

36 −90 72 −102 24 84 12 108

11–8
11–9 §11.4 THE 4-NODE QUADRILATERAL RING ELEMENT
 
232 −12 24 12 −24 −36 80 36
 −12 108 −24 84 −72 −102 −36 −90 
 24 −24 216 −120 72 −24 72 
 0 
 12 84 −120 −72 −282 36 −102 
= 
300
K(e)  −24 −72
3×3
 0 −72 216 120 24 24 

 −36 −102 72 −282 120 300 −12 84 
 80 −36 −24 36 24 −12 232 12

36 −90 72 −102 24 84 12 108
 
280 −12 24 12 −24 −36 104 36
 −12 108 −24 84 −72 −102 −36 −90 
 24 −24 216 −120 72 −24 72 
 0 
 12 −120 −72 −282 36 −102 
= 
84 300
K(e)  −24 −72
4×4
 0 −72 216 120 24 24 

 −36 −102 72 −282 120 300 −12 84 
 104 −36 −24 36 24 −12 280 12

36 −90 72 −102 24 84 12 108

As can be seen entries change substantially in going from p = 1 to p = 2. From then on only four
entries, associated with the r stiffness at nodes 1 and 4, change. The eigenvalues of these matrices
are:

Rule Eigenvalues of K(e) for varying integration rule

1×1 667.794 180.000 124.206 72.000 0 0 0 0
2×2 745.201 261.336 248.750 129.451 100.389 88.598 10.275 0
3×3 745.446 330.628 266.646 133.236 126.343 98.690 11.011 0
4×4 745.716 397.372 272.092 144.542 135.004 101.908 11.365 0 (11.6)

The stiffness matrix computed by the one-point rule is rank deficient by three. For p = 2 and up it has
the correct rank of 7. The eigenvalues do not change appreciably after p = 2. Because the nonzero
eigenvalues measure the internal energy taken up by the element in deformation eigenmodes, it can
be seen that raising the order of the integration stiffens the element.

§11.4.5. Consistent Node Forces for Body Force Field

The module Quad4IsoPRingBodyForces listed in Figure 11.8 computes the consistent force
associated with a body force field b = {bx , b y } given over a four-node iso-P quadrilateral in plane
stress. The field is specified per unit of volume in componentwise form.
For example if the element is subjected to a gravity acceleration field g due to self-weight in the
−z direction, br = 0 and bz = −ρg, where ρ is the mass density. If the element rotates at constant
angular velocity ω (radians per second) around z, thenbr = −ρω2r and bz = 0.
The arguments of the module are:
ncoor As in Quad4IsoPRingStiffness
Emat Not used: a dummy argument kept as placeholder.

11–9
Chapter 11: ISO-P QUADRILATERAL RING ELEMENTS 11–10

Quad4IsoPRingBodyForces[ncoor_,Emat_,Kfac_,options_,bfor_]:=
Module[{k,m,p=2,numer=False,qcoor,w,c,N1,N2,N3,N4,dNr,dNz,
Jdet,r1,r2,r3,r4,z1,z2,z3,z4,rk,br,bz,br1,bz1,br2,bz2,
br3,bz3,br4,bz4,brc,bzc,bk,fe=Table[0,{8}],
modname="Quad4IsoPRingBodyForces: "}, m=Length[bfor];
If [Length[options]==2,{numer,p}=options,{numer}=options];
If [p<1||p>4,Print[modname,"p out of range"];Return[Null]];
If [m==2,{br,bz}=bfor;br1=br2=br3=br4=br;bz1=bz2=bz3=bz4=bz];
If [m==4,{{br1,bz1},{br2,bz2},{br3,bz3},{br4,bz4}}=bfor];
{{r1,z1},{r2,z2},{r3,z3},{r4,z4}}=ncoor;
For [k=1,k<=p*p,k++,
{qcoor,w}= QuadGaussRuleInfo[{p,numer},k];
{{N1,N2,N3,N4},dNr,dNz,Jdet}=Quad4IsoPShapeFunDer[ncoor,qcoor];
rk=r1*N1+r2*N2+r3*N3+r4*N4; c=Kfac*w*Jdet*rk;
brk=br1*N1+br2*N2+br3*N3+br4*N4;
bzk=bz1*N1+bz2*N2+bz3*N3+bz4*N4;
bk={N1*brk,N1*bzk,N2*brk,N2*bzk,N3*brk,N3*bzk,N4*brk,N4*bzk};
If [numer,bk=N[bk]]; fe+=c*bk;
]; Return[Simplify[fe]]
];

Figure 11.8. Module for computation of consistent node forces from a given body force field.

Kfac As in Quad4IsoPRingStiffness
bfor Body forces per unit volume. Specified as a two-item one-dimensional list:
{ br,bz },
or as a four-entry two-dimensional list:
{ br1,bz1 },{ br2,bz2 }, { br3,bz3 },{ br4,bz4 }.
In the first form the body force field is taken to be constant over the element.
The second form assumes body forces to vary over the element and specified
by values at the four corners, from which the field is interpolated bilinearly.

ClearAll[Em,Ν,a,b,d,h,p,num];
a=4; b=2; d=0; Kfac=2*Pi; Kfac=1; p=2;
ncoor={{d,0},{a+d,0},{a+d,b},{d,b}}; num=False;
fe=Quad4IsoPRingBodyForces[ncoor,Emat,Kfac,{num,p},9*{br,bz}];
Print["fe =",fe];
fe=Quad4IsoPRingBodyForces[ncoor,Emat,Kfac,{num,p},
9*{{a*br,bz},{(a+d)*br,bz},{(a+d)*br,bz},{d*br,bz}}];
Print["fe =",fe];

Figure 11.9. Test statements for body force module.

The module returns fe as an 8×1 one dimensional array arranged { fr1,fz1,fr2,fz2,fr3,fz3,

fr4,fz4 } to represent the vector
f(e) = [ fr 1 f z1 fr 2 f z2 fr 3 f z3 fr 4 f z4 ]T . (11.7)

11–10
11–11 §11.4 THE 4-NODE QUADRILATERAL RING ELEMENT

Quad4IsoPRingStresses[ncoor_,Emat_,Kfac_,options_,ue_]:=
Module[{k,m,numer=False,qcoor,w,c,N1,N2,N3,N4,
dNr1,dNr2,dNr3,dNr4,dNz1,dNz2,dNz3,dNz4,Jdet,
r1,r2,r3,r4,z1,z2,z3,z4,rk,Be,recovery="D",
sigp,qctab={{-1,-1},{1,-1},{1,1},{-1,1}},
TG={{1+Sqrt[3]/2,-1/2, 1-Sqrt[3]/2,-1/2},
{-1/2, 1+Sqrt[3]/2,-1/2,1-Sqrt[3]/2},
{ 1-Sqrt[3]/2,-1/2,1+Sqrt[3]/2,-1/2},
{-1/2, 1-Sqrt[3]/2,-1/2,1+Sqrt[3]/2}},
sige=Table[0,{4},{3}]},
m=Length[options]; If [m>0,numer=options[[1]]];
If [m>1,recovery= options[[2]]];
If [recovery=="G", qctab=qctab/Sqrt[3]];
If [numer,qctab=N[qctab];TG=N[TG]];
{{r1,z1},{r2,z2},{r3,z3},{r4,z4}}=ncoor;
For [k=1,k<=Length[sige],k++, qcoor=qctab[[k]];
If [numer, qcoor=N[qcoor]];
{{N1,N2,N3,N4},{dNr1,dNr2,dNr3,dNr4},{dNz1,dNz2,
dNz3,dNz4},Jdet}=Quad4IsoPShapeFunDer[ncoor,qcoor];
rk=r1*N1+r2*N2+r3*N3+r4*N4;
Be={{dNr1, 0,dNr2, 0,dNr3, 0,dNr4, 0},
{ 0,dNz1, 0,dNz2, 0,dNz3, 0,dNz4},
{N1/rk, 0,N2/rk, 0,N3/rk, 0,N4/rk, 0},
{dNz1,dNr1,dNz2,dNr2,dNz3,dNr3,dNz4,dNr4}};
If [numer,Be=N[Be]];
sige[[k]]=Simplify[Emat.(Be.ue)] ];
If [recovery=="G", sigp=sige; sige=TG.sigp];
Return[Simplify[sige]];
];

Figure 11.10. Module for recovery of element corner stresses from displacements.

Running the test statements shown in Figure 11.9 (the factor of 9 in the supplied forces are to get
exact integer output) gives

f(e) = [ 24br 24bz 48br 48bz 48br 48bz 24br 24bz ]T , (11.8)

f(e) = [ 80br 24bz 176br 48bz 160br 48bz 64br 24bz ]T . (11.9)

§11.4.6. Recovery of Element Corner Stresses

A stress computation module is shown in Figure 11.10. The arguments of the module are:
The arguments of the module are:
ncoor As in Quad4IsoPRingStiffness
Emat As in Quad4IsoPRingStiffness
Kfac Not used: a dummy argument kept as placeholder.
ue The 8 corner displacements components. These may be specified as a 8-entry
one-dimensional list form:
{ ur1,uz1, ur2,uz2, ur3,ur3, ur4,uz4 },

11–11
Chapter 11: ISO-P QUADRILATERAL RING ELEMENTS 11–12

or as a 4-entry two-dimensional list:

{ ur1,uz1 },{ ur2,uz2 },{ ur3,uz3 },{ ur4,uz4 }.
The module returns the corner stresses stored in a 4-entry, two-dimensional list:
{ { sigrr1,sigzz1,sigtt1,sigrz1 }, { sigrr2,sigzz2,sigtt2,sigrz2 }, { sigrr3,
sigzz3,sigtt3,sigrz3 }, { sigrr4,sigzz4,sigtt4,sigrz4 } } to represent the array of cor-
ner stresses, stored row-wise:
 
σrr 1 σzz1 σθ θ 1 σr z1
σ σzz2 σθ θ 2 σr z2 
σ (e) =  rr 2  (11.10)
σrr 3 σzz3 σθ θ 3 σr z3
σrr 4 σzz4 σθ θ 4 σr z4

ClearAll[Em,ν,a,b,d,p,num,err,ezz,grz,ur,uz];
Em=100; ν=0; a=5; b=2; err=1/10; ezz=-5/100; grz=4/100;
ur[r_,z_]:=err*r; uz[r_,z_]:=ezz*z+grz*r;
ncoor={{d,0},{a+d,0},{a+d,b},{d,b}}; num=False;
Emat=Em/((1+ν)*(1-2*ν))*{{1+ν,ν,ν,0},{ν,1+ν,ν,0},
{ν,ν,1+ν,0},{0,0,0,1/2-ν}};
Print["Emat=",Emat//MatrixForm];
ue={}; For [n=1,n<=4,n++,{rn,zn}=ncoor[[n]]; urn=ur[rn,zn];
uzn=uz[rn,zn];AppendTo[ue,urn]; AppendTo[ue,uzn] ];
Print["ue=",ue];
sige=Quad4IsoPRingStresses[ncoor,Emat,Kfac,{num,"D"},ue];
Print["corner stresses=",sige//MatrixForm];

Figure 11.11. Test statements for stress recovery module.

100

0 0 0 








Emat

0 100 0 0 


 



 0 0 100 0 


 
 0 0 0 50 
d d 5d 5d 5d 1 5d d 1 d
ue
  ,  ,  ,  ,  ,     ,  ,     
10 25 10 25 10 10 25 10 10 25

10 5 10 2 

 


 5 


corner stresses

10 10 2 


 


 10 5 10 2 


 

 10 5 10 2 

Figure 11.12. Result from running the statements of Figure 11.11.

The module is tested by the statements listed in Figure 11.11. The technique used is to generate
the element node displacements by evaluating a test displacement field

u r (r, z) = err r, u z (r, z) = ezz r + γr z z (11.11)

11–12
11–13 §11.4 THE 4-NODE QUADRILATERAL RING ELEMENT

in which {err , ezz , γr z } are specified strains, assumed constant over the element. Note that the hoop
strain is eθ θ = u r /r = err .
The field (11.11) is evaluated at the corner nodes to construct ue, which is then fed to the stress
recovery module. For the data err = eθ θ = 1/10, ezz = −5/100, γr z = 4/100, E = 100 and ν = 0
the recovered stresses must be

σrr = σθ θ = 10, σrr = −5, σr z = 2, (11.12)

at the four corners. The test is run taking a = 5, b = 2 and leaving d as variable. The results,
shown in Figure 11.12, verify the stress values for both recovery methods, although only results for
recovery="D" are shown.

11–13
Chapter 11: ISO-P QUADRILATERAL RING ELEMENTS 11–14

§11.5. THE 9-NODE QUADRILATERAL RING ELEMENT

To be included.

§11.6. THE 16-NODE QUADRILATERAL RING ELEMENT

To be included.

11–14
11–15 Exercises

Homework exercises for Chatre 11

Iso-P Quadrilateral Ring Elements

EXERCISE 11.1

11–15