Numerical Analysis for Performance Evaluation of a Multi-Functional CO2 Heat Pump Water Heating System

In recent years, CO2 heat pump water heating systems have been developed, and their performances have been enhanced while their functions have been expanded. In a multi-functional system, used for both hot water supply and bath heating, hot water retrieved from the top of a storage tank is returned to its bottom or side after heat exchange for bath heating, which destroys the stratified temperature distribution in the storage tank and degrades the system performance. In this paper, the performance of a multi-functional CO2 heat pump water heating system has been evaluated by numerical simulation. A system model was created by combining component models for a CO2 heat pump, mixing valves, a storage tank, a heat exchanger, and a bathtub. Partly, perfect mixing by hot water return was assumed in the component model for the storage tank, and its validity was verified through experiments. A performance analysis has been conducted under daily repeated hot water and bath heating demands, and the system performance was evaluated at a periodically steady state. As a result, the system efficiency and the volume of unused hot water in the multi-functional system decreased by 4.9% and 16.3%, respectively, as compared to those in the uni-functional system, when hot water returned to the bottom of the storage tank. When the position for hot water return is heightened, the system efficiency becomes higher than that in the uni-functional system, while the volume of unused hot water decreases furthermore.


Introduction
In Japan, the energy consumption in the residential sector accounts for about 16% of the total one in all the sectors, and the energy consumption for hot water supply accounts for about 28% of that in the residential sector [1]. This is because hot water is used not only for washing, cooking, and showering, but also for baths. Thus, saving energy has been important for hot water supply in the residential sector. Under this situation, water heating systems (WHSs), each of which is composed of a heat pump (HP) using CO 2 as a natural refrigerant and a hot water storage tank, have been developed and commercialized widely [2]. In addition, their performances have been enhanced, while their functions have been expanded. Although HP WHSs have been introduced in other countries, their refrigerants are not CO 2 , but R410a and R134a. In addition, since the hot water consumption is smaller, sizes of storage tanks are also smaller. It seems that the situation in Japan is not true for other countries.
In earlier years, attention has been paid to analyzing and enhancing the performance of a CO 2 HP as the main component of a WHS. Some review papers on the performance analysis of CO 2 enhance the performance by extracting tepid water from a storage tank [36]. In addition, operation and design optimizations have been tried: Yokoyama et al. have estimated system performance values and have clarified the optimal operating conditions under daily change in simulated hot water demand [37]; Deng et al. have analyzed and optimized the performance of a combination of a solar heater and an HP WHS [38]. These results have been obtained for uni-functional systems with the function only of hot water supply. This is because hot water is retrieved from the top of the storage tank, and water is fed to its bottom, which keeps the temperature distribution in the storage tank stratified because of the dependence of water density on temperature, and thus, the stratified temperature distribution is expressed accurately by a one-dimensional model.
On the other hand, multi-functional CO 2 HP WHSs with the functions not only of hot water supply but also of bath heating have also been developed. In addition, the Japanese Industrial Standard "Residential Heat Pump Water Heaters" has been established, and the method of testing the performance of systems with the functions of hot water supply and bath heating has been prescribed [39]. It is necessary to distinguish these two functions for the performance analysis of a multi-functional system. This is because hot water retrieved from the top of the storage tank is returned to its bottom or side after heat exchange for bath heating, which destroys a stratified temperature distribution in the storage tank and causes three-dimensional convectional water flow, and the temperature distribution in the storage tank is essentially different from that of a uni-functional system. As a result, system performance values of the multi-functional system differ from those of the uni-functional system. Therefore, it is important to conduct the performance analysis of the multi-functional system in consideration of characteristics of the two functions.
In this paper, the performance of a multi-functional CO 2 HP WHS was analyzed by numerical simulation. First, a conventional component model for a storage tank was revised by assuming partly perfect mixing by hot water return based on the results obtained by experiments, which have been conducted to grasp the behavior of temperature distribution in the storage tank, in case that hot water retrieved from the top of the storage tank returned to its bottom or side after heat exchange for bath heating. Second, component models for a heat exchanger and a bathtub were newly created, and these models were combined with those for the other components to form that for the system. Finally, a performance analysis was conducted under daily repeated hot water and bath heating demands, and the system performance was evaluated at a periodically steady state after several days. The hourly changes in the temperature distribution in the storage tank and the HP coefficient of performance (COP) were investigated, and the daily HP COP, the storage and system efficiencies, and the volumes of stored and unused hot water were evaluated. The difference in the performance between uni-functional and multi-functional systems was also investigated. In addition, the influence of the position for hot water return on the system performance was examined. Figure 1 shows the configuration of the multi-functional CO 2 HP WHS with the functions of hot water supply and bath heating investigated in this paper. By the former function, the system supplies hot water for washing, cooking, shower and bath. By the latter function, the system keeps the temperature of hot water supplied to the bathtub. This is because hot water supplied to the bathtub once is usually used to family members who bath intermittently in Japan.

CO 2 HP WHS
This system is composed of a CO 2 HP, a hot water storage tank, mixing valves, a heat exchanger, a bathtub, pipes, and pumps. In the charging mode, the system heats water using the refrigeration cycle of the CO 2 HP and stores hot water in the storage tank. In the tapping mode, hot water stored in the storage tank is retrieved from the top of the storage tank, and is mixed with feed water. It is used for general and bath hot water supply. In addition, hot water stored in the storage tank is retrieved from the top of the storage tank, and is used to exchange heat for bath heating. The hot water is returned to the bottom or side of the storage tank after heat exchange.

Component Models
To conduct a numerical simulation, models of components which constitute the WHS were created. Mathematical equations of the component models are given in appendix A. Only brief descriptions of the component models are given here.
The component models for the CO2 HP and mixing valve were conventional ones which have been proposed previously [34,40]. A simplified static component model was adopted and expressed by algebraic equations for the CO2 HP. Mass and energy balance relationships, energy input-output relationship, and approximate relationships for the power consumption and HP COP were basic equations to be satisfied. A static component model was also adopted and expressed by algebraic equations for the mixing valve. Mass and energy balance relationships were basic equations to be satisfied.
The component model for the storage tank was revised from the one which has been proposed previously [31,34,40]. A one-dimensional dynamic component model was adopted and expressed by differential algebraic equations for the storage tank. It was revised partly for the control volume where hot water is returned by bath heating and the region where a stratified temperature distribution is destroyed by hot water return. Mass and energy balance relationships for the control volumes, which were created by dividing the storage tank vertically, were basic equations to be satisfied. Heat transfer by water flow was related with mass flow rates and temperatures. Outlet water temperatures were equalized with the temperatures at the corresponding positions in the storage tank. On the other hand, when hot water retrieved from the top of the storage tank is returned to its bottom or side after heat exchange for bath heating, a stratified temperature distribution is destroyed, and three-dimensional convectional water flow arises. To consider such a phenomenon, it was necessary to conduct a three-dimensional thermo-fluid numerical analysis. However, this takes a long computation time, and it does not meet the purpose of this paper. Here, a partly perfect mixing model was adopted for the region where a stratified temperature distribution was destroyed and three-dimensional convectional water flow arises. Namely, the region where the temperatures in lower control volumes become higher than those in higher ones is identified; it is assumed that water in the region was perfectly mixed, and that the temperature of water in the region became uniform instantaneously.
The component models for the heat exchanger and bathtub were created newly to consider the function of bath heating. The counter flow type was selected for the heat exchanger. Since the capacity of the heat exchanger was much smaller than that of the storage tank, a static component model was

Component Models
To conduct a numerical simulation, models of components which constitute the WHS were created. Mathematical equations of the component models are given in Appendix A. Only brief descriptions of the component models are given here.
The component models for the CO 2 HP and mixing valve were conventional ones which have been proposed previously [34,40]. A simplified static component model was adopted and expressed by algebraic equations for the CO 2 HP. Mass and energy balance relationships, energy input-output relationship, and approximate relationships for the power consumption and HP COP were basic equations to be satisfied. A static component model was also adopted and expressed by algebraic equations for the mixing valve. Mass and energy balance relationships were basic equations to be satisfied.
The component model for the storage tank was revised from the one which has been proposed previously [31,34,40]. A one-dimensional dynamic component model was adopted and expressed by differential algebraic equations for the storage tank. It was revised partly for the control volume where hot water is returned by bath heating and the region where a stratified temperature distribution is destroyed by hot water return. Mass and energy balance relationships for the control volumes, which were created by dividing the storage tank vertically, were basic equations to be satisfied. Heat transfer by water flow was related with mass flow rates and temperatures. Outlet water temperatures were equalized with the temperatures at the corresponding positions in the storage tank. On the other hand, when hot water retrieved from the top of the storage tank is returned to its bottom or side after heat exchange for bath heating, a stratified temperature distribution is destroyed, and three-dimensional convectional water flow arises. To consider such a phenomenon, it was necessary to conduct a three-dimensional thermo-fluid numerical analysis. However, this takes a long computation time, and it does not meet the purpose of this paper. Here, a partly perfect mixing model was adopted for the region where a stratified temperature distribution was destroyed and three-dimensional convectional water flow arises. Namely, the region where the temperatures in lower control volumes become higher than those in higher ones is identified; it is assumed that water in the region was perfectly mixed, and that the temperature of water in the region became uniform instantaneously.
The component models for the heat exchanger and bathtub were created newly to consider the function of bath heating. The counter flow type was selected for the heat exchanger. Since the capacity of the heat exchanger was much smaller than that of the storage tank, a static component model was adopted and expressed by algebraic equations. Mass and energy balance relationships for the flows with higher and lower temperatures, which correspond to the storage tank and bathtub sides, respectively, were basic equations to be satisfied. The exchanged heat was expressed as the product of overall heat transfer coefficient, heat transfer area, and log mean temperature difference. Since the capacity of the bathtub was as large as that of the storage tank, a dynamic component model was adopted and expressed by differential algebraic equations for it. Here, perfect mixing was assumed in the bathtub. Mass and energy balance relationships were basic equations to be satisfied. Outlet water temperatures were equalized with the temperature in the bathtub.

System Model
A model for the WHS was created by the aforementioned component models for the CO 2 HP, mixing valve, storage tank, heat exchanger, bathtub, a substance model for water, and other conditions.
The substance model gives properties of water which are commonly used in all the components. The other conditions were as follows. The connection conditions give equations which equalize the values of variables at the connection points between components, i.e., mass flow rates and temperatures between the CO 2 HP, mixing valves, storage tank, heat exchanger, and bathtub. The boundary conditions give the values variables at the boundaries of the system, i.e., the feed water temperature, the mass flow rate and temperature of hot water to the tapping site for both general and bath use, and the mass flow rate of hot water disposed of from the bathtub. The control conditions give the values of variables to be kept constant, i.e., the outlet water temperature of the CO 2 HP during its operation, and the mass flow rates of the heat exchanger on both the storage tank and bathtub sides for bath heating. The ambient condition gives the value of a parameter to be changed by the environment, i.e., the air temperature. The initial conditions give the values of variables at an initial time, i.e., the temperatures of water in the storage tank and bathtub, and the mass of water in the bathtub.

Numerical Solution
The aforementioned system model is expressed by the following set of nonlinear differential algebraic equations: f x(t), .
x(t), y(t), t = 0 where f is the vector for all the differential algebraic equations to be satisfied; x is the vector for variables with their derivatives, i.e. the temperatures of water in the storage tank and bathtub; .
x is the derivative of x; y is the vector for variables without their derivatives, i.e., all the other variables excluding the temperatures of water in the storage tank and bathtub; t 0 is the initial time; and .
x 0 is the initial values of x.
This set of equations was solved by combining numerical solution algorithms for nonlinear differential and algebraic equations hierarchically. At the upper level, the values of x(t) after the sampling time interval were calculated by solving nonlinear differential equations using the Runge-Kutta method at each sampling time. At the lower level, the values of .
x(t) and y(t) were calculated by solving nonlinear algebraic equations using the Newton-Raphson method. However, Equation (1) did not include the partly perfect mixing model of the storage tank. Thus, the partly perfect mixing model was incorporated into the numerical solution as shown in Figure 2. At each sampling time, the values of variables after the sampling time interval were obtained by the solution method without consideration of partly perfect mixing in the storage tank. Then, the region where the perfect mixing arised was identified, and the temperature distribution in the region was modified to form a uniform distribution.

Experimental Setup
The validity of the aforementioned numerical solution, which incorporates the partly perfect mixing model of the storage tank, should be investigated. Here, some experiments were conducted under the conditions with hot water return in the tapping mode, and simulated and experimental results are compared.
The schematic of the experimental setup is shown in Figure 3. The storage tank model was made of acryl, and its specification is shown in Table 1 [36]. The 5th inlet port from the top on the side, and the outlet port on the top were used for the experiments. Water temperatures in the storage tank were measured by a series of type-T thermocouples placed vertically with a span of 0.10 m. It was close to the inside wall of the tank, apart from the inlet ports on the side by a circumferential angle of 90°. Its measurement accuracy was ±1.0 °C. Water temperature in each port was measured by a type-K thermocouple inserted there. Its measurement accuracy was ±2.5 °C. These type-T and type-K thermocouples were shown by red and blue points, respectively, in Figure 3. The flow rates of water were adjusted by the valve and mass flow sensor. The measurement accuracy of the mass flow sensor was within ±3.0% of the full scale volumetric flow rate 0.067 L/s, or within ±2.0 × 10 −3 L/s, and its reproducibility accuracy was within ±0.5% of the same full scale volumetric flow rate, or within ±0.33 × 10 −3 L/s.

Experimental Setup
The validity of the aforementioned numerical solution, which incorporates the partly perfect mixing model of the storage tank, should be investigated. Here, some experiments were conducted under the conditions with hot water return in the tapping mode, and simulated and experimental results are compared.
The schematic of the experimental setup is shown in Figure 3. The storage tank model was made of acryl, and its specification is shown in Table 1 [36]. The 5th inlet port from the top on the side, and the outlet port on the top were used for the experiments. Water temperatures in the storage tank were measured by a series of type-T thermocouples placed vertically with a span of 0.10 m. It was close to the inside wall of the tank, apart from the inlet ports on the side by a circumferential angle of 90 • . Its measurement accuracy was ±1.0 • C. Water temperature in each port was measured by a type-K thermocouple inserted there. Its measurement accuracy was ±2.5 • C. These type-T and type-K thermocouples were shown by red and blue points, respectively, in Figure 3. The flow rates of water were adjusted by the valve and mass flow sensor. The measurement accuracy of the mass flow sensor was within ±3.0% of the full scale volumetric flow rate 0.067 L/s, or within ±2.0 × 10 −3 L/s, and its reproducibility accuracy was within ±0.5% of the same full scale volumetric flow rate, or within ±0.33 × 10 −3 L/s.    Numerical simulations were also conducted under the same conditions with those for the experiments. The number of control volumes for the storage tank model was set at 200, and the sampling time interval for the Runge-Kutta method was set at 10 s.

Case A
In case A, hot water was supplied to the inlet port on the side, and water was extracted from the outlet port on the top. The initial temperature distribution in tank 1 was realized by pouring city water with a temperature of about 15 • C into tank 1 through a piping connection between P 1 and P 2 in Figure 3. Although this temperature distribution does not arise during bath heating in an existing system, it was suitable for the purpose of investigating the change in the temperature distribution in an extreme case that hot water is supplied into water with a large temperature difference. Hot water in tank 2 heated up to about 65 • C by the electric heater, which supplied to the inlet port through a piping connection between P 1 and P 3 in Figure 3. Water extracted from the outlet port was exhausted into tank 3. The volumetric flow rates of hot water and water were adjusted by the valves and mass flow sensors to be set at 0.0305 L/s. The pump was operated for 400 s from the initial time. Figure 4 compares experimental and simulated results in terms of the hourly change in the temperature distribution in tank 1. According to both results, during the pump operation, the temperature distribution below the inlet port hardly changed, while that above the inlet port changed due to hot water return, and the temperature difference across the inlet port increased drastically. In addition, the gradient of temperature distribution near the inlet port remained large. This is because water mixed with hot water at the inlet port goes up by both convection and water flow. According to the experimental results, although the temperature distribution above the inlet port was not uniform, its gradient was small. According to the simulated results, on the other hand, the temperature distribution above the port was uniform because of the assumption of the partly perfect mixing. However, the simulated results coincide well with the experimental ones in terms of the aforementioned features of temperature distribution. After the pump operation, it was expected that the temperature difference across the inlet port would decrease quickly because of a large heat conduction, and that the simulated results coincide better with the experimental ones. Thus, the assumption of the partly perfect mixing was validated in the case that hot water was supplied into water. It is difficult to investigate with this experimental setup why there exists a slight gradient in the temperature distribution and how large the gradient is. It is necessary to understand the convection behavior for a short period by a three-dimensional analysis by experiment or simulation. Here, the performance of the overall system for a long period is focused on, and the detail of the convention behavior is not considered.
conduction, and that the simulated results coincide better with the experimental ones. Thus, the assumption of the partly perfect mixing was validated in the case that hot water was supplied into water. It is difficult to investigate with this experimental setup why there exists a slight gradient in the temperature distribution and how large the gradient is. It is necessary to understand the convection behavior for a short period by a three-dimensional analysis by experiment or simulation. Here, the performance of the overall system for a long period is focused on, and the detail of the convention behavior is not considered.

Case B
In case B, water was supplied to the inlet port on the side, and hot water was extracted from the outlet port on the top. The initial temperature distribution in tank 1 was realized by pouring hot water in tank 2 heated up to about 60 °C by the electric heater into tank 1 through a piping connection between P1 and P2, as seen in Figure 3. City water with a temperature of about 25 °C was supplied to the inlet port through a piping connection between P1 and P3, as seen in Figure 3. Although this situation does not arise during bath heating in an existing system, it was suitable for the purpose of investigating the change in the temperature distribution in the extreme case that water is supplied into hot water with a large temperature difference. Hot water extracted from the outlet port was exhausted into tank 3. The volumetric flow rates of water and hot water were adjusted by the valves and mass flow sensors to be set at 0.0305 L/s. The pump was operated for 400 s from the initial time. Figure 5 compares experimental and simulated results in terms of the hourly change in the temperature distribution in tank 1. According to both the results, during the pump operation, the temperature distribution above the inlet port hardly changed, except near the inlet port; while the area below the inlet port changed due to water return, the temperature difference across the inlet port

Case B
In case B, water was supplied to the inlet port on the side, and hot water was extracted from the outlet port on the top. The initial temperature distribution in tank 1 was realized by pouring hot water in tank 2 heated up to about 60 • C by the electric heater into tank 1 through a piping connection between P 1 and P 2 , as seen in Figure 3. City water with a temperature of about 25 • C was supplied to the inlet port through a piping connection between P 1 and P 3 , as seen in Figure 3. Although this situation does not arise during bath heating in an existing system, it was suitable for the purpose of investigating the change in the temperature distribution in the extreme case that water is supplied into hot water with a large temperature difference. Hot water extracted from the outlet port was exhausted into tank 3. The volumetric flow rates of water and hot water were adjusted by the valves and mass flow sensors to be set at 0.0305 L/s. The pump was operated for 400 s from the initial time. Figure 5 compares experimental and simulated results in terms of the hourly change in the temperature distribution in tank 1. According to both the results, during the pump operation, the temperature distribution above the inlet port hardly changed, except near the inlet port; while the area below the inlet port changed due to water return, the temperature difference across the inlet port increased drastically. However, the gradient of temperature distribution above and near the inlet port was smaller than that in case A. This is because hot water mixed with water at the inlet port goes down by convection but goes up by water flow. According to the experimental results, although the temperature distribution below the port was not uniform, its gradient was small. According to the simulated results, on the other hand, the temperature distribution below the port was uniform because of the assumption of the partly perfect mixing. However, the simulated results coincide well with the experimental ones in terms of the aforementioned features of temperature distribution. After the pump operation, it was expected that the temperature difference across the inlet port would decrease quickly because of a large heat conduction, and that the simulated results would coincide better with the experimental ones. Thus, the assumption of the partly perfect mixing was validated in the case that water was supplied into hot water. For the same reason, the detail of the convention behavior has not been considered here. increased drastically. However, the gradient of temperature distribution above and near the inlet port was smaller than that in case A. This is because hot water mixed with water at the inlet port goes down by convection but goes up by water flow. According to the experimental results, although the temperature distribution below the port was not uniform, its gradient was small. According to the simulated results, on the other hand, the temperature distribution below the port was uniform because of the assumption of the partly perfect mixing. However, the simulated results coincide well with the experimental ones in terms of the aforementioned features of temperature distribution. After the pump operation, it was expected that the temperature difference across the inlet port would decrease quickly because of a large heat conduction, and that the simulated results would coincide better with the experimental ones. Thus, the assumption of the partly perfect mixing was validated in the case that water was supplied into hot water. For the same reason, the detail of the convention behavior has not been considered here.

Conditions
To investigate the influence of bath heating on the system performance, the following two

Conditions
To investigate the influence of bath heating on the system performance, the following two systems were investigated: one is the multi-functional system with the functions of both hot water supply and bath heating (system M); the other is the uni-functional system with the function only of hot water supply (system U). Table 2 shows the specifications of the multi-functional CO 2 HP WHS investigated here. Their values are given based on an existing system. Figure 6 shows values measured using an existing system and approximate functions identified based on measured values for the power consumption, HP COP, and their product, or heat output of the CO 2 HP [34]. As an example, these are shown in relation only to the inlet water temperature for constant air and outlet water temperatures of 16 and 85 • C, respectively. Here, each value is relative to its rated one for the air and inlet/outlet water temperatures of 16, 17, and 65 • C, respectively. The reason why the performance of the CO 2 HP is expressed by relative values is as follows: the performance of the CO 2 HP was measured for an existing system; however, the performance has been enhanced significantly after that because components such as compressors, gas coolers, etc. have been improved technologically; thus, absolute values for the performance obtained by the measurement may mislead one to the performance of CO 2 HPs in the current market, and hence, they are not suitable. For the same reason, the performance of the WHS is also expressed by relative values.  The air temperature was set at 16 °C as the ambient condition, while the feed water temperature was set at 17 °C as one of the boundary conditions. These are prescribed for mid-season by the Japanese Industrial Standards [39]. Hot water and bath heating demands shown in Figure 7 were used, which were also prescribed as standard on the representative day by the Japanese Industrial Standards. Figure 7a shows the hourly change in the volumetric flow rate of the demand for general and bath hot water supply. The temperature of hot water supplied to the tapping site was set at 40 °C. Figure 7b shows the hourly change in the heat of the demand for bath heating. Here, the height of each vertical line corresponds to the total heat. Since the temperature of hot water at the top of the The air temperature was set at 16 • C as the ambient condition, while the feed water temperature was set at 17 • C as one of the boundary conditions. These are prescribed for mid-season by the Japanese Industrial Standards [39]. Hot water and bath heating demands shown in Figure 7 were used, which were also prescribed as standard on the representative day by the Japanese Industrial Standards. Figure 7a shows the hourly change in the volumetric flow rate of the demand for general and bath hot water supply. The temperature of hot water supplied to the tapping site was set at 40 • C. Figure 7b shows the hourly change in the heat of the demand for bath heating. Here, the height of each vertical line corresponds to the total heat. Since the temperature of hot water at the top of the storage tank changed with time, the temperature of hot water and duration could not be prescribed. Thus, when the total heat attained its prescribed value, the bath heating was considered to be finished. In addition, the volumetric flow rates of hot water with higher temperature on the storage tank side and lower temperature on the bathtub side were set at 0.06 and 0.15 L/s, respectively. On the other hand, in system U, the heat equivalent to the total heat was considered as an additional demand for bath hot water supply. Thus, the start time was common to both systems M and U, however, the stop time was independent in each system. storage tank changed with time, the temperature of hot water and duration could not be prescribed. Thus, when the total heat attained its prescribed value, the bath heating was considered to be finished. In addition, the volumetric flow rates of hot water with higher temperature on the storage tank side and lower temperature on the bathtub side were set at 0.06 and 0.15 L/s, respectively. On the other hand, in system U, the heat equivalent to the total heat was considered as an additional demand for bath hot water supply. Thus, the start time was common to both systems M and U, however, the stop time was independent in each system. At the initial time, the temperature of water in the storage tank was set at 17 °C. In addition, an infinitesimal positive value was given in place of zero to the mass of water in the bathtub, so that the set of equations could be solved; its temperature was set at 17 °C. The CO2 HP started up at 2:00 on each day, and shut down when the inlet water temperature attained its prescribed value. Here, the outlet water temperature during HP operation was set at 85 °C, and the inlet water temperature for HP shutdown was set at 30 °C.
The numerical simulation was conducted for eight days so that the system attained a periodically steady state; the system performance was evaluated on the 8th day. The number of control volumes for the storage tank was set at 200. The sampling time interval for the Runge-Kutta method was set at 10 s when there were water flows in the storage tank during HP operation, hot water supply, and bath heating, and was set at 180 s when there is no water flow in the storage tank. At the initial time, the temperature of water in the storage tank was set at 17 • C. In addition, an infinitesimal positive value was given in place of zero to the mass of water in the bathtub, so that the set of equations could be solved; its temperature was set at 17 • C. The CO 2 HP started up at 2:00 on each day, and shut down when the inlet water temperature attained its prescribed value. Here, the outlet water temperature during HP operation was set at 85 • C, and the inlet water temperature for HP shutdown was set at 30 • C.
The numerical simulation was conducted for eight days so that the system attained a periodically steady state; the system performance was evaluated on the 8th day. The number of control volumes for the storage tank was set at 200. The sampling time interval for the Runge-Kutta method was set at 10 s when there were water flows in the storage tank during HP operation, hot water supply, and bath heating, and was set at 180 s when there is no water flow in the storage tank.

Results and Discussion
First, system M was investigated for when hot water returned to the bottom of the storage tank after heat exchange for bath heating. The performance of system M was compared to that of system U. Figures 8 and 9 show the hourly changes in the temperature distribution in the storage tank for systems M and U, respectively. Figures 8a and 9a show the hourly changes during HP operation from 0:00 to 7:00, and Figures 8b and 9b show the hourly changes during hot water supply and bath heating, and hot water supply for systems M and U, respectively, from 7:00 to 24:00. In system U, since only the feed water flows into the bottom of the storage tank, the region with temperatures close to 17 • C expands. In system M, on the other hand, since only hot water demand occurs during 7:00 to 20:00, the temperature distribution at 20:00 was almost the same as that in system U. However, since bath heating demand occurred four times in addition to hot water demand during 20:00 to 22:00, the temperature from the bottom to the middle of the storage tank rose at 22:00. Afterwards, hot water demand occurred, and the temperature at the bottom of the storage tank dropped; it occurred once again from 22:00 to 24:00. As a result, since the temperature at the middle of the storage tank in system M was higher than that in system U, the gradient of temperature distribution in system M was slightly larger than that in system U, thus, the same heat was supplied in both systems.  Figure 10 shows the comparison in the hourly change in the HP COP between systems M and U. Here, each value is relative to its rated one for the air and inlet/outlet water temperatures of 16, 17, and 65 °C, respectively. It shows that the hourly change in the HP COP depends significantly on the temperature distribution from the bottom to the middle of the storage tank at 24:00 after daily hot water supply and bath heating, as shown in Figures 8b and 9b This is because water from the bottom  Figure 10 shows the comparison in the hourly change in the HP COP between systems M and U. Here, each value is relative to its rated one for the air and inlet/outlet water temperatures of 16, 17, and 65 °C, respectively. It shows that the hourly change in the HP COP depends significantly on the temperature distribution from the bottom to the middle of the storage tank at 24:00 after daily hot water supply and bath heating, as shown in Figures 8b and 9b This is because water from the bottom  Figure 10 shows the comparison in the hourly change in the HP COP between systems M and U. Here, each value is relative to its rated one for the air and inlet/outlet water temperatures of 16, 17, and 65 • C, respectively. It shows that the hourly change in the HP COP depends significantly on the temperature distribution from the bottom to the middle of the storage tank at 24:00 after daily hot water supply and bath heating, as shown in Figures 8b and 9b This is because water from the bottom to the middle of the storage tank enters the CO 2 HP during HP operation, and the performance of the CO 2 HP depends significantly on the inlet water temperature, as shown in Figure 6. Since system U had an almost uniform temperature distribution at the lowest quarter part of the storage tank, it had almost constant HP COP during the first half of HP operation. On the other hand, since system M had an increase in the temperature at the middle of the storage tank, it had an almost constant (although still lower) HP COP during the middle of HP operation. This difference resulted in a lower daily HP COP in system M as compared toto that in system U.  Table 3 shows the comparison in the following daily system performance values between systems M and U: HP COP, storage and system efficiencies, and volumes of stored and unused hot water. Here, the HP COP is defined as the ratio of the daily heat output to daily power consumption, the storage efficiency is defined as the ratio of daily heat supply to daily heat output, and the system efficiency is defined as the ratio of daily heat supply to daily power consumption, which is equal to the product of the HP COP and storage efficiency. The HP COP, and the storage and system efficiencies in system M are relative to those in system U. The volumes of stored and unused hot water are defined as the ones at 7:00 after HP operation and 24:00 after hot water supply and bath heating, respectively. They were evaluated as the volumes of hot water with a temperature of 40 °C obtained by mixing the hot water with temperatures higher than 40 °C and the feed water of 17 °C.
Since the temperature at the middle of the storage tank rose in system M, the HP COP in system M decreases by 4.7% as compared toto that in system U. In addition, since the average temperature in the storage tank in system M becomes higher than that in system U, the storage efficiency in system M also decreased as compared toto that in system U, but its decrease rate remains only 0.2%. As a whole, the system efficiency in system M decreases by 4.9% as compared to that in system U. Because of the difference in the temperature distribution in the storage tank between systems M and U, the volume of stored hot water in system M increased and that of unused hot water in system M decreased as compared to those in system U. Generally, there is a trade-off relationship between the system efficiency and the volume of unused hot water, depending on operating conditions in a system. However, both the system efficiency and the volume of unused hot water decreased in system M as compared to those in system U.   Table 3 shows the comparison in the following daily system performance values between systems M and U: HP COP, storage and system efficiencies, and volumes of stored and unused hot water. Here, the HP COP is defined as the ratio of the daily heat output to daily power consumption, the storage efficiency is defined as the ratio of daily heat supply to daily heat output, and the system efficiency is defined as the ratio of daily heat supply to daily power consumption, which is equal to the product of the HP COP and storage efficiency. The HP COP, and the storage and system efficiencies in system M are relative to those in system U. The volumes of stored and unused hot water are defined as the ones at 7:00 after HP operation and 24:00 after hot water supply and bath heating, respectively. They were evaluated as the volumes of hot water with a temperature of 40 • C obtained by mixing the hot water with temperatures higher than 40 • C and the feed water of 17 • C. Since the temperature at the middle of the storage tank rose in system M, the HP COP in system M decreases by 4.7% as compared toto that in system U. In addition, since the average temperature in the storage tank in system M becomes higher than that in system U, the storage efficiency in system M also decreased as compared toto that in system U, but its decrease rate remains only 0.2%. As a whole, the system efficiency in system M decreases by 4.9% as compared to that in system U. Because of the difference in the temperature distribution in the storage tank between systems M and U, the volume of stored hot water in system M increased and that of unused hot water in system M decreased as compared to those in system U. Generally, there is a trade-off relationship between the system efficiency and the volume of unused hot water, depending on operating conditions in a system. However, both the system efficiency and the volume of unused hot water decreased in system M as compared to those in system U. Second, system M was investigated when hot water returned to the side of the storage tank after heat exchange for bath heating. The side position was changed a parameter, and its influence on the system performance values was examined. Figure 11 shows the temperature distributions in the storage tank for systems M at 7:00 and 24:00 for five different control volumes for hot water return. The temperature distributions for the 160th control volume for hot water return were qualitatively similar to those for the 200th one (bottom). This is because the temperature, which was vertically constant by partly perfect mixing, was below 30 • C at 24:00. Additionally, the corresponding part disappeared at 7:00 by the inlet water temperature for HP shutdown of 30 • C. However, the temperature distributions for the 120th control volume for hot water return were qualitatively different from those for the 160th and 200th ones. This is because the temperature, which was vertically constant by partly perfect mixing, was above 30 • C at 24:00. It remained at 7:00, because of the inlet water temperature for HP shutdown of 30 • C. Thus, the volume of hot water stored at 7:00 was much smaller, and consequently, the volume of hot water unused at 24:00 was almost zero. The temperature distributions for the 40th and 80th control volumes for hot water return became close to those for the 160th and 200th ones. This is because the position for hot water return was heightened, and the temperature distributions at the lower part were not affected by the hot water return. is because the temperature, which was vertically constant by partly perfect mixing, was above 30 °C at 24:00. It remained at 7:00, because of the inlet water temperature for HP shutdown of 30 °C. Thus, the volume of hot water stored at 7:00 was much smaller, and consequently, the volume of hot water unused at 24:00 was almost zero. The temperature distributions for the 40th and 80th control volumes for hot water return became close to those for the 160th and 200th ones. This is because the position for hot water return was heightened, and the temperature distributions at the lower part were not affected by the hot water return. Figure 11. Influence of position for hot water return on temperature distributions. Figure 12 shows the system performance values for nine different control volumes for hot water return. Figure 12a shows the HP COP and storage and system efficiencies, and Figure 12b shows the volumes of stored and unused hot water. Each value is relative to that for the 200th control volume for hot water return. Since the deficit in hot water supply arised for the 140th control volume for hot water return, no data is shown in these figures. As aforementioned in relation to Figure 11, the volumes of stored and unused hot water decreased for the 120th control volume for hot water return. They increased for higher and lower control volumes for hot water return. The volumes of stored and unused hot water for lower control volumes for hot water return were slightly larger than those for higher ones. On the other hand, the storage efficiency increased for the 120th control volume for hot water return. This is because the average temperature in the storage tank was low. It decreased for higher and lower control volumes for hot water return. However, the storage efficiency for higher control volumes for hot water return was larger than that for lower ones. In addition, the HP COP for higher control volumes for hot water return was much larger than that for lower ones. This is because the temperature distributions at the lower part were not affected by hot water return. Thus, the system efficiency increases for the 120th control volume for hot water return. Although it decreased for higher and lower control volumes for hot water return, the system efficiency for higher control volumes was much larger than that for lower ones.  Figure 12 shows the system performance values for nine different control volumes for hot water return. Figure 12a shows the HP COP and storage and system efficiencies, and Figure 12b shows the volumes of stored and unused hot water. Each value is relative to that for the 200th control volume for hot water return. Since the deficit in hot water supply arised for the 140th control volume for hot water return, no data is shown in these figures. As aforementioned in relation to Figure 11, the volumes of stored and unused hot water decreased for the 120th control volume for hot water return. They increased for higher and lower control volumes for hot water return. The volumes of stored and unused hot water for lower control volumes for hot water return were slightly larger than those for higher ones. On the other hand, the storage efficiency increased for the 120th control volume for hot water return. This is because the average temperature in the storage tank was low. It decreased for higher and lower control volumes for hot water return. However, the storage efficiency for higher control volumes for hot water return was larger than that for lower ones. In addition, the HP COP for higher control volumes for hot water return was much larger than that for lower ones. This is because the temperature distributions at the lower part were not affected by hot water return. Thus, the system efficiency increases for the 120th control volume for hot water return. Although it decreased for higher and lower control volumes for hot water return, the system efficiency for higher control volumes was much larger than that for lower ones.
volumes of stored and unused hot water decreased for the 120th control volume for hot water return. They increased for higher and lower control volumes for hot water return. The volumes of stored and unused hot water for lower control volumes for hot water return were slightly larger than those for higher ones. On the other hand, the storage efficiency increased for the 120th control volume for hot water return. This is because the average temperature in the storage tank was low. It decreased for higher and lower control volumes for hot water return. However, the storage efficiency for higher control volumes for hot water return was larger than that for lower ones. In addition, the HP COP for higher control volumes for hot water return was much larger than that for lower ones. This is because the temperature distributions at the lower part were not affected by hot water return. Thus, the system efficiency increases for the 120th control volume for hot water return. Although it decreased for higher and lower control volumes for hot water return, the system efficiency for higher control volumes was much larger than that for lower ones.

Conclusions
The performance of a multi-functional CO 2 HP WHS with the functions of both hot water supply and bath heating has been analyzed by numerical simulation. First, a component model for a storage tank were revised by assuming partly perfect mixing based on the experimental results, in case that hot water retrieved from the top of the storage tank is returned to the bottom or side of the storage tank after heat exchange for bath heating. Second, component models for a heat exchanger and a bathtub have been newly created, and these models have been combined with those for the other components to form that for the system. Finally, a performance analysis has been conducted under daily repeated hot water and bath heating demands, and the system performance has been evaluated at a periodically steady state after several days. The hourly changes in the temperature distribution in the storage tank and HP COP, as well as the daily system performance values, have been investigated in comparison with those of a uni-functional system with the function only of hot water supply. In addition, the influence of the position for hot water return on the system performance has been examined. The following main results have been obtained:

•
In the multi-functional system, the temperature at the middle of the storage tank rises when hot water is returned to the bottom of the storage tank. As the position for hot water return is heightened, the part with a constant temperature at the middle of the storage tank is also heightened.

•
The multi-functional system has almost constant (although lower) HP COP during the middle of HP operation when hot water returned to the bottom of the storage tank. As a result, the daily HP COP in the multi-functional system decreased by 4.7% as compared to that in the uni-functional system. However, this decrease recovered when the position for hot water return heightened.

•
The average temperature in the storage tank in the multi-functional system became higher than that in the uni-functional system, when hot water returned to the bottom of the storage tank. Thus, the storage efficiency in the multi-functional system also decreased as compared to that in the uni-functional system, although its decrease rate remained only 0.2%. Since the average temperature in the storage tank decreased when the position for hot water return heightened, the storage efficiency became higher than that in the uni-functional system.
• As a whole, the system efficiency in the multi-functional system decreased by 4.9% as compared to that in the uni-functional system, when hot water returned to the bottom of the storage tank. When the position for hot water return heightened, the system efficiency also became higher than that in the uni-functional system. • The volume of unused hot water in the multi-functional system decreased by 16.3% as compared to that in the uni-functional system, when hot water returned to the bottom of the storage tank. Furthermore, when the position for hot water return heightened, the volume of unused hot water decreased. • When the constant temperature in the storage tank by the hot water return was slightly above the inlet water temperature for HP shutdown, unique temperature distributions arised. The storage and system efficiencies increased drastically, while the volume of unused hot water decreased drastically. Thus, a deficit in hot water supply may arise.
The validity of the performance, especially the temperature distribution of the storage tank, of the uni-functional system evaluated by numerical simulation was investigated by experiment in the previous paper [34]. However, the validity of the performance of the multi-functional system evaluated by numerical simulation has not been investigated by experiment yet. This should be done as future work.

Conflicts of Interest:
The authors declare no conflict of interest. The funder had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; and in the decision to publish the results. Mathematical equations of the component models are given in this appendix.