Sena Jeong

Eojin Kim

Kyungho Jeong

Doyoung Jeon

Yong Bok Lee

Effects of Residual Imbalance on the Rotordynamic Performance of Variable-Speed Turbo Blower

Keywords:

In the variable-speed turbo machinery, it is important to ensure the reliability in operation speed of the large area unlike the constant speed rotating machine. In particular, as a study on the performance maintenance and minimal vibration, the high-speed balancing and field balancing are necessary to remove vibrations due to the unbalanced mass. In using a dedicated balancing device, to add or remove the weight is often not obtained satisfactory results for expanded the range of variable rotating speed. Thus, in order to improve the performance and driving stability of variable speed machine, it has to make the residual imbalance to a minimum, through the field balancing.

This study presents the experimental results to apply the two-plane multi-speed balancing techniques, and ensure the driving stability by reducing the residual imbalance. This paper also shows rotordynamic simulation results of the rigid rotor and bearing system, including damped critical speed and imbalance responses. The operating maximum speed of the turbo blower is up to 20,000 rpm, usually it works normally in the 4,000 ~ 17,000 rpm as the variable-speed region. And this study was tried to reduce the residual imbalance vibration in all operating speed regions. First, by utilizing an influence coefficient method the rotor was adjusted to the ISO 1940 standard grade G1. Additional, as the least square method, the experiment was carried out to minimize the residual imbalance by applying the multi-speed balancing in the field. The first rigid mode was predicted in the near 17,500 rpm, the experimental results showed that rigid mode is consistent. Accordingly, the field balancing was performed at several points of the rigid mode speed region. This study was an effort to reduce excessive vibration due to the rigid mode by minimizing the residual imbalance. As a result, the vibration caused by the residual unbalance can be reduced to 70%, and experimentally verified that this can effectively improve the rotordynamic performance.

INTRODUCTION

According to high-speed rotating machine market increase, recently there is required a variable-speed turbo machine for energy saving and efficient operation [START_REF] Howard | Rotordynamic feasibility of a conceptual variable-speed power turbine propulsion system for large civil tilt-rotor applications[END_REF]. Basically the rotor bearing systems should be designed to avoid the critical mode from the operation region. In addition, performing the high-speed balancing is an essential requirement for satisfying durability and allowable vibration range of the turbo machine [START_REF] Rao | Rotor dynamics[END_REF]. Even though the design of the rotor completely, in actual industrial applications, imbalance vibration of the system can be increased due to the misalignment of the rotor assembly. Especially, the variable-speed turbo machine must be ensure the stability of the system with high precision balancing of the various speed. This balancing technology is largely divided into a modal technique and influence coefficient method. Modal technique has the effect of reducing the amplitude of the residual vibration of imbalance occurring in each mode by using the orthogonality of the modes [START_REF] Rieger | Balancing of Rigid and Flexible rotors, SVN 12, The shock and vibration center[END_REF][START_REF] Bishop | The Vibration of Rotating Shafts[END_REF][START_REF] Bishop | Second-Order Vibration of Flexible Shafts[END_REF]. However, there are some shortcomings to know the critical speed of the rotating shaft and the corresponding mode in advance. On the other hand, influence coefficient method does not require prior information of the rotor mode. It is a method of measuring magnitude and phase as installed a trial mass in sequence to each balancing plane. It can obtain the correction mass and angle from calculation of the residual imbalance at the desired speed [START_REF] Thearle | Dynamic Balancing of Rotating Machinery in the Field[END_REF]. The influence coefficient method is also sensitive to other system parameters which can result in an ill-conditioned influence coefficient matrix thereby producing erroneous corrections weights [START_REF] Darlow | The Identification and Elimination of Nonindepentent Balance Planes in Influence Coefficient Balancing[END_REF]. Further, when utilizing least square method proposed by Goodman [START_REF] Goodman | A Least-Squares Method for Computing Balance Corrections[END_REF], this approach allowed vibration at multiple locations and multiple speed conditions. Palazzolo [START_REF] Palazzolo | Multimass Flexible Rotor Balancing by the Least Squares Error Method[END_REF] also presents to be effective in reducing the overall vibration level by minimizing the sum of squares of the residual vibration. Moreover, multiple balance planes can be accommodated by expanding this equation which sets to zero the vibration magnitude at as many locations or speed conditions as balance planes are available. It is known as the exact point balancing procedure by Tessarzik and Badgley [START_REF] Tessarzik | Flexible Rotor Balancing by the Exact Point-Speed Influence Coefficient Method[END_REF][START_REF] Tessarzik | Flexible Rotor Balancing by the Exact PointSpeed Influence Coefficient Method[END_REF][START_REF] Tessarzik | Experimental Evaluation of the Exact Point-Speed and Least-Squares Procedures for Flexible Rotor Balancing by the Influence Coefficient Method[END_REF][START_REF] Badgley | Recent Developments in MultiPlane-MultiSpeed Balancing of Flexible Rotors in the United States[END_REF]. This method has been widely applied to balancing of large power generation turbines with many rotor critical speeds and bearing points to minimize vibration [START_REF] Lund | Analysis and Experiments on Multi-Plane Balancing of a Flexible Rotor[END_REF][START_REF] Kirk | Insights from Applied Field Balancing of Turbomachinery[END_REF]. It provides good results as long as the vibration effects at all points and speeds. As a case of practical study, Gunter [START_REF] Gunter | Field Balancing of 70 MW Gas Turbine Generators[END_REF] was provided the practical aspects of multi-plane field balancing 70 MW gas turbine generators. The optimum balance was determined by means of a constrained least squared error multi-plane program using vibration data at the various speeds and power levels.

This study describes the experimental results relating to remove residual imbalance over the entire speed range, including the operation speed. In order to satisfy the condition that minimizes the residual vibration in all driving regions of the variable-speed turbo machine, two-plane multi-speed balancing technology was applied in the field. Since the environment and condition of the turbo system applied to the different industrial field, and large complex turbo system is very difficult to assemble repeatedly. Then, it should be perform the high-speed precision balancing through trial and error in the field directly. The least square method was easily utilized in order to minimize the residual unbalance from the measured vibration. In particular, the speed selection is very important for vibration reduction desired and rotordynamic performance enhancing. Least square method can minimize the sum or average of the squared residual vibration data. However there is a limitation. When the outlier data was a mixed, it can make the wrong result. Thus, multi-speeds selection can be achieved experimentally for field balancing optimization. These can choose in some speed from the rigid mode vibration amplitude peak. By selecting within 10% from the peak of vibration amplitude, when the highest and low of two speed point, this experiment was effective in applying the least square method. Therefore, it could effectively reduced not only the vibration of the corresponding mode speed, also the vibration of the entire operating speed range. On the other hand, applying a lot of data, the balancing effect is decreased. Lastly, the rotordynamics simulation was performed with respect to the residual unbalance to obtain the test. By comparing the dynamic behavior of vibration analysis and balancing a result, it shows that can be minimized of the vibration response with residual imbalance of variable-speed turbo blower.

DESIGN AND MODE ANALYSIS OF THE RIGID ROTOR FOR VARIABLE-SPEED

The impeller of overhung structure is driven by a permanent magnet synchronous motor, both ends supported by the ball bearings as shown in Figure . 1. The total length of the rotor is 646.5 mm, weight is 21.25 kg. A finite element model of the rotor was performed base on the natural frequency analysis, the first bending mode is 1,068 Hz (64,071 rpm). In this study, since the maximum operation of the variable speed up to 20,000 rpm, the design margin is sufficient for the bending mode. Ball bearings in the 7200 series, the inner diameter is 40mm, allowable speed of oil-lubricated is 29,000 rpm. The stiffness of the bearing is 6x10 7 N/m. The quasi-static stiffness is calculated in Ref. [START_REF] Harris | Rolling Bearing Analysis[END_REF]. The rotordynamics analysis was performed by using these values, Figure 2 shows the campbell diagram for natural frequency analysis map. The results of nutural frequency analysis, backward is 17,566 rpm and forward is 19,467 rpm of the first rigid mode (cylindrical mode). The second rigid mode (conical mode) backward is 22,548 rpm and forward is 25,979 rpm. On the basis of this analysis result, during variable speed operation, the rotation stability may be concerned that excessive vibration of the bearings in the rigid mode of 17,000 ~ 20,000 rpm. In order to solve this problem, there are several ways. It is a method of raising the bearing mode to more than 20,000 rpm by increasing the bearing stiffness or axial preload. Alternatively, a method of re-design of the rotor bearing systems. But the redesign can be quite a waste of time and money, including the aerodynamic design or motor part. The last option, to apply a high-precision balancing method for reducing vibration of a desired operating speed range. Therefore, despite the limited condition, it is the best choice in a manner to prevent excessive vibration near the rigid modes.

INFLUENCE COEFFICIENT METHOD

The two plane balance corrections were calculated by peak to peak vibrations and phase measurements. This procedure is still being used with success and it is in essence the fundamental concept in which the influence coefficient method was derived for multi-speed rotor systems. An influence coefficient method with an optimal correction balancing algorithm is developed for balancing a high-speed rotor system.

The influence coefficient procedure [START_REF] Gunter | Balancing of multimass flexible rotors: Part Ⅰ theory and part Ⅱ experimental results[END_REF] for final balancing of a rotor at operating speed is as follows:

1. We choose two correction planes A and B, two measurement planes a and b. To obtain phase data, we need to install keyphasor. 2. Take the zero rotor vibration data, the rotor frequency vibration amplitudes and phase angles of all the planes for the original condition of the rotor. 3. Obtain the trial mass data of A plane and B plane.

To get the trial mass data, a trial mass is inserted in one balance plane, and the resulting rotor-frequency vibration amplitudes and phase angles of all the planes are measured. This procedure is repeated for the other balance plane.

1 [] bB bA BB aB aA AA ww B ww A                                (1)
Where, In addition such as least square balancing algorithm can be satisfied allowable residual vibration levels in overall speed ranges. The number of balancing plane is n, the number of desired balancing speed is l1, and the number of position measured vibration is l2. The number of measurement m is below.

12 m l l  ( 3
)
When m is larger than n, the correction mass vector {W} relates to the unbalance vibration vector {A} as below.

[

] { } { } {0} { } WA      (4)
Where [α] is the influence matrix of m x n, {0} is a zero vector of n th row and {ε} is the residual vector. As residual vector {ε} is lower, we can obtain the better performance of balancing. Therefore, the least square of residual vibration is as follow.

* * * 1 1 2 2 22 1 { } { } . ...... . {[ ( )] [ ( )] } . T mm m ii i Min Min Re Im Min                       (5)
Where {ε} T is a transpose of {ε} and composed of a conjugation complex numbers. By inserting (5) into (4), we can obtain the following equation.

[ ] [ ] { } [ ] { } {0} [ ] [ ] { } [ ] { } TT TT WA WA             (6)
As the equation ( 6) is n square matrix, it is possible to calculate inverse matrix. Thus, the equation ( 6) is converted to correction mass {W}.

1 { } {[ ] [ ]} [ ] { } TT WA      (7)

ZERO ROTOR DATA AND BALANCING SPEED CONDITION DETERMINATION

In order to apply influence coefficient the method, two of correction plane was determined as shown in figure 3. In addition, as part of a multi speed balancing technique, to reduce the residual unbalance of the rigid mode the speed point should determine from zero rotor data.

In the A and B plane, a distance of center and trial mass position is 120mm and 15mm respectively as shown in Figure 4. By using balancing dedicated machine, each residual imbalances of A and B plane are 172.44 g-mm at 167.66˚ and 92.25 g-mm at 178.07˚. The imbalance response simulation is performed with above data. This result shows in Figure 5. It indicates that a vibration increases caused by residual imbalance in high speed. The general speed of the variable-speed turbo blower is from 4,000 to 17,000 rpm, and a maximum operating speed is 20,000 rpm. The balancing speeds can choose in some speed from the rigid mode vibration amplitude peak. It was selected within 10% from the peak of vibration amplitude, the speed of highest peak is 19,500 rpm and the speed of low vibration amplitude is 17,000 rpm. In the critical speed range shown in Figure 5, the balancing speed point of four case that the following was optional determined.

EXPERIMENT RESULTS

Table 1 shows the data of inserted trial masses accordance with run type. Trial mass of A plane is 1.369 g at -11.25˚ and trial mass of B plane is 2.039 g at -22.5˚. In fact, A plane, trial mass were inserted into 0.697 g at 0˚ and -22.5˚, respectively. This can be seen to try 1.369 g at -11.25˚ by vector calculation.

Table 2 shows the measured imbalance vibration data of each rotating speed. In accordance with the trial mass inserted, both of the vibration and phase angle were measured. Hence, the residual imbalances can be calculated based on the previously determined multispeed case by utilizing least square method.

Table 3 shows the residual imbalance results of multi-speed balancing cases. The results was obtained the residual imbalance vibration in each multi-speed case by applying the influence coefficient and least square method. This study focused on the vibration of A plane, because the rotor has heavy load on one side, as overhung type. According to the analysis of the imbalance calculation for each multispeed case, the lowest residual imbalance result can be obtained in the case 2 that 17,000 and 19,500 rpm are selected. On the other hand, when comparing the cases 2 and 3, 18,500 rpm data is being added to the least square calculation. The residual unbalance of case 3 was not significantly reduced as much as case 2. This mean that the data of 18,500 rpm does not help in finding the optimal imbalance reduction. For multi-speed balancing, it is very important to determine the speed because an influence coefficient method is dependent in the speed. The least square method can minimize the sum or average of the squared residual vibration data. However there was an inappropriate data mixed, it cannot make the optimal result. Accordingly, this study was selected for optimum balancing speed as in case 2. In another cases, it did not appear less residual imbalance than case 2, Further, many multispeed points and influence coefficient are not getting good balancing effect. Table 4 shows the correction mass and angle calculated by the vibration and the phase measured in each case. This result can be carried out vibration response test by applied correction mass to the corresponding hole. On the other hand, correction's phase angle was not exactly be applied because it can only put mass into correction plane with intervals of 22.5˚. By a vector calculation, the experiment was tried to be as close as possible to apply correction results. The results of case 2 and case 3 has the smallest residual imbalance by the residual imbalance calculations. Thus, the operating test of case of 2 and 3 were performed by applying the correction mass, respectively. Figure 6 provides that simulation and experimental result of imbalance vibration responses are compared. There are three cases represented a set of experiments (zero rotor data, correction data of case 2, correction data of case 3). The vibration measurements were performed on the 'a' and 'b' plane as shown in Figure 3, it was expressed as front and rear, respectively.

Table 1. Data of inserted trial masses

Overall, experimental results and simulation trend have shown that good agreement. In particular, it can be confirmed that large vibration occurred in the rigid mode speed region. This imbalance responses were also found that it is considerably matching simulation and experiment results. Structurally of the rotor, the vibration in the rear just shows runout level due to assembling and production tolerances. Thus, the present research has focused on the front vibration level of the impeller position. Experimental and analytically, excessive vibrations of the near 17,000 rpm were successfully verified that approximately 70% decrease by applying the correction result of the case 2. Further, it showed the effect of the influence coefficient experimentally in the determinant speed, to obtain results that can significantly reduce the variation in the high-speed. It could be reduced not only the vibration of the corresponding mode speed, also the vibration of the entire operating speed range. By comparing the rotordynamic behavior of vibration analysis and balancing result, it shows that can be minimized of the vibration response with residual imbalance of variable-speed turbo blower.

Finally, Figure 7 is presented the orbital vibrations during the speed-down. From 20,000 to 16,000 rpm, vibration obit is shown the effect of residual imbalance on the operating performance in each balancing case. The case 2 can be proven that the good performance of the vibration, also the speed determination is the most important decisions in order to obtain an excellent balancing effect of influence coefficient.

CONCLUSION

This paper has provided a successful experiment and simulation results for verifying residual imbalance effect and rotordynamic performance with two-plane multi-speed balancing by applying influence coefficient method. A rotor bearing system application of a variable turbo blower, experiments were conducted to ensure rotational stability. Base on the least square balancing technique, the experiment was determined the multi-speed point. It was measured peak to peak vibration and phase angle at each speed, according to the multi-speed condition, calculated the residual unbalance and correction value. The rotational speed of the variable speed application turbo machine is up to 17,000 rpm, the maximum rotational speed is 20,000 rpm. Thus, the present study was tried to reduce the vibration near 17,000 rpm in field. This experiment showed the effects that can reduce the vibration more than 70%. In conclusion, for the most effective use of influence coefficient, it is very important to determine the balancing speed. Based on the balancing effect of the residual imbalance was demonstrated that to obtain a rotational stability of the variable-speed turbomachinery.

Figure 1 .

 1 Figure 1. Rigid rotor and mode shape analysis results.

Figure 2 .

 2 Figure 2. Campbell diagram for natural frequency analysis map.

4 .

 4 Calculate the response coefficients. The equations of influence coefficient (α) is bellow. In this study, we selected the balancing planes and measurement points as shown in Figure 3. A1, A2, A3 are measured vibration data of zero rotor, first trial data, second trial data in A plane. And B1, B2, B3 are measured vibration data of zero rotor, first trial data, second trial data in B plane. Also, TA and TB are first trial mass and second trial mass. The [α] and [w] are influence coefficient and corrected balance mass, respectively.

1

 1

5 .

 5 By using the influence coefficients to obtain a set of correction masses in the two balance planes that will reduce the residual vibration of the rotor.

Figure 3 .

 3 Figure 3. Indicated correction planes and vibration measurement points for balancing.

Figure 4 .

 4 Figure 4. Schematic view of correction A and B plane.

Figure 5 .

 5 Figure 5. Imbalance response simulation result at 'a' plane.

Figure 6 .

 6 Figure 6. Comparison of simulation and experimental results with balancing (during the speed-up).

Figure 7 .

 7 Figure 7. Orbital vibrations in impeller side during the rotor speed-down: (a) zero rotor data, (b) correction of case 2 and (c) correction of case 3.

•

  Case 1 -17,000 rpm, 18,500 rpm • Case 2 -17,000 rpm, 19,500 rpm • Case 3 -17,000 rpm, 18,500 rpm, 19,500 rpm • Case 4 -17,000 rpm, 19,000 rpm, 20,000 rpm

Table 2 .

 2 Measured imbalance vibration data

	Zero rotor	'a' plane	'b' plane
	Rotating speed [rpm]	Magnitude [p-p, mm]	Phase angle [degree]	Magnitude [p-p, mm]	Phase angle [degree]
	17,000	0.1750	-179.2	0.0091	134.9
	18,000	0.1838	-170.9	0.0313	137.3
	18,500	0.2472	179.4	0.0207	118.9
	19,000	0.4425	-179.6	0.0243	95.7
	19,500	0.6830	-160.7	0.0360	102.6
	20,000	0.6905	-150.6	0.0383	99.2
	Trial A plane	'a' plane	'b' plane
	Rotating speed [rpm]	Magnitude [p-p, mm]	Phase angle [degree]	Magnitude [p-p, mm]	Phase angle [degree]
	17,000	0.1527	-24.7	0.0245	-172.4
	18,000	0.1395	-0.1	0.0147	-130.7
	18,500	0.1285	-19.8	0.0121	179.8
	19,000	0.1882	-18.6	0.0192	178.2
	19,500	0.2847	18.2	0.0119	179.9
	20,000	0.2705	55.4	0.0234	138.7
	Trial B plane	'a' plane	'b' plane
	Rotating speed [rpm]	Magnitude [p-p, mm]	Phase angle [deg.]	Magnitude [p-p, mm]	Phase angle [deg.]
	17,000	0.3343	177.6	0.0308	42.0
	18,000	0.3668	-173.1	0.0384	81.6
	18,500	0.4536	-180.0	0.0385	73.2
	19,000	0.6125	-175.8	0.0469	68.7
	19,500	0.7713	-154.2	0.0493	79.0
	20,000	0.7650	-147.2	0.0601	71.3

Table 3 .

 3 Residual imbalance data of multi-speed balancing

		Speed [rpm]	A plane Magnitude Phase [p-p, mm] [deg.]	B plane Magnitude Phase [p-p, mm] [deg.]
	case1	17,000 18,500	0.2137 0.1881	-23.7 158.4	0.1184 0.0933	179.3 129.5
	case2	17,000 19,500	0.0124 0.0045	79.4 -38.8	0.0364 0.0867	-154.4 87.3
		17,000	0.2375	-25.9	0.092	-163.5
	case3	18,500	0.2016	157.4	0.0514	121.2
		19,500	0.0022	104.3	0.0843	114.8
		17,000	0.3141	21.9	0.1256	-113.1
	case4	19,000	0.3668	-137.5	0.0149	157.1
		20,000	0.1423	95.3	0.0932	89.0

Table 4 .

 4 Correction mass and phase angle data in each multi-speed balancing cases

		A plane		B plane
		Correction	Phase	Correction	Phase
		mass	angle	mass	angle
		[gram]	[deg.]	[gram]	[deg.]
	case1	0.937	-2.29	0.271	-17.80
	case2	1.045	-9.50	1.006	-49.92
	case3	1.015	-10.56	0.601	-61.15
	case4	1.080	-19.92	1.343	-90.75

ACKNOWLEDGMENTS

This material is based upon the works supported by the Knowledge Economy Technology Innovation Programs: "Driving efficiency maximizing technology development of pump/blower module for energy saving of buildings". The authors thank to researchers for their help in preparing the variable-speed turbo-blower test setup.


Zhenyue Ma

Qianqian Wu

Study on Dynamic Characteristics of Rotor-bearing Model in the Hydraulic-mechanical-electric Coupled System of Hydro Generator Unit

Keywords:

A hydraulic-mechanical-electric coupled model for hydroelectric shaft-bearing system was established, which contains the penstock model of water-carriage system, the model of governor system and the model of exciter system. A differential equations for coupled model, consisting of the key parameters from different systems (such as head, flow, speed and exciting current), were established and solved simultaneously. Further,a rotor-bearing system model is coupled to the whole system by User Programmable Features (UPFs).On this basis, a newly method which was focused on the study of rotor system dynamic nonlinear properties during the process of operation condition's changing was introduced. In this paper, the model was applied in the rotor dynamic analysis in aplicable start-up laws with different parameters. This analysis method could provide a reference to the stable and optimum operation of hydro generator units.

INTRODUCTION

The operation of hydro generator is a complex process with hydraulic-mechanical-electric coupled. A copuled hydropower system contains the hydro-turbine system,speed control system and generator system, which is involved in hydraulic transient model, mechanical transient and the electric-magnetic transient study, respectively.

The hydro-turbine system may be simplified as linear model and nonlinear model. For the linear model,more assumptions are built (such as non-elastic water,the proportional relationship between flow and gate positon and so on) , on the other hand, the nonlinear relationship of the turbine output, the turbine flow, the turbine head and the gate opening are considered in the study of nonlinear model.A linear model was used to study the factors' value for operating stablity with small load variation by Hovey [START_REF] Hovey | Optimum adjustment of hydro governors on Manitoba hydro system[END_REF]. With taking into consideration the compressibility of water and elastic of penstock, Wozniak and Fett [START_REF] Wozniak | Conduit representation in closed loop simulation of hydroelectric system[END_REF] establised a nonlinear model, and special emphasis was given to the pressure conduit in its properties which were mentioned above. On the basis of nonlinear model, Sanathanan [START_REF] Sanathanan | Accurate low order model for hydraulic turbine-penstock[END_REF] discussed the accuration of the low order model for hydro-turbine with long penstock. Similarly, Kishor et al. [START_REF] Nand Kishor | Simulation of reduce order hydro turbine models to study hydraulic transient charactersitic[END_REF] and Vournas [START_REF] Vournas | Second order hydraulic turbine models for multimachine stability studies[END_REF] too has proposed low order models for above study.In papers [START_REF]Dynamic models for steam and hydro turbines in power system studies[END_REF] [START_REF]Hydraulic turbine and turbine control models for system dynamic studies[END_REF], several hydropower plant models of linear and nonlinear were reported by IEEE group to research the effects of nonelastic and elastic water column on the plant.In recent years, many researchers give significant focus on the model....Over the past decades, many studies were carried out for the control theory of hydro-turbine governor system. A fundamental and pioneering study has been done by Paynter [START_REF] Paynter | A palimpsest on the electronic analog art[END_REF], hovey and Leum [START_REF] Leum | The development and field excitation of a transistor electric governor for hydro turbines[END_REF], while, their papers establised general guideling for hydro-turbine control theory. The permanent droop and generator damping which are neglected by above studied, are considered into the governor models in the papers of Chaudhry [START_REF] Chaudhry | Governing stability of a hydroelectric power plant[END_REF] and others researchers [START_REF] Agnew | The governing of Francis turbines[END_REF]. The PID control strategy is developed and studied for hydro-turbine control by many literatures unitl recent years. And with the delvelopment of computer technology, there is growing concern about the digital control system [START_REF] Findlay | Microprocessor based adaptive water turbine governor[END_REF][13] [START_REF] Murthy | Analysis and improvement of the stability of a hydro turbine generating unit with long penstock[END_REF] which is more sensitivity than conventional control system.In modern reseach, more control theory were deleloped such as optimal control, nonlinear control, adaptive control, robust control and so on [START_REF] Luquing | Variable structure and the time-varying parameter control for hydroelectric generating unit[END_REF][16] [START_REF] Herron | A state space pressure and speed sensing governor for hydro generators[END_REF][18] [START_REF] Luqing | An intelligent self-improving control strategy and its microprocessorbased implementation to a hydro-turbine governing system[END_REF].

The model of rotor-bearing system is the core component of hydropower generator unit. Mostly, the model is consisted of the shaft, the rotor, the turbine and the support structure. The dynamic characteristics of this model are the objective for reseachers to study,while the vibration of the rotor system is mainly affected by the source which copuled hydraulic,mechanical and electric factors. The self-excitation oscillation of water seal for turbine has discussed by many papers. The dynamic response of rotor system under the axial thrus from water is studied in paper [START_REF] Xueli An | Dynamic response of a rub-impact rotor system under axial thrust[END_REF]. The magnetic effect on the system commonly mannifested as a unblanced magnetic pull on the rotor centre. The eccentric force of rotor and turbine is the mechanical effect].

All of studies in the past of the transient (or coupled transient) and the rotor system are indepenednt. For the former the objective is mainly foucus on the stability of the control/power sytsem or condition's change on the basis of parameters optimization.The latter gives a significant attention to the structure dynamic response of the unit shaft system, and the loads effect on the rotor system is summarized form the vabration source system which are time-invariant. In these disscusion, the loads are changed with many factors except time. For the reason mentioned above, it is neccessory to carry out the research for the dynamic characteristic of rotor-bearing system coupled with the hydraulic system, mechanical system and electric system.

In this paper, a nonlinear copuled model of hydropower plant is development to study the structure's/rotor system dynamic characteristics during the coupeld transients process which are resulted in startup laws. First, a elastic penstock model is establised on the basis of the continuity and momentum equations describing the general behaviours of fluids in a pressure duct in terms of two variables, namely, H, piezometric head, and Q, fluid flow. A PID control strategy is used to model the governor system for hydro-turbine, which represented as a third order differential equation in mathematic. A third order model of synchronous generator which is depicted by three differential equations is employed. Then, a simultaneously differential equations with above sub-system equations were established for coupled system. Finally, a FEM model of classical rotorbearing system is developed.

Q i -discharge in penstock at the node i H i -pizeometric head in penstock at the node i H n -net water head of turbine D,A -the cross section diameter and area of penstock D 1 -the diameter of hydro-turbine H np -the operating head of hydro-turbine Q p ,Q 1 ′-hydroturbine's discharge and unit discharge a-water hammer wave velocity the superscript * represent the per-unit parameters n,n 1 ′-hydroturbine's mechanical speed and unit speed P t ,P 1 ′-power and unit power of hydroturbine ω m ,ω e -mechanical speed and electrical speed θ m ,θ e -mechanical/electrical rotation angle ω ms, ω es -mechanical/electrical synchronous speed f m ,f e -mechanical/electrical frenquency f ms ,f es -mechanical/electrical synchronous frenquency P t ,P e -mechanical/electrical power (active output) P c ,∆P e -certain/incremental generator ouput τ-gate openning η-turbine efficiency b p -permanent droop b t -temporary droop T d -reset time or dashpot constant T n -governor time constant T m turbine inertia time constant T w -water inertia time constant T y -main servo time constant T yB -pilot valve and servo motor time constant y d ,∆y,y 0 -turbine servomotor stroke and its deviation value,its initial value] J-inertia moment of the unit in direction of rotation (kg m2)

K P ,K I ,K d -proportional,

MODEL

A complex and coupled whole nonlinear system for hydopower plant is established as illustrated in Fig. 1. In order to develop the mathematical model, the system is decoupled into several modules and then dynamic nonlinear model is developed for each module. A new discussion can be carried out for the dynamic characteristics of the generator unit an palnt structure by using this numerical model.

Penstock Model

The the pressure conduit on the upstream of the plant is the connection between the hydro-turbine and reservoirs. The transition flow of the penstock can be described as Momentum Equation and Continuity Equation, as shown in Eq.( 1)~(2),

= 2 + + + DA Q fQ x H gA x Q A Q t Q ∂ ∂ ∂ ∂ ∂ ∂ (1) 0 = sin + + + 2 α A Q x Q Ag a x H A Q t H ∂ ∂ ∂ ∂ ∂ ∂ (2)
and, a Explicit Finite Differnetial Method (EFDM) is applied to disvretize the Eq.( 1) and Eq.( 2) at the junction nodes,as

) ( 2 ) + ( 2 1 = 1 1 + 1 1 + 1 + n i n i n i n i n i n i Q Q x ∆ A tQ ∆ Q Q Q - - - - - - -) ( 2 1 1 + n i n i H H x ∆ tg ∆ ) + ( ) + ( 8 1 1 + 1 1 + n i n i n i n i Q Q Q Q DA f - - (3)
) ( 2 ) + ( 2 1 = 1 1 + 1 1 + 1 + n i n i n i n i n i n i H H x ∆ A tQ ∆ H H H - - - - - -) ( 2 1 1 + 2 n i n i Q Q g a x ∆ A t ∆ α Q A t ∆ n i sin (4)
respectively, the subscripts i and superscripts n denote the position node i and time node n.

The hydro-turbine model consist of the relationship between H, Q, τ and n. It is the downstream boundary condition for penstock, the positive characteristic equation at time n is,

P a P P AH C A C Q - = (5) Where, a g C a = (6) ( )( )+ + = 2 L M M M M P Q Q Aa Q x ∆ A t ∆ A Q C - - ( )( ) ( ) t ∆ S S g H H Aa Q x ∆ A t ∆ H C f L M M M M a + ] + [ 0 - - (7)
The subsrcipts M and L denote the last two position nodes as shown in Fig. 2.

Hydro-turbine model

The Hydro-turbine converts the hydraulic power to the mechanical power which can supply the motivity to drive (as M t shown in Eq.( 8)) the shaft of rotor system,the rotation equation is, ( )

' , ' = 1 1 Q n f τ (9) ( ) ' , ' = 1 1 Q n g η (10
) Eq.( 9) to Eq.( 10) express the nonlinear relationship of the synthetic characteristic curves, which reflect the turbine's dynamic behaviors. With combining the H p , Q p which are solved from the penstock model as shown in Eq.(3)~Eq.( 5), the Eq.( 9) and Eq.( 10) and the law of hydraulic similitude, we can calculation the n, H n , Q p and M t for hydro-turbine at any operation point.

Governor model

The hydraulic governor system provide a reliable rotate speed regualtion of turbine with load variation in the power system. In this paper, a classical PID control strategy is applied in the coupled system on the basis of the predecessors' studies.

As shown in Fig. 3, the C f is ferenquency command. The transfer function from x to y in Fig. 3 can be written as, (12) Further,a equation of state is obtained from Eq.( 12) according to the Modern Control System [09], which is a first-order differential equations. The y can be sloved from the equations by Multistage Runge-Kutta algorithm.The gate opening can be obtained finally.

Rotor_bearing model

A FEM model of rotor system is establised in this paper by Ansys, which is a mature software for FEM methods. The beam188 element is employed to simulate the shaft of the rotor system. The rotor and the turbine are simplifed as mass21 element in the model. The combine14 element is used into simulate the guide bearings, and the dynamic characteristics of the bearings are equivalent to the coefficients of stiffness and damping (k ij, c ij ) of the element parameters.The rotor system model is mounted on a FEM model of powerhouse, as shown in Fig. 4.

l l l η h λ p h γ η p h η h ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ 6 = + 3 2 3 (14)
x=Rη l , λ=y/(L b /2), λ ∈ [-1,+1], the non-dimension parameters for the thickness of oil film is,

) cos( ) / ' 1 ( ) cos( + 1 = l b p l η β c c θ η ε h - - - - ) sin( / + l p η β b δ - (15)
Where, ε=e b /c b ,is eccentricity of the guide bearing, the b=c b /R b , is the ratio of clearance and radius. If the u and v is the circumferential velocity and radia velocity speed of the oil film, the Dynamic Reynolds Equation can be expressed as,

) cos + sin ( 6 + 6 = + 3 2 3 l l l l l η v η u η h λ p h γ η p h η h ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ (16)
Where,

2 + 4 3 = 2 2 l l η h h h η h a ∂ ∂ ∂ ∂
The oil pressure can be obtained from Eq.( 16) by FEM method, then the loads of the oil film can be generated through the interation of element solution [START_REF] Zhen-Yue | Dynamic Coefficients of Tilting Pad Guide Bearings of Hydraulic Turbine set[END_REF]. The loads' partial derivatives respect to the displacement and velocity are(i,s=x,y),

s f C s f K i is i is & ∂ ∂ ∂ ∂ = , = (17)
Which is the coefficents of the stiffness and damping for each pad.Then, the total coefficients (k ij ,c i j) of the bearing can be generated by combining the coefficients of pads [START_REF] Agnew | The governing of Francis turbines[END_REF].

Methods

The process of elecro-magnetic transient is not involved in the process of start-up, except the change of the field current during the process of voltage build up in the no-load operation, so, in this paper, the electrical system model is not taken into account. The rotor-bearing FEM model is built by Ansys, and other models are compiled in Fortran. With the UPFs of the Ansys, the programs of Fortran can be compiled into some commands from users, then, the operation process of hydro-turbine and generator's rotor can be simulated by the dynamic time-history analysis method with Ansys. On this basis, dynamic characteristics of the rotor system can be studied at any process during the operation of hydropower station.

The type of turbine is named HL180-LJ-410 in this paper, which is complied with the China's naming code. Its means that a francis turbine with runner type of 180, vertical arrangment,metallic spiral casing, and the diameters of runner is 410mm.

The mechanical eccentric forces of rotators and the unbalanced magnetic pull of rotor are considered for the kinetic rotor-bearing system [START_REF] Guo | The unbalanced magnetic pull and its effects on vibration in a three-phase generator with eccentric rotor[END_REF].

Result analysis

The datas of the whole model are listed in the Tab.1. During the process of start-up, the speed rised up smoothly with good dynamic response, short response time, small overshoot and fast attenuation. The gate opening fast reached the start opening at first, when the speed arrived at the rated point, the opening begins to attenuate with fluctuation with fluctuation until it stability at no-load opening. At the same time, in order to supply the power for growing speed, the turbine output increased rapidly. When the speed is stable, the output is reduced to 0 with a short fluctuation process.

When the gate is quickly opened at first, the head decreased and the flow increased quickly for the hydroturbine.Near the 12 seconds, the flow suddenly arrived The '*' represent the basic value of the per unit As the fig. 10 shown, the displacements of rotor and turbine are growing with the rotate speed. The trajectories are stable when the rotate speed at the rated point. When the process of voltage build up is started, the field current is increased. The rotor displace-ment is growing because of UMP's increasing as shown in Fig. 9, and with the UMP is balanced in fluctuation which is because the value of UMP is concerned with the rotor's displacement, the rotor's trajectories achieved a new stability again. As shown in Fig. 11, the axis trajectories' change of upper and lower guide bearings are similar to the rotor's because they provide the directly surpporting to the rotor.

(m) R a R b R r R i R o L p L r c b c r 0.9 2.12 5
Comparing Fig. 10 and Fig. 7 & 8, the vibration of rotor system has no influence on the hydraulic transient, because the head and discharge of turbine are stable when the rotor's displacement increased with the UMP.

Under the condition of this paper, the hydraulic force has not taken into consideration for the rotor system movement yet, but the hydraulic transient show its great influence on it. The generator power is zero during the process, so, when the turbine output is increased, the torque for the turbine is increased too, then the eccentric forces become larger with the higher rotate speed, finally, the rotor system vibrate with its ratotion.

Conclusions

In this paper, a whole model of the hydraulic power system was established. The model is used to studying the dynamic characteristics of rotor system when the coupled hydraulic-mechanical-electric factors is considered. The rotor-bearing system is the most important component in the model, meanwhile, the hydraulic, mechanical and electrical subsystems' are built and coupled together. The model is based on reliable mathematic and physical theories, and it gives a good simulation of the process of start-up as shown above.

On the basis of simulation, the results indicate that the hydraulic transient influence the rotor system 's rotation, it reflects the inner connection between the subsystems.

In the furture, the hdyraulic forces on the trubine and the impluse pressure on volute will be considered, a detial generator model will be coupled into the whole model, then, the more process of station operation can be simulated to study on the dynamic characteristics of rotor system.

Figure 1 .

 1 Figure 1. Sketch of Hydropower Plant System

Figure 2 .

 2 Figure 2. Boundary ConditonThe synthetic characteristic curves are applied into the numerical calculation which could consider the nonlinear variation sufficient during the ttransient process. The datas of the curves for whole operations can be obtained by fitting and interplation.

Figure 3 .

 3 Figure 3. The Transfer Function of the Frequency Regulation Mode

Figure 4 .Figure 5 .

 45 Figure 4. Rotor-bearing model A tilting pad guide bearing model is used in rotor system.The Reynolds Equation is used to describe dynamic oil film pressure fields around the pads in

Figure 6 .

 6 Figure 6. Simulation result of start-up transient

Figure 7 .

 7 Figure 7. The head change during start-up

Figure 8 .

 8 Figure 8. The discharge change during start-up at 0m 3 /s because of the closed gate, it caused a high head for the turbine, this phenomenon may create a great impulse pressure on the spiral case.

Figure 9 .Figure 10 .

 910 Figure 9. The change of UMP

Figure 11 .

 11 Figure 11. The trajectories of guide bearings' axis

  integral, derivative governor gain R a ,R r -radius of shaft and rotor R b -journal radius R i ,R o -radius of inner ring and outer ring for pad L p ,L r -height of pad, rotor length p-oil pressure e b ,e r -eccentricty of shaft axis, rotor centre c

b ,c r ,-clearance of bearing, rotor (air-gap length) h-thickness of oil film δ p -swing angle of pad α p -opening angle of pad η l -angle between the calculation location and y axis β-angle between a point of support for pad and y axis θ b -bearing axis deviation angle k ij ,c ij -stifness,damping cofficients,i=x,y;j=x,y M t ,M e -mechanical torque and electromechanical torque I f ,U f -field excitation current and voltage

Table 1 .

 1 Datas of models Nominal/no load value of Turbine

						Length/diameter of
	Head(m)	Weight (t)	Discharge(m 3 /s)	Power(kVA)	Speed(r/min)	penstock (m)
	117/124.9	420	86.65/14.7	88*	150*	495/8.5
			Time constants of governorn start-up (s)
	T m	T w	T d	T n	b p	b t	T y
	9	1.10	4T w	0.5T w	0.01	2(T w / T m )	0.1
			Parameter value of rotor-bearing system

ACKNOWLEDGMENTS

This research is supported by the National Natural Science Foundation of China (No.51379030).


Henrique C Zanúz
email: henrique.colao@gmail.com

Bruno Raffin
email: bruno.raffin@inria.fr

Omar A Mures

Emilio J Padrón

In-Transit Molecular Dynamics Analysis with Apache Flink

In this paper, an on-line parallel analytics framework is proposed to process and store in transit all the data being generated by a Molecular Dynamics (MD) simulation run using staging nodes in the same cluster executing the simulation. The implementation and deployment of such a parallel workflow with standard HPC tools, managing problems such as data partitioning and load balancing, can be a hard task for scientists. In this paper we propose to leverage Apache Flink, a scalable stream processing engine from the Big Data domain, in this HPC context. Flink enables to program analyses within a simple window based map/reduce model, while the runtime takes care of the deployment, load balancing and fault tolerance. We build a complete in transit analytics workflow, connecting an MD simulation to Apache Flink and to a distributed database, Apache HBase, to persist all the desired data. To demonstrate the expressivity of this programming model and its suitability for HPC scientific environments, two common analytics in the MD field have been implemented. We assessed the performance of this framework, concluding that it can handle simulations of sizes used in the literature while providing an effective and versatile tool for scientists to easily incorporate on-line parallel analytics in their current workflows.

INTRODUCTION

At the dawn of exascale computing, scientists are taking advantage more than ever of the evolving massively parallel computational resources available. Computational science simulations are increasingly generating more output data, due to the continuous improvements in the computational infrastructure and in the parallel numerical algorithms. Unfortunately, the I/O subsystem in HPC environments has not evolved at the same pace than processors and memories, so the storage and analysis of all that huge amount of data have started to become the real bottleneck in the scientist workflow.

Furthermore, there is an inherent need to simplify the creation of analysis code that can exploit all these heterogeneous resources. On the one hand, the need is to exceed the classic analytics pipeline based on a post-hoc approach to bypass the I/O bottleneck; and on the other hand, to get beyond a sequential analysis by leveraging the HPC infrastructure for enabling efficient analyses at large scale, but limiting the complexity of programming and deploying parallel algorithms at a relatively low level as required when relying on the classical HPC programming models (MPI/OpenMP, data partitioning and distribution. . . ).

, .

This paper explores an alternative approach based on Big Data frameworks. We use Apache Flink, a distributed streaming dataflow engine, to process in transit the data from the simulation. We leverage Flink high level stream processing programming model, and its runtime that takes care of the deployment, load balancing and fault tolerance. The analytics results are next written as soon as available in a scalable NoSQL database, using Apache HBase, a key-value database on top of HDFS. Both the results from the analytics and the raw data produced by the numerical simulation are stored in this database. Then, traditional post-hoc data processing based on batch analytics is also possible in our framework, as the output data from the simulation are stored after being analyzed in Flink. Besides, the same analysis kernels can be (re)used for the on-line in transit and the post-hoc data processing.

The numerical simulation used for the experiments is a Molecular Dynamics (MD) parallel mini-app, CoMD, that emulates the representative workloads of a real full scale MD simulation. At every given timestep, each CoMD process outputs the positions of all the atoms it simulates, giving what is commonly called an MD trajectory. Various analyses on these data are usually performed to extract meaningful knowledge.

The overall set-up consists in using multiple dedicated nodes of a cluster to run CoMD (simulation), Flink (analysis) and HBase (storage). Two common compute kernels for MD analytics were implemented to demonstrate the expressivity of Flink's programming paradigm: a position histogram of atoms and an atom neighbor computation based on a cutoff distance.

STATE OF THE ART

Traditional MD analysis tools associated with large MD simulation codes are mainly focused on post-hoc trajectory analysis. Gromacs [START_REF] Lindahl | GROMACS 3.0: A package for molecular simulation and trajectory analysis[END_REF] offers several trajectory analysis tools, for example. The included algorithms are mostly coded in C and only some of them are parallel. In terms of usability, each analysis is a different program. No framework for developing new analytics is provided, so whatever new analysis needed has to be coded from scratch. MD-Analysis [START_REF] Michaud-Agrawal | MDAnalysis: A toolkit for the analysis of molecular dynamics simulations[END_REF] or VMD [START_REF] Humphrey | VMD: Visual Molecular Dynamics[END_REF] are classical post-hoc trajectory analysis tools. They rely on MPI and/or OpenMP parallelization that are often complex to deploy at large scale and difficult to use efficiently by computational biologists. VMD has been extended to support in situ rendering [START_REF] Stone | High performance molecular visualization: In-situ and parallel rendering with EGL[END_REF].

Several in situ and in transit frameworks report experiments with MD simulations to demonstrate their scaling abilities. They all rely on HPC based approaches, mainly MPI+X, used in different ways [8-10, 22, 23]. To our knowledge these approaches mainly led to research prototypes and have not yet been adopted by scientists.

Computational biologists have considered the map/reduce model as an alternative to traditional HPC parallelization approaches to speed up the post-hoc MD trajectory analysis. Himach [START_REF] Tu | A scalable parallel framework for analyzing terascale molecular dynamics simulation trajectories[END_REF] was the first MD analysis framework to provide parallel execution capabilities inspired by Google's Map-Reduce. The authors initially considered Hadoop as a target, but due to its poor performance they developed a dedicated map/reduce framework based on MPI, Python and using VMD for the analysis kernels. Himach has a run-time responsible for assigning the tasks to the processors, coordinating the parallel I/O requests, storing and managing intermediate values (temporary key-value pairs) and orchestrating the data exchange between processes. Himach focuses on a temporal parallelization where each key-value is one simulation timestep with all atom positions. It shows a reasonable good performance and scalability, almost linearly until 32 MPI processes but degrading from that point due to I/O and communications.

On [START_REF] Paraskevakos | Taskparallel analysis of molecular dynamics trajectories[END_REF], the authors compare three general purpose task-parallel frameworks, Spark, Dask and RADICAL-Pilot, with respect to their ability to support post-hoc MD analytics on HPC environments. They also assess them in comparison to classical HPC MPI approaches. Their experiments show that Spark outperforms Dask when it comes to communication intensive tasks and iterative algorithms, due to the in-memory RDDs. Dask's low and high level APIs prove to be more versatile than Spark, although both have showed some limitations, requiring work-arounds to implement the analytics. Radical-Pilot proved to be more useful for coarse-grained task-level parallelism and when it is necessary to integrate other existing HPC analytics frameworks (such as MDAnalysis). Even though none of them outperformed their MPI counterpart, their easier programming paradigm and not so big performance gap still creates a strong case for using them for MD analytics.

If many tools have extended the map/reduce model for stream processing, very few works have attempted to combine this model with in situ analytics for large parallel simulations. The SMART [START_REF] Wang | Smart: A MapReducelike Framework for In-situ Scientific Analytics[END_REF] in situ framework proposes to rely on the map/reduce model to define in situ analysis. Implemented on top of MPI and OpenMP it outperforms Spark, but lacks features like support for in transit processing or fault tolerance. DataSpace [START_REF] Docan | DataSpaces: an Interaction and Coordination Framework forÂăCoupled Simulation Workflows[END_REF] is not a map/reduce framework but adopts the key/value storage concept to automatically index data produced by parallel simulations and store them in the memory of in transit nodes. Analytics query this in-memory database to get the needed data. DataSpace relies on MPI.

Many stream processing frameworks designed for internet and IoT applications exist, Spark, Flink and Storm being the most visible ones. For sake of conciseness refer to surveys like [START_REF] Chatterjee | Experimental Study on the Performance and Resource Utilization of Data Streaming Frameworks[END_REF] for comparative studies.

IN TRANSIT PARALLEL ANALYTICS FRAMEWORK FOR MD SIMULATIONS

The architecture of the proposed framework connects the CoMD simulation parallelized with MPI to Flink worker nodes using Ze-roMQ. Flink executes the analysis scripts in parallel and in-memory as soon as a new timestep is received, then injecting results to the HBase distributed database, that takes care of storing the results using its local disks (see an example of an 8 nodes configuration in Figure 1). An additional node is used to run both the Flink job manager and the HBase master.

Parallel Simulation

CoMD [START_REF] Exmatex | Exascale Co-Design Center for Materials in Extreme Environments[END_REF] is a mini-app for MD simulations, i.e. a simplified version of an actual MD application, sharing the same features and patterns in terms of operations, workload and work balance. We use the MPI version of CoMD, where each MPI updates the states of a fixed and different set of atoms at each timestep.

Distributed Stream Processing

Apache Flink [START_REF] Carbone | Apache Flink: Stream and Batch Processing in a Single Engine[END_REF] is an open source analytics engine supporting both stream and batch processing. Flink is a distributed, high-performing, highly availably and fault-tolerant framework. Flink is the result of the Stratosphere project [START_REF] Alexandrov | The Stratosphere Platform for Big Data Analytics[END_REF], an open source platform for massively parallel Big Data analytics. It includes various level of optimizations like Parallelization Contracts (PACTs) [START_REF] Alexandrov | MapReduce and PACT-Comparing Data Parallel Programming Models[END_REF]. To our knowledge Stratosphere/Flink has never been used for analyzing MD trajectories. Flink stream processing model relies on the dataflow model as defined in [START_REF] Akidau | The Dataflow Model: A Practical Approach to Balancing Correctness, Latency, and Cost in Massive-Scale, Unbounded, Out-of-Order Data Processing[END_REF]. It enables to continuously trigger the execution of map/reduce-like scripts [START_REF] Dean | MapReduce: Simplified Data Processing on Large Clusters[END_REF] on windows of data regularly extracted from data streams. If the window operates on keyed data, i.e. data that are partitioned according to a key defined by the user, the operations are executed in parallel.

A Flink script is transformed into a dataflow graph, taking into account the parallelization opportunities supported by the different operators being used. At execution, one TaskManager in each Flink worker node performs the operations over the input data. Each TaskManager provides one task-slot per core in the node, usually called sub-tasks. At specific points of the execution, the streams might be split into stream partitions that can be distributed or re-assigned to other sub-tasks. In opposite to a MPI based programming, the user does not need to explicitly control the distribution of data and computation.

Key-value Store

Apache HBase is an open-source distributed database based on Google proprietary BigTable. It is a non-relational database that runs on top of a distributed file system, commonly the Hadoop File System (HDFS). HBase can be seen as a distributed, sparse, persistent and multidimensional sorted map, indexed by a row key, a column key and a timestamp.

The worker nodes in HBase are called RegionServers. Each Re-gionServer contains an arbitrary number of regions. Each region is responsible for storing rows of one specific table, based on an interval of row-keys. The actual content of the rows are stored in HFiles on the underlying HDFS File System. An HBase master node coordinates the RegionServers and assigns their row-key intervals.

Putting the Pieces Together

Two connectors are needed in Flink to glue all the components together and get the whole system as depicted in Figure 1, with Flink receiving the data stream from CoMD, then processing the data and sending the raw data and the analysis results to the storage nodes. On the one hand, a source connector is used for the streaming data ingestion via ZeroMQ, and on the other one a sink connector is needed for dumping all data on the database.

The ZeroMQ [START_REF] Hintjens | ZeroMQ: Messaging for Many Applications[END_REF] (or 0MQ) asynchronous message library is used to connect the parallel simulation with Flink. ZeroMQ is datatype agnostic and provides its own sockets, which are able to send and receive atomic messages over different transport layers, such as TCP, multi-cast, inter-process and in-process communications. In this work, we rely on the 'Push-Pull' pattern.

Currently, Flink lacks an official ZeroMQ source connector (it was available in times of Stratosphere, but was not migrated to Flink due to license issues), so we have developed our own connector [START_REF] Mures | ZeroMQ connector for Apache Flink[END_REF]. The connector enables parallel connections: each CoMD process directly connects to a different thread on Flink nodes, these threads being distributed on the different Flink nodes and handling data to the local Flink sub-tasks. From Flink programming interface point of view, this is seen as a single stream that is actually distributed, effectively enabling parallel processing even before this data streams be keyed.

Flink includes an official HBase connector. So implementing the sink was more straightforward. The trajectory raw data are stored in one table where the rows are identified by a "MPI Rank_Timestep" key, so as to favor a partitioning of data based on their MPI rank. As most of the time queries need the data for a full timestep, this scheme reduces the risk for hotspots. The results from the analytics, way smaller than the raw trajectory, are stored in another table.

EXPERIMENTAL RESULTS

This whole environment runs on an homogeneous cluster with nodes having 2 processors Intel Xeon E5-2630 v3, with 8 Cores/CPU (i.e. a total of 16 cores per node), 128 GB, 2 x 558 GB HDD as storage and a 10 Gbps Ethernet network.

The first set of experiments used 8 simulation nodes (125 MPI CoMD processes, applying a 5 × 5 × 5 3D space partitioning), 1 Flink node with 16 task-slots and 1 HBase node. CoMD simulated a 32 million atoms model, which represents in number of atoms about half of the HIV virus molecular structure [START_REF] Perilla | Physical properties of the HIV-1 capsid from all-atom molecular dynamics simulations[END_REF].

The amount of streamed and processed data is controlled through the frequency of the timesteps output by CoMD: from every timestep down to one every 20 iterations. These trajectories are fully stored in HBase a couple of minutes after CoMD finishes. To limit data buffering throughout the workflow, low heap memory limits were imposed on Flink and HBase: 6 GB for Flink TaskManagers and 30 GB for HBase RegionServers.

These experimented saving frequencies are high compared to the standard of MD simulations. Real MD simulations iterate at high speed (250 it/s or more is common). Only a fraction of these timesteps are usually saved (1/500 or even less) to limit the trajectory size and the impact on the simulation performance [START_REF] Dreher | A Flexible Framework for Asynchronous in Situ and in Transit Analytics for Scientific Simulations[END_REF].

Three different on-line analysis scenarios were used for each experiment: a simple raw data injection, with no analysis at all; a low-demanding data processing that computes a position histogram of the atoms; and a high demanding analytics that computes the list of neighbors for each atom, i.e. at a distance lower than a given cutoff value. Histogram computation and neighbor identification are common analysis patterns for MD analytics [START_REF] Paraskevakos | Taskparallel analysis of molecular dynamics trajectories[END_REF]. The Flink pseudocodes for both algorithms are shown in Figure 2.

ZeroMQ ensures that messages stay in order and that none is lost. But we have to make sure that data processed together contain all and only the atoms from the same timestep. For that purpose we rely on Flink time event windows. The simulation timestep value contained in every message sent by CoMD is used as watermark. Flink splits the data streams according to tumbling windows of size one, i.e. once all messages for the next timestep have been received, ensuring the integrity of the data processed for each window.

For both algorithms we rely on Morton indexing (z-order curve), used to partition the space for the position histogram or as acceleration data structure when searching for neighbors.

For the position histogram, the data from the messages are first split (locally with a flatmap), associating with each atom the corresponding z-index. Because the data stream is parallel (one data source per MPI process) this operation is executed in parallel on Flink workers. Then a keyby shuffles the data so that all atoms belonging to the same z-cell are gathered together on the same Flink worker. This operation is combined with a reduction to count the number of atoms per cell. Again this shuffling/reduce operation is performed in parallel.

For the neighbor search, we use a z-indexing defined according to the cutoff radius. A flatmap associates with each atom a z-index for its z-cell but also for each z-cell it is neighbor of. This leads to duplicate the atom data 27 times for a 3D simulation like CoMD. Then the data are shuffled to gather all atoms with the same z-index together through the keyby, leading to gather all atoms belonging to the same neighborhood together. Then, a local operation (apply) computes the neighbor list for each atom.

The first experiment looked at analyzing the stability of the framework to cope with a state-of-the-art input pressure. Since MD simulations might run for extensive periods and generate huge trajectories, we want to be sure the system can reach and sustain a steady state, with non-increasing demand for buffering resources along time. Therefore, we tested our workflow with different frequencies of data outputs.

The maximum sustainable throughput is obtained at 71 MB/s (or 8.875 MB/s/node) when CoMD sends data every 11 iterations, for the direct data injection (Flink simply forwards the data to HBase) and the histogram computation. For both experiments CoMD produces Figure 3 shows some interesting data about the behavior of the system with the maximum experimented sustainable throughput (Figure 3a) and with the minimum unsustainable one (Figure 3b). In both subfigures, the graph on the left shows the cumulated number of messages processed throughout time at each point of the system. In this graph, three cumulative histograms are overlaid. The orange one counts the number of messages generated by CoMD; the blue one, the number of messages received by Flink and the green one, the number of messages handled by HBase. In the Figure 3a CoMD produces data every 11 iterations. The system is able to handle the data at the same rate it is produced by the simulation: CoMD generation, Flink processing and HBase ingestion histograms are almost perfectly superimposed. On the other hand, setting CoMD to produce one output every 10 iterations (Figure 3b) leads to a growing gap between the CoMD histogram and the almost matching Flink and HBase histograms. After about 200s Flink and HBase start to lag behind CoMD. After investigation, it appears that HBase is not able to cope with the data pressure. HBase imposes back pressure to Flink very quickly (both histograms almost match) that in turn imposes back pressure on ZeroMQ. The different buffers involved in this chain from HBase to ZeroMQ (ZeroMQ maintains buffers on the sender and receiver) get filled but are large enough to avoid impacting CoMD execution that proceeds at the same pace. After 750s CoMD finishes to produce messages while ZeroMQ, Flink and HBase keep on working for some extra time to empty the buffers. The right curves show that HBase accumulates data into buffers and flushes them at a regular pace.

The neighbor computing analysis leads to a significantly lower sustainable throughput of 1.10 MB/s. The bottleneck is clearly the overload of Flink as the data are replicated 27 times and need to be shuffled when keyed. We are exploring alternative algorithms that hopefully will achieve significantly better performance, addressing the issue of data duplication and optimizing the neighbor computing.

We also ran a weak scalability test, scaling both the workload and the number of processors, running the histogram computation on a 4 times bigger set-up, with 512 CoMD processes on 32 nodes and 64 cores for Flink and HBase (4 nodes each). The number of atoms in the simulation was also multiplied by 4, generating a trajectory of 256 GB. The maximum sustainable throughput achieved reached 100 MB/s, only 40% more than with a four time smaller configuration.

For situating these performance according to execution in a more traditional HPC context let consider the results of [START_REF] Dreher | A Flexible Framework for Asynchronous in Situ and in Transit Analytics for Scientific Simulations[END_REF], running a Gromacs MD simulation on a supercomputer. They test different in situ processing scenarios with FlowVR, saving simulation outputs to disk or performing some analytics. The closest scenario where all data produced are saved to disk, the simulation produces about 35 MB/s of data when running on 256 cores, and 75 MB/s when running on 512 cores. This is below the sustainable bandwidth we achieved with our set-up (this comparison should obviously be considered with care as it is not a direct comparison of the very same simulation running on the same machine).

CONCLUSION AND DISCUSSION

This paper investigates the HPC & BigData convergence for the on-line analysis of parallel simulation outputs. We presented early results using Flink to process in transit data produced by a parallel MD simulation and store the raw data from the simulation as well as the analytics results to HBase. We managed to achieve a sustainable throughput of 71 MB/s for 125 simulation processes (resp. 100 MB/s with 512 simulation processes) while saving 63 GB of raw results to HBase (resp. 256 GB) and computing a histogram. It also revealed performance issues when scaling or with the neighbor computations, but that should not be considered definitive as this early work leaves room for many possible improvements.

We are currently working on refining the analytics algorithms, using different Flink operators for increasing performance, moving to a different machine with InfiniBand network and SSD disks. We are also considering using a different storage. HBase was chosen as being standard in the Big Data domain, but other storage systems like Cassandra are known to have significantly better write performance.

This work also shows that Flink offers both a high level programming model easy to master, and an advanced runtime discharging the user of many time consuming aspects such as data communication, data partitioning, load balancing and fault tolerance. The programmer does not have to control data localization, i.e. to know explicitly the data each node keeps. He/She only has to control data aggregation by relying on data indexing with keys. Once the data are stored on HBase, the user can also query and process these data post-hoc without changing of programming environment. An HPC based approach would require more expertise and development time, often for the benefit of better performance. We discussed our approach with a few computational biologists that find it very appealing to trade some performance for saving on human expertise and time.

Figure 1 :

 1 Figure 1: Diagram of the proposed framework deployed on 8 nodes: the MD simulation (CoMD) running in 4 nodes; 2 nodes executing the on-line analytics with Flink, receiving the data stream through ZeroMQ; and other 2 nodes in charge of storing data (raw data from the simulation and results from the analytics) in a distributed database (HBase).

Figure 2 :

 2 Figure 2: Analysis kernels used in the experiments

  Sustainable throughput: 8.875 MB/s/node (for a total of 71 MB/s) Unsustainable throughput: 9.75 MB/s/node (for a total of 78 MB/s)

Figure 3 :

 3 Figure 3: Analyzing the stability of the framework for two different throughputs with the Histogram analytics

ACKNOWLEDGMENTS

Experiments ran on the Grid'5000 testbed, supported by a scientific interest group hosted by Inria and including CNRS, RENATER, several Universities and other organizations. This work is partially

A ARTIFACT DESCRIPTION APPENDIX: IN-TRANSIT MOLECULAR DYNAMICS ANALYSIS WITH APACHE FLINK A.1 Abstract

This section describes the artifacts of this work submitted to ISAV '18. We describe the environment and the context of the experiments presented in Section 4. We explain how we compile and execute CoMD and how we deploy and execute the in transit analytics pipeline proposed in the paper.

A.2 Description

A.2.1 Check-list (artifact meta information).

• Algorithms: Morton indexing (z-order curve), position histogram, neighbor identification based on a cutoff distance • Programs: C binary, Java jars, Python scripts, CoMD, ZeroMQ, JZMQ, Flink, HBase • Compilation: C, MPI, Java, Maven • Data set: Generated by CoMD • Run-time environment: Flink and HBase appropriately configured with access to JAVA_HOME. Java bindings for ZeroMQ installed (JZMQ). • Hardware: Grid 5000's paravance cluster; each node: 2 x Intel Xeon E5-2630 v3, 8 cores/CPU, 128 GB, 2 x 558 GB HDD, 2 x 10 Gbps • Publicly available?: Not yet A.2.2 How software can be obtained (if available). The framework we are presenting in this paper is not publicly available yet, but will be open-sourced in a near future. The whole software stack will be included in a git repository, including an appliance to ease the setup and deployment. For more information, feel free to send a request at emilio.padron@udc.gal.

A.2.3 Hardware dependencies.

The only hardware dependency is a cluster with nodes connected through an Ethernet network. Grid 5000 * (G5K), a testbed for experiments on high performance computing, distributed systems, big-data and cloud, was used in this work. G5K infrastructure is composed of several clusters on several sites throughout France and Luxembourg. In all the experiments, all the nodes used were equal and part of one single cluster, in order to avoid introducing another parameters that can alter the results. The specifications of the nodes used for all the experiments are: 2 x Intel Xeon E5-2630 v3, 8 cores/CPU, 128 GB, 2 x 558 GB HDD, 2 x 10 Gbps.

A.2.4 Software dependencies. Our software stack is composed of: (4) Execute the modified version of CoMD to generate the data stream to be analyzed (5) Collect experimental results (see A.5)

A.5 Evaluation and expected result

By performing the procedure above, the experiment is launched. Once CoMD finishes and all the data is persisted in HBase, it is possible to obtain (from the lines in HBase) the timestamps of creation, ingestion in Flink and storage in HBase for all the messages sent by CoMD. Computing a cumulation histogram counting the number of each type of event throughout the experiment execution results on the graphs we have shown in this paper. We consider a certain throughput value to be stable when no gap between the three curves (messages created, ingested by Flink and stored) is observable in any of the repetitions made. Otherwise, it is classified as unstable.

A.6 Experiment customization

By keeping the same scripts for the experiments, it is still possible to customize them by tuning a wide set of parameters (as shown in A.4), such as: the number of atoms computed; the number of nodes dedicated for each role (simulation, analysis, storage); the number of iterations the simulation computes; the output interval. . . These are some of the specific experiment invocations we have used in this paper:

A.7 Notes

To enable experiment reproducibility, we relied on different software tools that enabled us to control the full deployment and execution of experiments, from the OS up to the last lever software used. Due to the strong coupling that these tools has with G5K and its batch scheduler, OAR * * , we maintain a specific git repository with all the deployment stuff, in addition to the main repository with the Flink analysis kernels. This repository will be open-sourced as well, along with the main one. The use of these tools is not actually a strong requirement to deploy and run the experiments, but it is recommended mainly by two reasons:

• The first one is building the system image that will be used on the experiment nodes, allowing us to precisely control what is installed on them. This is done through the use of Kameleon Image Builder * * * , a tool that allows the creation on systems images based on scripts. • The second reason is to completely automatize the experiments workflow. For that a python script was made. This script uses Execo * * * * , a G5K python API that allows running remote commands over SSH and lower level commands, such as requesting and deploying nodes. In more details, this python script was used to set-up the configuration of the software stack, to configuring the clusters for Flink, HBase and CoMD, to launch the experiments and to recover and log the interesting data from the experiments in a result file. This script receives as input a YAML file which describes which experiments will be executed, their parameters and the amount of nodes dedicated for each role in this environment. Regarding the execution of the experiments, this python script was executed on a separate node from the ones used on the setup. This node runs the regular G5K system image, not the ones customized by us. After being launched, this script is the one in charge of reserving the master and slave nodes for the experiments. As explained on the paper, the slave nodes were used to run the CoMD MPI processes, the Flink's TaskManagers (workers) and HBase's RegionServers (workers). The master node is used to launch the MPI processes on the other nodes, run Flink's JobManager and HBase's Master process. Even though the master node was assigned with several tasks, those task were not compute-intensive. In addition, it has not showed any symptom of being overloaded during our experiments, staying with low memory and CPU usages.

A.8 References