IF YOU HAVE PROBLEM WITH DIFFERENTIAL MESH...........(2) -- 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:55:17 05/26/04 Wed
Author: Reza Khaksar Fard
Subject: IF YOU HAVE PROBLEM WITH DIFFERENTIAL MESH...........(2)

c.....
c..... Subroutine Set_Initial
c.....
subroutine Set_Initial
include 'param.inc'
include 'common.inc'
c.....
if(linearInterp) then
do I = 0,Imax
Xa = X(I,0)
Ya = Y(I,0)
Xb = X(I,Jmax)
Yb = Y(I,Jmax)
Xab = Xb-Xa
Yab = Yb-Ya
do J = 1,Jmax-1
Fraction = J*1./Jmax
X(I,J) = Xa+Xab*Fraction
Y(I,J) = Ya+Yab*Fraction
end do
end do
else
do I = 0,Imax
Xa = X(I,0)
Ya = Y(I,0)
Xb = X(I,Jmax)
Yb = Y(I,Jmax)
Xab = Xb-Xa
Yab = Yb-Ya
do J = 1,Jmax-1
Fraction = - 1./Sigma*log(1.-J*(1.-exp(-Sigma))/Jmax)
X(I,J) = Xa+Xab*Fraction
Y(I,J) = Ya+Yab*Fraction
end do
end do
end if
do J = 1,Jmax-1
X(Imax,J) = X(0,J)
Y(Imax,J) = Y(0,J)
end do
return
end

c.....
c..... Subroutine UnSt_LapSolver(Iter)
c.....
subroutine UnSt_LapSolver(Iter)
include 'param.inc'
include 'common.inc'
dimension UpperDiag(0:Imax),Diag(0:Imax),UnderDiag(0:Imax),
> Delta_X(0:Imax),Delta_Y(0:Imax),
> RHS_X(0:Imax),RHS_Y(0:Imax)
equivalence (RHS_X(0),RHS_Y(0))
equivalence (Delta_X(0),Delta_Y(0))
equivalence (RHS_X(0),Delta_X(0))
c.....
if(Iter.le.iter1) Omega = Omega1
if(Iter.gt.iter1.and.Iter.le.Iter2) Omega = Omega2
if(Iter.gt.iter2.and.Iter.le.Iter3) Omega = Omega3
if(Iter.gt.iter3) Omega = Omega4
call BoundCond
DelNorm = 0.0
c.....
c..... Calculation of X Coordinate
c.....
do J = 1,Jmax-1
do I = 0,Imax-1
dXdXi = (X(I+1,J)-X(I-1,J))*.5
dXdEta = (X(I,J+1)-X(I,J-1))*.5
dYdXi = (Y(I+1,J)-Y(I-1,J))*.5
dYdEta = (Y(I,J+1)-Y(I,J-1))*.5
c.....
Alpha = dXdEta**2+dYdEta**2
Beta = dXdXi*dXdEta+dYdXi*dYdEta
Gamma = dXdXi**2+dYdXi**2
c.....
UpperDiag(I) = Alpha
Diag(I) =-Alpha*2.-Gamma*2.
UnderDiag(I) = Alpha
RHS_X(I) =-Alpha*(X(I+1,J )-2.*X(I,J)+X(I-1,J ))
> -Gamma*(X(I ,J+1)-2.*X(I,J)+X(I ,J-1))
> +Beta *(X(I+1,J+1)-X(I+1,J-1)-
> X(I-1,J+1)+X(I-1,J-1))*.5
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
dXdXi = (X(I+1,J)-X(I-1,J))*.5
dXdEta = (X(I,J+1)-X(I,J-1))*.5
dYdXi = (Y(I+1,J)-Y(I-1,J))*.5
dYdEta = (Y(I,J+1)-Y(I,J-1))*.5
c.....
Alpha = dXdEta**2+dYdEta**2
Beta = dXdXi*dXdEta+dYdXi*dYdEta
Gamma = dXdXi**2+dYdXi**2
c.....
UpperDiag(I) = Alpha
Diag(I) =-Alpha*2.-Gamma*2.
UnderDiag(I) = Alpha
RHS_Y(I) =-Alpha*(Y(I+1,J )-2.*Y(I,J)+Y(I-1,J ))
> -Gamma*(Y(I ,J+1)-2.*Y(I,J)+Y(I ,J-1))
> +Beta *(Y(I+1,J+1)-Y(I+1,J-1)-
> Y(I-1,J+1)+Y(I-1,J-1))*.5
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 St_LapSolver(Iter)
c.....
subroutine St_LapSolver(Iter)
include 'param.inc'
include 'common.inc'
dimension UpperDiag(0:Imax),Diag(0:Imax),UnderDiag(0:Imax),
> Delta_X(0:Imax),Delta_Y(0:Imax),
> RHS_X(0:Imax),RHS_Y(0:Imax)
equivalence (RHS_X(0),RHS_Y(0))
equivalence (Delta_X(0),Delta_Y(0))
equivalence (RHS_X(0),Delta_x(0))
if(Iter.le.iter1) Omega = Omega1
if(Iter.gt.iter1.and.Iter.le.Iter2) Omega = Omega2
if(Iter.gt.iter2.and.Iter.le.Iter3) Omega = Omega3
if(Iter.gt.iter3) Omega = Omega4
call BoundCond
DelNorm = 0.0

CONTINUED.....

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