Material Non-linearity - Programming the Finite Element Method - page 302

eps=MATMUL(bee,eld); eps=eps-evpt(:,i,iel)
sigma=MATMUL(dee,eps); CALL invar(sigma,sigm,dsbar,lode_theta)
!-----------------------check whether yield is violated-------------------
CALL mocouf(phif,cf,sigm,dsbar,lode_theta,f); IF(f>fmax)fmax=f
IF(converged.OR.iters==limit)THEN; devp=sigma; ELSE
IF(f>=zero.OR.(converged.OR.iters==limit))THEN
CALL mocouq(psif,dsbar,lode_theta,dq1,dq2,dq3)
CALL formm(sigma,m1,m2,m3); flow=f*(m1*dq1+m2*dq2+m3*dq3)
erate=MATMUL(flow,sigma); evp=erate*dt
evpt(:,i,iel)=evpt(:,i,iel)+evp; devp=MATMUL(dee,evp)
END IF
END IF
IF(f>=zero.OR.(converged.OR.iters==limit))THEN
eload=MATMUL(devp,bee); bload=bload+eload*det*weights(i)
END IF
END DO gauss_points_2
!-----------------------compute the total bodyloads vector----------------
bdylds(g)=bdylds(g)+bload; bdylds(0)=zero
END DO elements_3
WRITE(*,'(A,F7.2,A,I4,A,F8.3)') &

" srf",srf(iy)," iteration",iters," F_max",fmax
IF(converged.OR.iters==limit)EXIT
END DO its; WRITE(11,'(F7.2,E12.4,I5)')srf(iy),MAXVAL(ABS(loads)),iters
IF(iters==limit)EXIT
END DO srf_trials
STOP
END PROGRAM p610

New scalar integers:
ifix sets boundary conditions,
(1 smooth-smooth, 2 rough-smooth, 3 rough-rough)

nze number of elements in z -direction

New scalar reals:
d1

depth of mesh in z -direction

This program demonstrates 3D plasticity analysis using 20-noded hexahedral elements
with “reduced” ( nip=8 ) integration. The example is of a simple 3D slope stability analysis,
and the program is very similar to its 2D counterpart Program 6.3. Only three additional
inputs are required as compared with the data for Program 6.3. The first is ifix which
fixes the front and back faces of the mesh (in the -direction) to either “rough” or “smooth”.

When ifix=1 , both boundaries are smooth, and if the slope is homogeneous the analysis
essentially reduces to plane strain, when ifix=2 the front face is rough and the back
face smooth, implying a line of symmetry along the centre of the embankment, and when
ifix=3 both boundaries are rough, enabling a full 3D analysis of a slope that may have
non-uniform and non-symmetric properties in the crest (

z

-) direction. The second new
input parameter is nze , which defines the number of slices of elements required in the
z

z

-direction.
Two new subroutines, emb 3d bc and emb 3d geom are introduced to generate,
respectively, the nodal freedom array nf , and the nodal coordinates and element node
numbering g coord and g num . The subroutines create a rather simple 3D geometry in

-direction, and the third is d1 which represents the depth of the slope in the z