Civil Engineering Reference
In-Depth Information
loads=exc_loads+bdylds; CALL spabac(kv,loads,kdiag)
!-----------------------check convergence---------------------------------
CALL checon(loads,oldis,tol,converged)
IF(iters==1)converged=.FALSE.
IF(converged.OR.iters==limit)THEN; bdylds=zero
DO k=1,nn; DO i=1,nodof
IF(lnf(i,k)/=0)tot_d(i,k)=tot_d(i,k)+loads(lnf(i,k))
END DO; END DO
END IF
!-----------------------go round the Gauss Points-------------------------
elements_4: DO iel=1,nels
phi=prop(1,etype(iel)); c=prop(2,etype(iel))
psi=prop(3,etype(iel)); e=prop(5,etype(iel))
v=prop(6,etype(iel)); IF(solid(iel)==0)e=zero; bload=zero
CALL deemat(dee,e,v); num=g_num(:,iel); CALL num_to_g(num,lnf,g)
coord=TRANSPOSE(g_coord(:,num)); eld=loads(g)
gauss_pts_4: DO i=1,nip
CALL bee8(bee,coord,points(i,1),points(i,2),det)
eps=MATMUL(bee,eld); eps=eps-evpt(:,i,iel)
stress=tensor(:,i,iel)+MATMUL(dee,eps)
!-----------------------air element stresses are zero---------------------
IF(solid(iel)==0)stress=zero
CALL invar(stress,sigm,dsbar,lode_theta)
!-----------------------check whether yield is violated-------------------
CALL mocouf(phi,c,sigm,dsbar,lode_theta,f)
IF(converged.OR.iters==limit)THEN; devp=stress; ELSE
IF(f>=zero)THEN; CALL mocouq(psi,dsbar,lode_theta,dq1,dq2,dq3)
CALL formm(stress,m1,m2,m3); flow=f*(m1*dq1+m2*dq2+m3*dq3)
erate=MATMUL(flow,stress); 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
!-----------------------if appropriate update the Gauss point stresses----
IF(converged.OR.iters==limit)tensor(:,i,iel)=stress
END DO gauss_pts_4
!-----------------------compute the total bodyloads vector ---------------
bdylds(g)=bdylds(g)+bload; bdylds(0)=zero
END DO elements_4; IF(converged.OR.iters==limit)EXIT
END DO its
WRITE(11,'(A,I3,A,I5,A)')" Increment",iy," took",iters,
&
" iterations to converge"
IF(iy==incs.OR.iters==limit)THEN
WRITE(11,'(A)') " Node x-disp y-disp"
DO i=1,nouts; WRITE(11,'(I5,2E12.4)')no(i),tot_d(:,no(i)); END DO
EXIT
END IF
END DO load_incs; IF(ii==layers.OR.iters==limit)EXIT
DEALLOCATE(kdiag,kv,exc_loads,bdylds,oldis,loads,exele)
END DO layer_number
loads(lnf(1,:))=tot_d(1,:); loads(lnf(2,:))=tot_d(2,:)
g_num(:,totex(:ntote))=0; CALL mesh(g_coord,g_num,12)
CALL dismsh(loads,lnf,0.1_iwp,g_coord,g_num,13)
CALL vecmsh(loads,lnf,0.1_iwp,0.1_iwp,g_coord,g_num,14)
STOP
END PROGRAM p68
