<<

. 48
( 54 .)



>>

xmu=props(2,im)/detf
xme=xmu-xlamb*log(detf)
call stress1(ndime,xmu,xme,btens,sigma)
call cisotp(ndime,xlamb,xme,ctens)



This constitutive equation is described in Section 5.4.3. The Cauchy
stresses are evaluated in terms of the b tensor from Equation (5.29), whereas
28 COMPUTER IMPLEMENTATION



the coe¬cients of the elasticity tensor are obtained in cisotp using Equa-
tion (5.37).

elemtk segment 5 “ Hyperelastic material in principal directions

else if(mat.eq.3) then
xlam=props(3,im)/detf
xmu=props(2,im)/detf
call jacobi(btens,stret,princ)
call stress3(ndime,xmu,xlam,detf,stret,princ,sigma,
& sprin)
call cprinc(ndime,xlam,xmu,stret,princ,sprin,ctens)



For this type of material the principal directions n± and principal stretches
are ¬rst computed by the routine jacobi and stored in princ and stret
respectively. Given these stretches and directions, the routine stress3 eval-
uates the Cartesian components of the Cauchy stress tensor with the help
of Equations (5.77) and (5.90). Finally, the routine cprinc uses Equa-
tions (5.86“7) and (5.91) to compute the corresponding elasticity tensor.

elemtk segment 6 “ Plane stress hyperelastic in principal direc-
tions

else if(mat.eq.4) then
det3d=detf**props(5,im)
xlam=props(3,im)/det3d
xmu=props(2,im)/det3d
call jacobi(btens,stret,princ)
call stress3(ndime,xmu,xlam,detf,stret,princ,sigma,
& sprin)
call cprinc(ndime,xlam,xmu,stret,princ,sprin,ctens)
vinc(ig)=vinc(ig)*props(4,im)*det3d/detf
stres(4,ig,ie)=props(4,im)*det3d/detf



This material is described in Section 5.6.7 and is implemented in a similar
way to the previous equation, the main di¬erence being the parameter γ,
which has been obtained in routine matprop using Equation (5.112) and is
stored as the ¬fth material parameter in props(5,im). Finally, the current
thickness is updated in terms of the volume ratio J, the area ratio j and
29
8.9 ROUTINE elemtk



the initial thickness stored in props(4,im) (see Exercise 3.3). The current
thickness is stored as the fourth stress for output purposes.

elemtk segment 7 “ Nearly incompressible neo-Hookean

else if(mat.eq.5) then
xmu=props(2,im)
call stress5(ndime,xmu,detf,btens,sigma)
call addpres(ndime,press,sigma)
call cdevia(ndime,xmu,detf,btens,ctens)
call cvolum(ndime,press,ctens)



This material is discussed in Sections 5.5.2 and 5.5.3. The deviatoric
Cauchy stresses are evaluated in stress5 using (5.51). Note that in the
mean dilatation method the pressure is constant over the element and there-
fore evaluated in Segment 2 outside the Gauss point loop. This pressure
value is now added to the deviatoric components in addpres. Once the
stresses are evaluated, routines cdevia and cvolum obtain the deviatoric and
volume components of the elasticity tensor by direct application of Equa-
tions (5.55a“b).

elemtk segment 8 “ Plane stress incompressible neo-Hookean

else if(mat.eq.6) then
xmu=props(2,im)
xme=xmu/(detf*detf)
xla=2.*xme
call stress6(ndime,xmu,detf,btens,sigma)
call cisotp(ndime,xla,xme,ctens)
vinc(ig)=vinc(ig)*props(4,im)/detf
stres(4,ig,ie)=props(4,im)/detf



This material is discussed in Exercise 5.1. Initially the e¬ective shear
coe¬cient µ is evaluated dividing µ by j 2 = det2—2 b. As shown in Ex-
ercise 5.1, an e¬ective lambda coe¬cient (see Equation (5.37)) emerges as
twice the e¬ective shear coe¬cient. The thickness is ¬nally computed as in
Material 4 with the exception that due to incompressibility det3d is equal
to 1.
30 COMPUTER IMPLEMENTATION




else if(mat.eq.7) then
xmu=props(2,im)/detf
elemtk segment 9 “ Nearly incompressible in principal directions
xlam=-2.*xmu/3.
call jacobi(btens,stret,princ)
call stress3(ndime,xmu,xlam,detf,stret,princ,
& sigma,sprin)
call addpres(ndime,press,sigma)
call cprinc(ndime,xlam,xmu,stret,princ, sprin,
& ctens)
call cvolum(ndime,press,ctens)



This material is discussed in Section 5.6.6. Stretches and principal di-
rections are ¬rst evaluated in Jacobi. Routine stress, which implemented
Equation (5.90), is reused with » being set to ’2µ/3 to give the deviatoric
Cauchy stresses in accordance with Equation (5.103). The internal hydro-
static pressure is the added in addpres. The deviatoric components of the
elasticity tensor are obtained using the same routine cprinc employed for
Material 3 by setting the argument xlam to ’2µ/3. Finally the volumetric
component is obtained as per Material 5.

elemtk segment 10 “ Plane stress incompressible in directions

else if(mat.eq.8) then
xmu=props(2,im)
xlam=2.*xmu
call jacobi(btens,stret,princ)
call stress3(ndime,xmu,xlam,detf,stret,princ,
& sigma,sprin)
call cprinc(ndime,xlam,xmu,stret,princ,sprin,
& ctens)
vinc(ig)=vinc(ig)*props(4,im)/detf
stres(4,ig,ie)=props(4,im)/detf
endif



This material is identical to Type 4 (see Segment 6), but because of
¯
incompressibility » ’ ∞ and hence γ = 0, J = 1, and » = 2µ.
31
8.10 ROUTINE ksigma




elemtk segment 11 “ Internal forces and tangent sti¬ness

call internal(ndime,nnode,lnods(1,ie),ldgof,
& sigma,elecd(1,1,ig),vinc(ig),
& resid,react)
call kconst(ndime,nnode,lnods(1,ie),ldgof,
& ctens,elecd(1,1,ig),vinc(ig),kprof,
& stifd,stifp)
call ksigma(ndime,nnode,lnods(1,ie),ldgof,
& sigma,elecd(1,1,ig),vinc(ig),kprof,
& stifd,stifp)
ks=1
do id=1,ndime
do jd=id,ndime
stres(ks,ig,ie)=sigma(id,jd)
ks=ks+1
enddo
enddo
enddo
enddo
return
end


Recalling that we remain within the Gauss loop, the Cauchy stress com-
ponents and the Cartesian derivatives of the shape functions are used in
routine internal to compute and assemble the equivalent nodal forces T
employing Equation (7.15). Subroutine kconst then evaluates and assem-
bles the constitutive component of the tangent matrix according to the in-
dicial Equation (7.35). The initial stress matrix is obtained and assembled
in ksigma using Equation (7.45). This routine will be described in detail
in the next section. Finally, the stresses are copied into the array stress
ready for the output routine.


8.10 ROUTINE ksigma

Although all sti¬ness matrix calculations are a direct implementation of the
relevant equations given in the text, it is instructive to consider the com-
putation and assembly of the initial stress matrix Kσ , which is of particular
signi¬cance in ¬nite deformation analysis.
32 COMPUTER IMPLEMENTATION




ksigma segment 1 “ Sti¬ness coe¬cient evaluation

subroutine ksigma(ndime,nnode,lnods,ldgof,sigma,
& gradn,vinc,kprof,stifd,stifp)
implicit double precision (a-h,o-z)
dimension lnods(nnode),ldgof(ndime,*),sigma(3,3),
& gradn(ndime,nnode),kprof(*),stifd(*),
& stifp(*)
do in=1,nnode
ip=lnods(in)
do jn=in,nnode
jp=lnods(jn)
sum=0.0
do kd=1,ndime
do ld=1,ndime
sum=sum+gradn(kd,in)*sigma(kd,ld)*
& gradn(ld,jn)*vinc
enddo
enddo



This routine is called from elemtk within a Gauss and element loop.
The array gradn contains the current Cartesian derivatives ‚Na /‚xi , and
the variable sum is computed according to Equation (7.45c).

ksigma segment 2 “ Assembly

do id=1,ndime
if=ldgof(id,ip)
jf=ldgof(id,jp)
if((if.le.0).or.(jf.le.0)) goto 10
if(jf.eq.if) then
stifd(if)=stifd(if)+sum
else
jc=max(jf,if)
ir=min(jf,if)
lp=kprof(jc)-jc+ir+1
stifp(lp)=stifp(lp)+sum
endif
10 enddo
33
8.11 ROUTINE bpress



enddo
enddo
return
end



Variables if and jf indicate the degree-of-freedom numbers correspond-
ing to the direction id of nodes in and jn. When both if and jf are active
and therefore positive, the value stored in sum needs to be assembled. For
the case where if is equal to jf, this value is added to the diagonal sti¬ness
vector stifp. Otherwise the correct position in stifp is obtained using
kprof as explained in Section 8.5.


8.11 ROUTINE bpress

bpress segment 1 “ Surface element loop

subroutine bpress(ndime,nnodb,ngaub,nbpel,xlamb,
& elbdb,lbnod,press,x,ldgof,
& tload,react,resid,kprof,stifd,

<<

. 48
( 54 .)



>>