diff --git a/resultsanalysis.m b/resultsanalysis.m
index d00710830cdd5c26daee872ff062dccd45975763..496b369dc46c323c266c40fc6660a68e0331c11f 100644
--- a/resultsanalysis.m
+++ b/resultsanalysis.m
@@ -44,70 +44,69 @@ function ResData = resultsanalysis(matFiles)
ResData = ResData_;
end
- ResData = getallphases(ResData, OFFSET_VALUES);
+ % Get all correct useful phases
+ ResData = calcallphases(ResData);
+ ResData = calcinertia(ResData);
+ ResData = calctrueforces(ResData);
+ ResData = calccoeff(ResData);
+ %
- ResData = getmeancoefficients(ResData);
+ plotcf(ResData);
+ plotclcd(ResData);
end
-function ResData = getallphases(ResData, offsetVal)
+function ResData = calcallphases(ResData)
% GETALLPHASES Determine the range for all phase offsets
for i = 1:size(ResData, 1)
for j = 1:size(ResData, 2)
+ % FIXME Refactor the rest to use this!
+ ResData(i, j).Front.trueFreq = getmainfrequency(ResData(i, j).Front.angles, ...
+ 1 / diff(ResData(i, j).time(1:2)));
+ ResData(i, j).Aft.trueFreq = getmainfrequency(ResData(i, j).Aft.angles, ...
+ 1 / diff(ResData(i, j).time(1:2)));
+
fprintf('\nFinding phase for element (%d,%d): %0.2fm, %0.2fHz, %0.1fm/s\n', ...
i, j, ...
ResData(i, j).Testcase.dx, ...
ResData(i, j).Testcase.freq, ...
ResData(i, j).Testcase.airspeed);
- for k = 1:length(offsetVal)
- ResData(i, j).Phase(k).phaseVal = offsetVal(k);
- ResData(i, j).Phase(k).idx = findwingphase(ResData(i, j), offsetVal(k));
- end
+ % Find all phase offsets between front and aft
+ [offsetsVal, signalParts] = findallphasesdiff(ResData(i, j));
+ ResData(i, j).AllPhases.offsets = offsetsVal;
+ ResData(i, j).AllPhases.parts = signalParts;
end
end
end
-function ResData = getmeancoefficients(ResData)
- % GETMEANCOEFFICIENTS Returns the true mean coefficients
+function ResData = calcinertia(ResData)
+ % CALCINERTIA Calcualte inertia forces
for i = 1:size(ResData, 1)
for j = 1:size(ResData, 2)
- % Remove inertia forces
+ % Index of the Wind-off data
woIdx = (j) - mod(j - 1, 4);
- ResData(i, j) = removeinertia(ResData(i, j), ResData(i, woIdx), 'Front');
- ResData(i, j) = removeinertia(ResData(i, j), ResData(i, woIdx), 'Aft');
-
- % Transform forces/moments into coefficients
- ResData(i, j) = calccoeff(ResData(i, j));
- end
- end
-end
+ ResDataWO = ResData(i, woIdx);
-function ResData = calccoeff(ResData)
+ for iPos = {'Front', 'Aft'}
+ pos = string(iPos);
- POSITIONS = {'Front', 'Aft'};
- CHORD = 0.05;
- SPAN = 2 * 0.2;
- C_TO_K = 273.15;
+ % Get first and last peaks for wind-off conditions
+ woBounds = getboundingpeaks(ResDataWO.(pos).angles, ResDataWO.sampling);
+ woIdx = woBounds(1):woBounds(2);
- for i = 1:length(ResData.Phase)
- if ~isempty(ResData.Phase(i).idx)
- for j = 1:length(POSITIONS)
- dens = airdensity(ResData.Testcase.temp + C_TO_K, ResData.Testcase.press);
- speed = ResData.Testcase.airspeed;
+ % Calculate mean force and moments due to inertia
+ meanForces = mean(ResDataWO.(pos).forces(woIdx));
+ meanMom = mean(ResDataWO.(pos).moments(woIdx));
- coeffF = forcetocoeff(ResData.Phase(i).(POSITIONS{j}).forces, ...
- dens, speed, CHORD, SPAN, 'force');
- coeffM = forcetocoeff(ResData.Phase(i).(POSITIONS{j}).moments, ...
- dens, speed, CHORD, SPAN, 'moment');
-
- ResData.Phase(i).(POSITIONS{j}).cF = coeffF;
- ResData.Phase(i).(POSITIONS{j}).cM = coeffM;
+ ResData(i, j).AllPhases.(pos).inertiaF = meanForces;
+ ResData(i, j).AllPhases.(pos).inertiaM = meanMom;
end
+
end
end
end
diff --git a/utils/analysis/findallphasesdiff.m b/utils/analysis/findallphasesdiff.m
new file mode 100644
index 0000000000000000000000000000000000000000..cad6ce71745f47c29f6fe0472a5b6c5b59932ceb
--- /dev/null
+++ b/utils/analysis/findallphasesdiff.m
@@ -0,0 +1,54 @@
+function [offsetsVal, signalParts] = findallphasesdiff(ResData)
+ % FINDALLPHASESDIFF Return all phases differences between the two angles and their range
+ % Todo
+
+ % ----------------------------------------------------------------------------------------------
+ % (c) Copyright 2022 University of Liege
+ % Author: Thomas Lambert <t.lambert@uliege.be>
+ % ULiege - Aeroelasticity and Experimental Aerodynamics
+ % MIT License
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ % Abbreviation
+ period = 1 / ResData.Front.trueFreq;
+ nPeriodIdx = period * ResData.sampling;
+
+ % Find peaks in the wing angles
+ peakLocs1 = findallpeaks(ResData.Front.angles, nPeriodIdx);
+ peakLocs2 = findallpeaks(ResData.Aft.angles, nPeriodIdx);
+ npeaks = min(length(peakLocs1), length(peakLocs2));
+
+ % Trim to same number of peaks
+ peakLocs1 = peakLocs1(1:npeaks);
+ peakLocs2 = peakLocs2(1:npeaks);
+
+ % Get all offsets and round them (1° precision on the offset is plenty enough)
+ offsets = getphaseoffsets(peakLocs1, peakLocs2, nPeriodIdx);
+ offsets = round(offsets);
+ % Remove first and last offsets (can't get one period before and one after the offset otherwise)
+ offsetsVal = offsets(2:end - 1);
+ peakLocs = peakLocs1(2:end - 1);
+
+ % singalPart is one period centered on the peak
+ % signalParts = [peakLocs(1:end - 2), peakLocs1(3:end)];
+ signalParts = round([peakLocs - nPeriodIdx / 2, peakLocs + nPeriodIdx / 2]);
+
+end
+
+function peakLocs = findallpeaks(signal, nPeriodIdx)
+ % FINDANGLEPEAKS Find all true peaks in the signal
+
+ maxAngle = max(signal);
+
+ [~, peakLocs] = findpeaks(signal, ...
+ 'MinPeakProminence', maxAngle / 2, ...
+ 'MinPeakDistance', nPeriodIdx / 2);
+end
+
+function offsets = getphaseoffsets(peaks1, peaks2, nPeriodIdx)
+ % GETPHASEOFFSETS Calculate phase offsets between every contiguous peaks
+ % Note: Offset is positive when peak2 is AFTER peak1
+
+ offsets = (peaks2 - peaks1) * 360 / nPeriodIdx;
+
+end
diff --git a/utils/analysis/plotcf.m b/utils/analysis/plotcf.m
new file mode 100644
index 0000000000000000000000000000000000000000..a3b54f1e197074a846e705c4818cba1b6fc25d18
--- /dev/null
+++ b/utils/analysis/plotcf.m
@@ -0,0 +1,93 @@
+function plotcf(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
+
+ ONLY_ERR_BARS = true;
+ OFFSET_STEPS = 20;
+
+ INTERESTING_PLOTS = [
+ ];
+
+ % for i = 1: size(ResData,1)
+ for i = 3
+ for j = 1:size(ResData, 2)
+ if ResData(i, j).Testcase.airspeed ~= 0 % No need to plot when cF = 0
+ expParam = sprintf('d = %0.2f, f = %0.2f, V = %0.1f', ...
+ ResData(i, j).Testcase.dx, ...
+ ResData(i, j).Testcase.freq, ...
+ ResData(i, j).Testcase.airspeed);
+
+ offsets = ResData(i, j).AllPhases.offsets;
+ cleanOffsets = -175:OFFSET_STEPS:175;
+ 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');
+ subplot(211);
+ hold on;
+ if ~ONLY_ERR_BARS
+ plot(ResData(i, j).AllPhases.offsets, ...
+ ResData(i, j).AllPhases.Front.cF(:, [1, 3]), 'o');
+ end
+ errorbar(cleanOffsets, meanF(:, 1), ynegF(:, 1), yposF(:, 1), xneg, xpos, 'o');
+ errorbar(cleanOffsets, meanF(:, 3), ynegF(:, 3), yposF(:, 3), xneg, xpos, 'o');
+ hold off;
+ title(['Front module -- ', expParam]);
+
+ setgca();
+
+ subplot(212);
+ hold on;
+ if ~ONLY_ERR_BARS
+ plot(ResData(i, j).AllPhases.offsets, ...
+ ResData(i, j).AllPhases.Aft.cF(:, [1, 3]), 'o');
+ end
+ errorbar(cleanOffsets, meanA(:, 1), ynegA(:, 1), yposA(:, 1), xneg, xpos, 'o');
+ errorbar(cleanOffsets, meanA(:, 3), ynegA(:, 3), yposA(:, 3), xneg, xpos, 'o');
+ hold off;
+ title(['Aft module -- ', expParam]);
+ setgca();
+ legend('Location', 'bestoutside', 'orientation', 'horizontal', 'Fx', 'Fz');
+
+ end
+ end
+ end
+
+end
+
+function setgca()
+ % SETGCA Sets the axes parameters
+ xlim([-180, 180]);
+ grid on;
+end
diff --git a/utils/analysis/plotclcd.m b/utils/analysis/plotclcd.m
new file mode 100644
index 0000000000000000000000000000000000000000..c0f1b65bcf1baa80ac7c9d2c1d6fd4557358d5f5
--- /dev/null
+++ b/utils/analysis/plotclcd.m
@@ -0,0 +1,84 @@
+function plotclcd(ResData)
+ % PLOTCLCD Plot CL/CD for each set of wing
+ % 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
+
+ % for i = 1: size(ResData,1)
+ for i = 2
+ for j = 1:size(ResData, 2)
+ if ResData(i, j).Testcase.airspeed ~= 0 % No need to plot when cF = 0
+ expParam = sprintf('d = %0.2f, f = %0.2f, V = %0.1f', ...
+ ResData(i, j).Testcase.dx, ...
+ ResData(i, j).Testcase.freq, ...
+ ResData(i, j).Testcase.airspeed);
+
+ clcdF = ResData(i, j).AllPhases.Front.cF(:, 3) ./ ...
+ ResData(i, j).AllPhases.Front.cF(:, 1);
+ clcdA = ResData(i, j).AllPhases.Aft.cF(:, 3) ./ ...
+ ResData(i, j).AllPhases.Aft.cF(:, 1);
+
+ offsets = ResData(i, j).AllPhases.offsets;
+ cleanOffsets = -175:OFFSET_STEPS:175;
+ 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));
+ meanClCdF(iOff) = mean(clcdF(offIdx), 1);
+ stdClCdF(iOff) = std(clcdF(offIdx), 0, 1);
+ meanClCdA(iOff) = mean(clcdA(offIdx), 1);
+ stdClCdA(iOff) = std(clcdA(offIdx), 0, 1);
+ end
+
+ ynegF = -stdClCdF;
+ yposF = stdClCdF;
+ ynegA = -stdClCdA;
+ yposA = stdClCdA;
+
+ xneg = cleanOffsets - meanOffset - stdOffset;
+ xpos = cleanOffsets - meanOffset + stdOffset;
+
+ figure('Name', 'errorbars');
+ subplot(211);
+ hold on;
+ plot(offsets, clcdF, 'o');
+ errorbar(cleanOffsets, meanClCdF, ynegF, yposF, xneg, xpos, 'o');
+ hold off;
+ title(['Front module -- ', expParam]);
+ setgca();
+
+ subplot(212);
+ hold on;
+ plot(offsets, clcdA, 'o');
+ errorbar(cleanOffsets, meanClCdA, ynegA, yposA, xneg, xpos, 'o');
+ hold off;
+ title(['Aft module -- ', expParam]);
+ setgca();
+
+ pause;
+
+ end
+ end
+ end
+
+end
+
+function setgca()
+ % SETGCA Sets the axes parameters
+ xlim([-180, 180]);
+ % ylim([-1.5, 1.5]);
+ grid on;
+end
diff --git a/utils/calccoeff.m b/utils/calccoeff.m
new file mode 100644
index 0000000000000000000000000000000000000000..75cfef0066f2ed5a33cd9cee629864f4b89d9d27
--- /dev/null
+++ b/utils/calccoeff.m
@@ -0,0 +1,38 @@
+function [ResData] = calccoeff(ResData)
+ % CALCCOEFF Calculates force and moment coefficients.
+ % TODO
+
+ % ----------------------------------------------------------------------------------------------
+ % (c) Copyright 2022 University of Liege
+ % Author: Thomas Lambert <t.lambert@uliege.be>
+ % ULiege - Aeroelasticity and Experimental Aerodynamics
+ % MIT License
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ CHORD = 0.05;
+ SPAN = 2 * 0.2;
+ C_TO_K = 273.15;
+
+ for i = 1:size(ResData, 1)
+ for j = 1:size(ResData, 2)
+ for iPos = {'Front', 'Aft'}
+ pos = string(iPos);
+ dens = airdensity(ResData(i, j).Testcase.temp + C_TO_K, ...
+ ResData(i, j).Testcase.press);
+ speed = ResData(i, j).Testcase.airspeed;
+
+ for iPart = 1:size(ResData(i, j).AllPhases(1).(pos).meanF, 1)
+
+ coeffF = forcetocoeff(ResData(i, j).AllPhases(1).(pos).meanF(iPart, :), ...
+ dens, speed, CHORD, SPAN, 'force');
+ coeffM = forcetocoeff(ResData(i, j).AllPhases(1).(pos).meanM(iPart, :), ...
+ dens, speed, CHORD, SPAN, 'moment');
+
+ ResData(i, j).AllPhases(1).(pos).cF(iPart, :) = coeffF;
+ ResData(i, j).AllPhases(1).(pos).cF(iPart, :) = coeffM;
+ end
+ end
+ end
+ end
+
+end
diff --git a/utils/calctrueforces.m b/utils/calctrueforces.m
new file mode 100644
index 0000000000000000000000000000000000000000..35c966765baf660dabcd3cc5377f66244851470c
--- /dev/null
+++ b/utils/calctrueforces.m
@@ -0,0 +1,43 @@
+function [ResData] = calctrueforces(ResData)
+ % CALCTRUEFORCES Calculates true forces by substracting inertia from mean raw measure.
+ % TODO
+
+ % ----------------------------------------------------------------------------------------------
+ % (c) Copyright 2022 University of Liege
+ % Author: Thomas Lambert <t.lambert@uliege.be>
+ % ULiege - Aeroelasticity and Experimental Aerodynamics
+ % MIT License
+ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+ for i = 1:size(ResData, 1)
+ for j = 1:size(ResData, 2)
+ for iPos = {'Front', 'Aft'}
+ pos = string(iPos);
+ indexes = ResData(i, j).AllPhases.parts;
+
+ for iPart = 1:length(indexes)
+
+ % Do not remove inertia from inertial case!
+ if ResData(i, j).Testcase.airspeed ~= 0
+ % Inertia
+ inF = ResData(i, j).AllPhases(1).(pos).inertiaF;
+ inM = ResData(i, j).AllPhases(1).(pos).inertiaM;
+ else
+ inF = 0;
+ inM = 0;
+ end
+ % Raw force
+ meanFraw = mean(ResData(i, j).(pos).forces(indexes(iPart, :), :));
+ meanMraw = mean(ResData(i, j).(pos).moments(indexes(iPart, :), :));
+
+ % True value
+ ResData(i, j).AllPhases(1).(pos).meanF(iPart, :) = meanFraw - inF;
+ ResData(i, j).AllPhases(1).(pos).meanM(iPart, :) = meanMraw - inM;
+
+ end
+
+ end
+ end
+ end
+
+end