diff --git a/src/solvers/bemt.m b/src/solvers/bemt.m index e596b1e64c26985f266f72d12210e2456414b477..fa64c77c125caef1f13387777f2e990a5545e604 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