NUMERICAL ANALYSIS OF NONLINEAR PlVEuMATIC STRUCTURES

J. T. Oden
Associate Professor of Engineering Mechanics

W. K. Kubitza
Professor of Engineering Mechanics
University of Alabama Research Institute
University of Alabama in Huntsville
Huntsville, Alabama
U. S. A.

GPO PRICE $

CFSTI PRICE(S1 $

Hard copy (HC)

Microfiche (MF)
ff 653 July 65

COLLOQUIUM ON €'NEUMATIC STRUCTURES

Stuttgart, Germany
1967
 EJUMERICAL AXALYSIS OF N O N L I W PNRlMLlTIC STRUCTURES

J. T. Ode$ and W. K. Kubitza*-x

Abstract. This paper presents a systematic numerical procedure for t h e

analysis of nonlinear behavior i n general pneumatic structures. Recent
advances i n t h e application of t h e f i n i t e element method t o t h e evaluation
of f i n i t e s t r a i n s and l a r g e displacements of e l a s t i c membranes a r e reviewed
and extensions of t h e method t o t h e analysis of l a r g e motions of reinforced
f a b r i c s , anisotropic metals, p l a s t i c s , v i s c o e l a s t i c , and nonlinearly elastic
materials a r e presented. Locai yielding of m e t a l l i c e l a s t o - p l a s t i c m a -
branes subjected t o e x t e r n a l pressure i s a l s o examined. By using l i n e a r
displacement approxbations and triangular f i n i t e elements, general non-
l i n e a r s t i f f n e s s r e l a t i o n s are derived. These lead t o systems of nonlinear
algebraic or ordinary d i f f e r e n t i a l equations i n t h e generalized displacements
and v e l o c i t i e s which are solved numerically. Numerical restilts are included
along with comparisons with available experinental data,

CONTENTS

NOTATION

1. INTRODUCTION
1.1 Opening Remarks.
1.2 Previous ,Related Work.
1.3 Scope.
2. GEOMETRIC AND KINEMATIC CONSIDEFATIONS
2.1 The d i s c r e t e Model.
2.2 S t r a i n s and Displacements.

. 3. THERMODYNAMICS OF FINITE ELEMENTS

3.1 The general Equation of Motion of a F i n i t e Element.
3.2 I n t e r n a l Energy.

4. NONLINEAR STIFFNESS RELATIONS

4.1 E l a s t i c MaAerials.
4.1.1 Natural and Synthetic Rubbers.
4.1.2 P l a s t i c s and Nonlinearly E l a s t i c Materials
4.1.3 Metals and Reinforced Fabrics.
4.2 Viscoelastic Makerials.
4.3 Elasto-Plastic Materials.

Associate Professor of Engineering Mechanics, University of A l a b a m a

++
Research I n s t i t u t e , Huntsville, Alabama, U. S. A.

++*Professor of Engineering Mechanics, University of Alabama Research

I n s t i t u t e , Huntsville, A l a b a m a , U..S, A.
 2

5. FORMSJLATION OF THE STRUCTU%L PROBIXM

5 . 1 Global Equations of Motion.
5.2 Ekternal Pressure.
5.3 Solution of Nonlinear Equations.
6. NUMERICAL AND EXR3RIMENTCIL RESULTS
6.1 S t r e s s Diffussion i n a Stiffened Panel.
6.2 Elasto-plastic Behavior of a Metallic Membrane.
6.3 F i n i t e Stretching of a Rubber Sheet.
6.4 I n f l a t i o n of an I n i t i a l l y F l a t Rubber Membrane.
 3

NOTATION

Indicia1 notation and t h e summation convention a r e used throughout

t h i s paper. Upper-case Latin indices indicate points i n space and lower-

case indices i n d i c a t e elements of an array. I n general, Greek indices are

associated with l o c a l coordinate systems and range from 1 t o 2. The

following symbols a r e used:

Constants i n displacement approximation

Area of deformed and undeformed element

Thermal load vector a t node M

Node displacement coefficients

Material constants

Components of rigid-body t r a n s l a t i o n

Deformation r a t e tensor

Element i d e n t i f i c a t i o n index

Total number of f i n i t e elements

E l a s t i c moduli

Multi-dimensional a r r a y of material constants

Yield surface

Body force per u n i t mass

A surface tensor

Shear moduli

H Heat input per u n i t mass

S t r a i n invariants

Consistent mass matrices

n Total number of nodes

Hydrostatic pressure

Generalized node forces of element e i n l o c a l coordinates

 4

Generalized node forces i n global coordinates

Ekternal pressure

Components of heat flux

Heat input

Element node forces due to q

Displacement components i n l o c a l coordinates

ui
Displacement of node N of element e i n xi d i r e c t i o n i n
l o c a l coordinates

U Total i n t e r n a l energy

'Ni
Displacement of node N i n Zi-direction i n global coodina,tes

v,v0 Volumes of deformed and undeformed elements

W S t r a i n energy per unit of undeformed volume

Xie Local coordinates of element e

X Local coordina,tes of node N of element e

Nie
Local coordina,tes of deformed element e
'ie
Local coordinates of node N of dement e a,fter deforma,tion
'Nie
'i Global coordinates

'Ni Global coordinates of node N

Orthogonal transformation ma,trix of element e

'ije
Lagrangian s t r a i n tensor
'i j

6ij Kronecker d e l t a
€ Two-dimensional permutation symbols

H Kinetic energy

x Ext ens ion rati o

Poisson's r a t i o s

I n t e r n a l energy per unit mass

Mass d e n s i t i e s
 5

0. S t r e s s tensor
a@
-0 Deviatoric s t r e s s 'censor
a@
R Power of e x t e r n a l forces

R Multi-dimensional array
iWe
 6

1 INTRODUCTION

1.1 Opening Remarks. U n t i l recent years, t h e behavior of the majority of

p r a c t i c a l s t r u c t u r e s could be adequately described by l i n e a r theory. The

deformations of most s t r u c t u r a l systems under working loads a r e usually

so small as t o be scarcely detectable with t h e unaided eye, and t h e s t r e s s -

s t r a i n r e l a t i o n s f o r such common materials as s t e e l , aluminum, and even

concrete can, f o r p r a c t i c a l purposes, be treated as l i n e a r . Solutions of

l i n e a r problems involving two and three-dimensional s t r u c t u r e s of general

shape w i t h complex boundary conditions, however, a r e often untractable by

c l a s s i c a l means and, even with t h e gross simplifications afforded by l i n e a r

t h e o r i e s , many important problems remain unsolved.

With t h e bulk of t h e available methods of analysis being applicable

t o only l i n e a r systems and with these methods being often inadequate i n

the face of complex geometries, t h e engineer must look upon the recent; .trend

toward t h e use of highly f l e x i b l e pneumatic s t r u c t u r e s w i t h some bewilderment.

The behavior of i n f l a t a b l e pneumatic s t r u c t u r e s i s inherently nonlinear:

such s t r u c t u r e s o f t e n acquire t h e i r primary load-carrying capacity a f t e r

undergoing deformations which, even under small pressures, may be so large

t h a t t h e o r i g i n a l undeformed shape i s unrecognizable. S t r a i n s appreciably

greater than u n i t y are not uncommon, and i n such cases Hooke's l a w i s not

applicable. Moreover, t h e materials used t o construct pneumatic s t r u c t u r e s

a r e o f t e n anisotropic and nonlinearly e l a s t i c and, t o further complicate

matters, t h e d i r e c t i o n s and magnitudes of t h e applied loads change w i t h t h e

deformation. To emphasize t h i s point, one need only r e f e r t o recent

experimental s t u d i e s on pneumatic s t r u c t u r e s wherein discrepancies on t h e

order of 400 per cent were encomtered when measured s t r e s s e s were compared

with those predicted by conventional linear theories.

 7

It i s c l e a r t h a t i n order t o accurately analyze general pnewmtic

s t r u c t u r e s one must t u r n t o the general nonlinear theories of s t r u c t u r a l

mechanics, many of which have long been regarded as only academic i n t e r e s t .

Although only a small number of exact solutions t o general nonlinear s t r u c t u r a l

problems a r e a v a i l a b l e i n t h e l i t e r a t u r e and although these, without excep-

t i o n , a r e concerned with only the most simple geometry, deformation, and

loading, a t l e a s t one can f i n d h e r & a rigorous and complete foundation on

on which t o base nonlinear analyses.

Fortunately, t h e trend toward the use of nonlinear s t r u c t u r a l systems

has been accompanied by s i g n i f i c a n t developments i n both large-scale d i g i t a l

computers and general methods of numerical analysis. W i t h t h e a i d of these

t o o l s , many of t h e nonlineax s t r u c t u r a l theories can be employed t o obtain

u s e f u l information concerning t h e nonlinear behavior of pneumatic s t r u c t u r e s b

Chief among the numerical schemes i s t h e so-called f i n i t e element method

wherein continuous s t r u c t u r a l systems a r e replaced by d i s c r e t e models whose

properties a r e consistent w i t h the general f i e l d equations defining t h e

behavior of t h e continuum. Notable progress has recently been made i n

applying t h i s method t o complex l i n e a r and nonlinear s t r u c t u r a l problems.

A comprehensive review of applications of the f i n i t e element method t o t h e

analysis of nonlinear behavior i n e l a s t i c membranes as w e l l as several

extensions of t h e method t o t h e analysis of l a r g e deformations of e l a s t i c ,

e l a s t o - p l a s t i c , and v i s c o e l a s t i c pneumatic s t r u c t u r e s i s t h e subject of

t h i s paper.

1.2 Previous Related Work. Since 1950, Lhe technical l i t e r a t u r e has contained

numerous applications of the f i n i t e element approach t o a wide v a r i e t y of

l i n e a r s t r u c t u r a l problems. Solutions t o complex plane-stress, p l a t e , and

s h e l l problems a r e available, along with sever solutions t o three-ditnens

e l a s t i c bodies. Applications t o nonlinear problems, however, have been

s i g n i f i c a n t l y l e s s extensive, It appears t h a t the f i r s t successful appli-

cation of t h e f i n i t e element concept t o the analysis of geometrically

nonlinear problems was presented b& Turner e t a1 c11. These authors

solved c e r t a i n nonlinear problems by dividing a large deformation i n t o a

number of steps. Within each s t e p t h e s t r u c t u r e i s assumed t o behave

l i n e a r l y and an instanteous ( l i n e a r ) s t i f f n e s s matrix i s computed i n the

deformed geometry. Argyris k,31, Gallagher and Padlog Csl, and Martin C5I

were among several investigators who l a t e r applied t h i s successive correction

technique t o large d e f l e c t i o n and s t a b i l i t y analyses. I n these papers, the

correction t o t h e l i n e a r s t i f f n e s s matrix i s often referred t o as the

11
geometric" s t i f f n e s s of the s t r u c t u r e , and it i s o r d i n a r i l y a function

of t h e s t r e s s e s associated w i t h some reference equilibrium s t a t e . A survey

of the liteera,ture using t h i s approach i s contained i n the paper by Martin

c'jl and a general formula f o r geometric s t i f f n e s s matrices was presented by

Oden C61. Following a d i f f e r e n t approach, Wissmann C7,8J obtained nonlinear

f i n i t e element forumlations f o r c e r t a i n problems involving l a r g e displace-

ments but small s t r a i n s of e l a s t i c structures.

&tensions of t h e f i n i t e element method t o t h e analysis of f i n i t e

deformations of e l a s t i c membranes and three-dimensional bodies were presented

by Oden C9,lOl and Oden and Sat0 [111. Using a somewhat d i f f e r e n t approach,

Becker h 2 1 obtained a numerical solution t o t h e problem of f i n i t e in-plane

deformations of rubber sheets subjected t o prescribed boundary displacements.

These nonlinear formulations contain t h e l i n e a r s t i f f n e s s matrices as

s p e c i a l cases and lead t o systems of nonlinear algebraic equations which

must be solved numerically. A t t h e same time, they a r e considerably more

general than t h e successive correction methods mentioned e a r l i e r i n that

 9

they contain c h a r a c t e r i s t i c s which are encountered only i n highly nonlinear

st r uctux a 1 behavior .
1.3 Scope. In t h e discussion t o follow, t h e basic philosophy of t h e f i n i t e

element representation of f l e x i b l e pneumatic structures i s presented.

General kinematic properties of t h i n membranes a r e c a s t i n t h e form of

algebraic equations defining t h e motion of an assembly of f i n i t e elements.

The f i r s t law of thermodynamics i s called upon t o provide general r e l a t i o n -

ships between kinematic and k i n e t i c variables associated with the behavior

of f i n i t e elements of a r b i t r a r y pneumatic structures. This leads t o t h e

general equation of motion of f i n i t e elements of t h i n membranes, and includes

such properties a s anisotropy, nonlinear v i s c o e l a s t i c i t y , thermoviscoelasticity,

nonhomogenuity, and p l a s t i c i t y with no r e s t r i c t i o n s on t h e magnitudes of the

deformations CEq. (20)I. In order t o obtain quantitative r e s u l t s , t h e

general formulation i s modified so t h a t it applies t o a number of important

s p e c i a l cases. These include a review of the aaalysis of f i n i t e deformations

of e l a s t i c membranes given i n Refs. 9, 10, and 11 and applications t o elasto-

p l a s t i c and v i s c o e l a s t i c structures. T h i s is followed by discussions of

procedures for assembling t h e f i n i t e elements, computing changes i n loading

due t o deformation, and solution o f t h e nonlinear equations generated i n t h e

analysis. Finally, t h e solutions of several representative problems a r e

presented and numerical values a r e compared with available experimental

data.

2. GEOMETRIC AND KINEMATIC COWSIDERATIONS

2.1 The Discrete Model. Classically, t h e analysis of continuous systems

begins with investigations of t h e properties of small d i f f e r e n t i a l elements

of the continuum under investigation. Relationships a r e established between

FIG. 1 F i n i t e element representation of a pneumatic s t r u c t w e .
 10

mean values of various q u a n t i t i e s associated with t h e i n f i n i t e s i m a l

elements and p a r t i a l d i f f e r e n t i a l equations governing t h e behavior of

t h e e n t i r e domain a r e obtained by allowing t h e dimensions of t h e elements t o

approach zero as t h e number of elements becomes i n f i n i t e l y large. In

contrast t o t h i s c l a s s i c a l approach, i n t h e present study we begin with

invcc:kigations of t h e properties of elements of f i n i t e dimensions. We

may employ equations of t h e continuous system i n order t o a r r i v e a t t h e

properties of these elements; but the dimensions of t h e elements remain

f i n i t e i n t h e analysis, integrations a r e replaced by f i n i t e summations, and

t h e d i f f e r e n t i a l equations of t h e continuous s t r u c t u r e a r e replaced by

systems of algebraic or ordinary d i f f e r e n t i a l equations. The continuous

system with i n f i n i t e l y many degrees of freedom i s thus represented by a

d i s c r e t e model which has f i n i t e ddgrees of freedom. Moreover, i f c e r t a i n

kinematic conditions a r e s a t i s f i e d , then, as t h e number of f i n i t e elements

i s increased and t h e i r dimensions a r e decreased, t h e behavior of t h e d i s c r e t e

system converges monotonically t o t h a t of t h e continuous system.

Consider, f o r example, t h e general pneumatic s t r u c t u r e shown i n Fig. l a

subjected t o a general system of applied loads. To define t h e i n i t i a l

geometry of t h i s system, a fixed rectangular Cartesian coordinate system

Z Z Z i s established which i s referred t o a s t h e global reference frame.

1, 2’ 3
I n general, an i n f i n i t e number of coordinates Zi are required t o completely

specify the i n i t i a l configuration of the membrane; but i n t h e present

.inalysis, we reduce t h i s continuously d i s t r i b u t e d system t o a d i s c r e t e one

by representing t h e s t r u c t u r e as an assembly of a f i n i t e number Ee of f l a t

triangular elements, as indicated i n Fig. l b . The v e r t i c e s of these

triangular a r e r e f e r r e d t o as t h e node points of t h e d i s c r e t e model. Thus,

if n denotes t h e t o t a l number of nodes i n t h e system and if ZNi(N=1,2,.r,n;

 22

/*3

FIG.2 Local and global coordinates of a typical f i n i t e element,

 11

i = l , 2 , 3 ) denote t h e global coordinates of a t y p i c a l node N, then the set

of numbers ZNi define the geometry of the d i s c r e t e system,

I n the f i n i t e element method, it i s convenient t o f i r s t describe the

behavior of each element independefitly i n terms of the displacements of i t s

nodes; t h e e n t i r e set of elements i s then connected together by establishing

c e r t a i n dependencies between appropriate node displacements. Toward t h i s

end, fixed l o c a l coordinate systems xie(i=1,2,3; e=l,2,. . ,E e ) a r e established

i n t h e neighborhoods of each f i n i t e element. "he q u a n t i t i e s xie a r e referred

t o as t h e l o c a l coordinates of element e. For simplicity, it i s assumed

t h a t t h e middle surface of each element e l i e s i n t h e ~ ~ - plane

~ of i, t s x

l o c a l coordinate system. Following a procedure similar t o t h a t used f o r

global coordinates, t h e l o c a l coordinates of a t y p i c a l node N of element e

are denoted xNie(N,i = 1,2,3; e=1,2,. ..,E e ) where, due %o our p a r t i c u l a r

choice of coordinates, x = 0. These coordinates are i l l u s t r a t e d i n Fig. 2,

We
A r i g i d r o t a t i o n of a t y p i c a l l o c a l system xie i n t o a l o c a l system
ie
whose coordinate l i n e s are p a r a l l e l t o the corresponding global coordinate

axes, i s accomplished by t h e orthogonal transformation

-x = 6.. x (no sum on e )

ie J X ~ je
where 6 , . i s the d i r e c t i o n cosine of t h e r7,ngle between Z (or x ) and x
-
Jle i ie je
The l o c a l coordinates of node N p a r a l l e l t o t h e global coordinates are

denoted %ie.
2.1 S t r a i n s and Displacements. I n order that the present discussion be self-

contained, we reproduce i n t h i s section the ' derivation of t h e general nonlinear

s t i f f n e s s r e l a t i o n s f o r f i n i t e e l e m n t s of membrane s t r u c t u r e s [g,lO,llj For

t h e present, we confine OUT a t t e n t i o n t o a ty-pical f i n i t e element and we

temporarily drop t h e element i d e n t i f i c a t i o n index e f o r c l a r i t y .

 12

Under a general deformation, l i n e elements o r i g i n a l l y straight i n t h e

undeformed structure become curved l i n e s i n the deformed structure. Hence,

an i n i t i a l l y f l a t element i s generally deformed i n t o a curved surface.

However, i f the node p i n t s a r e selected s u f f i c i e n t l y c l o s e of one ano-ther,

node l i n e s i n the deformed configuration are closely approximated by s t r a i g h t -

l i n e segments. Then plane elements remain plane a f t e r t h e deformation,

This i s equivalent t o assuming that t h e displacement f i e l d s within each

element a r e l i n e a r functions of t h e l o c a l coordinates of t h e element,

Assuming, for simplicity, that t h e element i s i n i t i a l l y i n t h e xL, x

2
plane and denoting by ui t h e components of displacement r e f e r r e d t o the

local coordinates of t h e element under consideration, it follows that

ui = di f aiaxa i = 1, 2, 3 ct = 1, 2 (2)

where d i a r e the rigid-body t r a n s l a t i o n s of t h e element and t h e aia a r e

undetermined constants. By evaluaking Eq. (2) a t each of t h e three nodes

of t h e element, we a r r i v e a t nine simultaneous equations i n the nine un-

knowns di, aia. Solving these we find t h a t

d. = k u
1 N Ni

and
a = c (4)
ia dN'Ni

where \i i s t h e displacement of node N i n t h e xi d i r e c t i o n ,

a
kl = -(
1
. x21x32 - X31x22)
k 2 = l ( x x - x x )
2Ao E 31 11 32

k = L ( x x - X x )
3 a0 11 22 12 31

and
 I n these equationq, xNa (N = 1, 2, 3; a = 1, 2) a r e t h e l o c a l coordinates

of node N and A. 5,s t h e area of t h e undeformed t r i a n g l e .

Substituting Eqs. (3) and (4) i n t o Eq. (2) gives

According t o Green and Adkins [131, t h e Lagrangian strain tensor for

a t h i n membrane i s given by

Ya3 = 0

- 1(h2 - 1)
y33 - 2

where h i s a s c a l a r function representing t h e exbension r a t i o a t t h e middle

. surface i n a d i r e c t i o n normal t o t h e surface. For very t h i n membranes, t h e

strains a r e e s s e n t i a l l y uniform over t h e thickness and

x = -h (9)
ho
where ho and h a r e t h e thicknesses of t h e membrane before and a f t e r deformation

respectively.

Introducing Eq. (7) i n t o Eqs. ( 8 ) , we f i n d for the straihs i n t h e d i s c r e t e

Ya3 = 0 y33 = $12 - 1)

a,@ = 1,2 M , N , i = 1,2,3
The s t r a i n components a r e thus constant throughout the f i n i t e element and,

since they a r e determined from prescribed displacement f i e l d s , they auto-

matically s a t i s f y t h e equations of compatibility throughout the element.

Note a l s o that t h e components of displacement along the boundaries of an

element a r e l i n e a r functions of t h e local coordinates, Since each edge

contains two nodes and since t h e displacements are l i n e a r along each edge,

it follows t h a t the displacements are continuous across element boundariesc

These properties of t h e approximation instwe monotonic convergence of the

solutions as t h e finite-element representation i s refined.

Three invariants can be formed from every symmetric second-order tensor.

I n t h e analysis of deformable membranes, it i s convenient t o form t h e

invariants of the deformation tensor (hij + 2Y. .), whew!

=J
is t h e fionecker
delta ('ij = 1 f o r i = j and hij t= 0 for i # j). These invariants are given '
by t h e formulas
I1
= X2 + 2(1 )
+,Y

I2
= - x4 + X2 I
1
1
+-I
A2 3

where

and cap, ExtJ. are the two-dimensional permutation symbols (El2 = 1, E21 = -1,

E:11 = EZ2 = 0).

3. THERMODYNAMICS OF FINITE EIXmmS

3.1 The General Equation of Motion of a F i n i t e Element. Having defined the

d i s c r e t e system i n t h e previous section, we now t u r n t o c e r t a i n ftmdamental

principles of mechanics i n order t o obtain general equations of motion f o r

a f i n i t e element. We begin with t h e f i r s t l a w of thermodynamics:

i+ir=n+g (J-3)
wlicre X i s t h e t i m e r a t e of change of k i n e t i c energy, U i s the time rate of

change of i n t e r n a l energy, CJ i s t h e power developed by t h e external forces,

and Q i s t h e heat input. By d e f i n i t i o n ,

U = I P S dv
V

= ‘F.G.Pdv + SiGi ds
J 1 1
V S

v
PH dv +
s
S
gini ds

I n these equations, P i s t h e mass density, dv i s t h e d i f f e r e n t i a l volume,

Gi a r e t h e v e l o c i t y components, 5 i s the i n t e r n a l energy per u n i t mass, Fi

a r e the body forces per u n i t mass, S i a r e t h e surface t r a c t i o n s per u n i t of

surface area s, H i s t h e heat input per u n i t mass, qi ape t h e components of

heat flux, and ni are t h e components of a u n i t vector normal t o the boundary

sa-faces of the elment.

Through arguments similar t o those used i n obtaining Eq. (71, it i s

found t h a t

and

where uNi i s t h e v e l o c i t y of node N i n d i r e c t i o n i and q i s t h e heat flux

Ni
a t node N i n d i r e c t i o n i, Moreover, i f Po, vo and P, v denote the mass
d e n s i t i e s and t h e volumes of t h e f i n i t e element i n the undeformed and t h e

deformed s t a t e s respectively, then according t o the p r b c i p l e of conservation

of mass,
 16

Povo = Pv (17)

Introducing Eqs. (15), (;t6),and (17) i n t o Eqs. (14) and simplifying,

we f i n d

where

The a r r a y i s the so-called consistent mass matrix for t h e f i n i t e element

P

and the q u a n t i t i e s pNi a r e t h e components of generalized force corresponding

'
to t h e generalized displacements uNi. Physically, pNi i s the generalized

force a t node N i n d i r e c t i o n i. The q u a n t i t i e s bNi are the generalized

thermal gradients usually r e f e r r e d t o as thermal loads. Note that the

heat input H i s regarded as a constant for each f i n i t e element.

Introducing Eqs. (18) i n t o Eq. (13), we obtain t h e general result

where %i
.. are t h e components of acceleration of the nodes. This r e l a t i o n

represents the most general d i s c r e t e representation of the equations of

motion f o r a f i n i t e element of a continuous media, It is applicable t o the

analysis of large deformations of nonlinearly e l a s t i c , p l a s t i c , v i s c o e l a s t i c ,

and thermoviscoelastic media since no r e s t r i c t i o n s have as yet been placed

QD the c ~ n s t i t ~ t i vequations
e r pp3;iteriaZ of.&&
f ~ the .Ithe element is

composed or t h e magnitude of t h e deformation, The second and t h i r d terms

on t h e right-hand s i d e of Eq. (20) represcnt thermal e f f e c t s on t h e motion

of t h e element.
.
The r a t e of change of i n t e r n a l energy 5 i s , i n general,

a l s o a function of t h e heat flux and. the temperature gradients. These e f f e c t s

will not be considered fwrther i n thls discussion; for our purposes, it

s u f f i c e s t o merely point out t h a t thermal e f f e c t s can be e a s i l y included

i n t h e analysis by r e t a i n i n g the terms bNiqNi and PovoH i n Eq, (20) and by

introducing, i n addition t o a c o n s t i t u t i v e equation involving the s t r e s s ,

a second c o n s t i t u t i v e equations which r e l a t e s heat f l u x to thermal gradients ,

deformation rates, s t r a i n s , e t c . Thus, w i t h thermal e f f e c t s omitted,

Eq. (20) reduces t o

3.2 I n t e r n a l Energx.
To apply equation to s p e c i f i c kmterials, it i s
.
necessary to obtain < a s a function of
this

generalized displacements and

the

t h e i r time derivatives. According t o Firingen [141, i f couple s t r e s s e s and

thermal e f f e c t s a r e omitted,

where oij i s t h e s t r e s s tensor and dij i s the deformation r a t e tensor, For

the f i n i t e element representation, it i s e a s i l y shown t h a t

0..d
13 i j
= 0 c
ij iN
(6k j + c j M k ) C m

We now assume t h a t t h e material i s homogeneous so that f i s no% a

 18

function of t h e l o c a l coordinates x . . I n view of Eq. (22), it i s a constant

1
w i t h respect t o x for the f i n i t e element and can therefore be factored
i
outside of t h e Q t e g r a l i n Eq. (21). Noting a l s o t h a t

P0 = J i 5 P

where 1 i s the t h i r d invariant i n Eq. (12), and introducing Eqs. (22), ( 2 3 ) ,

3
and (24) i n t o Eq. (21), we a r r i v e a t t h e equation

Since t h i s equation mst hold f o r a l l v e l o c i t i e s %i, we have

It i s now necessmy t o introduce i n t o t h i s equation t h e appropriate

relakion expressing t h e stress i n terms of t h e s t r a i n , s t r a i n r a t e s , higher-

order s t r a i n r a t e s , etc. i n order t o obtain the f i n i t e element representation

for s p e c i f i c materials.

4. NONLINEAR STIFFNESS RELATIOEJS

We now examine applications of Eqs. (21) and (26) t o various types of

mmbranes .
4.1 E l a s t i c Materials. I n t h e case of e l a s t i c materials, t h e s t r e s s a i j i s

derivable from a p o t e n t i a l M c t i o n W which represents the strain energy per

u n i t of undeformed volume. If each element i s homogeneous, then 5 and W a r e

not f'unctions of t h e s p a t i a l coordinates xi' It follows t h a t

spot dvo = G dvo = vo;

V V
where

t
Introducing these equations i n t o Eq. (21) and simplifying t h e results, we
obtain

This r e s u l t must be v a l i d f o r a r b i t r a r y node v e l o c i t i e s of t h e element;,

Therefore

Equation (30) i s t h e general equation of motion f o r f i n i t e forced oscillaticms

of an e l a s t i c membrane. For s p e c i f i c materials, t h e appropria'ke form of W

must be introduced i n t o t h i s equation,

The important problem of small o s c i l l a t i o n s about a s t a t e of l a r g e

deformation i s obtained as a s p e c i a l case of Eq. (30) by denoting

% i = %fi ?iN i (31)

where
-%i i s t h e prescribed l a r g e displacement a n d w is a small perturbation.
Ni
Then

%yfi + fu(cMi)wm= PNi

- (32)

where

Here f is known a known function of

-%i and i s independent of time. "he
NM
forces 5N i are known functions of t i m e and EqA (32) is l i n e a r i n t h e dependent

variables w and t h e i r derivatives.

Ni
If t h e membrane i s i n equilibriwn,Eq. (30) reduces t o [10,111

P (34)
Ni

For e l a s t i c membranes, t h e s t r e s s e s a r e calculated by means of the

formulas
 20

4..1.1 Natural and Synthetic Rubbers. I n t h e case of highly e l a s t i c

materials such as n a t u r a l and synthetic rubbers, several forms of t h e

s t r a i n energy function W a r e available, Ordinarily, such materials a r e

assumed t o be incompressible and, for i s o t r o p i c membranes, t h i s incom-

p r e s s i b i l i t y condition i s equivalent t o t h e condition

I = l
3

Then, W i s a function of only I1 and I and, instead of Eq, (35), t h e s t r e s s e s

2
a r e given by E131

q-J3 = 0 (37)

where h2 i s defined i n Eq. ( g ) , p i s the hydrostatic pressure, a n d

gap = 'a,@ 4- 2 E ~ o E i ~ ~ ~ ~ (38)

It i s assumed t h a t t h e membrane i s very t h i n so t h a t t h e s t r a i n s a r e

e s s e n t i a l l y uniform over the thickness a.nd 0 t h e s t r e s s normal t o t h e

33 '
deformed surface, i s negligible i n compa.rison w i t h (J Then p can be
&E3 e

determined from the condition (5 = 0:

33

I L

Thus /

Rivlin and S a n d e r s [15,161 v e r i f i e d experimentally t h a t t h e s t r a i n

energy function for m c t s t i s o t r o p i c , incompressible rubbers 5s of t h e form

w = C1(Il - 3) 4- $(I2- 3) (41)

where C1 i s a material constant and t h e function $ depends upon t h e ty-pe of

rubber. 1;: a n a l y t i c a l work, t h e most camon form of W i s t h e well-known

 21

Mooney form C17I:

where C and C are experimentally determined constants. Treloar C18I, using

1 2
a s t a f t i s t i c a l approach based on molecular theory, found for incompressible

rubbers

w = c(11 - 3) (43)

where C i s a constant. Rivlin h 9 1 refers t o such materials as neo-Hookean,

To obtain nonlinear s t i f f n e s s r e l a t i o n s f o r f i n i t e elements of rubber

membranes, f i r s t note t h a t for incompressi2de membranes

2
I1 = h 4- 2(1 -I- y010,)

and

(46)
and

gas = + EwE)$(Cp~u~h + hNI’N p c~JvflJ$-Niup.li) (47)

whereini, j , L , M , N , K = l , 2, 3 ; a , @ , h , ! J . = l , 2.
Introducing Eq. (42) and Eqs. (44)-(47) i n t o Eqs. (34) and ( b o ) , we

obtain t h e following equations for a f i n i b e element of a membrane of Mooney

 22

material C 113,

where

Similarly, for neo-Hookean membranes we find

governing the s t a t i c behavior of f i n i t e e1emen:ts of Mooney and neo-Hookean

membranes. Upon connecting t h e elements appropriately and applying boundary ’

conditions, these lead t o systems of nonlinear algebraic equations i n the

iiuae dicplzcemants . Once solved, the r e s u l t s are introduced i n t o Eqs (LO),

(48b), and (5Ob) t o obtain f i n a l s t r a i n s and s t r e s s e s i n t h e struc-bwce.

4.1.2 P l a s t i c s and Nonlinearly E l a s t i c Materials. Recent experiments on

p l a s t i c s and synthetic materials have attempted t o a r r i v e at; approximate

energy functions for such materials even through t h e i r behavior i s not always

perfectly elastic. These have led t o highly nonlinear forms for the

p o t e n t i a l function W. For example, experiments on a dimethyl siloxane

rubber by Hutchinson, Becker, and Landel E201 have indicated an energy

where B1, B2, B B4, bl,and b2 aye makerial constants.

3’
 23

I n t h i s case, t h e nonlinear s t i f f n e s s r e l a t i o n i s given by

Nk
- 3)1(S,p - 4
“ae>

- 4 4 2
2h - 21 Yw) f
baPll
a s before, s t r e s s e s a r e computed by means of Eq. (b).

14.1.3 Metals and Reinforced Fabrics. Deformations of metallic and reinforced

f a b r i c pneumatic s t r u c t u r e s a r e usually characterized by large displacements

accompanied by s t r a i n s which are small i n comparison with unity. In such

cases, the material i s assumed t o be homogeneous and e l a s t i c within each

f i n i t e element and t h e well-known Wookean form of t h e s t r a i n energy f’unction

i s applicable. Moreover, t h e s t r a i n normal t m the deformed surface can then *

be expressed as a l i n e a r f’unction of thc s t r a i n s i n the middle surface and

t h e s t a t e - o f - s t r a i n i s described by t h e two-dimansional tensor Yap. The

s t r a i n energy function i s therefore of t h e form

where E@hV i s a multi-dimensional a r r a y of c l a s t i c constants. The stress

tensor i s given by
I
‘a@ = Ea@&yXp= % ~ ~ ~ “ x N ( ‘ $ & l 5 “@M1-k%i 1 (544

and t h e nonlinear s t i f f n e s s r e l a t i o n €or the f i n i t e element i s obtained by

introducing Eq. (53) i n t o Eq. (34):
 24

and, f o r i s o t r o p i c materials, it i s given by

where E i s Young's modulus and V i s Poisson's r a t i o .

The form of Eapw i n Eq. (56) i s applicable t o most engineering metals

( s t e e l , aluminum, and tungsten, e t c . ) and t o i s o t r o p i c f a b r i c s . However,

most of t h e reinforced f a b r i c and composite materials are e i t h e r completely

aniso4xopic or transversely i s o t r o p i c only with respect t o a normal t o the

middle surface a

A t y p i c a l example of such a reinforced f a b r i c i s indicated i n Fig. 3.

Here a f a b r i c core material, which i s assumed t o be i s o t r o p i c and t o have

a modulus E and Poisson's rate V, i s reinforced by a network of fibers of

modulus and Poisson's r a t i o 5. The fibers form an angle a with respect

t o t h e yl-axis, as shown, When t h e f i b e r spacing i s r e l a t i v e l y small, it

i s convenient t o introduce mean values for t h e a3propriate material

constants of t h e composite s t r u c t u r e . Thus, we introduce t h e following

constants h.1:
El = kE + (1-k)E
EZ
E =
2 kE + (1-112
GE (57)
G =
C kG + (1-k)G

where k i s t h e r a t i o of the cross-sectional a r e a of t h e fibers t o t h e t o t a l

cross-sectional area of the element, G and 5 are t h e respective shear moduli

of the core fabric and t h e f i b e r s , El and E2 are respectively the effective

mean moduli of e l a s t i c i t y i n t h e y1 and y2 d i r e c t i o n s , Gc i s the mean shear

modulus, and Vc i s t h e mean Poisson's r a t i o of the composite material. We

FIG. 3 Fiber reinforcement p a t t e r n i n composite fabric sheet.
then have
a cos4a + aZ2sin4a + 2(a12 + 2GC)sin% cos2a
E1lll=11

E1122= E2211
-
= (au + a22 bGc)sin2a cos 2U + a (sin4a
12
+ cosk)
E2222= allsin 4a f a 2 2 a~ +~2(a12
~ + 2Gc)sin2a cos2a
4 (58)
- - (aLE + a - 2aE)sin2a cos%
Ek7m= Ea221 - = 5323. 22
+ ~,(cos% - s i n k ) 2
- - - - =o
E1112= E1121 - E2212 - E2221 - E1211- E1222

where

- El .- (59)
a22 - 1 - V 2 E 'E,
c l ,

4.2 Viscoelastic Materials. The equati'ons of motion of a f i n i t e element of

a v i s c o e l a s t i c material a r e obtained by introducing the appropriate c o n s t i t u t i v e

equation, w r i t t e n e x p l i c i t l y i n terms of the stress tensor aij, i n t o Eq. (26),

We confine our a t t e n t i o n t o Lsotropic materials f o r which the stress i s

dependent on the displacement and v e l o c i t y gradients. I n general, t h e

c o n s t i t u t i v e equations for such materials can be w r i t t e n ' i n the form

T= -02

w h e r e a i s a tensor functional of t h e indfcated

of motion for a f i n i t e element are

 26

Specific forms of these r e l a t i o n s for given materin.1.s a r e obtained by

introducing t h e appropriate expanded form of a i n Eq. (61).

Although a d e t a i l e d discussion of f i n i t e element formulations f o r

nonlinearly v i s c o e l a s t i c materials i s outside t h e scope of t h e present

investigation, t h e general procedure i s amply demonstrated by a simple

exmple. Consider t h e case of plane s t r e s s i n a t h i n membrane undergoing

large displacements but strains which are small i n comparison with unity,

and, for simplicity, assume t h a t the membrane i s constructed of a simple

l i n e a r l y v i s c o e l a s t i c material of Voigt type:

Here AM^ and Bapxp a r e arrays of material parameters and may be f'unctions

of time. Then Eq. (26) becomes

%'Mi
.. + caNcxI(si@ + C@~uMi)C&@)$.(dwm + $w',Jm)u&B@lp(6pm
+

We obtain a system of linear d i f f e r e n t i a l equations describing t h e in-plane

motion of f l a t v i s c o e l a s t i c elements by deleting products and squares of t h e

displacements and v e l o c i t i e s i n Eq. (63) and l j m i t i n g t h e ranges of lower-

case indices t o 2:

It i s seen t h a t t h e above procedure provides a systematic and r a t i o n a l means

for identifying t h e material damping c o e f f i c i e n t s f o r any material f o r which

t h e stress i s given e x p l i c i t l y as a function of s t r a i n s , s t r a i n r a t e s , and

other kinematic variables.

4.3 Elasto-plastic Materials. It i s not d i f f i c u l t t o modify f i n i t e element

 ~'ormiilationsso as t o account for yielding and p l a s t i c deformation of

m e t a l l i c membranes. Since, i n view of Eqs. (lo), t h e s t r a i n components

a r e uniform throughout each t r i a n g u l a r element, so a l s o i n t h e s t r e s s

[Eq. (54)l. Thus, i f l o c a l yielding i s imminent, it i s characterized by

uniform yielding and strain hardening of t h e associated l o c a l f i n i t e elements.

This means t h a t e l a s t i c - p l a s t i c boundaries cannot move continuously during

an incremental loading process; but t h i s c h a r a c t e r i s t i c of t h e d i s c r e t e

b

model need not lead to divergent or inconsistent results.

I n t h i s investigation, t h e approach suggested by Pope [22,231 i s

extended so as t o apply t o l a r g e displacements of e l a s t o - p l a s t i c membranes.

The e l a s t o - p l a s t i c behavior of a t y p i c a l f i n i t e element i s analyzed through

an incremental loading process. During each increment, t h e material responds

l i n e a r l y , b u t t h e o v e r a l l response obtained by summing t h e incremental

values may be highly nonlinear.

0 0
Let @,: ' PNiY and %i denote %?heknown values of t h e stress, s t r a i n ,
node forces, and node displacements a t some ref'erence s t a t e o i n 8. t;flical

f i n i t e element and l e t boa@, hap, 6pNi, and 811s-i denote small imxements i n

these quantities. If these increments a r e s u f f i c i e n t l y small, it can be

e a s i l y shown t h a t

e
where 8YaB i s t h e e l a s t i c s t r a i n increment and 8Yp i s t h e p l a s t i c s t r a i n
aa
increment. The term (bpi+ c @ ~ u & i) n Eq. (65c) represents t h e inflxence of

l a r g e displacements on t h e r e l a t i o n s h i p between s t r e s s e s and t h e node forces,

 28

The term cpMGi an be neglected i n the case of' sma.11 displacements.

The e l a s t i c s t r a i n increment i s r e l a t e d t o t h e s t r e s s increment according

to

0
i n which the arra-y Eo may, i n general, be a flmction of Yap.
aSw
The yield condition may be represented by a convex y i e l d surface i n

s t r e s s space which i s given by [243

f(Oij) = 0 (67)

The y i e l d function f i s symmetrical with respect t o G and 0.. and depends

13 J1
upon t h e s t r a i n h i s t o r y of t h e membrane. Since it i s asswned t h a t ' p l a s t i c

deformation i s independent of t h e hydrostatic s t r e s s , we can rewrite (67) i n

t h e form

where Tij i s t h e stress devia.tor:

When t h e s t r e s s increment bap((J3a = 033 = 0) does not have a positive

component i n the d i r e c t i o n of an inward normal t o t h e y i e l d surface and when

t h e stress point Oap l i e s on t h e yield surface, t h e material yields and t h e

p l a s t i e s t r a i n increment i s given by

The f a c t o r P depends upon t h e s t r a i n h i s t o r y and i s independent of a l l

components of boap except t h a t normal t o t h e f i e l d surface.

A t i n i t i a l yielding, t h e y i e l d condition (68) i s satisfied. After an

additional increment of p l a s t i c s t r a i n , t h e y i e l d condition i s given by

where taE3
describes t h e strain-hardening properties of t h e material.

Introducing Eqs. ( 6 6 ) , (681, and (70) i n t o (71) and sirnplifylng t h e r e s u l t s ,

we f i n d t h a t

Thus

(73 )

where

(74 1

Equation (77) represents t h e s t i f f n e s s r e l a t i o n between t h e incremental

loads 6pNi and t h e i r cooresponding incremental displacements. These equations

a r e inverted for each load increment; t h e solutions 6u a r e introduced i n t o

Ni
e
Eq. (75) t o obtain t h e associated s t r e s s increment. Incremental strains 6Y
aP
and 6Yp a r e then calculated by means of Eqs. (66) and (73). Once the
a@
i r ~ c r ~ ~ evalues
nta~ 6 and 6Y are determined, they a r e added
%ifd c r a f i y a$
algebraically to those of t h e reference state (i.e., Gi' (so a
,YO+ to

obtain a new reference state (uo' =

Ni
Gi+ 6%iy
a*' =
a@
'
'
0
up
+ boas, etc.)

and the process i s repeated f o r a new load increment. Following the

procedure indicated by Pope [221, during each cycle the f a c t o r P of each

element i s examined. If P 0, t h e element i s added t o the e l a s t i c region

of t h e membrane and t h e analysis i s repeated; i f t h e mean stress over a

load increment i s on or outside t h e yield surface, t h e analysis i s repeated

*with t h e element permitted t o deform p l a s t i c a l l y . Accuracy i s improved by

choosing t h e load increments such t h a t one element, a t t h e most, yields

during each load increment. Other d e t a i l s of t h e procedure a r e i d e n t i c a l

to those of the small displacement case and can be found In references f-221
and f-231.

5. FORMULBTION, OF THE STRUCTURAL PROBLEM

5.1 Global Equations of Motion. The nonlinear equations derived i n the

previous section describe the behavior of a s i n g l e f i n i t e membrane within

i t s l o c a l reference frame; these r e l a t i o n s are independen6 of t h e loading

on t h e membrane, t h e boundary conditions, or t h e location of the element i n

t h e assembled system. It i s now necessary t o connect t h e elements a t
appropriate node points and t o sum t h e i r properties so as t o represent a

pneumatic s t r u c t u r e of specified shape with specified boundary conditions.

To accomplish t h i s , it i s convenient to first rotate t h e node forces,

displacements, v e l o c i t i e s , accelerations and l o c a l cookdinates associated

w i t h each element so t h a t they a r e p a r a l l e l t o the global reference &me Zie

This i s accomplished through the transformations

 3 1.

where t h e underscore (-
e ) indicates t h a t no sum i s t o be taken on the repeated

index e .
-
I n these equations, pMie, uMie,
- . ..
and %ie
-ke, a r e respectively,

t h e node forces, displa,cements, v e l o c i t i e s , and a.ccel~erationsof node M of

element e i n t h e d i r e c t i o n of Zi,
-xie a r e t h e rotated local coordinates

defined i n Eq. (I), and Pkie i s t h e cosine of t h e angle between Xie and xke
The ranges of t h e indices i n these equa,tions a r e M, i, k, = 1, 2 , 3 and

e = 1, 2, ..., E, where Ee i s t h e t o t a l number of f i n i t e elements.

It was pointed out e a r l i e r t h a t the s e t of numbers Z N i ( i = 1,2,3;

q

N = 1,2, ..., n) describes t h e geometry of t h e assembled (connected)

-
system whereas xNie(N, i 1,2,3; e 1,2, ..., Ee) decribes t h a t of t h e
= =

individual elements. The connectivity of t h e system i s estaklished by

r e l a t i n g t h e members of t h e set ZNi t o those of TNie by the tro.xisformtion

-
X -- 0 z (M = 1, 2, 3; N = 1, 2, ..., n) (79)
Mie We N i

where

1 if node M of element e i s i d e n t i c a l t o node N i n

t h e assembled system
hl
MNe (80)
0 i f otherwise

The transformation indicated i n Eq. (79) defines a mapping of points i n t h e

-
set ZNi i n t o points i n xme and, i n e f f e c t , assembles t h e elements i n t o a

single unit. The process i s i l l u s t r a t e d symbolically far the case E, = 4,

n = 5 i n Fig. 4.
 I

LOCAL SYSTEMS GLOBAL SYSTEM

L

X 2
Mie Ni
e = 1,2,3,4 (ELEMENT INDEX) N = I ,2,3,4,6 (GLOBAL NODE INDEX)
M = 1,2,3 (LOCAL NODE INDEX) i 8 l,2,3 (COORDINATE INDEX)
i = 1,2,3 (COOROINATE INDEX 1

-x = a z
M ie MNe Mi

FIG. 4 Assembly of elements throii~;'htrnnsformatjons of points In Elobal s e t

i n t o local s e t s TMie.
'Ni
 32

Similarly, i f PNk, U%,

..
Urn, and UNk denote t h e values of node f o r c e s ,
displacements, velocities,and accelerations i n t h e assembled system, it

can be shown that

In t h i s case, t h e repeated indices M and e are summed throughout t h e i r

e n t i r e ranges: N= 1, 2, ..., n; e = 1, 2 , ...)Eea
Application of Eqs, (78) through (81) assembles t h e f i n i t e elements

i n t o a single d i s c r e t e system. When t h e l o c a l variables appearing i n l o c a l

eqwtions of motion [such as Eqs. (21) or (26)I a r e transformed i n accordance

with Eqs. (78) and (81), t h e r e s u l t i n g r e l a t i o n i s r e f e r r e d to as a g l o b a l

equation of motion, I n p a r t i c u l a r , Eq. (21) becomes

.. .
P 6
Ni Ni
= %U M iUN i + Po<dVo
VO

where Ee

e=l
I n t h i s case N, M = 1, 2 ..., n; R , S = 1, 2, 3 and t h e Lntegration i s

carried out over t h e e n t i r e volume Vo of t h e undeformed s t r u c t u r e . Equation

(82) i s t h e general global equation of motion of a f i n i t e element.

Global equations f o r t h e case of s t a t i c behavior are of s p e c i a l i n t e r e s t .

I n t h i s case,the l o c a l equations of motion reduce t o nonlinear s t i f f n e s s

relakions of t h e form
 33

where k(ume) i s t h e appropriate nonlinear function of t h e node displacements.

For example, t h e function k(% ) f o r synthetic rubbers, nonlinear3y e l a s t i c

e
materials, and Hookean metals are defined by Eqs. (48a), (52), and (54)

respectively. When t h e components of node forces and displacements a r e

' rotated so t h a t they a r e para,llel t o t h e coordinates xie' Eq. (84) becomes

J?
Mie
(85 I

Finally, t h e global s t i f f n e s s r e l a t i o n s a r e obtained through t h e transformations

imlicated i n Eq. (81):

where

, Boundazy conditions are applied by prescribing generalized (global) d i s -

placements a t appropriate boundary nodes. Then Eqs. (86) reduce t o a system

of independent nonlinear algebraic equations i n t h e unknown node displacements.

5.2 External Pressure. Up t o t h i s point, t h e r e l a t i o n s derived a r e applicable

only t o s i t u a t i o n s i n which t h e loads do not change i n direction as t h e

s t r u c t u r e deforms. Since, f o r pneumatic s t r u c t u r e s , t h i s i s obviously a

severe r e s t r i c t i o n , a procedure [7,111 for accbunting for changes i n t h e

external loading due t o deformation i s now examined.

It i s f i r s t assumed t h a t t h e dimensions of each f i n i t e element are

s u f f i c i e n t l y small t h a t t h e pressure q can be regarded as uniform over

t h e surfaces of each f i n i t e element. Then, i f A denotes t h e area of t h e

deformed element, the t o t a l force exerted normal t o t h e plane of the element

is
 34

L e t ni denote t h e components of a n outward u n i t vector normal t o A,

Then t h e components of the pressure force < are given by

N
qi = "$9

To determine the components ni, t h e coordinates of nodes of t h e deformed

element are denoted

For convenience i n writing, the o r i g i n of t h e reference frme y i s transferred

i
t o node 3 of t h e deformed element. If t h e r e s u l t i n g coordinate system i s

denoted zi, it follows t h a t

- (91)
Z N ~- Y N ~- Y3i

Now consider t w o unit vectors emanating from t h e o r i g i h of' the coordinates

z
i (node 3 ) .
The components ni of the u n i t normal are obtained by forming

the vector product of these two vectors:

ni =aL' i j k Z l j Z 2 k (92)

where Eijk i s the permutation symbol. Thus, equation (89) can be w r i t t e n

The net external f o r c e a t each node i s obtained by simply representing q by

t&ee forces, one a t each node, whose components are

Introducing Eq. (91) gives

This result defines t h e generalized external force i n the deformed elemen%

produced by external pressure. Note that no node i d e n t i f i c a t i o n index i s

needed since Qi i s t h e same f o r each node of t h e elementb

 35

To complete t h e analysis, these forces a r e now transformed i n t o components

-Qi p a r a l l e l to t h e l o c a l reference frame Zi. The r e s u l t f o r element e is
3

wherein t h e underscore again indicates t h a t no sum i s t o be taken on t h e

repeated index e.

In t h e case of pressure loadings, the component qie take t h e place of the

node forces 5Nie of Eqs. (78) and (85). It i s seldom necessary t o transform

these components i n t o the global system, however, since it i s more convenient

t o f i r s t transform displacements i n t o t h e zi system w i t h t h e a i d of Eq. (1)

and then t o transform t h e r e s u l t i n g forces i n t o t h e global system.

5.3 Solution of Nonlinear Equations, In t h e case of time-dependent phenomena,.

f i n i t e element formulations lead t o systems of simultaneous nonlinear

d i f f e r e n t i a l equations of t h e form indicated i n Eqs. (26), (63), and (82).

The solution of such systems of equations i s a formidable task, even with

t h e a i d of modern d i g i t a l computers. Generalizations of the well-known

Runga-Kutta techniques may lead t o acceptable r e s u l t s i n some cases; but

general procedures for solving such large systems of coupled nonlinear

d i f f e r e n t i a l equations are, a t b e s t , s t i l l i n t h e e a r l y stages of development.

It i s important t o note, however, t h a t numerical procedures a r e available for

t h e solution of nonlinear algebraic equations; and by incrementing t h e time

variable t , t h e o r i g i n a l s e t of d i f f e r e n t i a l equations reduce t o a system of

nonlinear algebraic equations for each time increment. Moreover, f i n i t e

element representations of s t a t i c behavior i n pneumatic s t r u c t u r e s also lead

t o systems of nonlinear algebraic equations. In view of t h i s , the solution

of large systems of nonlinear d i f f e r e n t i a l equ&tions i s not considered f u r t h e r

 36

i n t h i s study, Rather, consideration i s given t o procedures f o r solving

systems of nonlinear algebraic equations, it being understood t h a t , a t

t h e cost of g r e a t l y increasing t h e computing time, these procedures can

also be applied t o certain systems of nonlinear d i f f e r e n t i a l equations.

Several numerical schemes f o r solving simultaneous non1inea.r

algebraic equations are a v a i l a b l e i n t h e l i t e r a t u r e ; but not a l l of t h e s e

are s u i t a b l e f o r systems of equations a s l a r g e as those encountered i n t h e

present investigation. A comprehensive review and comparison of numerical

procedures for equations of t h i s type was r e c e n t l y contributed by Remmler,

Cawood , Stanton, and H i l l [251, wherein numerical experimentat ion showed

t h a t t h e c l a s s i c a l Newton-Raphson method and t h e Fletcher-Powell method are

among t h e most e f f i c i e n t and r e l i a b l e techniques available. To these may

be added t h e method of incremental loading, which i s somewhat r e l a t e d t o

the Newton-Raphson method, except t h a t the loading i s assumed t o b e applied

small increments during each of which the. s t r u c t u r e responds l i n e a r l y . This

l a t t e r technique i s p a r t i c u l a r l y well-suited f o r t h e analysis of s t a b i l i t y

and p l a s t i c behavior. The numerical r e s u l t s t o be presented subsequently

were obtained using v a r i a t i o n s of the Newton-Raphson method. Thus, f o r t h e

present discussion, it s u f f i c i e s t o merely o u t l i n e t h i s procedure. Details

of t h i s and other numerical procedures can be found i n t h e report by Remmler

e t a1 [ 2 5 ] and i n t h e papers by Sprang [26I and Brooks [271.

Consider a system of nonlinear s t i f f n e s s r e l a t i o n s of t h e form i n

Eq. ( 8 6 ) . Assuming t h a t a f t e r appropriate boundary conditions have been

applied t h e r e remain r unprescribed components of node displacements, t h i s

system represents a s e t of r independent nonlinear equations i n t h e unknown

displacements UNk' For simplicity, suppose t h a t these equations a r e represented

i n matrix form a s
 37

H ( U ) =O (97)

where H i s a r x 1 column matrix, each row of which represents an independent

nonlinear s t i f f n e s s equation, and u i s the solution vector. To solve these

equations, we expand H i n a Taylor s e r i e s about an a r b i t r a r y point uo which

represents an i n i t i a l estimate of the solution U . The vector u0 may,

f o r example, correspmd to t h e linearized solution. Taking only two terms,

we f i n d

H(U) = H( Uo) + J ( #io)(Li - 0')

where i s t h e jacobian matrix

Equation (98) i s l i n e a r i n u I( Solving t h i s equation, we f i n d

-1
U= Uo - J ( Uo) H ( Uo) (100)

The corrected solution U serves as t h e i n i t i a l estimate i n a second cycle,

and the process i s continued until a desired degree of accuracy i s obtained.

6. NUMERICAL AND EXPERIMENTAL RESULTS

I n t h i s section, we examine numerical r e s u l t s obtained by applying t h e

theory developed i n the preceding sections t o several representative

problems. Whenever possible, these r e s u l t s are compared with available

experimental or a n a l y t i c a l data.

6.1 S t r e s s Diffussion i n a Stiffened Panel, Ordinarily, stresses obtained

from f i n i t e element analyses based on approximate displacement f i e l d s are

less accurate than displacements obtained fmm such models, To get an

indication of t h e accuracy of s t r e s s e s derived from a r a t h e r coarse f i n i t e

element representation, and a l s o to examine an el-ementary problem involving

 38

composite bar and p l a t e elements, t h e problem of plane stress i n a reinforced

panel w a s considered as a simple first example, I n t h i s case, deformations

a r e assumed t o be small and e l a s t i c and t h e material i s assumed t o be homogeneous

and isotropic. Equations (54) and (56) a r e applicable except t h a t only

' in-plane deformtions a r e considered and products and squares of displacements

are neglected i n comparison with t h e displacements themselves.

The s t i f f e n e d panel shown i n Fig. 5a was analyzed using t h e f i n i t e

element representation i n Fig. 5b. Here a rectangular panel 0.127 mrp~thick,

s t i f f e n e d by longitudinal rods of area 1.613 em2 on t h e outside and 0,807 cm2

along t h e centerline, i s subjected t o concentrated forces of 1,587 kg, a s i s

indicated. An e l a s t i c modulus of 703,000 kg/cm 2 and a Poisson's r a t i o of 0.3

were used. The comwted v a r i a t i o n of t h e ilorma,l s t r e s s i n t h e exterior

longitudinal s t i f f e n e r s with t h e distance x from t h e fixed edge i s shown i n

Fig. 5c compared with r e s u l t s obtained fran an approximate theory developed

by Kuhn E281 We observe t h a t t h e coarse finite-element representation

yielded s t r e s s e s , i n t h i s case, which a r e i n close agreement with those

predicted by KUhn's theory.

6.2 Elasto-plastic Behavior of a Metallic Membrane. The f i n i t e element

formulation described i n Section 4.3 was used i n t h e analysts of p l a s t i c

behavior of a square aluminum membrane subjected t o external pressure. A

b i l i n e a r stress-strain l a w of t h e form indicated i n Fig. 6 was assumed with

(3 = 2,514 kg/cm 2 Yy = 0.0034, Ee = 20E = 740,000 kg/cm 2 e These properties

Y P
.correspond t o t h e aluminum a l l o y 2014-73 and agree closely with those used

i n experiments on rectangular s h e l l plating by NeuberL and Sommer C291, A

t h i n metallic sheet, 60 an. square and 0.14 cm. t h i c k , i s subjected t o a

uniformly d i s t r i b u t e d hydrostatic pressure. As t h e pressure i s slowly increased,

d,r:egion near t h e center of t h e p l a t e y i e l d s and p l a s t i c flow i s i n i t i a t e d .
 1,587 kg

IO00

900

8 00

@J
- 700

3E 600
tT
Y
u 500
in
2
a
400 -
A=+ FINITE ELEMENT
KUHN

300

200

IO0
4

IO 20 30 40

(C)
FIG. 5 Plane s t r e s s In a s t i f f e n e d panel.
FIG. 6 Bilinear stress-strain law.

60cm

FIG, 7 Finite-element representation of a square metallic membrane.

 33

This behavior w a s ana.l.y-zed numerically by using t h e f i n i t e element

representation Shawn i n Fig. 7. Figure 8 shows the computed variation i n the

center displacement of t h e sheet with external pressure compared with t h e

experimental r e s u l t s of Neubert and Sommer [291 and w i t h r e s u l t s obtained

using approximate t h e o r i e s proposed by Foppl' C301 and Hencky c311. Again

*
note t h a t t h e rather coarse network was adequate t o obtain displacements i n

excellent agreement with experimental data.

6.3 F i n i t e Stretching of a Rubber Sheet. I n order t o indicate t h e r a t e of

convergence of t h e r e s u l t s a s t h e finite-element network i s refined, we

reproduce here r e s u l t s similar t o those obtained e a r l i e r by Oden and Sato h11.

I n t h i s example, an i n i t i a l l y square rubber sheet, 0.127 cm. thick, i s

stretched i n i t s plane t o twice i t s o r i g i n a l length. The material i s assumed
*
t o be of the Mooney type, with material constants C1 and C2 of 1.75 and 0.15
, kg/cm 2 respectively. Thus, Eqs. (48) a r e applicable.
Various f i n i t e element representatiorls of t h e sheet a r e shown i n Fig. 9

along w i t h t h e v a r i a t i o n s i n t h e t o t a l edge force w i t h the t o t a l number of

f i n i t e elements. I n t h i s example, the edge force converged monotonically t o

approximately 16.32 kg .
6.4 I n f l a t i o n of an I n i t i a l l y F l a t Rubber Membrane. I n a recent paper, Hart-

Smith and Crisp c321 presented experimental data on the i n f l a t i o n of t h i n

rubber membranes. Although these investigators used an exponential form of

the s t r a i n energy function, s u f f i c i e n t in.formation was given t o deduce equivalent

Mooney constants f o r t h e material used. Specifically, we consider t h e i n f l a t i o n

of an i n i t i a l l y f l a t , c i r c u l a r , synthetic rubber membrane subjected t o uniform

external pressure. The membrane i s i n i t i a l l y 50.8 mm.in diameter and 0.2 mm.

thick and i s held fixed around i t s edges i n a metal clamp.

In t h e f i n i t e element analysis of t h i s membrane, it was assumed t h a t t h e

 50

40

E
E
U

2 30
W
25
w
0
4
1
a
v)
23 20
a
w
I-
z 0 EXPERIMENTS, NEUBERT zf SOMMER
W
0 -
--I
FINITE ELEMENTS
FOPPL'S THEORY
--- HENCKY'S THEORY

I I I I I 1 I I I I

I 2 3 4 5 6 7 8 9 10
PRESSURE ( ATMOSPHERES 1

FIG. 8 Variation i n center displacement of a square plate w i t h external pressure,

 t
20cm
1

E =4 16 t

17.0

cn 16.5
X
U

w
0
LT
0
LL
16.0
W

12 24 36 40 I 72
NUMBER OF ELEMENTS
FIG. 9 F i n i t e stretching of a square rubber shest.
 40

rubber possessed a s t r a i n energy function of t h e Mooney form I: Eq. (41)I SO

t h a t t h e nonlinear s t i f f n e s s r e l a t i o n s i n Eqs. (48) were applicable. Values

of t h e Mooney constants of C1 = 9.5C2 = 1.75 kg/m 2 were derived from t h e

data given i n C32I. The case considered i s t h a t i n which t h e membrane i s

subjected t o a uniform pressure of 0.097 kg/cm2. According t o the experimental

data, t h i s corresponds t o an extension r a t i o a t t h e crown of 5.5.

It i s important t o note t h a t t h e inflaking pressure i s a highly nonlinear

flmction of the extension r a t i o a t t h e crown and, consequently, of t h e displace-

ments. Thus, more than one equilibrium configuration can exist f o r a given

pressure. No provisions f o r determining a l l possible equilibrium states f o r

a specified pressure were incorporated i n the present analysis, and t h e

particu1a.r configuration obtained depends upon t h e choice of i n i t i a l values

employed i n t h e i t e r a t i v e solution of t h e nonlinear s t i f f n e s s r e l a t i o n s .

Several f i n i t e element networks were used i n t h e analysis, beginning

with a s i n g l e 30 degree element and eventually using t h e 10-element represent-

a t i o n shown i n Fig. 10. For a given f i n i t e element network, t h e r a t e of

convergence of t h e Newton-Raphson method depends on the choice of i n i t i a l

'
values of the displacements. Convergence r a t e s a.re considerably higher f o r

plane problems (such a s t h a t i n Fig. 9 ) than i n t h e case of large out-of-plane

deformations. I n t h e present example, r a t e s of convergence were increased 'by

f i r s t analyzing a coarse finite-element representation of t h e membrane using
a small number of i t e r a t i o n s . The r e s u l t s were then used as s t a r t i n g values

f o r a more refined representation, t h e displacements of t h e added node points

being obtained through l i n e a r interpolation.

Figure 11 shows t h e computed p r o f i l e of t h e i n f l a t e d sheet compared with

. t h e p r o f i l e s obtained experimentally and t h e o r e t i c a l l y by Hart-Smith and

Crisp. We observe t h a t t h e agreement i s q u i t e good, t h e m a x i m h difference

 r--50.8 mm 1-

FIG, 10 (a) Finite-elemegt represcritation of a circular membrane and

(b) a t y p i c a l 30 segment.
 -0- FINITE ELEMENTS

-
-I- MEASURE0 HART-SMITH
THEORY 1 8 CRISP

1 I.'

T-"

FIG. 11 Comparison of inflated profiles calculated by the finite-element

method with profiles obtained experimentally and theoretically
by Hart-Smith and Crisp c 3 d .
 41
between t h e displacements computed by the finite-element analysis and t h e

experimental values being approximakely six per cent.

6.5 Experiments on Rubber Membranes. A s a f i n a l example, we consider

b r i e f l y t h e r e s u l t s of experiments performed a t t h e Structures and Ma,terials

Laboratory of t h e University of Alabama Research I n s t i t u t e on t h i n n a t u r a l

rubber membranes. I n these experiments, c i r c u l a r disks, 0.0068 inc. (0.0173

cm.) thick and 15.0 i n , (38.1 cm.) i n diameter, of pure gum natural rubber sheet

were clamped around t h e i r edges i n a metal clamp, The disks were marked a t

t h e center and on t e n equally spaced concentric c i r c l e s , The disks were then

i n f l a t e d i n stages of pressures of approximately 100 mm. of water, which

corresponded t o a polar extension r a t i o of around A= 5. After r e s i s t i n g

maximum pressures f o r 45 minutes, the specimens were then deflated i n stages

u n t i l a l l applied pressure was removed.

Figure 12 shows t h e ex-perimental apparatus and a t y p i c a l i n f l a t e d c i r c u l a r

membrane. Figure 13 indicates t h e variation i n polar extension r a t i o with

pressure. It i s seen t h a t t h e behavior i s highly nonlinear and t h a t some

energy i s dissipated i n t h e unloading process. A residual extension a t t h e

center of a,pproximately J, = 1.25 was experienced, which was completely

recovered w i t h i n 24 hours af%er unloading.

A f i n i t e element representation with 96 elements, four nodes

'
along 30° radia,l l i n e s was used t o determine t h e i n f l a t e d p r o f i l e of a t y p i c a l

specimen subjected t o a pressure of 61 mm. of water. Again, t h e material was

2
assumed t o be of t h e Mooney type, with constants C1 = 1 . 1 4 kg/cm and

~ c2 = 0.14 kg/cm2 determined by t h e method of Hart-Smith and Crisp C32I and data i n
Figure 13. No attempt was made t o predict t h e obvious viscoelastic character

of t h e behavior. Results of these calculations a r e given i n Fig. 14,

 100
r lNFLATiNG

6 80
k3
LI
0
E
U
E 60
W
U
3
v,
v,
2e 40

20

a I I I * I

2 3 4 5 6

POLAR EXTENSION RATIO

FIG. 13 Variation i n polar extension r a t l o with i n t e r n a l pressure.

Fig. 12 I n f l a t e d membrane used i n experiments.

-..-
 - MEMBRANE
-01. FINITE ELEMENT

h AXIS OF SYMMETRY
19.05 cm

FIG. 14 Comparison of measured profile of an inflated rubber disk with

p r o f i l e obtained from finite element analysis.
 ACKNOWLEDGMENT

It i s a pleasure t o acknowledge t h e assistance of M r . T, .Sato,

who developed t h e computer programs used i n a11 of -the calculaLion$,

This work was supported by %he National Aeronautics and Space

Administration through General Research Grant NsG-381, and t h e

synthetic rubber materials used i n t h e experiment were donated by the

Materials Laboratory of NASA's.Marshal1 Space FlighL Center*

 43

REFERENCES

1. TURNER, M. J., DILL, E. H., MARTIN, H. C., and MELOSH, R. J., " h r g e
,
Deflections of Structures Subjected t o Heating and External Loads I' Journal
of t h e Aero/Space Sciences, V o l . 27, February 1960, pp. 97-106.

2. ARGYRIS, J. H., "Recent Advances i n Matrix Methods of S t r u c t u r a l Analysis, I t

Progress i n Aeronautical Sciences , Pergamon Press, Oxford, England, 1964,
PP. 115-1330

3. ARGYRIS, J. H., "Matrix Analysis of Three-Dimensional E l a s t i c Media-

Small and Large Displacements, ALAA Journal, V o l . 3, No. 1, January lN5,
?l

PP* 45-51.

4. GALLAGHER, R. H. and PADMG, J., "Discrete Element Approach t o S t r u c t u r a l

I n s t a b i l i t y Analysis," AIAA Journal, V o l . 1, No. 6, June 1963, pp. 1437-1439.

5. MARTIN, H. C., "On t h e Derivation of S t i f f n e s s Matrices f o r t h e Analysis

of Large Deflection and S t a b i l i t y Problems, " Conference on Matrix Methods
i n S t r u c t u r a l Mechanics, Wright-Patterson AFB, Dayton, Ohio, October 1965.
6. ODEN, J. T., "Calculation of Geometric S t i f f n e s s Matrices for Complex
Structures," AIAA Journal, Vole 4, No. 6, A u g u s t 1966, pp. 1480-1482.
'I (1
7. WISSMAIVN, J. W., Nwnerische Berechnung nichtlinear elastischer Koerper,
Dissertation, Hannover, 1963.

8. WISSMANN, J. W. , "Nonlinear Structural Analysis; Tensor Approach,''

Conference on Matrix Methods i n Structural Mechanics, Wright-Patterson AEB,
Dayton, Ohio, October 1965.

9- ODEN, J. T., "Analysis of Large Deformations of E l a s t i c Membranes by t h e

1 F i n i t e Element Method," Proceedings, IASS Congress on,..Large Span Shells,
,
Leningrad, September 1966.

10. ODEN, J. T . , "Numerical Formulation of Nonlinear E l a s t i c i t y Problems, ''

Journal of t h e Engineering Mechanics Division, ASCE, June 1967.

11. ODEN, J. T. and SATO, T., "Finite Strains and Displacements of E l a s t i c

Membranes by t h e F i n i t e Element Method,'' International Journal of Solids
and Structures, June 1967.

12. BECICER, E. B. , "A Numerical Solution of a Class of Problems of F i n i t e

E l a s t i c Deformation," Dissertation, University of California, Berkeley,
Ca,lifornia, 1966.

13. GREEN, A. E. and fYDKINS, J. E., Large E l a s t i c Deformations and Non-linear

Continuum Mechanics, Oxford University Press, London, 1960.

* 14. ERINGEN, A. C . , Nonlinear Theory of Continuous Media, McGraw H i l l Rook

Co., New York, 1962.
15. R I V L I N , R. S. and SAUNDERS, D. W., "Large E l a s t i c Deformations of
of Isotropic Materials. V I 1 Experiments on t h e Deformation of Rubber,
Philosophy Transactions of t h e Royal Society, No. 243, London, 1951,
pp. 251-288.

16. RIVLIN, R. S. and SAUNDERS, D. W., "The Free Energy of Deformation

for Vulcanized Rubbers, It Transactions. of t h e Faraday Societe, No. 48,
1952, pp. 200-206.
It
17. MOONEY, M., A Theory of Large E l a s t i c Deformations," Journal of
Applied Physics, No. 11, 1940, pp. 582-592.

18. TRELOAR, L. R. G., The Physics of Rubber E l a s t i c i t y , 2nd Edition,

Oxford University Press, London, 1958.

19. RIVLIN, R. S. , "Large E l a s t i c Deformations of Isotropic Materials. I

Fundamental Concepts, 11 Philosophy Transactions of the Royal Society,
No. 240, London, 1948, pp. 459-490.

20. HUTCHINSON, W. D., BECKER, G . W., and LANDETJ, R. F., ( 1Determination

of t h e Stored Energy Function of Rubber-like Materials I t Bulletin, Fourth
; vier , 1965.
21. BORESI, A. P., LANGHAAR, H. L. and MILL;ER, It. E., "Buckling of Axially-
Compr e s s ed Bilayer ed Fiber -Reinforced E l a s t i c Cylindr ica 1 She11," Develop-
ments i n Theoretical and Applied Mechanics. Vol. 2, Proceedings of t h e
Second Southeastern Conference on Theoretical and Applied Mechanics, W. A.
SHAW Ed., Pergamon Press , London, 1965, pp. 95-115.

22. POPE, G. G., "The Application of t h e Matrix Displacement Method i n

Plane Elasto-Plastic Problems, " Conference on Makrix Methods i n S t r u c t u r a l
Mechanics, Wright-Patterson AE'B, Dayton, Ohio, October 1965.

23. POPE, G. G. , "A Discrete Element Method f o r t h e Analysis of Plane

E l a s t o p l a s t i c S t r e s s Problems ," The Aeronautical Quarterly, Februa.ry 1966.

24. HILL, R., The Mathematical Theory of P l a s t i c i t & Oxford University Press,
London, 1950.

25. REMMLER, K. L., CAWOOD, D. W., STANTON, J. A., and HILL, R . , "Solutions
of Systems of Nonlinear Equations, I' Lockheed MSC/HREC , NASA 8-20178,
October 1966,
11
26. SPRANG, H. A,, "A Review of Minimization Techniques for Nonlinear Functions,
SLAM Review, Vol. 4, No. 4, October 1962, pp. 343-364.
27. BROOKS, S. H., "A Comparison of Maximum Seeking Methods ,
" Journal of
Operations Resea-rch, V o l . 7, 1959, pp. 430-457.

28. KUHN, P., Stresses i n A i r c r a f t S t r u c t w e s , McGraw-Hill Book Co., New

York, 1956, pp. 105-110.
29. FJEUBERT, M. and SOMMER, A., "Rechteckige Blechhaut unter g l e i c d s s i g
verteiltem Flusigkeitsdruck," Luftfahrtforschung, Vol. 17, No. 7, July 20,
1940, pp. 207-210 ( a l s o published as NACA Technical Memorandu, 965, December
1940).
30. FOPPL, A. and FOPPI;, L., Drang und Zwanq, V o l . 1, Oldenbourg, 1924.

Biegungsteif igkei't ,''

..
31. HENCKY, H., "Die Berechnung &inner rechteckiger Platten m i t verschwindender
9
(ZAMM), Bd. 1, Heft

32. HART-SMITH, L. J. and CRISP, J, D. C., "Large E l a s t i c Deformations of

Thin Rubber Membranes, If I n t e r n a t i o n a l Journal of Engineering Science, V o l . 5,
No. 1, January 1967, p-~?. 1-24.