260 lines
8.1 KiB
Text
260 lines
8.1 KiB
Text
|
c 2D XFEM Corrosion Element
|
||
|
|
|
||
|
|
subroutine uel(rhs,amatrx,svars,energy,ndofel,nrhs,nsvars,props,
|
||
|
|
1 nprops,coords,mcrd,nnode,u,du,v,a,jtype,time,dtime,kstep,kinc,
|
||
|
|
2 jelem,params,ndload,jdltyp,adlmag,predef,npredf,lflags,
|
||
|
|
3 mlvarx,ddlmag,mdload,pnewdt,jprops,njprop,period)
|
||
|
|
c
|
||
|
|
include 'aba_param.inc'
|
||
|
|
c
|
||
|
|
dimension rhs(mlvarx,*),amatrx(ndofel,ndofel),svars(nsvars),
|
||
|
|
1 energy(8),props(*),coords(mcrd,nnode),
|
||
|
|
2 u(ndofel),du(mlvarx,*),v(ndofel),a(ndofel),time(2),
|
||
|
|
3 params(3),jdltyp(mdload,*),adlmag(mdload,*),
|
||
|
|
4 ddlmag(mdload,*),predef(2,npredf,nnode),lflags(*),jprops(*)
|
||
|
|
c
|
||
|
|
dimension phig(8),phih(8),phi(8),phix(8),phiy(8)
|
||
|
|
dimension crdnx(4),crdny(4),dndg(4),dndh(4)
|
||
|
|
dimension theta(4),rjac(2,2),rjaci(2,2)
|
||
|
|
dimension gx(100,4),hx(100,4),xi(2),yi(2),gi(2),hi(2)
|
||
|
|
dimension c(2),gp(2,8),gm2(8,8)
|
||
|
|
c
|
||
|
|
parameter(zero=0.d0,one=1.d0)
|
||
|
|
c material property definition
|
||
|
|
thick = 1.
|
||
|
|
rho = 1.
|
||
|
|
beta=0.
|
||
|
|
dpos=0.6
|
||
|
|
c initialization (nrhs=1)
|
||
|
|
do k1=1,ndofel
|
||
|
|
rhs(k1,nrhs)=zero
|
||
|
|
do k2=1,ndofel
|
||
|
|
amatrx(k2,k1)=zero
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
if (lflags(1).eq.33) then
|
||
|
|
do icrd=1,4
|
||
|
|
crdnx(icrd)=coords(1,icrd)
|
||
|
|
crdny(icrd)=coords(2,icrd)
|
||
|
|
theta(icrd)=abs(crdnx(icrd)-dpos)*
|
||
|
|
1 sign(1.,crdnx(icrd)-dpos)
|
||
|
|
enddo
|
||
|
|
c if (sign(1.,theta(1))/=sign(1.,theta(2)))then
|
||
|
|
if (2==1)then
|
||
|
|
c possible enriched element
|
||
|
|
npart=10
|
||
|
|
rpart=float(npart)
|
||
|
|
ienr=npart*npart
|
||
|
|
do isdx=1,npart
|
||
|
|
do isdy=1,npart
|
||
|
|
rmidx=-1.-1./rpart+(2./rpart)*float(isdx)
|
||
|
|
rmidy=-1.-1./rpart+(2./rpart)*float(isdy)
|
||
|
|
isubindex=npart*(isdy-1)+isdx
|
||
|
|
gpos=1./(sqrt(3.)*rpart)
|
||
|
|
gx(isubindex,1)=rmidx-gpos
|
||
|
|
gx(isubindex,2)=rmidx+gpos
|
||
|
|
gx(isubindex,3)=rmidx+gpos
|
||
|
|
gx(isubindex,4)=rmidx-gpos
|
||
|
|
hx(isubindex,1)=rmidy-gpos
|
||
|
|
hx(isubindex,2)=rmidy-gpos
|
||
|
|
hx(isubindex,3)=rmidy+gpos
|
||
|
|
hx(isubindex,4)=rmidy+gpos
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
c check if int points are on different sides of front
|
||
|
|
icheck=0
|
||
|
|
do i=1,ienr
|
||
|
|
do j=1,4
|
||
|
|
g=gx(i,j)
|
||
|
|
h=hx(i,j)
|
||
|
|
phi(1)=0.25*(1.-g)*(1.-h)
|
||
|
|
phi(3)=0.25*(1.+g)*(1.-h)
|
||
|
|
phi(5)=0.25*(1.+g)*(1.+h)
|
||
|
|
phi(7)=0.25*(1.-g)*(1.+h)
|
||
|
|
rLS=theta(1)*phi(1)+theta(2)*phi(3)
|
||
|
|
1 +theta(3)*phi(5)+theta(4)*phi(7)
|
||
|
|
if (i==1 .and. j==1)then
|
||
|
|
sgn=sign(1.,rLS)
|
||
|
|
else
|
||
|
|
if (sign(1.,rLS)/=sgn)then
|
||
|
|
icheck=1
|
||
|
|
endif
|
||
|
|
endif
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
if (icheck==0)then
|
||
|
|
c regular element - fix extra dofs
|
||
|
|
ienr=1
|
||
|
|
gpos=1/sqrt(3.)
|
||
|
|
gx(1,1)=-gpos
|
||
|
|
gx(1,2)=gpos
|
||
|
|
gx(1,3)=gpos
|
||
|
|
gx(1,4)=-gpos
|
||
|
|
hx(1,1)=-gpos
|
||
|
|
hx(1,2)=-gpos
|
||
|
|
hx(1,3)=gpos
|
||
|
|
hx(1,4)=gpos
|
||
|
|
endif
|
||
|
|
else
|
||
|
|
c Normal Shp Funs
|
||
|
|
ienr=1
|
||
|
|
gpos=1/sqrt(3.)
|
||
|
|
gx(1,1)=-gpos
|
||
|
|
gx(1,2)=gpos
|
||
|
|
gx(1,3)=gpos
|
||
|
|
gx(1,4)=-gpos
|
||
|
|
hx(1,1)=-gpos
|
||
|
|
hx(1,2)=-gpos
|
||
|
|
hx(1,3)=gpos
|
||
|
|
hx(1,4)=gpos
|
||
|
|
endif
|
||
|
|
c assemble amatrx and rhs
|
||
|
|
do i=1,ienr
|
||
|
|
do j=1,4
|
||
|
|
g=gx(i,j)
|
||
|
|
h=hx(i,j)
|
||
|
|
phi(1)=0.25*(1.-g)*(1.-h)
|
||
|
|
phi(3)=0.25*(1.+g)*(1.-h)
|
||
|
|
phi(5)=0.25*(1.+g)*(1.+h)
|
||
|
|
phi(7)=0.25*(1.-g)*(1.+h)
|
||
|
|
rLS=theta(1)*phi(1)+theta(2)*phi(3)
|
||
|
|
1 +theta(3)*phi(5)+theta(4)*phi(7)
|
||
|
|
cond=1.
|
||
|
|
spec=1.
|
||
|
|
do iter=1,4
|
||
|
|
phi(2*iter)=phi(2*iter-1)*
|
||
|
|
1 (abs(rLS)-abs(theta(iter)))
|
||
|
|
enddo
|
||
|
|
phig(1)=0.25*-(1.-h)
|
||
|
|
phig(3)=0.25*(1.-h)
|
||
|
|
phig(5)=0.25*(1.+h)
|
||
|
|
phig(7)=0.25*-(1.+h)
|
||
|
|
phih(1)=0.25*-(1.-g)
|
||
|
|
phih(3)=0.25*-(1.+g)
|
||
|
|
phih(5)=0.25*(1.+g)
|
||
|
|
phih(7)=0.25*(1.-g)
|
||
|
|
diLSg=sign(1.,rLS)*(phig(1)*theta(1)+phig(3)*
|
||
|
|
1 theta(2)+phig(5)*theta(3)+phig(7)*theta(4))
|
||
|
|
diLSh=sign(1.,rLS)*(phih(1)*theta(1)+phih(3)*
|
||
|
|
1 theta(2)+phih(5)*theta(3)+phih(7)*theta(4))
|
||
|
|
do iter=1,4
|
||
|
|
phig(2*iter)=phig(2*iter-1)*(abs(rLS)-
|
||
|
|
1 abs(theta(iter)))+phi(2*iter-1)*diLSg
|
||
|
|
phih(2*iter)=phih(2*iter-1)*(abs(rLS)-
|
||
|
|
1 abs(theta(iter)))+phi(2*iter-1)*diLSh
|
||
|
|
enddo
|
||
|
|
rjac=zero
|
||
|
|
do iter=1,4
|
||
|
|
rjac(1,1)=rjac(1,1)+phig(2*iter-1)*crdnx(iter)
|
||
|
|
rjac(1,2)=rjac(1,2)+phig(2*iter-1)*crdny(iter)
|
||
|
|
rjac(2,1)=rjac(2,1)+phih(2*iter-1)*crdnx(iter)
|
||
|
|
rjac(2,2)=rjac(2,2)+phih(2*iter-1)*crdny(iter)
|
||
|
|
enddo
|
||
|
|
djac=rjac(1,1)*rjac(2,2)-rjac(1,2)*rjac(2,1)
|
||
|
|
rjaci(1,1)= rjac(2,2)/djac
|
||
|
|
rjaci(2,2)= rjac(1,1)/djac
|
||
|
|
rjaci(1,2)=-rjac(1,2)/djac
|
||
|
|
rjaci(2,1)=-rjac(2,1)/djac
|
||
|
|
do iter=1,8
|
||
|
|
phix(iter)=rjaci(1,1)*phig(iter)+
|
||
|
|
1 rjaci(1,2)*phih(iter)
|
||
|
|
phiy(iter)=rjaci(2,1)*phig(iter)+
|
||
|
|
1 rjaci(2,2)*phih(iter)
|
||
|
|
enddo
|
||
|
|
dtdx=zero
|
||
|
|
dtdy=zero
|
||
|
|
t =zero
|
||
|
|
told=zero
|
||
|
|
do iter=1,8
|
||
|
|
dtdx=u(iter)*phix(iter)+dtdx
|
||
|
|
dtdy=u(iter)*phiy(iter)+dtdy
|
||
|
|
t=u(iter)*phi(iter)+t
|
||
|
|
told=(u(iter)-du(iter,nrhs))*phi(iter)+told
|
||
|
|
end do
|
||
|
|
dtdt=(t-told)/dtime
|
||
|
|
we=djac
|
||
|
|
do ki=1,8
|
||
|
|
c loop over nodes
|
||
|
|
rhs(ki,nrhs) = rhs(ki,nrhs) -
|
||
|
|
1 (we/float(ienr))*(phi(ki)*rho*spec*dtdt+
|
||
|
|
2 cond*(phix(ki)*dtdx + phiy(ki)*dtdy))
|
||
|
|
do kj=1,8
|
||
|
|
amatrx(ki,kj)=amatrx(ki,kj)+(we/float(ienr))
|
||
|
|
1 *(phi(ki)*phi(kj)*rho*spec/dtime+cond
|
||
|
|
2 *(phix(ki)*phix(kj)+phiy(ki)*phiy(kj)))
|
||
|
|
end do
|
||
|
|
end do
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
c if interface is in the element a penalty term is needed
|
||
|
|
if(ienr>1)then
|
||
|
|
icount=0
|
||
|
|
if (sign(1.,theta(1))/=sign(1.,theta(2)))then
|
||
|
|
icount=icount+1
|
||
|
|
f=abs(theta(1))/(abs(theta(1))+abs(theta(2)))
|
||
|
|
xi(icount)=f*(crdnx(2)-crdnx(1))+crdnx(1)
|
||
|
|
yi(icount)=f*(crdny(2)-crdny(1))+crdny(1)
|
||
|
|
gi(icount)=(2.*xi(icount)-(crdnx(1)+crdnx(2)))
|
||
|
|
1 /(-crdnx(1)+crdnx(2))
|
||
|
|
hi(icount)=-1.
|
||
|
|
endif
|
||
|
|
if (sign(1.,theta(2))/=sign(1.,theta(3)))then
|
||
|
|
icount=icount+1
|
||
|
|
f=abs(theta(2))/(abs(theta(2))+abs(theta(3)))
|
||
|
|
xi(icount)=f*(crdnx(3)-crdnx(2))+crdnx(2)
|
||
|
|
yi(icount)=f*(crdny(3)-crdny(2))+crdny(2)
|
||
|
|
gi(icount)=1.
|
||
|
|
hi(icount)=(2.*yi(icount)-(crdny(2)+crdny(3)))
|
||
|
|
1 /(-crdny(2)+crdny(3))
|
||
|
|
endif
|
||
|
|
if (sign(1.,theta(3))/=sign(1.,theta(4)))then
|
||
|
|
icount=icount+1
|
||
|
|
f=abs(theta(3))/(abs(theta(3))+abs(theta(4)))
|
||
|
|
xi(icount)=f*(crdnx(4)-crdnx(3))+crdnx(3)
|
||
|
|
yi(icount)=f*(crdny(4)-crdny(3))+crdny(3)
|
||
|
|
gi(icount)=(2.*xi(icount)-(crdnx(4)+crdnx(3)))
|
||
|
|
1 /(-crdnx(4)+crdnx(3))
|
||
|
|
hi(icount)=1.
|
||
|
|
endif
|
||
|
|
if (sign(1.,theta(1))/=sign(1.,theta(4)))then
|
||
|
|
icount=icount+1
|
||
|
|
f=abs(theta(1))/(abs(theta(1))+abs(theta(4)))
|
||
|
|
xi(icount)=f*(crdnx(4)-crdnx(1))+crdnx(1)
|
||
|
|
yi(icount)=f*(crdny(4)-crdny(1))+crdny(1)
|
||
|
|
gi(icount)=-1.
|
||
|
|
hi(icount)=(2.*yi(icount)-(crdny(1)+crdny(4)))
|
||
|
|
1 /(-crdny(4)+crdny(1))
|
||
|
|
endif
|
||
|
|
c(1)=1.
|
||
|
|
c(2)=1.
|
||
|
|
do iter=1,2
|
||
|
|
Gp(iter,1)=0.25*(1.-gi(iter))*(1.-hi(iter))
|
||
|
|
Gp(iter,3)=0.25*(1.+gi(iter))*(1.-hi(iter))
|
||
|
|
Gp(iter,5)=0.25*(1.+gi(iter))*(1.+hi(iter))
|
||
|
|
Gp(iter,7)=0.25*(1.-gi(iter))*(1.+hi(iter))
|
||
|
|
Gp(iter,2)=-Gp(iter,1)*abs(theta(1))
|
||
|
|
Gp(iter,4)=-Gp(iter,3)*abs(theta(2))
|
||
|
|
Gp(iter,6)=-Gp(iter,5)*abs(theta(3))
|
||
|
|
Gp(iter,8)=-Gp(iter,7)*abs(theta(4))
|
||
|
|
enddo
|
||
|
|
do i=1,8
|
||
|
|
rhs(i,nrhs)=rhs(i,nrhs)+
|
||
|
|
1 beta*(Gp(1,i)*c(1)+Gp(2,i)*c(2))
|
||
|
|
enddo
|
||
|
|
c find GtG
|
||
|
|
gm2=0.
|
||
|
|
do i=1,8
|
||
|
|
do j=1,8
|
||
|
|
gm2(i,j)=gp(1,i)*gp(1,j)+gp(2,i)*gp(2,j)
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
c add penalty stiffness
|
||
|
|
do i=1,8
|
||
|
|
do j=1,8
|
||
|
|
amatrx(i,j)=amatrx(i,j)+beta*gm2(i,j)
|
||
|
|
enddo
|
||
|
|
enddo
|
||
|
|
endif
|
||
|
|
end if
|
||
|
|
return
|
||
|
|
end
|