Skip to content
Snippets Groups Projects
Verified Commit 11501f87 authored by Thomas Lambert's avatar Thomas Lambert :helicopter:
Browse files

feat(plots): add plots for the presentation

parent a336d6f2
No related branches found
Tags v0.0.0
No related merge requests found
Hide whitespace changes
Inline Side-by-side
utils/analysis/plotcf.m
+ 0
− 2
View file @ 11501f87
......@@ -26,7 +26,6 @@ function plotcf(ResData)
0.03, 3.5, 4.6
0.03, 3, 7.7
0.10, 3, 7.7
0.10, 3, 4.6
% 0.03, 2, 7.7
% 0.03, 2, 7.7
......@@ -326,7 +325,6 @@ function setgca(varargin)
% xlim([-180, 180]);
xlim([0, 360]);
ylim([-0.3, 0.6]);
xlabel('Phase offset [deg]');
ylabel('\cl and \cd [-]');
......
utils/analysis/plotcfpres.m 0 → 100644
+ 176
− 0
View file @ 11501f87
function plotcfpres(ResData)
% PLOTCF Plot force coefficients
% Todo
% ----------------------------------------------------------------------------------------------
% (c) Copyright 2022 University of Liege
% Author: Thomas Lambert <t.lambert@uliege.be>
% ULiege - Aeroelasticity and Experimental Aerodynamics
% MIT License
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% X-axis = flapping frequency
% Y-axis = Cf
% One figure per dx
% Plot airspeeds in different colors
set(0, 'defaulttextinterpreter', 'latex');
ONLY_ERR_BARS = true;
ERR_CAPSIZE = 3;
OFFSET_STEPS = 20;
CHORD = 0.05;
INTERESTING_PLOTS = [
0.03, 3, 4.6
0.03, 3.5, 4.6
0.10, 3, 4.6
];
SAVE_TIKZ_PARAM = []; % Points to save ([0,0,0] to make sur nothing is save)
for j = 1:size(ResData, 2)
for i = 1:size(ResData, 1)
if isinteresting(ResData(i, j), INTERESTING_PLOTS) && ...
~(ResData(i, j).Testcase.dx == 0.05 && ResData(i, j).Testcase.freq <= 2) && ...
~(ResData(i, j).Testcase.dx == 0.03 && ResData(i, j).Testcase.freq == 1.25)
trueF = ResData(i, j).Front.trueFreq;
redFreq = redfreq(trueF, CHORD, ResData(i, j).Testcase.airspeed);
expParam = sprintf('d = %0.2f, f = %0.2f, V = %0.1f (k = %0.02f)', ...
ResData(i, j).Testcase.dx, ...
ResData(i, j).Front.trueFreq, ...
ResData(i, j).Testcase.airspeed, ...
redFreq);
offsets = ResData(i, j).AllPhases.offsets;
offsets = wrapTo360(offsets); % Wrap offsets in [0,360];
cleanOffsets = OFFSET_STEPS / 2:OFFSET_STEPS:360 - OFFSET_STEPS / 2;
meanOffset = nan(size(cleanOffsets));
stdOffset = nan(size(cleanOffsets));
for iOff = 1:length(cleanOffsets)
curOffsets = cleanOffsets(iOff);
offIdx = find(abs(offsets - curOffsets) < diff(cleanOffsets(1:2)) / 2);
meanOffset(iOff) = mean(offsets(offIdx));
stdOffset(iOff) = std(offsets(offIdx));
meanF(iOff, :) = mean(ResData(i, j).AllPhases.Front.cF(offIdx, :), 1);
stdF(iOff, :) = std(ResData(i, j).AllPhases.Front.cF(offIdx, :), 0, 1);
meanA(iOff, :) = mean(ResData(i, j).AllPhases.Aft.cF(offIdx, :), 1);
stdA(iOff, :) = std(ResData(i, j).AllPhases.Aft.cF(offIdx, :), 0, 1);
end
ynegF = -stdF;
yposF = stdF;
ynegA = -stdA;
yposA = stdA;
xneg = cleanOffsets - meanOffset - stdOffset;
xpos = cleanOffsets - meanOffset + stdOffset;
figure('Name', 'myfig');
ax1 = subplot(211);
setcolormap();
hold on;
if ~ONLY_ERR_BARS
plot(ResData(i, j).AllPhases.offsets, ...
ResData(i, j).AllPhases.Front.cF(:, [1, 3]), 'o');
end
e1 = errorbar(cleanOffsets, meanF(:, 3), yposF(:, 3), 's', 'Capsize', ERR_CAPSIZE);
e2 = errorbar(cleanOffsets, meanF(:, 1), yposF(:, 1), 'o', 'Capsize', ERR_CAPSIZE);
e1.MarkerFaceColor = e1.Color;
e2.MarkerFaceColor = e2.Color;
hold off;
setgca('XTickLabel', [], ...
'XLabel', []);
% legend({'\cl', '\cd'}, 'Location', 'SouthWest');
% legend('boxoff');
% ax1.Position(3) = 1.25 * ax1.Position(3);
ax2 = subplot(212);
setcolormap();
hold on;
if ~ONLY_ERR_BARS
plot(ResData(i, j).AllPhases.offsets, ...
ResData(i, j).AllPhases.Aft.cF(:, [1, 3]), 'o');
end
e3 = errorbar(cleanOffsets, meanA(:, 3), yposA(:, 3), 's', 'Capsize', ERR_CAPSIZE);
e4 = errorbar(cleanOffsets, meanA(:, 1), yposA(:, 1), 'o', 'Capsize', ERR_CAPSIZE);
e3.MarkerFaceColor = e3.Color;
e4.MarkerFaceColor = e4.Color;
hold off;
setgca();
linkaxes([ax1 ax2], 'xy');
ax2.Position(3) = ax1.Position(3);
% Add annotations
% Increase axes Y lim by 10% of the data range before writing text
FACTOR = 0.1;
yDiff = diff(ax1.YLim);
ax1.YLim = [ax1.YLim(1), ax1.YLim(2) + FACTOR * yDiff];
if isinteresting(ResData(i, j), SAVE_TIKZ_PARAM)
pause(1);
figdir = 'figures/results/';
filename = sprintf('flapWindOn-Dx%0.2f_F%0.2f_V%0.1f_PRES.tex', ...
ResData(i, j).Testcase.dx, ...
ResData(i, j).Testcase.freq, ...
ResData(i, j).Testcase.airspeed);
save2tikz([figdir, filename], '\small');
end
% Update titles after having saved
title(ax1, ['Front - ', expParam]);
title(ax2, ['Aft - ', expParam]);
end
end
end
end
function setgca(varargin)
% xlim([-180, 180]);
xlim([0, 360]);
ylim([-0.3, 0.6]);
xlabel('$\Delta\phi$ [deg]');
% ylabel('\cl and \cd [-]');
% xticks([-180:90:180]);
xticks([0:90:360]);
yticks([-0.3:0.3:0.6]);
set(gca, ...
'Box', 'off', ...
'TickDir', 'out', ...
'TickLength', [.02 .02], ...
'XMinorTick', 'off', ...
'YMinorTick', 'off', ...
'YGrid', 'on', ...
'XGrid', 'on', ...
'XColor', 'k', ...
'YColor', 'k', ...
'GridLineStyle', ':', ...
'GridColor', 'k', ...
'GridAlpha', 0.25, ...
'LineWidth', 1, ...
'FontName', 'Helvetica', ...
'Fontsize', 14);
if nargin > 0
set(gca, varargin{:});
end
end
function addtxt(ax, txt, ypos)
text(ax, 10, ypos, ...
txt, ...
'HorizontalAlignment', 'left', 'VerticalAlignment', 'top', 'FontWeight', 'bold');
end
utils/plotsemiflappres.m 0 → 100644
+ 378
− 0
View file @ 11501f87
function ResData = plotsemiflap(type)
% PLOTWINDONON Plots the forces in wind on and wind off conditions
% Todo
% ----------------------------------------------------------------------------------------------
% (c) Copyright 2022 University of Liege
% Author: Thomas Lambert <t.lambert@uliege.be>
% ULiege - Aeroelasticity and Experimental Aerodynamics
% MIT License
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
set(0, 'defaulttextinterpreter', 'latex');
POSITION = 'Front';
INTERESTING_POINTS = [
];
PLOT_TIME = true;
SAVE_TIKZ_PARAM = []; % INTERESTING_POINTS(1:2, :);
MYGREY = [186, 216, 221] / 255;
ALLPLOTS = false;
SAMPLING = 1000;
ERR_CAPSIZE = 3;
CHORD = 0.05;
SPAN = 2 * 0.2;
C_TO_K = 273.15;
switch type
case 'Front' % Front wing is flapping
dat3cm = '2022-3cm-flapFront.mat';
dat5cm = '2022-5cm-flapFront.mat';
notType = 'Aft';
case 'Aft' % Aft wing is flapping
dat3cm = '2022-3cm-flapAft.mat';
dat5cm = '2022-5cm-flapAft.mat';
notType = 'Front';
end
Tmp = load(['data/cleaned_MAT/2022/', dat3cm]);
ResData(1, :) = Tmp.ExpData;
Tmp = load(['data/cleaned_MAT/2022/', dat5cm]);
ResData(2, :) = Tmp.ExpData;
% Remove useless data
ResData = rmfield(ResData, {'arduStart', 'tunnelStart'}); % Useless
% FIXME: remove deprecated
ResData = rmfield(ResData, {'arduFile', 'tunnelFile', ...
'dx', 'airspeed', 'freq'});
for i = 1:size(ResData, 1)
% for i = 1
for j = 1:size(ResData, 2)
dens = airdensity(ResData(i, j).Testcase.temp + C_TO_K, ...
ResData(i, j).Testcase.press);
speed = ResData(i, j).Testcase.airspeed;
% Calculate main frequency
ResData(i, j).(type).trueFreq = getmainfrequency(ResData(i, j).(type).angles, ...
1 / diff(ResData(i, j).time(1:2)));
ResData(i, j).(notType).trueFreq = 0;
% Calculate inertia
woIdx = (j) - mod(j - 1, 4); % Corresponding wind-off index
[ResData(i, j).(type).inertiaF, ...
ResData(i, j).(type).inertiaM] = calcinertia(ResData(i, woIdx), type);
ResData(i, j).(notType).inertiaF = zeros(size(ResData(i, j).(type).inertiaF));
ResData(i, j).(notType).inertiaM = zeros(size(ResData(i, j).(type).inertiaM));
% Average everything over a mean period for the active wing
[avSignal, peakIdx] = averageperiods(ResData(i, j).(type).angles, ...
SAMPLING, ...
ResData(i, j).(type).trueFreq);
ResData(i, j).(type).avAngles = avSignal;
[ResData(i, j).(type).avFz, ResData(i, j).(type).stdFz] = getavsignal(ResData(i, j).(type).forces(:, 3), peakIdx); % mh:ignore_style (WONTFIX)
[ResData(i, j).(type).avFx, ResData(i, j).(type).stdFx] = getavsignal(ResData(i, j).(type).forces(:, 1), peakIdx); % mh:ignore_style (WONTFIX)
[ResData(i, j).(notType).avFz, ResData(i, j).(notType).stdFz] = getavsignal(ResData(i, j).(notType).forces(:, 3), peakIdx); % mh:ignore_style (WONTFIX)
[ResData(i, j).(notType).avFx, ResData(i, j).(notType).stdFx] = getavsignal(ResData(i, j).(notType).forces(:, 1), peakIdx); % mh:ignore_style (WONTFIX)
% Remove inertia and invert FX so it is oriented alongside drag diraction
ResData(i, j).(type).avFz = ResData(i, j).(type).avFz - ...
ResData(i, j).(type).inertiaF(3);
ResData(i, j).(type).avFx = -(ResData(i, j).(type).avFx - ...
ResData(i, j).(type).inertiaF(1));
% Calculate coefficient
for iPos = {'Aft', 'Front'}
pos = string(iPos);
ResData(i, j).(pos).cFz = forcetocoeff(ResData(i, j).(pos).avFz, ...
dens, speed, CHORD, SPAN, 'force');
ResData(i, j).(pos).cFx = forcetocoeff(ResData(i, j).(pos).avFx, ...
dens, speed, CHORD, SPAN, 'force');
ResData(i, j).(pos).stdFz = forcetocoeff(ResData(i, j).(pos).stdFz, ...
dens, speed, CHORD, SPAN, 'force');
ResData(i, j).(pos).stdFx = forcetocoeff(ResData(i, j).(pos).stdFx, ...
dens, speed, CHORD, SPAN, 'force');
% Calculate mean CL and mean CD over the average period
ResData(i, j).(pos).cFz_mean = mean(ResData(i, j).(pos).cFz);
ResData(i, j).(pos).cFx_mean = mean(ResData(i, j).(pos).cFx);
% Get mean STD over the period
ResData(i, j).(pos).cFz_std = mean(ResData(i, j).(pos).stdFz);
ResData(i, j).(pos).cFx_std = mean(ResData(i, j).(pos).stdFx);
end
% -----------------------------------------------
% Plot each test case
if ALLPLOTS
redFreq = redfreq(ResData(i, j).(type).trueFreq, ...
CHORD, ResData(i, j).Testcase.airspeed);
expParam = sprintf(['Aft angle = %0.1f - d = %0.2f, f = %0.2f, '...
'V = %0.1f (k = %0.02f)'], ...
ResData(i, j).(notType).angles(1), ...
ResData(i, j).Testcase.dx, ...
ResData(i, j).(type).trueFreq, ...
ResData(i, j).Testcase.airspeed, ...
redFreq);
[~, donwStrokeEndIdx] = min(ResData(i, j).(type).avAngles);
nPeriodIdx = 1 / ResData(i, j).(type).trueFreq * SAMPLING;
time = (0:length(ResData(i, j).(type).avAngles) - 1) / nPeriodIdx;
patchX = [0, donwStrokeEndIdx / nPeriodIdx, donwStrokeEndIdx / nPeriodIdx, 0];
patchY = [-5, -5, 10, 10];
figure('Name', 'HalfFlap');
ax1 = subplot(211);
setcolormap();
hold on;
patch(patchX, patchY, [186, 216, 221] / 255, 'EdgeColor', 'none', 'FaceAlpha', 0.5);
h1 = plot(time, ResData(i, j).Front.avFz);
h2 = plot(time, ResData(i, j).Front.avFx);
hold off;
setgca('XTickLabel', [], ...
'XLabel', []);
legend([h1, h2], {'\fz', '\fx'}, 'Location', 'NorthEast');
legend('boxoff');
% ax1.Position(3) = 1.25 * ax1.Position(3);
ax2 = subplot(212);
setcolormap();
hold on;
patch(patchX, patchY, [186, 216, 221] / 255, 'EdgeColor', 'none', 'FaceAlpha', 0.5);
h1 = plot(time, ResData(i, j).Aft.avFz);
h2 = plot(time, ResData(i, j).Aft.avFx);
hold off;
setgca();
linkaxes([ax1 ax2], 'xy');
ax2.Position(3) = ax1.Position(3);
% Increase axes Y lim by 10% of the data range before writing text
FACTOR = 0.1;
yDiff = diff(ax1.YLim);
ax1.YLim = [ax1.YLim(1), ax1.YLim(2) + FACTOR * yDiff];
addtxt(ax1, '\bf Front', ax1.YLim(2));
addtxt(ax2, '\bf Aft', ax2.YLim(2));
text(ax2, (donwStrokeEndIdx / nPeriodIdx) / 2, 0.5, ...
'Downstroke', ...
'HorizontalAlignment', 'center');
text(ax2, donwStrokeEndIdx / nPeriodIdx + ...
(1 - (donwStrokeEndIdx / nPeriodIdx)) / 2, 0.5, ...
'Upstroke', ...
'HorizontalAlignment', 'center');
%
% if isinteresting(ResData(i, j), SAVE_TIKZ_PARAM)
% pause(1);
% figdir = 'figures/results/';
% filename = sprintf('flap%s-Dx%0.2f_F%0.2f_V%0.1f.tex', ...
% type, ...
% ResData(i, j).Testcase.dx, ...
% ResData(i, j).Testcase.freq, ...
% ResData(i, j).Testcase.airspeed);
% save2tikz([figdir, filename], '\small');
% end
% Update titles after having saved
title(ax1, ['Front - ', expParam]);
title(ax2, ['Aft - ', expParam]);
end
end
end
% Loop on all cases and create useful vectors
nSpeeds = 4; % Number of airspeeds tested
for iDx = 1:size(ResData, 1)
for j = 2:nSpeeds
speedIdx = j:4:size(ResData, 2);
for k = 1:length(speedIdx)
speeds(j, k) = ResData(iDx, speedIdx(k)).Testcase.airspeed;
freq(j, k) = ResData(iDx, speedIdx(k)).Testcase.freq;
dx(j, k) = ResData(iDx, speedIdx(k)).Testcase.dx;
ang2(j, k) = mean(ResData(iDx, speedIdx(k)).Aft.angles);
cfzF(j, k) = ResData(iDx, speedIdx(k)).Front.cFz_mean;
cfxF(j, k) = ResData(iDx, speedIdx(k)).Front.cFx_mean;
cfzA(j, k) = ResData(iDx, speedIdx(k)).Aft.cFz_mean;
cfxA(j, k) = ResData(iDx, speedIdx(k)).Aft.cFx_mean;
stdcfzF(j, k) = ResData(iDx, speedIdx(k)).Front.cFz_std;
stdcfxF(j, k) = ResData(iDx, speedIdx(k)).Front.cFx_std;
stdcfzA(j, k) = ResData(iDx, speedIdx(k)).Aft.cFz_std;
stdcfxA(j, k) = ResData(iDx, speedIdx(k)).Aft.cFx_std;
end
figure('Name', 'Dx');
ax1 = subplot(211);
setcolormap();
hold on;
e11 = errorbar(ang2(j, 1:3), cfzF(j, 4:end), ...
stdcfzF(j, 1:3), '-s', 'Capsize', ERR_CAPSIZE);
e12 = errorbar(ang2(j, 1:3), cfxF(j, 4:end), ...
stdcfxF(j, 1:3), '-o', 'Capsize', ERR_CAPSIZE);
e11.MarkerFaceColor = e11.Color;
e12.MarkerFaceColor = e12.Color;
hold off;
setgca2('XTickLabel', [], ...
'XLabel', []);
ax2 = subplot(212);
setcolormap();
hold on;
e21 = errorbar(ang2(j, 4:end), cfzA(j, 4:end), ...
stdcfzF(j, 4:end), '-s', 'Capsize', ERR_CAPSIZE);
e22 = errorbar(ang2(j, 4:end), cfxA(j, 4:end), ...
stdcfxF(j, 4:end), '-o', 'Capsize', ERR_CAPSIZE);
e21.MarkerFaceColor = e21.Color;
e22.MarkerFaceColor = e22.Color;
hold off;
setgca2();
linkaxes([ax1 ax2], 'xy');
% Add annotations
% Increase axes Y lim by 10% of the data range before writing text
FACTOR = 0.1;
yDiff = diff(ax1.YLim);
ax1.YLim = [ax1.YLim(1), ax1.YLim(2) + FACTOR * yDiff];
pause(1);
figdir = 'figures/results/';
filename = sprintf('half%s-AllF-Dx%0.2f_V%0.1f_PRES.tex', ...
type, ...
dx(j, 1), ...
speeds(j, 1));
save2tikz([figdir, filename], '\small');
expParam = sprintf('Dx = %0.02f, U = %0.2f', ...
dx(j, 1), ...
speeds(j, 1));
title(ax1, expParam);
end
end
end
function addtxt(ax, txt, ypos)
text(ax, -10, ypos, ...
txt, ...
'HorizontalAlignment', 'left', 'VerticalAlignment', 'top', 'FontWeight', 'bold');
end
function [inertiaF, inertiaM] = calcinertia(ResDataWO, type)
% CALCINERTIA Calcualte inertia forces
SAMPLING = 1000;
% Get first and last peaks for wind-off conditions
woBounds = getboundingpeaks(ResDataWO.(type).angles, SAMPLING);
woIdx = woBounds(1):woBounds(2);
% Calculate mean force and moments due to inertia
inertiaF = mean(ResDataWO.(type).forces(woIdx, :));
inertiaM = mean(ResDataWO.(type).moments(woIdx, :));
end
function [avSignal, peakIdx] = averageperiods(signal, sampling, freq)
% AVERAGEPERIOD Average all periods of the signal and return the average signal
% Find all peaks location
period = 1 / freq;
nPeriodIdx = period * sampling;
peakLocs = getallpeaks(signal, sampling);
% Remove first and last peaks
peakLocs = peakLocs(2:end - 1);
% Split signal in a matrix whith rows for the peak number and cols to the INDEXES of values in
for iPeak = length(peakLocs):-1:1
peakIdx(iPeak, :) = peakLocs(iPeak):peakLocs(iPeak) + round(nPeriodIdx);
peakSignal(iPeak, :) = signal(peakIdx(iPeak, :));
end
% Take the mean value of the signal formed with all these peaks
avSignal = mean(peakSignal);
end
function [avSignal, stdSignal] = getavsignal(signal, peakIdx)
for iPeak = size(peakIdx, 1):-1:1
peakSignal(iPeak, :) = signal(peakIdx(iPeak, :));
end
% Take the mean value of the signal formed with all these peaks
avSignal = mean(peakSignal);
stdSignal = std(peakSignal);
end
function setgca(varargin)
ylim([-1.5, 1.5]);
xlim([0, 1]);
xlabel('t/T [-]');
ylabel('\fz and \fx [N]');
set(gca, ...
'Box', 'off', ...
'TickDir', 'out', ...
'TickLength', [.02 .02], ...
'XMinorTick', 'off', ...
'YMinorTick', 'off', ...
'YGrid', 'off', ...
'XGrid', 'off', ...
'XColor', 'k', ...
'YColor', 'k', ...
'GridLineStyle', ':', ...
'GridColor', 'k', ...
'GridAlpha', 0.25, ...
'LineWidth', 1, ...
'FontName', 'Helvetica', ...
'Fontsize', 14);
if nargin > 0
set(gca, varargin{:});
end
end
function setgca2(varargin)
xlim([-13, 40]);
xticks([-10:10:40]);
xlabel('$\Lambda_\text{aft}$ [deg]');
set(gca, ...
'Box', 'off', ...
'TickDir', 'out', ...
'TickLength', [.02 .02], ...
'XMinorTick', 'off', ...
'YMinorTick', 'off', ...
'YGrid', 'off', ...
'XGrid', 'off', ...
'XColor', 'k', ...
'YColor', 'k', ...
'GridLineStyle', ':', ...
'GridColor', 'k', ...
'GridAlpha', 0.25, ...
'LineWidth', 1, ...
'FontName', 'Helvetica', ...
'Fontsize', 14);
if nargin > 0
set(gca, varargin{:});
end
end
utils/plotwindonoff.m
+ 3
− 0
View file @ 11501f87
......@@ -303,6 +303,9 @@ end
function setgcpres(varargin)
xlabel('t');
ylabel('\fz');
set(gca, ...
'Box', 'off', ...
'TickDir', 'out', ...
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment