<< Предыдущая стр. 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 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 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,
implicit double precision (a-h,o-z)
dimension lnods(nnode),ldgof(ndime,*),sigma(3,3),
& 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
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,