From 61bb6acba7df0ecd6c1ceb5088a787a312fa9ddb Mon Sep 17 00:00:00 2001
From: Robin Tamburrini <robin.tamburrini@student.uliege.be>
Date: Wed, 21 Feb 2024 10:54:16 +0100
Subject: [PATCH] fix: the code was no longer adapted for a single rotor case
---
src/solvers/bemt.m | 101 ++++++++++++++++++++++-----------------------
1 file changed, 50 insertions(+), 51 deletions(-)
diff --git a/src/solvers/bemt.m b/src/solvers/bemt.m
index e596b1e..fa64c77 100644
--- a/src/solvers/bemt.m
+++ b/src/solvers/bemt.m
@@ -38,64 +38,63 @@ 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
+ for i = 2:length(OpRot)
- % 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;
+ % Deriving parameters of the second rotor to neutralize the overall torque
+ if Mod.coax.verif % 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
- % Update the bounds based on the torque
- if OpRot(2).torque < OpRot(1).torque
- OpRot_l = OpRot(2).Op.(strategy);
+ 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
- OpRot_u = OpRot(2).Op.(strategy);
+ % 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
-
- 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;
+
+ else % If there are more than 2 blades, the function can not support to cancel the total torque
+
+ 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
-
- 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
- end
-
end
end
--
GitLab