Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • am-dept/matlab_airfoil_toolbox
1 result
Show changes
Commits on Source (13)
Hide whitespace changes
Inline Side-by-side
+af_tools/+utils/spacedvector.m
View file @ 0995fa2a
function vect = spacedvector(spacing, nPoints)
% SPACEDVECTOR Returns a vector based on chosen spacing and number of points
%
% Spacing:
% - linear : points spaced linearily
% - halfcosine : more dense at beginning and end of vector
% - cosine : more dense at the beginning and more spaced out at the end
% - sine : more dense at the end and more spaced out at the beginning
%
% -----
% (c) Copyright 2022 University of Liege
% (c) Copyright 2022-2023 University of Liege
% Author: Thomas Lambert <t.lambert@uliege.be>
% ULiege - Aeroelasticity and Experimental Aerodynamics
% Apache 2.0 License
......@@ -18,6 +25,9 @@ function vect = spacedvector(spacing, nPoints)
case 'cosine'
beta = linspace(0, pi / 2, nPoints);
vect = 1 - cos(beta);
case 'sine'
beta = linspace(0, pi / 2, nPoints);
vect = sin(beta);
end
end
+af_tools/@Airfoil/Airfoil.m
View file @ 0995fa2a
......@@ -99,4 +99,9 @@ classdef Airfoil < handle
end
end
methods (Static)
% Create an airfoil by interpolating between two different Airofil objects
self = interpairfoil(Airfoil1, Airfoil2, weightAirfoil2)
end
end
+af_tools/@Airfoil/interpairfoil.m 0 → 100644
View file @ 0995fa2a
function self = interpairfoil(Airfoil1, Airfoil2, weightAirfoil2)
% INTERPAIRFOIL Create a new Airfoil object by interpolation of two other Airfoils.
% This static method can be used in place of a normal constructor for the Airfoil class.
% This method can be used when the user wants to create a non-standard airfoil made from
% the interpolation bewteen two known airfoils.
% Consider that the airfoil at x=0 is a NACA 0012, and the one at x=1 is a CLARY-Y.
% Rather than assuming that all segments in [0,1] are NACA 0012 and all segments after 1 are
% CLARK-Y, this function allows to create polars for the [0,1] interval that are a mix between
% the two bounds. If we are at 0.1, the newly created Airfoil object, will then be 90% of NACA
% 0012 and 10% of CLARK-Y.
% Note:
% The two input Airfoils should contain Polars of similar structure
% (same Reynolds, Mach and nCrit).
% -----
%
% Usage:
% MixedAirfoil = Airfoil.INTERPPOLAR(Airfoil1, Airfoil2, weightAirfoil2) create an new Airofil
% object that is a linear interpolation between Airfoil1 and Airfoil2, with a given weight for
% Airfoil2.
%
% Inputs:
% Airfoil1 : Airfoil object for one part of the interpolation.
% Airfoil2 : Airfoil object for the second part of the interpolation.
% Airfoil2.Polar must have the same Re, Mach and nCrit values as
% Airfoil1.Polar.
% weightAirfoil2 : Scalar in [0-1] to skew the interpolation properly between Airfoil1 and
% Airfoil2.
%
% Output:
% self : Airfoil object that results from the interpolation.
%
% Example:
% MixedAirfoil = af_tools.interpairfoil(AfNACA0012, AfClarkY, 0.6); Creates an interpolated
% Airfoil object MixedAirfoil that results from the interpolation of AfNACA0012
% and AfClarkY, asssuming that the MixedPolar is made of 60% ClarkY and 40%
% NACA0012.
%
% See also: af_tools.Airfoil, af_tools.Polar.
%
% -----
% <a href="https://gitlab.uliege.be/am-dept/matlab_airfoil_toolbox">Documentation (online)</a>
% ----------------------------------------------------------------------------------------------
% (c) Copyright 2022-2024 University of Liege
% Author: Thomas Lambert <t.lambert@uliege.be>
% ULiege - Aeroelasticity and Experimental Aerodynamics
% MIT License
% Repo: https://gitlab.uliege.be/am-dept/matlab_airfoil_toolbox
% Issues: https://gitlab.uliege.be/am-dept/matlab_airfoil_toolbox/-/issues
% ----------------------------------------------------------------------------------------------
% --- Defaults and constants
DEBUG = false;
% --- Input checks
% Check if both Airfoils are Airfoil objects and weightAirfoil2 is in [0,1]
validateattributes(Airfoil1, {'af_tools.Airfoil'}, {'scalar'}, mfilename(), 'Airfoil1', 1);
validateattributes(Airfoil2, {'af_tools.Airfoil'}, {'scalar'}, mfilename(), 'Airfoil2', 2);
validateattributes(weightAirfoil2, {'double'}, {'scalar', 'nonnegative', '<=', 1}, ...
mfilename(), 'weightAirfoil2', 3);
% --- Interpolate the Airfoils
self = af_tools.Airfoil;
if Airfoil1 == Airfoil2
self.name = Airfoil1.name;
self.coord = Airfoil1.coord;
self.upper = Airfoil1.upper;
self.lower = Airfoil1.lower;
self.Polar = Airfoil1.Polar;
return
end
% Private properties to prevent accidental re-assignment
wAf2 = weightAirfoil2;
wAf1 = 1 - wAf2;
% Set base Airfoil information
self.name = sprintf('MIX-%d%%%s-%d%%%s', round(wAf1 * 100), Airfoil1.name, ...
round(wAf2 * 100), Airfoil2.name);
% Interpolate coordinates and Polars
self.upper = interpsurf(wAf1, Airfoil1.upper, wAf2, Airfoil2.upper);
self.lower = interpsurf(wAf1, Airfoil1.lower, wAf2, Airfoil2.lower);
self.coord = [flipud(self.upper); self.lower(2:end, :)];
self.Polar = af_tools.Polar.interppolar(Airfoil1.Polar, Airfoil2.Polar, wAf2);
% --- DEBUGGING
if DEBUG
% Plot the original Airfoil Coordinates as well as the interpolated one
figure('Name', 'Interpolated Airfoil');
hold on;
plot(Airfoil1.coord(:, 1), Airfoil1.coord(:, 2));
plot(Airfoil2.coord(:, 1), Airfoil2.coord(:, 2));
plot(self.coord(:, 1), self.coord(:, 2));
hold off;
legend(Airfoil1.name, Airfoil2.name, self.name, 'Location', 'NorthWest');
grid on;
pause;
end
end
function newSurf = interpsurf(weight1, surfCoords1, weight2, surfCoords2)
% INTERPSURF Interpolate airfoil surfaces (ONE SURFACE AT A TIME, NOT FULL COORDS)
% Reinterpolate original data based on the longest chord vector
if numel(surfCoords1(:, 1)) > numel(surfCoords2(:, 1))
newX = surfCoords1(:, 1);
newY1 = surfCoords1(:, 2);
newY2 = interp1(surfCoords2(:, 1), surfCoords2(:, 2), newX);
else
newX = surfCoords2(:, 1);
newY1 = interp1(surfCoords1(:, 1), surfCoords1(:, 2), newX);
newY2 = surfCoords2(:, 2);
end
newY = weight1 * newY1 + weight2 * newY2;
newSurf = [newX, newY];
end
+af_tools/@Polar/Polar.m
View file @ 0995fa2a
......@@ -145,7 +145,7 @@ classdef Polar < handle
%
% See also: af_tools.findcllinearrange.
[self.Lin.clRange, self.Lin.aoaRange, self.Lin.slope] = ...
[self.Lin.aoaRange, self.Lin.clRange, self.Lin.slope] = ...
af_tools.findcllinearrange(self);
end
......@@ -180,6 +180,8 @@ classdef Polar < handle
self = loadpolar(self, polarFile) % Load polar from file
self = polypolar(polyCl, polyCd) % Create polar using polynomial coefficients
% Create a polar by interpolating between two different polar objects
self = interppolar(Polar1, Polar2, weightPolar2)
end
end
+af_tools/@Polar/interppolar.m 0 → 100644
View file @ 0995fa2a
function self = interppolar(Polar1, Polar2, weightPolar2)
% INTERPPOLAR Create a new polar object by interpolation of two other polars.
% This static method can be used in place of a normal constructor for the Polar class.
% This method can be used when the user wants to interpolate the polars bewteen two known
% airfoils. Consider that the airfoil at x=0 is a NACA 0012, and the one at x=1 is a CLARY-Y.
% Rather than assuming that all segments in [0,1] are NACA 0012 and all segments after 1 are
% CLARK-Y, this function allows to create polars for the [0,1] interval that are a mix between
% the two bounds. If we are at 0.1, the newly created Polar object, will then be 90% of NACA
% 0012 and 10% of CLARK-Y.
% Note:
% The two input Polars should have the same range of Reynolds, Mach, and nCrit.
% -----
%
% Usage:
% MiedPolar = Polar.INTERPPOLAR(Polar1, Polar2, weightPolar2) create an new Polar object that
% is a linear interpolation between Polar1 and Polar2, with a given weight for Polar2.
%
% Inputs:
% Polar1 : Polar object for one part of the interpolation.
% Polar2 : Polar object for the second part of the interpolation. Muse have the same Re, Mach
% and nCrit values as Polar1.
% weightPolar2: Scalar in [0-1] to skew the interpolation properly between Polar1 and Polar2.
%
% Output:
% self : Polar object that results from the interpolation.
%
% Example:
% MixedPolar = af_tools.interppolar(PolarNACA0012, PolarClarkY, 0.6); Creates an interpolated
% Polar object MixedPolar that results from the interpolation of PolarNACA0012
% and PolarClarkY, assuming that the MixedPolar is made of 60% ClarkY and 40%
% NACA0012.
%
% See also: af_tools.Polar.
%
% -----
% <a href="https://gitlab.uliege.be/am-dept/matlab_airfoil_toolbox">Documentation (online)</a>
% ----------------------------------------------------------------------------------------------
% (c) Copyright 2022-2024 University of Liege
% Author: Thomas Lambert <t.lambert@uliege.be>
% ULiege - Aeroelasticity and Experimental Aerodynamics
% MIT License
% Repo: https://gitlab.uliege.be/am-dept/matlab_airfoil_toolbox
% Issues: https://gitlab.uliege.be/am-dept/matlab_airfoil_toolbox/-/issues
% ----------------------------------------------------------------------------------------------
% --- Defaults and constants
DEBUG = false;
% --- Input checks
% Check if both Polars are Polars objects and weightPolar2 is in [0,1]
validateattributes(Polar1, {'af_tools.Polar'}, {'scalar'}, mfilename(), 'Polar1', 1);
validateattributes(Polar2, {'af_tools.Polar'}, {'scalar'}, mfilename(), 'Polar2', 2);
validateattributes(weightPolar2, {'double'}, {'scalar', 'nonnegative', '<=', 1}, ...
mfilename(), 'weightPolar2', 3);
% Check that both Polars have the Reynolds, Mach and nCrit
equalfieldscheck(Polar1, Polar2, 'reynolds');
equalfieldscheck(Polar1, Polar2, 'mach');
equalfieldscheck(Polar1, Polar2, 'nCrit');
% --- Interpolate the polars
self = af_tools.Polar;
if Polar1 == Polar2
self.properties = self.properties;
self.airfoil = self.airfoil;
self.reynolds = self.reynolds;
self.mach = self.mach;
self.nCrit = self.nCrit;
self.aoa = self.aoa;
self.cl = self. cl;
self.cd = self.cd;
self.cm = self.cm;
self.Ext = self.Ext;
self.Stall = self.Stall;
self.Zero = self.Zero;
self.Lin = self.Lin;
self.extrapMethod = self.extrapMethod;
return
end
wPol2 = weightPolar2;
wPol1 = 1 - wPol2;
nPolars = numel(Polar1.reynolds);
% Set base polar information
afName = {sprintf('MIX-%d%%%s-%d%%%s', round(wPol1 * 100), Polar1.airfoil{1}, ...
round(wPol2 * 100), Polar2.airfoil{2})};
self.airfoil = repmat(afName, 1, nPolars);
self.reynolds = Polar1.reynolds;
self.mach = Polar1.mach;
self.nCrit = Polar1.nCrit;
% Interpolate the angle of attack, cl, cd and cm using the base polars
if numel(Polar1.aoa(:, 1)) ~= numel(Polar2.aoa(:, 1)) || ...
any(Polar1.aoa(:, 1) ~= Polar2.aoa(:, 1))
% Reinterpolate the original polars so they share the same range of AOA.
aoaMin = min(min(Polar1.aoa(:, 1)), min(Polar2.aoa(:, 1)));
aoaMax = max(max(Polar1.aoa(:, 1)), max(Polar2.aoa(:, 1)));
aoaStep = min(diff(Polar1.aoa(1:2, 1)), diff(Polar2.aoa(1:2, 1)));
aoaVect = aoaMin:aoaStep:aoaMax;
aoa1 = repmat(aoaVect', 1, nPolars);
aoa2 = repmat(aoaVect', 1, nPolars);
cl1 = interp1(Polar1.aoa(:, 1), Polar1.cl, aoaVect');
cl2 = interp1(Polar2.aoa(:, 1), Polar2.cl, aoaVect');
cd1 = interp1(Polar1.aoa(:, 1), Polar1.cd, aoaVect');
cd2 = interp1(Polar2.aoa(:, 1), Polar2.cd, aoaVect');
cm1 = interp1(Polar1.aoa(:, 1), Polar1.cm, aoaVect');
cm2 = interp1(Polar2.aoa(:, 1), Polar2.cm, aoaVect');
else
% Use direclty the polar data
aoa1 = Polar1.aoa;
aoa2 = Polar2.aoa;
cl1 = Polar1.cl;
cl2 = Polar2.cl;
cd1 = Polar1.cd;
cd2 = Polar2.cd;
cm1 = Polar1.cm;
cm2 = Polar2.cm;
end
self.aoa = lininterparray(wPol1, aoa1, wPol2, aoa2);
self.cl = lininterparray(wPol1, cl1, wPol2, cl2);
self.cd = lininterparray(wPol1, cd1, wPol2, cd2);
self.cm = lininterparray(wPol1, cm1, wPol2, cm2);
% --- DEBUGGING
if DEBUG
% Plot the original polars for each Reynolds as well as the interpolated one
for j = 1:nPolars
figure('Name', 'Combined Polars');
hold on;
plot(Polar1.aoa(:, j), Polar1.cl(:, j));
plot(Polar2.aoa(:, j), Polar2.cl(:, j));
plot(self.aoa(:, j), self.cl(:, j));
hold off;
legend(Polar1.airfoil{1}, Polar2.airfoil{1}, self.airfoil{1}, 'Location', 'NorthWest');
grid on;
pause;
end
end
end
function interpArray = lininterparray(weight1, array1, weight2, array2)
% LININTERPARRAY Linear interpolation between two arrays
interpArray = weight1 * array1 + weight2 * array2;
end
function equalfieldscheck(Polar1, Polar2, field)
% FIELDCHECK Check if some fields are equal in both Polar objects
if Polar1.(field) ~= Polar2.(field)
error('Polar:interppolar:different%s', ...
['Impossible to interpolate between the two Polar objects as they have '...
'different %s.\n'...
'Please use two Polar objects with identical AOA, Reynolds, Mach and nCrit'], ...
field, field);
end
end
+af_tools/findstall.m
View file @ 0995fa2a
function [alpha_s, c_l_s, c_d_s] = findstall(varargin)
% FINDSTALL Finds the stall point and returns the corresponding AOA, CL and CD.
% This function uses Matlab's findpeak in order to determine the position of the stall point.
% This function looks for the highest local maximm in the CL curve to determine the position
% of the stall point.
%
% FINDSTALL accepts inputs in the form of a Polar structure (generated with XF2MAT) or values
% for angles of attack and their associated cl and cd. While alpha and cl are required for
......@@ -99,22 +100,11 @@ function [alpha_s, c_l_s, c_d_s] = stallfinder(alpha, c_l, c_d)
% Loop for every polar data (every column of cl)
for i = 1:size(c_l, 2)
% Find stall for each polar using findpeaks
maxFound = false;
minWidth = 0;
% Iterate in order to find only one peak (largest one should indicate the stall)
while ~maxFound && minWidth <= length(alpha(:, i))
[~, locs] = findpeaks(c_l(:, i), 'MinPeakWidth', minWidth);
if numel(locs) == 1
maxFound = true;
else
minWidth = minWidth + 1;
end
end
% Find all local maxima in c_l, then assume the highest one is the stall point
locMax = islocalmax(c_l(:, i));
[~, absMaxIdx] = max(locMax);
if ~maxFound
% If findpeaks does not work
if ~any(locMax)
if PLOT_POLAR_IF_ERROR
figure('Name', 'Debug: Polar with no stall');
plot(rad2deg(alpha(:, i)), c_l(:, i));
......@@ -124,15 +114,15 @@ function [alpha_s, c_l_s, c_d_s] = stallfinder(alpha, c_l, c_d)
grid on;
end
error('MATLAB:findstall:noStallFound', ...
['findpeaks did not found any stall point the polar. '...
['The stall point could not be determined as no local maximum was detected. '...
'Please provide a polar that goes a bit after the stall.\n' ...
'See the attached plot for more details']);
end
alpha_s(i) = alpha(locs, i);
c_l_s(i) = c_l(locs, i);
alpha_s(i) = alpha(absMaxIdx, i);
c_l_s(i) = c_l(absMaxIdx, i);
if ~isempty(c_d)
c_d_s(i) = c_d(locs, i);
c_d_s(i) = c_d(absMaxIdx, i);
end
end
......
+af_tools/nacacamber.m
View file @ 0995fa2a
......@@ -6,13 +6,13 @@ function [coord, gradY] = nacacamber(varargin)
% -----
%
% Usage:
% [X, Y] = NACACAMBER prompts the user for the airfoil digits, then calculates the coordinates
% of the camberline using 100 points spaced with the half-cosine rule.
% [coord, gradY] = NACACAMBER prompts the user for the airfoil digits, then calculates the
% coordinates of the camberline using 100 points spaced with the half-cosine rule.
%
% [X, Y] = NACACAMBER(DIGITS) calculates the coordinates of the camberline for the NACA
% airfoil specified by DIGITS, without further input.
% [coord, gradY] = NACACAMBER(DIGITS) calculates the coordinates of the camberline for the
% NACA airfoil specified by DIGITS, without further input.
%
% [X, Y] = NACACAMBER(DIGITS, NPOINTS) calculates the coordinates of the NACA airfoil
% [coord, gradY] = NACACAMBER(DIGITS, NPOINTS) calculates the coordinates of the NACA airfoil
% specified by DIGITS, using NPOINTS points spaced with the half-cosine rule.
%
% [...] = NACACAMBER(..., 'spacing', SP) uses the defined spacing rule to determine to
......@@ -20,14 +20,14 @@ function [coord, gradY] = nacacamber(varargin)
% 'linear' - linear
% 'halfcosine' - (default) finer at leading and trailing edges
% 'cosine' - finer at leading edge, coarser at the trailing edge
% 'sine' - finer at trailing edge, coarser at the leading edge
%
% Inputs:
% digits : Character vector representing the NACA airfoil (ex: '0012', '24120', '24012')
% nPoints : Number of output coordinates
%
% Output:
% xc : Camberline X coordinates (0 <= X <= 1)
% yc : Camberline Y coordinates
% coord : Camberline coordinates [nPts x 2]
% gradY : Gradient of Y (required to get the whole airfoil coordinates)
%
% Example:
......@@ -41,7 +41,7 @@ function [coord, gradY] = nacacamber(varargin)
% <a href="https://gitlab.uliege.be/am-dept/matlab_airfoil_toolbox">Documentation (online)</a>
% -------------------------------------------
% (c) Copyright 2022 University of Liege
% (c) Copyright 2022-2023 University of Liege
% Author: Thomas Lambert <t.lambert@uliege.be>
% ULiege - Aeroelasticity and Experimental Aerodynamics
% MIT License
......
CHANGELOG.md
View file @ 0995fa2a
......@@ -17,6 +17,43 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
### Fixed
## [4.4.2] - 2024-03-13
### Fixed
- Wrong coordinates array for interpolated Airfoil
## [4.4.1] - 2024-03-13
### Fixed
- Issue with interpolation of polars between polars with different range of AOA
## [4.4.0] - 2024-03-13
### Added
- Creation of interpolated Polar between two reference Polars
- Creation of interpolated Airfoil between two reference Airfoils
## [4.3.0] - 2023-12-14
### Changed
- **findstall**: use `islocalmax` rather than `findpeaks` to find stall point.
## [4.2.1] - 2023-02-17
### Fixed
- **filesinput**: issue when single fileinput
## [4.2.0] - 2023-02-13
### Added
- **Polar**: add `extrapMethod` to `getcoeff`
## [4.1.0] - 2022-11-25
### Added
......@@ -127,7 +164,13 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
- Initial release
[Unreleased]: https://gitlab.uliege.be/am-dept/matlab_airfoil_toolbox/-/compare/v4.1.0...master
[Unreleased]: https://gitlab.uliege.be/am-dept/matlab_airfoil_toolbox/-/compare/v4.4.2...master
[4.4.2]: https://gitlab.uliege.be/am-dept/matlab_airfoil_toolbox/-/compare/4.4.1...v4.4.2
[4.4.1]: https://gitlab.uliege.be/am-dept/matlab_airfoil_toolbox/-/compare/4.4.0...v4.4.1
[4.4.0]: https://gitlab.uliege.be/am-dept/matlab_airfoil_toolbox/-/compare/4.3.0...v4.4.0
[4.3.0]: https://gitlab.uliege.be/am-dept/matlab_airfoil_toolbox/-/compare/4.2.1...v4.3.0
[4.2.1]: https://gitlab.uliege.be/am-dept/matlab_airfoil_toolbox/-/compare/4.2.0...v4.2.1
[4.2.0]: https://gitlab.uliege.be/am-dept/matlab_airfoil_toolbox/-/compare/4.1.0...v4.2.0
[4.1.0]: https://gitlab.uliege.be/am-dept/matlab_airfoil_toolbox/-/compare/4.0.0...v4.1.0
[4.0.0]: https://gitlab.uliege.be/am-dept/matlab_airfoil_toolbox/-/compare/3.0.0...4.0.0
[3.0.0]: https://gitlab.uliege.be/am-dept/matlab_airfoil_toolbox/-/compare/v2.0.1...v3.0.0
......