# A algorithm problem meets in the FTC3D code

 Register Blogs Members List Search Today's Posts Mark Forums Read

September 14, 2015, 10:24
A algorithm problem meets in the FTC3D code
#1
New Member

Liao Bin
Join Date: Jun 2012
Posts: 1
Rep Power: 0
Hi all,
Recently,I download the FTC(Front-Tracking code)3D code from the internet. But the FTC3D had been modified already. Specially, I think the subroutine about updating the density field is not exact. Then I found the FTC2D code from Dr. Tryggvason's personal web:http://www3.nd.edu/~gtryggva/. And I rewrite the subroutine of the FTC3D record to the algorithm of the FTC2D, but I meet a problem!
The subroutine in the FTC2D is:
subroutine dens2(r,rtmp,tmp1,tmp2,nxp2,nyp2,iii,idim,pt,ptcon ,icp,elcon,
elprop,maxpt,maxel)

dimension tmp1(nxp2+2,nyp2+2), tmp2(nxp2+2,nyp2+2)
dimension rtmp(nxp2,nyp2)
dimension r(nxp2,nyp2)
dimension iii(idim,2)

dimension pt(maxpt,2), ptcon(maxpt)
dimension icp(maxel,2),elcon(maxel),elprop(maxel,3)

integer ptcon,elcon,np,ffp,lfp,fep, ne,ffe,lfe,fee

common/grid/hxi,hyi,nxp1,nyp1,nx,ny,xl,yl
common/bounds/ibdry,jbdry,nxlast,nylast
common/front/np,ffp,lfp,fep, ne,ffe,lfe,fee,amax,amin
common/fparam/xmv,mv,fxl,fyl
common/flprop/r1,r2,m1,m2,gx,gy,rro

xx=(xmv+0.5)*xl
yy=(xmv+0.5)*yl

sc=1./2.
pi=4.*ATAN(1.)
pd2=0.5*pi
cnstx=(.25*hxi)**2
cnsty=(.25*hyi)**2

do 10 i=1,nxp2+2
do 10 j=1,nyp2+2
tmp1(i,j)=0.0
tmp2(i,j)=0.0
10 continue

k=ffe

do 30 ll=1,ne

p0x=pt(icp(k,1),1)
p0y=pt(icp(k,1),2)

p1x= p0x+amod((pt(icp(k,2),1)-p0x+xx),xl)-0.5*xl
p1y= p0y+amod((pt(icp(k,2),2)-p0y+yy),yl)-0.5*yl

xc=0.5*(p1x+p0x)
yc=0.5*(p1y+p0y)

x1=p1x-p0x
y1=p1y-p0y

c calculate components of n.dA

rx=-y1 *elprop(k,1)
ry= x1 *elprop(k,1)

c generate a grid value

xt=amod((xc+0.5/hxi)+xmv*xl,xl)
yt=amod((yc+0.5/hyi)+xmv*yl,yl)

ir=1+int(nx*xt/xl)
jr=1+int(ny*yt/yl)

do 40 i1=1,4
do 40 j1=1,4

ii=ir-2+i1
jj=jr-2+j1

drx = 1. + cos((xt*hxi - float(ii-1))*pd2)
dry = 1. + cos((yt*hyi - float(jj-1))*pd2)

tmp1(ii+1,jj+1)=tmp1(ii+1,jj+1)+rx*drx*dry*cnstx
tmp2(ii+1,jj+1)=tmp2(ii+1,jj+1)+ry*drx*dry*cnsty

40 continue

k=elcon(k)
30 continue

c correct boundaries:

call gbd1(tmp1,tmp2,ibdry,jbdry,nxp2,nyp2)

do 900 i=1,nxp2
do 900 j=1,nyp2
900 rtmp(i,j)=0.0

do 80 i=2,nxp1
do 80 j=2,nyp1
rtmp(i,j)= 0.5*(hxi*(tmp1(i+1+1,j+1)-tmp1(i-1+1,j+1))+
1 hyi*(tmp2(i+1,j+1+1)-tmp2(i+1,j-1+1)))
80 continue

c set connection
node=0
do 100 i=2,nxp1
do 100 j=2,nyp1
cc=0.0
do 105 io=1,3
do 105 jo=1,3
105 cc=cc+rtmp(i+io-2,j+jo-2)
if(abs(cc) .gt. 1.0e-03)then
node=node+1
iii(node,1)=i
iii(node,2)=j
end if
c open(13,file='cc',status='new')
c write(13,*) cc,i,j,node
100 continue
c pause

numnd=node

c then iterate

c-----------the following should be read in
itmax=1000
bet=1.2
errmax=1.0e-05
c------------------------------------------
cf=0.5/(hxi*hxi+hyi*hyi)
do 200 iter=1,itmax

err=0.0
do 210 node=1,numnd
i=iii(node,1)
j=iii(node,2)
c do 210 i=2,nxp1
c do 210 j=2,nyp1
rold=r(i,j)
r(i,j)=(1.0-bet)*r(i,j)+bet*cf*
1 (hxi*hxi*(r(i+1,j)+r(i-1,j))+
2 hyi*hyi*(r(i,j+1)+r(i,j-1))-rtmp(i,j))
err=err+abs(rold-r(i,j))
210 continue

c set boundary terms...
if(ibdry .eq. 0)then
do 110 j=1,nyp2
r(nxp2,j)=r(2,j)
r(1,j)=r(nxp1,j)
110 continue
else
do 115 j=1,nyp2
r(nxp2,j)=r(nxp1,j)
r(1,j)=r(2,j)
115 continue
end if

if(jbdry .eq. 0)then
do 120 i=1,nxp2
r(i,nyp2)=r(i,2)
r(i,1)=r(i,nyp1)
120 continue
else
do 125 i=1,nxp2
r(i,nyp2)=r(i,nyp1)
r(i,1)=r(i,2)
125 continue
end if

c check convergence......

if(err.lt.(errmax*float(numnd)))go to 300
200 continue
300 continue

write(*,*)' DENSITY: ',err,' after ',iter,' iterations',
1 ' using ',numnd,' nodes'

return
end

And I rewrite the 3D subroutine:
subroutine creater(r,rtmp,tmp1,tmp2,tmp3,nxp2,nyp2,nzp2,iii,i dim,pt,ptcon,
icp,elcon,maxpt,maxel)

dimension tmp1(nxp2,nyp2,nzp2),tmp2(nxp2,nyp2,nzp2)
dimension tmp3(nxp2,nyp2,nzp2),rtmp(nxp2,nyp2,nzp2)
dimension r(nxp2,nyp2,nzp2)
dimension iii(idim,3)

dimension pt(maxpt,3), ptcon(maxpt)
dimension icp(maxel,3),elcon(maxel)

dimension work(10000),bd(1,1)

integer ptcon,elcon,np,ffp,lfp,fep, ne,ffe,lfe,fee
integer p1,p2,p3

common/grid/hxi,hyi,hzi,nxp1,nyp1,nzp1,nx,ny,nz,xl,yl,zl
common/front/np,ffp,lfp,fep, ne,ffe,lfe,fee,amax,amin,aspmax
common/flprop/r1,r2,m1,m2,gx,gy,gz,st,rro

sc=1./3.
pd2=0.5*pimach(1.0)

cnstx=(.25*hxi)**3
cnsty=(.25*hyi)**3
cnstz=(.25*hzi)**3

dxh=0.5/hxi
dyh=0.5/hyi
dzh=0.5/hzi

do 10 i=1,nxp2
do 10 j=1,nyp2
do 10 k=1,nzp2

tmp1(i,j,k)=0.0
tmp2(i,j,k)=0.0
tmp3(i,j,k)=0.0

10 continue

k=ffe

do 30 ll=1,ne

p1=icp(k,1)
p2=icp(k,2)
p3=icp(k,3)

xc=(pt(p1,1)+pt(p2,1)+pt(p3,1))*sc
yc=(pt(p1,2)+pt(p2,2)+pt(p3,2))*sc
zc=(pt(p1,3)+pt(p2,3)+pt(p3,3))*sc

x1=pt(p2,1)-pt(p1,1)
y1=pt(p2,2)-pt(p1,2)
z1=pt(p2,3)-pt(p1,3)

x2=pt(p3,1)-pt(p1,1)
y2=pt(p3,2)-pt(p1,2)
z2=pt(p3,3)-pt(p1,3)

c calculate components of n.dA

rx=(y1*z2-y2*z1)
ry=(z1*x2-z2*x1)
rz=(x1*y2-x2*y1)

c generate a grid value

xt=amod((xc+dxh)+5.0*xl,xl)
yt=amod((yc+dyh)+5.0*yl,yl)
zt=amod((zc+dzh)+5.0*zl,zl)

ir=1+int(nx*xt/xl)
jr=1+int(ny*yt/yl)
kr=1+int(nz*zt/zl)

do 40 i1=1,4
do 40 j1=1,4
do 40 k1=1,4

ii=ir-2+i1
jj=jr-2+j1
kk=kr-2+k1

drx = 1. + cos((xt*hxi - float(ii-1))*pd2)
dry = 1. + cos((yt*hyi - float(jj-1))*pd2)
drz = 1. + cos((zt*hzi - float(kk-1))*pd2)

ii=mod((ii-1+5*nx),nx)+1
jj=mod((jj-1+5*ny),ny)+1
kk=mod((kk-1+5*nz),nz)+1

tmp1(ii,jj,kk)=tmp1(ii,jj,kk)+rx*drx*dry*drz*cnstx
tmp2(ii,jj,kk)=tmp2(ii,jj,kk)+ry*drx*dry*drz*cnsty
tmp3(ii,jj,kk)=tmp3(ii,jj,kk)+rz*drx*dry*drz*cnstz

40 continue

k=elcon(k)

30 continue

c correct the ends of the box

call gbdry2(tmp1,nxp2,nyp2,nzp2)
call gbdry2(tmp2,nxp2,nyp2,nzp2)
call gbdry2(tmp3,nxp2,nyp2,nzp2)

do 900 i=1,nxp2
do 900 j=1,nyp2
do 900 k=1,nzp2
rtmp(i,j,k)=0.0
900 continue

do 80 i=2,nxp1
do 80 j=2,nyp1
do 80 k=2,nzp1

rtmp(i,j,k)=0.5*(hxi*hxi*(tmp1(i+1,j,k)-tmp1(i-1,j,k))+
1 hyi*hyi*(tmp2(i,j+1,k)-tmp2(i,j-1,k))+
2 hzi*hzi*(tmp3(i,j,k+1)-tmp3(i,j,k-1)))

80 continue

node=0

do 100 i=2,nxp1
do 100 j=2,nyp1
do 100 k=2,nzp1

cc=0.0

do 105 io=1,3
do 105 jo=1,3
do 105 ko=1,3

cc=cc+rtmp(i+io-2,j+jo-2,k+ko-2)

105 continue

if(abs(cc) .gt. 1.0e-3)then

node=node+1
iii(node,1)=i
iii(node,2)=j
iii(node,3)=k

end if

100 continue

numnd=node

itmax=1000
errmax=1.0e-3
beta=1.2

cf1=0.5/(hxi*hxi*hxi+hyi*hyi*hyi+hzi*hzi*hzi)

do 200 iter=1,itmax

err1=0.0

do 210 node=1,numnd

i=iii(node,1)
j=iii(node,2)
k=iii(node,3)

rold=r(i,j,k)

r(i,j,k)=(1.0-beta)*r(i,j,k)+beta*cf1*
1 (hxi*hxi*hxi*(r(i+1,j,k)+r(i-1,j,k))+
2 hyi*hyi*hyi*(r(i,j+1,k)+r(i,j-1,k))+
3 hzi*hzi*hzi*(r(i,j,k+1)+r(i,j,k-1))-rtmp(i,j,k))

err1=err1+abs(rold-r(i,j,k))

210 continue

call gbdry(r,nxp2,nyp2,nzp2)

if(err1.lt.errmax) go to 300

200 continue
300 continue

write(*,*)' DENSITY: ',err1,' after ',iter,' iterations',
1 ' using ',numnd,' nodes'

return
end

The results of the 2D and 3D code are in the attachment.
I readed the article of Dr.Tryggvason "A Front-Tracking Method for the Computations of Multiphase Flow'', the section 3.5 "Updating the Material Properties" explain the algorithm of Updating the density field. But I can't understand. So I need some help to modify the code. Thanks!
Yours'
Liaobin
Attached Images
 density field picture.jpg (121.9 KB, 12 views)

 Thread Tools Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post CuriousBoy Main CFD Forum 1 April 28, 2016 01:37 Hooman Main CFD Forum 9 August 10, 2011 17:05 Esto Main CFD Forum 3 June 15, 2006 11:02 ali OpenFOAM Running, Solving & CFD 3 December 6, 2005 05:43 Lans Main CFD Forum 13 October 27, 1998 11:20

All times are GMT -4. The time now is 22:27.