. 43
( 54 .)


I,J,K,L=1 V
and explain the meaning of CIJKL .
4. With the help of Exercise 6.3, derive a two-dimensional equation for the
external pressure component of the tangent matrix K p for a line element
along an enclosed boundary of a two-dimensional body under uniform
pressure p.
5. Recalling the line search method discussed in Section 7.6.2, show that
minimizing Π(·) = Π(xk + ·u) with respect to · gives the orthogonality

R(·) = uT R(xk + ·u) = 0



We have seen in the previous chapters that the solution to the nonlinear
equilibrium equations is basically achieved using the Newton“Raphson it-
erative method. In addition, in a ¬nite element context it is advisable to
apply the external forces in a series of increments. This has the advantage of
enhancing the converging properties of the solution and, where appropriate,
provides possible intermediate equilibrium states. Furthermore, it is clear
that the two fundamental quantities that facilitate the Newton“Raphson
solution are the evaluation of the residual force and the tangent matrix. In
this chapter we shall describe the FORTRAN implementation of the solu-
tion procedure in the teaching program FLagSHyP (Finite element LArGe
Strain HYperelasticity Program).
It is expected that the reader already has some familiarity with the com-
puter implementation of the ¬nite element method in the linear context.
Consequently this chapter will emphasize those aspects of the implemen-
tation that are of particular relevance in the nonlinear ¬nite deformation
context. In this respect, it is essential to understand two crucial routines.
Firstly, the master routine that controls the overall organization of the pro-
gram and, secondly, the subroutine elemtk. This latter routine computes
the equivalent nodal forces due to internal stress and the main components
of the tangent sti¬ness matrix. Consequently, it provides a vehicle for exam-
ining those aspects of the computation that are particular to ¬nite deforma-
tion analysis. In addition to the master routine and elemtk, the subroutine
bpress, which evaluates the equivalent nodal forces due to surface pres-
sure, and the corresponding tangent matrix components Kp and ksigma,


which computes and assembles the initial stress matrix Kσ , will also be de-
scribed. The remaining routines in the program are either relatively similar
to standard ¬nite element elasticity codes or are a direct implementation of
equations contained in the book.
The program description includes user instructions, a dictionary of vari-
ables and subroutines, and sample input and output for a few typical exam-
ples. The program can deal with a number of two-dimensional and three-
dimensional elements together with a variety of compressible and nearly
incompressible hyperelastic constitutive equations. It can be obtained from
the publishers via accessing the WWW site http://www.cup.cam.ac.uk.
The source program can be accessed via the ftp site: ftp.swansea.ac.uk
use “anonymous” user name and go to directory pub/civeng and get
Alternatively access the WWW site address:
http://www.swansea.ac.uk/civeng/Research/Software/ flagshyp/ to
obtain the program and updates. Alternatively, it can be obtained by
request to the authors or
e-mail j.bonet@swansea.ac.uk
r.d.wood@swansea.ac.uk. Finally, as a last resort, the program can be
obtained by sending a PC-formatted 3.5-inch disk together with a self-
addressed envelope to the authors™ address.


The input ¬le required by FLagSHyP is described below. The ¬le is free-
formatted, so items within a line are simply separated by commas or spaces.
Input Lines Comments
1-title(1:80) 1 Problem title
2-eltype(1:6) 1 Element type (see note 1)
3-npoin 1 Number of mesh points
4-ip,icode(ip), Coordinate data:
ip: Node number
icode(ip): Boundary code
(see note 2)
x(i,ip): x, y, z Coordinates
ndime: No. of dimensions
5-nelem 1 Number of elements
6-ie,matno(ie), Topology data:
element number
(lnods(i,ie), i=1, nnode) ie:
material number
nodal connectivities
No. of nodes per

7-nmats 1 Number of di¬erent materials
8-im,matyp(im) 2—nmats Material data:
im: material number
matyp(im): constitutive equation
props(i,im): properties (see note 3)
9-nplds,nprs,nbpel, 1 Loading data:
nplds: No. of loaded nodes
nprs: No. of non-zero
individual prescribed
nbpel: No. of line or surface
elements with applied
gravt: gravity vector
10-ip,(force(i),i=1,ndime) Point loads:
ip: node number
force: force vector
11-ip,id,presc Prescribed displacements:
ip: node number
id: spatial direction
presc: total nominal presc. displ.
12-je,(lbnod(i,ie),i=1, Pressure loads:
ie: surface element number
nbnod), press
lbnod(i,ie): nodal connectivities
nbnod: No. of nodes per element
press: nominal pressure
(see note 4)
13-nincr,xlmax,dlamb, 1 Solution control parameters:
nincr: number of load/
displacement increments
xlmax: max. value of load
scaling parameter
dlamb: load parameter increment
miter: maximum allowed no. of
iteration per increment
cnorm: convergence tolerance
searc: line search parameter
(if 0.0 not in use)
arcln: arc length parameter
(if 0.0 not in use)

Note 1: The following element types are recognized (see also Section 8.4
for more details about these elements):

tria3: 3-noded linear triangle;
tria6: 6-noded quadratic triangle;
quad4: 4-noded bilinear quadrilateral;
tetr4: 4-noded linear tetrahedron;
tetr10: 10-noded quadratic tetrahedron;
hexa8: 8-noded trilinear hexahedron.
Di¬erent element types cannot be mixed in a mesh. Given the el-
ement name the program automatically sets the number of nodes
per element (nnode), the number of dimensions (ndime) and the
number of Gauss points (ngaus). It also identi¬es the associated
type of surface or line element for pressure loads and sets the
corresponding number of nodes per element (nbnod) and Gauss
points (ngaub).

Note 2: The boundary codes are as follows:
0: free;
1: x prescribed;
2: y prescribed;
3: x, y prescribed;
4: z prescribed;
5: x, z prescribed;
6: y, z prescribed;
7: x, y, z prescribed.
Prescribed degrees of freedom are assumed to be ¬xed (that is,
no displacement) unless otherwise prescribed to be di¬erent from
zero in input items 9 and 11. If a displacement is imposed in
input items 9 and 11, the corresponding degree of freedom must
have already been indicated as prescribed in input item 4.

Note 3: The following constitutive equations have been implemented (see
also Section 8.6):
1: plane strain or three-dimensional compressible neo-Hookean;
2: not de¬ned;
3: plane strain or three-dimensional hyperelastic in principal
4: plane stress hyperelastic in principal directions;
5: plane strain or three-dimensional nearly incompressible neo-
6: plane stress incompressible neo-Hookean;
7: plane strain or three-dimensional nearly incompressible hy-
perelasticity in principal directions;

8: plane stress incompressible hyperelasticity in principal direc-
The corresponding list of properties to be read in Item 8 are
shown in the following table in which ρ is the density, » and µ
are the Lame coe¬cients, κ is the bulk modulus, and h is the
thickness for plane stress cases:

Type props(1) props(2) props(3) props(4)
1 ρ µ » ”
3 ρ µ » ”
4 ρ µ » h
5 ρ µ κ ”
6 ρ µ h ”
7 ρ µ κ ”
8 ρ µ h ”

Note 4: All the loads and prescribed displacements are nominal and will
be multiplied by the load scaling parameter » shown in Equa-
tion (7.74). The value of » is controlled indirectly by the arc-
length method or directly by the control parameters xlmax and
dlamb. For surface elements the direction of positive pressure is
given by the right-hand screw rule following the surface element
numbering. For line elements the positive pressure is at 90 de-
grees turning anticlockwise to the direction of increasing node

The somewhat contrived example shown in Figure 8.1 has been chosen
to illustrate as many diverse features of these input instructions as possible.
The required input ¬le is listed in Box 8.1. Note that point, gravity and
pressure loads, and prescribed displacements are all subject to the same
incrementation in the solution procedure.


The program FLagSHyP produces a largely unannotated output ¬le that is
intended to be an input ¬le for a postprocessor to be supplied by the user.
The contents and structure of this ¬le are shown in the table below.

y 3.4
1 8 9
Type 6
2 3 4
Type 4
4 5 6


. 43
( 54 .)