CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Main CFD Forum

A algorithm problem meets in the FTC3D code

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

Reply
 
LinkBack Thread Tools Display Modes
Old   September 14, 2015, 10:24
Default A algorithm problem meets in the FTC3D code
  #1
New Member
 
Liao Bin
Join Date: Jun 2012
Posts: 1
Rep Power: 0
liaoyun0126 is on a distinguished road
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

c generate gradient vector

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

c generate gradient vector

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
File Type: jpg density field picture.jpg (121.9 KB, 7 views)
liaoyun0126 is offline   Reply With Quote

Reply

Thread Tools
Display Modes

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 Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Fortran 90 code for bicgstab algorithm CuriousBoy Main CFD Forum 1 April 28, 2016 01:37
Unsteady flow code - Problem with space loop Hooman Main CFD Forum 9 August 10, 2011 17:05
Looking for a open FEM code to solve a FSI problem Esto Main CFD Forum 3 June 15, 2006 11:02
Problem convertng a code from 11 to 12 ali OpenFOAM Running, Solving & CFD 3 December 6, 2005 05:43
What kind of Cmmercial CFD code you feel well? Lans Main CFD Forum 13 October 27, 1998 11:20


All times are GMT -4. The time now is 11:12.