[lehoanganh07]

Dear CFD user,
I am using fortran code to make the grid for my model: 2D nozzle.
The method is structure meshing: algibra method with stretching grid near the wall. The mesh look good, but when I ran my simulation code, after 10000 calculation step, the result become NaN. And I think the error cause of the grid.

This is my grid code (simple)

C************************************************* *******************
C FILE NAME is tecy.f
C************************************************* *******************
PARAMETER(NI1=81,NI2=201,NI3=201,NI=483,NJ=81)
IMPLICIT REAL*8 (A-H,K,L,O-Z)
C
DIMENSION X(NI,NJ),
& Y(NI,NJ),
& A(NI,NJ),
& X1(NI1,NJ),
& Y1(NI1,NJ),
& X2(NI2,NJ),
& Y2(NI2,NJ),
& X3(NI3,NJ),
& Y3(NI3,NJ),
& HD1(NI1,NJ),
& HD2(NI2,NJ),
& HD3(NI3,NJ),
& YL1(NI1,NJ),
& YL2(NI2,NJ),
& YL3(NI2,NJ)
C
OPEN(2,FILE="Mesh_Nozzle_Rectang")
C OPEN(3,FILE="X")
C OPEN(4,FILE="Y")
C
H = 32./1000
H1 = 4./1000
L1 = 150./1000
L3 = 16./1000
L4 = 300./1000
H2 =(H-H1)/2.
ET2=.5
ET1=H2/H+0.2*H1/H
ET3=(H-H2-0.2*H1)/H
ANPHA =3.
ANPHAX=3.
ANPHAY=3.
DXI1 =1./(NI1-1)
DXI2 =1./(NI2-1)
DXI3 =1./(NI3-1)
DET =1./(NJ-1)
C ********************BOUNDARY WALL BLOCK1****************
DO I=1,NI1
Y1(I,1)=0.
Y1(I,NJ)=H
E=Y1(I,NJ)-Y1(I,1)
YL1(I,1)=Y1(I,1)
HD1(I,1)=E
LXI=DXI1*(I-1)
K =EXP(ANPHAX1)-1.
K1=EXP(ANPHAX1)-EXP(ANPHAX1*(1.-LXI))
X1(I,1)=L1*K1/K
X1(I,NJ)=X1(I,1)
ENDDO
C -----------ALB GRID----------------------------
DO J=1,NJ
DO I=1,NI1
LET=DET*(J-1)
IF(LET.LE.ET2)THEN
Y1(I,J)= HD1(I,1)*ET2*(EXP(ANPHAY*LET/ET2)-1.)/(EXP(ANPHAY)-1.)
& +YL1(I,1)
ELSEIF(LET.GE.ET2)THEN
Y1(I,J)= HD1(I,1)*(1.-(1.-ET2)*(EXP(ANPHAY*(1.-LET)/(1.-ET2))-1.)
& /(EXP(ANPHAY)-1.)) +YL1(I,1)
ENDIF
LXI=DXI1*(I-1)
K =EXP(ANPHAX1)-1.
K1=EXP(ANPHAX1)-EXP(ANPHAX1*(1.-LXI))
X1(I,J)=L1*K1/K
ENDDO
ENDDO
C
C ********************BOUNDARY WALL BLOCK2****************
DO I=1,NI2
Y2(I,1)=H2
Y2(I,NJ)=H2+H1
E=Y2(I,NJ)-Y2(I,1)
YL2(I,1)=Y2(I,1)
HD2(I,1)=E
LXI=DXI2*(I-1)
X2(I,1)=LXI*L3
X2(I,NJ)=X2(I,1)
ENDDO

C -----------ALB GRID----------------------------
DO J=1,NJ
DO I=1,NI2
LET=DET*(J-1)
IF(LET.LE.ET2)THEN
Y2(I,J)= HD2(I,1)*ET2*(EXP(ANPHAY*LET/ET2)-1.)/(EXP(ANPHAY)-1.)
& +YL2(I,1)
ELSEIF(LET.GE.ET2)THEN
Y2(I,J)= HD2(I,1)*(1.-(1.-ET2)*(EXP(ANPHAY*(1.-LET)/(1.-ET2))-1.)
& /(EXP(ANPHAY)-1.)) +YL2(I,1)
ENDIF
LXI=DXI2*(I-1)
X2(I,J)=LXI*L3
ENDDO
ENDDO
C
C ********************BOUNDARY WALL BLOCK3****************
DO I=1,NI3
Y3(I,1)=0.
Y3(I,NJ)=H
E=Y3(I,NJ)-Y3(I,1)
YL3(I,1)=Y3(I,1)
HD3(I,1)=E
LXI=DXI3*(I-1)
K =EXP(ANPHAX2)-1.
K1=EXP(ANPHAX2*LXI)-1.
X3(I,1)=L4*K1/K
X3(I,NJ)=X3(I,1)
ENDDO
C -----------ALB GRID----------------------------
DO J=1,NJ
DO I=1,NI3
LET=DET*(J-1)
IF(LET.LE.ET2)THEN
Y3(I,J)= HD3(I,1)*ET2*(EXP(ANPHAY*LET/ET2)-1.)/(EXP(ANPHAY)-1.)
& +YL3(I,1)
ELSEIF(LET.GE.ET2)THEN
Y3(I,J)= HD3(I,1)*(1.-(1.-ET2)*(EXP(ANPHAY*(1.-LET)/(1.-ET2))-1.)
& /(EXP(ANPHAY)-1.)) +YL3(I,1)
ENDIF
LXI=DXI3*(I-1)
K =EXP(ANPHAX2)-1.
K1=EXP(ANPHAX2*LXI)-1.
X3(I,J)=L4*K1/K
ENDDO
ENDDO
C
C ****************MODEL DOMAIN*************************
DO I=1,NI1
DO J=1,NJ
X(I,J)=X1(I,J)
Y(I,J)=Y1(I,J)
ENDDO
ENDDO
DO J=1,NJ
A(NI1,J)=X(NI1,J)
ENDDO
DO I=1,NI2
DO J=1,NJ
X(I+NI1,J)=X2(I,J)+A(NI1,J)
Y(I+NI1,J)=Y2(I,J)
ENDDO
ENDDO
DO J=1,NJ
A(NI1+NI2,J)=X(NI2+NI1,J)
ENDDO
DO I=1,NI3
DO J=1,NJ
X(I+NI1+NI2,J)=X3(I,J)+A(NI1+NI2,J)
Y(I+NI1+NI2,J)=Y3(I,J)
ENDDO
ENDDO
C
DO J=1,NJ
DO I=1,NI
C WRITE(3,10) X(I,1)
C10 FORMAT(F18.16)
C WRITE(4,11) Y(I,1)
C11 FORMAT(F18.16)
WRITE(2,10) X(I,J),Y(I,J)
10 FORMAT(F18.16,2X,F18.16)
ENDDO
ENDDO
C
END
[/SIZE][/I][/SIZE]

My model is like this, I devided Nozzle into 3 block: block1, block2 (Nozzle hole) and block 3 and make grid for each block. After that I combine into one model as in the code.
http://img69.imageshack.us/img69/4205/wkgt.png

Could anyone help me to solve this problem?, is that my code correct?