<<

. 47
( 54 .)



>>

the previous increment by the routine force. This prevents errors in the
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

<<

. 47
( 54 .)



>>