вход по аккаунту



код для вставкиСкачать
Numerical and Experimental Investigation of Wind
Turbine Wakes
X. Huang∗1 , S. Vey†2 , M. Meinke1 , W. Schröder1 , G. Pechlivanoglou2 , C. N. Nayeri2 , and
C. O. Paschereit2
Institute of Aerodynamics, RWTH Aachen University, Wüllnerstr. 5a, 52062 Aachen, Germany
Institute of Fluid Dynamics and Technical Acoustics, TU Berlin, Müller-Breslau-Str. 8, 10623 Berlin
Numerical and experimental investigations of the flow characteristics around the wind
turbine blade and the wake flow field are performed. The numerical method is based on
a large-eddy simulation formulated in a rotating reference frame on a structured bodyfitted grid to predict the flow field around a complete rotating wind turbine blade at an
average blade Reynolds number of 300,000. Periodic boundary conditions are applied in the
circumferential direction, such that the flow over only one of the three blades of the wind
turbine is computed. The numerical results show a flow separation on the suction side on
the blade with a subsequent transition to turbulence and flow reattachment. Wind tunnel
measurements are performed on a 3 m diameter wind turbine installed in the TU Berlin
large wind tunnel. Initial experimental data were measured using the data acquisition
system installed inside the rotating hub of the wind turbine. First results of the inflow
angle and blade root strain gauge measurements at varying inflow angles are presented.
To efficiently exploit available installation space and to reduce maintenance and installation costs, wind
turbines are clustered in large wind farms. Therefore, a substantial number of wind turbines are located in
the wake of upstream wind turbines. The velocity deficit in the wake will lead to a reduction of the power
production and the high turbulence intensity will increase dynamic loads on the downstream wind turbine.
Depending on the location and wind direction, the power loss of one downstream wind turbine can reach
40%. If the data are averaged over the various wind directions, approx. 8% power reduction is observed for
the onshore and 12% for the offshore installations [1].
The wake of a wind turbine is usually divided into two regions, the near wake and the far wake. In the
near wake up to one or two diameters downstream of the rotor, the influence of the rotor geometry can be
directly identified in the form of the tip and root vortex. The tip and root vortex can cause large velocity
gradients and turbulence intensities and even form shear layers in the case of high tip-speed ratios. In the
far wake, the influence of the detailed rotor geometry is no more detectable, but the tip and root or hub
vortices still have a significant influence on the turbulence intensity. A better understanding of the tip and
root vortex dynamics can help to accelerate the recovery of the velocity deficit and decrease the overall
turbulence intensity in the far wake.
To investigate the development of the wake downstream of a turbine blade including the vortex system,
numerical simulations and wind tunnel experiments have been conducted. Numerical methods, based on
Reynolds averaged Navier-Stokes (RANS) equations including the k-ω turbulence model, have been widely
used for wind turbine simulations [2, 11, 18] because of its robustness and accuracy for most turbulent flows.
It is well known, however, that RANS solutions for the tip and root vortex are strongly dependent on the
turbulence models. Compared to RANS, large-eddy simulation methods resolve most of the turbulence energy
containing structures and the modeling approach only focuses on the subgrid scale of the turbulent spectrum.
∗ PhD Student, Institute of Aerodynamics, RWTH Aachen University, Wüllnerstr. 5a, 52062 Aachen,, Tel.: +49 241 80 90404
† Dr.-Ing, Institute of Fluid Dynamics and Technical Acoustics, TU Berlin, Müller-Breslau-Str.
8, 10623 Berlin,, Tel.: +49 30 314 22693
1 of 12
American Institute of Aeronautics and Astronautics
LES has been proven to be a successful approach to simulate unsteady flows around airfoils and there are
applications of LES coupled with actuator line techniques to wind turbine wake simulations [5, 6, 10, 17].
In this paper, the turbulent flow around the wind turbine blade and hub at a Reynolds number of
Rec = 300, 000 based on the mean chord will be investigated by LES. The flow field around the blade and
the near wake is predicted on a block-structured body-fitted mesh, to ensure an adequate resolution of the
boundary layer on the blade. The emphasis of the investigation is put on the accurate prediction of the tip
and hub vortices in the wake which cannot be accurately captured by RANS methods.
The experimental data are obtained on the Berlin Research Turbine (BeRT). The BeRT is a load control
research turbine that is equipped with sensor and actuator systems installed in the rotating system. Chordwise and spanwise pressure distributions, blade root bending moment, tip acceleration, and inflow angle and
velocity are captured at a sample rate of 10 kHz. The dynamic flow effects on the rotor blade can be studied
under yawed inflow angle of −30◦ ≤ Θyaw ≤ 30◦ .
Governing Equations
The Navier-Stokes equations are transformed into a rotating frame of reference, where the relative velocity
vr is defined by the difference of the absolute velocity →
va and the rotational velocity →
ve ,
vr = →
va − →
For simplicity, the x-coordinate direction is assumed to be the rotation axis, hence, the Navier-Stokes equations in the relative frame of reference for the unsteady compressible flow can be written as
dV + F · ndA = SdV
where n is the normal vector of the surface dA and Q = [ρ, ρu, ρE] is the vector of conservative variables.
The quantities ρ, u, E denote the density, the velocity vector in the relative frame and the relative total
energy. The relative total energy is defined as
E =e+
vr 2
ve 2
u2 + v 2 + w2
ω2 →
where ω is the rotation speed and r is the radius of the rotation axis. The flux vector F contains an inviscid
Fi and a viscous part Fv
1 
F = Fi + Fv =  ρuU + p · n  +
τu + q
The Reynolds number is defined by the fluid properties at rest state and p is the static pressure. The relative
total enthalpy in a steadily rotating frame of reference is defined as
p →
vr 2
ve 2
U = nx u + ny v + nz w,
H =e+
U is the contravariant velocity
where nx , ny , nz represent components of the normal vector of the surface A. The viscous stress τ can be
expressed by
τ = η (∇ · u) I − 2η ∇u + (∇u)
where η is the dynamic viscosity computed by Sutherland’s law. The vector of heat conduction q is
Pr (γ − 1)
2 of 12
American Institute of Aeronautics and Astronautics
The quantity k is the thermal conductivity and the Prandtl number is defined as Pr = η∞ cp /k∞ , where cp
is the specific heat at constant pressure and γ the ratio of specific heats.
The source term S contains the Coriolis and centrifugal force
ρω(yω + 2w) + ρfy
ρω(zω − 2v) + ρfz
ρf vr + q̇h
with f denoting the body force and q̇h being the time rate of heat transfer. The system of equations is closed
by the equation of state for a perfect gas
ρ (γ − 1)
Numerical Methods
The Navier-Stokes equations are spatially filtered assuming an implicit grid filter and discretized using a
finite-volume technique on a structured body-fitted grid. The monotonic integrated large-eddy simulation
(MILES) approach proposed in [3] is used, where the dissipative part of the truncation error is assumed
to mimic the effect of the dissipation of the unresolved subgrid scales. Various difference stencils are used
to discretize the convective and viscous fluxes. The convective fluxes are formulated by a low dissipation
variant of the AUSM scheme [8] and a MUSCL interpolation at second-order accuracy. The discretization of
the viscous fluxes is performed by a second-order accurate central difference scheme [9]. An explicit five-step
Runge-Kutta formulation also at second-order accuracy is used for the temporal integration. This solution
method has been extensively validated for various turbulent flow problems, see e.g. [7, 13, 14].
Computational Setup and Grid Strategy
The wind turbine geometry is identical to that used in the wind tunnel of TU Berlin. The rotor uses
the downscaled geometry from the previous European project MEXICO. The Clark-Y profile is used for the
blade from the root to the tip. Experimental measurements are performed in the wind tunnel of TU Berlin
with a downscaled rotor diameter of approx. 3m and the hub height of approx. 2.1m. The Reynolds number,
based on the freestream velocity (6.4m/s) and radius of the rotor (1.5m), is 640,000 and the tip-speed ratio
is equal to 4.6, which matches the experimental value.
Figure 1. Domain size and boundary conditions for the LES of the flow around the wind turbine blade
3 of 12
American Institute of Aeronautics and Astronautics
Fig. 1 shows the geometric setup and boundary conditions in the CFD model of the wind turbine blade.
To reduce the computational effort, the simulation domain around the wind turbine will be restricted to a 120
degree section in the circumferential direction such that the flow around one blade has to be resolved. Periodic
boundary conditions are used in the circumferential direction as shown in the Fig. 1. The computational
domain has a length of 10R, where R is the radius of the rotor, upstream and downstream of the turbine
blade to avoid an influence of the approximate boundary conditions. At the inflow boundary the mean inflow
velocity is prescribed.
Figure 2. Mesh and topology near the blade
Figure 3. Mesh and topology for the periodic planes in the circumferential direction
A structured body-fitted grid with 127 million mesh points around the blade and hub was created with
a commercial mesh generation software. The mesh resolution required by the LES of wall-bounded flows
(∆y + ≈ 2, ∆z + ≈ 30, ∆x+ ≈ 50) is reached in the blade-tip region. The grid system consists of 48 blocks
covering one blade and one third of the hub. To obtain a better grid orthogonality near the blade surface, a
C-type topology is used along the span. Due to the periodic boundary conditions, the number of cells and
the distributions upstream and downstream of the blade must be identical such that the special topology
depicted in Fig. 2 was designed to meet the requirements. Besides, additional blocks, indicated in blue in
Fig. 3, are inserted as buffer zone between the blade blocks and periodic surfaces, which can reduce the mesh
skewness near the blade and improve the mesh resolution of the near wake.
4 of 12
American Institute of Aeronautics and Astronautics
Numerical Results
In Fig. 4, the contours of the second invariant of the rate of strain tensor is used to illustrate the vortical
structures around the wind turbine blade. A stable tip vortex and turbulent structures in the boundary
layer and the wake can be seen. The axial velocity distribution in the wake flow on a vertical plane through
the blade is shown in Fig. 5. The higher velocity in the blade tip region and the tip vortex are evident. The
wake deficit downstream of the wind turbine blade is visible in the lower values of the axial velocity.
(a) Pressure side
(b) Suction side
Figure 4. Q-Criterion of vortical and turbulent structures on the pressure side (a) and suction side (b) of the
Figure 5. Axial velocity contours on a cut plane parallel to the rotation axis and through the turbine blade.
The instantaneous and mean pressure contours on the suction and pressure surfaces of the blade are
shown in Fig. 6. A plateau and sudden decline of the static pressure is visible on the suction side of the
blade indicating flow separation. In Fig. 6, I, II and III denote the slices at 63%, 80%, and 90% span. Figs.
7 and 8 show the averaged chordwise pressure coefficient and the skin-friction coefficient against φ at 63%,
80%, and 90% span of the blade (I, II, and III as shown in Fig.6). The quantity φ is defined as the arcsin
(L/r), where L is the distance between the point on the blade and XY plane and r the distance to the X-axis.
The pressure coefficient and skin-friction coefficient are computed by
5 of 12
American Institute of Aeronautics and Astronautics
Figure 6. (a) Instantaneous pressure contours on the suction and pressure surface; (b) time-averaged pressure
contours on the suction and pressure surface.
Figure 7. Pressure coefficient around the blade at 63%(a), 80%(b), and 90%(c) span of the turbine blade
Figure 8. Skin-friction coefficient around the blade at 63%(a), 80%(b), and 90%(c) span of the turbine blade
Cp =
(P − P∞ )
+ (ωr)2 ])
2 ρ[U∞
6 of 12
American Institute of Aeronautics and Astronautics
Figure 9. Mach number contours on the slices at 63%, 80%, and 90% span of the turbine blade
Figure 10. Contours of turbulent kinetic energy on the slices at 63%, 80%, and 90% span of the turbine blade
Cf =
2 ρ[U∞
+ (ωr)2 ])
In Figs. 7 and 8, the pressure distribution on the suction side also show the plateau, while the skin-friction
coefficient reaches a negative value at the same position showing the flow to separate and reattach further
downstream. The instantaneous contours of the relative velocity magnitude and the turbulent kinetic energy
around the blade section at 63%, 80%, and 90% span is shown in Fig. 9 and 10. It can be seen that the
separated flow region is small and a turbulent boundary layer develops downstream of the reattachment.
Research Turbine and Experimental Setup
A 3 m diameter modular load control research wind turbine was designed, built, and installed in the TU
Berlin large wind tunnel wind energy measurement section. The ”Berlin Research Turbine”, BeRT is the
first German load control research turbine. A further unique feature for this turbine size is the installation of
data acquisition and real time load control processor within the rotating system. With this setup a multitude
of sensor signals can be read simultaneously at a sampling frequency of up to 10 kHz. The rotor blades can
be easily interchanged so that various load control principles can be tested experimentally.
The setup in the wind energy measurement box (WOX) is illustrated in Fig. 11. The existing wind
tunnel’s settling chamber was extended in two stages (WOX 1 and WOX 2, respectively) by a total of 5 m.
This new test section has a cross-section of 4.2 m x 4.2 m. It is mainly used for measurements incorporating
the BeRT and features a 1.1 m diameter turntable with eccentric mounting for the turbine tower so that
yaw movements do not interfere with blade tip – side wall distance symmetry. The turbine can be set to
yaw angles in the range −30◦ ≤ Θyaw ≤ 30◦ .
The baseline turbines power curve is shown in Fig. 12. Maximum power output is reached at a tip speed
ratio of TSR=4. At a maximum inflow velocity of U∞ = 8 m/s the blade Reynolds number is in the range
250 000 ≤ Re ≤ 350 000. Note that the measurements were performed with only WOX 1 installed. This
means that the wind tunnel nozzle is located 0.5 rotor diameters downstream of the rotor-plane and therefore
7 of 12
American Institute of Aeronautics and Astronautics
an interaction with the wake-development of the turbine cannot be ruled out. WOX 2 is now installed, so
that future measurements will be performed with the nozzle further downstream so that the interaction
between wake and nozzle will be reduced.
The blade is equipped with a Clark Y airfoil along the full span – a cylindrical root section is omitted.
This airfoil was chosen because of its suitability for flaps, while still being relatively insensitive to Reynolds
number effects. Furthermore, it provides the necessary thickness to house the required sensors and actuators
inside the blade.
Figure 12. Power curve of BeRT.
Figure 11. Setup in wind tunnel.
The initial set of three baseline blades was used in the development and validation of the data acquisition
and load control system. The baseline blades were therefore equipped with stick-on pressure sensors, threehole probes, blade-root strain gauges, and a tip accelerometer. These sensors are representative of the sensors
that will be integrated into the active blade. With these preliminary measurements some problems could be
identified and fixed at an early-on stage, such as an insufficient stiffness of the turbine foundations. The fully
sensor- and actuator-equipped rotor blade has just been commissioned and measurements will be performed
in the near future. A picture of the BeRT with the sensor- and actuator-equipped rotor-blade is shown in
Fig. 13.
Figure 13. BeRT with smart rotorblade.
Figure 14. View of BeRT with flow tufts.
Measurement Techniques and Results
Quantitative Tuft Flow Visualization
The BeRT presented the unique opportunity for the development of a novel measurement technique where
a quantitative tuft flow visualization technique was synchronized with time resolved pressure and vibration
8 of 12
American Institute of Aeronautics and Astronautics
measurements. With this technique, an arbitrary measured variable (e.g. pressure) can be linked to an
instantaneous surface flow field on the rotor blades. It is possible to capture the complete rotor in one
image, with all three blades equipped with flow tufts and image registration markers. A dedicated in-house
code then distinguishes between the three blades and computes instantaneous tuft-angles for each tuft in
each image. The resulting data-set is then processed to characterize the flow regime at an individual tuft
location. In combination with e.g. the relative angle of attack reading from the 3-hole probe and the root
strain gauge a deeper insight into the response of the turbine to dynamic inflow conditions can be gained.
Variation of Yaw Angle
Some interesting observations can be made from the development of the measured variables under various
yawed inflow cases. An example of which is shown for the strain gauge (top) and surface pressure (bottom)
signals in Fig. 15. The surface pressure signal was measured using a stick-on surface pressure port at
x/c=75% and a radial position of 75% blade length. It shows the expected behavior, namely that signal
amplitude is increased for larger deviations from Θyaw = 0◦ . In contrast to the pressure signals, the strain
gauge signal shows no material change in signal amplitude when the yaw angle is altered. This behavior
was not expected and the source was identified to be a natural frequency overlap caused by an insufficient
stiffness of the wind turbine foundation. This causes turbine vibrations which are directly transmitted to the
turbine blade root. Once the source of the problem was identified the base plate was replaced with a stiffer
construction that incorporates a central stiffening element. Initial qualitative vibration tests were promising
and the improvement will be quantified in the next measurement campaign.
yaw angle
blade root strain gage
10% c pressure
75% radial position
Figure 15. Signals from blade root strain gauge and 10 % chord, 75 % span stick-on surface pressure port.
Note the change in pressure signal amplitude over yaw angle while the strain gauge amplitude remains the
Three-Hole Probe Measurements
One of the tested load control schemes requires the accurate and time resolved measurement of the relative
blade section inflow angle. This is accomplished using three hole probes mounted to the respective blade
section. Shown in Fig. 16 is the orientation of the probe itself and the pressure signal from the three probe
holes for various yaw-cases. These pressure signals will be used directly as an input to a machine learned
control scheme (MLC). In the future, a calibration function will be applied and both the relative inflow
speed and direction (i.e. the angle of attack) can be derived from the raw signals. Furthermore, the number
of required pressure transducers can be reduced from currently three to two per three hole probe. This is
accomplished by measuring against the central three hole probe tube.
Results from Quantitative Tuft Flow Visualization Technique
It was already mentioned above that a quantitative tuft flow visualization technique was combined and
synchronized with the time resolved measurement system. The setup is shown in Fig. 17. The turbine was
set to a yaw angle of 30◦ and operated at a tip speed ratio of T SR = 4. All three blades were equipped with
flow tufts and image registration markers (see also the photo in Fig. 14). A high resolution digital camera
was mounted downstream of the turbine in the high speed test section of the wind tunnel. With the high
camera-resolution it was possible to capture the whole rotor plane in a single image.
9 of 12
American Institute of Aeronautics and Astronautics
Figure 16. Raw signals from three hole probe. Signals are used as an input to the AFC control loop.
flow tufts
Θ yaw
strain gage
flash detector
2 kW flash-lamp
Figure 17. Experimental setup for tuft flow visualization measurements.
Figure 18. Surface flow patterns obtained with
phase averaged quantitative tuft flow visualization.
A high-energy, high-speed flash unit provided the necessary lighting. The light pulse was detected by
the flash detector, whose signal was recorded simultaneously with the other measured variables by the
data acquisition system. In a post processing step the individual images can then be associated with the
instantaneous measured variables.
The raw tuft images were then post-processed using the in-house FlowViz code. The resulting per-blade
vector fields were summarized to yield per-tuft statistics. These can then be further processed using tools
and techniques known from PIV data post-processing such as phase- and ensemble-averaging and POD. This
technique has already been applied to utility scale wind turbines [15], and full scale test cars and trucks.
The FlowViz-technique is under continuous development because surprisingly the evaluation of wind
tunnel FlowViz-data poses the larger effort compared to the free-field measurements. The reason is the
better lighting conditions in the free field measurements. It is a non trivial task to create acceptable lighting
conditions in a 4.2 m by 4.2 m cross section with black walls and a rotating object with curved surfaces.
Some preliminary results of this exciting new measurement possibility are shown in Fig. 18. The blade’s
vector fields were phase averaged for Θwing =270◦ , 0◦ , and 90◦ . The top part of Fig. 18 shows the average
10 of 12
American Institute of Aeronautics and Astronautics
Figure 19. Mean surface flow angle over phase angle. Data extracted from full flowviz-dataset.
deviation ∆Φ between the chord-wise and the local flow vectors. The lower half of the figure shows the
streamlines visualized using a line integral convolution (lic) [4] with a visual notion of the streamline deviation. The sequence shows the periodic behavior of the crossflows of the BeRT under Θyaw = 30◦ . The data
is again summarized in Fig. 19. Here, each data point represents the mean flow direction of a complete blade
at a certain rotor azimuthal angle Θwing . In the azimuthal rotor angle range of −200◦ ≤ Θwing ≤ −125◦ the
turbine’s tower is blocking sight. This data can now be correlated to the pressure or strain gauge signals
shown in Fig. 15.
Conclusion and Outlook
This study presents a high fidelity flow simulation for a wind turbine blade. The approach uses a
large-eddy simulation solver to investigate the aerodynamic characteristics around wind turbine blade. The
turbulence around the blade and the tip vortex in the near wake is well captured by the LES.
For the numerical part, the future work will focus on a detailed comparison of the numerical and experimental data and the analysis of the stability of the tip vortex system in the wind turbine wake [12, 16]. Two
oscillating flaps will be installed close to the wind turbine blade tip and the mesh deformation method will
be used to mimic the motion of the flaps. The influence of the perturbation frequency of the flaps on the
stability of the vortex system as well as the aerodynamic loads on the blade will be investigated.
Initial results from an experimental wind tunnel study on the Berlin Research Turbine (BeRT) were
presented for varying yaw angles. A quantitative tuft flow visualization technique was conducted and the
variations in surface flow patterns on the rotor blades under a yaw angle of Θyaw = 30◦ were shown. In
future experiments active trailing edge flaps will be used in an attempt to excite tip vortex instabilities. The
tip vortex structure and dynamics will then be studied using stereoscopic particle image velocimetry.
This numerical and experimental research is conducted in the DFG funded project Active control of wind
turbine wakes PAK 780. All computations were performed at the High Performance Computing Center
Stuttgart (HLRS).
11 of 12
American Institute of Aeronautics and Astronautics
1 R. Barthelmie and L. Jensen. Evaluation of wind farm efficiency and wind turbine wakes at the nysted offshore wind
farm. Wind Energy, 13(6):573–586, 2010.
2 A. Bechmann and N. Sorensen. Cfd simulation of the mexico rotor wake. In European wind energy conference and
exhibition, Marseille, 2009.
3 J. P. Boris, F. F. Grinstein, E. S. Oran, and R. L. Kolbe. New insights into large eddy simulation. Fluid Dynamics
Research, 10:199–228.
4 B. Cabral and L.C. Leedom. Imaging vector fields using line integral convolution. In Proceedings of ACM SIGGRAPH
’93, pages 263–270, 1993.
5 S. Ivanell. Numerical analysis of the tip and root vortex position in the wake of a wind turbine. 2007. Journal of Physics:
conference series, 75-012035.
6 A. Jimenez. Application of a les technique to characterize the wake deflection of a wind turbine in yaw. Wind Energy,
13:559–572, 2010.
7 M. Konopka, M. Meinke, and W. Schrder. Large-eddy simulation of shock-cooling-film interaction at helium and hydrogen
injection. Phys.Fluids, 25(10),106101, 2013.
8 M. S. Liou and C. J. Steffen-Jr. A new flux splitting scheme. J. Comput. Phys., 107:23–39, 1993.
9 M. Meinke, W. Schröder, E. Krause, and Th. Rister. A comparison of second and sixth-order methods for large-eddy
simulations. Comp. Fluids, 31:695–718, 2002.
10 J Meyers and C Meneveau. Large eddy simulations of large wind-turbine arrays in the atmospheric boundary layer. In
48th AIAA Aerospace Sciences Meeting, 2010. AIAA 2010-827.
11 N. Sorensen and M.O.L. Hansen. Rotor performance predictions using navier-stockes method. In 36th AIAA Aerospace
Sciences Meeting and Exhibit, 1998. AIAA-98-0025.
12 V. L. Okulov and J. N. Sorensen. Stability of helical tip vortices in a rotor far wake. J. Fluid Mech., 576:1–25, 2007.
13 P. Renze, W. Schröder, and M. Meinke. Large-eddy simulation of film cooling flows at density gradients. International
Journal of Heat and Fluid Flow, 29:18–34, 2008.
14 V. Statnikov, T. Sayadi, M. Meinke, P. Schmid, and W. Schrder. Analysis of pressure perturbation sources on a generic
space launcher after-body in supersonic flow using zonal rans/les and dynamic mode decomposition. Phys.Fluids, 27,106101,
15 S. Vey, H.M. Lang, C.N. Nayeri, C.O. Paschereit, and G. Pechlivanoglou. Extracting quantitative data from tuft flow
visualizations on utility scale wind turbines. In The Science of Making Torque from Wind 2014, Journal of Physics: Conference
Series 524, 2014. 012011.
16 S.E. Widnall. The stability of a helical filament. J. Fluid Mech., 54(4):641–663, 1973.
17 Y Wu and Porte-Agel F. Large-eddy simulation of wind turbine wakes: evaluation of turbine parameterization. BoundaryLayer Meteorology 2011, 138.
18 F. Zahle and N. Sorensen. Wind turbine rotor-tower interaction using an incompressible overset grid method. Wind
Energy, 12:594–619, 2009.
12 of 12
American Institute of Aeronautics and Astronautics
Без категории
Размер файла
8 566 Кб
2015, 2310
Пожаловаться на содержимое документа