implicit double precision (a-h,o-z)

dimension elbdb(ndime,nnodb+1,*),lbnod(nnodb,*),

& press(*),x(ndime,*),

& ldgof(ndime,*),tload(*),

& react(*),resid(*),kprof(*),

& stifd(*),stifp(*),dxis(3,2)

data ((dxis(i,j),i=1,3),j=1,2)

& /0.,0.,0.,0.,0.,0./

if(nbpel.eq.0) return

do ie=1,nbpel

epres=press(ie)

do ig=1,ngaub

34 COMPUTER IMPLEMENTATION

For cases where there are surface or line elements with pressure applied,

this routine ¬rst loops over these elements and then over their Gauss points.

bpress segment 2 “ Equivalent forces and sti¬ness component

if(ndime.eq.2) dxis(3,2)=-1.

do id=1,ndime

do jd=1,ndime-1

dxis(id,jd)=0.0

do in=1,nnodb

ip=lbnod(in,ie)

dxis(id,jd)=dxis(id,jd)+x(id,ip)*

& elbdb(jd+1,in,ig)

enddo

enddo

enddo

call pforce(ndime,nnodb,epres,xlamb,elbdb(1,1,ig),

& lbnod(1,ie),ldgof,tload,resid,react,

& dxis)

apres=epres*xlamb

call kpress(ndime,nnodb,apres,elbdb(1,1,ig),

& lbnod(1,ie),ldgof,dxis,kprof,stifd,

& stifp)

enddo

enddo

return

end

For each Gauss point, ¬rst the tangent vectors ‚x/‚ξ and ‚x/‚· are

evaluated. For 2-D cases the vector ‚x/‚· is set to [0, 0, ’1]T , so that when

the cross product with ‚x/‚ξ is taken the correct normal vector is obtained.

These tangent vectors are stored in the array dxis. Using these vectors

and the pressure value epress, the routine pface evaluates and assembles

into tload in a standard manner the equivalent nodal force components

corresponding to this Gauss point using Equation (7.15c). At the same time

these nodal forces are added to the residual force vector and to the reaction

vector. Finally, the routine kpress is called to evaluate the initial pressure

sti¬ness component.

35

8.12 EXAMPLES

y

1

1

1/4 1/4

1/4

3/4

1

1/2 1 x

2

1

FIGURE 8.5 Patch test.

8.12 EXAMPLES

A limited number of examples are described in this section in order to show

the capabilities of FLagSHyP and illustrate some of the di¬culties that may

arise in the analysis of highly nonlinear problems.

8.12.1 SIMPLE PATCH TEST

As a simple check on the code the nonlinear plane strain patch test example

shown in Figure 8.5 is studied. Two irregular six-noded elements making up

a square are employed and boundary displacements are prescribed as shown

in the ¬gure in order to obtain a uniform deformation gradient tensor given

by,

20

F= ; J = 3/2

0 3/4

Assuming Material 1 with values » = 100 and µ = 100, the program gives the

correct state of uniform stress that can be evaluated from Equation (5.29)

as,

µ » 227.03 0 40

b = FFT =

σ= (b ’ I) + (ln J)I = ;

0 ’2.136 0 9/16

J J

8.12.2 NONLINEAR TRUSS

The nonlinear truss example already described in Chapter 1, Section 1.3.2,

can now be re-examined. For this purpose, a single four-noded element is

36 COMPUTER IMPLEMENTATION

99 1

0.4

down

0.3

up

0.2

0.1

F/EA

y

99 0

’0.1

(not to scale) ’0.2

’0.4

1 x

’4 ’3 ’2 ’1 0 1 2 3 4

x/L

(a) (b)

FIGURE 8.6 Truss example: (a) geometry; (b) load-displacement curve.

used with the support conditions shown in Figure 8.6. In order to mini-

mize the two-dimensional e¬ects a large length-to-thickness ratio of approx-

imately 100 is used. The initial angle is 45 √ degrees as in Figure 1.5 of

Chapter 1 and the initial thickness is set to 1/ 2 so that the initial area of

the truss is 1. Material 8 is used with µ = 1/3, which for ν = 0.5 implies

E = 1. The displacement of the top node is prescribed either upward or

downward in small increments of 0.5. The resulting vertical force is shown

in Figure 8.6(b) and is practically identical to that obtained in Section 1.3.2

using the logarithmic strain equation and shown in Figure 1.5.

It is worth noting that when the truss is pushed downward and reaches

the ’45 degree position the program appears to fail to converge. The reason

for this is simply that at this position the stresses in the truss, and hence

all the reactions, are zero to within machine accuracy. Consequently, when

the convergence norm is evaluated by dividing the norm of the residuals by

the norm of the reactions a practically random number is obtained. In order

to escape from this impasse, an arti¬cially large convergence norm must be

chosen for the increment that corresponds to this position (0.98 has been

used by the authors). The computation can then be continued using normal

convergence norms for the rest of the increments. Note that the restarting

facilities in FLagSHyP enable this arti¬ce to be easily performed.

8.12.3 STRIP WITH A HOLE

This is a well-known plane stress hyperelasticity example where an initial

6.5 — 6.5 — 0.079 mm3 strip with a hole 0.5 mm in diameter is stretched in

37

8.12 EXAMPLES

6

Dimensionless force F/Aµ

5

6.5 mm

4

3

F

6.5 mm

F

0.5 mm

2

1

0

1 2 3 4 5 6

Longitudinal stretch »

(b)

(a)

(c)

FIGURE 8.7 Strip with hole: (a) geometry; (b) load-stretch curve; (c) ¬nal mesh.

the horizontal direction and constrained in the vertical direction as shown in

Figure 8.7. Given the two planes of symmetry, 100 four-noded quadrilateral

elements are used to describe a quarter of the problem. A plane stress incom-

pressible neo-Hookean material model is used with µ = 0.4225 N/mm2 . The

strip is stretched to six times its horizontal length using only ¬ve increments.

The load-displacement curve is shown in Figure 8.7(b), and Figure 8.7(c)

displays the ¬nal mesh.

8.12.4 PLANE STRAIN NEARLY INCOMPRESSIBLE STRIP

This well-known example has been included in order to illustrate the di¬cul-

ties that can be encountered using the penalty type of nearly incompressible

plane strain or three-dimensional materials implemented in FLagSHyP in

conjunction with a displacement control process. In this case a 20 — 20 mm2

strip is clamped and stretched as shown in Figure 8.8. Because of the

symmetry only a quarter of the problem is actually modelled using 256

four-noded mean dilatation quadrilateral elements. Material 5 is used with

µ = 0.4225 N/mm2 and κ = 5 N/mm2 . The ¬nal mesh is shown in Fig-

ure 8.8(b), where for a total horizontal stretch of three a vertical stretch of

38 COMPUTER IMPLEMENTATION

20 mm

20 mm

(a)

(b)

FIGURE 8.8 Incompressible strip: (a) geometry; (b) ¬nal mesh.

0.3711 is observed (smaller values are obtained as κ is increased to reach the

incompressible limit.)

Although the material is slightly compressible (ν ≈ 0.46) the solution

requires 200 increments in order to ensure that the Newton“Raphson it-

erative process reaches convergence in a reasonable number of steps. In

fact if the value of κ is increased in order to approach the incompressible

limit (to say 50 or 500), an even larger number of increments is needed.

This is in clear contrast with the previous example, which could be run in

¬ve increments. Unfortunately, using large increments in the present case

leads to large changes in the volumes of those elements containing nodes

with prescribed displacements, which, because of the relatively large bulk

modulus, give extremely large internal pressures at these elements. These

unrealistic pressure values lead in turn to very large residual forces from

which the Newton“Raphson process fails to converge. This problem is typi-

cally solved using the augmented Lagrangian technique whereby equilibrium

is ¬rst reached using small initial values of κ, which are then progressively

increased, always maintaining equilibrium by further iterations if necessary,

until the desired κ/µ ratio is reached. This technique enables very large

¬nal values of κ/µ to be used without increasing the number of increments

needed. In order to keep the code as simple as possible, however, this tech-