Civil Engineering Reference
In-Depth Information
ress=stress-(tensor(:,i,iel)+sigma-caflow)
CALL fmacat(vmfl,acat); acat=acat/dsbar
acatc=MATMUL(dee,acat); qmat=acatc*dlam
DO k=1,4; qmat(k,k)=qmat(k,k)+one; END DO; CALL invert(qmat)
vmtemp(1,:)=vmfl; vmtemp=MATMUL(vmtemp,qmat)
vmflq=vmtemp(1,:); top=DOT_PRODUCT(vmflq,ress)
vmtemp=MATMUL(vmtemp,dee); vmfla=vmtemp(1,:)
bot=DOT_PRODUCT(vmfla,vmfl); dslam=(fnew-top)/bot
qinvr=MATMUL(qmat,ress); qinva=MATMUL(MATMUL(qmat,dee),vmfl)
dsigma=-qinvr-qinva*dslam; stress=stress+dsigma
CALL invar(stress,sigm,dsbar,lode_theta)
fnew=dsbar-SQRT(d3)*prop(1,etype(iel)); dlam=dlam+dslam
IF(fnew<tol)EXIT
END DO iterate_on_fnew
dl(i,iel)=dlam; elso=tensor(:,i,iel)+sigma-stress
eload=MATMUL(elso,bee); bload=bload+eload*det*weights(i)
CALL vmflow(stress,dsbar,vmfl)
CALL fmrmat(vmfl,dsbar,dlam,dee,rmat); caflow=MATMUL(rmat,vmfl)
bot=DOT_PRODUCT(vmfl,caflow); CALL formaa(vmfl,rmat,daatd)
dee=rmat-daatd/bot
END IF
IF(fstiff<zero)
&
CALL deemat(dee,prop(2,etype(iel)),prop(3,etype(iel)))
km=km+MATMUL(MATMUL(TRANSPOSE(bee),dee),bee)*det*weights(i)
!-----------------------update the Gauss Point stresses-------------------
tensor(:,i,iel)=tensor(:,i,iel)+sigma-elso; stress=tensor(:,i,iel)
eload=MATMUL(stress,bee); dload=dload+eload*det*weights(i)
END DO gauss_points_2
!-----------------------compute the total bodyloads vector----------------
bdylds(g)=bdylds(g)+bload; bdylds(0)=zero; ddylds(g)=ddylds(g)+dload
ddylds(0)=zero; storkm(:,:,iel)=km
DO k=1,ndof; diag_precon(g(k))=diag_precon(g(k))+km(k,k); END DO
END DO elements_4; diag_precon(1:neq)=one/diag_precon(1:neq)
diag_precon(0)=zero; tloads=SUM(bdylds); IF(iters==1)converged=.FALSE.
IF(iters/=1.AND.tloads<ltol)converged=.TRUE.; totd=totd+loads
IF(converged.OR.iters==limit)EXIT
END DO plastic_iters; totd=totd+loads
WRITE(11,'(I5,2E12.4,I5,F17.2)')
&
iy,ptot,totd(nf(2,no(1))),iters,REAL(cg_tot)/REAL(iters)
IF(iters==limit)EXIT
END DO load_increments
CALL dismsh(totd,nf,0.05_iwp,g_coord,g_num,13)
CALL vecmsh(totd,nf,0.05_iwp,0.1_iwp,g_coord,g_num,14)
STOP
END PROGRAM p66
The “mesh-free” version of the constant stiffness method, Program 6.2, was inefficient
(at least on a scalar computer) due to the large number of repeated equation solutions.
Program 6.5 has shown that a tangent stiffness approach, with consistent return, has reduced
the number of equation solutions to 4 per load increment. With significantly fewer equations
to be solved, there is less benefit to be gained from factorisation, implying that a “mesh-free”
approach using an iterative preconditioned conjugate gradient solver could be appropriate.
Program 6.6 should be seen as a merging of Programs 6.2 and 6.5. The data for the problem
solved, again the original one of Figure 6.9, are shown in Figure 6.27. Extra information,
