Akiko Imura

Dr M R Yeadon

Mechanics of the Fouetté turn

Keywords: Turn, Simulation, Angular Momentum, Ballet

The Fouetté turn in classical ballet is performed repeatedly on one leg with swinging of the free limbs, producing a continued sequence of turns with one turn leading into the next. The purpose of this study was to determine the possible time history profiles of the twisting torque between the supporting leg and the remainder of the body that will allow continued performances of the Fouetté turn. Simulations were performed using a model which comprised the supporting leg and the remainder of the body to find torque profiles that maintain the initial angular velocity so that the state after one revolution is the same as the initial state. The solution space of torque profiles was determined for various rotation times and coefficients of friction between foot and floor. As the time for one revolution became shorter the solution space became smaller and for a given turn time there was a lower limit on the coefficient of friction. As the frictional coefficient became smaller the solution space became smaller and for a given coefficient there was a lower limit on the turn time. Turns of a given tempo can be performed on floors with different friction by modifying the twisting torque profile. When a turn is completed with a net change in angular velocity this can be compensated for in the next turn by adjusting the twisting torque profile.

Introduction

Skilled ballet dancers can continuously perform repeated Fouetté turns (Fig. 1) and while there are various performance styles, such as Italian and Russian Fouetté turns according to the ballet style of the dancer [START_REF] Warren | Classical ballet technique[END_REF], the basic technique is the same. The turn is started from one or two revolutions of the pirouette which is initiated with both feet in contact with the floor to produce the initial angular momentum. The dancer then keeps turning to music, swinging the arms and the free leg while the supporting foot is stationary in full contact with the floor (Figs. 1 A-D).

After the swinging, the dancer adopts the pirouette position during which the foot slips (Fig. 1 E-I) before starting to swing the free leg again (Fig. 1 J-K). The dancer regains the angular momentum lost due to friction during the slipping phase by swinging the free limbs when the foot is stationary, which enables the floor to exert a large frictional torque T F on the foot in the same direction as the swinging [START_REF] Laws | The physics of dance[END_REF][START_REF] Laws | Momentum transfer in dance movement[END_REF][START_REF] Imura | Biomechanics of the continuity and speed change during one revolution of the Fouetté turn[END_REF]. The dancer can keep turning for more than 30 revolutions by repeating these movements.

The frictional torque T F is the only external torque during the Fouetté turn and consequently determines the changes in the angular momentum of the whole body.

This frictional torque is dependent on the limiting frictional torque (limiting T F ) and the twisting torque T used to swing the free limbs. The supporting foot during the Fouetté turn is essentially on tiptoe during the slipping or in full contact with the floor when the foot stops turning. *** insert Fig. 1 here *** Dancers have to perform the Fouetté turn in time to the music in the choreography, facing the front at the same position for the aesthetics of ballet [START_REF] Laws | The physics of dance[END_REF]. However, they sometimes turn to music tempo that is too fast or too slow and struggle with performing successful turns. The friction coefficient between the shoes and the floor of the performance stage may be different to that of the practice studio and this will require technique to be modified.

The purpose of this study was to determine the possible time history profiles of the twisting torque between the supporting leg and the remainder of the body that will allow continued performances of the Fouetté turn. Techniques for coping with changes in tempo and friction were also investigated.

Methods

A computer simulation model of the Fouetté turn was used to investigate the solution space of twisting torque profiles that permitted performances of successful turns for various coefficients of friction between foot and floor and various time periods of turn. The body (mass 49.5 kg) of a typical dancer who participated in a previous study [START_REF] Imura | Biomechanics of the continuity and speed change during one revolution of the Fouetté turn[END_REF] was modeled as two cylinders (Fig. 2): the supporting leg L and the remainder of the body B whose moment of inertia I B about a vertical axis changed according to the positions of the free leg and arms. The time profiles of the foot radius r, the moment of inertia I B , and the normal ground reaction force N were based on experimental data from the study of [START_REF] Imura | Biomechanics of the continuity and speed change during one revolution of the Fouetté turn[END_REF] and were represented by joining adjacent maximum and minimum values using monotonic quintic functions with zero first and second time derivatives at the endpoints (Fig. 3). *** insert Fig. 2 here ***

The radius r was taken to be 0.12 m at the maximum when the supporting foot is fully in contact with the floor and to be 0.05 m at the minimum when the dancer stands on tiptoe (Fig. 3a). These bounding values were determined using the distance between the toe and the center of the pressure (CoP) from the experimental data of [START_REF] Imura | Biomechanics of the continuity and speed change during one revolution of the Fouetté turn[END_REF]. The foot radius time profile was matched to the experimental data with time normalized to the turn time, recognizing that CoP locations were unreliable when the normal reaction force was small.

The time profile of I B (Fig. 3b) was calculated using the theorem of parallel axes by scaling inertia data [START_REF] Ae | Estimation of inertia properties of the body segments in Japanese athletes (in Japanese)[END_REF]) of a subject with similar body mass as the model and using arm and leg positions based upon those exhibited in the experimental study of [START_REF] Imura | Biomechanics of the continuity and speed change during one revolution of the Fouetté turn[END_REF]. The maximum and minimum I B were calculated to be 2.67 and 1.06 kg.m 2 , respectively. The moment of inertia of the leg I L was assumed to be constant during the turn and was calculated to be 0.085 kg.m 2 . The normal ground reaction force N was determined by the following four values based on experimental data from the study of [START_REF] Imura | Biomechanics of the continuity and speed change during one revolution of the Fouetté turn[END_REF]: 2.68 body weights at the time of full foot contact, 0.0 body weights once slipping had started, 1.24 and 0.82 body weights during slipping on tiptoe. It was assumed that the fitted profile was symmetrical about the mid-turn time (Fig. 3c). The average N was one body weight during one turn.

The friction coefficient µ between the shoes and the floor was estimated to be 0.2, calculated from the slipping phase of the experimental data of [START_REF] Imura | Biomechanics of the continuity and speed change during one revolution of the Fouetté turn[END_REF] using the equation = T F / Nr. The leg L was rotated further than B (Fig. 2) by 0.22 radians at the start of the turn according to the experimental data of [START_REF] Imura | Biomechanics of the continuity and speed change during one revolution of the Fouetté turn[END_REF].

Simulations were performed for 1 s starting with the supporting foot in full contact with the floor (Fig. 1A), with the requirement that the dancer rotates 2 in 1 s using an appropriate twisting torque T. Three variables were used to define the time profile of T: T max the initial (positive) torque to swing the free limbs and trunk, T min the (negative) torque which is maintained during the sliding to reverse the swinging, and the time t 1 at which T becomes T min . T max , and T min were joined by monotonic quintic functions, assuming the profile of T to be symmetric since the dancer swings the limbs at the end of the turn as at the start (Fig. 3d). *** insert Fig. 3 here ***

The frictional torque T F while the dancer is slipping can be assumed to be the limiting T F , which is the product of the friction coefficient µ, the normal ground reaction force N, and the radius r of the foot contact area with the floor. When the foot is stationary, the frictional torque T F acting on the supporting foot is equal to T. The frictional torque T F was defined as T F = T when T µ Nr and T F = µ Nr once L slipped. Values for µ , N, r, I B , I L , and T were input into the equations of motion.

The angular momentum of the body B is given by h = I B B φ where B φ is the angle turned by the body B. The torque T applied to the body B by the supporting leg L is equal to the rate of change of h. Thus:

B B B B I I T φ + φ = from which: B B B B I / ) I T ( φ - = φ (1)
The net torque acting on the leg L in the direction of the turn is T F -T and so:

L L F I T T φ = - from which: L F L I / ) T T ( - = φ
(2)

The angles and angular velocities of B and L ( B B , φ φ and L L , φ φ

) were calculated from the accelerations derived in ( 1), (2) using stepwise integration.

A grid search was made for T max between 10 and 30 Nm, T min between -1 and -10 Nm and t 1 between 0.1 and 0.4 s to find the possible time profiles of T for turns which satisfied all of the following conditions. These were: (a) the leg L rotates 2 radians and stops by the end of the simulation, (b) the body B rotates 2 radians and the angular velocity of B at the end of the simulation is the same as at the start, (c) the foot does not slip more than 0.13 radians in the direction opposite to that of the turn (as for the experimental data of [START_REF] Imura | Biomechanics of the continuity and speed change during one revolution of the Fouetté turn[END_REF]. The coefficient of friction µ = 0.2 and time for one turn was 1.0 s. After determining the bounding time profiles of minimum and maximum T max from the solution space of T, additional cases were considered in which the friction coefficient ranged from 0.1 to 0.3 (t end = 1.0 s) and the time for one turn ranged from 0.7 to 1.0 s ( µ = 0.2), spanning the experimentally determined values of µ = 0.2 and t end = 0.85 s.

Results

For each combination of frictional coefficient µ and turn time t end there were maximum and minimum values for the initial (maximum) value T max of the twisting torque T (Tables 1234). Each pair of max-min solutions gave rotation angle time histories similar to experimental data and had similar angular velocity time histories (Figs. 4 and5). For each value of T max lying between the maximum and minimum values there existed unique values of T min and t 1 for which the body and supporting leg each rotated one revolution and any counter-slipping was less than 0.13 radians as shown in the example (Table 2, Fig. 6). With µ fixed at 0.2, smaller values of t end lead to larger values of minimum T max and consequently to a smaller range of solutions (Table 1). For t end = 0.76 there was a narrow range of solutions with T max lying between 28.1 Nm and 28.3 Nm and for t end = 0.75, or smaller, there were no solutions. For decreasing t end the maximum angle difference d φ between the body and the supporting leg increased (Table 1). T min (Table 3). As µ decreased the maximum and minimum values of T max , T min , and t 1 became closer (Table 3), giving a narrower range of solutions and for µ < 0.12 there were no solutions. The initial angular velocity B φ remained essentially constant for the various values of µ .

For each solution with given values of µ and t end , there existed solutions with equivalent torque profiles for different combinations of µ and t end . For example if µ increased from 0.20 to 0.25, t end decreased from 1.0 to 0.9 in the equivalent solution (Table 4). The torque profiles of these equivalent solutions have the same t 1 /t end values and merely have different scaling factors for torque and time.

*** insert Table 3 here

Discussion

The purpose of this study was to determine the possible time history profiles of T for various tempos and friction. For each combination of time of turn and frictional coefficient there is a range of solutions for T (Tables 1234). The solutions for T satisfy the requirement that the net change in the angular velocity of the body B is zero and that the leg L and body B each rotate 2 at t end . Because the body B rotates under the action of only T, the time integral of T should be zero after one revolution so that there is no net change in angular momentum. The rotation of the leg L is dependent on the net torque T F -T, so the foot should rotate 2 during slipping as a consequence of sufficient acceleration produced by the net torque. Again the time integral of T F -T must be zero since the start and end velocities are zero. Hence, T max , T min , and t 1 are such that the integral of T is zero and T F -T rotates L just one revolution during slipping. For a given value of T max these two constraints give a unique solution for the remaining two degrees of freedom (T min and t 1 ). Thus for a given set of conditions ( µ , t end ) there is a set of torque profiles defined by T max , T min , and t 1 each of which lie within the bounds shown in Tables 1234. While there is a range of solutions for a given set of conditions ( µ , t end ), the tight specification of any individual solution may make the Fouetté turn difficult to perform in a steady state since precise timing would be required by the dancer. On the other hand since there are various time histories within the general profiles which produce one revolution of the Fouetté turn (Fig. 6), such time profiles may represent the different styles from different schools.

*** insert Fig. 6 here *** Initially the twisting torque is positive and less than the limiting frictional torque, accelerating the upper body in the direction of twist while the supporting leg remains stationary with the twisting torque and frictional torque cancelling each other (Figs. 4a and5a). If the twisting torque exceeds the limiting frictional torque in this phase (T > 0) the foot will slip in the direction opposite to that of the twist. Soon after the upper body starts to decelerate relative to the supporting leg and T becomes negative, the magnitude of T exceeds that of limiting T F and the foot slips in the direction of twist (Figs. 4a and5a). Once the magnitude of T falls below that of limiting T F , the angular velocity of the foot decreases until the foot stops rotating at which time T = T F .

For a given coefficient of friction there is a lower bound on the time of the turn since faster turns require larger torques (Table 1) and these are bounded by the limiting frictional torque. This explains why the range of solutions narrows for faster turns (Table 1). There will also be an upper bound on the turn time arising from the assumed time profile of the vertical reaction force N (Fig. 3) since the supporting leg must flex in order to reduce N below one body weight and there is a limit to the amount of flexion possible. The average angular velocity of the body will be inversely proportional to the turn time, and since the angular velocity profiles of the body are similar in different solutions (Figs. 4c and5c) the initial angular velocity B φ will be approximately inversely proportional to t end (Table 1). Thus if a dancer completes a turn with a net decrease in B φ this could be compensated for in the next turn by choosing to continue to turn at the new angular velocity and selecting a twisting torque profile with the corresponding t end . Alternatively, a larger T max for the same T min could be used to produce a net gain in angular momentum in the next turn while keeping the same turn time. Real time adjustments for B φ and L φ could also be made by modifying the time profiles of I B and r.

As the friction coefficient becomes smaller, the magnitude of T also decreases but so does the range of possible profiles for T so that for a given t end there is a lower limit for µ below which there are no solutions. This is a consequence of not having sufficient T F -T to produce the required rotation of the foot. For values of µ below this limit there are solutions with longer turn times since for a given solution there are equivalent solutions with smaller µ and larger t end (Table 4). For larger values of µ , solutions require larger values of T (Table 3), and so there will be an upper limit on µ imposed by the ability of the dancer to exert large torques. For floors with different µ it will be possible to turn at a given tempo providing µ lies within a certain range.

For µ below the lower bound slower turns will be possible and for µ above the upper bound faster turns will be possible up to a limit.

For a given pair of values of µ and t end and a particular solution there are other corresponding pairs of values with an equivalent solution for which t 1 /t end is the same (Table 4). In comparing these solutions a change of t end by a scaling factor k will correspond to a change in each of B I , B φ , and L φ by a factor of 1/k and a change in each of B φ and L φ by a factor 1/k 2 . As a consequence of equations ( 1) and ( 2) T and T F will change by a factor 1/k 2 .

There are a number of simplifications associated with the model. The moment of inertia about the longitudinal axis of the supporting leg is assumed to be constant. In an actual performance the knee bends and extends in order to stand on tiptoe and so the moment of inertia will vary. However knee flexion occurs primarily during full foot contact and changing the moment of inertia in this phase would have minimal effect on a simulation since the foot slips very little during this phase. During the majority of the slipping phase the leg is straight and the assumption of constant moment of inertia is reasonable. Although the time profiles of the variables r, I B , and N are simplifications they were based on the experimental data of a dancer and the calculated leg and body rotation angles were similar to performance data (Figs. 4 and5). The coefficient of static friction has been assumed to be the same as that of dynamic friction rather than a little larger. The effect of this will have been to reduce the range of possible solutions slightly. The limit of 0.13 radians on foot slippage in the direction opposite to the turn constrains the maximum angle difference between body and support leg to anatomically feasible values of around 1 radian or less in the simulations (Tables 1234).

If this constraint is removed the solution space is much larger but includes simulations with large relative rotations between body and support leg which are beyond anatomical limits.

This simple model has been used to describe the solution space of the possible time profiles of the twisting torque T that produce the required rotation about the longitudinal axis in a Fouetté turn. A change of floor to one with increased friction will require a larger twisting torque to turn at the same tempo. A net reduction in the angular velocity of the body B after one turn can be compensated for by increasing T max for the same T min . The model could be applied not only to the Fouetté turn but also to other ballet turns such as the pirouette and Grand Fouetté Italien. φ is the maximum difference between leg and body angles, B φ is the initial angular velocity of the body B Fig. 6. An example of an intermediate twisting torque profile (dashed line) between maximum T max (thick line) and minimum T max (thin line) corresponding to Table 2.

List of Figure Captions

Fig. 1 .

 1 Fig. 1. Sequential view of one revolution of Fouetté turn. Each picture is shown every

Fig. 2 .

 2 Fig. 2. The model comprises the supporting leg L and the remainder of the body B.

Fig. 4 .

 4 Fig. 4. Time profiles of a successful Fouetté turn ( µ = 0.2) for which the twisting

Fig. 5 .

 5 Fig. 5. Time profiles of a successful Fouetté turn ( µ = 0.2) for which the twisting

Fig

  Fig. 1.

Table 1 .

 1 Ranges of the twisting torque parameters T max , T min , t 1 for which the dancer can turn in steady state ( µ = 0.2)

	T	t end	T max -T min t 1 /t end	φ	d	φ	B
	max 1.00 28.3	4.9	0.15 0.66 4.49
	min 1.00 18.4	5.5	0.23 0.76 4.49
	max 0.90 28.3	5.4	0.16 0.68 4.97
	min 0.90 20.2	6.0	0.23 0.80 4.97
	max 0.80 28.3	6.0	0.18 0.75 5.57
	min 0.80 23.9	6.7	0.22 0.91 5.58
	Notes: d φ is the maximum difference between leg and body angles,
	B φ is the initial angular velocity of the body B

Table 3 .

 3 Ranges of the twisting torque parameters T max , T min , and t 1 for which the dancer can turn in steady state (t end = 1.0 s)

	T	µ T max -T min t 1	φ	d	φ	B
	max	0.15	21.2	4.2	0.16	0.69	4.46
	min	0.15	15.9	4.7	0.23	0.86	4.47
	max	0.20	28.3	4.9	0.15	0.66	4.49
	min	0.20	18.4	5.5	0.23	0.76	4.49
	max	0.25	35.3	5.7	0.14	0.65	4.51
	min	0.25	20.1	6.8	0.25	0.89	4.52
	max	0.30	42.4	6.5	0.13	0.65	4.52
	min	0.30	22.2	8.1	0.27	1.05	4.55
	Notes: d φ is the maximum difference between leg and body angles,
	B φ is the initial angular velocity of the body B


Nan Liang

Shao Shuangquan

Changqing Tian

Y Y Yan

Two phase flow instabilities in horizontal straight tube evaporator

It is essential to ensure the stability of a refrigeration system if the oscillation in evaporation process is the primary cause for the whole system instability. This paper is concerned with an experimental investigation of two phase flow instabilities in a horizontal straight tube evaporator of a refrigeration system. The relationship between pressure drop and mass flow with constant heat flux and evaporation pressure is measured and determined. It is found that there is a negative slope section in the middle of positive slope pressure drop versus mass flow velocity thus making the velocity a multi-valued function of pressure drop.

Three types of dynamic instabilities including the density wave instability, pressure drop instability and thermal instability are found under the conditions of heat flux from 5 to 17 kW/m 2 , mass velocity from 150 to 1500 kg/(m 2 •s), and evaporating pressure from 0.5 to 0.8 MPa. The density wave oscillation occurs at almost all mass velocities; its period is about 1~3 seconds and its amplitude is the lowest. The pressure drop oscillation takes place in a negative slope section; its period is about 10 seconds. The thermal oscillation can be induced at a high mass flow velocity, and its period is of 60 seconds and the amplitude is the highest among the

Introduction

The stability which is an essential condition for a successful and high-efficiency refrigeration system has been discussed by researchers recently [START_REF] Zahn | A visual study of two-phase flow while evaporating in horizontal tubes[END_REF][START_REF] Wedekind | Theoretical model for predicting the transient response of the mixture-vapor transition point in horizontal evaporating flow[END_REF][START_REF] Mithraratne | An experimental and numerical study of the dynamic behaviour of a counter-flow evaporator[END_REF][START_REF] Ibrahim | Effect of sudden changes in evaporator external parameters on a refrigeration system with an evaporator controlled by a thermostatic expansion valve[END_REF][START_REF] Liang | Dynamic simulation of variable capacity refrigeration systems under abnormal conditions[END_REF]. The cause of instability may be explained as the influence of the inherent characteristics of two-phase flow in evaporator and also the impact of the system control characteristics [START_REF] Chen | A study on the operational stability of a refrigeration system having a variable speed compressor[END_REF]. The former emphasizes particularly on the oscillations of two-phase flow of the refrigerant in evaporation process; its result is regarded as the basis of the system instability.

The instability induced by the oscillations of two-phase flow boiling, which conventionally refers mainly water system such as steam generators, boiling water reactors and reboiliers, can be classified as either the static instability or the dynamic instability. The static instability, caused by the steady-state performance of a boiling process, mainly includes Ledinegg

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

instability and flow pattern transition instability; the dynamic instability, initiated by the transient changes in boiling process, mainly covers the density-wave oscillations, thermal oscillations and pressure drop oscillations [START_REF] Boure | Review of two-phase flow instability[END_REF][START_REF] Lahey | An assessment of the literature related to LWR instability mode[END_REF].

Indeed, the results of water system have provided necessary experience and methodology for the researches of refrigeration loop. However, some differences between the two-phase flow instabilities in refrigeration system and that in water system can be identified. For example, water systems are often open circulation and the evaporating pressure is governable as a fixed value with a surge tank or some compressible containers, but the similar containers in refrigeration systems do not have such functions; the heat flux of evaporators tube in water systems is often much higher than those in refrigeration systems; and moreover, the working condition in a refrigeration system is also different from that in water system, such as throttling device and the quality at the inlet and exit. These differences may cause the absence of certain types of oscillations. In addition, they are also obstructions for applying the water two-phase flow instabilities theories and the results into refrigeration systems. Some preceding investigations focused on the thermal character of refrigerant with several types of tubes, but the system setup is the same as the water system. Among these results, only Veziroglu et al. [START_REF] Veziroglu | Fundamentals of two-phase oscillations and experiments in single-channel systems[END_REF] investigated the stability boundaries of sustained density-wave oscillations in an electrically heated single channel, up-flow system with R-11 as the test fluid.

Early experimental observations for two-phase flow instabilities of a refrigeration system were reported by Zahn [START_REF] Zahn | A visual study of two-phase flow while evaporating in horizontal tubes[END_REF]. Wedekind and Stoecker [START_REF] Wedekind | Theoretical model for predicting the transient response of the mixture-vapor transition point in horizontal evaporating flow[END_REF] have found the oscillatory motion of mixture-vapour transition point in evaporator and the motion was considered as an oscillatory nature, even for steady flow condition. Based on experimental and theoretical

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
analysis, several researchers including Zahn [START_REF] Zahn | A visual study of two-phase flow while evaporating in horizontal tubes[END_REF], Wedekind & Beck [START_REF] Liang | Dynamic simulation of variable capacity refrigeration systems under abnormal conditions[END_REF], and Barnhart & Peters [START_REF] Barnhart | An experimental investigation of entrained liquid carry-over from a serpentine evaporator[END_REF] suggested that the occurrence of slug flow at the inlet of evaporator should have resulted in the oscillation; and actually the slug flow should be a kind of thermal oscillations or the result of density-wave oscillations. Meanwhile, Wedekind and Stoecker [START_REF] Wedekind | Theoretical model for predicting the transient response of the mixture-vapor transition point in horizontal evaporating flow[END_REF] believed that it was related to the density-wave oscillation. Moreover, Gruhle and Isermann [START_REF] Gruhle | Modeling and control of a refrigerant evaporator[END_REF] proposed a probability that it may be influenced by the strong nonlinear character of heat transfer coefficient. However, the variation of heat transfer coefficient may actually be an integrated result of density-wave oscillation and other oscillations like pressure drop oscillation. In general, since only these a few experimental studies have focused on evaporations in refrigeration system, there is not a general explanation and experimental proof for the instabilities in refrigeration system.

This paper presents the results of an experimental study on two phase flow oscillations of R-22 in a refrigeration system with a horizontal straight tube evaporator.

Experimental measurement

As shown in Fig. 1, the experimental refrigeration system is consist of four main components, namely the compressor, the water-cooled condenser, the evaporator and the electronic expansion valve (EEV), and other accessorial components. As a main part of the experimental system, the evaporator is a horizontal straight copper tube of 3 m long, 8 mm inner diameter and 9.52 mm outer diameter. The tube is wrapped by a thin zonal heater of 3 kW heat input capacity. The heat flux of the evaporator is controlled by a voltage regulator.

All temperatures of the evaporator are measured using copper-constantan thermocouples of displayed and recorded by the data acquisition system and can be processed in the computer.

Detailed information for test devices and their precision is listed in Table 1.

The experiment is performed as follows:

1) At first, open all valves and turn on the compressor, R-22 is transferred through the whole system, and the evaporation pressure is adjusted to a fixed value by changing power frequency of the compressor and the EEV opening.

2) The heat flux is controlled by setting values with the voltage regulator. The evaporation pressure of R22 will change with the variation of heat flux, and it can be adjusted to a fixed value by changing power frequency of the compressor and the EEV opening again. With constant heat flux and evaporation pressure, change the refrigerant mass flow rate, and observe and measure the parameters of two-phase flow in evaporator.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

4) Change the evaporating pressure, and repeat the step 2 and step 3.

Results and discussion

Static instability

Fig. 2 show the results of measured pressure drop between inlet and outlet of the evaporator changing with mass flow velocity in the experiment system in which the pressure of evaporation is of 0.7MPa and the heat flux is of 17 kW/m 2 . It can be seen from Fig. 2 that there is a negative slope region (point 2 to 4) in addition to the two positive slope regions (point 1 to 2, point 4 to 5). In general, the pressure drop rises with the increase of mass velocity. This kind of characteristic is called positive slope. The negative slope in two-phase flow of evaporator is formed in the conditions that the pressure increase by the evaporation exceeds the pressure drop of friction. The negative slope region means that there must be several mass velocities (point 1, 3 and 5) for a given pressure drop, which may cause the oscillation due to the mass velocity leap (between point 1 and 3, or between point 3 and 5) at a fixed pressure drop.

The relationship of pressure drop changes with mass velocity of the evaporator as shown in Fig. 2 is similar to the results of water system observed by Ledinegg [START_REF] Ledinegg | Instability of flow during natural and forced circulation[END_REF], which is called Ledinegg instability and belongs to the static instability. 2.

The density wave oscillation is the results of fluid waves of alternately higher and lower density mixtures due to multiple regenerative feedbacks among the flow rate, vapour generation rate and pressure drop [START_REF] Boure | Review of two-phase flow instability[END_REF]. In the present experiment, the density wave oscillation is observed from very small mass flow rate to large mass flow rate. Figs. second, and it is proportional to the transit time of refrigerant particle through the evaporator.

And the oscillation amplitude of evaporating temperature is of 0.6 o C, and the according mass velocity is about 60 kg/ (m 2 •s).

It has been pointed out that the pressure drop oscillations occur only when the pressure

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
drop decreases with the increase of mass flow rate (or mass velocity), and the oscillation period is governed by the volume and the compressibility of vapour in the water system [START_REF] Lahey | An assessment of the literature related to LWR instability mode[END_REF].

This means that the necessary condition for the occurrence of the pressure drop type oscillation is the pressure drop curve changing with mass velocity with negative slope (point 2 to 4 in Fig. 2) and the existence of a compressible volume (e.g. surge tank) in the flow circuit. In refrigeration system, most compressible volume usually comes from the accumulator and receiver. In addition, the amplitude for evaporating temperature is 1 o C; the amplitude for mass velocity is about 150 kg/ (m 2 •s).

The term of 'thermal oscillation' is often used to describe the situations of temperature fluctuations occurring in solid interacting with fluid. The thermal oscillation used here is related to the instability of boiling heat transfer [START_REF] Lahey | An assessment of the literature related to LWR instability mode[END_REF], and it is characterized by large amplitude fluctuations in the heater wall temperature. The typical character of thermal oscillations is its long period and great amplitude.

In the present experimental study, thermal oscillations have been found at relatively high mass velocity. For example the oscillations are found, respectively, when mass velocity is in the range of 637-750 kg/ (m 2 •s) and heat flux is at 5 kW/m 2 ,; and also when mass velocity is in the range of 959-1500 kg/ (m The onset of oscillation is concerned with the inlet subcooling, system pressure, heat flux at inlet and exit, liquid temperature, cross-section of the channel, entrance throttling and the inlet velocity [START_REF] Saha | An experimental investigation of the thermally induced two-phase flow oscillations two-phase systems[END_REF][START_REF] Aritomi | Instabilities in parallel channel of forced-convection boiling upflow system, (II) Experimental results[END_REF]. The stability thresholds and the correlations for density-wave oscillations The density wave oscillations (DWO) first appear at very small mass velocity (e.g., 150

kg/(m 2 •s)) and it exists the whole mass velocity range (150 to 1500 kg/(m 2 •s)) in our experiments. From experimental results, it can be concluded that the density wave oscillation takes place when the heat flux is high enough for a given mass velocity.

The pressure drop oscillation (PDO) occurs in the region between the broken line a and broken line b, where the slope is negative for the pressure drop versus mass velocity. When the heat flux increases, the boundary between the density wave oscillations and pressure drop oscillations moves to higher mass velocity, this means that the region of density wave oscillations enlarges with the increase of the heat flux.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

The thermal oscillations (TO) occur in the region of right hand side from broken line b and it happens when the heat flux is relatively low at the given mass velocity since the mass velocity is so high that the liquid R-22 cannot be evaporated adequately at a low heat flux.

In refrigeration system, the users always focus on the temperature at evaporator outer wall while the wall temperature is more stable than the actual evaporating temperature. As a result, density wave oscillations can be seen as a kind of oscillations with little harm to system operation because of its low amplitude and short period. However, the pressure drop oscillations and thermal oscillations should be paid more attention for their high amplitude and long period, and the temperature at evaporator outer wall will then appear to be unstable.

Mathematical modelling of instability boundary

The boundaries of the given state were obtained from the experimental results. As the onset conditions of a given oscillations always change along with the working conditions, it is necessary to develop a mathematical model to determine the boundaries at different conditions.

In the present study, a homogeneous phase model is proposed to simplify the governing equations. To develop the mathematical model, the following assumptions have been made :

1) the flow inside the tube is one dimensional homogeneous phase flow;

2) there is a thermodynamic equilibrium between two phases;

3) the fluid is incompressible;

4) the dissipated heat by friction, kinetic energy and potential energy is negligible in the energy balance equation.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

5) the heat flux is homogeneous.

6) the axial heat conduction of tube wall and the heat stored in tube wall is negligible.

Based on the physical phenomena described, the governing equations of the flow in terms of continuity, momentum and energy equations can be written in Eqs. ( 1) to (3):

Continuity: d 0 d G z = (1)
Momentum:

2 2 2 2 2 f a 1 d d d 2 d (1 ) d d d d ( 1
)
g l p p p fG x x G z z z D z φ ρ ρ α ρ α     -     = + = - - +           -         (2)
Energy: d( ) d

Gh DQ z π = (3)
Considering the above homogeneous phase model, Eqs. ( 1) to ( 3) can be transferred into integral equations along the tube length. Then the equations can be converted to dimensionless format as

2 2 l G C gDρ = (4) 0 1 2 2 2 2 2 2 2 2 2 (1 ) 2 d 1 (1 ) z l z l l l g p G G x x f z gD gD D gD ρ φ ρ ρ ρ α ρ α     ∆ - = + + -         -     ∫ (5) 0 0 i fg r h Q x Gz h G π ∆ + = (6)
It is supposed that,

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT ( , , ) g l a f x ρ ρ = (7) 2 ( , ) g l f x ρ φ ρ = (8)
So that the equations can be expressed as:

2 2 l G C gDρ = , (9) 2 2 fg
, , , , ,

g i l l l h p G L f x f gD gD D h ρ ρ ρ ρ   ∆ ∆ =       , ( 10
) 0 0 ( , , ) i fg fg fg Q Dz h h L f x h Gh h D π ∆ ∆ = = . (11)
Then, with these equations, the non-dimensional equations can be further expressed as: ,,,,, )

g i fg fg l l l h Q L G p f f Gh h D gD gD ρ ρ ρ ρ ∆ ∆ = . ( 12
)
For a given system, the correspondingly change of f and L/D can be omitted. So that the final mathematical expression for the oscillation formation condition is: ,, )

g i fg fg l l l h Q G p f Gh h gD gD ρ ρ ρ ρ ∆ ∆ = . ( 13
)
With the experiment results of evaporating pressure 0.7MPa, empirical correlation for the

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

onset of pressure drop oscillations and that of thermal oscillation are obtained as Eqs.( 14) and ( 15) respectively.

g i fg l fg l l h Q G p Gh gD h gD ρ ρ ρ ρ         ∆ ∆ =                   (14)
2.243 0.852 1.174 0.329 2 2 1.649

g i fg l fg l l h Q G p Gh gD h gD ρ ρ ρ ρ         ∆ ∆ =                   (15)
A comparison between the calculated boundaries and measured boundaries is shown in 11 that the trend of calculated boundaries is similar to those measured ones. However, there is still an error between the simulation and the measurement, and the corresponding error is about 30% ± while the root mean square error is about 20% ± . The corresponding comparison between the calculated boundaries and the measured ones is given in Table 3.

Based on Eqs. ( 5), ( 6), ( 14) and [START_REF] France | Measurement and analysis of dynamic instabilities in fluid-heated two-phase flow[END_REF], the boundary between density-wave oscillations and pressure drop oscillations and the boundary between pressure drop oscillations and thermal oscillations can be calculated for a given refrigeration system when the mass velocity, heat flux, and evaporating pressure are given.

In practical design, a system should be tested before operation. Generally speaking, there are several parameters which can be measured for evaporation, such as the degree of refrigerant subcooling at the expansion valve inlet, the pressure of condensation, the pressure

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
of evaporation, the heat flux (heat load and heat transfer area), the degree of superheat at evaporator outlet, etc. It is also convenient to record these parameters when the system oscillations occurred. If the refrigerant at expansion valve inlet is subcooled, the period and amplitude of the oscillations will be measured and recorded. Then the oscillations will be classified by the mass velocity, heat flux and evaporating pressure. For the working range (of evaporating pressure), it is possible to draw out such onset boundaries. If possible, it is better to design an evaporator tube with no negative slope in the curve of pressure drop versus mass velocity to avoid pressure drop oscillations. If the evaporator is designed with a negative slope in the curve of pressure drop versus mass velocity, the working range should be set to the region on the left hand side of the pressure drop oscillation onset boundary. If the working range crosses the unstable region (the region on the right hand side of the pressure drop oscillation onset boundary), the system control strategies should be prepared on how to adjust the system parameters to make the system stable when the boundaries are determined.

Conclusions

In this paper, the static and dynamic instabilities of two phase flow in a horizontal straight tube evaporator of refrigeration system are investigated.

The static instability of evaporator can be expressed by pressure drop versus flow rate. The characteristic curve has a negative slope region in addition to two positive slope regions, thus making the mass velocity a multi-valued function of the pressure drop. The multi-valued characteristic provides the potentiality for static oscillations and dynamic oscillations.

Three types of dynamic oscillations have been found and observed, namely, the density

List of table captions

  oscillations. The boundary as the onset of pressure drop oscillation and that of thermal oscillation are obtained experimentally. Finally, the empirical equations for the boundaries between three types of oscillations are obtained according to the mathematical model and test data.Keywords: evaporator; instability; refrigeration; two phase flowNomenclature

Fig. 3

 3 Fig.3shows the influence of heat flux on pressure drop when the pressure of evaporation

4 and 5

 5 show the measured transient behaviour of evaporation temperature and mass velocity, respectively, with the heat flux of 17 kW/m 2 when the density wave oscillation happens. The pressure of evaporation is referred to the measured value from the pressure sensors in the middle of evaporator, and the temperature of evaporation is the corresponding saturated temperature to the pressure of evaporation. The oscillation period of density wave oscillations is about 4

Figs. 6

 6 Figs. 6 and 7 show the measured transient behaviour of evaporating temperature and

  M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPTmeasured transient behaviour of evaporating temperature and mass velocity respectively during the period of thermal oscillation which occurs when heat flux is of 17 kW/m 2 and mass velocity is of 900-1000 kg/ (m 2 •s). The period of thermal oscillations is about 60 second, and the amplitude of evaporating temperature is of 3 o C. There are many short period density wave oscillations in a cycle of the long period thermal oscillation, and the oscillations shown in Figs.8 and 9can be seen as an integrated results of thermal oscillations and density wave oscillations.3.3 Instability boundaryIt is necessary to determine the boundaries of the regions where the density wave oscillations, pressure drop oscillations and thermal oscillations occur. From a mass of the test data, the boundary as the onset of pressure drop oscillation and that of thermal oscillation can be obtained, and the empirical equations for the boundaries among density wave oscillations, pressure drop oscillations and thermal oscillations are developed according to the mathematical model and test data.

[ 15 -

 15 20], pressure drop oscillations [21] and thermal oscillations[22] in water system were given and the effects of various parameters on stability boundaries and oscillation periods were studied. The equations of instabilities have been reviewed by Kakac and Bon[22].However, it can be identified that they reviewed equations are aimed only at a certain type ofM A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPToscillations; and the aim of the present researches is to simulate the period and amplitude of oscillations. Such theoretical simulations and analyses have the advantages of understanding the mechanism of a certain type of oscillations, but the results are lack of applicability related to practical designing, testing and controlling of refrigeration system. In this paper, a mathematical model is built on the basis of measured instability boundaries to describe the relations between the oscillation boundaries and these parameters which can be measured.The boundary between density-wave oscillations and pressure drop oscillations is the onset of pressure drop oscillations, and the pressure drop oscillations can occur at the right hand side of the boundary. Similarly, the boundary between pressure drop oscillations and thermal oscillations is the onset of thermal oscillation, and the thermal oscillations can take place at the right hand side of the boundary. The boundary as the onset of pressure drop oscillation (broken line a) and that of thermal oscillation (broken line b) were obtained according to large numbers of the test data (see Fig.10).

Fig. 11 .

 11 Fig.11. The experimental data is obtained from the test under the conditions of heat flux from

  pressure drop oscillations and thermal oscillations. The density wave oscillation, which takes place at almost all the mass velocities, has the smallest periods and the lowest amplitude among three types. Pressure drop oscillations are also incidental and the period is about 10 seconds. The thermal oscillation occurs with longest periods and the highest amplitude among three types.With the experiment results, the experimental onset conditions of pressure drop oscillations and thermal oscillations are obtained. The empirical correlations of the onset of pressure drop oscillations and that of thermal oscillation are also obtained with the homogeneous phase model. The empirical correlations can be used to analyze the onset conditions for pressure drop oscillations and thermal oscillations in a certain R-22 system. The method can be used to guide practical design through some simple tests.The results in this paper provide a further understanding for the mechanism of two-phase flow instabilities of refrigeration system. These results can also be used to provide guidance for design and control strategy of refrigeration system in order to avoid serious oscillation problems. a convective boiling upflow system, International Journal of Heat and Fluid Flow 15 (3) (1994) 241-246.[19] Z. Gao, J. Li, L. Qian, Analysis of stability in a low pressure, low quality and two phase flow with natural circulation, in: Proceedings of the 2nd International Symposium on Multiphase Flow and Heat Transfer, Xian, China, Hemisphere, Washington, DC, 1990. [20] H. Lai, Y.Y. Yan, The effect of choosing dependent variables and cell-face velocities on convergence of the SIMPLE algorithm using non-orthogonal grids. International Journal of Numerical Methods for Heat & Fluid Flow, 11(5/6) (2001) 524-546. [21] L. Guo, Z.P. Fen, X.J. Chen, Transient convective heat transfer of steam-water two-phase flow in a helical tube under pressure drop type oscillations, International Journal of Heat and Mass Transfer 43 (2002) 533-542. [22] S. Kakac, B. Bon, A Review of two-phase flow dynamic instabilities in tube boiling systems, International Journal of Heat and Mass Transfer 51 (2008) 399-433.

Fig. 2 Fig. 3 Fig . 4 Fig. 5 Fig. 6 Fig. 7 Fig. 8 Fig. 9 17 kW/m 2 )Fig. 10 Fig. 1 Fig. 2 Fig. 3 Fig. 4 Fig. 5 Fig. 6 Fig. 7 Fig. 8 Fig. 9 Fig. 10

 234567891721012345678910 Fig. 1 Experimental system

Table 1 :

 1 Detailed information for the test devices

Table 2 :

 2 Experiment data of different heat flux

Table 3 :

 3 Comparison of calculated boundaries point and experimental boundaries point

	M A N U S C R I P T

Table 1 :

 1 Detailed information for the test devices

	Components	Model	Remark
	Compressor	YZB-21R	Variable frequency from 30 to 120Hz
	EEV		DPF 1.6	Variable pulse number from 0 to 550
	Evaporator		Self-designed	Horizontal copper tube with 3m long, inner
				diameter of 8 mm and outer diameter of 9.52
				mm
	Condenser		SS-075GT
	Refrigerant	Mass 2100	R22 flow rate: 0-1500 kg/h
	mass meter	flow		Precision: 0.1kg/h ±
	Water	flow	LZB-10	Flow rate: 0.36-3.6 m 3 /h
	meter			Precision: 0.5% ±

Table 2 :

 2 Experiment data of different heat flux

	Heat flux	Oscillation	Period	Average Amplitude ( o C or kg/ (m 2 •s))
	(kW/m 2 )	Type	(s)	Temperature	Mass velocity
	5	Density wave	0.8	0.9	110
		Pressure drop	9.5	2.1	228
		Thermal	55.4	4.7	1053
	9	Density wave	1.0	0.8	100
		Pressure drop	9.9	1.9	195
		Thermal	57.6	4.2	1028
	11	Density wave	1.1	0.7	93
		Pressure drop	10.2	1.7	185
		Thermal	58.7	4.0	1008
	13	Density wave	1.3	0.7	88
		Pressure drop	10.5	1.6	170
		Thermal	59.3	3.5	968
	17	Density wave	1.6	0.6	75
		Pressure drop	10.7	1.4	150
		Thermal	60.4	3.2	950

Table 3 :

 3 Comparison of calculated boundaries point and experimental boundaries point

	Heat flux (kW/m 2 )	Evaporating Pressure (MPa)	Calculated boundaries point (kg/ (m 2 •s))	Experimental boundaries point (kg/ (m 2 •s))	Error (%)
			PDO	TO	PDO	TO	PDO	TO
		0.5			397.4	545.2	7.6	7.5
		0.7	425.4	589.7	470.0	637.5	10.5	8.1
		0.8			524.7	720.4	23.3	22.1
		0.5			484.7	621.5	11.5	13.8
	9.0	0.7	547.5	712.8	517.3	686.4	5.6	3.8
		0.8			621.8	812.1	13.6	13.9
		0.5			511.6	655.6	13.4	13.2
	11.1	0.7	590.3	755.2	591.4	751.1	0.2	0.5
		0.8			674.5	885.7	14.3	17.3
		0.5			567.5	697.9	11.5	14.7
	13.4	0.7	641.6	817.9	624.3	806.1	2.7	1.5
		0.8			726.2	935.9	13.2	14.4
		0.5			657.5	864.6	7.3	3.9
	16.9	0.7	709.1	899.6	734.2	958.8	3.5	6.5
		0.8			800.2 1012.4	12.8	12.5

Acknowledgement

The study was supported by Natural Science Foundation of China (Grant No. 50676099).


N Pardo

Á Montero

A Sala

J Martos

J F Urchueguía

Efficiency improvement of a ground coupled heat pump system from energy management

Keywords: Ground coupled heat pump, Energy efficiency, Numerical simulation

The installed capacity of an air conditioning system is usually higher than the average cooling or heating demand along the year. So, most of the time, the system is working under its actual capacity. In this contribution, we study the way to improve the efficiency of a ground coupled heat pump air conditioning system by adapting its produced thermal energy to the actual thermal demand. For this purpose, an air conditioning system composed by a ground coupled heat pump and a central fan-coil linked to an office located in a cooling dominated area was simulated, and a new management strategy aiming to diminish electrical consumption was developed under the basic constraint that comfort requirements are kept. This strategy takes advantage of the possibility of managing the air flow in the fan, the water mass flows in the internal and external hydraulic systems, and the set point temperature in the heat pump to achieve this objective. The electrical consumption of the system is calculated for the new management strategy and compared with the results obtained for a conventional one, resulting in estimated energy savings around 30%.

-INTRODUCTION

It is a well known fact that building equipments are responsible of around the 40% of all greenhouse gas emissions produced in developed countries, being approximately 60% produced by cooling and heating systems [START_REF] Orme | Estimates of the energy impact of the ventilation and associated financial experience[END_REF][2][START_REF] Champsour | Recensement de la Population 1999 -Des logements Plus Grands et Plus Comfortables[END_REF]. In this context, the development of air conditioning (HVAC) systems should be oriented to improve its energy efficiency and to reduce its environmental impact.

Ground Coupled Heat Pumps (GCHP) systems represent a technically viable technology for heating, cooling and domestic hot water systems in buildings [START_REF] Lund | World-wide direct uses of geothermal[END_REF][START_REF] Lund | World-wide direct uses of geothermal energy[END_REF][START_REF] Lund | Geothermal heat pumps -an overview[END_REF][START_REF] Sanner | Current status of ground source heat pumps and underground thermal storage in Europe[END_REF][START_REF] Urchueguía | Experimental comparative analysis of a ground coupled heat pump system versus a conventional air-to-air heat pump in typical conditions of the European Mediterranean coast[END_REF][START_REF] Omer | Ground-source heat pump systems and applications[END_REF]. They offer several interesting characteristics for the potential user, such as a lower electrical demand and maintenance requirements than conventional air to water heat pump systems and, therefore, lower annual cost [START_REF] Catan | An optimized ground coupled heat pump systems design for northern climate applications[END_REF][START_REF] Martin | A design and economic sensitivity study of single pipe horizontal ground coupled heat pump systems[END_REF][START_REF] Lund | Ground Source (geothermal) heat pumps[END_REF][START_REF] Urchueguía | Comparison between the energy performance of a group coupled water to water heat pump system and an air to water heat pump system for heating and cooling in typical conditions of the European Mediterranean Coast[END_REF]. From a non economical point of view, they offer competitive levels of comfort compared with standard technologies, reduced noise levels and visual impact, savings of greenhouse gas emissions and reasonable environmental safety. Furthermore, the Environmental Protection Agency (EPA) recognizes ground source systems among the most efficient and comfortable heating and cooling systems available today [14]. A proof of the goodness of this technology is the amount of installed units worldwide, estimated in 1.1 million [START_REF] Spitler | Ground-Source Heat Pump System Research -Past, Present and Future[END_REF].

In the standard design of an air conditioning system, the references taken to estimate the heating and cooling capacity of the system to be installed are usually based on the coldest and the warmest day along the year. A ground coupled heat pump air conditioning system can be designed to meet the whole load or just a part of it, using an auxiliary system to achieve comfort conditions in days in which thermal demand is higher than the capacity of the heat pump system. In both cases, the thermal energy required by the thermal load is under the actual capacity of the system during most part of the time. In this context, the development of strategies for the operation of the air conditioning system and, particularly, for the system based on the ground coupled heat pumps, allowing to adapt the thermal energy generated by the system with the thermal load is a good way to improve the system energy efficiency while satisfying the thermal comfort [START_REF] Fraisse | An Analysis of the Performance of Different Intermittent Heating Controllers and an Evaluation of Comfort and Energy Consumption[END_REF][START_REF] Liao | An experimental study on an inferential control scheme for optimising the control of boiler in multi-zone heating system[END_REF][START_REF] Fong | HVAC system optimization for energy management by evolutionary programming[END_REF][START_REF] Fraisse | A numerical comparison of different methods for optimizing heatingrestart time in intermittently occupied buildings[END_REF]. An approach to achieve this objective based in the combination of the ground coupled system with other production system and in the decoupling of energy production from energy distribution was presented in [START_REF] Pardo | Optimization of hybrid -ground coupled and air source -heat pump systems in combination with thermal storage[END_REF]. This contribution tries to achieve the same objective just taking advantage of the possibility of managing the operational variables of the ground coupled heat pump system. Following this idea, the aim of this work is to evaluate the energy savings that a new management strategy can produce in a HVAC system composed by a ground coupled heat pump and a central fan coil linked to a standard office space. In this new management strategy, the air mass flow in the fan, the water mass flow in the internal and external hydraulic system and the set point temperature in the heat pump, usually fixed in conventional strategies, have the possibility of a continuous regulation that allows us to design a more efficient way to achieve the desired thermal comfort.

This new management strategy is based on five capacity levels developed from the total electrical power equation of the HVAC system. In our particular case, this equation indicates that to achieve energy savings is desirable to work with low water flows. For this reason, the new management strategy tries to achieve a steady state in which the water mass flow is maintained at the lower level that guarantees the production of enough thermal energy to satisfy the thermal demand.

An office space in a cooling dominated area is modelled to evaluate in these conditions the energy performance of the ground coupled heat pump HVAC system when it is managed by the new management strategy and by a conventional one. The thermal comfort criterion which has to be satisfied is the Predicted Mean Vote index (PMV). This comfort index predicts the mean value of votes of a large group of people in the Thermal Sensation Scale and it is defined by the ISO7730-1994 standard [START_REF] Iso | Moderate thermal environments-determination of the PMV and PPD indices and specification of the conditions for thermal comfort[END_REF].

We compare the annual electrical energy consumption of the HVAC system when it is managed by our new management strategy and by a conventional one and an evaluation of its energy efficiency is presented. We also analyze the factors which allow improving the energy efficiency of the air conditioning system by the new management strategy. m & rated,air : Maximum mass flow rate of air that can pass through the fan m & con,fc : Flow rate of condensate draining from the coil h air,in,fc : Enthalpy of air entering the coil h air,out,fc : Enthalpy of air exiting the coil h con,fc : Enthalpy of condensate draining from the coil

Nomenclature

-SIMULATED SYSTEM

In this section, we present a description of the simulated office area as well as the PMV index which is the criterion used to evaluate the thermal comfort. Afterwards, we describe the ground coupled heat pump system linked to the office area, the energy model of the employed devices, the new management strategy and the conventional one. Finally, these systems are simulated using TRNSYS software package [START_REF]TRNSYS 16.1: Transient System Simulation Program[END_REF].

-Simulated Office Area

The simulated office area comprises 108 m² (12m x 9m) with two windows in the façades north and south and three in the façades east and west. To perform the simulation the office area is characterized by its building materials, its dimensions, distribution and orientation. There are four different kinds of construction elements: external walls, floor, roof and window glasses. External walls are defined as ventilated façades composed by four elements: perforated brick, 5 cm of insulation, air chamber and a Naturex plate cover; its global conductivity is 0.51 W/m 2 K. The floor and the roof are built with hollow blocks with 5cm of insulation; its global conductivity is 0.51 W/m 2 K. Finally, the window is composed by a glass, with solar radiation transmissivity equal to 0.837 and conductivity equal to 5.74 W/m 2 K, and a window frame with conductivity equal to 0.588 W/m 2 K. The windows size is 1.5 m 2 , dedicating a 15% of this area to the frame surface. The internal and external shadow factor for these windows is estimated in 0.7.

These parameters are included in the TRNSYS software package through the TRNBuild tool specifically designed to simulate the thermal behaviour in a multi-zone building area. Finally, the weather data base of the Spanish city of Valencia is used to characterize the Mediterranean coast weather.

-Comfort criteria

The Predicted Mean Vote (PMV) index is one of the criteria to estimate the comfort state proposed by the ISO7730-1994 standard. We choose this parameter to evaluate the comfort state in our simulation. The PMV index for the air conditioned space is calculated at each simulation step time. To do it, three parameters, activity, thermal resistance of clothing and mean air velocity has to be estimated. We follow ISO7730-1994 recommendations for this purpose. The value for the activity parameter is 1.2 met for moderated office activity. The thermal resistance of clothing is 1.0 clo, when the considered working clothes are shirts, trousers, jackets, and shoes. Finally, the value for the mean air velocity is estimated in 0.1 m/s. The other three parameters, air temperature, mean radiant temperature and partial water vapour pressure are computed at each time step.

-Air Conditioning System and Electric Power Consumptions

The air conditioning system linked to the simulated office area is based on a ground coupled heat pump. This system is composed by two parts: the external circuit composed by the external water pump (EWP) and the geothermal heat exchanger (GHE), and the internal circuit, composed by the internal water pump (IWP) and a central fan coil (FC) linked to the office area. The water to water heat pump (WWHP) is a reversible heat pump and it is the connecting element between both circuits. In figure 1 a schematic diagram showing the layout of these devices is included.

FIGURE 1

The only difference between the hydraulic circuit for the new management strategy and the conventional management strategy is in the active devices (WWHP, EWP, IWP, EMF) employed. For the new one the active devices have a continuous regulation, whereas for the conventional one these devices can be only switched on or switched off.

The following paragraphs detail the model describing the behaviour of each component. We also derive the way to compute the electrical power consumption of each one, for given working conditions. Notice that the purpose of this study is to improve the energy efficiency of the system keeping the comfort requirements.

Therefore, the knowledge of its electrical power consumption behaviour (see below eq. ( 9)) allows us to develop a suitable management strategy to achieve this objective. [START_REF] Hellström | Duct Ground Heat Storage Model[END_REF] is used to simulate our GHE, composed by three U-tube vertical boreholes of fifty meters depth placed in parallel with a distance between them of 5.3 meters. This model is implemented in TRNSYS library as type557. To estimate the size of the ground heat exchanger the method proposed in reference [START_REF]Ground-source heat pump project analysis[END_REF] was employed. U-tube pipes are characterized by an outer radius equal to 0.01664 m., an inner radius equal to 0.01372 m., a center to center half distance equal to 0.0254 m., and a pipe thermal conductivity equal to 0.42 W/m K. Boreholes of radius equal to 0.1016 m. are filled with bentonite (with thermal conductivity equal to 1.3 W/m K). Finally, ground thermal properties are characterized by a storage thermal conductivity equal to 2 W/m K, and a storage heat capacity equal to 2016 kJ/m3/K.

Ground Heat Exchanger (GHE). 'Duct Ground Heat Storage Model'

To calculate the electrical energy consumption of the water to water heat pump the outlet temperature of the ground heat exchanger (T GHE,out ) is needed. This temperature is computed by the 'Duct Ground Heat Storage Model' as:

a in GHE out GHE T T T ⋅ - + ⋅ = ) 1 ( , , β β (1)
Where β is a dumping factor defined as follows:

        ⋅ - =         ⋅ - = f f p p f f v Q C L Q C V α α β exp exp (2)
Water to Water Heat Pump (WWHP). The model employed for the water to water heat pump (type 666 from TRNSYS library) calculates its power consumption from the operational variables. The heat pump power consumption (P HP ) is calculated as the quotient between its capacity at current conditions (G) and its coefficient of performance at current conditions (COP), all multiplied by the fraction of full load power, f flp :

flp HP f COP G P = (3)
The fraction of full load power is a function of the heat pump Partial Load Ratio (PLR), the ratio between the thermal load and heat pump capacity at current conditions. These three quantities needed to calculate P HP are related on different ways with the operational variables. The heat pump capacity, G, and the coefficient of performance, COP, are calculated using the following expressions:

] , [ ,out GHE set ratio rate T T G G G = (4) ] , [ ,.out GHE set ratio rate T T COP COP COP = (5)
G rate and COP rate are the heat pump capacity and the coefficient of performance both in nominal conditions. In table 1, the values for G rate and COP rate used in our simulations are included. G ratio and COP ratio are dimensionless coefficients calculated from the values of the operational variables T set , set point heat pump temperature, and T GHE,out , temperature of the water returning from the GHE. These quantities are read from two data files (standard files from TRNSYS library) and characterize the behaviour of the heat pump for the different values of the operational conditions. Another data file (also available together with type 666 from TRNSYS library) is employed to relate the fraction of full load power, f flp, and the partial load ratio, PLR.

Fan Coil (FC).

The fan coil model (type 774 from TRNSYS library) considers that the final air temperature is the average temperature of the fluid in the coil. Note that the inlet air temperature to the coil is the ambient temperature plus the heating effect due to the heat losses of the fan electric motor. The heat absorbed by the fluid is the same extracted from the air minus the losses due to condensation as can be read from equation [START_REF] Lund | Geothermal heat pumps -an overview[END_REF].

fc cond fc con fc out air fc in air fc air fc fluid h m h h m Q , , , , , , , , ) ( & & - - = (6)
The electrical consumption of this device comes from the fan electric motor, and it is computed as follows:

where P rated,fan and m & rated,air are the rated fan power consumption and the rated air mass flow when the fan is operating at full-speed. In table 1, the values for P rated,fan and m & rated,air used in our simulations are included.

Finally, m & air is the air mass flow through the fan in each time step.

Fan power consumption is calculated in common software tools as a function of the percent design airflow using a second order equation [START_REF] Zheng | Development and validation of a fan performance model for typical packaged HVAC systems[END_REF]. The choice of equation ( 7) is done to be representative of a generic case.

Other choices in which the second order equation slightly differs from equation ( 7) produces a very low impact in energy performance results, confirming that this choice can be considered as representative of a generic case.

Water Pump (IWP, EWP). We use variable speed water pumps for the internal and external hydraulic circuit. The pump speed will be adjusted by a pressure regulator so that the pressure drop keeps constant in the hydraulic circuit. As a consequence the electric pumping consumption is linear with respect to the water mass flow through the pump. This linear relation is defined in the following expression:

water rated water pump rated pump m m P P , , & & = ( 8
)
where P rated,pump and m & rated,water are the rated pump power consumption and the rated water mass flow when the pump is operating at full capacity. In table 1, the values for P rated,pump and m & rated,water used in our simulations are included. Finally, m & water is the water mass flow through the pump in each time step.

TABLE 1

Electrical consumption of the air conditioning system. The addition of the consumption from the different devices is the total electric power consumption of our air conditioning system. Equation ( 9) allows us to calculate this total power consumption at every time step. In this expression, the constant coefficients C FC , C IWP , C EWP and C WWHP can be obtained by looking in the previous electrical consumption expressions for each device, equations ( 3), ( 7) and (8

m T T T f T T COP T T G C m C m C m C P & & & & + + + = (9)
This equation shows the relation between the electrical consumption of the air conditioning system and six A proper understanding of the behaviour of this power consumption expression is the key to design a management strategy that will improve the energy efficiency of the system.

The consumptions of the electric motor of the fan coil and the water pumps only depend on the air mass flow and the water mass flow respectively, and they are independent of the electric consumptions of the other devices. In contrast, the consumption of the heat pump depends on its set point temperature and the internal water mass flow and this last variable is included in the consumption equation of the internal water pump.

Let study this expression for given conditions of T GHE,out and T load,in . The non trivial behaviour of equation ( 9), corresponding to the heat pump plus water pumping power consumption is shown in figure 2. This figure

shows the overall electric consumption of these two devices as a function of the internal water mass flow and the set point for the outgoing temperature in the cold side when delivering 5 kW in cooling mode with a heat pump inlet temperature from the ground heat exchanger of 24ºC.

FIGURE 2.

From figure 2 we see that to achieve energy savings it is desirable to work with low water flows and high set point temperatures. Notice that the change in electric power consumption is higher in the direction of the internal water mass flow than in the set point temperature of the heat pump. Therefore, if it is necessary to supply more energy it is more convenient to modify first the set point temperature keeping low water mass flows. Once the set point temperature achieves its minimum allowed value the water pump should increase its water mass flow.

In addition to the conclusions derived from figure 2, it would be desirable to have similar flows in both sides of the heat pump. Therefore, the ratio between internal and external water mass flow must be the same with time. Finally, notice that the electric motor of the fan is the unique device which can be switched on with independence of the heat pump. And the internal and external water pumps along with the heat pump have to be switched on all together.

-Management Strategies

In this section we present both management strategies: the new one based on the behaviour of equation ( 9)

and the conventional one.

For both systems, we have fixed two aspects: the working period and the heating and cooling seasons. The working period in the office is fixed from 9:00 to 18:00, which is the period when is switched on the air conditioning system. The heating season is considered from November to December and from January to March and the cooling season from April to October. Finally, there are not holiday periods, as a consequence the air conditioning system works every day.

-New management strategy

The objective of the new management strategy is to improve the energy efficiency of the air conditioning system by adapting properly its generated energy to the actual thermal demand in the office. The new management strategy is detailed in the following paragraphs.

To maintain neutral comfort conditions in the office the PMV value should be zero in the Thermal Comfort Scale. The deviation of the PMV from zero indicates a variation of the thermal demand in the office area.

Notice that, to satisfy this thermal demand, there are several configurations of the management variables which are suitable to compensate this PMV deviation. Our management strategy is based on the choice of a particular configuration of these variables. This choice classifies the air conditioning system capacity in five capacity levels given by the steady state value of the management variables. In figure 3, we show a diagram illustrating this classification. We now describe each capacity level.

FIGURE 3.

• First capacity level; in steady state conditions the fan capacity is between the 0% and the 50% of its maximum capacity and the other devices are switched off. In this level the blown air by the fan modify the convective factor and homogenize the temperature in the office to achieve the thermal comfort.

• Second capacity level; the fan, the hydraulic pumps and the heat pump are switched on. In steady state conditions the air blown by the fan is fixed to the 50% of its maximum capacity, the water mass flows of the internal and external hydraulic pumps are between the 0% and the 10% of its maximum allowed value, and the set point temperature is 0%, which indicates the lowest or highest set point temperature for heating or cooling mode respectively.

• Third capacity level; in steady state conditions the air blown by the fan is fixed to the 50% of its maximum capacity, the water mass flows of the internal and external hydraulic pumps are fixed to the 10% of its maximum allowed value and the set point temperature is between 0% and 100%, meaning that the set point temperature can be any value of its range in heating or cooling mode.

• Fourth capacity level; in steady state conditions the air blown by the fan is fixed to the 50% of its maximum capacity, the water mass flows of the internal and external hydraulic pumps are between the 10% and the 100% of its maximum allowed value and the set point temperature is 100% which indicates the highest or lowest set point temperature for heating or cooling mode respectively.

• Fifth capacity level; in steady state conditions the fan capacity is between the 50% and the 100% of its maximum capacity, the water mass flows of the internal and external hydraulic pumps are fixed to the 100% of its maximum allowed value and the set point temperature is 100%.

If the energy supplied by the air conditioning system when all the active devices are given the 100% of its capacities is not enough to maintain the thermal comfort conditions, the PMV diverts from zero and the thermal demand is not satisfied.

The choice of these five capacity levels is based on the behaviour of equation ( 9). This equation, explained in previous section, indicates that to achieve energy savings is desirable to work with low water flows in the air conditioning system. For this reason, the new management strategy tries to achieve a steady state in which the water mass flow is maintained at the lower level that guarantees the production of enough thermal energy to satisfy the thermal demand.

We also want to point out that in actual conditions a high air flow blown by the fan can produce an excessive noise due to the vibration in the fan coil, to avoid this effect, the new management strategy tries to maintain in steady state conditions the air mass flow blown by the fan around half of its maximum capacity.

allows to quickly achieve the thermal comfort state and, after this, drive the management variables to the steady state of the capacity levels described before without comfort penalties. In the simulations, the action of the controllers has place every 28.8 seconds, which is fast enough in comparison with the building dynamic behaviour.

-Conventional management strategy

In the following paragraphs, we are going to describe the conventional management strategy for the air conditioning system.

This strategy uses on/off regulators to manage the air conditioning system. This kind of regulators has only two possible operation points. They give its maximum capacity when switched on and nothing when switched off. The on/off regulators are installed in active elements: the heat pump, the electric motor of the fan coil and the internal and external water pumps.

In this strategy, three aspects are important. First, the difference between the set point temperature and the inlet temperature of fluid in the heat pump from the office area indicates the connection or disconnection of the generator energy system, composed by the heat pump, the external water pump and the ground heat exchanger. Therefore, the internal water pump is always switched on to have a measure of this difference during the working period. Second, to maintain the comfort conditions, the PMV is located in a comfort band between 0.5 and -0.5, in order to try to avoid an excessive number of connections and disconnections of the system which, in an actual situation, could damage the actuators of the fan coil. Finally, the set point temperature in the heat pump is constant, for the heating season is fixed to 45ºC and for the cooling season to 7ºC.

In heating mode the conventional management strategy works as follows. When the value of the PMV variable is below the lower limit of the comfort band, PMV equal to -0.5, the electric motor of the fan is switched on. Then, the air goes through the coil, where a constant water mass flow is pumped from the internal circulation pump. This heat exchange produces a variation of the temperature of the water in the internal circuit which is detected by the heat pump; immediately this device and the external water pump are switched on to provide the necessary energy to fix the value of the PMV to the upper limit of the comfort band, PMV equal to 0.5. After achieving this value, the heat pump, the external water pump and the electric motor of the fan are switched off.

In cooling mode the conventional management strategy is the same as the previous one except that the references of the comfort band are inverted. The electric motor of the fan is connected when the PMV is above the upper limit, PMV equal to 0.5, and the air conditioning system stops supplying energy when the PMV arrives to -0.5.

-SIMULATION RESULTS AND DISCUSSION

In the following paragraphs, we present and discuss the results obtained from the simulations. First, we evaluate the accuracy of the obtained results. And second, we present the energy consumption of the air conditioning system for both management strategies, evaluating the energy savings achieved by the new one.

To evaluate the accuracy of the simulation results a set of simulations are carried out for several time steps.

The obtained results have to show a tendency towards a value independent of the time step employed to solve the set of differential equations programmed in our software package. Simulations are performed at the following solver time steps: 14.40 seconds, 7.20 seconds, 3.6 seconds, 1.8 seconds and 1.44 seconds.

The results obtained for the conventional management strategy are almost independent of the solver time step chosen, with variations among them always smaller than 1%. Nevertheless, the results from the simulation of the new one show a higher dispersion. This behaviour is due to the fast dynamic introduced by the controllers which are not possible to simulate correctly when the solver time step is not small enough. Indeed, it is a wellknown fact that simulating multiple levels of control with widely different time constants leads to a system of stiff ordinary differential equations, for which the numerical integration routines in TRNSYS may not be specially adapted: this is the reason of needing a time step around 1 second in order to reliably simulate a whole year. In spite of this, the tendency of the electrical energy consumptions in all the simulations is the same through the year. In the cases where is used the smallest time steps, 1.80 seconds and 1.44 seconds, the obtained results do not show significant variation between them (always smaller than 3%) and they point out to the same continuum values. Now, we present the energy consumption of the air conditioning system for both management strategies, and we evaluate the energy savings achieved by the new one. In the simulation conditions, the weather database employed models the Mediterranean coast weather, which is characterized to have hot summers and warm winters and the working period coincides when the external ambient temperature and the solar radiation is the highest throughout the day. The combination of these two factors produces that the energy demand in cooling mode is much higher than in heating mode.

In figure 4 we present the monthly electrical energy consumptions of the air conditioning system for the two management strategies. These data correspond to the simulations employing a solver time step equal to 1.44 seconds. From this figure we can see that in all months, except April, the system consumes less electrical energy when it is managed by the new management strategy.

FIGURE 4

In winter season, our simulation results show that the influence of the external temperature and the solar radiation in the office area during the working period is enough to maintain the thermal comfort in it during most days of the winter season, therefore the energy provided by the air conditioning system is very low. A remarkable difference between both management strategies exists because the conventional strategy employs the temperature change of the water in the internal hydraulic circuit to communicate the energy demand in the office room to the generator system. As a consequence, the internal water pump is always switched on during the working time and, therefore, consuming electrical energy independently of the air conditioning necessities. Whereas, in the new management strategy, the air conditioning system only consumes electrical energy according with the thermal demand in the office area. Then, the unnecessary electrical consumptions, which can suppose a great energy loss throughout the year, are avoided.

In summer season, the influence of the external temperature and the solar radiation which were an advantage during the winter season are now a disadvantage. Therefore, both management strategies are more active to provide the necessary high cooling power to achieve the thermal comfort. In these conditions is when the new management strategy has more chances to manage the air conditioning system to reduce its electrical consumption. In figure 4, the electrical consumption for both strategies grows up in function of the cooling demand, being the maximum electrical consumptions during the warmest months, July and August. In all months during the cooling season, except April, the electrical consumptions of the air conditioning system when is employed the new one are significantly smaller than the ones achieved when using the conventional one. Furthermore, the new strategy improves the efficiency of the air conditioning system as the cooling demand increases.

In April, the electrical consumption of the air conditioning system is smaller when the conventional management strategy is used. This behaviour is because in the first day of this month the air conditioning system changes from heating mode to cooling mode. When the conventional management strategy manages the air conditioning system the value of the PMV index is between 0.0 and 0.5 most of the time belonging to the working periods. In these conditions, this strategy keeps switched off the fan, the heat pump and the external water pump. Meanwhile, the new strategy activates these three devices to supply energy to the office to compensate the small deviation of the PMV from the neutral thermal comfort state. This shows us the difficulty to choose the suitable dates to change from heating mode to cooling mode and vice versa which are particulars for each system. The previous situation can be avoided by delaying the date to activate the cooling mode of the heat pump.

Finally, the annual electrical consumption of the air conditioning system when it is managed by the conventional management strategy is estimated in 3359 kWh, and in 2289 kWh when managed by the new one (data for solver time step equal to 1.44 s.). Then, the annual energy savings achieved are estimated in a value around the 30%.

-CONCLUSIONS

In this paper, a new management strategy for a HVAC system composed by a ground coupled heat pump and a central fan coil linked to an office area is presented. The air mass flow in the fan, the water mass flow in the internal and external hydraulic systems and the set point temperature in the heat pump are managed to obtain the desired comfort in the office while reducing the electrical consumption of the air conditioning system.

This new management strategy is based on five capacity levels developed from the total electrical power consumption equation of the HVAC system. This equation indicates that to achieve energy savings is desirable to work with low water flows in the air conditioning system. For this reason, the new management strategy tries to achieve a steady state in which the water mass flow is maintained at the lower level that guarantees the production of enough thermal energy to satisfy the thermal demand.

The total power consumption of the HVAC system when it is managed by the new strategy and by a conventional one are compared, obtaining that the annual electrical consumption is reduced in a value around the 30% when the air conditioning system is managed by the new management strategy.

To conclude, the energy savings obtained in this study are particular of the chosen conditions of our simulations. Nevertheless, we think that this improvement in the energy efficiency of the HVAC systems driven by ground coupled heat pumps is an estimation of the effect when these kinds of systems are managed by the new management strategy. The presented new management strategy can be easily adapted to other areas with similar characteristics and conditions and the energy savings achieved should be of the order of the savings obtained in this work.

P

  pump : Circulation pump electrical power consumption P rated,pump : Rated circulation pump electrical power consumption P fan : Fan electrical power consumption P rated,fan : Rated fan electrical power consumption Q load : Load met by the heat pump P HP : Heat pump electrical power consumption at current conditions G : Heat pump capacity at current conditions G ratio : Heat pump capacity at current conditions divided by the rated capacity COP : Heat pump coefficient of performance at current conditions COP rate : Heat pump rated coefficient of performance at current conditions COP ratio : Heat pump COP at current conditions divided by the rated COP f flp : Heat pump fraction of full load power PLR : Heat pump partial load ratio T set : Set point temperature -Desired outlet temperature of fluid in the heat pump fluid stream T load,in : Inlet temperature of fluid in the heat pump load stream T GHE,in : Inlet temperature of fluid to the ground heat exchanger T GHE,out : Outlet temperature of fluid from the ground heat exchanger T a : Temperature of the surrounding ground β: dumping factor α v : heat transfer coefficient between the fluid and the surrounding ground per unit of volume α p : heat transfer coefficient between the fluid and the surrounding ground per unit length of pipe V: Volume of storage L p : total pipe length in the storage volume C f : Volumetric heat capacity of the fluid Q f : Total fluid flow rate in the ground heat exchanger m & water : Mass flow rate of fluid passing through the circulation pump m & rated,water : Maximum mass flow rate of fluid that can pass through the circulation pump m & water,IWP : Mass flow rate of fluid passing through the internal circulation pump m & water,EWP : Mass flow rate of fluid passing through the external circulation pump Cp water : Specific heat of fluid m & air : Mass flow rate of air passing through the fan

  system variables: air mass flow ( m & air ), internal water mass flow ( m & IWP ), external water mass flow ( m & EWP ), heat pump set point temperature (T set ), heat pump inlet temperature from the ground heat exchanger (T GHE,out ) and heat pump inlet temperature from the load (T load,in ). Out of these six variables the first four are management variables ( m & air , m & IWP , m & EWP , T set ) .

  This work has been supported by the Spanish Government under projects "Modelado y simulación de sistemas energéticos complejos" (2005 Ramón y Cajal program), "Modelado, simulación y validación experimental de la transferencia de calor en el entorno de la edificación" (ENE2008-0059/CON). A. Sala is grateful to the financial support of grants DPI2008-06731-c02-01 (Spanish Government), and Generalitat Valenciana prometeo/2008/088.

Figure 1 .

 1 Figure 1. Air conditioning system: the geothermal heat exchanger (GHE), the external water pump (EWP), the internal water pump (IWP), the water to water heat pump (WWHP) and the fan coil (FC) with the electric motor of the fan (EMF).

Figure 2 .

 2 Figure 2. Electrical power consumption for the heat pump and the internal water pump when the heat pump supplies 5 kW in cooling mode with a heat pump inlet temperature from the ground heat exchanger of 24ºC.

Figure 3 .

 3 Figure 3. Diagram illustrating the classification in capacity levels of the system capacity given by the steady state values of the management variables. For a given thermal demand there is a unique choice for the operational point of the air conditioning system.

  ).

	total	FC	2 air	IWP	IWP	EWP	EWP	WWHP	[ ratio ratio	, set set [	] out out , GHE , GHE ,	]	flp	[	load	, in	,	GHE	,	out	,	set	,	IWP	]