[dinhanh]

Hi everyone,

I tried to make the c-grid over NACA foil by using the Poisson eq.

First, the distribution of grid point on boundary is specified and fixed by stretching function.

Next, internal grid point is specified by tranfinte interpolation.

Next, the control function in boudary is calculated based on the the fixed grid at boundary and of course they are fixed value.

The control function in internal point is calculated based on the linear interpolation of the control function at boundary as paper of Thompson.

Then the Poisson eq. is solved

However, the grid that I did con not force the grid close to the foil wall. I dont know why it is happened.

In here, the grid point is stretched to the nuy=0 and xi=0.5 (nuy and xi is the diection in computation coordinate). But it seem to large grid space from the wall as in the pictures

Anyone has experiency on this proplem please help me??

Thank you.

Below is my code:

!c
!C************************************************ ***************
!C************************************************ ********************
PARAMETER(ni=301,nj=70)
IMPLICIT REAL*8 (a-h,k,l,o-z)
!c
dimension x(ni,nj),y(ni,nj),&
x1(ni,nj),y1(ni,nj),&
phi(ni,nj),pha(ni,nj)
!c
pai=4.*atan(1.)
! Naca 4-rigid constant
tf = 0.15
ac0= 0.2969
ac1=-0.126
ac2=-0.3516
ac3= 0.2843
ac4=-0.1015
c = 30./1000
! Number of point in each domain
ni1= 101
ni2= 51
! Domain Boundary Length
ac= 5.*c
bc= 3.*c
oc= 5.*c
! Stretch constant
alf = 5.
alfy= 5.
! Boundary in X-direction
! domain 2
do i=1,ni2
dxi=1./(ni2-1)
xi =dxi*(i-1)
a1=exp(alf*xi)-1 !!stretching grid at boudary
a2=exp(alf)-1.
x(ni1+i-1, 1)=oc*a1/a2+ac
y(ni1+i-1, 1)=0.
x(ni1+i-1,nj)=x(ni1+i-1,1)
y(ni1+i-1,nj)=bc
end do
! domain 1
do i=1,ni1
dxi=1./(ni1-1)
xi =dxi*(i-1)
xi0=0.0
a1=exp(alf*xi)-1
a2=exp(alf)-1.
xc=c*a1/a2
y(i,1)=c*tf/0.2*(ac0*(xc/c)**0.5+ac1*(xc/c)+ac2*(xc/c)**2&
+ac3*(xc/c)**3+ac4*(xc/c)**4)
x(i,1)=xc+ac-c
!
a1=exp(alf*xi)-1
a2=exp(alf)-1.
xc=ac*a1/a2
x(i,nj)=xc
y(i,nj)=bc*sqrt(1.-((ac-xc)/ac)**2)
enddo
!
do i=1,ni1+ni2-1
x1(i, 1)=x(ni1+ni2-1-i+1, 1)
y1(i, 1)=y(ni1+ni2-1-i+1, 1)
x1(i,nj)=x(ni1+ni2-1-i+1,nj)
y1(i,nj)=y(ni1+ni2-1-i+1,nj)
end do
! Full domain in X-direction
ini=ni1+ni2-1
do i=1,ini
x1(ini+i-1, 1)= x(i, 1)
y1(ini+i-1, 1)=-y(i, 1)
x1(ini+i-1,nj)= x(i,nj)
y1(ini+i-1,nj)=-y(i,nj)
end do
do i=1,ni
x(i, 1)=x1(i, 1)
y(i, 1)=y1(i, 1)
x(i,nj)=x1(i,nj)
y(i,nj)=y1(i,nj)
end do
! boundary eta direction
do j=1,nj
eta=1./(nj-1)*(j-1)
a1=exp(alfy*eta)-1.
a2=exp(alfy)-1.
x( 1,j)=ac+oc
x(ni,j)=ac+oc
y( 1,j)=bc*a1/a2
y(ni,j)=-y(1,j)
end do
! Tran Finite Inter.
do j=2,nj-1
do i=2,ni-1
eta=1./(nj-1)*(j-1)
xi =1./(ni-1)*(i-1)
x(i,j)=(1.-xi)*x(1,j)+xi*x(ni,j)+(1.-eta)*x(i,1)+eta*x(i,nj)&
-(1.-xi)*(1.-eta)*x(1,1)-(1.-xi)*eta*x(1,nj)&
-xi*(1.-eta)*x(ni,1)-xi*eta*x(ni,nj)
y(i,j)=(1.-xi)*y(1,j)+xi*y(ni,j)+(1.-eta)*y(i,1)+eta*y(i,nj)&
-(1.-xi)*(1.-eta)*y(1,1)-(1.-xi)*eta*y(1,nj)&
-xi*(1.-eta)*y(ni,1)-xi*eta*y(ni,nj)
end do
end do
!C ******************** ELIPTIC GRID ******************
ca1=0.06
ca2=-1000.
cc1=2.
cc2=10.36
dxi=1./(ni-1)
det=1./(nj-1)
xtemp=0.
ytemp=0.
omega=0.27
residual=1.

do i=2,ni-1
xxi=.5*(x(i+1,1)-x(i-1,1))/dxi
yxi=.5*(y(i+1,1)-y(i-1,1))/dxi
xxxi=(x(i+1,1)-2.*x(i,1)+x(i-1,1))/dxi**2
yxxi=(y(i+1,1)-2.*y(i,1)+y(i-1,1))/dxi**2
phi(i,1)=-(xxi*xxxi+yxi*yxxi)/(xxi**2+yxi**2)
xxi=.5*(x(i+1,nj)-x(i-1,nj))/dxi
yxi=.5*(y(i+1,nj)-y(i-1,nj))/dxi
xxxi=(x(i+1,nj)-2.*x(i,nj)+x(i-1,nj))/dxi**2
yxxi=(y(i+1,nj)-2.*y(i,nj)+y(i-1,nj))/dxi**2
phi(i,nj)=-(xxi*xxxi+yxi*yxxi)/(xxi**2+yxi**2)
end do

do j=2,nj-1
xet=.5*(x(1,j+1)-x(1,j-1))/det
yet=.5*(y(1,j+1)-y(1,j-1))/det
xeet=(x(1,j+1)-2.*x(1,j)+x(1,j-1))/det**2
yeet=(y(1,j+1)-2.*y(1,j)+y(1,j-1))/det**2
pha(1,j)=-(xet*xeet+yet*yeet)/(xet**2+yet**2)
xet=.5*(x(ni,j+1)-x(ni,j-1))/det
yet=.5*(y(ni,j+1)-y(ni,j-1))/det
xeet=(x(ni,j+1)-2.*x(ni,j)+x(ni,j-1))/det**2
yeet=(y(ni,j+1)-2.*y(ni,j)+y(ni,j-1))/det**2
pha(ni,j)=-(xet*xeet+yet*yeet)/(xet**2+yet**2)
end do
1001 continue
do j=2,nj-1
do i=2,ni-1
residual=-10000.
xxi=.5*(x(i+1,j)-x(i-1,j))/dxi
xet=.5*(x(i,j+1)-x(i,j-1))/det
yxi=.5*(y(i+1,j)-y(i-1,j))/dxi
yet=.5*(y(i,j+1)-y(i,j-1))/det
aj = xxi*yet-xet*yxi

alpha = xet**2+yet**2
beta = xxi*xet+yxi*yet
beta=0.
gama = xxi**2+yxi**2

phi(i,j)=phi(i,1)+(phi(i,nj)-phi(i,1))*(j-1)&
/(nj-1)
pp=phi(i,j)*(xix**2+xiy**2)
pha(i,j)=pha(1,j)+(pha(ni,j)-pha(1,j))*(i-1)&
/(ni-1)
qq=pha(i,j)*(etx**2+ety**2)
! x equation
xtemp = (dxi*det)**2/(2.*(alpha*det*det+gama*dxi*dxi))&
* (alpha/(dxi*dxi)*(x(i+1,j)+x(i-1,j))&
+ gama/(det*det)*(x(i,j+1)+x(i,j-1))&
- beta/(2.*dxi*det)*(x(i+1,j+1)+x(i-1,j-1)&
- x(i-1,j+1)-x(i+1,j-1))+(aj*aj)*(xxi*pp+xet*qq))
ytemp = (dxi*det)**2/(2.*(alpha*det*det+gama*dxi*dxi))&
* (alpha/(dxi*dxi)*(y(i+1,j)+y(i-1,j))&
+ gama/(det*det)*(y(i,j+1)+y(i,j-1))&
- beta/(2.*dxi*det)*(y(i+1,j+1)+y(i-1,j-1)&
- y(i-1,j+1)-y(i+1,j-1))+(aj*aj)*(yxi*pp+yet*qq))

xtemp = omega*xtemp + (1.-omega)*x(i,j)
ytemp = omega*ytemp + (1.-omega)*y(i,j)
residual =max(residual,sqrt((x(i,j)-xtemp)**2)+sqrt((y(i,j)-ytemp)**2))
x(i,j)=xtemp
y(i,j)=ytemp
enddo
!c x(1,j)=x(2,j)
!c x(ni,j)=x(ni-1,j)
enddo
if(residual.gt.(1.d-6)) goto 1001
if(residual.le.(1.d-6)) continue

!c ************************ MAKE TECPLOTFILE ************************
OPEN(13,FILE='test.tec')
write(13,*) "VARIABLES = ", "X ", "Y "
write(13,*) "ZONE I=", ni, ",J=",nj, ", F=POINT"
DO j=1,nj
DO i=1,ni
write(13,1000) x(i,j),&
y(i,j)
ENDDO
ENDDO
CLOSE(UNIT=13)
1000 FORMAT(2E24.10)
100 FORMAT(I4.4)
END