The power conversion system is responsible for the electricity production of the nuclear power plant. The system takes in heat from the nuclear reactor through a heat exchanger. Following the heat exchanger, the fluid passes through a turbine, turning its shaft through a generator to produce electricity. For this subsystem, a Rankine cycle design with a single turbine was chosen. This cycle has 4 key components: first, the fluid flows from the reactor heat exchanger through the turbine. From there, it flows through a condenser where it expels waste heat to the cooling tower, then flows through a pump before going through the reactor heat exchanger, also referred to below as the evaporator.
The overall goal of this subsystem is to maximize the ratio of electrical energy output for a given heat flow rate into the system. In order to do this, the system aims to maximize the work output of the turbine while minimizing the work input from the pump and the amount of waste heat expelled through the condenser.
As stated above, the power conversion system is based on a Rankine cycle. The state variables are determined based on thermodynamic and fluid dynamic equations. In many cases, fluid property libraries are used to determine specific values.
Given that the Rankine cycle is a well-understood thermodynamic cycle, many of the state properties can be determined beforehand. As such, there are only three design variables for this system. They are the temperature at the exit of the evaporator, the temperature drop across the turbine, and the pressure inside the evaporator.
There are many more variables tracked at different stages of the system. There are 4 key states in the model. State 1 is the state between the turbine and the condenser. State 2 is the state between the condenser and the pump. State 3 is between the pump and the evaporator. Finally, state 4 is the state between the evaporator and the turbine. At each state, the pressure, temperature, specific enthalpy, and specific entropy are calculated. These values are used to determine the value at next states, determine work and heat flow rates, and verify that thermodynamic laws are followed.
The objective of this sub-system is to maximize the energy output of electricity per unit of energy that flows into the system as heat. The optimization of this system is done using negative-null formatting, which gives the following equation:
Minimize 1/ηconverter where ηconverter = Welec / Qreactor
The condenser and evaporator heat exchanges are reversible thermodynamic processes. As such, the specific entropy between states 1 and 2, along with states 3 and 4, does not change.
The inefficiencies from the pump and the turbine are converted to heat, with a subsequent entropy generation due to this additional heat.
Friction losses in the system are negligible.
The pressure in the evaporator and condenser does not change.
The temperature in the condenser stays constant at the saturation temperature.
Pump Efficiency (ηpump=0.92)
Turbine Efficiency (ηturb=0.92)
Coolant Type = water
Condenser Output Quality: Q = 0
State values of water and determined using Refprop, a fluid table library with a MATLAB library
Evaporator Output temperature in Kelvin (Th)
Temperature drop across the turbine in Kelvin (Tdrop)
Pressure in the evaporator (Pevap)
Boundary Constraints:
Evap Output Temp (K) 350 < Th < 500
Turbine Temp Drop (K) 20 < Tdrop < 120
Evaporator Pressure (kPa) 100 < Pevap < 10000
Linear Constraints:
Condenser Inlet Temperature > 293K
Non-linear Constraints:
Reactor Heat Flow (Qreactor) = 4500MW
Mass flow rate (Mdot) > 16.7 kg/s
Mass flow rate (Mdot) < 2000 kg/s
Heat flow to the condenser (Qcool)
Electrical Work Produced (Wturbine)
Pump Work (Wpump)
Mass flow rate (Mdot)
In order to optimize the design variable, the fmincon MATLAB function was utilized. The function used the interior point algorithm to deal with constraint violations. Originally, the SQP algorithm was used, but this resulted in poor optimization results.
Optimized Solution
In the most optimized solution (trial 1), we saw a peak temperature of 412, a temperature difference of about 120, and pressure around 357. This resulted in a work output of 1100 MW of power and a 24.46% efficiency. Key aspects to note about the optimized parameters include the temperature drop. The temperature drop was 119.7K, which is almost equal to the upper bound of this parameter at 120K. This makes sense as the ideal efficiency of a Rankine cycle is equal to (1-Tc/Th). This means that we would expect higher efficiencies from greater temperature differences.
With the equation listed above, the ideal efficiency at these operating temperatures is 29%, which is a little bit higher than our final solution, as expected.
Differences across Solutions
The 5 trials run gave very different solutions throughout. The actual efficiency across trials varied treating with a peak efficiency of 24.46% and a bottom efficiency of 13.50%. This indicates that there are likely many local minima throughout the design space.
All the trials except for 1 saw a temperature drop very close to 120K. As discussed above, this is expected to help maximize efficiency. There was one trial that only saw a drop of 66K, but this resulted in the lowest efficiency of any optimized solution.
The pressure seemed to have less of an impact on the final result. There were dramatic swings in pressure from as low as 357 kPa up to 2300 kPa. There was no clear trend between the difference in pressure and the efficiency.
Model Verification and Validation
There were many methods used to verify and validate the model.
First, the efficiency was compared to the ideal efficiency, which showed it was slightly lower than expected.
The net energy of the system was calculated and was deemed to be approximately 0, as it was 16 orders of ten smaller than the cooling rate and electrical work rate. This ensured that the model properly tracked all energy flows in and out of the system.
Finally, the entropy generation rate was calculated, and this value was positive as expected, indicating that the efficiencies were properly integrated into the system.
Compared to a standard Rankine cycle, which has an efficiency of approximately 35%, this cycle is quite inefficient. This could be do to several issues, including the simplicity of the model. Other Rankine cycles have multiple turbines, which allow for a greater enthalpy drop across the turbines. They also have reheating elements to further maximize work output.
function [Work_turb_spec, Work_pump_spec, Qreac_spec, Qcool_spec] = RakineCycle(Th, Tdrop, Pevap)
fluid = 'water';
first_law_pump = .92; %first law efficiency of the pump
first_law_turbine = .92; %first law efficiency of the turbine
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Stage 4 Values
%%%%%%%%%%%%%%
temp_4 = Th;
pressure_4 = Pevap;
enthal_4 = refpropm('H', 'T', temp_4, 'P', pressure_4, fluid);
entropy_4 = refpropm('S', 'T', temp_4, 'P', pressure_4, fluid);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Stage 1 Values
%%%%%%%%%%%%%%%
temp_1 = Th-Tdrop;
enthal_1_prime = refpropm('H', 'T', temp_1, 'S', entropy_4, fluid);
enthal_1 = enthal_4-(enthal_4-enthal_1_prime)*first_law_turbine;
entropy_1 = refpropm('S', 'T', temp_1, 'H', enthal_1, fluid);
pressure_1 = refpropm('P', 'T', temp_1, 'S', entropy_1, fluid);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Stage 2 Values
%%%%%%%%%%%%%%%
temp_2 = temp_1;
quality_2 = 0;
pressure_2 = pressure_1;
entropy_2 = refpropm('S', 'T', temp_2, 'Q', 0, fluid);
enthal_2 = refpropm('H', 'T', temp_2, 'Q', 0, fluid);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Stage 3 Ideal Values
%%%%%%%%%%%%%%%
pressure_3_prime = Pevap;
entropy_3_prime = entropy_2;
temp_3_prime = refpropm('T', 'P', pressure_3_prime, 'S', entropy_3_prime, fluid);
enthal_3_prime = refpropm('H', 'T', temp_3_prime, 'P', pressure_3_prime, fluid);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Stage 3 Values
%%%%%%%%%%%%%%
enthal_3 = enthal_2 + (enthal_3_prime-enthal_2)*first_law_pump;
pressure_3 = pressure_3_prime;
temp_3 = refpropm('T','P', pressure_3, 'H', enthal_3, fluid);
entropy_3 = refpropm('S','P', pressure_3, 'H', enthal_3, fluid);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Key Value Calculations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Work_turb_spec = (enthal_4-enthal_1);
Work_pump_spec = -(enthal_3-enthal_2);
Qreac_spec = (enthal_4-enthal_3);
Qcool_spec = (enthal_2-enthal_1);
%Sgen = (entropy_3 - entropy_3_prime + entropy_1 - entropy_1_prime)*mdot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
function fmin = Performance(x)
Th = x(1);
Tdrop = x(2);
Pevap = x(3);
Qreac = 4500;
[Work_turb_spec, Work_pump_spec, Qreac_spec, Qcool_spec] = RakineCycle(Th, Tdrop, Pevap);
mdot = Qreac/Qreac_spec;
Work_turb = Work_turb_spec*mdot;
Work_pump = Work_pump_spec*mdot;
firstLawEff = (Work_turb+Work_pump)/Qreac;
fmin = -firstLawEff;
end
function [c, ceq] = constraints(x)
Th = x(1);
Tdrop = x(2);
Pevap = x(3);
Qreac = 4500e6;
max_flow = 125000/60; %kg/sec
min_flow = 1000/60; %kg/sec
[Work_turb_spec, Work_pump_spec, Qreac_spec, Qcool_spec] = RakineCycle(Th, Tdrop, Pevap);
c(1) = -(Qreac_spec)+Qreac/max_flow;
c(2) = (Qreac_spec)-Qreac/min_flow;
ceq = [];
end
min_temp_high = 350;
max_temp_high = 500;
min_temp_drop = 20;
max_temp_drop = 120;
min_pevap = 100;
max_pevap = 10000;
lb = [min_temp_high, min_temp_drop, min_pevap]; % Lower bounds
ub = [max_temp_high, max_temp_drop, max_pevap]; % Upper bounds
% Initial guesses (5 different starting points)
x0_set = [
375, 100, 1000;
457, 30, 900;
400, 103, 4000;
365, 35, 9000;
435, 65, 500;
]
results = zeros(5, 4);
% Run optimization
% Optimization options
options = optimoptions('fmincon', 'Algorithm', 'interior-point', 'Display', 'iter');
% Store results
for i = 1:5
% Select initial guess
x0 = x0_set(i, :);
% Run optimization
[x_opt, fval] = fmincon(@Performance, x0, [-1,1,0], [-293], [], [], lb, ub, @constraints, options);
% Store results
results(i, :) = [x_opt, fval];
end