Trying to solve a system of heat transfer equations more efficiently

In one operating mode of my code, there is a solid zone annulus and a coolant zone (think a rod with a hole drilled through it), and it is split up into a 1-D grid with an axial step size of 0.01 meters, or 1 centimeter nodes.

There are a huge number of parameters going into this function, as it evaluates power generation according to the initial time step's thermal state of all the current zone's solid material and coolant, as well as other zone thermal parameters, but this is a simplified mode of the code that only accounts for two regions at this time.

FM is fuel material, FC is fuel coolant,

q_dot_FM is the power input to the solid material and is evaluated once against thermal parameters and other items for the OLD time and extrapolated into an axial profile using a power profile (W/m^3) to be a vector of equal length to the grid.

Cp and K are thermally sensitive for each material, and are computed for the OLD time step's corresponding temperature value for that material

the uknowns in this problem are the convective heat transfer into the coolant (q_conv4), new solid fuel temp (T1_FM), new enthalpy of the coolant (enth_FC), and the conductive heat transfer into the downstream node (q_cond6)

the code iterates through the entire 1-D grid with the same dx defined above, and vpasolve is used to solve for those four unknowns

Then, as a way to circumvent me not being able to set this up to solve in matrix form, i have to redefine the initial temperature of the downstream solid material temperature (T0_FM(i+1)) by using the specific heat capacity and that q_cond6 term,

and likewise, I have to redefine the coolant thermal properties by using the enthalpy found from vpasolve to get the new coolant temperature to serve as the inlet to the next node so i can use it to find the new heat transfer coefficient from my core_thermal_hydraulics function

Is there any way to reorient this code to make it run against an array of initial temperature, conductivity, specific heat capacity, power input, heat transfer coefficient (although that one is more tricky because it is evaluated against the coolant conditions as it changes through the grid)? Running long time step spans is causing code runs to take hours

BOL_distance = axial_step_size*i*100; %cm
q_dot_FM(i) = (Fuel_Power_Profile_Function(1)*BOL_distance^2 + Fuel_Power_Profile_Function(2)*BOL_distance + Fuel_Power_Profile_Function(3))*q_dot_FM_temp;   

total_Q = (q_dot_FM(i) * Power_deposition_fuel * FM_vol);

cp_FM = (FM_CP_C0 + FM_CP_C1*(T0_FM(i)/1000) + FM_CP_C2*(T0_FM(i)/1000)^2 + FM_CP_C3/((T0_FM(i)/1000)^2))*1000;          % to transform J/gK to J/kgK            
k_FM   = FM_k_C0 + FM_k_C1*(T0_FM(i)/1000) + FM_k_C2*(T0_FM(i)/1000)^2 + FM_k_C3*(T0_FM(i)/1000)^2;            % W/mK

[HTC_FC,P0_new, enth0_FC] = core_thermal_hydraulics(m_dot_annulus,FC_d_hyd,T_FC(i),P_FC(i),axial_step_size,fluid, total_Q,FC_roughness,cp_FM);

%             syms q_conv4 T1_FM enth_FC  q_cond6
eqn1 = (q_dot_FM(i) * Power_deposition_fuel * FM_vol) - q_conv4 == Fuel_material_density * FM_vol * cp_FM * (T1_FM-T0_FM(i)) / dt;
eqn5 = q_conv4 == HTC_FC* pi * (2*FC_r_o) * axial_step_size * (T1_FM - T_FC(i));
eqn6 = q_conv4 == m_dot_annulus*(enth_FC - enth0_FC);
eqn10 = q_cond6 == k_FM * pi*((2*FM_r_o)^2 - (2*FM_r_i)^2) *(T1_FM - T0_FM(i+1)) / axial_step_size;

sol = vpasolve(eqn1,eqn5,eqn6,eqn10);

enth_FC_outlet = sol.enth_FC;
T1_FM_outlet(i) = sol.T1_FM;
q_cond6_outlet = sol.q_cond6;

%multiply watts by time to get joules, divide by mass and divide by time to get K increase and add to old temp
T0_FM(i+1) = T0_FM(i+1) + (q_cond6_outlet*dt)/(FM_vol*Fuel_material_density)/cp_FM;

% evaluate new coolant temperature
T_FC(i) = (enth_FC_outlet - enth0_FC)/Cp0 + T_FC(i);
T_FC(i+1) = T_FC(i);

% assign inlet pressure for next node
P_FC(i+1) = P0_new;

[Sam Chak]

@Anthony,

Can you fix the unrecognized function?

axial_step_size = 1;

BOL_distance = axial_step_size*i*100; %cm
q_dot_FM(i) = (Fuel_Power_Profile_Function(1)*BOL_distance^2 + Fuel_Power_Profile_Function(2)*BOL_distance + Fuel_Power_Profile_Function(3))*q_dot_FM_temp;   

Unrecognized function or variable 'Fuel_Power_Profile_Function'.

[Anthony]

Fuel_Power_Profile_Function = [-4.00E-04 4.50E-02 1.97E-01];

I did not include values for these parameters since there is just so much going on here, and you will likely run into the same problem when you arrive at the core_thermal_hydraulics function. I was hoping there would be a conceptual answer upon visualizing how the code runs and what type of arrays are printed. That whole code excerpt above is just a portion of the function evaluating heat transfer. I can't include much more than this as it is proprietary to my project.

The code below is more representative of a larger portion of the code (rearranged to exclude the other modes of operation). There is an if statement for if the code has reached the final node (last slice of the rod) so it does not evaluate conduction downstream to the next grid node because there is no downstream node. I could not figure out how to circumvent that when trying to set up a solver to take the whole array of q_dot_FM, T0_FM, enth0_FC, etc. Hypothetically, I could get arrays for heat transfer coefficient (one of my prior noted concerns) by feeding it the entire old thermal array with coolant temp, pressure, density. So, all the known parameters would either have known vectors or be scalar values, like FM_vol or Power_deposition_fuel, etc.

for i = 1:(length_Reactor/axial_step_size)
    if annulus_solver_mode == 0 %one fuel region, no insulator, NO CONDUCTION only convection
        if i < (length_Reactor/axial_step_size)
            
            if i == 1
                T_FC(i) = T_FC0;
                P_FC(i) = P_FC0;
            end
            
            Cp0=PropsSI('Cpmass','T',T_FC(i),'P',P_FC(i)*P_atm,fluid);
            D_FC(i) = PropsSI('D','T',T_FC(i),'P',P_FC(i)*P_atm,fluid); % this is the Density for the coolant using temperature and pressure
            
            % q_dot below is the volumetric power generation that is a function of the temperature of the node before heating (either initialized or the old time step's final temperature, only evaluated once)
            if i == 1
                if isempty(previous_reactivity) == 0
                    if previous_reactivity == 99999;
                        q_dot_FM_temp = 0;
                        final_reactivity = [];
                        Reactivity_Array = []; reactivity_contributions = [];
                    else                        
                        [q_dot_FM_temp,final_reactivity,reactivity_contributions] =  fuel_power_generation(general_reactor_properties,neutron_decay_constants,neutron_decay_fractions, reactivity_curves, ...
                            mean(T0_FM), mean(T0_MM), mean(T0_RM), mean(old_D_FC), mean(D_MC), ...
                            control_drum_rotation_array, current_time_step, previous_time_step,previous_reactivity);                 
                    end
                else
                    [q_dot_FM_temp,final_reactivity,reactivity_contributions] =  fuel_power_generation(general_reactor_properties,neutron_decay_constants,neutron_decay_fractions, reactivity_curves, ...
                        mean(T0_FM), mean(T0_MM), mean(T0_RM), mean(old_D_FC), mean(D_MC), ...
                        control_drum_rotation_array, current_time_step, previous_time_step,previous_reactivity);                
                end
            end
            % axial power peaking profile used to generate true array of heating power over reactor length
            BOL_distance = axial_step_size*i*100; %cm
            q_dot_FM(i) = (Fuel_Power_Profile_Function(1)*BOL_distance^2 + Fuel_Power_Profile_Function(2)*BOL_distance + Fuel_Power_Profile_Function(3))*q_dot_FM_temp;   
            
            total_Q = (q_dot_FM(i) * Power_deposition_fuel * FM_vol);
            
            cp_FM = (FM_CP_C0 + FM_CP_C1*(T0_FM(i)/1000) + FM_CP_C2*(T0_FM(i)/1000)^2 + FM_CP_C3/((T0_FM(i)/1000)^2))*1000;          % to transform J/gK to J/kgK            
            k_FM   = FM_k_C0 + FM_k_C1*(T0_FM(i)/1000) + FM_k_C2*(T0_FM(i)/1000)^2 + FM_k_C3*(T0_FM(i)/1000)^2;            % W/mK
            
            [HTC_FC,P0_new, enth0_FC] = core_thermal_hydraulics(m_dot_annulus,FC_d_hyd,T_FC(i),P_FC(i),axial_step_size,fluid, total_Q,FC_roughness,cp_FM);
            
            %             syms q_conv4 T1_FM enth_FC  q_cond6
            eqn1 = (q_dot_FM(i) * Power_deposition_fuel * FM_vol) - q_conv4 == Fuel_material_density * FM_vol * cp_FM * (T1_FM-T0_FM(i)) / dt;
            eqn5 = q_conv4 == HTC_FC* pi * (2*FC_r_o) * axial_step_size * (T1_FM - T_FC(i));
            eqn6 = q_conv4 == m_dot_annulus*(enth_FC - enth0_FC);
            eqn10 = q_cond6 == k_FM * pi*((2*FM_r_o)^2 - (2*FM_r_i)^2) *(T1_FM - T0_FM(i+1)) / axial_step_size;
            
            sol = vpasolve(eqn1,eqn5,eqn6,eqn10);
            
            enth_FC_outlet = sol.enth_FC;
            T1_FM_outlet(i) = sol.T1_FM;
            q_cond6_outlet = sol.q_cond6;
            
            %multiply watts by time to get joules, divide by mass and divide by time to get K increase and add to old temp
            T0_FM(i+1) = T0_FM(i+1) + (q_cond6_outlet*dt)/(FM_vol*Fuel_material_density)/cp_FM;
            
            % evaluate new coolant temperature
            T_FC(i) = (enth_FC_outlet - enth0_FC)/Cp0 + T_FC(i);
            T_FC(i+1) = T_FC(i);
            
            % assign inlet pressure for next node
            P_FC(i+1) = P0_new;
            
        elseif i == (length_Reactor/axial_step_size)
            
            Cp0=PropsSI('Cpmass','T',T_FC(i),'P',P_FC(i)*P_atm,fluid);
            D_FC(i) = PropsSI('D','T',T_FC(i),'P',P_FC(i)*P_atm,fluid); % this is the Density for the coolant using temperature and pressure
            
            BOL_distance = axial_step_size*i*100; %cm
            q_dot_FM(i) = (Fuel_Power_Profile_Function(1)*BOL_distance^2 + Fuel_Power_Profile_Function(2)*BOL_distance + Fuel_Power_Profile_Function(3))*q_dot_FM_temp;   
            
            total_Q = (q_dot_FM(i) * Power_deposition_fuel * FM_vol);
            
            cp_FM = (FM_CP_C0 + FM_CP_C1*(T0_FM(i)/1000) + FM_CP_C2*(T0_FM(i)/1000)^2 + FM_CP_C3/((T0_FM(i)/1000)^2))*1000;            % J/kgK
            
            [HTC_FC,P0_new, enth0_FC] = core_thermal_hydraulics(m_dot_annulus,FC_d_hyd,T_FC(i),P_FC(i),axial_step_size,fluid, total_Q,FC_roughness,cp_FM);
            
            %             syms q_conv4 T1_FM enth_FC 
            eqn1 = (q_dot_FM(i) * Power_deposition_fuel * FM_vol) - q_conv4 == Fuel_material_density * FM_vol * cp_FM * (T1_FM-T0_FM(i)) / dt;
            eqn5 = q_conv4 ==HTC_FC* pi * (2*FC_r_o) * axial_step_size * (T1_FM - T_FC(i));
            eqn6 = q_conv4 == m_dot_annulus*(enth_FC - enth0_FC);
            
            sol = vpasolve(eqn1,eqn5,eqn6);
            
            enth_FC_outlet = sol.enth_FC;
            T1_FM_outlet(i) = sol.T1_FM;
            
            % evaluate new coolant temperature
            T_FC(i) = (enth_FC_outlet - enth0_FC)/Cp0 + T_FC(i);
            
            % assign inlet pressure for final node
            P_FC(i) = P0_new;            
            
            Fuel_thermal_array = double([q_dot_FM', T1_FM_outlet', T_FC', P_FC',D_FC']);
            Reactivity_Array = [final_reactivity];
            Outlet = [m_dot, T_FC(end), P_FC(end)];
            
        end     
        
        
    end
    
end

Answer 1

[Nipun]

Hi Anthony,

I understand that you are trying to optimize the code to run it against an array of initial temperature, conductivity, specific heat capacity, power input, heat transfer coefficient given the system of equations. I assume that you are working with the code that was shared in the comments.

To improve the performance of your code and make it more scalable, you can consider vectorizing the calculations and using matrix operations instead of loops and symbolic solvers. This approach is generally more efficient in MATLAB.

Here's a high-level overview of how you might reorganize your code:

1. Preallocate Arrays: Preallocate arrays to store the results, which can significantly improve performance in MATLAB.
2. Vectorized Power Profile Calculation: Instead of looping to calculate the axial power profile, use vectorized operations. You can create a vector of distances and calculate the power profile for all nodes at once.
3. Use Array Operations for Equations: Replace symbolic solvers with array operations. You can create arrays for all the equations and solve them in one go using array operations.

Here's a simplified example to give you an idea:

% Preallocate arrays
T_FC = zeros(1, length_Reactor/axial_step_size + 1);
P_FC = zeros(1, length_Reactor/axial_step_size + 1);
% ... other arrays

% Vectorized distance calculation
distances = axial_step_size * (1:length(T_FC)) * 100; % Assuming axial_step_size is a scalar

% Vectorized power profile calculation
q_dot_FM = (Fuel_Power_Profile_Function(1) * distances.^2 + Fuel_Power_Profile_Function(2) * distances + Fuel_Power_Profile_Function(3)) * q_dot_FM_temp;

% Vectorized property calculations
cp_FM = (FM_CP_C0 + FM_CP_C1*(T0_FM/1000) + FM_CP_C2*(T0_FM/1000).^2 + FM_CP_C3./((T0_FM/1000).^2)) * 1000;
k_FM = FM_k_C0 + FM_k_C1*(T0_FM/1000) + FM_k_C2*(T0_FM/1000).^2 + FM_k_C3*(T0_FM/1000).^2;

% Vectorized thermal hydraulics
[HTC_FC, P0_new, enth0_FC] = core_thermal_hydraulics_vectorized(m_dot_annulus, FC_d_hyd, T_FC, P_FC, axial_step_size, fluid, q_dot_FM, FC_roughness, cp_FM);

% Vectorized equations
eqn1 = q_dot_FM * Power_deposition_fuel * FM_vol - q_conv4 == Fuel_material_density * FM_vol * cp_FM .* (T1_FM - T0_FM) / dt;
eqn5 = q_conv4 == HTC_FC .* pi .* (2*FC_r_o) .* axial_step_size .* (T1_FM - T_FC);
eqn6 = q_conv4 == m_dot_annulus .* (enth_FC - enth0_FC);

% Solve the system of equations
sol = solve(eqn1, eqn5, eqn6);

% Extract results from the solution
enth_FC_outlet = sol.enth_FC;
T1_FM_outlet = sol.T1_FM;
q_cond6_outlet = k_FM .* pi .* ((2*FM_r_o)^2 - (2*FM_r_i)^2) .* (T1_FM - T0_FM(2:end)) / axial_step_size;

% Update temperatures
T0_FM(2:end) = T0_FM(2:end) + (q_cond6_outlet * dt) ./ (FM_vol * Fuel_material_density) ./ cp_FM(2:end);

% Update coolant temperatures
T_FC = (enth_FC_outlet - enth0_FC) ./ Cp0 + T_FC;

% Update pressures
P_FC(2:end) = P0_new;

Adjust it according to the initial conditions and boundary conditions.

Hope this helps.

Regards,

Nipun