Transient Simulation of an Air-source Heat Pump under Cycling of Frosting and Reverse-cycle Defrosting

Frost accumulation is a common but undesired phenomenon for air-source heat pump (ASHP) systems in winter operations. Reverse-cycle defrosting that applies heat to the outdoor coil by reversing the thermodynamic cycle, is one of the predominant means for periodic removal of the accumulated frost. This paper presents a dynamic modeling framework for ASHPs under cycling of frosting and reverse-cycle defrosting operations. A uniform model structure was applied to frost formation and melting models, which were incorporated into the evapo-rator model, without a need for reinitializing the system when the operating mode switches between heating and defrosting. The developed model was simulated to predict transients of a residential ASHP under multiple cycles of frosting and defrosting operations, and results yielded good agreement with measurements.


Introduction
Frost accumulation on evaporator coil surfaces can significantly degrade the performance of air-source heat pump (ASHP) systems in winter operations.The continued buildup of frost eventually necessitates a defrosting mode to remove the accumulated frost and return the system to its normal operating characteristics.Among various approaches currently employed on ASHP systems, reversecycle defrosting (RCD) that applies heat to the outdoor coil by reversing the thermodynamic cycle, is one of the predominant means because of its easy implementation without additional heat sources.The nature of frosting-defrosting cycling operations imposes significant challenges for ASHP control designs, optimization, fault detection and diagnostics (FDD).Therefore, a simulation tool capable of capturing the system dynamics under mode-switching operations is extremely useful in development and evaluation of improved control and FDD algorithms.
Due to the complicated underlying physical behaviors of frosting and defrosting of ASHPs and significant computational complexities, a very limited number of systemlevel modeling efforts can be found in the open literature.Many studies solely focus on performance of heat ex-changers either with frost formation or during RCD, e.g., (Dopazo et al. 2010;Breque and Nemer 2017;Padhmanabhan et al. 2011).Within the context of complete cycle modeling considering the dynamics on both the refrigerant side and air side, Qiao, Aute, and Radermacher (2017) integrated a one-dimensional frost growth model into an evaporator model coupled to other components of a twostage flash tank vapor injection heat pump system.Transients of the system going through a start-up period, a stable frosting stage, followed by an unstable hunting stage due to the performance degradation by frost accumulation were simulated.Comparisons of the simulation results against experimental data indicate that the developed model can reasonably predict the system responses under frosting conditions.Steiner and Rieberer (2013) simulated the defrosting process of a CO2 heat pump system used for an electric vehicle by integrating a lumped frost melting model into the heat exchanger model.Predictions of the refrigerant dynamics and the compressor power showed reasonable agreements with the measurements.Then the model was used to conduct a parametric analysis on different expansion valve openings during defrosting.It was found that an optimal valve opening existed regarding defrost time and efficiency.However, only the system operation after the refrigerant pressure levels were equalized and the compressor was turned on again in defrost cycle, was simulated omitting the reverse flow transients.Han et al. (2022) carried out a simulation study of a residential heat pump system during defrosting cycles.A multistage frost melting model was coupled to a finite volume heat exchanger model.Results were compared with experimental data and yielded good agreement.Similar to Steiner and Rieberer (2013), both studies concentrated on system dynamics during a short time period of the defrosting cycle, and excluded the refrigerant pressure equalization stage triggering flow reversal.In a subsequent simulation study, Qiao, Aute, and Radermacher (2018) incorporated a multi-stage frost melting model in the evaporator model.With a reversing valve model, the system model was able to capture transients when the refrigerant flow was reversed during initiation and termination of RCD.No experimental validation was provided.
To address the need for a general simulation tool that is directly applicable to control and FDD, this paper presents a dynamic modeling framework of ASHP incorporated with frost formation and melting.The developed cycle model was validated using experimental data collected from a residential ASHP unit operating under multiple cycles of frosting and RCD.Section 2 presents major component models.Implementation of the developed models for simulating the system dynamics under cycling of frosting and defrosting is then described in Section 3. Section 4 reports simulation results and comparisons against the measurements, followed by conclusions of the present work summarized in Section 5.

Heat exchanger model
Much of the modeling efforts have focused on heat exchanger models since the dominant dynamics of a general vapor compression system reside in two-phase heat exchangers.In the present work, a heat exchanger model is integrated with detailed frost formation and melting models to investigate the overall impact of frost on the system performance.To describe dynamics of the refrigerant in heat exchanges, a number of assumptions are required to simplify the model construction of two-phase flow: • The refrigerant flow is one-dimensional with uniform fluid properties at cross sections.
• Changes in kinetic energy and potential energy are negligible; viscous dissipation is negligible; axial heat conduction along the refrigerant flow direction is negligible.
• The liquid and vapor of the two-phase region are in thermodynamic equilibrium.
• Body forces are neglected in the momentum balance.
The finite volume approach, that segments a heat exchanger into an arbitrary number of equally sized control volumes, is utilized for discretization.The refrigerant pressure and enthalpy as a pair of independent thermodynamic properties are chosen as state variables to express the discretized governing equations.A staggered grid scheme is adopted to solve the balance equations in those control volumes (CVs) (Laughman et al. 2015).As shown in Figure 1, the mass and energy balances are solved in the upper grid, referred to as volume cells, where thermodynamic properties are determined, while the momentum balances are solved in the lower grid, referred to as flow cells, where dynamics of the mass flow rate are evaluated for neighboring volume cells.
With a discretization of n control volumes for a heat exchanger, n − 1 momentum balances are formed, where each of the two boundary flow cells has an extended length of half the cell.In this way, the momentum balances solves interface mass flow rates on the n − 1 inner edges of the volume cells, leaving mass flow rates on the outer edges ( ṁ1 , ṁn+1 ) as boundary conditions to the grid.Meanwhile, the first and last volume cells expose the thermodynamic states.As a result, connection with other flow devices avoids the solution of large nonlinear systems for algebraic pressures (Franke, Casella, Sielemann, et al. 2009).The discretized governing equations of refrigerant mass, momentum and energy balances are given respectively in (1)-( 3), where V i denotes the volume of each volume cell, L i denotes the length of a flow cell, vi denotes the average velocity of a volume cell, Qi is the convective heat transfer rate with the metal wall, H i denotes the enthalpy flow rate computed at the edge of a volume cell.The upwind difference scheme is used to approximate thermodynamic quantities at the volume cell interfaces.Therefore, the enthalpy flow rate under the nominal flow direction can be calculated as In determining quantities transported via convection such as the specific enthalpy, it is important to note that upstream and downstream are concerned with a certain flow direction.Since this work deals with reverse-flow modeling, flow dependent values need to be evaluated as the mass flow rate crossing zero.In Modelica-based modeling, this is handled by using stream variables that can describe the bi-directional flow in a numerically reliable way when singularity arises at zero mass flow (Franke, Casella, Otter, et al. 2009).Friction force across a flow cell is evaluated by the equivalent frictional pressure drop ∆p f ,i : where A s,i is the surface area of a control volume.Note that a sign function of the mass flow rate is applied since the friction force is always in the direction opposite to the flow.
Assuming a uniform temperature of the tube wall and associated fins, conservation of energy for the coil metal structure can be derived as where Q f ,i denotes heat conduction with the frost layer.
The air temperature and humidity profile can be derived from one-dimensional, quasi-steady-state mass and energy balances with a uniform surface temperature and inlet air conditions, T ai,out = T as,i + (T a,in − T as,i )e −Ntu (7) where Ntu is the number of transfer units for sensible heat transfer, ω as,i is the saturated humidity ratio evaluated at the surface temperature T as,i .The heat and mass transfer analogy through the Lewis number Le 2/3 = 0.9 is adopted to correlate convection coefficients of heat and mass transfer, which is valid in both cases of condensation and sublimation (Bergman et al. 2011;Leoni et al. 2017).The total heat transfer consisting of sensible and latent parts can be obtained by Qi = ṁa,i c p,a (T ai,out − T a,in ) + ṁa,i (ω ai,out − ω a,in )∆h lat (10) where ∆h lat represents the latent heat of condensation or sublimation according to the surface temperature.When the coil fan is off, heat transfer between the ambient air and coil through natural convection becomes dominant.The heat transfer rate is calculated assuming a uniform air temperature around the coil and neglecting latent heat transfer, The value of the Grashof number over the square of the Reynolds number Gr/Re 2 measures effects of buoyancy forces relative to inertial forces, which can be applied to indicate the dominant effect when the coil fan is energized and shut off.In this work, ( 7) -( 10) are used when Gr/Re 2 ≤ 0.1, otherwise the air-side heat transfer rate is calculated from (11).
The prediction of frost formation on the outdoor coil is essential to characterize the performance of an ASHP system under frosting conditions due to blockage of the air flow passage and additional thermal resistance.To make a numerical model tractable, it is often assumed that frost growth normal to the cold surface is of primary interest.The frost layer growth and densification rates can be determined by solving the heat and mass diffusion equations.As shown in Figure 2, at the frost-air interface, the total heat flux from the air stream q ′′ tot is composed of sensible and latent heat due to the temperature and humidity differences, where α h denotes the heat transfer coefficient, α m denotes the mass transfer coefficient, T f s denotes the frost surface temperature, ω f s denotes the saturated humidity ratio at the frost surface, ∆h sg denotes the water sublimation heat.
Note that a heat transfer rate obtained from (12) represents the same heat transfer rate as (10) under frosting conditions.Indeed, the heat transfer rate will be implemented as a connection variable of the air flow model and the frost formation model, such that a equation system can be formed to solve for intermediate variables (e.g., frost surface temperature).It also gives the amount of heat conducted to the metal structure ( Q f ,i ) in ( 6), assuming that the frost growth process is quasi-steady-state.When the frost layer is of thickness δ f ,i and lumped density ρ f ,i , mass conservation of the frost layer can be formed as which can rewritten as (14) indicates that the total water mass flux ṁ′′ a,i is divided into two components: ṁ′′ δ contributes to an increment of the frost thickness and ṁ′′ ρ , which is diffused into the frost layer, contributes to its densification.Consider a differential control volume of thickness dx within the frost layer, the heat and mass diffusion can be formulated and then solved associated with boundary conditions at the wall surface and the frost-air interface.For the sake of brevity, the solution procedure is not included in the present paper.Refer to (Qiao, Aute, and Radermacher 2017) for derivations to obtain the mass fluxes in ( 14).
Since the model is derived following a quasi-steadystate assumption, the frost thickness and density are updated at a fixed time step ∆t which enables the overall frost growth process to remain transient over time.
The rather stochastic frost melting process is idealized such that the defrost process progresses through several predictable stages.The present work applies a five-stage defrost model subdividing the overall process into preheating, melting start, melting, vaporizing, and dry heating (Qiao, Aute, and Radermacher 2018).The governing equations of each stage are presented in Appendix.The model is basically established by a lumped-capacitance analysis along with energy and mass conservation at interfaces, and switches between different stages according to the current temperature profile and phase presence.The overall melting process is illustrated in Figure 3.
A switching algorithm based on the Fuzzy logic is implemented (Kim et al. 2021).Since transitions between different stages are triggered by values of the wall temperature T w,i , frost and a water film thickness δ f ,i and δ water,i , three linguistic variables T w,i − T melt , δ water,i /δ water,max and δ f ,i are defined.T melt is the melting temperature of water (e.g., 273.15 K). δ water,max is a constant maximum thickness of melted water retained on the coil due to surface tension.Fuzzy numbers and their associated membership functions are constructed for each linguistic variable.For example, Figure 4 shows Fuzzy numbers N, P and membership functions for T w,i − T melt .Denote the governing equations of each stage as f j (•) ( j = {1, 2, 3, 4, 5}).The system dynamics can be determined based on linguistic descriptions of IF-THEN rules (Rule i, i={1,2,3,4,5}), which are subsequently converted into numerical outputs (known as the defuzzification process).Take the rule for the first stage (preheating) as an example, where x is a state vector consisting of temperatures and thicknesses of different phase presences.Then a numerical output associated with the rule can be obtained which is used to compute a weighted combination of dy-namics from all the stages as described in (18).
In this way, the IF-THEN conditional rules can be eliminated, and the overall dynamics are evaluated by weighting the dynamics of all stages.The proposed switching algorithm systemically covers all possible switches and avoids discontinuities in switches.Consequently, the frost melting model can be used to evaluate dynamics of the frost density and thickness during RCD, where f ρ,melt (•) and f δ ,melt (•) are the corresponding governing equations.

Component models
A quasi-static model is developed to describe performance of a variable-speed compressor using efficiency maps.The refrigerant mass flow rate is determined by where V s is a fixed displacement, N is the number of rotations per minute.The volumetric efficiency η v is regressed by a polynomial of pressure ratio and compressor speed, where a 0 − a 3 are constant coefficients.The power consumption is estimated using an isentropic efficiency which is fitted by The refrigerant discharge state can be determined by forming an energy balance for the compressor, where the heat loss ratio has a strong dependence on the operating conditions and the ambient temperature T amb ,  The expansion process is assumed to be isenthalpic.The mass flow rate is determined by where A v is the varying valve opening area, and C d is the discharge coefficient that accounts for corrections to the mass flow rate at different operating conditions.Note that the valve is modeled as a bi-directional device that the inlet and outlet can switch when flow is reversed.The opening area is adjusted based on superheat control.An empirical correlation based on a power law is utilized to estimate the discharge coefficient (Liu, Cai, and Kim 2022): where T sc is the subcooling at the valve inlet, T c is the critical temperature of the refrigerant, and φ is the normalized valve opening.An accumulator is modeled assuming that vapor and liquid inside are saturated and in thermal equilibrium, which lead to the mass and energy conservation equations as shown in ( 29)-( 30), Note that these saturated properties are solely dependent on the refrigerant pressure, thus the lumped pressure can be selected as a state variable.Given the constituent relation that the total volume of the tank is occupied by the saturated liquid and vapor volumes V f +V g = V acc , either one of the volumes can be selected as another state variable.The exit flow is assumed to be saturated vapor to account for the fact that the superheated vapor from the evaporator mixes with the liquid refrigerant stored inside the accumulator.
The coil fan is described via a polynomial characteristic curve: where V is the total volume flow rate.Since the coil pressure drop is assumed to be solely a result of friction, a hydraulic equilibrium is established between the fan and the coil, which leads to the system of equations below: where ∆p a,i is the air side pressure drop of the ith control volume.Along with the mass balance where ṁa,i represents the air flow distribution of each control volume, a non-linear algebraic equation system is formed to solve the air flow maldistribution under nonuniform frost formation.
A reversing valve model is essential to switch the heat pump system model between heating and defrosting modes.It is challenging to model the exact physical process of a reversing valve, since it mainly involves mechanical movements.In the current work, two pairs of check valves are used to represent the reversing valve model (Qiao, Aute, and Radermacher 2015).As shown in Figure 5, each port of the reversing valve has a pair of check valves to control the flow direction.One is kept fully opened, while the other is closed.When the on-off statuses of them are switched, the refrigerant flow direction is reversed.A check valve is modeled as where A v is the fully opened area, φ is valve opening degree normalized between 0 and 1.When reverse flow needs to be triggered, opening and closing a check valve is achieved by changing the valve opening.However, this process is not instantaneous during mode switching.For numerical stability concerns, a first-order filter is applied to smooth the transition, where φ SP is the opening setpoint, which is either 1 or 0 depending on the on/off status of each check valve.The heat loss from the high pressure side is primarily driven by the temperature difference between the discharge gas and the suction gas, and can be characterized by heat transfer loss/gain coefficients (Damasceno, Rooke, and Goldschmidt 1991).

Model implementation
To simulate the cycling behavior of non-uniform frost growth and melting on the outdoor coil, the frost formation and melting models are integrated into each discretized control volume of the heat exchanger model.four elements can be identified across the predominant heat flow direction, which is perpendicular to the refrigerant flow direction: the refrigerant flow, the coil metal structure, the frost layer, and the air flow.Note that the frost layer in Figure 6 represents simultaneously running the frost growth and melting models.These frost models are implemented as generic control volumes to maintain a consistent model structure, thus the frost layer dynamics coupled to the heat pump system dynamics can switch between frost growth and melting, according to the system operating mode.The standard HeatPort extended from Modelica.Thermal.HeatTransfer is employed at the interface between the frost and air flow.However, at the interface between the frost and metal wall, a heat source is utilized for the metal wall model whose input switches between computed heat flow rates from the frost growth and melting models.Furthermore, the lumped metal wall temperature is considered as an input to the frost models.
Since the frost density and thickness characterize its dynamics, they are selected as internal states.That means, the current states of frost density and thickness at each time step are treated as known and inputs to the frost models, while their dynamics are obtained as outputs for predictions over the next time interval.
After that, the overall frost dynamics are evaluated as a weighted combination of dynamics obtained from these two models.Nonetheless, due to different assumptions applied for derivations, it is not straightforward to implement the scheme before having a uniform structure among them.Since the frost formation model is derived in a quasi-steady-state fashion, and updates the frost density and thickness in a fixed time step as shown in (15, 16), the progress made over time can be considered as reference values.Then a first-order filter is applied to track the reference values where ρ ref,i and δ ref,i are quantities evaluated in (15, 16), ρ f ,i and δ f ,i are states representing the frost behavior, τ is a time constant, which should be selected such that the signal tracking is faster than update of the reference values.In this way, dynamics of the frost formation can be described in a continuous-time domain.Note that these reference values evaluated in the frost formation model are non-decreasing.This is because throughout the derivation the mass flux is assumed to be non-negative, meaning that no melting behavior is considered in the frost formation model.Therefore, in case of a mode switch from defrosting to frosting (heating), the reference values should be reset to the states at the termination of frost melting as initial conditions for the next frosting period.The overall frost dynamics are evaluated by where f ρ,melt (ρ f ,i , δ f ,i ) and f δ ,melt (ρ f ,i , δ f ,i ) are the dynamics obtained from the frost melting model, φ is a variable that switches the frost dynamics between frost formation and melting.The reversing valve opening in ( 35) can be used for this purpose.The scheme allows the frost formation and melting models to run simultaneously during online simulations, which enables modeling the cycling of frosting and defrosting without interruptions.Consequently, component models are interconnected through the standard fluid connectors to form a cycle model.Refrigerant properties are evaluated using a regression approach based on Artificial Neural Networks (Ma, Kim, and Braun 2018).

Simulation results and validation
This section presents simulation results and experimental of the developed cycle model for capturing transient characteristics of a residential heat pump system under multiple cycles of frosting and reverse-cycle defrosting (RCD).The heat pump unit is a commercial variable-speed system having a nominal capacity of 2 tons that utilizes R410A as the working fluid.The split system was installed in a pair of independently controlled psychrometric chambers, where the indoor unit was ducted to a nozzle box for air flow rate measurements.The experimental setup was fully instrumented with T-type thermocouples, pressure transducers, and a coriolis flow meter on the refrigerant loop.T-type thermocouple grids and chilled mirror dew point meters were installed for air-side measurements.Power meters were installed to measure power of the indoor and outdoor units as well as the outdoor fan.The temperature set points of the indoor and outdoor conditions were 291.5 K and 271 K, respectively.The relative humidity was set to 40% for the indoor and 85% for the outdoor.The unit was turned on in heating mode after these conditions were met.At low ambient temperature, a system built-in controller determined operating modes between heating and defrosting based on a fixed-interval defrost logic.During defrost cycles, the compressor is operated at a lower speed than the heating mode but never shut off.The EXV opening is regulated to control the evaporator outlet superheat in the heating mode but is fully opened in the defrosting as well as cooling modes.The indoor and outdoor coil fan speeds are controlled according to the compressor speed, however, the outdoor fan is off in defrost cycles.The defrost interval, representing the compressor run time between two defrost cycles, was set to 120 minutes as the longest.However, upon startup of the unit, the first defrost interval was defaulted to 30 minutes.Regarding this defrost control pattern, a test that lasted three hours ran through a start-up with minor frost accumulation for 30 minutes, followed by the first defrost cycle, then normal heating operation for 120 minutes with considerable frost buildup, followed by a second defrost cycle, and then another heating operation period of about 10 minutes.Even though the inlet air conditions were kept constant in both chambers, various durations of each frosting cycle and different defrost cycle durations yielded different cycling transients of the system, since frost buildup and removal are inherently transient processes.
Each heat exchanger was discretized into 30 control volumes based on a preliminary sensitivity analysis that the number of control volumes is determined until the solution is invariant to further increases of control volumes.The heat pump cycle model was initialized using the pressure and temperature measurements across various components.The initial refrigerant pressure of each component was almost at the same level, though the temperature differed depending on the room temperature.As a result, initial conditions of most components reflect gas phase refrigerant presence, while the refrigerant residing in the outdoor coil was initialized as two-phase.Without access to the liquid level of the suction accumulator, steady-state initialization was performed using the cycle model to determine the vapor and liquid volumes.Simulations were carried out to predict the complete cycling operations of frosting and reverse-cycle defrosting lasting three hours.and suction pressures.As the compressor is powered on, its speed increases to the desired setting regulated by the control board in 3 minutes as shown in Figure 8, which results in the rapid rise of the discharge pressure.The EXV is fixed to a preset position during this period before the superheat control is enabled.After that, rise of the discharge pressure slows down since the compressor is running at a pre-defined speed.As the EXV is closing to achieve the superheat set point temperature, the resulting mass flow imbalance between the compressor and EXV leads to the continuing pressure increase.As mentioned before, the first defrost cycle is defaulted to 30 minutes after the initial power-up of the compressor.When the system initiates defrost operation, the EXV is fully opened, which leads to a dramatic pressure change due to the pressure equalization across the EXV.Shortly the reversing valve is energized, which further brings down the discharge pressure, and slightly increases the suction pressure due to the connection to the high pressure side when the valve is switching port positions.After the flow is reversed, the discharge gas is pumped into the outdoor coil, which operates as a condenser, and the liquid refrigerant residing in the indoor coil, which operates as an evaporator, vaporizes due to the abrupt pressure drop and flows towards the vapor line.As a consequence, the discharge and suction pressures increase until termination of the defrost mode.The first defrost cycle lasts about 5 minutes.
Then the reversing valve is energized again to switch the flow directions, that reduces the pressure difference between the high pressure side and the low pressure side.Meanwhile, the EXV opening is reset to an initial position, and then adjusted following superheat control.A similar trend of the discharge pressure can be observed as during the start-up, a rapid rise due to the compressor speed regulation, followed by a gradual increase from 40 min to around 100 min as the refrigerant mass flow is being balanced between the vapor-line and liquid-line.Though frost starts to form shortly after the system switches back to heating mode, the refrigerant dynamics during this period are dominated by the mass flow imbalance without noticeable impact from the frost accumulation.However, performance degradation due to frost formation becomes more evident since 100 min until initiation of the next defrost cycle (160 min).The reduced air flow due to frost blockage results in a slight decline of the evaporating pressure (see zoom-in plot in Figure 7), yet a significant decline of the discharge pressure which is more sensitive to changes on the suction side.The second defrost cycle is initiated as the accumulated compressor run time reaches the preset defrost interval.The refrigerant pressures exhibit a similar trend as the previous defrost cycle, triggered by the same consecutive actuation of the EXV and reversing valve, although the pressure rises are relatively smooth because a large amount of heat is taken to melt frost that slows down the refrigerant energy storage.This defrost cycle lasts 10 minutes, after which the system returns to heating mode.The simulation results agree sufficiently well with the measurements, which demonstrate that the developed cycle model is able to capture the complicated system dynamics under the frosting-defrosting cycling operations.The discrepancy between the predictions and measurements during defrost cycles is due to model simplification and unknown response times of actuators including EXV and reversing valve that trigger the mode switches.For instance, rises in the predicted suction pressure around 35 min and 170 min during mode switches are not observed in the measurements.This is due to implementation of the reversing valve model using a system of check valves.When check valves are closed to switch flow direction, the rapid mass flow decrease almost results in a pressure equalization of the discharge and suction pressures.
Figure 9 shows the refrigerant mass flow rate prediction in the indoor liquid-line.It is important to note that the refrigerant flow can be bi-directional in most of the components depending on the operating mode.Following the sign convention that inflow quantities are positive and outflow quantities are negative, a positive mass flow rate at the inlet of a component indicates the nominal flow direction (e.g., heating mode), while a negative flow rate indicates reverse flow (e.g., defrosting mode).The sign convention applies to both mass flow and heat flow in the following context.In heating mode, the refrigerant mass flow rate undergoes abrupt changes when the compressor speed is regulated, otherwise fluctuates according to the discharge and suction properties.From 87 min until initiation of the second defrost cycle at 160 min, the mass flow slightly declines during the stable heating operation due to an insignificant decline of the pressure ratio.When the system switches to defrost mode at 30 min and at 160 min, a spike of the liquid-line flow rate occurs as the EXV is fully opened.The reversing valve is then engaged to trigger the flow reversal.The mass flow rate steadily increases during the defrost cycle since expansion valves are kept fully open as the reversed refrigerant flow develops, until the flow is reversed at the end of each cycle.Note that the output signal of the mass flow meter saturates to zero when the flow direction is opposite to the configured direction in the test.That means, the flow rate measurement during defrost cycle is unavailable.The comparison between the simulated and measured mass flow shows that the model can capture major flow transients in heating operations and the flow reversal between mode switches.The discrepancy during each start-up in heating mode can be attributed to prediction errors of quasi-static models of the compressor and EXV that were correlated using steady-state data.As can be seen from 80 min to 160 min, the nearly steady-state mass flow after start-up transients is well captured.Simulation results of the indoor unit air-side capacity and the compressor power are revealed in Figure 10.The indoor unit fan is continuously running throughout the test.The air flow rate climbs after the fan motor is powered, and reaches a constant value based on the measurements.The air-side capacity evolves corresponding to the refrigerant dynamics.Since the compressor operates at a nearly fixed speed, a decline of the heating capacity can be observed during 100 min to 160 min, when the refrigerant condensing pressure drops due to the frosting operation.In defrost mode, the indoor coil operates as an evaporator and the heat flow direction is reversed compared to heating mode.During this period, heat is removed from the indoor space, which can potentially result in thermal discomfort.The model captures dynamics of the indoor capacities under both condenser and evaporator modes well.
The change of compressor power consumption follows regulation of the compressor speed, while also varying due to the refrigerant dynamics.During heating cycles, the power consumption changes mainly due to transients of the discharge pressure, which in turn dominates the pressure ratio.The lower power consumption of defrost cycles can be attributed to the lower compressor speed (Figure 8) and the smaller pressure ratio.Predictions of the compressor power correspond well with the measurements since the model accurately captures refrigerant pressures.flow direction in heating mode, with CV1 representing the inlet connected to the EXV and CV30 representing the exit connected to the reversing valve.A significant amount of frost is formed during the 2-hour frosting cycle.Clearly more frost accumulates close to the liquid-line of the coil, where the refrigerant temperature is lower at the beginning of each frosting cycle.Correspondingly, the frost tends to be denser close to the vapor-line.During defrost cycles, the frost is considered to be completely melted when the thickness is smaller than 0.01 mm.Then the frost density is reset to 25 kg/m 3 to provide a reasonable initial condition for the next frosting cycle.Figure 13 shows the simulated free airflow passage of the outdoor coil due to frost blockage.It can be noted that 23% of the total flow passage is blocked after two hours of frosting operation.Relatively large values of NER in the indoor air-side capacity and compressor power can be attributed to prediction errors in the refrigerant mass flow rate associated with performance maps.However, results of NER all lie below 0.025, which are sufficiently good for transient simulations.Additionally, predictions of the evolution of frost properties over time provides insights into the non-uniform frost formation and melting phenomena, which are typically difficult to measure and characterize experimentally at a system level.The proposed models are attractive for development and evaluation of improved control and fault detection designs for ASHPs with frosting and defrosting taken into considerations.interface remains at 273.15 K.The governing equations of this stage are (41), ( 43), ( 46)-( 50) where T water,s,i is the water layer surface temperature that appears in this stage as an intermediate variable to quantify heat dissipation to the ambient, ρ amb and ρ water,s,i are water vapor densities evaluated at the ambient temperature and the water layer surface temperature, respectively.c v is the vaporization coefficient (Qiao, Aute, and Radermacher 2018).
In the dry heating stage, heat is directly dissipated from the metal wall to the ambient air.Variables of the thicknesses and temperatures in the frost melting model are inactive.

Figure 1 .
Figure 1.Staggered grid for discretization of balance equations.

Figure 4 .
Figure 4. Fuzzy numbers and membership functions for linguistic variable T w,i − T melt .

Figure 5 .
Figure 5. Model representation of a reversing valve.
Figure 6 depicts the cross-flow heat exchanger model composed of flow and heat&mass transfer elements.The equivalent refrigerant flow length is divided into n equally sized control volumes (CVs).Of each discretized length,

Figure 6 .
Figure 6.Segments of finite-volume heat exchanger model with frost incorporated.

Figure 7
Figure 7 reports validations of the refrigerant discharge

Figure 7 .
Figure 7. Validations of the refrigerant discharge and suction pressures.

Figure 10 .
Figure 10.Validations of the indoor unit air-side capacity and compressor power.

Figure 11 and
Figure11and Figure12show the simulated frost thickness and density trajectories associated with selected control volumes.Note that the CVs are indexed following the

Figure 13 .
Figure 13.Simulated free airflow passage of the outdoor coil.

Figure 14 .
Figure 14.Normalized error residuals (NER) for: discharge pressure p dis , suction pressure p suc , indoor air-side capacity Qa and compressor power P c .