IF YOU HAVE PROBLEM WITH DIFFERENTIAL MESH...........(3) -- CFD

[ Show ]

Support VoyForums

[ Shrink ]

VoyForums Announcement: Programming and providing support for this service has been a labor of love since 1997. We are one of the few services online who values our users' privacy, and have never sold your information. We have even fought hard to defend your privacy in legal cases; however, we've done it with almost no financial support -- paying out of pocket to continue providing the service. Due to the issues imposed on us by advertisers, we also stopped hosting most ads on the forums many years ago. We hope you appreciate our efforts.

Show your support by donating any amount. (Note: We are still technically a for-profit company, so your contribution is not tax-deductible.) PayPal Acct: Feedback:

Donate to VoyForums (PayPal):

[ Login ] [ Contact Forum Admin ] [ Main index ] [ Post a new message ] [ Search | Check update time ]

[ Next Thread | Previous Thread | Next Message | Previous Message ]

Date Posted: 09:56:40 05/26/04 Wed
Author: Reza Khaksar Fard
Subject: IF YOU HAVE PROBLEM WITH DIFFERENTIAL MESH...........(3)

c.....
c..... Calculation of X Coordinate
c.....
do J = 1,Jmax-1
dEtaP = Eta(J+1)-Eta(J )
dEtaM = Eta(J )-Eta(J-1)
do I = 0,Imax-1
dXiP = Xi(I+1)-Xi(I )
dXiM = Xi(I )-Xi(I-1)
c.....
S1 = dXiM/dXiP/(dXip+dXiM)
S2 = (dXiP-dXiM)/dXiP/dXiM
S3 = dXiP/dXiM/(dXip+dXiM)
S4 = dEtaM/dEtaP/(dEtaP+dEtaM)
S5 = (dEtaP-dEtaM)/dEtaP/dEtaM
S6 = dEtaP/dEtaM/(dEtaP+dEtaM)
S7 = 2./dXiP/(dXip+dXiM)
S8 = 2./dXiP/dXiM
S9 = 2./dXiM/(dXiP+dXiM)
S10 = 2./dEtaP/(dEtaP+dEtaM)
S11 = 2./dEtaP/dEtaM
S12 = 2./dEtaM/(dEtaP+dEtaM)
c.....
dXdXiIJ = S1 *X(I+1,J )+S2 *X(I ,J )-S3 *X(I-1,J )
dYdXiIJ = S1 *Y(I+1,J )+S2 *Y(I ,J )-S3 *Y(I-1,J )
dXdXiIJP1 = S1 *X(I+1,J+1)+S2 *X(I ,J+1)-S3 *X(I-1,J+1)
dXdXiIJM1 = S1 *X(I+1,J-1)+S2 *X(I ,J-1)-S3 *X(I-1,J-1)
dXdEtaIJ = S4 *X(I ,J+1)+S5 *X(I ,J )-S6 *X(I ,J-1)
dYdEtaIJ = S4 *Y(I ,J+1)+S5 *Y(I ,J )-S6 *Y(I ,J-1)
d2XdXi2 = S7 *X(I+1,J )-S8 *X(I ,J )+S9 *X(I-1,J )
d2XdEta2 = S10*X(I ,J+1)-S11*X(I ,J )+S12*X(I ,J-1)
d2XdXiEta = S4 *dXdXiIJP1 +S5 *dXdXiIJ -S6 *dXdXiIJM1
c.....
Alpha = dXdEtaIJ**2+dYdEtaIJ**2
Beta = dXdXiIJ*dXdEtaIJ+dYdXiIJ*dYdEtaIJ
Gamma = dXdXiIJ**2 +dYdXiIJ **2
c.....
UpperDiag(I) = S7*Alpha
Diag(I) =-S8*Alpha-S11*Gamma-2.*S2*S5*Beta
UnderDiag(I) = S9*Alpha
RHS_X(I) =-Alpha*d2XdXi2+2.*Beta*d2XdXiEta-Gamma*d2XdEta2
end do
Il = 0
Iu = Imax-1
call Thomas(Il,Iu,UnderDiag,Diag,UpperDiag,RHS_X,Delta_X)
do I = 0,Imax-1
X(I,J) = X(I,J)+Omega*Delta_X(I)
c DelNorm = DelNorm+abs(Delta_X(I))
DelNorm = DelNorm+Delta_X(I)**2
end do
c.....
c..... Calculation of Y Coordinate
c.....
do I = 0,Imax-1
dXiP = Xi(I+1)-Xi(I )
dXiM = Xi(I )-Xi(I-1)
c.....
S1 = dXiM/dXiP/(dXip+dXiM)
S2 = (dXiP-dXiM)/dXiP/dXiM
S3 = dXiP/dXiM/(dXip+dXiM)
S4 = dEtaM/dEtaP/(dEtaP+dEtaM)
S5 = (dEtaP-dEtaM)/dEtaP/dEtaM
S6 = dEtaP/dEtaM/(dEtaP+dEtaM)
S7 = 2./dXiP/(dXip+dXiM)
S8 = 2./dXiP/dXiM
S9 = 2./dXiM/(dXiP+dXiM)
S10 = 2./dEtaP/(dEtaP+dEtaM)
S11 = 2./dEtaP/dEtaM
S12 = 2./dEtaM/(dEtaP+dEtaM)
c.....
dXdXiIJ = S1 *X(I+1,J )+S2 *X(I ,J )-S3 *X(I-1,J )
dYdXiIJ = S1 *Y(I+1,J )+S2 *Y(I ,J )-S3 *Y(I-1,J )
dYdXiIJP1 = S1 *Y(I+1,J+1)+S2 *Y(I ,J+1)-S3 *Y(I-1,J+1)
dYdXiIJM1 = S1 *Y(I+1,J-1)+S2 *Y(I ,J-1)-S3 *Y(I-1,J-1)
dXdEtaIJ = S4 *X(I ,J+1)+S5 *X(I ,J )-S6 *X(I ,J-1)
dYdEtaIJ = S4 *Y(I ,J+1)+S5 *Y(I ,J )-S6 *Y(I ,J-1)
d2YdXi2 = S7 *Y(I+1,J )-S8 *Y(I ,J )+S9 *Y(I-1,J )
d2YdEta2 = S10*Y(I ,J+1)-S11*Y(I ,J )+S12*Y(I ,J-1)
d2ydXiEta = S4 *dYdXiIJP1 +S5 *dYdXiIJ -S6 *dYdXiIJM1
c.....
Alpha = dXdEtaIJ**2+dYdEtaIJ**2
Beta = dXdXiIJ*dXdEtaIJ+dYdXiIJ*dYdEtaIJ
Gamma = dXdXiIJ**2 +dYdXiIJ **2
c.....
UpperDiag(I) = S7*Alpha
Diag(I) =-S8*Alpha-S11*Gamma-2.*S2*S5*Beta
UnderDiag(I) = S9*Alpha
RHS_Y(I) =-Alpha*d2YdXi2+2.*Beta*d2YdXiEta-Gamma*d2YdEta2
end do
Il = 0
Iu = Imax-1
call Thomas(Il,Iu,UnderDiag,Diag,UpperDiag,RHS_Y,Delta_Y)
do I = 0,Imax-1
Y(I,J) = Y(I,J)+Omega*Delta_Y(I)
c DelNorm = DelNorm+abs(Delta_Y(I))
DelNorm = DelNorm+Delta_Y(I)**2
end do
end do
c DelNorm = DelNorm/(Imax*Jmax)
DelNorm = sqrt(DelNorm/(Imax*Jmax))
return
end
c.....
c..... Subroutine Stretch_XiEta
c.....
subroutine Stretch_XiEta
include 'param.inc'
include 'common.inc'
do I = 0,Imax
Xi(I) = Real(I)/Imax
end do
do J = 0,Jmax
Eta(J) = (Real(J)/Jmax)**P_Eta
end do
Xi(-1) = -Xi(1)
return
end
c.....
c..... BoundCond
c.....
subroutine BoundCond
include 'param.inc'
include 'common.inc'
do J = 0,Jmax
X(-1 ,J) = X(Imax-1,J)
Y(-1 ,J) = Y(Imax-1,J)
X(Imax,J) = X(0,J)
Y(Imax,J) = Y(0,J)
end do
return
end
c.....
c..... Subroutine Thomas
c.....
subroutine Thomas(Il,Iu,Bb,Dd,Aa,Cc,Ee)
include 'param.inc'
dimension Aa(0:imax),Bb(0:imax),Cc(0:imax),Dd(0:imax),Ee(0:imax)
Lp = Il+1
do I = Lp,Iu
R = Bb(I)/Dd(I-1)
Dd(I) = Dd(I)-R*Aa(I-1)
Cc(I) = Cc(I)-R*Cc(I-1)
end do
Ee(Iu) = Cc(Iu)/Dd(Iu)
do I = Lp,Iu
J = Iu-I+Il
Ee(J) = (Cc(J)-Aa(J)*Ee(J+1))/Dd(J)
end do
return
end

CONTINUED.....

[ Next Thread | Previous Thread | Next Message | Previous Message ]