From 55dff5cacede17b4f11e8272622e7687b7b09733 Mon Sep 17 00:00:00 2001 From: Thomas Lambert <t.lambert@uliege.be> Date: Mon, 12 Dec 2022 15:00:40 +0100 Subject: [PATCH] feat(plots): add static and semi flap plots --- utils/analysis/plotstatictests.m | 128 ++++++++++ utils/plotsemiflap.m | 418 +++++++++++++++++++++++++++++++ 2 files changed, 546 insertions(+) create mode 100644 utils/analysis/plotstatictests.m create mode 100644 utils/plotsemiflap.m diff --git a/utils/analysis/plotstatictests.m b/utils/analysis/plotstatictests.m new file mode 100644 index 0000000..12c5142 --- /dev/null +++ b/utils/analysis/plotstatictests.m @@ -0,0 +1,128 @@ +function ParsedDat = plotstatictests(StatData) + % PLOTSTATICTESTS Plot CL and CD for each set of wing w.r.t front and aft dihedral + % Todo + + % ---------------------------------------------------------------------------------------------- + % (c) Copyright 2022 University of Liege + % Author: Thomas Lambert <t.lambert@uliege.be> + % ULiege - Aeroelasticity and Experimental Aerodynamics + % MIT License + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + close all; + + % Constants and defaults + CHORD = 0.05; + SPAN = 2 * 0.2; + C_TO_K = 273.15; + AIRSPEEDS = [5, 8, 10, 14]; + + % Create useful vectors for each sets of front/aft position + % Vectors dimensions = number of speeds + + % First we get all unique angles for front and aft wings that were tested + for i = numel(StatData):-1:1 + angleF(i) = StatData(i).Front.angles(1); + angleA(i) = StatData(i).Aft.angles(1); + end + angleF = unique(angleF); + angleA = unique(angleA); + nVeloc = numel(StatData) / (length(angleF) * length(angleA)); % Number of velocities per config + + % Create base mesh grid for all tests cases + [gridF, gridA] = meshgrid(angleF, angleA); + + % Initialize force coefficients in X and Y for Front and Aft + cfzF = zeros([size(gridF), nVeloc]); + cfxF = zeros([size(gridF), nVeloc]); + cfzA = zeros([size(gridF), nVeloc]); + cfxA = zeros([size(gridF), nVeloc]); + stdfz = zeros([size(gridF), nVeloc]); + + for i = 1:numel(StatData) + + idx1 = gridF == StatData(i).Front.angles(1) & gridA == StatData(i).Aft.angles(1); + [~, idx2] = min(abs(AIRSPEEDS - StatData(i).Testcase.airspeed)); + + for iPos = {'Front', 'Aft'} + pos = string(iPos); + + fz = mean(StatData(i).(pos).forces(:, 3)); + fx = mean(StatData(i).(pos).forces(:, 1)); + % stdfz = std(StatData(i).(pos).forces(:, 3)); + % stdfx = std(StatData(i).(pos).forces(:, 1)); + + dens = airdensity(StatData(i).Testcase.temp + C_TO_K, ... + StatData(i).Testcase.press); + speed = StatData(i).Testcase.airspeed; + + cfz = forcetocoeff(fz, dens, speed, CHORD, SPAN, 'force'); + cfx = forcetocoeff(fx, dens, speed, CHORD, SPAN, 'force'); + + % Nondim with projected area + cfz = cfz / sind(90 + StatData(i).(pos).angles(1)); + cfx = cfx / sind(90 + StatData(i).(pos).angles(1)); + + if strcmp(pos, 'Front') + cfzF(:, :, idx2) = cfzF(:, :, idx2) + idx1 * fz; + cfxF(:, :, idx2) = cfxF(:, :, idx2) + idx1 * fx; + cfzF(:, :, idx2) = cfzF(:, :, idx2) + idx1 * fz; + cfxF(:, :, idx2) = cfxF(:, :, idx2) + idx1 * fx; + else + cfzA(:, :, idx2) = cfzA(:, :, idx2) + idx1 * fz; + cfxA(:, :, idx2) = cfxA(:, :, idx2) + idx1 * fx; + end + + end + + end + + ParsedDat.cfzF = cfzF; + ParsedDat.cfxF = cfxF; + ParsedDat.cfzA = cfzA; + ParsedDat.cfxA = cfxA; + + for i = 1:nVeloc + figure('Name', ['V=', i]); + ax1 = subplot(221); % Front CL + surf(gridF, gridA, cfzF(:, :, i)); + setgca(); + + ax2 = subplot(222); % Front CD + surf(gridF, gridA, cfxF(:, :, i)); + setgca(); + + ax3 = subplot(223); % Aft CL + surf(gridF, gridA, cfzA(:, :, i)); + setgca(); + + ax4 = subplot(224); % Aft CD + surf(gridF, gridA, cfxA(:, :, i)); + setgca(); + + % linkprop([ax1 ax3], 'Zlim'); % CL plots + % linkprop([ax2 ax4], 'Zlim'); % CD plots + end + + for i = 1:length(angleF) + figure; + size(angleA); + size(cfzF(i, :, 4)); + hold on; + title([angleF(i)]); + errobar(angleA, cfzF(:, i, 3), stdfzF(:, i, 3), '-o'); + plot(angleA, cfzF(:, i, 3), '-o'); + plot(angleA, cfzA(:, i, 3), '-o'); + plot(angleA, cfzF(:, i, 3) + cfzA(:, i, 3), '--x'); + hold on; + end + +end + +function setgca() + % SETGCA Sets the axes parameters + + xlabel('Front'); + ylabel('Aft'); + grid on; +end diff --git a/utils/plotsemiflap.m b/utils/plotsemiflap.m new file mode 100644 index 0000000..499862d --- /dev/null +++ b/utils/plotsemiflap.m @@ -0,0 +1,418 @@ +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(221); + setcolormap(); + hold on; + e11 = errorbar(ang2(j, 1:3), cfzF(j, 1:3), ... + stdcfzF(j, 1:3), '-s', 'Capsize', ERR_CAPSIZE); + e12 = errorbar(ang2(j, 1:3), cfxF(j, 1:3), ... + stdcfxF(j, 1:3), '-o', 'Capsize', ERR_CAPSIZE); + e11.MarkerFaceColor = e11.Color; + e12.MarkerFaceColor = e12.Color; + hold off; + setgca2('XTickLabel', [], ... + 'XLabel', []); + + legend({'\cl', '\cd'}, 'Location', 'SouthWest'); + legend('boxoff'); + + ax2 = subplot(222); + setcolormap(); + hold on; + e21 = errorbar(ang2(j, 4:end), cfzF(j, 4:end), ... + stdcfzF(j, 4:end), '-s', 'Capsize', ERR_CAPSIZE); + e22 = errorbar(ang2(j, 4:end), cfxF(j, 4:end), ... + stdcfxF(j, 4:end), '-o', 'Capsize', ERR_CAPSIZE); + e21.MarkerFaceColor = e21.Color; + e22.MarkerFaceColor = e22.Color; + hold off; + setgca2('XTickLabel', [], ... + 'XLabel', [], ... + 'YTickLabel', [], ... + 'YLabel', []); + ax3 = subplot(223); + setcolormap(); + hold on; + cosd(ang2(j, 4:end)); + e31 = errorbar(ang2(j, 1:3), cfzA(j, 1:3) ./ cosd(ang2(j, 4:end)), ... + stdcfzA(j, 1:3), '-s', 'Capsize', ERR_CAPSIZE); + e32 = errorbar(ang2(j, 1:3), cfxA(j, 1:3) ./ cosd(ang2(j, 4:end)), ... + stdcfxA(j, 1:3), '-o', 'Capsize', ERR_CAPSIZE); + e31.MarkerFaceColor = e31.Color; + e32.MarkerFaceColor = e32.Color; + hold off; + setgca2(); + ylabel('$\cl^*$ and $\cd^*$ [-]'); + + ax4 = subplot(224); + setcolormap(); + hold on; + e41 = errorbar(ang2(j, 4:end), cfzA(j, 4:end) ./ cosd(ang2(j, 4:end)), ... + stdcfzA(j, 4:end), '-s', 'Capsize', ERR_CAPSIZE); + e42 = errorbar(ang2(j, 4:end), cfxA(j, 4:end) ./ cosd(ang2(j, 4:end)), ... + stdcfxA(j, 4:end), '-o', 'Capsize', ERR_CAPSIZE); + e41.MarkerFaceColor = e41.Color; + e42.MarkerFaceColor = e42.Color; + hold off; + setgca2('YTickLabel', [], ... + 'YLabel', []); + + linkaxes([ax1 ax2 ax3 ax4], '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]; + addtxt(ax1, '$\mathbf{\freq}$ = 1.25 Hz (Front)', ax1.YLim(2)); + addtxt(ax2, '$\mathbf{\freq}$ = 2.00 Hz (Front)', ax2.YLim(2)); + addtxt(ax3, '$\mathbf{\freq}$ = 1.25Hz (Aft)', ax3.YLim(2)); + addtxt(ax4, '$\mathbf{\freq}$ = 2.00 Hz (Aft)', ax4.YLim(2)); + pause(1); + figdir = 'figures/results/'; + filename = sprintf('half%s-AllF-Dx%0.2f_V%0.1f.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); + title(ax3, 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) + + % ylim([-1.5, 1.5]) + xlim([-13, 40]); + xticks([-10:10:40]); + + % + xlabel('Aft wing dihedral [deg]'); + ylabel('\cl and \cd [-]'); + + 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 -- GitLab