Civil Engineering Reference
In-Depth Information
END DO elements_1; CALL mesh(g_coord,g_num,12)
DO i=2,neq; kdiag(i)=kdiag(i)+kdiag(i-1); END DO; ALLOCATE(kv(kdiag(neq)))
WRITE(11,'(2(A,I7))')
&
" There are",neq," equations and the skyline storage is",kdiag(neq)
!-----------------------add fluid bulk modulus effective dee matrix-------
CALL deemat(dee,e,v); pi=ACOS(-one)
DO i=1,nst; DO k=1,nst; IF(i/=3.AND.k/=3)dee(i,k)=dee(i,k)+bulk
END DO; END DO; snph=SIN(phi*pi/d180)
dt=d4*(one+ v)*(one-two*v)/(e*(one-two*v+snph*snph))
CALL sample(element,points,weights); kv=zero; tensor=zero; etensor=zero
!-----------------------element stiffness integration and assembly--------
elements_2: DO iel=1,nels
num=g_num(:,iel); coord=TRANSPOSE(g_coord(:,num)); g=g_g(:,iel); km=zero
gauss_pts_1: DO i=1,nip
CALL shape_fun(fun,points,i)
CALL bee8(bee,coord,points(i,1),points(i,2),det)
gc=MATMUL(fun,coord); bee(4,1:ndof-1:2)=fun(:)/gc(1)
km=km+MATMUL(MATMUL(TRANSPOSE(bee),dee),bee)*det*weights(i)*gc(1)
tensor(1:2,i,iel)=cons; tensor(4,i,iel)=cons
END DO gauss_pts_1; CALL fsparv(kv,km,g,kdiag)
END DO elements_2
!-----------------------read displacement data and factorise equations----
READ(10,*)fixed_freedoms
IF(fixed_freedoms/=0)THEN
ALLOCATE(node(fixed_freedoms),sense(fixed_freedoms),no(fixed_freedoms),&
storkv(fixed_freedoms))
READ(10,*)(node(i),sense(i),i=1,fixed_freedoms)
DO i=1,fixed_freedoms; no(i)=nf(sense(i),node(i)); END DO
kv(kdiag(no))=kv(kdiag(no))+penalty; storkv=kv(kdiag(no))
END IF; CALL sparin(kv,kdiag); CALL deemat(dee,e,v)
!-----------------------displacement increment loop-----------------------
READ(10,*)tol,limit,incs,presc
WRITE(11,'(/A)')" step
disp
dev stress pore press iters"
oldis=zero; totd=zero
disp_incs: DO iy=1,incs
iters=0; bdylds=zero; evpt=zero
!-----------------------plastic iteration loop----------------------------
its: DO
iters=iters+1
WRITE(*,'(A,E11.3,A,I4)')" displacement",iy*presc," iteration",iters
loads=zero; loads(no)=storkv*presc; loads=loads+bdylds
CALL spabac(kv,loads,kdiag)
!-----------------------check plastic convergence-------------------------
CALL checon(loads,oldis,tol,converged); IF(iters==1)converged=.FALSE.
!-----------------------go round the Gauss points ------------------------
elements_3: DO iel=1,nels
bload=zero; num=g_num(:,iel); coord=TRANSPOSE(g_coord(:,num))
g=g_g(:,iel); eld=loads(g)
gauss_pts_2: DO i=1,nip
CALL shape_fun(fun,points,i)
CALL bee8(bee,coord,points(i,1),points(i,2),det)
gc=MATMUL(fun,coord); bee(4,1:ndof-1:2)=fun(:)/gc(1)
eps=MATMUL(bee,eld); eps=eps-evpt(:,i,iel); sigma=MATMUL(dee,eps)
stress=sigma+tensor(:,i,iel)
CALL invar(stress,sigm,dsbar,lode_theta)
!-----------------------check whether yield is violated-------------------
CALL mocouf(phi,c,sigm,dsbar,lode_theta,f)
IF(f>=zero)THEN
