A Low Order Model for Vertical Axis Wind Turbines
Isaac M. Asher?, Mark Drela?and Jaime Peraire?
Massachusetts Institute of Technology, Cambridge, MA 02139, U.S.A.
A new computational model for initial sizing and performance prediction of vertical axis wind turbines
is presented. The model uses a 2D hybrid dynamic vortex and blade element momentum approach. Each
airfoil is modeled as a single vortex of time varying strength with an analytical model for the influence of the
shed vorticity. The vortex strengths are calculated by imposing a flow tangency condition at the three-quarter
chord location on each airfoil, modified in the case of stall. The total blade forces and the momentum-based
streamtube deceleration are then obtained using pre-computed cd and cm 2D blade profile characteristics.
Model fidelity is improved over previous models because flow curvature, dynamic vortices, blade interactions,
static stall, and streamtube changes are all taken into account. Fast convergence is obtained for a large range
of solidity and tip speed ratio, which allows optimization of various parameters, including blade pitch angle
variation.
Nomenclature
Letters
A frontal area, taken as 2R (unit span).
A0, A1 pitch angle modulation coefficients.
b = 2/?, wake vorticity calibration constant.
c blade airfoil chord.
C(?) stall model leakage velocity.
cd coefficient of drag.
cl coefficient of lift computed using ?.
?cl = 0.05, stall onset parameter.
clmin,max stall thresholds.
cm coefficient of moment.
CP = NbP? /
1
2
?V 3?A, average power coefficient.
CT = NbT? /
1
2
?V 2?A, average thrust coefficient.
Note that CT > 0 corresponds to windmilling.
~D = 1
2
?~V 2c/4ccd, drag force (along ~Vc/4).
h = 1
2
, lift curve slope parameter.
Ks = 40, stall ?leakage? velocity factor.
Kcd drag coefficient factor.
~L = ??(~Vc/4 ??z?) + ?c??(~Vc/4 ??z?)/Vc/4,
lift force.
? stream-tube parameter.
~M = 1
2
?(~V 2c/4 ? ??)c2cm,
moment, positive counter clockwise.
~Mtot = ~M + ~rc/4 ? (~L + ~D),
moment due to ~M , drag, and lift.
Nb number of blades.
~n unit vector normal to airfoil camber line at ~rcp.
P = ~Mtot ? ~?, instantaneous power.
R radius of turbine (origin to leading edge of blade).
~rc/4 quarter chord point.
~rcp control point, the three-quarter chord.
~rle vector from x, y origin to leading edge of foil.
Re chord Reynolds number.
T thrust (component of ~L and ~D in x?-direction).
TSR = ?R/V?, tip speed ratio.
~V? external freestream velocity far upstream.
~V ??(~r) external freestream velocity evaluated at ~r
using stream tube deceleration model.
~Vrot = ?~?? ~r velocity due to rotation ~?,
~Vcp relative velocity at the control point
~Vc/4 relative velocity at the quarter chord
used for calculating aerodynamic forces,
includes vortex-induced velocities.
~V? = ??~n/?c, velocity induced at ~rcp
by the vortex at ~rc/4.
~Vw = [?b??/~Vcp ? s?]~n, wake-induced velocity.
Symbols
?0 angle between ??? and chord line s?,
positive clockwise about the leading edge.
? strength of vortex located at ~rc/4,
positive clockwise.
? rotation angle of blade,
measured counter-clockwise from +y-axis.
? density of air.
?Research Assistant, Dept. of Aeronautics and Astronautics, 77 Massachusetts Avenue Room 37-442, Cambridge, MA 02139,
isaaca@alum.mit.edu, AIAA Student Member.
?Terry J. Kohler Professor of Fluid Dynamics, Dept. of Aeronautics and Astronautics, AIAA Fellow.
?Professor, Dept. of Aeronautics and Astronautics, AIAA Associate Fellow.
1 of 9
American Institute of Aeronautics and Astronautics
? = c/R, solidity.
? pitch angle modulation phase offset.
~? = ??z? angular velocity of blade rotation,
positive counterclockwise.
?0 = ??0, pitch rate of blade.
Coordinate systems
x, y non-rotating, inertial, absolute frame.
r, ? frame centered at the leading edge of the airfoil
with r? away from the origin and ?? counterclockwise.
s, t frame centered at the leading edge of the airfoil
with s? along the chord line and t? perpendicular.
z perpendicular to the plane (spanwise direction).
I. Introduction
Darrieus-type vertical axis wind turbines (VAWTs) have advantages in omni-directionality and structural simplicity
over traditional horizontal axis wind turbines (HAWTs), but they are more complex aerodynamically. Current VAWT
designs are in general not as efficient as HAWTs, and one contributing factor is a deficiency in design tools. Most
simple, low order models cannot capture enough physics to reliably predict performance. More accurate models tend
to be too slow and prevent fast design iterations. The tool developed in this project is computationally efficient and
incorporates many of the salient features of the accurate tools, thus giving the designer the ability to search the design
space effectively.
The model uses the Blade Element/Momentum method introduced by Templin1. The entire turbine is assumed to
be a single actuator disc. The blade interactions captured by the more intricate multiple and double-multiple streamtube
models of Strickland2 and Paraschivoiu3 are here captured with the airfoil vortices. Variations in the streamtube
deceleration function form (constant, linear, and inverse tangent) had little effect on the resulting CP , so it is assumed
that the single actuator disc with linear deceleration is sufficiently accurate.
The airfoils are modeled by a vortex at the quarter chord of varying strength as Strickland4 did. The shed vorticity
(due to varying vortex strength) is set to a 2D analytic approximation. The model does not track the shed vorticity,
since its influence is generally small far from the airfoil. Strickland extended his vortex method to include static and
dynamic stall and ?flow curvature?. Static stall and flow curvature are both included in this model.
The model is a hybrid between the Blade Element/Momentum and Vortex methods, and is accurate and fast enough
to be wrapped in an optimization loop and give good results. One shortcoming is that solidity must not be allowed
to grow too large during optimization, because there is little penalty in the model for very high solidities (blade
interactions are approximate).
This model does not include 3D effects that are common to previous models, although wind shear and tip vortices
could easily be included. At the moment, the model is useful for assessing the performance of straight-bladed turbines
with high blade aspect ratios, moderate to low solidities, and moderate to high tip speed ratios.
II. Assumptions
We take the simple case of a 2D airfoil rotating in the plane at constant angular velocity ? (see Fig. 1 for geometry
setup). The freestream is equal to ~V? far from the airfoil. Near the airfoil, the freestream decelerates linearly with the
net deceleration determined by the average thrust produced by the airfoil?. The air flow around the airfoil is assumed to
be incompressible and to have constant Reynolds and Mach numbers. The precise velocity and pressure distributions
are approximated by replacing the airfoils with a discrete set of vortices. One main vortex is located at the quarter
chords of the airfoils, and a secondary wake vortex is located at the trailing edges?. The net induced velocity from the
freestream, kinematics (rotation and pitch angle changes), and the main vortices of other airfoils (though not the wake
vortices of the other airfoils) is used to calculate the forces on a given vortex. Losses (drag) and moments are captured
with locally 2D airfoil profile drag characteristics?. The chosen reference units are R = 1, ~V? = 1x?, and ? = 1. We
enforce a flow tangency condition and the thrust coefficient equation and solve for ?(?) and CT .
?The distance over which the flow decelerates is set a-priori.
?The wake vorticity is proportional to the shed vorticity ??.
?Taken at the nominal Reynolds number. Given the vortex strengths, we can compute cl and use cd(cl) data. For cl beyond stall, a linear
extrapolation is used. The same is done for the moment, except that cm is assumed zero beyond stall.
2 of 9
American Institute of Aeronautics and Astronautics
x?y?
?
~?
??
R
r?
?0
s?
t?
~rcp~rc/4
?
~n
~Vcp
s?
t?
?
~Vc/4
~L
~D
~M
Forces and moments on
? due to net induced
velocity at ~rc/4
Figure 1. Turbine geometry and model setup
III. Model Algorithm
A. Streamtube Deceleration
For the streamtube model only, the turbine is assumed to be a single actuator disc, and ~V ?? is calculated using the
stream tube deceleration implied by CT . A linear deceleration over a distance 2? is assumed?,
~V ?? =
?????
????
~V?, x < ??,
~V?
[
1 +
?
1?CT?1
2
(
1 + x?
)]
, ?? ? x ? ?,
~V?
(?
1 ? CT
)
, x > ?.
An example ~V ?? distribution for ? = 6 and CT = 0.8 is show in Fig. 2.
B. Blade-relative Velocities
The various contributing velocities at ~rcp are calculated. The freestream velocity is calculated with the streamtube
deceleration model. The kinematic velocity, or the blade-observed velocity due to its movement (varying ? and ?0) is
~Vrot = ?
[
~? ? ~rle + c3
4
??s
]
.
?? is chosen a-priori
3 of 9
American Institute of Aeronautics and Astronautics
The total freestream velocity, which corresponds to the terms not multiplying ? in Eq. (2), is then
~Vcp = ~V
?
? + ~Vrot = ~V
?
?(~rcp) ? ~?? ~rle ? c
3
4
??s.
The velocity at ~rcp due to the main vortex (which becomes the term multiplying ? in Eq. (2)) is
~V? = ?? ~n
?c
.
?10 ?5 0 5 10
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
1.1
x
F
re
es
tr
ea
m
ve
lo
ci
ty
Streamtube Deceleration Model, ? = 6, CT = 0.8
Figure 2. Example of streamtube velocity distribution
Note that the direction is taken as ~n, which is not
strictly true if the airfoil has camber. However, it
is important numerically to preserve the magnitude
of ~V?.
The influence of the main vortices of the other
airfoils are also calculated. The effect of the main
vortex of foil i on the control point for foil j is
~V?i,j = ?
~rcp,j ? ~rc/4,i
2?|~rcp,j ? ~rc/4,i|2
.
The wake velocity is based on a calibrated ap-
proximation to the velocity induced by the trailing
wake vortex sheet5,
~Vw =
?b
~Vcp(~rcp) ? s?
??~n.
This introduces a ?? term, which turns Eq. (2) into
a rate equation for ?.
C. Flow Tangency and Stall Model
The model assumes that the airfoil can be represented by a vortex at the quarter chord of strength ?(t). The standard
practice is to compute a ?(?) to enforce flow tangency at a control point on the airfoil. That is, a ? is found such that
~V (~rcp) ? ~n = 0. no-stall (1)
Here ~V includes contributions from the rotation (~Vrot), the external flow (~V ??), the vortex (~V?), and the wake (~Vw).
In this model, we add a simple stall model by setting the right hand side of Eq. (1) to some non-zero velocity ?leakage?
at high cl. The equation becomes
~V (~rcp) ? ~n = Vcp
4?h
Ks?cl log
1 + exp [(cl ? clmax)/?cl]
1 + exp [(clmin ? cl)/?cl] , stall (2)
where we take h = 1/2, Ks = 40, and ?cl = 0.05. clmax and clmin are found manually. The leakage velocity is
plotted in Fig. 3 versus cl.
The control point is chosen to be the 3
4
-chord point? to be consistent with thin airfoil theory. The normal vector is
perpendicular to the airfoil camber line at the control point (at x = 3
4
).
Once a ? is found to satisfy Eq. (2), outputs such as lift, drag, moment, power out, and power coefficient may be
computed. Lift is computed using the vortex strength, and drag and moment are computed using the lift and 2D airfoil
section polars. Power coefficient is computed at each ? and integrated numerically to obtain CP .
?Note that ~rcp here always lies along the chord line, rather than on the camber line or aligning with the local freestream. This is done to simplify
the computation.
4 of 9
American Institute of Aeronautics and Astronautics
D. Circulation Evolution Equation
Now Eq. (2) is of the form (
?~n
?
? + ~Vcp +
?b
~Vcp ? s?
??~n
)
? ~n = C(?), (3)
where C is the right hand side that constitutes the stall model (which depends on cl and therefore ?). Note that CT
is an input (needed to calculate ~Vcp) and an output (based on forces generated), so both ? and CT must be computed
iteratively. This is done by Newton iterations to drive the residual of the above equation to zero. In order to obtain an
expression for CT , we need to calculate the forces on the blades.
E. Force and Moment Calculation
Given a guess of ? and CT , the velocity at any point in the resulting flow field can be found. Forces and moments are
evaluated using the ?infinity? velocity at the vortex location
~Vc/4 = ~V (~rc/4) = ~V
?
?(~rc/4) ? ~?? ~rle ? c
1
4
??s? b??~n
~Vcp(~rcp) ? s?
+ ~V?i .
The ?? term accounts for the blade?s own near-wake shed vorticity, and ~V?i term accounts for the bound vortices of the
other Nb ? 1 blades. The near-wake shed vorticity of the other blades is ignored. ~L is calculated directly from ?
~L = ??(~Vc/4 ??z?) + ?c??(~Vc/4 ??z?)/Vc/4,
and a corresponding cl is found
cl =
2(? + c??/Vc/4)
cVc/4
.
Note that during iterations, the ?? used here is the derivative of the guess ?. 2D airfoil section data are used to find cd
and cm corresponding to the cl (see below). The drag and moment are calculated
~D = 1
2
?~V 2c/4ccd,
~M = 1
2
?(~V 2c/4 ? ??)c2cm.
Note that the moment is positive nose down, so positive moment reinforces the motion.
F. Drag and Moment Coefficient Model
For each ?, the associated cl is mapped to a cd and cm using blade profile characteristics6, augmented by the stall
model. For clmin < cl < clmax, the polar file data are interpolated using a 10th order polynomial, thereby obtaining cd
and cm. The stall behavior implies a change in the cd curve for cl outside this range. The cd(cl) imported from the
polar files is augmented by a cdstall :
cd =
??????
?????
cd(cl), clmin < cl < clmax
cd(clmax) +
dcd
dcl
????
clmax
(cl ? clmax), cl > clmax
cd(clmin) +
dcd
dcl
????
clmin
(cl ? clmin), cl < clmin
. (4)
This defines how the drag polar should be extended through stall, which is a simple linear function cd(cl), as shown
in Fig. 4.
Note that, during iterations, the cl used here is calculated from the input ?. Since the polynomials are quite invalid
beyond the input data range, any cl outside the range of the polar file will give cm = 0.
5 of 9
American Institute of Aeronautics and Astronautics
?2 ?1 0 1 2
?15
?10
?5
0
5
10
15
Lift coefficient
L
ea
ka
g
e
V
el
o
ci
ty
Stall model: leakage velocity
? Clmin
? Clmax
Figure 3. Control point leakage velocity for stall model.
0 0.02 0.04 0.06 0.08 0.1
?2
?1.5
?1
?0.5
0
0.5
1
1.5
2
Drag coefficient
L
if
t
co
effi
ci
en
t
Stall model: drag polar extension
? Clmin
? Clmax
Figure 4. 2D profile drag polar
G. Thrust and Power Coefficients
The thrust, net moment about the origin, and power are calculated for one foil from the force, with positive power
corresponding to the power extracted from the air,
T = (~L + ~D) ? x?, P = ~Mtot ? ~? =
(
~M + ~rc/4 ? (~L + ~D)
)
? ~?.
The thrust and power are averaged over the cycle (using the Midpoint rule for the integration), and the thrust and power
coefficients for the entire turbine are then calculated
CT =
TNb
1
2
?V 2?A
, CP =
PNb
1
2
?V 3?A
. (5)
The thrust is normalized by the freestream dynamic pressure and turbine area, and the power extracted is normal-
ized by the total power available in the external flow through the turbine region of area A.
IV. Solution Method
All quantities are expressed in terms of ? and are implicitly assumed periodic. We discretize ? into m points?
?i = 2?i/m, i = 0..m ? 1. The unknowns that we solve for are then ?i = ?(?i) and the averaged quantity CT .
Assuming that ? is sufficiently smooth, a spectrally accurate finite differentiation matrix7 is used to calculate ?? from
?,
??j ?
m
2
?1?
i=?m
2
+1
di?j+i, di =
{
1
2
(?1)i+1 cot (?im ) , m 6= 0
0, m = 0
.
Initial guesses to ? and CT are supplied by a simpler model (no stall or ??, which only requires iterations on CT ), or
assumed zero. Newton iterations to drive point-wise residuals of Eqs. (2) and (5) to zero. The residuals are
R? =
(
?~n?? + ~Vcp + ?b~Vcp?s? ??~n
)
? ~n? C(?),
RCT = CT ? TNb1
2
?V 2
?
A
.
?m is taken as an integer multiple of the number of blades, we use m = 16 for the results below.
6 of 9
American Institute of Aeronautics and Astronautics
The residuals are checked for convergence (using a predefined tolerance). A Jacobian matrix is constructed using
analyitcal derivatives of the terms in the residuals??,
J =
[
?R?/?? ?R?/?CT
?RCT /?? ?RCT /?CT
]
.
And finally the solutions are updated: [
?
CT
]
=
[
?
CT
]
? J?1
[
R?
RCT
]
.
Because of the square root in the ~V ?? equation, a CT > 1 will result in imaginary coefficients. The actual Newton step
is scaled down in this case (under-relaxed Newton). This generally occurs with high solidity and high TSR, since this
can result in CT very close to 1.
V. Optimization Results
The model was wrapped in an optimization loop in order to study the isolated affects of TSR and ?0(?) on
efficiency (CP ). The pitch variation is assumed sinusoidal (one-per-cycle variation) so
?0 = A0 + A1 sin(? + ?),
where A0 is the pitch offset, A1 is the oscillation amplitude, and ? is the phase. In addition, the affect of ?dirty
blades? was simulated with increasing factors multiplying cd. In all of the cases below, we report CP versus TSR,
with varying ?0 parameters or cd factors. All designs have three blades with NACA0012 airfoils, a reference Re=
3 ? 105, a streamtube deceleration factor? of ? = 6R, and optimized ?. In addition, the baseline case is optimized
for all parameters (with the cd factor Kcd = 1), and each figure represents sensitivity of the baseline to variations
in a single parameter. The baseline (optimized) case has ? = 0.29, TSR= 2.17, and a pitch modulation function
?0 = 7.6 + 6.2 sin(? + 21.5) degrees, which yields CP = 0.534.
1 2 3 4 5 6
0.1
0.2
0.3
0.4
0.5
0.6
Tip Speed Ratio
E
ffi
ci
en
cy
Varying pitch offset
A0=5
A0=7
A0=9
A0=11
A0=13
Figure 5. Efficiency vs tip speed ratio and pitch offset A0
1 2 3 4 5 6
4
5
6
7
8
9
10
11
Tip Speed Ratio
O
p
ti
m
a
l
p
it
ch
o
ff
se
t
(d
eg
)
Optimal pitch offset over TSR
Figure 6. Variation of optimal pitch offset A0 with tip speed ratio
First, we vary the pitch offset A0, while keeping the sinusoidal component constant and Kcd = 1. Fig. 5 shows that
there is an optimal A0, since the peak efficiency increases and then decreases with A0. In addition, as A0 increases,
the dependence on TSR increases. This is essentially the increased sensitivity of operating the turbine closer to stall
conditions. In addition, the peak efficiency shifts slightly toward lower TSR, where the pitch variation takes advantage
of the larger variation in blade-relative velocities. These trends can be seen in the variation of the optimal A0 with
TSR, as shown in Fig. 6.
??The effects of one main vortex on another are small and not included in calculating the derivatives due to complexity. This does not appear to
affect the convergence of the residuals.
?Chosen to agree with experimental results: D. W. Erickson, J. J. Wallace and J. Peraire (In preparation)
7 of 9
American Institute of Aeronautics and Astronautics
In Fig. 7, the pitch oscillation amplitude A1 is varied. Larger A1 tends to decrease the maximum efficiency point
(excessive stall) and shift it toward smaller TSR. Again, the dependence on TSR is stronger for higher A1 due to
operating closer to the stall limits.
1 2 3 4 5 6 7
?0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
Tip Speed Ratio
E
ffi
ci
en
cy
Varying pitch angle amplitude
A1=3
A1=5
A1=7
A1=9
A1=11
Figure 7. Efficiency vs tip speed ratio and pitch modulation ampli-
tude A1
1 2 3 4 5 6 7
?0.1
0
0.1
0.2
0.3
0.4
0.5
0.6
Tip Speed Ratio
E
ffi
ci
en
cy
Varying pitch angle phase
?=0
?=20
?=40
?=60
?=80
Figure 8. Efficiency vs tip speed ratio and pitch modulation phase
angle ?
In Fig. 8, the phase ? is varied. As can be seen, the efficiency is somewhat sensitive to the phase angle. Phase
lead reduces the peak efficiency and reduces off-design performance (i.e. efficiency at higher TSR). Interestingly, this
effect is very weak at low TSR.
Finally, variation of the cd factor is shown in Fig. 9. Clearly, larger cd reduces efficiency. In addition, larger Kcd
shifts the optimal TSR downward, since operating at a higher TSR is more sensitive to drag.
Note that in all of the above cases, ? = c/R (solitidy) is optimized. This is done because there is a relationship
between ? and TSR in the optimal case, specifically ?TSR ?constant. Since TSR is a ratio of velocities and ? is
a ratio of length scales, ?TSR is a ratio of time scales between the turbine (?) and the freestream (c/~V?). In an
optimized design, the velocity, length, and time scales are all matched. Fig. 10 shows this trend for the baseline case
(only varying TSR).
1 2 3 4 5 6
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
Tip Speed Ratio
E
ffi
ci
en
cy
Varying Cd scale
K=1
K=1.5
K=2
Figure 9. Efficiency vs tip speed ratio and drag coefficient factor
Kcd
1 2 3 4 5 6
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
Tip Speed Ratio
O
p
ti
m
a
l
so
li
d
it
y
Optimal solidity over TSR
Figure 10. Variation of optimal solitidy ? with tip speed ratio
8 of 9
American Institute of Aeronautics and Astronautics
VI. Future Work
Further improvements to the model would include accounting for dynamic stall and Reynolds number variation
over the angle ?. The model converges well for moderate to low ?, moderate to high TSR, and moderate ?0. Dynamic
stall modeling may give better accuracy at high ? and low TSR, but the model has been shown to be a good approxima-
tion to a high fidelity 2D Navier-Stokes solution with high solidity and low tip speed ratio. The high fidelity solution
has vortices shedding at the airfoil leading edge (dynamic stall), which produces fluctuations in the resulting ~Ftot(?)
and ~Mtot(?). The low order model predicts forces and moments that are close to the high fidelity results in regions
without dynamic stall, and the overall CP is similar. Reynolds number variation over ? would have to be implemented
as an empirical cd scaling and would not change the model results much for moderate to high TSR.
This model is particularly easy to extend to various pitch control mechanisms. First of all, an optimal ?0(?)
variation can be computed by wrapping the model in an optimization loop. This can be compared to forced one-per-
cycle pitch variation, where ?0(?) = A0 + A1 sin(? + ?). One could also include passive pitch control mechanisms
such as a spring or an aerodynamic spring. This is simply done by introducing new state variables ?0(?) and ??0(?)
and calculating the appropriate additional velocities and forces.
References
1Templin, R. J., ?Aerodynamic Performance Theory for the NRC Vertical-Axis Wind Turbine,? NASA STI/Recon Technical Report N, Vol. 76,
June 1974, pp. 16618?+.
2Strickland, J. H., ?A Performance Prediction Model for the Darrieus Turbine,? International Symposium on Wind Energy Systems., British
Hydrodynamics Research Association, Cranfield, Beds., England, 1977, pp. 3?+.
3Paraschivoiu, I., ?Double-Multiple Streamtube Model for Studying Vertical-Axis Wind Turbines,? Journal of Propulsion Power, Vol. 4, Aug.
1988, pp. 370?377.
4Strickland, J. H., Webster, B. T., and Nguyen, T., ?Vortex Model of the Darrieus Turbine: An Analytical and Experimental Study,? NASA
STI/Recon Technical Report N, Vol. 80, Feb. 1980, pp. 25887?+.
5Drela, M., ?Integrated Simulation Model for Preliminary Aerodynamic, Structural, and Control-Law Design of Aircraft,?
AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference and Exhibit, Vol. 3, AIAA, Washington, DC, 1999.
6Drela, M., ?XFOIL: An Analysis and Design System for Low Reynolds Number Airfoils,? Conference on Low Reynolds Number Airfoil
Aerodynamics, University of Notre Dame, June 1989.
7Gopinath, A. K. and Jameson, A., ?Time Spectral Method for Periodic Unsteady Computations over Two- and Three- Dimensional Bodies,?
43rd AIAA Aerospace Sciences Meeting and Exhibit, AIAA, Washington, DC, Jan. 2005.
9 of 9
American Institute of Aeronautics and Astronautics