22 COMPUTER IMPLEMENTATION

converged solution (for instance due to a large tolerance value cnorm) from

propagating throughout the solution process. The further addition to the

residual forces due to an increment in applied surface pressure is evaluated by

calling bpress with the ¬fth argument set to the load parameter increment

dlamb. At the same time bpress evaluates the addition to the tangent

matrix due to the change in its initial pressure component that results from

the increment in the magnitude of the pressure.

If prescribed displacements have been imposed, then the current incre-

ment in their value will immediately change the current geometry and neces-

sitate a complete re-evaluation of equivalent internal and surface pressure

forces and the tangent matrix (this is e¬ectively a geometry update). In

this case, subroutine prescx will reset the current geometry for those nodes

with prescribed displacements based on their initial coordinates, the nominal

value of the prescribed displacement, and the current load scaling param-

eter xlamb. Subroutine initrk initializes the tangent matrix to zero and

the residuals to the current value of external and point loads. Subsequently,

subroutine bpress, now called with the ¬fth argument set to xlamb, adds

the total value of the nodal forces due to surface pressure and obtains the

corresponding initial surface pressure tangent matrix component. Finally,

routine elemtk evaluates the current equivalent nodal forces due to stresses,

subtracts these from the residual vector, and computes the constitutive and

initial stress components of the tangent matrix.

¬‚agshyp segment 6 “ Newton“Raphson loop

niter=0

rnorm=2.*cnorm

do while((rnorm.gt.cnorm).and.(niter.lt.miter))

niter=niter+1

call datri(stifp,stifp,stifd,kprof,

& ndgof,.false.,6)

call dasol(stifp,stifp,stifd,resid,displ,

& kprof,ndgof,6,rtu0)

if(arcln.ne.0.0) then

call dasol(stifp,stifp,stifd,tload,resid,

& kprof,ndgof,6,r)

call arclen(ndgof,niter,arcln,displ,resid,

& xincr,xlamb,dlamb)

endif

23

8.8 MAIN ROUTINE flagshyp

The do while loop controls the Newton“Raphson iteration process. This

continues until the residual norm rnorm is smaller than the tolerance cnorm

and, of course, while the iteration number niter is smaller than the maxi-

mum allowed miter. (Note that rnorm is initialized in such a way that this

loop is completed at least once.) The ¬rst step of the Newton“Raphson pro-

cess is the factorization of the current tangent matrix performed in dasol.

Note that for the ¬rst iteration of the ¬rst increment, the tangent matrix

is available either from initel in Segment 3 or from elemtk in Segment

5. Subsequently, the current matrix is that evaluated at the previous iter-

ation unless nonzero displacements have been prescribed, in which case it

is re-evaluated in Segment 5. Factorization is followed by the backsubsti-

tution routine dasol from which the incremental displacements displ are

obtained. If the arc-length method is in use, dasol is called again to evalu-

ate the auxiliary vector uF (see Equation 7.80c), which is temporarily stored

in the array resid. The arc-length manipulation is completed by evaluating

a new value of the scaling parameter xlamb in the routine arclen.

¬‚agshyp segment 7 “ Solution update and equilibrium check

loop

eta0=0.

eta=1.

nsear=0

rtu=rtu0*searc*2.

do while((abs(rtu).gt.abs(rtu0*searc)).and.

& (nsear.lt.5))

nsear=nsear+1

call update(ndime,npoin,ldgof,x,displ,eta-eta0)

call initrk(ndgof,nprof,negdf,xlamb,eload,tload,

& resid,react,stifd,stifp)

call bpress(ndime,nnodb,ngaub,nbpel,xlamb,elbdb,

& lbnod,press,x,ldgof,tload,react,resid,

& kprof,stifd,stifp)

call elemtk(ndime,nnode,ngaus,nstrs,nelem,x,x0,

& eledb,lnods,matno,matyp,props,ldgof,

& stres,resid,kprof,stifd,stifp,react,

& vinc,elecd,elacd,vol0)

call search(ndgof,resid,displ,eta0,eta,rtu0,rtu)

enddo

call checkr(incrm,niter,ndgof,negdf,xlamb,resid,

24 COMPUTER IMPLEMENTATION

& tload,react,rnorm)

enddo

If the line search parameter searc is input as 0.0, the program resets the

value in the routine incontr to 105 . This simple device ensures that the do

while loop is performed at least once per Newton“Raphson iteration in order

to update the solution and check for equilibrium. If searc is not 0.0, the line

search procedure is active and this inner line search loop may be performed

up to ¬ve times. The routine update adds the displacements evaluated in

Segment 6, scaled by the line search factor · (initially set to 1), to the

current nodal positions. Given this new geometry, the routines initrk,

bpress, and elemtk serve the same purpose as described in Segment 5. The

routine search evaluates a new line search factor · if required. Finally,

checkr computes the residual norm rnorm given the residual force vector

resid and the total external force vector tload.

¬‚agshyp segment 8 “ Retry or output results

if(niter.ge.miter) then

write(6,100)

call restar1(title,eltyp,ndime,npoin,nnode,

& ngaus,nstrs,nelem,ndgof,negdf,

& nprs,nprof,nmats,incrm,xlamb,

& nbpel,nnodb,ngaub)

call restar2(ndime,npoin,nnode,ngaus,nelem,

& ndgof,negdf,nprof,nmats,nnodb,

& ngaub,nbpel,matyp,props,matno,

& lnods,x,x0,kprof,stifd,stifp,

& resid,eload,ldgof,icode,eledb,

& pdisp,vol0,elbdb,lbnod,press)

call incontr(nincr,xlmax,dlamb,miter,cnorm,

& searc,arcln,5)

else

call output(ndime,nnode,ngaus,nstrs,npoin,

& nelem,eltyp,title,icode,incrm,

& xlamb,x,lnods,ldgof,matno,stres,

& tload,react)

call dump(title,eltyp,ndime,npoin,nnode,

& ngaus,nstrs,nelem,ndgof,negdf,

25

8.9 ROUTINE elemtk

& nprs,nprof,nmats,incrm,xlamb,

& nbpel,nnodb,ngaub,matyp,props,

& matno,lnods,x,x0,kprof,stifd,

& stifp,resid,eload,ldgof,icode,

& eledb,pdisp,vol0,elbdb,lbnod,

& press)

endif

enddo

stop

100 format(™Solution not converged, restarting from

previous step™)

end

In the event that the Newton“Raphson iteration does not converge, the

user has the opportunity of interacting with the program to restart from the

previous converged increment with revised control parameters. Note that

interactively setting the number of increments to any value smaller than

the current increment number will result in the termination of the program.

Subroutine output writes out the output ¬le described in Section 8.3 while

dump creates a binary ¬le with enough information to enable future restarts.

8.9 ROUTINE elemtk

This routine evaluates equivalent nodal forces from element stresses and, at

the same time, the main components of the tangent matrix. The coe¬cients

of both these items are immediately assembled upon computation. Observe

that in the presence of applied surface pressure, the initial pressure com-

ponent of the tangent matrix has previously been evaluated and assembled

by the routine bpress. This routine will be described in more detail in

Section 8.11. Again recall that the concatenation of the following segments

constitutes the complete routine.

elemtk segment 1 “ Element loop

subroutine elemtk(ndime,nnode,ngaus,nstrs,nelem,x,x0,

& eledb,lnods,matno,matyp,props,ldgof,

& stres,resid,kprof,stifd,stifp,react,

& vinc,elecd,elacd,vol0)

26 COMPUTER IMPLEMENTATION

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

dimension lnods(nnode,*),props(8,*),ldgof(ndime,*),

& kprof(*),stifd(*),stifp(*),x(ndime,*),

& x0(ndime,*),matno(*),

& eledb(ndime+1,nnode+1,ngaus),

& elacd(ndime,*),resid(*),ctens(3,3,3,3),

& sigma(3,3),btens(3,3),matyp(*),react(*),

& stres(nstrs,ngaus,*),stret(3),princ(3,3),

& sprin(3),vol0(*),elecd(ndime,nnode,*),

& vinc(*)

data ((btens(i,j),i=1,3),j=1,3)

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

do ie=1,nelem

im=matno(ie)

mat=matyp(im)

call gradsh(ie,ndime,nnode,ngaus,eledb,lnods,x,

& elecd,vinc)

Elemtk loops ¬rst over all the elements. For each element the stan-

dard routine gradsh computes current Cartesian derivatives ‚Na /‚x of the

shape functions using Equation 7.11. For convenience gradsh computes

these derivatives at all Gauss points of the element and stores them in elecd.

Additionally, gradsh returns the weighted Jacobian per Gauss point.

elemtk segment 2 “ Mean dilatation implementation

if(mat.eq.5) then

call getjbar(ngaus,vol0(ie),vinc,barj)

xkapp=props(4,im)

press=xkapp*(barj-1.)

xkapp=xkapp*barj

call kvolume(ndime,nnode,ngaus,xkapp,vinc,elecd,

& elacd,lnods(1,ie),ldgof,kprof,stifd,

& stifp)

else if (mat.eq.7) then

call gebarj(ngaus,vol0(ie),vinc,barj)

xkapp=props(4,im)

press=xkapp*log(barj)/barj

xkapp=(xkapp/barj)-press

27

8.9 ROUTINE elemtk

call kvolume(ndime,nnode,ngaus,xkapp,vinc,elecd,

& elacd,lnods(1,ie),ldgof,kprof,stifd,

& stifp)

endif

For nearly incompressible Materials 5 and 7, the mean dilatation tech-

nique described in Chapters 6 and 7 is used to obtain ¬rst the internal

¯

element pressure p (press) from the volume ratio J (barj) using Equa-

tion 6.50a. The dilatational component Kκ of the tangent matrix is com-

puted and assembled in routine kvolume using Equation 7.60.

elemtk segment 3 “ Gauss loop and evaluation of b

do ig=1,ngaus

call leftcg(ndime,nnode,lnods(1,ie),x0,

& elecd(1,1,ig),detf,btens)

For each element the routine now loops over all Gauss points. Given

that the current rather than initial Cartesian derivatives are known, routine

leftcg ¬rst evaluates the inverse deformation gradient as F ’1 = ‚X/‚x =

Na , then obtains F , and ¬nally b = F F T . This routine also

a Xa —

evaluates the pointwise volume ratio J = det F (or area ratio in the case of

plane stress). Using this left Cauchy“Green tensor, the seven segments that

follow will obtain the Cauchy stresses and elasticity tensor for each of the

seven material types implemented.

elemtk segment 4 “ Compressible neo-Hookean material

if(mat.eq.1) then

xlamb=props(3,im)/detf