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

YOU CAN FIND BOTH THE UNSTRETCHED AND STRETCHED LAPLACE IN THE FOLLOWING CODE.

/=========================FORTRAN======================/

program GridProject1
include 'param.inc'
include 'common.inc'
include 'fig.inc'
logical Converged
open(unit=1,file='airf.dat')
open(unit=2,file='init.dat')
open(unit=3,file='final.dat')
open(unit=4,file='Res.dat')
open(unit=100,file='Input.dat')
c.....
call ReadInput
call Set_AlgebBound
Printcode=1
call PrintResults(Printcode)
call Set_Initial
Printcode=2
call PrintResults(Printcode)
Converged = .false.
Delmax = 0.
Res = 0.
if(StretchMesh) then
call Stretch_XiEta
do Iter = 1,60
call St_LapSolver(Iter)
if(Iter.eq.1) Delmax = DelNorm
Res = alog10(DelNorm/Delmax)
write(5,*) Iter,Res
if(Iter-Iter/10*10.eq.0) then
write(4,*) Iter,Res
end if
call CheckConverg(Converged)
if(iter.lt.50) converged=.false.
if(Converged) goto 100
end do
else
do Iter = 1,100000
call UnSt_LapSolver(Iter)
if(Iter.eq.1) Delmax = DelNorm
Res = alog10(DelNorm/Delmax)
write(5,*) Iter,Res
if(Iter-Iter/10*10.eq.0) then
write(4,*) Iter,Res
end if
call CheckConverg(Converged)
if(Converged) goto 100
end do
end if
100 continue
Printcode=3
call PrintResults(Printcode)
close(1)
close(2)
close(3)
close(4)
end
c.....
c..... Subroutine ReadInput
c.....
Subroutine ReadInput
include 'param.inc'
include 'common.inc'
namelist
>/Input/Naca,Error,Sigma,P_Eta,OuterRadius,Omega1,Omega2,
> Omega3,Omega4,Iter1,Iter2,Iter3,
> LinearInterp,StretchMesh,
> IbackColor,ImeshColor,IairfoilColor
read(100,Input)
close(100)
IMFactor = Naca/1000
MFactor = IMFactor/100.
IP = Naca/100
IPFactor = (IP-IP/10*10)
PFactor = IPFactor/10.
ITauFactor = (Naca-Naca/100*100)
TauFactor = ITauFactor/100.
return
end
c.....
c..... Subroutine Set-AlgebBound
c.....
subroutine Set_AlgebBound
include 'param.inc'
include 'common.inc'
Pi2 = 2.*acos(-1.)
ImaxD2 = Imax/2
c..... Internal Boundary
do I = 0,Imax
Theta = I*Pi2/Imax
X_IntBound = .5*(1.+cos(Theta))
if(X_IntBound.lt.PFactor) then
YCamber = MFactor/PFactor**2
> *(2.0*PFactor*X_IntBound-X_IntBound**2)
else
YCamber = MFactor/(1.-PFactor)**2
> *((1.-2.*PFactor)+2.*PFactor*X_IntBound-X_IntBound**2)
end if
YThickD2 = 5.*TauFactor*(.2969*sqrt(X_IntBound)-0.126
> * X_IntBound-.3516*X_IntBound**2+.2843*X_IntBound**3
> -.1015*X_IntBound**4)
if(I.le.ImaxD2) then
Y_IntBound = YCamber + YThickD2
else
Y_IntBound = YCamber - YThickD2
end if
X(I,0) = X_IntBound
Y(I,0) = Y_IntBound
c..... External Boundary
X(I,Jmax) = OuterRadius*cos(Theta)+.5
Y(I,Jmax) = OuterRadius*sin(Theta)
end do
X(Imax,0) = 1.
X(0 ,0) = 1.
Y(Imax,0) = 0.
Y(0 ,0) = 0.
X(Imax,Jmax) = OuterRadius+.5
X(0 ,Jmax) = OuterRadius+.5
Y(Imax,Jmax) = 0.
Y(0 ,Jmax) = 0.
return
end

CONTINUED...

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