Supermarket refrigeration system (Muscod)

The differential equations for the Supermarket refrigeration system in Muscod code read as follows

// number of compressors
#define NVALVES 2
 
// constants
#define M_GOODS 200.0
#define C_P_GOODS 1000.0
#define UA_GOODS_AIR 300.0
#define M_WALL 260.0
#define C_P_WALL 385.0
#define UA_AIR_WALL 500.0
#define M_AIR 50.0
#define C_P_AIR 1000.0
#define UA_WALL_REF_MAX 4000.0
#define M_REF_MAX 1.0
#define TAU_FILL 40.0
#define T_SH 10.0
#define V_SUC 5.0
#define V_SL 0.08 // 2 display cases - 2 compressors
// #define V_SL 0.095 // 3 display cases - 3 compressors
#define ETA_VOL 0.81
 
// disturbances - day scenario
#define Q_AIRLOAD 3000.0
#define M_REF_CONST 0.2
 
// disturbances - night scenario
// #define Q_AIRLOAD 1800.0
// #define M_REF_CONST 0.0
 
double delta_h = (0.0217*xd[0]*xd[0] - 0.1704*xd[0] + 2.2988)*1e5;
double T_e = -4.3544*xd[0]*xd[0] + 29.224*xd[0] - 51.2005;
double rho_suc = 4.6073*xd[0] + 0.3798;
double rho_suc__P_suc = -0.0329*xd[0]*xd[0]*xd[0] + 0.2161*xd[0]*xd[0] - 0.4742*xd[0] + 5.4817;
double f = (0.0265*xd[0]*xd[0]*xd[0] -0.4346*xd[0]*xd[0] + 2.4923*xd[0] + 1.2189)*1e5;
 
double Q_goods_air[NVALVES];
double Q_air_wall[NVALVES];
double UA_wall_ref[NVALVES];
double Q_e[NVALVES];
double m[NVALVES];
 
double m_in_suc = 0.0;
 
int i;
 
for (i=0; i<NVALVES; i++){
    Q_goods_air[i] = UA_GOODS_AIR*(xd[1 + i*4] - xd[3 + i*4]);
    Q_air_wall[i] = UA_AIR_WALL*(xd[3 + i*4] - xd[2 + i*4]);
    UA_wall_ref[i] = UA_WALL_REF_MAX * xd[4 + 4*i]/M_REF_MAX;
    Q_e[i] = UA_wall_ref[i]*(xd[2 + 4*i] - T_e);
 
    m[i] = Q_e[i]/delta_h;
    m_in_suc += m[i];
}
 
double V_comp = 0.0;
double comp_scale = (double) 1.0/NCOMPS;
V_comp = comp_scale*u[NVALVES]*ETA_VOL*V_SL;
 
// suction pressure
rhs[0] = (m_in_suc + M_REF_CONST - V_comp*rho_suc) / (V_SUC*rho_suc__P_suc);
 
// for each display/valve
for (i=0; i<NVALVES; i++){
 
    // temperatures:
 
    // goods
    rhs[1 + i*4] = - Q_goods_air[i]/(M_GOODS*C_P_GOODS);
 
    // wall
    rhs[2 + i*4] = (Q_air_wall[i] - Q_e[i])/(M_WALL*C_P_WALL);
 
    // air
    rhs[3 + i*4] = ((Q_goods_air[i] + Q_AIRLOAD - Q_air_wall[i]) /(M_AIR*C_P_AIR));
 
    // mass of liquefied refrigerant:
    rhs[4 + i*4] = ((M_REF_MAX - xd[4 + 4*i])/TAU_FILL * u[i] - m[i] * (1 - u[i]));
}