fix: improved robustness for multi-cooling system configurations including spray cooling

Author: agonzale

GenerateSimulinkThermalModel.mlx

@@ -66,4 +66,4 @@ The license is available in the License file within this repository.
 ## Community Support
 [MATLAB Central](https://www.mathworks.com/matlabcentral)
 
-Copyright 2022-2024 The MathWorks, Inc.
+Copyright 2022-2025 The MathWorks, Inc.

README.md

@@ -887,7 +887,7 @@ function updateLabLossTables(obj)
 
             obj.turnOffLossDependenceWithTemperatureOrSpeed();
 
-            % 1) COMPUTE ROM DATA
+            % 1) COMPUTE ROM DATA  ---------------------------------------
 
             % Enable the cooling systems and get the flowrate and inlet
             % temperature property names.
@@ -912,6 +912,7 @@ function updateLabLossTables(obj)
                 end
             end
             obj.updateModel();
+            AdjacencyMat = obj.AdjacencyMat;
 
             % Get ambient temperature and initial node temperatures
             Tambient_degC =  obj.Tambient_degC;
@@ -942,18 +943,47 @@ function updateLabLossTables(obj)
             samplingGridStruct = getSamplingGridFromBkpts(BkptsStruct);
             StateSpaceND.SamplingGrid = samplingGridStruct;
 
-            % Get coolant data
-            AdjacencyMat = obj.AdjacencyMat;
+            % Build inlet/outlet index rows that match coolingSystemsEnabled one-to-one
+            numCoolSys   = numel(coolingSystemsEnabled);
+            inletCells   = cell(1,numCoolSys);
+            outletCells  = cell(1,numCoolSys);           
+            for kSys = 1:numCoolSys
+                csName = coolingSystemsEnabled{kSys};
+                % identify which rows in the object belong to csName
+                if any(startsWith(obj.MultiNozzleSprayCoolingSystems, string(csName))) % is multi-nozzle spray cooling (front and rear)
+                    % Multi-nozzle aggregate: collect all spray rows
+                    rowMask = ismember(obj.CoolingSystemNamesAndMcadIdxes(:,1), ...
+                                       obj.MultiNozzleSprayCoolingSystems);
+                else
+                    % Normal case: exact match
+                    rowMask = strcmp(obj.CoolingSystemNamesAndMcadIdxes(:,1), csName);
+                end
+                rows = find(rowMask);
+                if isempty(rows)
+                    error(['Cooling system "%s" is not enabled in Motor-CAD (or its ' ...
+                           'name does not match SupportedCoolingSystemsMcadNames).'], csName);
+                end
+                % concatenate indices from all matching rows
+                inIdx  = [];
+                outIdx = [];
+                for r = rows(:)'
+                    inIdx  = [inIdx  obj.InletArrayIdxs{r}];           %#ok<AGROW>
+                    tmpOut = obj.OutletArrayIdxs{r};
+                    if isempty(tmpOut)          % per-element fallback
+                        tmpOut = obj.InletArrayIdxs{r};
+                    end
+                    outIdx = [outIdx tmpOut];   %#ok<AGROW>
+                end
+                inletCells{kSys}  = inIdx;
+                outletCells{kSys} = outIdx;
+            end
+            obj.InletArrayIdxs = inletCells;
+            obj.OutletArrayIdxs = outletCells;
+
             % Simulink does not accept cell arrays as parameters, need to
             % transform them to arrays:
             InletCoolIdxs = adaptInletOutletArrays(obj.InletArrayIdxs); 
             OutletCoolIdxs = adaptInletOutletArrays(obj.OutletArrayIdxs);
-            % Some cooling systems (multinozzle spray) have no outlet node in Motor-CAD.
-            if isempty(OutletCoolIdxs)
-                % FixMe: temporary workaround, assume outlet=inlet (no heat transfer)
-                OutletCoolIdxs = InletCoolIdxs;
-                % TO DO: estimate outlet coolant temperature from the heat transfer
-            end
 
             % Calculate state-space model at each breakpoint -------
             idxsCell = getNestedForLoopIdxs(sizeMatND(1:end-2));
@@ -1036,7 +1066,8 @@ function updateLabLossTables(obj)
             xCoolantArrayIdxs = obj.CoolantArrayIdxs;            
             xCapMat = obj.CapMat;
 
-            % 2) GENERATE SIMULINK MODEL
+            % 2) GENERATE SIMULINK MODEL ---------------------------------
+            
             disp("Creating Simulink model...")
             % create and open the model
             open_system(new_system(modelName));

ValidateSimulinkThermalModel.mlx

@@ -16,20 +16,23 @@
     % is the SROTM node temperatures and Tmcad is the Motor-CAD node
     % temperature, in the transient simulation.
 
-    % Copyright 2022 The MathWorks, Inc.
+    % Copyright 2022-2025 The MathWorks, Inc.
 
     mdlWks = get_param(modelName,'ModelWorkspace');
     StateSpaceND = mdlWks.getVariable('StateSpaceND');
     samplingGrid = StateSpaceND.SamplingGrid;
 
-    if length(bkptIdxComb)==3 % one cooling system
+    % Obtain Simulink results ------------------------------------------
+    if length(bkptIdxComb)==3 % one cooling system - e5_IM_HWJ case
         wBkpts = squeeze(samplingGrid.w(:,1,1))';
         fr1Bkpts = squeeze(samplingGrid.fr1(1,:,1));
         Tin1Bkpts = squeeze(samplingGrid.Tin1(1,1,:))';
         rpmVal = wBkpts(bkptIdxComb(1));
         fr1Val = fr1Bkpts(bkptIdxComb(2));
-        Tin1Val = Tin1Bkpts(bkptIdxComb(3));        
-    elseif length(bkptIdxComb)==5 % two cooling systems
+        Tin1Val = Tin1Bkpts(bkptIdxComb(3));
+        set_param(strcat(modelName, '/HousingWaterJacket_Flowrate_lpm'), 'Value', num2str(fr1Val));
+        set_param(strcat(modelName, '/HousingWaterJacket_InletTemp_degC'), 'Value', num2str(Tin1Val));
+    elseif length(bkptIdxComb)==5 % two cooling systems - e8_IPMSM_HWJandVent case
         wBkpts = squeeze(samplingGrid.w(:,1,1,1,1))';
         fr1Bkpts = squeeze(samplingGrid.fr1(1,:,1,1,1));
         fr2Bkpts = squeeze(samplingGrid.fr2(1,1,:,1,1))';
@@ -40,46 +43,63 @@
         fr2Val = fr2Bkpts(bkptIdxComb(3));
         Tin1Val = Tin1Bkpts(bkptIdxComb(4));
         Tin2Val = Tin2Bkpts(bkptIdxComb(5));
-    end
-
-    set_param(strcat(modelName, '/TorqueNm'), 'Value', num2str(torqueVal));
-    set_param(strcat(modelName, '/SpeedRPM'), 'Value', num2str(rpmVal));
-    % set HWJ for e8 model
-    try
         set_param(strcat(modelName, '/HousingWaterJacket_Flowrate_lpm'), 'Value', num2str(fr1Val));
         set_param(strcat(modelName, '/HousingWaterJacket_InletTemp_degC'), 'Value', num2str(Tin1Val));
-    catch
-        % just continue
-    end
-    % set spray cooling for e8 model
-    try
+        set_param(strcat(modelName, '/Ventilated_Flowrate_lpm'), 'Value', num2str(fr2Val));
+        set_param(strcat(modelName, '/Ventilated_InletTemp_degC'), 'Value', num2str(Tin2Val));
+    elseif length(bkptIdxComb)==7 % three cooling systems - e8_IPMSM_SprayMultiNozzle case
+        wBkpts = squeeze(samplingGrid.w(:,1,1,1,1,1,1))';
+        fr1Bkpts = squeeze(samplingGrid.fr1(1,:,1,1,1,1,1)); 
+        Tin1Bkpts = squeeze(samplingGrid.Tin1(1,1,:,1,1,1,1))'; 
+        fr2Bkpts = squeeze(samplingGrid.fr2(1,1,1,:,1,1,1))'; 
+        Tin2Bkpts = squeeze(samplingGrid.Tin2(1,1,1,1,:,1,1))'; 
+        fr3Bkpts = squeeze(samplingGrid.fr3(1,1,1,1,1,:,1))'; 
+        Tin3Bkpts = squeeze(samplingGrid.Tin3(1,1,1,1,1,1,:))'; 
+        rpmVal = wBkpts(bkptIdxComb(1));
+        fr1Val = fr1Bkpts(bkptIdxComb(2));       
+        fr2Val = fr2Bkpts(bkptIdxComb(3));
+        fr3Val = fr3Bkpts(bkptIdxComb(4));
+        Tin1Val = Tin1Bkpts(bkptIdxComb(5));
+        Tin2Val = Tin2Bkpts(bkptIdxComb(6));
+        Tin3Val = Tin3Bkpts(bkptIdxComb(7));
         set_param(strcat(modelName, '/Spray_RadialHousing_Flowrate_lpm'), 'Value', num2str(fr1Val));
         set_param(strcat(modelName, '/Spray_RadialHousing_InletTemp_degC'), 'Value', num2str(Tin1Val));
-    catch
-        % just continue
+        set_param(strcat(modelName, '/HousingWaterJacket_Flowrate_lpm'), 'Value', num2str(fr2Val));
+        set_param(strcat(modelName, '/HousingWaterJacket_InletTemp_degC'), 'Value', num2str(Tin2Val));
+        set_param(strcat(modelName, '/Spray_RadialRotor_Flowrate_lpm'), 'Value', num2str(fr3Val));
+        set_param(strcat(modelName, '/Spray_RadialRotor_InletTemp_degC'), 'Value', num2str(Tin3Val));
     end
 
-    if length(bkptIdxComb)==5
-        set_param(strcat(modelName, '/Ventilated_Flowrate_lpm'), 'Value', num2str(fr2Val));
-        set_param(strcat(modelName, '/Ventilated_InletTemp_degC'), 'Value', num2str(Tin2Val));
-    end
+    set_param(strcat(modelName, '/TorqueNm'), 'Value', num2str(torqueVal));
+    set_param(strcat(modelName, '/SpeedRPM'), 'Value', num2str(rpmVal));
 
     stopTime = 1000; % s
     out = sim(modelName, 'StopTime', num2str(stopTime));
     TnodesSeries = out.yout{1}.Values;
-
     TnodesInit = mdlWks.getVariable('TnodesInit');
+
+    % Obtain Motor-CAD baseline ------------------------------------------
     mcadIntf.Tambient_degC = TnodesInit(1);
     mcadIntf.Shaft_Speed_RPM = rpmVal;
-    mcadIntf.HousingWaterJacket_FlowRate_m3ps = fr1Val/60/1000; % lpm to m3ps
-    mcadIntf.HousingWaterJacket_InletTemperature_degC = Tin1Val;
-    mcadIntf.Spray_RadialHousing_FlowRate_m3ps = fr1Val/60/1000; % lpm to m3ps
-    mcadIntf.Spray_RadialHousing_InletTemperature_F_degC = Tin1Val;
-    mcadIntf.Spray_RadialHousing_InletTemperature_R_degC = Tin1Val;
-    if length(bkptIdxComb)==5
+    if length(bkptIdxComb)==3 % one cooling system - e5_IM_HWJ case
+        mcadIntf.HousingWaterJacket_FlowRate_m3ps = fr1Val/60/1000; % lpm to m3ps
+        mcadIntf.HousingWaterJacket_InletTemperature_degC = Tin1Val;
+    elseif length(bkptIdxComb)==5 % two cooling systems - e8_IPMSM_HWJandVent case
+        mcadIntf.HousingWaterJacket_FlowRate_m3ps = fr1Val/60/1000; % lpm to m3ps
+        mcadIntf.HousingWaterJacket_InletTemperature_degC = Tin1Val;
         mcadIntf.Ventilated_FlowRate_m3ps = fr2Val/60/1000; % lpm to m3ps
         mcadIntf.Ventilated_InletTemperature_degC = Tin2Val;
+    elseif length(bkptIdxComb)==7 % three cooling systems - e8_IPMSM_SprayMultiNozzle case
+        mcadIntf.Spray_RadialHousing_FlowRate_m3ps = fr1Val/60/1000; % lpm to m3ps
+        mcadIntf.Spray_RadialHousing_InletTemperature_F_degC = Tin1Val;
+        mcadIntf.Spray_RadialHousing_InletTemperature_R_degC = Tin1Val;
+        mcadIntf.HousingWaterJacket_FlowRate_m3ps = fr2Val/60/1000; % lpm to m3ps
+        mcadIntf.HousingWaterJacket_InletTemperature_degC = Tin2Val;
+        mcadIntf.Spray_RadialRotor_FlowRate_m3ps = fr3Val/60/1000; % lpm to m3ps
+        mcadIntf.Spray_RadialRotor_InletTemperature_F_degC = Tin3Val;
+        mcadIntf.Spray_RadialRotor_InletTemperature_R_degC = Tin3Val;
     end
+
     mcadIntf.EnableStatorTempCoeffRes = 1;
     mcadIntf.EnableRotorTempCoeffRes = 1;
     mcadIntf.updateModel();    
@@ -88,6 +108,7 @@
     allMcadIdxs = [mcadIntf.NodeNamesAndMcadIdx{:,2}];
     [tVecMcad, TnodesMcad] = mcadIntf.getTransientTemperatureForNodeMcadIdxs(allMcadIdxs);
 
+    % Plot results ------------------------------------------------------
     TnodesMcadSeries = timeseries(TnodesMcad', tVecMcad);
     [TnodesMcadSeries,TnodesSeries] = synchronize(TnodesMcadSeries,TnodesSeries,'Union');
     TnodesMcadSeries.Name = 'Motor-CAD';