AER E 361 Computational Techniques for Aerospace Design
Aerodynamic Characteristics of a NACA 23012 Airfoil ------------------------------------------------------------------------------------------------------------------------------------------------------Computation and Analysis
Group Members:
Jonathan DeVries
Ray Anaya
Alex Bakken
Abstract Aerodynamic analysis is a subject rooted with a rich history of genius contributions th
made by some of the greatest minds of the 20 century. For decades, these minds poured over wind tunnels, computers, and the classic pen and paper with the goal of discovering analysis techniques than those that were previously offered. We explore three methods commonly used in aerodynamic analysis: thin airfoil theory, the discrete vortex panel panel method, and a full full airfoil vortex panel method. This analysis is performed on an NACA 23012, in inviscid, incompressible, steady flow. The bulk of this report covers the use of a vortex panel method which takes the airfoil’s
entire contour into account; account; which, as expected, gave stunning stunning accuracy. For instance: the NACA 23012 has a design lift coefficient of 0.3 (the first digit of an NACA 5-series multiplied by 0.15 gives the coefficient), while our numerical method provides a l ift coefficient of 0.298 at zero angle of attack. While there are no real implications of our results, the comparison of our computational results to theoretical and experimental values provides a good lesson on when and how the various techniques should be used. Thus, this project was an excellent learning learning tool for numerical methods in aerospace design.
2
Contents List of Symbols .............................................................................................................................................. 4 1. Introduction .............................................................................................................................................. 5 2. Theory ....................................................................................................................................................... 6 2.1 Kutta condition conditio n................................... .................. .................................. ................................... ................................... ................................... .................................... .......................... ........ 6 2.2 Thin Airfoil Theory ................................. ................ .................................. ................................... ................................... .................................. ................................... ........................ ...... 6 2.3 Thin Airfoil Theory for Cambered Airfoils .................... ............................. ................... ................... .................. .................. .................. .................. ................ ....... 8 2.4 Discrete Vortex Panel Method .......................................................................................................... 10 2.5 Vortex Panel Method for Full Airfoil ................................................................................................. 10 3. Analysis ................................................................................................................................................... 14 4. Results ..................................................................................................................................................... 17 5. Conclusion ............................................................................................................................................... 19 References .................................................................................................................................................. 20 Appendix ..................................................................................................................................................... 21
3
List of Symbols
Influence coefficient of bound vortex j on control point i (Dimensionless) Fourier coefficient Influence coefficient of wake vortex ( Dimensionless) Chord Length (dimensionless) Section lift coefficient Moment coefficient about leading edge Moment coefficient about quarter-chord point Pressure coefficient Control Point index Bound Vortex and vortex panel index Lift per unit span Vortex strength (dimensionless) Moment about leading edge Number of Panels used to simulate the airfoil Distance vector from bound vortex j to control point i Trailing Edge Velocity tangential to panel Flow velocity at infinite distance from the airfoil (dimensionless) Coordinate of collocation point Location of center of pressure Airfoil angle of attack Constant vortex Circulation Panel j length Angle between influence radius vectors Variable used for integral transformations Instantaneous deflection of control surface Dummy variable Density (dimensionless)
4
1. Introduction In this report, we discuss three different methods of airfoil analysis: Thin Airfoil Theory, Discretized Vortex Panel method, and the normal 2-dimensional panel method. Thin airfoil theory is a great “fly “fl y by night” technique for estimating trends in aerodynamic characteristics,
however, it make several assumptions which leave the result much to be desired when used for any practical application. The discretized vortex panel method is a more accurate method method for analyzing an airfoil; however, it still lacks the accuracy accuracy necessary for rigorous design. While the first two methods may seem completely undesirable for any sort of practical use, they work great for very thin airfoils. The advantage of of the normal vortex 2-dimensional 2-dimensional method, however, is its ability to analyze an airfoil with noticeable thickness while taking into account the entire contour of the airfoil. The goal of this analysis was to use computational methods to determine the aerodynamic characteristics characteristics of an NACA 23012 airfoil and to compare these values with theoretical and experimental results.
5
2. Theory 2.1 Kutta condition Simply stated, the Kutta condition describes nature’s tendency to adopt a certain value of circulation such that the fluid flow leaves the trailing edge of the airfoil smoothly. One result of this is that, as shown in Figure 2.1, the velocities of f low leaving the upper and lower surfaces of the trailing edge must be both equal to zero for steady flow. However, adding a cusp to the airfoil provides such a small angle between the two surfaces at the leading edge that, due to the fact that both flows are orthogonal to each other, the velocities need not be equal to zero; yet they still must be equal to each other to satisfy this condition. Because the velocities are equal at the trailing edge of any airfoil, the vorticity at the trailing edge will always be zero, or
2.2 Thin Airfoil Theory
(2.1)
Thin airfoil theory is a basic tool for calculating fundamental aerodynamic characteristics characteristics of thin, symmetrical airfoils. By placing a vortex sheet along the camber camber line, the airfoil can be simulated. Following the derivation outlined in Anderson(4.7), one can arrive at the the fundamental equation of thin airfoil theory
∫
(2.2)
which states that that the camber line is a streamline of of the flow. Further, in order to satisfy the the Kutta condition Eq.(2.2) needs to be satisfied for Anderson(324). The solution to (2.2) is
. This can be verified by referring to
(2.3)
⁄
The fundamental fundamental equation of thin airfoil theory can theory can be used to develop basic aerodynamic characteristics of a thin, symmetrical airfoil. Because the airfoil is thin, thin, the term The total circulation around the airfoil is
(2.4)
∫
Using the following transformations,
6
.
∫
(2.5) (2.6)
And substituting the result into Eq.(2.3),
(2.7)
[1]
From the Kutta-Jukowski theorem ,
(2.8)
and from fundamental aerodynamics knowledge,
we arrive at
(2.9)
(2.10)
which is an important result of thin airfoil theory. This relation shows that the theoretical theoretical lift coefficient of an airfoil varies linearly l inearly with angle of attack, which can be verified experimentally
α
for any symmetric airfoil. Note that lift-coefficient computations computations must use units for angle of attack
in radians.
Thin airfoil theory may also be expanded to resolve theoretical values for the moment coefficient of a given airfoil. The equation for the theoretical leading edge edge moment is found to be
∫
(2.11)
Using Eqs(2.5,2.6) and performing the integration, (2.11) becomes (2.12)
From fundamental aerodynamics knowledge
and from Eq.(2.10)
7
(2.13)
(2.14)
the moment coefficient about the leading edge is found to be
(2.15)
More importantly, the moment coefficient about the quarter-chord point can be found from Eq.(2.15). Since
and substituting (2.15) into (2.16)
⁄ ⁄
(2.16)
(2.17)
This result is important because it demonstrates that the theoretical center of pressure is at the quarter-chord point for symmetric airfoils. A plot of the theoretical lift coefficient can be found in Figure(4.1).
2.3 Thin Airfoil Theory for Cambered Airfoils Thin airfoil theory can be further expanded for cambered cambered airfoils. For a cambered airfoil, dz/dx is dz/dx is no longer zero, and a finite result for the derivative needs to be obtained. Solving Eq.(2.2), we find that
∑ ⁄ ∫ ∑ ∫ ∫ ⁄ ∑ ∫
Our goal is to use Eq.(2.18) to arrive at values for (2.2)
and
(2.18)
. Substituting (2.18) into
(2.19)
and evaluating the integrals
and reducing (2.19) using the results obtained herein, a value for
(2.20)
is found: (2.21)
A Fourier analysis, covered in Hildebrand(217), gives the final two values sought:
and 8
(2.22)
∫
(2.23)
To obtain aerodynamic coefficients, the the circulation must be determined. Substituting (2.18) into (2.4) and using the transformations (2.5) and (2.6)
[ ∫ ∑ ∫ ] ∫ * + ⁄ * +
(2.24)
Solving (2.24) gives
(2.25)
and the lift per unit span is
(2.26)
from which the lift coefficient becomes
This result shows that the lift coefficient slope for any shape airfoil is
(2.27)
. Referring to
Figure(4.9) in Anderson, the geometry shows explicitly
(2.28)
where the zero-lift angle of attack is found to be
(2.29)
By extension, the equation for the moment coefficient about the leading edge is (2.30)
Substituting (2.30) into (2.16) gives the moment coefficient about the quarter chord: (2.31)
The location of the center of pressure is known from fundamental aerodynamics knowledge: (2.32)
Substituting (2.30) into (2.32) gives the center of pressure as a function of the lift coefficient:
9
(2.33)
While the center of pressure locations and leading edge moments depend on the l ift coefficient, and thus angle of attack, attack, the quarter chord moment does not. Thus, this moment stays around a constant value for each airfoil. This extension of thin airfoil theory is important because it can be used for any type of thin airfoil, even symmetric airfoils. Evaluating these these equations for a symmetric symmetric airfoil would result in everything reducing to their symmetric thin airfoil theory counterparts.
2.4 Discrete Vortex Panel Method The vortex panel method has greater capabilities than thin airfoil theory, where thickness is no longer a restriction. It is analogous to the source panel method, with which it
shares many characteristics. The method is the same, with the exception of introducing a vortex. By adding the vortex, the airfoil now has a circulation
the quantity of which is the
total imbalance of pressure distribution, which effectively causes lift. To start the Vortex panel method, a mean camber line can be used to approximate the aerodynamic characteristics, just like in thin airfoil theory. The mean camber line however is discretized into N number of panels, each a straight line. Source/sinks are placed, yielding best results with as few panels as possible at the ¼ chord point of each panel. The strengths of the singularities affect each panel at its control point, or the ¾ chord of each panel. However, assumptions assumptions for this simplified version must be made. No Kutta condition is involved since no circulation is introduced, which gives a system of N of N equations with N unknowns. Then, the use of the geometry geometry and influence coefficients gains access to the singularity strengths. Since this is a simplified version of the normal panel method, we can think of the mean camber line as essentially the top half of an airfoil.
2.5 Vortex Panel Method for Full Airfoil The discrete vortex panel method can be expanded to use not only the mean camber line, but the shape of the airfoil itself. This is the pinnacle pinnacle of this analysis: an accurate way way of numerically determining aerodynamic characteristics of an airfoil. To begin, the airfoil shape is divided into
number of straight straight panels. This can be
accomplished in a number of ways, although this group chose to use a method discussed in the analysis section below. These panels are then used to determine five necessary characteristics characteristics of each panel. First, a fictitious control point must be chosen for each panel, which is determined by
Second, the length of each panel needs to be determined:
10
.
(2.34)
√
.
(2.35)
Third, the angles between each panel and the horizontal need to be computed: (2.36)
Fourth, the influence radii between each panel:
(( ) ( ) (()()()()())(()()()()())
(2.37)
Finally, the angles between each influence radius:
when
when
and
(2.38a)
(2.38b)
In order to satisfy the fact that there is i s no velocity within the airfoil body, a no-penetration condition must be applied. Applying this condition yields:
for
∑ ( ) ( ) ∑ ( ) ( ) ∑
(2.39)
, and where
(2.40)
and
(2.41)
and
(2.42)
The Kutta condition must must also be satisfied. Applying this yields:
where 11
(2.43)
∑ ( ) ( ) ∑ ∑ ( ) ( )
(2.44)
and
(2.45)
and
The
terms form what is known as the aerodynamic aerodynamic influence coefficient matrix, coefficient matrix, or
This matrix is used to form a system of
where
(2.46)
is the vortex strength,
(2.47)
equations, in the form of
is the vorticity, and
is the column vector
formed by (2.42) and (2.46).
Using the solutions to this system of equations, pressure coefficients and lift coefficients can be obtained. The pressure pressure coefficient is derived from Bernoulli’s Equation as
(2.48)
where
( )( ) ( ) ( ) (2.49) [5]
The lift coefficient is then found using
12
∑
(2.50)
This more accurate lift coefficient will in turn give more accurate values of moment coefficients and locations of center of pressure, using the equations developed in the thin airfoil theory discussion.
13
3. Analysis The first step to the vortex panel airfoil analysis was to obtain coordinates for the [3]
airfoil’s profile from the University of Illinois at Urbana Champaign’s airfoil database . For an
analysis using twelve panels, data points were removed from the coordinate profile until only every fifth point remained, which became the nodes of each panel. panel. This method can be repeated for any number of panels, up to the maximum amount of data points in the coordinate file. A plot showing the airfoil contour, panels, and mean camber line is shown in Figure(3.1).
Figure(3.1) Airfoil contour showing discretized panels After the previous step was completed, a Fortran Fortran code was written for the actual analysis. analysis. It starts by taking the nodes as input and passing them through Eq.(2.34) to determine the collocation points on each panel. These node points are also used in Eq(2.36) to determine the angles between each panel. Note that this analysis was was completed using a unit free-stream velocity. The collocation points are then used to calculate influence radii between each panel and their corresponding angles. These quantities are used in Eqs.(2.37-2.46) Eqs.(2.37-2.46) to determine each element of the aerodynamic coefficient matrix, Eq(2.47). Eq(2.47). Figure(3.2) shows the resulting influence influence matrix.
14
Figure(3.2). AIC Matrix. The vortex strengths
and vorticity
equations, where
are found as solutions to (2.47) to form a system of
is the unknown, and
. This system is solved through through
Gaussian elimination with with partial pivoting. These values are presented in Figure(3.3). Figure(3.3). x coll. 0.966635 0.84197 0.625935 0.37624 0.16163 0.035975 0.031015 0.15667 0.37624 0.625935 0.84197 0.966635
x node 1.00003 0.93324 0.7507 0.50117 0.25131 0.07195 0 0.06203 0.25131 0.50117 0.7507 0.93324
theta(i) -3.010591 -3.0472641 -3.0783181 -3.1377904 3.0472593 2.7965539 0.714649 0.1164926 -0.04799 -0.1071899 -0.1380448 -0.1742877
m(i) 1.1027468 0.9915361 0.9598768 0.9893461 1.0603056 1.2512424 0.2681776 -0.8212711 -1.076278 -1.1661928 -1.2429978 -1.4052664
gamma 0.1396484 0.1396484 0.1396484 0.1396484 0.1396484 0.1396484 0.1396484 0.1396484 0.1396484 0.1396484 0.1396484 0.1396484
delta(s) 0.0673673 0.1833551 0.2500304 0.2498618 0.180161 0.0764561 0.0821238 0.1905716 0.250148 0.2509704 0.1842932 0.0678175
cp 0.100001 0.010001 -0.03007 -0.0066 0.09863 0.549998 -1.10001 -0.9495 -0.55002 -0.30002 0.000107 0.081664
vt -0.94868 -0.99499 -1.01492 -1.0033 -0.94941 -0.67082 1.449142 1.396246 1.244999 1.140184 0.999946 0.958298
Figure(3.3) Tabulated output data for angle of attack of 4 degrees The resulting values for vortex strength and vorticity are used in Eqs.(2.50,2.59) to find a computational lift coefficient and and tangential velocities of each panel. The tangential velocities are further used to find center of pressure distribution around the the airfoil, using Eq(2.48). The pressure distribution is shown in Figure(4.2), for an angle of attack of four degrees. The vortex panel analysis was compared to theoretical, experimental, and previous computational coefficients for this airfoil. The lift coefficient, computed here as 0.568 at four degrees angle of attack, does not differ [4]
much from other methods. methods. The theoretical value gives 0.438, the experimental value gives 0.5, and the previous computational computational method provides 0.53. More important, however, is the comparison of lift coefficients at zero angle of attack. attack. As shown in section two of this report, the theoretical lift coefficient from thin airfoil theory gives an angle of of attack of zero. The experimental value is nearly 0.1, and the previous computed computed value is 0.15. For the full vortex panel method, the lift coefficient is computed as 0.298, which is negligibly smaller than the
15
design lift coefficient of 0.3. Thus, the vortex panel panel method for full airfoils is found found to be an extremely accurate approximation, approximation, even when using only twelve panels for analysis. The moment coefficient about the quarter chord point did not need values obtained from the vortex panel method. method. Thus, there are only two values values of concern: those obtained from both types of thin airfoil theory. The value from symmetrical thin thin airfoil theory always returns returns a value of zero. However, a cambered airfoil produces produces additional lift that is independent of the angle of attack; it acts through a point point on the airfoil which depends on the airfoil’s shape. This
value is calculated as -0.01283, and is not dependent on angle of attack. The pressure coefficient distribution, as computed by the full vortex panel method, was compared to the theoretical theoretical value as being nearly spot-on with theoretical values. An analysis using more than twelve panels would increase the accuracy of these data, however.
16
4. Results This section provides figures and tables of data that were achieved using methods described in this report.
Figure(4.1). Comparison of lift coefficients obtained using various methods.
Figure(4.2). Pressure distribution of NACA 23012 at 4 degrees angle of attack.
17
Figure(4.3). Center of pressure location as function of angle of attack.
18
5. Conclusion There are many ways to to analyze the aerodynamic characteristics characteristics of an an airfoil. The method which is chosen, though, depends on how much accuracy the analyst wants to achieve. Some applications may only need the amount of accuracy achieved by thin airfoil theory alone, however, thick and complex airfoils such as the 23012 b enefit greatly from using the full vortex panel method. Many times in engineering analysis level of accuracy during during analysis has a strong correlation with cost cost of design. Therefore, it is beneficial if the engineer engineer is able to differentiate between these numerical methods, and apply each one to situations where they are most appropriate. That said, this analysis provided a great perspective on what is produced with increasing accuracy. The vortex panel method produced produced values that coincide with what is to be expected for this airfoil. The lift coefficient at zero angle angle of attack is almost almost exactly at the design design lift coefficient value.
19
References [1] Anderson, John D. Fundamentals of Aerodynamics. Aerodynamics . 5th ed. New York: McGraw-Hill, 2007. Print. [2] Hildebrand, F.B.: Advanced F.B.: Advanced Calculus for Applications , 2d ed., Prentice-Hall Inc., Englewood Cliffs, N.J., 1976. [3] "NACA23012.dat." "NACA23012.dat." UIUC Airfoil Coordinates Database. Database . UIUC Applied Aerodynamics Group, 2011. Web. 26 March 2012. . se.html>. [4] "NACA 23012 12%." Airfoil 12%." Airfoil Investigation I nvestigation Database. Database . UIUC Airfoil Coordinates Database, 25 Feb. 2012. Web. 31 Mar. 2012. . . [5] “Panel Method.” Methods and Tools for Aerospace Aerospace Design. Ambar K. Mitra, Iowa Iowa State
university Dept. of Aerospace Aerospace Engineering, 2010. 2010. Web. 31 Mar 2012. . web/index.html>. [6] Rajagopalan, Ganesh. AerE 361 Lecture. Pearson 1115, Ames, IA. Spring 2012. Lecture.
20
Appendix vortexpanel panels coeff partialpiv
PROGRAM USE USE USE
USE write USE
aerodynamics
IMPLICIT NONE REAL , DIMENSION (20)
:: x, y, xc4, yc4, theta, b, m, dels :: r, eta, a, ln, vt, cp INTEGER, DIMENSION(20) :: INDX REAL , ALLOCATABLE :: ans(:) INTEGER :: n, i REAL :: alpha, vinf, delsum, cl WRITE(*,*) "what is the number of panels?" READ(*,*) n ALLOCATE (ANS(20)) WRITE(*,*) "what is alpha in radians?" read(*,*) alpha CALL panelstuff(x,y,xc4,yc4,r,eta,theta,n,ln,dels,delsum) CALL AIC(ln,eta,theta,n,a,b,alpha,vinf) CALL gaussianelim(n+1,a,b,ans) !cl = ((4*ans(13))/vinf)*delsum CALL aero(vt,cp,n,ans,theta,ln,eta,vinf,alpha,delsum) CALL output(n,x,xc4,theta,dels,ans,cp,vt) END PROGRAM vortexpanel REAL , DIMENSION (20,20)
21
MODULE
panels
CONTAINS
SUBROUTINE
mcl()
IMPLICIT NONE
:: pi, deltatheta :: xmcl, ymcl, x, y INTEGER :: n, i, error, j REAL REAL
pi = 3.141593 deltatheta = 30 * pi / 180 OPEN (unit
= 12, file = "naca23012.dat" ) OPEN (unit = 14, file = "nacatest2.dat" ) n = 10000
i = 1, 61 x, y IF (x <= 1 .AND. x > 0.2025) THEN ymcl = 0.02208*(1-x) 0.02208*(1-x) ELSE IF (x <= 0.2025 .AND. x >= 0) THEN ymcl = 2.6595*(x**3-0.6075*x** 2.6595*(x**3-0.6075*x**2+0.1147*x) 2+0.1147*x)
DO
READ(12,*)
END IF IF
(x == 0)
EXIT
WRITE(14,*)
x, ymcl deltatheta = deltatheta + (30*pi/180) END DO END SUBROUTINE SUBROUTINE
mcl
panelstuff(x,y,xc4,yc4,r,eta,theta,n,ln,dels,delsum)
IMPLICIT NONE REAL ,DIMENSION(20), INTENT (OUT)
:: x, y, xc4, yc4,theta, dels :: r, eta, ln
REAL ,DIMENSION(20,20), INTENT (OUT)
:: j, error, i, jp1, ip1 :: n REAL :: pi REAL , INTENT (OUT) :: delsum INTEGER
INTEGER, INTENT(IN)
pi = 4.*ATAN 4.*ATAN(1.) (1.) OPEN (unit
= 12, file = 'panels2.dat' 'panels2.dat') ) OPEN (unit = 13, file = 'collocations.dat' ) OPEN (unit = 14, file = 'r.dat' 'r.dat') ) OPEN (unit = 15, file = 'eta.dat' 'eta.dat') ) 'theta.dat') ) OPEN (unit = 16, file = 'theta.dat' OPEN (UNIT = 17, FILE = 'lnn.dat' 'lnn.dat') ) OPEN (UNIT = 18, FILE = 'lnmat.dat' 'lnmat.dat') ) 'xyc4.dat') ) OPEN (UNIT = 19, FILE = 'xyc4.dat' DO j = 1, n+1 READ(12,*,iostat=error) x(j),y(j) IF (j == 13) EXIT END DO
22
delsum = 0 dels(1) = SQRT SQRT((x(1)-x(2))**2+(y(1)-y(2))**2) ((x(1)-x(2))**2+(y(1)-y(2))**2) DO i = 1, n ip1 = i+1 IF (i==n) THEN ip1 = 1 END IF
(i >= 2) THEN dels(i) = SQRT((x(i)-x(ip1))**2+(y(i)-y(ip1))**2) SQRT((x(i)-x(ip1))**2+(y(i)-y(ip1))**2)
IF
END IF
delsum = dels(i) + delsum END DO
j = 1,n jp1 = j+1 ! xc4(j) = 3*(x(jp1)-x(j))/4 + x(j) ! yc4(j) = 3*(y(jp1)-y(j))/4 + y(j) xc4(j) = (x(j)+x(jp1))/2 yc4(j) = (y(j)+y(jp1))/2 WRITE(19,*) xc4(j), yc4(j) DO
END DO
xc4(12) = xc4(1) DO j = 1,n WRITE(13,*) xc4(j), yc4(j) END DO
i = 1,n ip1 = i+1 IF (i == n) ip1 = 1
DO
THEN
END IF WRITE(17,*)
"i=",i "i=",i theta(i) = ATAN2 ATAN2(y(ip1)-y(i),x(ip1)-x(i)) (y(ip1)-y(i),x(ip1)-x(i)) WRITE(14,*) "i= ",i ",i WRITE(15,*) "i= ",i ",i WRITE(16,*) "i= ",i ",i WRITE(16,*) theta(i) DO j = 1,n jp1 = j+1 r(i,j) = SQRT SQRT((xc4(i)-x(j))**2+(yc4(i)-y(j))**2) ((xc4(i)-x(j))**2+(yc4(i)-y(j))**2) r(i,jp1) = SQRT SQRT((xc4(i)-x(jp1))**2+(yc4(i)-y(jp1))**2) ((xc4(i)-x(jp1))**2+(yc4(i)-y(jp1))**2) eta(i,j) = ATAN2 ATAN2((yc4(i)-y(jp1))*(xc4(i)-x(j))-(xc4(i)-x(jp1))*(yc4(i)-y(j)),(yc4(i)((yc4(i)-y(jp1))*(xc4(i)-x(j))-(xc4(i)-x(jp1))*(yc4(i)-y(j)),(yc4(i)y(jp1))*(yc4(i)-y(j))+(xc4(i)-x(jp1))*(xc4(i)-x(j))) IF (i == j) THEN eta(i,j) = pi END IF
ln(i,j) = LOG LOG(r(i,jp1)/r(i,j)) (r(i,jp1)/r(i,j)) WRITE(17,*) ln(i,j) WRITE(15,*) eta(i,j) WRITE(14,*) r(i,j), r(i,jp1) END DO END DO DO
i = 1, n
(18,'(12(f10.7,1X))' ) WRITE(18,'(12(f10.7,1X))'
(ln(i,j),j=1,n)
END DO
panelstuff panels
END SUBROUTINE END MODULE
23
MODULE
coeff
CONTAINS SUBROUTINE
AIC(ln,eta,theta,n,a,b,alpha,vinf)
IMPLICIT NONE REAL , DIMENSION (20,20), INTENT(IN)
:: eta, ln REAL , DIMENSION (20), INTENT (IN) :: theta REAL , DIMENSION (20,20), INTENT(OUT) :: a REAL , DIMENSION (20,20) :: sumavortextan, lnterm, r REAL , DIMENSION (20), INTENT (OUT) :: b REAL , DIMENSION (20,20) :: avortexnor, avortextan, asourcenor, asourcetan INTEGER :: i, j, k, np1, jp1, ip1 INTEGER, INTENT(IN) :: n REAL :: pi REAL , INTENT (OUT) :: Vinf REAL , INTENT (IN) :: alpha 'atest.dat') ) OPEN (UNIT = 12, FILE = 'atest.dat' OPEN (UNIT = 13, FILE = 'ln.dat' 'ln.dat') ) OPEN (UNIT = 14, FILE = 'atest2.dat' 'atest2.dat') ) OPEN (UNIT = 15, FILE = 'bvector.dat' 'bvector.dat') ) OPEN (UNIT = 18, FILE = 'eta2.dat' 'eta2.dat') ) pi = 4.*ATAN 4.*ATAN(1.) (1.) Vinf = 1 np1 = n+1 i=1, n ip1 = i+1 b(i) = Vinf*sin Vinf* sin(theta(i)-alpha) (theta(i)-alpha) DO j=1, n jp1 = j+1 a(i,j) = (sin ( sin(theta(i)-theta(j))*ln(i,j))/(2*pi)+(eta(i,j)* (theta(i)-theta(j))*ln(i,j))/(2*pi)+(eta(i,j)* cos cos(theta(i)(theta(i)theta(j)))/(2*pi) DO
END DO END DO
b(13) = -vinf*cos -vinf* cos(theta(1)-alpha)-vinf* (theta(1)-alpha)-vinf* cos cos(theta(12)-alpha) (theta(12)-alpha) DO i = 1, n (13,'(12(f10.7,1X))' ) (ln(i,j),j=1,n) WRITE(13,'(12(f10.7,1X))' END DO
i = 1, n j = 1, n WRITE(12,*) a(i,j)
DO
DO
END DO END DO
a(0,np1) = 0 DO i = 1, n DO j = 1, n avortexnor(i,np1) avortexnor(i,np1) = ( cos cos(theta(i)-theta(j))*ln(i,j))/(2*pi)-(eta(i,j)* (theta(i)-theta(j))*ln(i,j))/(2*pi)-(eta(i,j)* sin sin(theta(i)(theta(i)theta(j))/(2*pi)) a(i,np1) = a(i,np1) + avortexnor(i,np1) END DO END DO
a(np1,0) = 0 DO j = 1, n DO k = 1, n
24
asourcetan(np1,j) = (eta(k,j)*sin (eta(k,j)* sin(theta(k)-theta(j)))/(2*pi)-( (theta(k)-theta(j)))/(2*pi)-(cos cos(theta(k)(theta(k)theta(j))*ln(k,j))/(2*pi) a(np1,j) = a(np1,j)+asourcetan(np1,j) END DO END DO
avortextan(i,np1) avortextan(i,np1) = 0 DO i = 1, n DO j = 1, n sumavortextan(i,np1) = (eta(i,j)* cos cos(theta(i)-theta(j)))/(2*pi)+( (theta(i)-theta(j)))/(2*pi)+(sin sin(theta(i)(theta(i)theta(j))*ln(i,j))/(2*pi) avortextan(i,np1) = avortextan(i,np1) + sumavortextan(i,np1) END DO END DO
a(np1,np1) = 0 DO i = 1, n a(np1,np1) = a(np1,np1) + avortextan(i,np1) END DO DO
i = 1, np1
WRITE(14,'(13(f10.7,1X))' (14,'(13(f10.7,1X))' ) WRITE(15,*)
(a(i,j),j=1,np1)
b(i)
END DO
AIC coeff
END SUBROUTINE END MODULE
25
MODULE
partialpiv
CONTAINS SUBROUTINE
gaussianelim(n, aa, bb, xx)
IMPLICIT NONE INTEGER, INTENT(IN)
:: n :: aa(:,:), bb(:) REAL , ALLOCATABLE , INTENT (OUT) :: xx(:) INTEGER :: i, j, k, m, prow REAL(KIND=8) REAL (KIND=8) :: btemp, factor, pivot, det, detmin, sum REAL(KIND=8), REAL (KIND=8), ALLOCATABLE :: a(:,:), b(:), x(:), atemp(:) REAL , INTENT (IN)
ALLOCATE (a(n,n), b(n), x(n), xx(n), atemp(n)) a = DBLE DBLE(aa) (aa) b = DBLE(bb) DBLE(bb)
k=1, n-1 !~~Partial Pivoting~~! pivot = DABS DABS(a(k,k)) (a(k,k)) prow = k DO m=k+1, n IF (DABS DABS(a(m,k)) (a(m,k)) > pivot) pivot = DABS DABS(a(m,k)) (a(m,k)) prow = m
DO
THEN
END IF END DO
atemp(:) = a(k,:) a(k,:) = a(prow,:) a(prow,:) = atemp(:) btemp = b(k) b(k) = b(prow) b(prow) = btemp !~~Gaussian Elimination~~! DO i=k+1, n factor = a(i,k)/a(k,k) DO j=k, n a(i,j) = a(i,j) - factor*a(k,j) factor*a(k,j) END DO
b(i) = b(i) - factor*b(k) END DO END DO
!~~Determine if solution exists~~! det = 1.d0 detmin = a(1,1) DO i=1, n det = det*a(i,i) IF (DABS DABS(a(i,i)) (a(i,i)) < DABS DABS(detmin)) (detmin)) detmin = a(i,i) END DO
detmin = DABS DABS(detmin)*1.d-6 (detmin)*1.d-6 DABS(det) (det) < detmin) THEN IF (DABS WRITE (*,*) "The determinant of the reduced matrix is 0" WRITE (*,*) "The system of equations has no real solution." x(:) = 0. END IF
!~~Backward Substitution~~! DABS(det) (det) >= detmin) THEN IF (DABS 26
x(n) = b(n)/a(n,n) DO i=n-1, 1, -1 sum = 0.d0 DO j=i+1, n, 1 sum = sum + (a(i,j)*x(j)) END DO
x(i) = (b(i) - sum)/a(i,i) END DO END IF
xx = SNGL SNGL(x) (x) RETURN END SUBROUTINE END MODULE
MODULE
gaussianelim
partialpiv
aerodynamics 27
CONTAINS SUBROUTINE
aero(vt,cp,n,ans,theta,ln,eta,vinf,alpha,delsum)
IMPLICIT NONE REAL , DIMENSION (20), INTENT (IN)
:: ans, theta!, vt, cp, vsum1, vsum2,
vsumfactor,vsumfactor2 REAL , DIMENSION (20,20), INTENT(IN) :: eta, ln !REAL, DIMENSION(20,20), INTENT(OUT) :: vt, cp REAL , DIMENSION (20) :: vsum1, vsum2, vsumfactor, vsumfactor2 REAL , DIMENSION (20), INTENT (OUT) :: vt, cp !REAL, DIMENSION(20,20) :: vsum1, vsum2, vsumfactor, vsumfactor2 REAL , INTENT (IN) :: vinf, alpha, delsum REAL :: pi, cl, cmle, xcp, cl_theor, cmle_theor, cmc4_theor, xcp_theor, alphaL0_theor, A0, A1, A2 INTEGER, INTENT(IN) :: n INTEGER :: dummy, i, j pi = 4.*ATAN 4.*ATAN(1.) (1.) WRITE(*,*) delsum, ans(13) cl = ((2*ans(13))/vinf)*delsum dummy = 1 DO
i = 1, n+1 ans(i)
WRITE(*,*) END DO
vt(i) = 0 DO i = 1, n DO j = 1, n vsumfactor(i) vsumfactor(i) = vinf*cos vinf*cos(theta(i)-alpha)+(ans(i)/(2*pi))*(eta(i,j)* (theta(i)-alpha)+(ans(i)/(2*pi))*(eta(i,j)* sin sin(theta(i)(theta(i)theta(j))-cos theta(j))- cos(theta(i)-theta(j))*ln(i,j))+(ans(13)/(2*pi))*( (theta(i)-theta(j))*ln(i,j))+(ans(13)/(2*pi))*( sin sin(theta(i)(theta(i)theta(j))*ln(i,j)+eta(i,j)* cos cos(theta(i)-theta(j))) (theta(i)-theta(j))) vt(i) = vt(i)+vsumfactor(i) vt(i)+vsumfactor(i) END DO END DO DO
i = 1, n
WRITE(*,*)
vt(i)
END DO
i = 1, n cp(i) = 1.-(vt(i)/vinf)**2 1.-(vt(i)/vinf)**2 ! WRITE(*,*) cp(i,dummy) DO
END DO
cmle = -cl/4. xcp = -cmle/cl
!Coefficients calculated by hand A0 = .041155768 A1 = .0954840633 A2 = .0791502466 cl_theor = pi*(2.*A0+A1) cmle_theor = -.25*(cl_theor + pi*(A1-A2)) cmc4_theor = .25*pi*(A2-A1) xcp_theor = .25*(1. + pi/cl_theor*(A1-A2)) 28
alphaL0_theor = alpha - (A0+.5*A2) file= 'theoretical.dat' ) OPEN (unit=20, file='theoretical.dat' WRITE (20,*) cl_theor WRITE (20,*) cmle_theor WRITE (20,*) cmc4_theor WRITE (20,*) xcp_theor WRITE (20,*) alphaL0_theor alphaL0_theor
WRITE(*,*)
"Cl=",cl, "Cl=" ,cl,"CmLE=" "CmLE=",cmle, ,cmle,"xcp=" "xcp=",xcp ,xcp aero aerodynamics
END SUBROUTINE END MODULE
29
MODULE write CONTAINS SUBROUTINE
output(n,x,xc4,theta,dels,ans,cp,vt)
IMPLICIT NONE REAL , DIMENSION (20), INTENT (IN)
:: x, xc4, theta, ans, dels :: cp, vt
REAL , DIMENSION (20,20), INTENT(IN) INTEGER, INTENT(IN)
:: n :: i, j, dummy OPEN (UNIT = 12, FILE = 'tabulation.dat' ) dummy = 1 DO i = 1, n (12,'(8(f13.7,1X))' ) xc4(i), x(i), theta(i), ans(i), ans(13), dels(i), cp(i,dummy), WRITE(12,'(8(f13.7,1X))' vt(i,dummy) INTEGER
END DO END SUBROUTINE
output
END MODULE write
30
