Model an Ice Maker by Using Mode Charts in a Custom Simscape Block

Design the Ice Maker

The ice maker freezes water in 0.9 L batches. Each batch starts at 298.15 K and freezes inside the compartment at 255 K. The ambient air outside the ice maker stays at 298.15 K.

Specify Ice Compartment Geometry and Thermal Characteristics

The ice maker consists of a refrigeration cycle and an ice compartment. The freezing water is in an ice tray within the ice compartment. Specify parameters that define the ice compartment, tray geometry, and their thermal properties:

h_tray = 0.03; % [m] Height of ice tray
w_tray = 0.3; % [m] Width of ice tray
d_tray = 0.1; % [m] Depth of ice tray
t_tray = 0.002; % [m] Thickness of ice tray wall
V_water = h_tray * w_tray * d_tray; % [m^3] Volume of freezing water

k_tray = 0.2; % [W/m*K] Thermal resistivity of ice tray wall
Cp_tray = 20; % [kJ/kg*K] Specific heat capacity of ice tray wall
rho_tray = 1000; % [kg/m^3] Density of ice tray wall

h_compartment = 0.5; % [m] Height of ice compartment
w_compartment = 0.3; % [m] Width of ice compartment
d_compartment = 0.4; % [m] Depth of ice compartment
t_compartment = 0.05; % [m] Thickness of ice compartment wall

k_compartment = 0.02; % [W/m*K] Thermal conductivity of ice compartment wall
Cp_compartment = 1500; % [kJ/kg*K] Specific heat capacity of ice compartment wall
rho_compartment = 30; % [kg/m^3] Density of ice compartment wall

h_stillAir = 5; % [W/m^2*K] Heat transfer coefficient of still air
h_blownAir = 30; % [W/m^2*K] Heat transfer coefficient of blown air

Calculate Refrigeration Requirements

Next, calculate the heat that the refrigeration cycle must remove to freeze the water. This value depends on the thermodynamic properties and mass of each batch of water. First, define the relevant thermophysical properties of the water:

rho_liquid = 997; % [kg/m^3} Density of water
Cp_liquid = 4.184; % [kJ/kg*K] Specific heat capacity of liquid water
dH_fusion = 334; % [kJ/kg] Latent heat of fusion of water
T_initial = 298.15; % [K] Initial temperature of water
T_freezing = 273.15; % [K] Freezing temperature of water at 1 atm
T_environment = 298.15; % [K] Temperature of ambient environment outside ice maker
T_set = 255; % [K] Ice compartment temperature set point

Calculate the mass of water in each batch and the temperature difference between the initial and freezing temperatures of the water.

m_water = V_water * rho_liquid; % [kg] Mass of water
dT_cooling = T_freezing - T_initial; % [K] Temperature change to freezing

Next, compute the total heat the ice maker must remove to freeze each batch of water. This value includes the sensible heat and the latent heat of fusion,

$$Q_{freeze}=\Delta H_{freeze}=m_{water}\left(C_{p,l}\Delta T_{cooling}-\Delta H_f\right)$$,

where:

dH_freeze = m_water * (Cp_liquid * dT_cooling - dH_fusion); % [kJ] Enthalpy change to freeze one batch
Q_freeze = dH_freeze; % [kJ] Heat of freezing

Calculate the cooling rate,

$${\underset{}{\overset{\dot{}}{Q}}}_{freeze}=\frac{Q_{freeze}}{\Delta t}$$,

where a dot represents a rate per time.

dt_freeze = 7200; % [s] Estimated time to freeze one batch of ice
Qdot_freeze = Q_freeze / dt_freeze % [kW] Cooling rate for freezing

Qdot_freeze = 
-0.0547

Design the refrigeration cycle

Next, use the cooling requirement to design the refrigeration cycle. Start by defining the subcooling and superheat margins.

dT_subcooling = 10; % [K] Subcooling
dT_superheating = 10; % [K] Superheating

To understand how the refrigerant behaves, use a pressure-enthalpy (P-H) diagram for R-134a. To view a P-H diagram, you can open an empty model, add a Two-Phase Fluid Predefined Properties (2P) block, and double-click on the block to open the block dialog. Set Two-phase fluid to R-134a. In the Plots section, click the Plot button for the Two-phase fluid properties (contours). In the plot window, in the top-left corner, click the Enthalpy-x-axis button.

Using the design specifications, select four key points on the diagram:

These points define the thermodynamic cycle. For more information about how you select these four points, see Model a Refrigeration Cycle. Use the enthalpy values at these points to calculate the coefficient of performance (COP),

$$COP=\frac{h_1-h_4}{h_2-h_1}$$,

where $$h_1$$, $$h_2$$, $$h_3$$, and $$h_4$$ are the enthalpies of the refrigerant at points 1 through 4.

h1 = 395.7; % [kJ/kg] Enthalpy at point 1
h2 = 443.7; % [kJ/kg] Enthalpy at point 2
h3 = 227.5; % [kJ/kg] Enthalpy at point 3
h4 = 227.5; % [kJ/kg] Enthalpy at point 4
COP = (h1 - h4) / (h2 - h1) % Coefficient of Performance

COP = 
3.5042

The evaporator removes heat at a rate equal to the cooling rate of freezing water when the ice maker is perfectly insulated:

$${\underset{}{\overset{\dot{}}{Q}}}_{evap}={\underset{}{\overset{\dot{}}{Q}}}_{freeze}$$,

$$\underset{}{\overset{\dot{}}{W}}=\frac{{\underset{}{\overset{\dot{}}{Q}}}_{evap}}{COP}$$,

$${\underset{}{\overset{\dot{}}{Q}}}_{cond}={\underset{}{\overset{\dot{}}{Q}}}_{evap}+\underset{}{\overset{\dot{}}{W}}$$,

where:

Qdot_evap = Qdot_freeze; % [kJ/hr] Heat rate of evaporator
Wdot = Qdot_evap / COP % [kJ/hr] External work requirement to refrigeration cycle

Wdot = 
-0.0156

Qdot_cond = Qdot_evap + Wdot % [kJ/hr] Heat rate of condenser

Qdot_cond = 
-0.0703

You can use these heat rates to calculate the mass flow rates of air through the condenser and evaporator,

$$\underset{}{\overset{\dot{}}{m}}=\frac{\underset{}{\overset{\dot{}}{Q}}}{C_{p,air}\Delta T}$$,

where:

Cp_air = 1.00; % [kJ/kg*K] Specific heat capacity of air
dT_evap = 10; % [K] Temperature drop of air across evaporator
dT_cond = 10; % [K] Temperature drop of air across compressor
mdot_evap = Qdot_evap / (Cp_air * dT_evap) % [kg/hr] Mass flow rate of air through evaporator

mdot_evap = 
-0.0055

mdot_cond = Qdot_cond / (Cp_air * dT_cond) % [kg/hr] Mass flow rate of air through compressor

mdot_cond = 
-0.0070

Build the Model

Next, build a model to represent the ice maker. The model includes three domains:

It also includes a custom block that models freezing water and the effect of the water on the rest of the ice maker.

Model the Refrigeration Cycle

To ensure the condenser expels heat into the ambient environment, verify that the value of the Condenser External Fluid Nominal inlet temperature is less than the condenser refrigerant outlet temperature. The condenser refrigerant outlet temperature is equal to the value of the Nominal condensing (saturation) temperature minus the Nominal condenser subcooling. To ensure the evaporator can absorb heat, verify that the value of the Evaporator External Fluid Nominal inlet temperature is greater than the evaporator refrigerant outlet temperature. The evaporator refrigerant outlet temperature is equal to the value of the Nominal evaporating (saturation) temperature plus the Nominal (static + opening) evaporator superheat.

Model the Ambient Environment

Build the Ice Compartment

Use the moist air and thermal domains to model the ice compartment. Create a custom Simscape block to model the freezing water and its effects on the ice maker.

Next, build the ice maker thermal network and connect it to port H of the Constant Volume Chamber (MA) block. The thermal network models heat transfer from the ambient environment to the ice compartment. Include an unconnected physical signal in the model to connect to the custom Simscape block. In this system, the:

Create a Custom Simscape Block

You can use the Simscape language to create a custom block that models the discontinuous thermodynamic properties of freezing water. For more information about creating custom Simscape blocks, see Creating Custom Components. For this example, open the freezingWater.ssc file:

open('freezingWater.ssc')

This file uses the modecharts construct to represent freezing and melting water:

Complete the thermal network by adding the custom Simscape block using the Simscape Component block, Temperature Sensor blocks to measure the temperature of the ice compartment and freezing water, and Outport and PS-Simulink Converter blocks to pass Simulink signals out of the subsystem. Name the subsystem Thermal Network.

Add Control Feedback and Scopes

Open the IceMaker model to see the complete model.

Observe the Results

Simulate the IceMaker model and plot the temperature of the ice compartment and freezing water.

mdl = 'IceMaker';
open_system(mdl);
sim(mdl);
t_model = tout/3600; % [hr] Model time
T_ice_compartment = simlog_IceMaker.Ice_Compartment.T_I.series.values; % [K] Ice Compartment Temperature
T_freezing_water = simlog_IceMaker.Thermal_Network.Freezing_Water.T.series.values; % [K] Freezing water temperature
plot(t_model, T_freezing_water);
hold on;
plot(t_model, T_ice_compartment);
legend('T_{water}', 'T_{compartment}');
xlabel('Time (hr)');
ylabel('Temperature (K)');
hold off;

Both the ice compartment and water get colder for the first 40 minutes. The water reaches its freezing point around 40 minutes and both the ice compartment and water temperatures stay constant as the water freezes. After the water is completely freezes, the temperatures start to drop until the system turns off the compressor and fans. The refrigeration cycle turns on and off to keep the ice compartment at the temperature setpoint for the rest of the simulation.

Next, plot the phases of the water.

phase = simlog_IceMaker.Thermal_Network.Freezing_Water.internal_mode_var_mc1__.series.values;
plot(t_model, phase, 'LineWidth', 2);
ylim([0, 4]);
yticks([1, 2, 3]);
yticklabels({'Solid', 'Liquid', 'Solid <-> Liquid'});
xlabel('Time (hr)');
ylabel('Phase');

The freezing water block is in the liquid phase for approximately 40 minutes. The water then transitions between liquid and solid states before becoming solid. These phases and their timing are consistent with the discontinuous temperature changes in the freezing water block.