feat: add a loop on the lower rotor operation parameter to cancel the torque

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;