diff --git a/+PropulsionPkg/Contents.m b/+PropulsionPkg/Contents.m index 3b00c8f..7311193 100644 --- a/+PropulsionPkg/Contents.m +++ b/+PropulsionPkg/Contents.m @@ -16,4 +16,5 @@ % PropulsionSizing - PropulsionPkg.PropulsionSizing is a function. % RecomputeSplits - PropulsionPkg.RecomputeSplits is a function. % TestCreatePropArch - PropulsionPkg.TestCreatePropArch is a function. +% TestEngineMotorRegressions - PropulsionPkg.TestEngineMotorRegressions is a function. % TestPowerAvailable - PropulsionPkg.TestPowerAvailable is a function. diff --git a/+PropulsionPkg/PowerSupplementCheck.m b/+PropulsionPkg/PowerSupplementCheck.m index feb5e3d..6e99fb3 100644 --- a/+PropulsionPkg/PowerSupplementCheck.m +++ b/+PropulsionPkg/PowerSupplementCheck.m @@ -1,8 +1,8 @@ function [Psupp] = PowerSupplementCheck(Preq, Arch, Lambda, Eta, TrnType, EtaFan) % -% [Psupp] = PowerSupplementCheck(Preq, Arch, Lambda, Eta, TrnType, EtaFan, itrn) +% [Psupp] = PowerSupplementCheck(Preq, Arch, Lambda, Eta, TrnType, EtaFan) % written by Paul Mokotoff, prmoko@umich.edu -% last updated: 04 aug 2025 +% last updated: 16 sep 2026 % % In the propulsion architecture, check if any components are either % suppling/siphoning power from the gas-turbine engines. If a component is @@ -100,23 +100,36 @@ % find the gas-turbine engine being supplemented Driving = find((Arch(:, icomp) > 0)' & (TrnType == 1)); - % check if there are multiple (or none) "driving" gas-turbines - if (length(Driving) ~= 1) - - % don't know how much each motor is powering each gas-turbine + % a shared target needs both an engine and a motor contribution + Helping = find((Arch(:, icomp) > 0)' & (TrnType == 0)); + if (isempty(Driving) || isempty(Helping)) continue; - end - - % find the electric motors that are supplementing - Helping = find((Arch(:, icomp) > 0)' & (TrnType == 0)); - - % add the power supplement, accounting for the fan efficiency - Psupp(:, Driving) = Psupp(:, Driving) + sum(Preq(:, Helping), 2) .* EtaFan; %#ok, ignore warning about "find" ... easier to read this way + + % Downstream splits describe fractions of the target demand. Motor + % power at the target already includes path losses in Preq, so using + % motor input power here would apply the same split a second time. + MotorPower = Preq(:, icomp) .* sum(Lambda(icomp, Helping)); + + % Only a fan target needs the engine fan-efficiency conversion. + TargetEfficiency = 1; + if (TrnType(icomp) == 2) + TargetEfficiency = EtaFan; + end + + % When engines share the target, allocate the motor credit in + % proportion to their own contributions to that target. + EngineShares = Lambda(icomp, Driving); + if (sum(EngineShares) <= 0) + EngineShares = ones(size(EngineShares)); + end + Psupp(:, Driving) = Psupp(:, Driving) + ... + MotorPower .* TargetEfficiency .* ... + (EngineShares ./ sum(EngineShares)); end end % ---------------------------------------------------------- -end \ No newline at end of file +end diff --git a/+PropulsionPkg/ProcessPropArch.m b/+PropulsionPkg/ProcessPropArch.m index 4a90ced..5d08d79 100644 --- a/+PropulsionPkg/ProcessPropArch.m +++ b/+PropulsionPkg/ProcessPropArch.m @@ -2,7 +2,7 @@ % % [Aircraft] = ProcessPropArch(Aircraft) % written by Paul Mokotoff, prmoko@umich.edu -% last updated: 10 mar 2026 +% last updated: 16 sep 2026 % % given a propulsion architecture, find how each gas turbine engine is % connected to a propeller. the power required at the propeller is @@ -55,7 +55,7 @@ EtaDwn = ones(ncomp, ncomp); % allocate memory for indexed transmitter arrays -WhichProp = zeros(1, ntrn); +WhichProp = repmat({[]}, 1, ntrn); HEcoeff = zeros(1, ntrn); % create logical array transmitters @@ -75,53 +75,15 @@ Preq = zeros(ncomp, 1); PSLS = zeros(ncomp, 1); -% loop through all propellers to find indirect gas turbine engines +% Propagate each propeller's demand through every upstream path. A direct +% engine connection must not hide a second engine reached through a cable. for iprop = PropIdx - - % remember a copy of the index - jprop = iprop; - - % remember the index - krow = jprop; - - % loop until gas turbine engine is found - while (~isempty(krow)) - - % check if any are gas turbine engines - AnyGTE = logical(sum(Arch(:, krow) > 0 & ID == 1 & itrn, 2)); - - if (any(AnyGTE)) % assume only 1 GTE per fan - - % get its index - igte = find(AnyGTE); - - % perturb the power vector - PowerVector(iprop) = SLSPower(iprop - nsrc); - - % propagate power downstream - Pout = PropulsionPkg.PowerFlow(PowerVector, Arch', LamDwn, EtaDwn, -1, 1.0e-06); - - % remember the power output by the engine and the total power - % required - Preq(igte) = Preq(igte) + Pout(igte); - PSLS(igte) = PSLS(igte) + SLSPower(iprop - nsrc); - - % remove the perturbation - PowerVector(iprop) = 0; - - % break out of the loop - break; - - else - - % remember the indices - jrow = krow; - - % search a level deeper - [krow, ~] = find(Arch(:, jrow)); - - end - end + PowerVector(iprop) = SLSPower(iprop - nsrc); + Pout = PropulsionPkg.PowerFlow(PowerVector, Arch', LamDwn, EtaDwn, -1, 1.0e-06); + ActiveEngines = EngIdx(Pout(EngIdx) > 0); + Preq(ActiveEngines) = Preq(ActiveEngines) + Pout(ActiveEngines); + PSLS(ActiveEngines) = PSLS(ActiveEngines) + SLSPower(iprop - nsrc); + PowerVector(iprop) = 0; end % loop through all propellers to find directly connected gas turbine engines @@ -134,7 +96,9 @@ if (~isempty(GTEIdx)) % remember the index - WhichProp(GTEIdx - nsrc) = iprop; + for igte = GTEIdx(:)' + WhichProp{igte - nsrc}(end + 1) = iprop; + end end @@ -150,4 +114,4 @@ Aircraft.Specs.Propulsion.Engine.HEcoeff = HEcoeff; Aircraft.Specs.Propulsion.PropArch.WhichProp = WhichProp; -end \ No newline at end of file +end diff --git a/+PropulsionPkg/PropAnalysis.m b/+PropulsionPkg/PropAnalysis.m index 24c236f..7f0b71a 100644 --- a/+PropulsionPkg/PropAnalysis.m +++ b/+PropulsionPkg/PropAnalysis.m @@ -486,13 +486,13 @@ icol = HasEng(ieng) + nsrc; % check if it has a propeller - iprop = WhichProp(HasEng(ieng)); + iprop = WhichProp{HasEng(ieng)}; % check if the engine is connected to a propeller - if (iprop ~= 0) + if (~isempty(iprop)) % get the thrust requirement from the propeller - TEng = Tout(ibeg:iend, iprop); + TEng = sum(Tout(ibeg:iend, iprop), 2); else @@ -712,4 +712,4 @@ % ---------------------------------------------------------- -end \ No newline at end of file +end diff --git a/+PropulsionPkg/PropulsionSizing.m b/+PropulsionPkg/PropulsionSizing.m index ee3f97d..426168c 100644 --- a/+PropulsionPkg/PropulsionSizing.m +++ b/+PropulsionPkg/PropulsionSizing.m @@ -195,7 +195,7 @@ if (Aircraft.Specs.Propulsion.Engine.DesignThrust > 1.0e-06) % size the engine - Engine = EngineModelPkg.TurbofanNonlinearSizing(Aircraft.Specs.Propulsion.Engine, Psupp(ieng(1))); + Engine = EngineModelPkg.TurbofanNonlinearSizing(Aircraft.Specs.Propulsion.Engine, Psupp(ieng(jeng))); % turn off engine sizing Engine.Specs.Sizing = 0; @@ -207,7 +207,7 @@ InletArea = pi * Engine.FanDiam ^ 2 / 4; % find the fan connected to the engine - ifan = find((Arch(jeng+nsrc, idx) == 1) & Prop); + ifan = find((Arch(ieng(jeng)+nsrc, idx) == 1) & Prop); % remember the inlet area Aircraft.Specs.Propulsion.InletArea(ifan) = InletArea; @@ -319,4 +319,4 @@ % ---------------------------------------------------------- -end \ No newline at end of file +end diff --git a/+PropulsionPkg/RecomputeSplits.m b/+PropulsionPkg/RecomputeSplits.m index e5daf8c..f1d8dc0 100644 --- a/+PropulsionPkg/RecomputeSplits.m +++ b/+PropulsionPkg/RecomputeSplits.m @@ -2,14 +2,13 @@ % % [Aircraft] = RecomputeSplits(Aircraft, SegBeg, SegEnd) % written by Paul Mokotoff, prmoko@umich.edu -% last updated: 03 sep 2025 +% last updated: 16 sep 2026 % % Re-compute the operational power splits for a "full throttle" setting % during the mission. % -% this function previously worked for two elements connected in parallel. -% it has now been generalized to work for any number of elements connected -% in parallel. +% Each parallel target is recomputed from the power actually delivered to +% that target. Separate targets may share an engine without sharing a split. % % INPUTS: % Aircraft - structure with information about the aircraft and mission @@ -53,11 +52,9 @@ % get the propulsion architecture Arch = Aircraft.Specs.Propulsion.PropArch.Arch; -% get the number of parallel connections -npar = length(ParIndx); - % get the number of sources and transmitters nsrc = length(Aircraft.Specs.Propulsion.PropArch.SrcType); +TrnType = Aircraft.Specs.Propulsion.PropArch.TrnType; % get the power available (equal to power output for "full throttle" case) Pav = Aircraft.Mission.History.SI.Power.Pav(SegBeg:SegEnd, :); @@ -72,56 +69,64 @@ % get the number of downstream splits nsplit = length(Aircraft.Specs.Power.LamDwn.SLS); -% get a temporary power split -TmpSplit = LamDwn(1, :); +OperUps = Aircraft.Specs.Propulsion.PropArch.OperUps; +OperDwn = Aircraft.Specs.Propulsion.PropArch.OperDwn; +EtaUps = Aircraft.Specs.Propulsion.PropArch.EtaUps; + +% Collect each parallel target's incoming edges. An engine can contribute +% to several targets, so its helpers cannot be treated as one group. +Source = []; +Target = []; +for itarget = nsrc + (1:length(TrnType)) + Parents = find(Arch(:, itarget))'; + Transmitters = Parents(Parents > nsrc & Parents <= nsrc + length(TrnType)); + Engine = Transmitters(TrnType(Transmitters - nsrc) == 1); + Motor = Transmitters(TrnType(Transmitters - nsrc) == 0); + if (~isempty(Engine) && ~isempty(Motor)) + Source = [Source, Parents]; + Target = [Target, repmat(itarget, 1, length(Parents))]; + end +end -% get the original downstream matrix -OperDwn = PropulsionPkg.EvalSplit(Aircraft.Specs.Propulsion.PropArch.OperDwn, TmpSplit); +if (isempty(Source) || nsplit == 0) + return +end -% loop through each power split -for ipar = 1:npar - - % get the index of the main connection - imain = ParIndx(ipar); - - % get the supplemental connection(s) - isupp = ParConns{imain}; - - % account for the source indices - imain = imain + nsrc; - - % get all indices - jdx = [imain, isupp]; - - % find the upstream split - iups = find(sum(Arch(jdx, :), 1) == length(isupp) + 1); - - % get the total power output at any given time from those sources - Pout = sum(Pav(idx, jdx), 2); - - % find the downstream split contributing - for kdx = jdx - for isplit = 1:nsplit - - % perturb a split - TmpSplit(isplit) = TmpSplit(isplit) + 0.01; - - % get the new matrix - OperNew = PropulsionPkg.EvalSplit(Aircraft.Specs.Propulsion.PropArch.OperDwn, TmpSplit); - - % check if the matrices are different - if (abs(OperDwn(iups, kdx) - OperNew(iups, kdx)) > 1.0e-06) - - % recompute the power split - LamDwn(idx, isplit) = Pav(idx, kdx) ./ Pout; - - end - - % remove the perturbation - TmpSplit(isplit) = TmpSplit(isplit) - 0.01; - - end - end +ncomp = size(Arch, 1); +DownIndex = sub2ind([ncomp, ncomp], Target, Source); +UpIndex = sub2ind([ncomp, ncomp], Source, Target); +for ipoint = find(idx)' + Active = Pav(ipoint, Target) > 0; + if (~any(Active)) + continue + end + Current = LamDwn(ipoint, :); + UpMatrix = PropulsionPkg.EvalSplit(OperUps, LamUps(ipoint, :)); + DownMatrix = PropulsionPkg.EvalSplit(OperDwn, Current); + EdgePower = Pav(ipoint, Source(Active))'; + TargetPower = Pav(ipoint, Target(Active))'; + UpFraction = UpMatrix(UpIndex(Active)); + UpEfficiency = EtaUps(UpIndex(Active)); + Desired = EdgePower .* UpFraction(:) .* UpEfficiency(:) ./ TargetPower; + Baseline = DownMatrix(DownIndex(Active)); + Baseline = Baseline(:); + Sensitivity = zeros(sum(Active), nsplit); + for isplit = 1:nsplit + Perturbed = Current; + Perturbed(isplit) = Perturbed(isplit) + 0.01; + NewMatrix = PropulsionPkg.EvalSplit(OperDwn, Perturbed); + NewFraction = NewMatrix(DownIndex(Active)); + Sensitivity(:, isplit) = (NewFraction(:) - Baseline) / 0.01; + end + Change = Sensitivity \ (Desired(:) - Baseline(:)); + Updated = Current + Change'; + NewMatrix = PropulsionPkg.EvalSplit(OperDwn, Updated); + Actual = NewMatrix(DownIndex(Active)); + if (any(~isfinite(Updated)) || any(abs(Actual(:) - Desired(:)) > 1.0e-06)) + error("FAST:EdgePowerInconsistent", ... + "Edge power is not consistent. Your power splits are over-constrained."); + end + LamDwn(ipoint, :) = Updated; end % if any are NaN, return 0 (assume it's from 0 power available) @@ -132,4 +137,4 @@ % ---------------------------------------------------------- -end \ No newline at end of file +end diff --git a/+PropulsionPkg/TestEngineMotorRegressions.m b/+PropulsionPkg/TestEngineMotorRegressions.m new file mode 100644 index 0000000..052b05e --- /dev/null +++ b/+PropulsionPkg/TestEngineMotorRegressions.m @@ -0,0 +1,744 @@ +function [Success] = TestEngineMotorRegressions() +% +% [Success] = TestEngineMotorRegressions() +% written by Triet Ho +% last updated: 16 sep 2026 +% +% Check gas-turbine and electric-motor accounting on controlled propulsion +% architectures. Include physically consistent split and efficiency cases, +% since arbitrary test matrices can conceal edge-power attribution errors. +% +% INPUTS: +% none +% +% OUTPUTS: +% Success - 1 when all 33 checks pass, otherwise 0. +% size/type/units: 1-by-1 / int / [] +% + +%% TEST CASE SETUP %% +%%%%%%%%%%%%%%%%%%%%% + +% Keep one result per case so a failure identifies the affected behavior. +Pass = false(33, 1); + +%% ENGINE TO PROPELLER CONNECTIONS %% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% one engine driving one propeller +aircraft = PropulsionPkg.ProcessPropArch(makeAircraft([1])); +Pass(1) = isequal(propellerConnections(aircraft, 1), 3); + +% one engine driving two propellers must retain both connections +aircraft = PropulsionPkg.ProcessPropArch(makeAircraft([1, 1])); +Pass(2) = isequal(propellerConnections(aircraft, 1), [3, 4]); + +% the connection list must also retain a third propeller +aircraft = PropulsionPkg.ProcessPropArch(makeAircraft([1, 1, 1])); +Pass(3) = isequal(propellerConnections(aircraft, 1), [3, 4, 5]); + +% separate engines must not acquire each other's propeller +aircraft = PropulsionPkg.ProcessPropArch(makeAircraft([1, 2])); +Pass(4) = isequal(propellerConnections(aircraft, 1), 4) && ... + isequal(propellerConnections(aircraft, 2), 5); + +%% PER ENGINE MOTOR SUPPLEMENTS %% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% different motor shares must produce distinct per-engine supplements +Pass(5) = CheckEngineSupplements(0.2, 0.4); + +% reversing the shares must not reuse the first engine's supplement +Pass(6) = CheckEngineSupplements(0.4, 0.2); + +%% TARGET EFFICIENCY AND MOTOR CREDIT %% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% fan target applies fan efficiency to the motor contribution +Pass(7) = CheckValue(parallelSupplement([1, 0, 2], [0, 10, 0], 0.8), 8); + +% generator and cable targets do not incur a fan-efficiency loss +Pass(8) = CheckValue(parallelSupplement([1, 0, 3], [0, 10, 0], 0.8), 10); +Pass(9) = CheckValue(parallelSupplement([1, 0, 4], [0, 10, 0], 0.65), 10); + +% contributions from two motors must both be counted at a shared target +Pass(10) = CheckValue(parallelSupplement([1, 0, 0, 3], [0, 7, 11, 0], 0.8), 18); + +% unity fan efficiency leaves the motor contribution unchanged +Pass(11) = CheckValue(parallelSupplement([1, 0, 3], [0, 10, 0], 1), 10); + +% only one quarter of the first fan's demand comes from the motor +architecture = zeros(4); +architecture(1, 3) = 1; +architecture(2, 3:4) = 1; +splits = zeros(4); +splits(3, 1) = 0.75; +splits(3, 2) = 0.25; +splits(4, 2) = 1; +actual = PropulsionPkg.PowerSupplementCheck([15, 20, 20, 15], architecture, splits, ones(4), [1, 0, 2, 2], 0.8); +Pass(12) = CheckValue(actual(1), 4); + +% a half-efficient motor path delivers half of its available power +architecture = zeros(3); +architecture(1:2, 3) = 1; +splits = zeros(3); +splits(3, 1:2) = 0.5; +efficiencies = ones(3); +efficiencies(3, 2) = 0.5; +actual = PropulsionPkg.PowerSupplementCheck([10, 20, 20], architecture, splits, efficiencies, [1, 0, 2], 0.8); +Pass(13) = CheckValue(actual(1), 8); + +% shared targets must credit each split motor contribution only once +architecture = zeros(4); +architecture(1:2, 3:4) = 1; +splits = zeros(4); +splits(3:4, 1:2) = 0.5; +actual = PropulsionPkg.PowerSupplementCheck([20, 20, 20, 20], architecture, splits, ones(4), [1, 0, 2, 2], 0.8); +Pass(14) = CheckValue(actual(1), 16); + +%% EDGE-CONSISTENT PARALLEL HYBRIDS %% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% A 50/50 fan demand has 50 W of delivered motor power, or 49.5 W after +% fan efficiency. A downstream split is normalized by fan demand, not by +% motor power; applying it to both quantities would count the split twice. +Pass(15) = CheckValue(consistentParallelSupplement(0.99, 0.99), 49.5); + +% The graph's fan-edge efficiency can differ from the engine model's +% EtaFan. Only the latter reduces the delivered motor power again. +Pass(16) = CheckValue(consistentParallelSupplement(0.8, 0.99), 49.5); + +% One motor sends 25 W to a 100 W fan and 75 W to a 200 W fan. Their +% downstream shares differ from the motor's 25/75 outgoing allocation. +architecture = zeros(4); +architecture(1:2, 3:4) = 1; +splits = zeros(4); +splits(3, 1:2) = [0.75, 0.25]; +splits(4, 1:2) = [0.625, 0.375]; +efficiencies = ones(4); +efficiencies(3:4, 1:2) = 0.8; +actual = PropulsionPkg.PowerSupplementCheck( ... + [250, 125, 100, 200], architecture, splits, efficiencies, [1, 0, 2, 2], 0.99); +Pass(17) = CheckValue(actual(1), 99); + +% If two TSEs and one motor feed one fan, silently skipping the motor +% loses 49.5 W of assistance after fan efficiency. The total credit +% should conserve this contribution regardless of its engine allocation. +architecture = zeros(4); +architecture(1:3, 4) = 1; +splits = zeros(4); +splits(4, 1:3) = [0.25, 0.25, 0.5]; +efficiencies = ones(4); +efficiencies(4, 1:3) = 0.8; +actual = PropulsionPkg.PowerSupplementCheck( ... + [31.25, 31.25, 62.5, 100], architecture, splits, efficiencies, [1, 1, 0, 2], 0.99); +Pass(18) = CheckValue(sum(actual(1:2)), 49.5); + +%% GRAPH INDEXING AND FULL-THROTTLE SPLITS %% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% User-supplied graphs need not place engines before motors. Each sized +% engine must assign its inlet area to its own fan in an interleaved graph. +Pass(19) = CheckInterleavedEngineInletAreas(); + +% One engine can share two different fans with two different motors. +% Full-throttle split recomputation must not assume one common target. +Pass(20) = CheckSeparateParallelTargets(); + +% A direct engine and an engine connected through a cable can both feed +% one propeller. Finding the direct engine must not hide the indirect one. +Pass(21) = CheckMixedDepthEngineConnections(); + +%% MULTI-ENGINE ARCHITECTURES %% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +% Exercise more than two engines at increasing graph sizes. Each count +% checks distinct propellers, one shared fan, interleaved sizing, and +% direct plus indirect paths to a common propeller. +EngineCounts = [3, 5, 8]; +for icount = 1:length(EngineCounts) + engineCount = EngineCounts(icount); + caseIndex = 22 + 4 * (icount - 1); + Pass(caseIndex) = CheckManyEngineConnections(engineCount); + Pass(caseIndex + 1) = CheckManyEngineMotorCredit(engineCount); + Pass(caseIndex + 2) = CheckManyInterleavedEngineInletAreas(engineCount); + Pass(caseIndex + 3) = CheckManyMixedDepthEngineConnections(engineCount); +end + +%% CHECK THE TEST RESULTS %% +%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +Success = all(Pass); +if (Success) + fprintf(1, "EngineMotorRegressions tests passed!\n"); +else + fprintf(1, "EngineMotorRegressions tests failed:\n"); + fprintf(1, " Test %d\n", find(~Pass)); +end + +end + +function [Pass] = CheckEngineSupplements(firstMotorShare, secondMotorShare) +% +% [Pass] = CheckEngineSupplements(firstMotorShare, secondMotorShare) +% written by Triet Ho +% last updated: 16 sep 2026 +% +% Confirm that sizing gives each engine its own electric-motor supplement. +% The supplements must differ enough to expose accidental reuse. +% +% INPUTS: +% firstMotorShare - fraction of the first engine target fed by its motor. +% size/type/units: 1-by-1 / double / [] +% secondMotorShare - corresponding fraction for the second engine. +% size/type/units: 1-by-1 / double / [] +% +% OUTPUTS: +% Pass - true if both sized engines receive their own load. +% size/type/units: 1-by-1 / logical / [] +% + +aircraft = makeTwoEngineSizingAircraft(firstMotorShare, secondMotorShare); +aircraft = PropulsionPkg.PropulsionSizing(aircraft); +expected = aircraft.Specs.Propulsion.PowerSupp(1:2); +actual = arrayfun(@(engine) engine.FanSysObject.ElecWork, aircraft.Specs.Propulsion.SizedEngine); +actual = actual(:)'; +Pass = abs(diff(expected)) > 1e5 && ... + isequal(size(actual), size(expected)) && ... + all(isfinite(actual)) && all(isfinite(expected)) && ... + all(abs(actual - expected) <= 1e-8 * max(1, abs(expected))); +end + +function [Pass] = CheckValue(actual, expected) +% +% [Pass] = CheckValue(actual, expected) +% written by Triet Ho +% last updated: 16 sep 2026 +% +% Check a finite scalar power contribution within a 1e-10 W tolerance. +% +% INPUTS: +% actual - calculated supplemental power. +% size/type/units: 1-by-1 / double / [W] +% expected - expected supplemental power. +% size/type/units: 1-by-1 / double / [W] +% +% OUTPUTS: +% Pass - true if the values agree within tolerance. +% size/type/units: 1-by-1 / logical / [] +% + +Pass = isscalar(actual) && isfinite(actual) && abs(actual - expected) <= 1e-10; +end + +function aircraft = makeTwoEngineSizingAircraft(firstMotorShare, secondMotorShare) +% +% aircraft = makeTwoEngineSizingAircraft(firstMotorShare, secondMotorShare) +% written by Triet Ho +% last updated: 16 sep 2026 +% +% Build two otherwise matched engines with independently adjustable motor +% loads. The graph includes two engines, two motors, two fans, and a sink. +% +% INPUTS: +% firstMotorShare - fraction of the first fan power from its motor. +% size/type/units: 1-by-1 / double / [] +% secondMotorShare - fraction of the second fan power from its motor. +% size/type/units: 1-by-1 / double / [] +% +% OUTPUTS: +% aircraft - ERJ175LR structure with the test architecture. +% size/type/units: 1-by-1 / struct / [] +% + +aircraft = DataStructPkg.PreSpecProcessing(AircraftSpecsPkg.ERJ175LR()); +aircraft = DataStructPkg.SpecProcessing(aircraft); +architecture = zeros(9); +architecture(1, 3:4) = 1; +architecture(2, 5:6) = 1; +architecture(3, 7) = 1; +architecture(5, 7) = 1; +architecture(4, 8) = 1; +architecture(6, 8) = 1; +architecture(7:8, 9) = 1; +splits = ones(9); +splits(7, 3) = 1 - firstMotorShare; +splits(7, 5) = firstMotorShare; +splits(8, 4) = 1 - secondMotorShare; +splits(8, 6) = secondMotorShare; +splits(9, 7:8) = 0.5; +aircraft.Specs.Propulsion.PropArch.Arch = architecture; +aircraft.Specs.Propulsion.PropArch.SrcType = [1, 0]; +aircraft.Specs.Propulsion.PropArch.TrnType = [1, 1, 0, 0, 2, 2]; +aircraft.Specs.Propulsion.PropArch.OperDwn = @() splits; +aircraft.Specs.Propulsion.PropArch.EtaDwn = ones(9); +aircraft.Specs.Power.LamDwn.SLS = []; +end + +function aircraft = makeAircraft(propellerEngine) +% +% aircraft = makeAircraft(propellerEngine) +% written by Triet Ho +% last updated: 16 sep 2026 +% +% Build a fuel-to-engine-to-propeller graph with one sink. Each vector +% entry identifies the engine that drives the corresponding propeller. +% +% INPUTS: +% propellerEngine - one-based engine index for each propeller. +% size/type/units: 1-by-n / integer / [] +% +% OUTPUTS: +% aircraft - aircraft structure with the test graph and powers. +% size/type/units: 1-by-1 / struct / [] +% + +engineCount = max(propellerEngine); +propellerCount = numel(propellerEngine); +componentCount = 1 + engineCount + propellerCount + 1; +architecture = zeros(componentCount); +for engine = 1:engineCount + architecture(1, 1 + engine) = 1; +end +for propeller = 1:propellerCount + propellerIndex = 1 + engineCount + propeller; + architecture(1 + propellerEngine(propeller), propellerIndex) = 1; + architecture(propellerIndex, componentCount) = 1; +end +aircraft.Specs.Propulsion.PropArch.Arch = architecture; +aircraft.Specs.Propulsion.PropArch.SrcType = 1; +aircraft.Specs.Propulsion.PropArch.TrnType = [ones(1, engineCount), 2 * ones(1, propellerCount)]; +aircraft.Specs.Propulsion.PropArch.OperDwn = @() ones(componentCount); +aircraft.Specs.Propulsion.SLSPower = [100 * ones(1, engineCount), 40 * ones(1, propellerCount)]; +aircraft.Specs.Power.LamDwn.SLS = []; +end + +function connections = propellerConnections(aircraft, engine) +% +% connections = propellerConnections(aircraft, engine) +% written by Triet Ho +% last updated: 16 sep 2026 +% +% Return sorted propeller indices, accepting both scalar and cell-array +% representations of the engine's connection list. +% +% INPUTS: +% aircraft - processed aircraft structure with propeller links. +% size/type/units: 1-by-1 / struct / [] +% engine - one-based transmitter index of the engine. +% size/type/units: 1-by-1 / integer / [] +% +% OUTPUTS: +% connections - sorted, one-based propeller component indices. +% size/type/units: 1-by-n / integer / [] +% + +entry = aircraft.Specs.Propulsion.PropArch.WhichProp(engine); +if iscell(entry) + entry = entry{1}; +end +connections = sort(entry(:)'); +end + +function supplement = parallelSupplement(types, requiredPower, fanEfficiency) +% +% supplement = parallelSupplement(types, requiredPower, fanEfficiency) +% written by Triet Ho +% last updated: 16 sep 2026 +% +% Make one engine and each motor feed a common target. The input powers +% determine matching target demand and downstream fractions. +% +% INPUTS: +% types - transmitter type for each graph component. +% size/type/units: 1-by-n / integer / [] +% requiredPower - required power of each component. +% size/type/units: 1-by-n / double / [W] +% fanEfficiency - fan efficiency applied at fan targets. +% size/type/units: 1-by-1 / double / [] +% +% OUTPUTS: +% supplement - motor power credited to the engine. +% size/type/units: 1-by-1 / double / [W] +% + +componentCount = numel(types); +architecture = zeros(componentCount); +architecture(1:end-1, end) = 1; +targetPower = sum(requiredPower(1:end-1)); +requiredPower(end) = targetPower; +splits = zeros(componentCount); +splits(end, 1:end-1) = requiredPower(1:end-1) ./ targetPower; +efficiencies = ones(componentCount); +output = PropulsionPkg.PowerSupplementCheck(requiredPower, architecture, splits, efficiencies, types, fanEfficiency); +supplement = output(1); +end + +function supplement = consistentParallelSupplement(edgeEfficiency, fanEfficiency) +% +% supplement = consistentParallelSupplement(edgeEfficiency, fanEfficiency) +% written by Triet Ho +% last updated: 16 sep 2026 +% +% Return motor assistance for a power-consistent, equal-share fan. The fan +% requires 100 W; each parent delivers 50 W after its edge loss. +% +% INPUTS: +% edgeEfficiency - efficiency of both parent-to-fan connections. +% size/type/units: 1-by-1 / double / [] +% fanEfficiency - fan efficiency used for motor credit. +% size/type/units: 1-by-1 / double / [] +% +% OUTPUTS: +% supplement - motor assistance credited to the engine. +% size/type/units: 1-by-1 / double / [W] +% + +architecture = zeros(3); +architecture(1:2, 3) = 1; +splits = zeros(3); +splits(3, 1:2) = 0.5; +efficiencies = ones(3); +efficiencies(3, 1:2) = edgeEfficiency; +parentPower = 50 / edgeEfficiency; +output = PropulsionPkg.PowerSupplementCheck( ... + [parentPower, parentPower, 100], architecture, splits, efficiencies, ... + [1, 0, 2], fanEfficiency); +supplement = output(1); +end + +function pass = CheckInterleavedEngineInletAreas() +% +% pass = CheckInterleavedEngineInletAreas() +% written by Triet Ho +% last updated: 16 sep 2026 +% +% Reorder the two-engine sizing graph so motors precede engines, then check +% that each sized engine assigns its inlet area to its connected fan. +% +% INPUTS: +% none +% +% OUTPUTS: +% pass - true if both fan slots receive finite inlet areas. +% size/type/units: 1-by-1 / logical / [] +% + +aircraft = makeTwoEngineSizingAircraft(0.2, 0.4); +propArch = aircraft.Specs.Propulsion.PropArch; +order = [1, 2, 5, 3, 6, 4, 7, 8, 9]; +splits = propArch.OperDwn(); +propArch.Arch = propArch.Arch(order, order); +propArch.OperDwn = @() splits(order, order); +propArch.TrnType = propArch.TrnType(order(3:8) - 2); +aircraft.Specs.Propulsion.PropArch = propArch; +aircraft.Specs.Propulsion.InletArea = NaN(1, 6); +aircraft = PropulsionPkg.PropulsionSizing(aircraft); +areas = aircraft.Specs.Propulsion.InletArea; +pass = numel(areas) == 6 && all(isfinite(areas(5:6))) && ... + all(isnan(areas([1, 3]))); +if (~pass) + fprintf(1, "Interleaved engine inlet areas: %s\n", mat2str(areas)); +end +end + +function pass = CheckSeparateParallelTargets() +% +% pass = CheckSeparateParallelTargets() +% written by Triet Ho +% last updated: 16 sep 2026 +% +% Check full-throttle recomputation when one engine has a different motor +% helper at each of two fans. Expected splits use each fan's own output. +% +% INPUTS: +% none +% +% OUTPUTS: +% pass - true if both fan-specific splits are recomputed correctly. +% size/type/units: 1-by-1 / logical / [] +% + +aircraft.Specs.Propulsion.PropArch.Arch = zeros(5); +aircraft.Specs.Propulsion.PropArch.Arch(1, 4:5) = 1; +aircraft.Specs.Propulsion.PropArch.Arch(2, 4) = 1; +aircraft.Specs.Propulsion.PropArch.Arch(3, 5) = 1; +aircraft.Specs.Propulsion.PropArch.SrcType = []; +aircraft.Specs.Propulsion.PropArch.TrnType = [1, 0, 0, 2, 2]; +aircraft.Specs.Propulsion.PropArch.OperUps = @separateParallelUps; +aircraft.Specs.Propulsion.PropArch.OperDwn = @separateParallelSplits; +aircraft.Specs.Propulsion.PropArch.EtaUps = ones(5); +aircraft.Specs.Power.LamDwn.SLS = [0.5, 0.5]; +aircraft.Mission.History.SI.Power.Pav = [100, 20, 30, 70, 80]; +aircraft.Mission.History.SI.Power.LamUps = 1; +aircraft.Mission.History.SI.Power.LamDwn = [0.5, 0.5]; +aircraft = PropulsionPkg.PropArchConnections(aircraft); +try + output = PropulsionPkg.RecomputeSplits(aircraft, 1, 1); + actual = output.Mission.History.SI.Power.LamDwn; + expected = [20 / 70, 30 / 80]; + pass = isequal(size(actual), size(expected)) && ... + all(abs(actual - expected) < 1e-10); +catch error + fprintf(1, "Separate parallel targets: %s\n", error.message); + pass = false; +end +end + +function splits = separateParallelSplits(firstMotorShare, secondMotorShare) +% +% splits = separateParallelSplits(firstMotorShare, secondMotorShare) +% written by Triet Ho +% last updated: 16 sep 2026 +% +% Evaluate independent motor shares for two different fans. +% +% INPUTS: +% firstMotorShare - motor fraction at the first fan. +% size/type/units: 1-by-1 / double / [] +% secondMotorShare - motor fraction at the second fan. +% size/type/units: 1-by-1 / double / [] +% +% OUTPUTS: +% splits - downstream operational matrix. +% size/type/units: 5-by-5 / double / [] +% + +splits = zeros(5); +splits(4, 1:2) = [1 - firstMotorShare, firstMotorShare]; +splits(5, [1, 3]) = [1 - secondMotorShare, secondMotorShare]; +end + +function splits = separateParallelUps(~) +% +% splits = separateParallelUps(~) +% written by Triet Ho +% last updated: 16 sep 2026 +% +% Route half the engine output to each fan and each motor to its own fan. +% +% INPUTS: +% unused - full-throttle upstream split marker. +% size/type/units: 1-by-1 / double / [] +% +% OUTPUTS: +% splits - upstream operational matrix. +% size/type/units: 5-by-5 / double / [] +% + +splits = zeros(5); +splits(1, 4:5) = 0.5; +splits(2, 4) = 1; +splits(3, 5) = 1; +end + +function pass = CheckMixedDepthEngineConnections() +% +% pass = CheckMixedDepthEngineConnections() +% written by Triet Ho +% last updated: 16 sep 2026 +% +% Check a propeller receiving half its power from a direct engine and half +% from an engine reached through a cable. Both engines must be accounted for. +% +% INPUTS: +% none +% +% OUTPUTS: +% pass - true if both engines receive their expected hybrid coefficient. +% size/type/units: 1-by-1 / logical / [] +% + +aircraft.Specs.Propulsion.PropArch.Arch = zeros(6); +aircraft.Specs.Propulsion.PropArch.Arch(1, 2:3) = 1; +aircraft.Specs.Propulsion.PropArch.Arch(2, 5) = 1; +aircraft.Specs.Propulsion.PropArch.Arch(3, 4) = 1; +aircraft.Specs.Propulsion.PropArch.Arch(4, 5) = 1; +aircraft.Specs.Propulsion.PropArch.Arch(5, 6) = 1; +aircraft.Specs.Propulsion.PropArch.SrcType = 1; +aircraft.Specs.Propulsion.PropArch.TrnType = [1, 1, 4, 2]; +splits = zeros(6); +splits(2:3, 1) = 1; +splits(4, 3) = 1; +splits(5, [2, 4]) = 0.5; +splits(6, 5) = 1; +aircraft.Specs.Propulsion.PropArch.OperDwn = @() splits; +aircraft.Specs.Propulsion.SLSPower = [100, 100, 100, 100]; +aircraft.Specs.Power.LamDwn.SLS = []; +aircraft = PropulsionPkg.ProcessPropArch(aircraft); +actual = aircraft.Specs.Propulsion.Engine.HEcoeff(1:2); +pass = all(abs(actual - [1.5, 1.5]) < 1e-10); +if (~pass) + fprintf(1, "Mixed-depth engine coefficients: %s\n", mat2str(actual)); +end +end + +function pass = CheckManyEngineConnections(engineCount) +% +% pass = CheckManyEngineConnections(engineCount) +% written by Triet Ho +% last updated: 16 sep 2026 +% +% Check that each engine retains its own direct propeller connection. +% +% INPUTS: +% engineCount - number of TSEs and propellers in the graph. +% size/type/units: 1-by-1 / integer / [] +% +% OUTPUTS: +% pass - true if every engine has its expected propeller. +% size/type/units: 1-by-1 / logical / [] +% + +aircraft = PropulsionPkg.ProcessPropArch(makeAircraft(1:engineCount)); +pass = true; +for engine = 1:engineCount + pass = pass && isequal(propellerConnections(aircraft, engine), ... + 1 + engineCount + engine); +end +end + +function pass = CheckManyEngineMotorCredit(engineCount) +% +% pass = CheckManyEngineMotorCredit(engineCount) +% written by Triet Ho +% last updated: 16 sep 2026 +% +% Check proportional motor credit when many TSEs share one fan. The +% downstream shares and parent powers describe the same edge flows. +% +% INPUTS: +% engineCount - number of TSEs feeding the fan. +% size/type/units: 1-by-1 / integer / [] +% +% OUTPUTS: +% pass - true if each TSE receives its expected motor credit. +% size/type/units: 1-by-1 / logical / [] +% + +fan = engineCount + 2; +architecture = zeros(fan); +architecture(1:fan-1, fan) = 1; +engineShares = 0.6 * (1:engineCount) / sum(1:engineCount); +splits = zeros(fan); +splits(fan, :) = [engineShares, 0.4, 0]; +efficiencies = ones(fan); +efficiencies(fan, 1:fan-1) = 0.8; +parentPower = 100 * splits(fan, 1:fan-1) / 0.8; +actual = PropulsionPkg.PowerSupplementCheck( ... + [parentPower, 100], architecture, splits, efficiencies, ... + [ones(1, engineCount), 0, 2], 0.99); +expected = 40 * 0.99 * engineShares / sum(engineShares); +pass = numel(actual) == fan && ... + all(abs(actual(1:engineCount) - expected) < 1e-10); +end + +function pass = CheckManyInterleavedEngineInletAreas(engineCount) +% +% pass = CheckManyInterleavedEngineInletAreas(engineCount) +% written by Triet Ho +% last updated: 16 sep 2026 +% +% Size independently supplemented engines interleaved with motors. Each +% engine must assign its inlet area and motor load to its own fan. +% +% INPUTS: +% engineCount - number of TSE, motor, and fan triples. +% size/type/units: 1-by-1 / integer / [] +% +% OUTPUTS: +% pass - true when all fan areas and supplements are assigned. +% size/type/units: 1-by-1 / logical / [] +% + +aircraft = makeTwoEngineSizingAircraft(0.2, 0.4); +ncomp = 3 * engineCount + 3; +sink = ncomp; +fans = 2 * engineCount + 3 : 3 * engineCount + 2; +motors = 3:2:2 * engineCount + 1; +engines = 4:2:2 * engineCount + 2; +architecture = zeros(ncomp); +architecture(1, engines) = 1; +architecture(2, motors) = 1; +splits = ones(ncomp); +for engine = 1:engineCount + fan = fans(engine); + architecture([motors(engine), engines(engine)], fan) = 1; + motorShare = 0.1 + 0.7 * engine / (engineCount + 1); + splits(fan, [motors(engine), engines(engine)]) = ... + [motorShare, 1 - motorShare]; +end +architecture(fans, sink) = 1; +splits(sink, fans) = 1 / engineCount; +aircraft.Specs.Propulsion.PropArch.Arch = architecture; +aircraft.Specs.Propulsion.PropArch.TrnType = ... + [repmat([0, 1], 1, engineCount), 2 * ones(1, engineCount)]; +aircraft.Specs.Propulsion.PropArch.OperDwn = @() splits; +aircraft.Specs.Propulsion.PropArch.EtaDwn = ones(ncomp); +aircraft.Specs.Propulsion.InletArea = NaN(1, 3 * engineCount); +aircraft = PropulsionPkg.PropulsionSizing(aircraft); +areas = aircraft.Specs.Propulsion.InletArea; +supplements = aircraft.Specs.Propulsion.PowerSupp(2:2:2 * engineCount); +actual = arrayfun(@(engine) engine.FanSysObject.ElecWork, ... + aircraft.Specs.Propulsion.SizedEngine); +pass = numel(areas) == 3 * engineCount && ... + all(isnan(areas(1:2 * engineCount))) && ... + all(isfinite(areas(2 * engineCount + 1:end))) && ... + numel(actual) == engineCount && all(isfinite(actual)) && ... + all(abs(actual(:)' - supplements) <= 1e-8 * max(1, abs(supplements))); +if (~pass) + fprintf(1, "%d interleaved engine areas: %s\n", engineCount, mat2str(areas)); +end +end + +function pass = CheckManyMixedDepthEngineConnections(engineCount) +% +% pass = CheckManyMixedDepthEngineConnections(engineCount) +% written by Triet Ho +% last updated: 16 sep 2026 +% +% One TSE connects directly to a propeller; the others connect through +% separate cables. Each should deliver an equal share of the power. +% +% INPUTS: +% engineCount - number of TSEs feeding the propeller. +% size/type/units: 1-by-1 / integer / [] +% +% OUTPUTS: +% pass - true when every hybrid coefficient is correct. +% size/type/units: 1-by-1 / logical / [] +% + +propeller = 2 * engineCount + 1; +sink = propeller + 1; +engines = 2:engineCount + 1; +cables = engineCount + 2:2 * engineCount; +architecture = zeros(sink); +architecture(1, engines) = 1; +architecture(engines(1), propeller) = 1; +splits = zeros(sink); +splits(engines, 1) = 1; +for engine = 2:engineCount + cable = cables(engine - 1); + architecture(engines(engine), cable) = 1; + architecture(cable, propeller) = 1; + splits(cable, engines(engine)) = 1; +end +architecture(propeller, sink) = 1; +splits(propeller, [engines(1), cables]) = 1 / engineCount; +splits(sink, propeller) = 1; +aircraft.Specs.Propulsion.PropArch.Arch = architecture; +aircraft.Specs.Propulsion.PropArch.SrcType = 1; +aircraft.Specs.Propulsion.PropArch.TrnType = ... + [ones(1, engineCount), 4 * ones(1, engineCount - 1), 2]; +aircraft.Specs.Propulsion.PropArch.OperDwn = @() splits; +aircraft.Specs.Propulsion.SLSPower = [100 * ones(1, 2 * engineCount - 1), 90]; +aircraft.Specs.Power.LamDwn.SLS = []; +aircraft = PropulsionPkg.ProcessPropArch(aircraft); +actual = aircraft.Specs.Propulsion.Engine.HEcoeff(1:engineCount); +pass = all(abs(actual - (2 - 1 / engineCount)) < 1e-10); +if (~pass) + fprintf(1, "%d mixed-depth engine coefficients: %s\n", ... + engineCount, mat2str(actual)); +end +end diff --git a/TestFAST.m b/TestFAST.m index 5424d03..7d9a8f5 100644 --- a/TestFAST.m +++ b/TestFAST.m @@ -30,6 +30,7 @@ % run test cases and print results PropulsionPkg.TestCreatePropArch(); PropulsionPkg.TestPowerAvailable(); +PropulsionPkg.TestEngineMotorRegressions(); %% TEST UNIT CONVERSION PACKAGE %% @@ -51,4 +52,4 @@ fprintf(1, "\nAll FAST tests completed! Check above to see if any have failed.\n"); -end \ No newline at end of file +end