sym

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

28 DISCRETIZATION AND SOLUTION

condition,

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

CHAPTER EIGHT

COMPUTER IMPLEMENTATION

8.1 INTRODUCTION

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,

1

2 COMPUTER IMPLEMENTATION

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:

flagshyp.f.

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.

8.2 USER INSTRUCTIONS

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:

npoin

ip: Node number

(x(i,ip),i=1,ndime)

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:

nelem

element number

(lnods(i,ie), i=1, nnode) ie:

material number

matno(ie):

nodal connectivities

lnods(i,ie):

No. of nodes per

nnode:

element

3

8.2 USER INSTRUCTIONS

7-nmats 1 Number of di¬erent materials

8-im,matyp(im) 2—nmats Material data:

im: material number

matyp(im): constitutive equation

type

props(i,im): properties (see note 3)

props(1,im),

props(2,im),...

9-nplds,nprs,nbpel, 1 Loading data:

nplds: No. of loaded nodes

(gravt(i),i=1,ndime)

nprs: No. of non-zero

individual prescribed

displacements

nbpel: No. of line or surface

elements with applied

pressure

gravt: gravity vector

10-ip,(force(i),i=1,ndime) Point loads:

nplds

ip: node number

force: force vector

11-ip,id,presc Prescribed displacements:

nprs

ip: node number

id: spatial direction

presc: total nominal presc. displ.

12-je,(lbnod(i,ie),i=1, Pressure loads:

nbpel

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/

miter,cnorm,searc,

displacement increments

arcln

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):

4 COMPUTER IMPLEMENTATION

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

directions;

4: plane stress hyperelastic in principal directions;

5: plane strain or three-dimensional nearly incompressible neo-

Hookean;

6: plane stress incompressible neo-Hookean;

7: plane strain or three-dimensional nearly incompressible hy-

perelasticity in principal directions;

5

8.3 OUTPUT FILE DESCRIPTION

8: plane stress incompressible hyperelasticity in principal direc-

tions.

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

numbering.

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.

8.3 OUTPUT FILE DESCRIPTION

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.

6 COMPUTER IMPLEMENTATION

y 3.4

0.25

1.2

1 8 9

7

Type 6

2 3 4

Type 4

4 5 6