Civil Engineering Reference
In-Depth Information
In the limit as
β
t
→∞
, this expression tends to,
3
σ
3
(K
p
−
1
)
D
f
=
(6.85)
(K
p
+
2
)
although for numerical purposes, undrained behaviour is essentially captured for
β
t
>
20
Program 6.9 Axisymmetric “undrained” strain of an elastic-plastic (Mohr-Cou-
lomb) solid using 8-node rectangular quadrilaterals. Viscoplastic strain method.
PROGRAM p69
!-------------------------------------------------------------------------
! Program 6.9 Axisymmetric 'undrained' strain of an elastic-plastic
! (Mohr-Coulomb) solid using 8-node rectangular
! quadrilaterals. Viscoplastic strain method.
!-------------------------------------------------------------------------
USE main; USE geom; IMPLICIT NONE
INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
INTEGER::fixed_freedoms,i,iel,incs,iters,iy,k,limit,ndim=2,ndof=16,nels, &
neq,nip=4,nn,nod=8,nodof=2,nr,nst=4,nxe,nye
REAL(iwp)::bulk,c,cons,det,dq1,dq2,dq3,dsbar,dt,d4=4.0_iwp, &
d180=180.0_iwp,e,f,lode_theta,one=1.0_iwp,phi,pi,presc,psi,sigm,snph, &
penalty=1.e20_iwp,tol,two=2.0_iwp,v,zero=0.0_iwp
CHARACTER(LEN=15)::element='quadrilateral'; LOGICAL::converged
!-----------------------dynamic arrays------------------------------------
INTEGER,ALLOCATABLE::g(:),g_g(:,:),g_num(:,:),kdiag(:),nf(:,:),no(:),
&
node(:),num(:),sense(:)
REAL(iwp),ALLOCATABLE::bdylds(:),bee(:,:),bload(:),coord(:,:),dee(:,:), &
der(:,:),deriv(:,:),devp(:),eld(:),eload(:),eps(:),erate(:), &
etensor(:,:,:),evp(:),evpt(:,:,:),flow(:,:),fun(:),gc(:),g_coord(:,:), &
jac(:,:),km(:,:),kv(:),loads(:),m1(:,:),m2(:,:),m3(:,:),oldis(:),
&
points(:,:),pore(:,:),sigma(:),storkv(:),stress(:),tensor(:,:,:),
&
totd(:),weights(:),x_coords(:),y_coords(:)
!-----------------------input and initialisation--------------------------
OPEN(10,FILE='fe95.dat'); OPEN(11,FILE='fe95.res')
READ(10,*)nxe,nye,phi,c,psi,e,v,bulk,cons
CALL mesh_size(element,nod,nels,nn,nxe,nye)
ALLOCATE(nf(nodof,nn),points(nip,ndim),weights(nip),g_coord(ndim,nn), &
x_coords(nxe+1),y_coords(nye+1),num(nod),evpt(nst,nip,nels), &
coord(nod,ndim),g_g(ndof,nels),tensor(nst,nip,nels),fun(nod), &
etensor(nst,nip,nels),dee(nst,nst),pore(nip,nels),stress(nst), &
jac(ndim,ndim),der(ndim,nod),deriv(ndim,nod),g_num(nod,nels), &
bee(nst,ndof),km(ndof,ndof),eld(ndof),eps(nst),sigma(nst),bload(ndof), &
eload(ndof),erate(nst),evp(nst),devp(nst),g(ndof),m1(nst,nst),
&
m2(nst,nst),m3(nst,nst),flow(nst,nst),gc(ndim))
READ(10,*)x_coords,y_coords
nf=1; READ(10,*)nr,(k,nf(:,k),i=1,nr); CALL formnf(nf); neq=MAXVAL(nf)
ALLOCATE(kdiag(neq),loads(0:neq),bdylds(0:neq),oldis(0:neq),totd(0:neq))
!-----------------------loop the elements to find global arrays sizes-----
kdiag=0
elements_1: DO iel=1,nels
CALL geom_rect(element,iel,x_coords,y_coords,coord,num,'y')
CALL num_to_g(num,nf,g); g_num(:,iel)=num
g_coord(:,num)=TRANSPOSE(coord); g_g(:,iel)=g; CALL fkdiag(kdiag,g)
