|
|
发表于 2003-4-23 01:42:42
|
显示全部楼层
谁有面元法的求解程序？
program main
c
parameter( npmx = 101 )
c
c---- All common blocks listed here.
c
common/bod / nlower, nupper, nodtot, np1, ss,
& x(npmx+1), y(npmx+1), xm(npmx), ym(npmx),
& costhe(npmx), sinthe(npmx)
common/cof / a(npmx,npmx+10), kutta
common/cpd / cp(npmx), xp
common/inf1 / aan(npmx,npmx), bbn(npmx,npmx)
common/sing / q(npmx), gamma
common/forces/ cl, cm
character*10 argi
c
pi = acos(-1.0)
c
c---- Retrieve parameters from the command line.
c
call getarg(1,argi)
read(argi,*) naca
call getarg(2,argi)
read(argi,*) nodtot
call getarg(3,argi)
read(argi,*) alpi
c
c---- Valid parameter check.
c
if (( naca .gt. 25999 ).or.( naca .lt. 1 )) naca = 12
if ( nodtot .lt. 20 ) nodtot = 20
if ( nodtot .gt. 100 ) nodtot = 100
if ( alpi .gt. 90.0 ) alpi = 90.0
if ( alpi .lt. -90.0 ) alpi = -90.0
c
nlower = nodtot / 2
nupper = nodtot - nlower
np1 = nodtot + 1
c
alpi = alpi*pi/180
c
c---- Generate the airfoil coordinates.
c
call airfoil(naca)
c
c---- Compute a steady solution at AOA = alpi.
c
call steady(alpi)
c
c---- Write the data.
c
call print_cp(naca,alpi)
c
stop
end
c------------------------------------------------------------------
subroutine airfoil(naca)
c
parameter( npmx = 101 )
c
common/bod / nlower, nupper, nodtot, np1, ss,
& x(npmx+1), y(npmx+1), xm(npmx), ym(npmx),
& costhe(npmx), sinthe(npmx)
c
c---- Compute the airfoil coordinates and panal angles.
c
pi = acos( -1.0 )
c
c---- Decompose the NACA number to determine airfoil coefficients.
c
ieps = naca / 1000
iptmax = naca / 100 - 10 * ieps
itau = naca - 1000 * ieps - 100 * iptmax
c
c---- Compute the coefficients.
c
epsmax = ieps * 0.01
ptmax = iptmax * 0.1
tau = itau * 0.01
c
c---- Error correction for bogus NACA numbers.
c
if (( naca .le. 9999 ) .and. ( epsmax .gt. 0 ) .and.
& ( ptmax .eq. 0 )) ptmax = .1
c
c---- If NACA 5 digit coding is used, make neccessary changes.
c
if ( ieps .ge. 10 ) then
if ( ieps .eq. 21 ) then
ptmax = 0.0580
ak1 = 361.4
elseif ( ieps .eq. 22 ) then
ptmax = 0.1260
ak1 = 51.64
elseif ( ieps .eq. 23 ) then
ptmax = 0.2025
ak1 = 15.957
elseif ( ieps .eq. 24 ) then
ptmax = 0.2900
ak1 = 6.643
elseif ( ieps .eq. 25 ) then
ptmax = 0.3910
ak1 = 3.230
endif
epsmax = ak1 * ptmax**3 / 6
endif
c
c---- initialize indexing for lower surface.
c
npoint = nlower
nstart = 0
c
c---- Loop over lower surface.
c
do n = 1, npoint
z = ( 1 + cos( pi * ( n-1 ) / npoint )) / 2
i = nstart + n
call naca45(naca,tau,epsmax,ptmax,z,thick,camber,beta)
x(i) = body_x(thick,beta,z,-1.0)
y(i) = body_y(thick,camber,beta,-1.0)
enddo
c
c---- Reinitialize indexing for upper surface
c
npoint = nupper
nstart = nlower
c
c---- Loop over upper surface.
c
do n = 1, npoint
z = ( 1 - cos( pi * ( n-1 ) / npoint )) / 2
i = nstart + n
call naca45(naca,tau,epsmax,ptmax,z,thick,camber,beta)
x(i) = body_x(thick,beta,z,1.0)
y(i) = body_y(thick,camber,beta,1.0)
enddo
c
c---- Load final point.
c
x(np1) = x(1)
y(np1) = y(1)
c
c---- Compute slopes of panals and arc length of airfoil skin.
c
ss = 0.0
do i = 1, nodtot
c
c---- Control points.
c
xm(i) = ( x(i+1) + x(i)) / 2
ym(i) = ( y(i+1) + y(i)) / 2
c
c---- Arc length.
c
dx = x(i+1) - x(i)
dy = y(i+1) - y(i)
dist = sqrt( dx**2 + dy**2 )
ss = ss + dist
c
c---- Slope.
c
sinthe(i) = dy / dist
costhe(i) = dx / dist
enddo
c
return
end
c------------------------------------------------------------------
function body_x(thick,beta,z,sign)
c
c---- Compute the x coordinate of an airfoil point.
c
body_x = z - sign * thick * sin( beta )
c
return
end
c------------------------------------------------------------------
function body_y(thick,camber,beta,sign)
c
c---- Compute the y coordinate of an airfoil point.
c
body_y = camber + sign * thick * cos( beta )
c
return
end
c------------------------------------------------------------------
subroutine cofish(alpha)
c
parameter( npmx = 101 )
c
common/bod / nlower, nupper, nodtot, np1, ss,
& x(npmx+1), y(npmx+1), xm(npmx), ym(npmx),
& costhe(npmx), sinthe(npmx)
common/cof / a(npmx,npmx+10), kutta
common/inf1 / aan(npmx,npmx), bbn(npmx,npmx)
c
pi = acos( -1.0 )
pi2inv = 0.5 / pi
kutta = nodtot + 1
cosalf = cos( alpha )
sinalf = sin( alpha )
c
c---- Initialize coefs.
c
do 10 j = 1, kutta
a(kutta,j) = 0.0
10 continue
c
c---- Set vn = 0 at midpoint of ith panel.
c
do 30 i = 1, nodtot
a(i,kutta) = 0.0
c
c---- Find contribution of jth panel.
c
do 20 j = 1, nodtot
if ( j .eq. i ) then
flog = 0.0
ftan = pi
else
dxj = xm(i) - x(j)
dxjp = xm(i) - x(j+1)
dyj = ym(i) - y(j)
dyjp = ym(i) - y(j+1)
flog = alog( ( dxjp**2 + dyjp**2 )
& / ( dxj**2 + dyj**2 )) / 2
ftan = atan2( dyjp * dxj - dxjp * dyj,
& dxjp * dxj + dyjp * dyj )
endif
ctimtj = costhe(i) * costhe(j) + sinthe(i) * sinthe(j)
stimtj = sinthe(i) * costhe(j) - costhe(i) * sinthe(j)
a(i,j) = pi2inv * ( ftan * ctimtj + flog * stimtj )
bbn(i,j) = pi2inv * ( flog * ctimtj - ftan * stimtj )
a(i,kutta) = a(i,kutta) + bbn(i,j)
if (( i .eq. 1 ) .or. ( i .eq. nodtot )) then
a(kutta,j) = a(kutta,j) - bbn(i,j)
a(kutta,kutta) = a(kutta,kutta) + a(i,j)
endif
c
c---- Load aan and bbn values now so that they don't have to be
c---- recomputed in infl.
c
aan(i,j) = a(i,j)
20 continue
c
c---- Fill in known sides.
c
a(i,kutta+1) = sinthe(i) * cosalf - costhe(i) * sinalf
30 continue
a(kutta,kutta+1) = -( costhe(1) + costhe( nodtot )) * cosalf
& -( sinthe(1) + sinthe( nodtot )) * sinalf
c
return
end
c------------------------------------------------------------------
subroutine formom(alpha)
c
parameter( npmx = 101 )
c
common/bod / nlower, nupper, nodtot, np1, ss,
& x(npmx+1), y(npmx+1), xm(npmx), ym(npmx),
& costhe(npmx), sinthe(npmx)
common/cpd / cp(npmx), xp
common/forces/ cl, cm
c
cosalf = cos( alpha )
sinalf = sin( alpha )
cfx = 0.0
cfy = 0.0
cm = 0.0
c
do i = 1, nodtot
c
c---- Moment coefficient is computed about the elastic axis.
c
xmid = xm(i) - xp
ymid = ym(i)
dx = x(i+1) - x(i)
dy = y(i+1) - y(i)
cfx = cfx + cp(i) * dy
cfy = cfy - cp(i) * dx
cm = cm + cp(i) * ( dx * xmid + dy * ymid )
enddo
c
c---- Compute lift.
c
cl = cfy * cosalf - cfx * sinalf
c
return
end
c------------------------------------------------------------------
subroutine gauss(nrhs,m,nitr)
c
parameter( npmx = 101 )
c
common/cof / a(npmx,npmx+10), neqns
c
c---- Performs Gaussian elimination of matrix a.
c
np = neqns + 1
ntot = neqns + nrhs
c
if (( m .le. 1 ) .and. ( nitr .eq. 0 )) then
c
c---- Do full matrix elimination sequence.
c
do 10 i = 2, neqns
im = i - 1
c
c---- Eliminate the (i-1)th unknown from i-th through
c---- neqns-th equations.
c
do 10 j = i, neqns
r = a(j,im) / a(im,im)
do 10 k = i, ntot
a(j,k) = a(j,k) - r * a(im,k)
10 continue
else
c
c---- Elimination on RHS only.
c
do 20 i = 2, neqns
im = i - 1
do 20 j = i, neqns
r = a(j,im) / a(im,im)
do 20 k = np, ntot
a(j,k) = a(j,k) - r * a(im,k)
20 continue
endif
c
c---- Back subtitution.
c
do 40 k = np, ntot
a(neqns,k) = a(neqns,k) / a(neqns,neqns)
do 40 l = 2, neqns
i = neqns + 1 - l
do 30 j = i+1, neqns
a(i,k) = a(i,k) - a(i,j) * a(j,k)
30 continue
a(i,k) = a(i,k) / a(i,i)
40 continue
c
return
end
c------------------------------------------------------------------
subroutine naca45(naca,tau,epsmax,ptmax,z,thick,camber,beta)
c
c---- Compute the thickness, camber, and angular location of an
c---- airfoil point.
c
c---- Thickness, corrected when z is very small.
c
if ( z .lt. 1.0e-10 ) then
thick = 0.0
else
thick = tau * 5 * ( 0.2969 * SQRT(Z)
& - Z * ( 0.1260
& + Z * ( 0.3537
& - Z * ( 0.2843
& - Z * 0.1015))))
endif
c
if ( epsmax .eq. 0.0 ) then
c
c---- For NACA 4-digit symmetrical arfoils.
c
camber = 0.0
beta = 0.0
else
if ( naca .gt. 9999 ) then
c
c---- For NACA 5 digit numbers.
c
c---- Ptmax = m and epsmax = (k_1*m^3)/6 from Abbott and Doenhoff.
c
if ( z .gt. ptmax ) then
camber = epsmax * ( 1.0 - z )
dcamdx = - epsmax
else
w = z / ptmax
camber = epsmax * ( w**3 - 3*w**2 +(3.-ptmax)*w)
dcamdx = epsmax/ptmax*(3*w**2 - 6*w + ( 3.0-ptmax))
endif
else
c
c---- For NACA 4 digit airfoils.
c
if ( z .gt. ptmax ) then
camber = epsmax / ( 1.0 - ptmax )**2
& * ( 1. + z - ptmax * 2 ) * ( 1. - z )
dcamdx = epsmax * 2 / ( 1.0 - ptmax )**2
& * ( ptmax - z )
else
camber = epsmax / ptmax**2 * ( ptmax*2 - z ) * z
dcamdx = epsmax * 2 / ptmax**2 * ( ptmax - z )
endif
endif
c
beta = atan( dcamdx )
endif
c
return
end
c------------------------------------------------------------------
subroutine print_cp(naca,alpi)
c
parameter( npmx = 101 )
c
common/bod / nlower, nupper, nodtot, np1, ss,
& x(npmx+1), y(npmx+1), xm(npmx), ym(npmx),
& costhe(npmx), sinthe(npmx)
common/cpd / cp(npmx), xp
c
c---- Write out data.
c
open(8,file='cp.fmt',form='formatted',status='unknown')
do n=1,nodtot
write(8,*) xm(n), ym(n), cp(n)
enddo
close(8)
c
return
end
c------------------------------------------------------------------
subroutine steady(alpha)
c
c---- alpha is angle of attack in radians, positive clockwise.
c
call cofish(alpha)
call gauss(1,0,0)
call veldis(alpha)
call formom(alpha)
c
return
end
c------------------------------------------------------------------
subroutine veldis(alpha)
c
parameter( npmx = 101 )
c
common/bod / nlower, nupper, nodtot, np1, ss,
& x(npmx+1), y(npmx+1), xm(npmx), ym(npmx),
& costhe(npmx), sinthe(npmx)
common/cof / a(npmx,npmx+10), kutta
common/cpd / cp(npmx), xp
common/inf1 / aan(npmx,npmx), bbn(npmx,npmx)
common/sing / q(npmx), gamma
c
cosalf = cos( alpha )
sinalf = sin( alpha )
c
c---- Retrieve solution from a matrix.
c
sum = 0.0
do i = 1, nodtot
q(i) = a(i,kutta+1)
ds = sqrt((x(i+1)-x(i))**2 + (y(i+1)-y(i))**2)
enddo
gamma = a(kutta,kutta+1)
c
c---- Find vt and cp at mid-point of i-th panel.
c
do i = 1, nodtot
vtfree = cosalf * costhe(i) + sinalf * sinthe(i)
vtang = vtfree
c
c---- Add contribution of j-th panel.
c
do j = 1, nodtot
vtang = vtang - bbn(i,j) * q(j) + gamma * aan(i,j)
enddo
cp(i) = 1.0 - vtang**2
enddo
c
return
end |
|
IP卡
狗仔卡
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
显身卡
楼主
