From d8f31c79347f7d7e1e6cd222e77c5489017580ff Mon Sep 17 00:00:00 2001
From: Robin Tamburrini <robin.tamburrini@student.uliege.be>
Date: Thu, 15 Feb 2024 16:31:59 +0100
Subject: [PATCH] feat: add a loop on the lower rotor operation parameter to
cancel the torque
---
src/solvers/bemt.m | 64 +++++++++++++++++++++++++++++++++++++++++++---
1 file changed, 60 insertions(+), 4 deletions(-)
diff --git a/src/solvers/bemt.m b/src/solvers/bemt.m
index 8924ea4..e596b1e 100644
--- a/src/solvers/bemt.m
+++ b/src/solvers/bemt.m
@@ -37,12 +37,66 @@ function bemt(OpRot, Mod)
% First rotor always calculated as if isolated
bemtsinglerot(OpRot(1), Mod);
+
+ % Deriving parameters of the second rotor to neutralize the overall torque
+ if length(OpRot) == 2 % The following torque cancellation method is only valid for N_ROTORS = 2
+
+ % Solve iteratively
+ isConverged = 0;
+ loopCount = 0;
+
+ % Parameters to vary to make both rotor torques coincide
+ strategy = Mod.coax.strategy;
+ if ~strcmp(strategy, 'rpm') && ~strcmp(strategy, 'coll')
+ error('bemt: the strategy must be''rpm'' or''coll''.');
+ end
+
+ % Lower and Upper bounds for bissection method
+ if strcmp(strategy, 'rpm')
+ OpRot_l = 1/10*OpRot(1).Op.(strategy);
+ OpRot_u = 10*OpRot(1).Op.(strategy);
+ else
+ OpRot_l = -10*OpRot(1).Op.(strategy);
+ OpRot_u = 10*OpRot(1).Op.(strategy);
+ end
+
+ while (~isConverged) && loopCount <= Mod.Num.maxIter
+
+ OpRot(2).Op.(strategy) = (OpRot_l + OpRot_u)/2; %Midpoint
+
+ OpRot(2).ElPerf.updateupstreamvel(OpRot(1), Mod); % Update upstream velocity
+ bemtsinglerot(OpRot(2), Mod); % Calculate BEMT as usual, with new upstream vel distribution
+
+ % Convergence criterion
+ if abs(OpRot(1).torque - OpRot(2).torque) < Mod.Num.convCrit
+ isConverged = 1;
+ else
+ % Update the bounds based on the torque
+ if OpRot(2).torque < OpRot(1).torque
+ OpRot_l = OpRot(2).Op.(strategy);
+ else
+ OpRot_u = OpRot(2).Op.(strategy);
+ end
+ end
+
+ if loopCount == Mod.Num.maxIter
+ warning('Rotare:Solvers:bemt:NotConverged', ...
+ ['Solution not converged for the coaxial case after reaching the '...
+ 'maximum number of iterations (%d iter).\n'], loopCount);
+ end
+ loopCount = loopCount + 1;
+ end
+
+ else % If there are more than 2 blades, the function can not support to cancel the total torque
+
+ warning('bemt: the total number of rotors in coaxial configuration must be equal to 2.')
- for i = 2:length(OpRot)
- OpRot(i).ElPerf.updateupstreamvel(OpRot(i - 1), Mod); % Update upstream velocity
- bemtsinglerot(OpRot(i), Mod); % Calculate BEMT as usual, with new upstream vel distribution
+ for i = 2:length(OpRot)
+ OpRot(i).ElPerf.updateupstreamvel(OpRot(i - 1), Mod); % Update upstream velocity
+ bemtsinglerot(OpRot(i), Mod); % Calculate BEMT as usual, with new upstream vel distribution
+ end
+
end
-
end
function bemtsinglerot(OpRot, Mod)
@@ -53,9 +107,11 @@ function bemtsinglerot(OpRot, Mod)
% Determine angles, velocities and forces for each element using the proper solver
solver = str2func(Mod.solver);
solver(OpRot, Mod);
+
% Calculate elemental forces, torque and power
OpRot.ElPerf.calcforces;
+
% Calculate overall rotor performance
OpRot.calcperf;
--
GitLab