Civil Engineering Reference
In-Depth Information
terms
in (6.33). In the viscoplastic algorithm however, only the potential function derivatives
{
∂Q/∂
σ
}
be approximated in this way due to the inclusion of both
{
∂Q/∂
σ
}
and
{
∂F/∂
σ
}
are not needed by equation (6.19).
All the programs in this chapter for solving two-dimensional problems have been based
on the 8-node quadrilateral element, together with reduced integration (four Gauss points
per element). This particular combination has been chosen for its simplicity, and also its
well-known ability to compute collapse loads accurately (e.g. Zienkiewicz
et al
., 1975;
Griffiths, 1980, 1982). Of course, other element types could be used if required, by making
similar changes to those described in Chapter 5. To alleviate “locking” problems with
lower order elements, “reduced integration” could be used selectively on the volumetric
components of the stiffness matrix (see e.g. Hughes, 1987; Griffiths and Mustoe, 1995).
are approximated since
{
∂F/∂
σ
}
Program 6.1 Plane strain bearing capacity analysis of an elastic-plastic (von Mises)
material using 8-node rectangular quadrilaterals. Viscoplastic strain method.
PROGRAM p61
!-------------------------------------------------------------------------
! Program 6.1 Plane strain bearing capacity analysis of an elastic-plastic
! (von Mises) material using 8-node rectangular
! quadrilaterals. Viscoplastic strain method.
!-------------------------------------------------------------------------
USE main; USE geom; IMPLICIT NONE
INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
INTEGER::i,iel,incs,iters,iy,k,limit,loaded_nodes,ndim=2,ndof=16,nels,
&
neq,nip=4,nn,nod=8,nodof=2,nprops=3,np_types,nr,nst=4,nxe,nye
REAL(iwp)::ddt,det,dq1,dq2,dq3,dsbar,dt=1.0e15_iwp,d3=3.0_iwp,d4=4.0_iwp,&
f,lode_theta,one=1.0_iwp,ptot,sigm,tol,two=2.0_iwp,zero=0.0_iwp
CHARACTER(LEN=15)::element='quadrilateral'; LOGICAL::converged
!-----------------------dynamic arrays------------------------------------
INTEGER,ALLOCATABLE::etype(:),g(:),g_g(:,:),g_num(:,:),kdiag(:),nf(:,:), &
node(:),num(:)
REAL(iwp),ALLOCATABLE::bdylds(:),bee(:,:),bload(:),coord(:,:),dee(:,:), &
der(:,:),deriv(:,:),devp(:),eld(:),eload(:),eps(:),erate(:),evp(:),
&
evpt(:,:,:),flow(:,:),g_coord(:,:),jac(:,:),km(:,:),kv(:),loads(:),
&
m1(:,:),m2(:,:),m3(:,:),oldis(:),points(:,:),prop(:,:),qinc(:),
&
sigma(:),stress(:),tensor(:,:,:),totd(:),val(:,:),weights(:),
&
x_coords(:),y_coords(:)
!-----------------------input and initialisation--------------------------
OPEN(10,FILE='fe95.dat'); OPEN(11,FILE='fe95.res')
READ(10,*)nxe,nye,np_types; 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),dee(nst,nst),g_g(ndof,nels), &
prop(nprops,np_types),etype(nels),evpt(nst,nip,nels),stress(nst), &
tensor(nst,nip,nels),coord(nod,ndim),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))
READ(10,*)prop; etype=1; IF(np_types>1)READ(10,*)etype
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')
