Edson Luiz Padoin
email: padoin@unijui.edu.br

Philippe O A Navaux
email: navaux@inf.ufrgs.br

Jean-François Méhaut
email: jean-francois.mehaut@imag.fr

Using Power Demand and Residual Load Imbalance in the Load Balancing to Save Energy of Parallel Systems

Keywords: Energy Consumption, Load Balancing, DVFS, Power Demand, Energy Saving

The power consumption of the High Performance Computing (HPC) systems is an increasing concern as large-scale systems grow in size and, consequently, consume more energy. In response to this challenge, we have develop and evaluate new energy-aware load balancers to reduce the average power demand and save energy of parallel systems when scientific applications with imbalanced load are executed. Our load balancers combine dynamic load balancing with DVFS techniques in order to reduce the clock frequency of underloaded computing cores which experience some residual imbalance even after tasks are remapped. The results show that our load balancers present power reductions of 7.5% in average with the fine-grained variant that performs per-core DVFS, and of 18.75% with the coarse-grained variant that performs per-chip DVFS over real applications.

Introduction

The main focus of HPC systems has been performance, and nowadays the HPC community works toward building Exascale systems, which will provide unprecedented computational power, allowing to solve even larger scientific problems. However, if we conceive Exascale supercomputers by scaling the current technology, these systems could demand over a Gigawatt [START_REF] Hsu | Towards efficient supercomputing: A quest for the right metric[END_REF][START_REF] Feng | The green500 list: Encouraging sustainable supercomputing[END_REF]. In this way, a global research effort has risen to try to break this barrier while avoiding such high power demands.

Reductions in the total execution time of real applications are also relevant in the perspective of energy consumption, as energy is saved when hardware resources are used for a shorter time. However, it is possible to achieve even greater energy savings if the runtime system is able to exploit the residual imbalances to fine tune the voltage and frequency of cores accordingly. In this case, a challenge lies in reducing the energy consumption of the application while maintaining a similar performance.

Scientific applications generally have an irregular computational load and tasks with different processing demands, which makes the load distribution among the processors more difficult. In this context, different approaches have been used to mitigate these imbalanced workloads and achieve an efficient use of all available parallel resources. However, most of these load balancing strategies focus only on reducing the execution time using only computational load, neglecting power demand and total energy consumption.

Generally, approaches to reduce the power demand are used separately of load balancing research for increasing application's performance. Several proposal have used DVFS, however this technique may cause performance degradation and an increase of the total execution time of parallel applications. Other proposals have used load balancing strategies to mitigate the imbalanced workloads and achieve an efficient use of all parallel resources, which reduce the overall execution time, and, consequently, save energy. Our developed approaches make the case that both strategies, load balancing and DVFS, could be used together over real applications to save energy. They identify the possibility of reducing the processors clock during the runtime to achieve better gains over other load balancer algorithms. In this form, energy improvements are achieved due the reduction of average power demand during the execution and also through of reducing of the execution time by reducing the amount of tasks migrated.

In this paper we propose and validate over real applications our energy-aware load balancers, which minimize the total energy consumption when applied over iterative applications with imbalanced loads. We show that, combining load balancing with DVFS, our energy-aware load balancers are able to reduce the energy spent with load balancing when they are used over other state-of-the-art algorithms.

The remaining sections of this paper are organized as follows. Section 2 discusses some of the related works on DVFS and energy-aware load balancing. The main concepts of our energy-ware load balancer (EnergyLB) as well as its two variants and implementation details are discussed in Section 3. In Section 4 we present the evaluation methodology and the real applications used in the conducted experiments. In addition, in Section 5 we address the results obtained from the experiments. Finally, the Section 6 emphasizes the scientific contribution of the work and notes several challenges that we can address in the future.

Related Work

Different approaches are used to reduce the power demand of the HPC systems aim improve their energy efficiency. In this way Dynamic voltage and frequency scaling (DVFS) has been used in different contexts as a means to save energy.

Spiliopoulos et al. [START_REF] Spiliopoulos | Introducing DVFS-management in a full-system simulator[END_REF] extended the Gem5 simulator to support full-system DVFS modeling and to study the behavior of different DVFS governors. They concluded that the interactive governor is faster than on-demand to adapt to the workload changes and thus achieves better performance at about the same energy consumption. Gerards et al. [START_REF] Marco | Analytic clock frequency selection for global DVFS[END_REF] analyze the use of global DVFS in the context of multi-core processors. They proposed a theoretical method to transform the problem of finding an optimal clock frequency on global DVFS systems to a single core problem by using the amount of parallelism of applications. Their main goal is to minimize the energy consumption of nontrivial real-time applications. Kin et al. [START_REF] Kim | Energy-centric DVFS controlling method for multi-core platforms[END_REF] proposed a realistic DVFS performance prediction method and a practical DVFS control policy (eDVFS) that aims to minimize total energy consumption in multi-core processors. Their experimental results show that eDVFS can save a substantial amount of energy compared with Linux ondemand. Isci et al. [START_REF] Isci | An analysis of efficient multi-core global power management policies: Maximizing performance for a given power budget[END_REF] proposed to fine-tune the processor's clock frequency by using workload characteristics aim to maintain a chip-level power below a specified budget without degrading the performance significantly.

Other challenging problem that has been studied extensively to improve the performance of parallel applications is Energy-Aware Load balancing [START_REF] Kalé | Load balancing in parallel molecular dynamics[END_REF][START_REF] Zheng | Periodic hierarchical load balancing for large supercomputers[END_REF]. However, few works have made some efforts to further improve the energy consumption. Aupy [START_REF] Aupy | Energy-aware scheduling under reliability and makespan constraints[END_REF] proposed energy-aware scheduling models to schedule tasks under reliability and makespan constraints. They designed and evaluated them using simulations with different heuristics based on the failure probability, the task weights, and the processor speeds. Sarood et al. [START_REF] Sarood | A 'cool' way of improving the reliability of HPC machines[END_REF] proposed a load balancing strategy that limits the processors' temperatures to reduce the energy spent in cooling and to prevent hot spots. Their results achieved energy savings of up to 63%, with a timing penalty from 2% to 23%. Goel et al. [START_REF] Goel | Portable, scalable, per-core power estimation for intelligent resource management[END_REF] proposed a model that uses CPU performance counters and CPU temperature to generate accurate per-core power estimates in real-time. They showed that the model can be used to guide scheduling decisions in power-aware resource managers. Hartog et al. [START_REF] Hartog | Mapreduce framework energy adaptation via temperature awareness[END_REF] studied the relationship between CPU temperature and energy consumption in clusters and provided a method of estimating the power consumption of the system. This method was then used to implement a MapReduce framework that can evaluate the current status of each node and dynamically react to estimated power usage without having to rely on readings from expensive power monitoring hardware affixed to each node in the cluster.

As opposed to these works, our approach performs load balancing along with DVFS to improve the performance and to reduce the energy consumption by exploiting residual imbalances of parallel applications. In addition, we also reduce the cost of task migrations, since we only migrate tasks between processors when necessary. The performance and energy efficiency of our approach were also analyzed on a set of real-world application running on top of a real platform without the need of simulations.

3 Energy-Aware Load Balancer (EnergyLB) Scientific applications generally have tasks with different processing demands. Load balancing approaches have been used to mitigate these imbalanced load by making a better load distribution among the available processors in the system. Load balancing algorithms try to mitigate load imbalance by moving tasks between processors, which incurs in overheads to data collection, decision making and mainly tasks migration, and may affect the application runtime and consequently the total of energy spent. However, the most of the strategies aim only on reducing the execution time. Only a few recent strategies began to use the power demand and energy consumption information in their decision making process.

Thus, when imbalanced applications are executed on parallel platforms, the program only finishes its execution when the most loaded processor completes all its tasks. This way, the total execution time of an application is generally defined by the most loaded processor. This means that all other processors or cores in charge of lightweight tasks finish first and remain consuming energy without performing any actual work for the application.

In the case of iterative applications, energy savings can be achieved if the clock frequency of these underloaded processors or cores is decreased through DVFS in such a way that they can still and their tasks before or at the same time of the most loaded ones. Load balancers can be applied to reduce such imbalances but some residual imbalance may still remain after task migrations, making room for further energy improvements.

In principle, DVFS can be performed mainly at two levels of granularity. Systems featuring per-chip DVFS use the same power delivery network to reach every core, and consequently, bind each core to the same DVFS schedule. On the other hand, systems that allow per-core DVFS control use a separate voltage regulator for each core and therefore allow every core to have an independent DVFS schedule. Taking these two levels of granularity in mind, we propose two variants of Energy-Aware Load Balancer (Figure 1). The first one, is a centralized approach, called Fine-Grained Ener-gyLB ( FG-EnergyLB). It is suitable for platforms composed of few tens of cores that allow per-core DVFS. It considers a flat view of the underlying platform, balancing the load among all cores and performing DVFS on each individual core. When the load imbalance is considerably high, we rely on currently available load balancers to better distribute the tasks among the cores. On the other hand, when migrating tasks is not beneficial, we adjust the clock frequency of underloaded cores individually to save energy. Although this approach allows for a finer tune of the load and frequency, it becomes subject to considerable overheads on platforms having hundreds of cores. To overcome this issue, the second, called Coarse-Grained EnergyLB ( CG-EnergyLB), uses a hierarchical approach suitable for current HPC platforms, where the underlying platform is represented as a two-level tree. In this scheme, the platform (root) is composed of several multi-core processors (leaves). Our hierarchical approach performs load balancing to distribute the load among the processor's cores in each leaf. Then, in the root level, it adjusts the frequency of each multi-core processor according to its load to save energy.

Beyond the load balancers, we implemented two other modules. The first, called of Energy Daemon is responsible for gathering power and energy measurements through the Intel's Running Average Power Limit (RAPL). And the second is a Lib Energy that offers an interface to get system information such as the maximum frequency that can be set to a core or processor, to set the clock frequency of individual processors/cores or to set a specific Linux governor.

Our strategies take three input parameters in their execution. The first one, is the load balancer that is used to migrate tasks when the imbalance is high. The second one, is the maximum frequency available by processors that can be set to a core, and the last one, is a threshold value, used to decide whether call the load balancer or perform DVFS strategy.

In this way, in each call of the load balancer, is verified if the weighted load of each processor exceeds or not the threshold, make decisions to adjust the frequencies (determining so that the frequency will be decreased or increased) or invoke other load balancer to migrate tasks.

We used the Charm++ parallel runtime system [START_REF] Laxmikant V Kalé | Programming petascale applications with CHARM++ and AMPI[END_REF] to implement the two variants of EnergyLB. This framework provides a simple API that can be used to implement load balancers without changing the source code of the applications. The instrumented data load can be dynamically obtained during the application execution through the use of the API.

4 Evaluation Methodology The platform runs an unmodified SUSE Linux Enterprise Server operating system with kernel 3.0.101-0.29 installed. All applications as well as the Charm++ programming model were compiled with GCC 4.8.2. The Charm++ version used in our experiments was linux64-6.5.1. The results presented in Section 5 are the average of at least 10 runs. The relative error was less than 5% using a 95% statistical confidence by Student's t-distribution.

Applications

To evaluate the runtime, power demand and total energy consumption of our Energy-Aware Load Balancers, we selected three real-world applications. They were chosen due to their varied range of communication patterns and workload characteristics. The description of the applications is given below:

• Ondes3D is a seismic wave propagation simulator employed to estimate the damage in future earthquake scenarios [START_REF] Dupros | Exploiting intensive multithreading for the efficient simulation of 3d seismic wave propagation[END_REF]. In Ondes3D, seismic waves are modeled as a set of elastodynamics equations. These equations are then solved by applying a finite difference method. The application is over decomposed into multiple virtual MPI processes per core.

Ondes3D presents load irregularity due to the boundary conditions producing additional work, and load dynamicity from the simulation of waves spreading through space;

• Lulesh simulates a variety of science and engineering problems require modeling hydrodynamics, which describes the motion of materials relative to each other when subject to forces. The Livermore Unstructured Lagrange Explicit Shock Hydrodynamics (LULESH) application was originally developed as one of the five challenge problems in the DARPA Ubiquitous High Performance Computing (UHPC) program. Lulesh solves one octant of the spherical Sedov problem using Lagrange hydrodynamics [START_REF] Ss Dosanjh | Exascale design space exploration and co-design[END_REF]; and

• Lassen is a mini-application that explores applications with varying load balance. It simulates the propagation of a wave through as it travels around an unstructured mesh and it works on a 2D and 3D version. The computational load of the mesh is sub-divided into domains, which are assigned to processors. Therefore, this becomes a challenging problem to effectively parallelize since the workload is constantly changing [START_REF] Mccandless | Lassen mini-app[END_REF].

Input parameters

Table 1 summarizes the characteristics of the applications and parameters used in our experiments. Different load balancing frequencies have been chosen for different applications in order to strike a balance between the benefits of remapping tasks and the overheads of moving tasks and computing a new task mapping.

Load Balancers

We applied our energy-aware load balancer along with different standard load balancers availin Charm++. The rationale behind that is to evaluate the energy savings that can be obtained by exploiting residual imbalances of these load balancers. To compare the power demand and energy consumption results, we selected the load balancers GreedyLB (G) and GreedyCommLB (GC) that have a greedy algorithm, RefineLB (R) and RefineCommLB (RC) that use a limit of task migrations, RandCentLB (Ra) that randomly assigns objects to processors and ScotchLB (S) that map based on the recursive bi partitioning using both the source process graph and the target architecture graph.

Experimental Results

Real applications as Ondes3D, Lulesh and Lassen have different initial imbalance and different dynamicity in their computational load, which leaving space for the use of load balancing algorithms. In this context, some energy improvements were achieved.

Thus, this section provides an analysis of total energy spend with load balancing over total energy consumption and the percentage of energy spend with load balancing when FG-EnergyLB and CG-EnergyLB are used over these real applications. The application runtime depends on several issues, among them, the number of parallel tasks and their load, the duration of each timestep, and the selected load balancing strategy. The impact of load balancing is directly related to the load balancing frequency once load balancing overhead can overcome the gains achieved with load balancing.

Fine-Grained EnergyLB Evaluation

The propagation of the wave of this application introduce a dynamism to the load. In this way, Ondes3D has both irregularity and dynamicity in its execution. Applying FG-EnergyLB alone, it reduces the clock frequency of cores used during the runtime, which reduces the power demand in 15% on average, from 35.5 W to 30.2 W. By reducing the processors power according their computational load, FG-EnergyLB achieve saving energy of 13.9%.

Similar to Ondes3D, Lulesh and Lassen have imbalance and dynamicity of load. Considering this, FG-EnergyLB make several DVFS to adjust the clock frequency of the cores according their relative load, reducing the average power of 35.06 W to 33.64 W, to Lulesh and of 34.8 W to 30.8 W (11.48%).

When FG-EnergyLB was applied over the residual imbalance of other load balancing algorithms the average power demand is reduced in up to 5.29%, 4.05% and 5.23% respectively.

Considering this reductions in average power demand, we select the energy spent with load balancing when load balancers are used alone (LB) and the energy spent with load balancing when load balancers are used with FG-EnergyLB (LB+FG) of the total energy spent in the application execution (App), as depicted in the Figure 3. The energy spent with load balancing comes from of the time overhead of every load balancer call. In the test with Ondes3D, the load balancers are responsible in average for spending of 7.87% of the total energy consumed during execution of the application (red column in Figure 3(a)). When FG-EnergyLB was used together with load balancers (green column), it is able to saving energy of own load balancing. The energy spent with load balancing was reduced from 7.87% to 2.07%, on average. This occurs due the reduction of the total execution time and mainly by reductions on average power demand.

Over Lulesh and Lassen applications, the load balancing is responsible for 4.91% and 3.57% respectively of the total energy consumed during execution of the application. Using FG-EnergyLB together with other load balancers, the total energy spent with load balancing is reduced to 3.52% and 3.13% respectively, as shown in the Figure 3(b). In this way, using FG-EnergyLB together with other load balancers, the average of energy spent with load balancing is reduced from 5.45% to 2.90%.

Coarse-Grained EnergyLB Evaluation

This section provides a CG-EnergyLB evaluation over real-world applications. Each one these scientific applications present a different and dynamic behavior, the load of its tasks change through the iterations, which provides a more challenging scenario for energy aware load balancing. Since all the 192 cores will be used, different parameters are used in the evaluation of the CG-EnergyLB, as shown in the Table 1.

These applications has irregularity in the creation of its tasks and load dynamicity during execution. So, when load balancing algorithms are applied over real applications, different improvements are achieved. By updating the frequency of the processors according to their weighted loads, When CG-EnergyLB is employed alone, it reduces the average power demand of the system while the application is running. Using this approach, reductions of 27.4%, 25.5% and 9.17% on average in power demand are achieved, as depicts the Figure 4. Once applied load balancers, residual imbalances still are present. In this way, when CG-EnergyLB is employed to mitigate this residual imbalance left by the load balancing algorithm, it identifies the possibility of reducing the processors clock to achieve better gains over these algorithms. For theses applications, the average power demand is reduced in 21%, 20.6% and 8.83% respectively.

Considering the different power reductions, we analyse the load balancing energy when load balancers are used alone (LB) and the load balancing energy when load balancers are used with FG-EnergyLB (LB+FG). These energy value are depicted in the Figure 5 compared to total energy spent in the application execution (App).

We realize that our CG-EnergyLB has load balancing overhead lower than other stateof-the-art load balancers when tested with real applications. So, when the application call CG-EnergyLB, it spend less time to adjust the clock frequency, which results in a reduction in total execution time and an equivalent saving energy.

In the tests with the Ondes3D, Lulesh and Lassen applications, the load balancers are responsible in average by spent of 28.87%, 9.87% and 6.80% of the total energy consumed during execution of the applications. When CG-EnergyLB was used together with load balancers, it reduces the energy spent with load balancing on average to 24.55%, 8.46% and 5.71% respectively.

Conclusions

The exponential increase in power consumption related to a linear increase in the clock frequency and a higher complexity involved in the processors' design changed the course of development of new processors. Power consumption has become a critical aspect to the development of both large and small scale systems. This concern is now enough to warrant the research on techniques to improve the energy efficiency of parallel applications running on top of HPC platforms.

In this paper, we focused on reducing the energy consumption and power demand of parallel applications without considerably degrading their overall performance. Our main idea was to exploit the existence of residual imbalances on iterative applications to adjust the clock frequency of underloaded cores/processors through DVFS. We proposed two variants of Ener-gyLB, an energy-aware load balancer and validate both variants in Charm++. Both strategies can be used together with any other load balancer to improve energy efficiency of real world applications.

Experimental results present improvements in energy consumption when our EnergyLB load balancers were used with state-of-the-art algorithms over real-world applications. FG-EnergyLB and CG-EnergyLB are able to reduce the load balancing energy in average of 32.6% and 11.1% respectively. This work can be extended in different directions. One possibility would be to develop a new load balancer that performs load balancing and DVFS at the same time in each load balancing step. For that, it would be necessary to create a heuristic that take into account the cost of task migrations between cores/processors that operate in different clock frequencies. Another possibility would be to develop a hierarchical energy-aware load balancer that performs task migrations between cores of the same processor and only migrate tasks between processors when needed. In this scheme, only the processors involved in task migrations would need their clock frequencies to be adjusted, reducing overhead of performing DVFS on all processors at each load balancing step. Finally, we also intend to evaluate the benefits of FG-EnergyLB and CG-EnergyLB on other real-world applications and platforms.

Figure 1 :

 1 Figure 1: EnergyLB Implementation

Figure 2 :

 2 Figure 2: Average Power Demand reductions for each real application with FG-EnergyLB.

Figure 3 :

 3 Figure 3: Energy spend with load balancing over total energy consumption for each real application with FG-EnergyLB.

Figure 4 :

 4 Figure 4: Average Power Demand reductions for each real application with CG-EnergyLB.

Figure 5 :

 5 Figure 5: Energy spend with load balancing over total energy consumption for each real application with FG-EnergyLB.

  Our experiments were conducted on an Altix UV 2000 platform designed by SGI. The platform is composed of 24 NUMA nodes. Each node has an Intel Xeon E5-4640 Sandy Bridge-EP x86-64 processor with 8 physical cores running at 2.40 GHz. There are 14 clock frequency levels available in this processor, allowing us to vary the clock frequency of the processor from 1.2 GHz up to 2.4 GHz.Each core of the Intel Xeon E5-4640 has 32 KB instruction and 32 KB data L1 caches and 256 KB of L2 cache. Each node has 32 GB DDR3 memory, which is shared with other nodes in a cc-NUMA fashion through SGI's proprietary NUMAlink6. Overall, this platform has 192 physical cores and 768 GB DDR3 memory.

	4.1 Experimental Environment

Table 1 :

 1 Summary of the input parameters of applications.

	Application FG-EnergyLB CG-EnergyLB Iterations LB Frequency
		Tasks	Tasks
	Ondes3D	128	1024	500	20
	Lulesh	729	5832	1000	50
	Lassen	256	512	550	50

Acknowledgements

This work was supported by CNPq, CAPES, FAPERGS and FINEP. The research has received funding from the EU H2020 Programme and from MCTI/RNP-Brazil under the HPC4E Project, grant agreement number 689772. Work developed on the context of the associated international laboratory between UFRGS and Université de Grenoble -LICIA.


Sudhakar Gantasala
email: sudhakar.gantasala@ltu.se

Jean-Claude Luneno

Jan-Olov Aidanpää

Michel Cervantes

Aeroelastic simulations of wind turbine using 13 DOF rigid beam model

Keywords:

The vibration behavior of wind turbine substructures is mainly dominated by their first few vibration modes because wind turbines operate at low rotational speeds. In this study, 13 degrees of freedom (DOF) model of a wind turbine is derived considering fundamental vibration modes of the tower and blades which are modelled as rigid beams with torsional springs attached at their root. Linear equations of motion (EOM) governing the structural behavior of wind turbines are derived by assuming small amplitude vibrations. This model is used to study the coupling between the structural and aerodynamic behavior of NREL 5 MW model wind turbine. Aeroelastic natural frequencies of the current model are compared with the results obtained from the finite element model of this wind turbine. Quasi-steady aerodynamic loads are calculated considering wind velocity changes due to height and tower shadow effects. In this study, vibration responses are simulated at various wind velocities. The derived 13 DOF simplified model of the wind turbine enables to simulate the influence of change in parameters and operating conditions on vibration behavior with less computational effort. Besides that, the results of the simplified models can be interpreted with much ease.

INTRODUCTION

Wind turbine structures consist of both rotating and nonrotating substructures manufactured in different materials. In order to perform a dynamic analysis of wind turbines, detailed structural model along with loads predicted from the aerodynamic analysis of the blades are needed. Dynamic changes in the wind and rotational effects force flexible blades to vibrate. Rotation and vibrations of the blades change the effective wind velocity generating aerodynamic loads. Thus the structural and aerodynamic analyses of the wind turbine are coupled to each other, dynamic behavior of the wind turbine structures can be accurately predicted considering this coupling. Finite element methods (FEM) can be used to model structural behavior and computational fluid dynamics (CFD) techniques are used for predicting aerodynamic loads. Coupled structural and aerodynamic analysis using detailed FEM and CFD models is computationally expensive for parametric studies in the design stage. To determine the influence of this coupling on the vibration behavior, simple models for structural and aerodynamic calculations are used in this study.

In literature, rigid and flexible beam models are used to analyze the dynamics of the isolated rotating blades considering pitch, flap, and lead-lag vibrations. Chopra [START_REF] Chopra | Nonlinear dynamic response of wind turbine rotors[END_REF][START_REF] Chopra | Non-linear dynamic response of a wind turbine blade[END_REF] considered rigid wind turbine blades with root hinges and analyzed its linear and nonlinear dynamic behavior with and without external loadings. Kim and Lee [START_REF] Kim | Stability analysis of large-scale wind turbines considering nonsymmetric configuration of rotor blades[END_REF] considered a simple structural model for wind turbine with only 10 degrees of freedom (DOF). They analyzed the effect of asymmetry in the blade stiffnesses on the system stability over a range of rotational speeds. Ramakrishnan and Feeny [START_REF] Ramakrishnan | In-plane nonlinear dynamics of wind turbine blades[END_REF] modeled wind turbine blade as a non-linear flexible beam and reduced its equation of motion to a nonlinear Mathieu equation to study super-harmonic resonances. Wind turbine blades are made of aerofoil cross sections which are twisted and tapered along the blade length. The vibration behavior in bending and torsion are therefore coupled. Kim and Lee [START_REF] Kim | Coupled bending and torsional vibration analysis of flexible wind turbine blades by using assumed modes method[END_REF] derived the equations of motion for coupled flapping, lead-lag and torsional vibrations of pre-twisted wind turbine blades. They calculated modal frequencies of the blade using the assumed modeshape method.

Fundamental vibration modes of all the substructures contribute more to the dynamic response of the wind turbine than the other higher frequency vibration modes because wind turbines operate at low rotational speeds. In this study, 13 DOF model of the wind turbine structure is derived considering fundamental vibration modes of all the substructures. Tower and blades are modeled as rigid beams with equivalent inertias and supported by springs at their roots. This configuration duplicates first bending vibration modes of the substructures through the deformation of the torsional spring (as shown in Fig. 1). Lagrange's equations are used to derive the equations of motion (EOM) for a three blade horizontal axis wind turbine. MATLAB symbolic math toolbox is used to derive the EOM. Turbine components like the hub and nacelle are assumed to be rigid bodies. Tower and blades kinetic energies are calculated considering their distributed mass properties. Small amplitude vibrations are assumed for all the DOF and EOM are linearized by retaining only the first order terms. Structural modeling of the wind turbine stationary and rotating substructures together will yield coupled nonlinear differential equations with periodic coefficients irrespective of the type of the reference coordinate frame. Periodic coefficients of the linear differential equations for symmetric blades can be eliminated by using multi-blade coordinate transformation [START_REF] Kim | Stability analysis of large-scale wind turbines considering nonsymmetric configuration of rotor blades[END_REF][START_REF] Bir | User's guide to MBC3: Multi-blade coordinate transformation code for 3-bladed wind turbine[END_REF] and eigenvalue analysis can be used to calculate natural frequencies of the structure. MBC transformation models collective blades motion in the stationary frame of reference. For symmetric blades, MBC transformation yields constant coefficients in the EOM if the gravity effect of the rotating blades is ignored. The influence of the gravity on natural frequencies can be ignored [START_REF] Kim | Coupled bending and torsional vibration analysis of flexible wind turbine blades by using assumed modes method[END_REF]. Blade element momentum (BEM) theory [START_REF] Martin | Aerodynamic of wind turbines[END_REF] is used for calculating quasi-steady aerodynamic loads acting on the blades. Wind shear is modeled using wind profile power law using mean wind velocity at nacelle height. Tower shadow is modeled using potential theory of flow around a cylinder for calculating wind velocity near the tower [START_REF] Martin | Aerodynamic of wind turbines[END_REF]. NREL 5 MW model wind turbine data [START_REF] Jonkman | Definition of a 5-MW reference wind turbine for offshore system development[END_REF][START_REF] Bir | Modal dynamics of large wind turbines with different support structures[END_REF] is considered as an example for this study. Aeroelastic natural frequencies of this 13 DOF rigid beam model in the operating wind velocity range of 3-25 m/s are compared with the results predicted from the FEM model of this wind turbine [START_REF] Bir | Aeroelastic instabilities of large offshore and onshore wind turbines[END_REF]. Vibration responses of the 13 DOF model are simulated at various wind velocities using the operating parameters like rotating speeds, pitch angles of the NREL 5 MW model wind turbine [START_REF] Jonkman | Definition of a 5-MW reference wind turbine for offshore system development[END_REF] and considering aerodynamic coupling with structural behavior. In this study, variations in the wind velocity due to wind shear and tower shadow are only considered. Thus the spectral vibration responses of the 13 DOF model consist of frequencies 1P, 3P (P refers to the rotational frequency) and their harmonics.

STRUCTURAL MODEL

Geometric model of the wind turbine with the considered degrees of freedom of the substructures is shown in Fig. 2(a). Reference coordinate system OXYZ is fixed at the base of the tower which is used to derive the EOM. The tower can oscillate about X and Y axes at the origin O with DOF's α and β respectively, which are known as tower fore-aft and side-to-side vibrations. Hub and nacelle are considered as two rigid masses with their center of gravity (C.G.) locations C and B at distances l h and l n respectively from the yaw axis (point A). These two bodies are fixed on top of the tower and can rotate relative to the tower about the yaw axis with a DOF γ . The nacelle is a stationary part which encloses the generator and other systems. The generator has both stationary and rotating parts. The generator rotor is connected to the hub through a gearbox. In this model, inertias of the generator rotor and the hub are lumped to a single inertia and attached to an equivalent drive shaft. Torsional vibration of the equivalent drive shaft is denoted asψ . Blade vibrates with inplane (lead-lag), out-of-plane (flap) and torsional DOF i i i , , ε ϕ θ (i refers to the blade number) respectively. These motions are defined about a local coordinate system

x bi y bi z bi fixed at the root of the blade as shown in Fig. 2(a).

Position vectors are defined initially and then energy equations for all the substructures are obtained to derive the EOM using Lagrange's equations. Kinetic energy is denoted by KE and potential energy is denoted by PE.

Tower:

As the tower cross section changes along the length, kinetic energy of a small element of length dy located at a distance y from the origin O along the axis Y is defined initially and it is integrated over the full length of the tower to calculate total kinetic energy of the tower. Position vector of a small element of length dy is defined below.

( ) ( )           = 0 0 y T T q Z X tower β α  (1)
where T X (α) and T Z (β) are the transformation matrices about X and Z axis (given in Appendix).

Kinetic, potential energies of the tower and generalized dissipative forces due to damping are given below.

( ) β α ρ β α ρ β β α α β α     C Q , C Q dy q g ) A ( K K PE dy dt q d ) A ( KE L tower tower tower L tower tower tower Y = = + + = = ∫ ∫ 0 2 2 0 2 2 1 2 1 2 1 (2)
where KE tower and PE tower are the kinetic and potential energies of the tower; L is the length of the tower; Q α and Q β are the generalized forces; ρA tower is mass per unit length of the tower;

( ) Y tower q 
is the Y coordinate of the position vector of the tower; g is acceleration due to gravity; K α , K β and C α , C β are stiffness, damping coefficients of the springs and dampers restricting tower vibrations α, β.

Hub and Nacelle:

A hub is a joint that connects turbine blades to the rotating shaft and transfers energy to the generator through the gearbox. In this model, hub rotating on the low speed shaft and generator rotating on the high speed shaft are modeled as a single inertia. This equivalent inertia mass rotates with an angular velocity of ω and undergoes torsional oscillations on an equivalent drive shaft whose torsional stiffness and damping coefficients are known. As this equivalent inertia mass is located on top of the tower, gyroscopic moments are generated in the system due to the change in orientation of its angular velocity vector due to the tower vibrations and yaw motion. Rotational kinetic energy of the equivalent inertia mass accounts for the gyroscopic moments which are calculated using inertia tensor defined about its C.G. and angular velocity vector defined in the local coordinate system R 1 :x 1 y 1 z 1 as shown in Fig. 2(b). The order of rotations followed to get the position vector is ( )

ψ ω γ β α + t , , ,
respectively and the same rotation order is used to transform the velocity vectors

γ β α    , ,
into the R 1 :x 1 y 1 z 1 coordinate system. Position, angular velocity vectors, energy expressions of the rigid hub and nacelle masses and generalized dissipative forces due to damping are given in Eq. ( 3) & (4). L 1 is the distance of hub rotational axis above the ground along the tower centerline;

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ψ γ ψ γ α β γ ψ ω β γ ψ ω γ ψ ω ψ ω γ β α ψ ψ γ γ ψ γ             C Q , C Q q g m K K PE IP ID ID dt q d m KE T T t T T t T t T l L T T T q Y h h hub h h h T h h h hub T Z T Y T Z T Y T Z T Z h h Y Z X h = = + + = Ω           Ω + =           + +           + +           + +           + = Ω           = 2 2 1 2 1 2 1 2 1 0 0 0 0 0 0 2 1 2 1 0 0 0 0 0 0 0 0 0 (3) ( ) ( ) ( ) ( ) ( ) ( ) ( ) Y n n nacelle n n n n T n n n nacelle T Z T Y T Y n n Y Z X n q g m PE IP ID ID dt q d m KE T T T l L T T T q          = Ω           Ω + =           +           +           = Ω           - = 0 0 0 0 0 0 2 1 2 1 0 0 0 0 0 0 0 2 1 α β γ β γ γ γ β α (4) where h q  , n q  and h Ω  , n Ω 
( ) ( ) Y Y n h q , q   are the Y coordinates of

Blades:

Blade in-plane (lead-lag) vibrations defined about z bi axis of the rotating coordinate system placed at the root of the blade are denoted by i θ , out-of-plane (flap) vibrations defined about y bi axis are denoted by i ϕ and torsional vibrations defined about x bi axis are denoted by i ε (refer Fig. 2(a)),

where i = 1,2,3 denote the blade number. The kinetic energy of the blade is obtained by integrating the kinetic energy of a small element of length dr located at a radial distance of r from the blade root, over the full length of the blade (l b ).

Position vector of the small element dr is given below.

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )                       + +           + +           = z y i X i Y i Z Z h Z h Y Z X b cg cg r T T T t T r t T l L T T T q i ε ϕ θ ψ ω ψ ω γ β α 0 0 0 1  (5)
where

( ) ( ) ( ) i X i Y i Z T , T , T ε ϕ θ
are the transformation matrices given in the Appendix; (r,cg y ,cg z ) are the coordinates of the C.G. of the blade section in local coordinate system x bi y bi z bi . Energy expressions for the i th blade are given below.

( )

i b i b i b l b blade i i i b l b blade b i i i i i i b Y i i i i i b i i C Q C Q C Q dr q g ) A ( K K K PE dr dt q d ) A ( KE ε ϕ θ ρ ε ϕ θ ρ ε ϕ θ ε ϕ θ ε ϕ θ      = = = + + + = = ∫ ∫ 0 2 2 2 0 2 2 1 2 1 2 1 2 1 (6)
where

i i i i i b b b b b Q , Q , Q , PE , KE ε ϕ θ
are the kinetic, potential energies and generalized forces of the i th blade; ρA blade is the mass per unit length of the blade;

( ) Y i b q
 is the Y coordinate of the position vector of the i th blade;

i i i K , K , K ε ϕ θ and i i i C , C , C ε ϕ θ
are stiffness, damping coefficients of the springs and dampers restricting i th blade in-plane, out-of-plane, torsional vibrations Equations of motion are derived from the total energy expressions of the system using Lagrange's equations given in Eq. ( 8) for all the generalized coordinates.

( ) ( ) ( )

0 = + ∂ ∂ + ∂ ∂ -         ∂ ∂ j q j j j Q q PE q KE q KE dt d   (8)
where j q  generalized coordinates

i i i , , , , , , ε ϕ θ ψ γ β α (i =1,2,3)
Equations of motion obtained from the Eq. ( 8) are coupled nonlinear differential equations in the generalized coordinates. In the current work, small amplitude vibrations are assumed and linear EOM are obtained after retaining first order terms. Linear EOM in the matrix form is given in Eq. ( 9) and whose matrices contain periodic terms.

{ } { } { } { } ) ( ) ( ) ( ) ( ) ( ) ( ) ( t f t q t K t q t C t q t M = + + ω ω ω    (9)
where M(ωt), C(ωt), K(ωt) are the system matrices, ,,,,,,,,,,(t) q

{ } [ ] T , ,
3 2 1 3 2 1 3 2 1 ε ε ε ϕ ϕ ϕ θ θ θ ψ γ β α = and
{f(t)} is the force vector Rotating blade's DOF are transformed into multi blade coordinates [START_REF] Bir | User's guide to MBC3: Multi-blade coordinate transformation code for 3-bladed wind turbine[END_REF] using the transformation matrix given in Eq. ( 10).

+ = + = = t , t , t
and I 4 is a 4x4 identity matrix which retains first four nonrotating DOF as it is in the final transformed coordinate vector. The relationship between { } ) (t q and multi-blade coordinate's vector { }

) (t q B
is given in Eq. ( 11).

{ } { }

) t ( q T ) t ( q B = (11) where { } [ ] T s c s c s c B , , , , , , , , , , , , (t) q
ε ε ε ϕ ϕ ϕ θ θ θ ψ γ β α 0 0 0 = Difference between { } ) (t q and { } ) (t q B
is only in the blades in-plane, out-of-plane and torsional vibrations. Subscripts 0, c, s refer to collective, progressive and regressive modes of the blade assembly which can excite tower vibration modes. In collective mode, the vibrations of the blades are in phase with each other, whereas in the progressive and regressive modes, phase difference exists between blade vibrations as defined in the matrix T B (ωt).

Detailed explanation of these coordinates and MBC transformations are given in [START_REF] Bir | User's guide to MBC3: Multi-blade coordinate transformation code for 3-bladed wind turbine[END_REF]. The EOM obtained after MBC transformation in the matrix form are given in Eq. ( 12).

( ) ( ) { } ( ) ( ) { } ( ) ( ) { } ( ) { } t f t q t K t q t C t q t M B B B B B B B = + + ω ω ω    (12) where ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) t f T t f , T t M T T t C T T t K T t K , T t M T T t C T t C , T t M T t M T B T T T B T T B T B = + + = + = =     ω ω ω ω ω ω ω ω ω 2
Gravity effect changes stiffness over the azimuthal rotation of the blade [START_REF] Kim | Coupled bending and torsional vibration analysis of flexible wind turbine blades by using assumed modes method[END_REF] which makes K B matrix still periodic in time even after MBC transformation. Neglecting gravity term makes the system matrices time invariant for symmetric blades, and therefore eigenvalue analysis can be used to calculate natural frequencies. However, for asymmetric blades, time invariant system matrices cannot be obtained using MBC transformation even after neglecting gravity terms.

Data required to simulate the derived 13 DOF model are computed from the structural details of NREL 5 MW model wind turbine. The inertia properties of all the subsystems can be calculated from the data given in [START_REF] Jonkman | Definition of a 5-MW reference wind turbine for offshore system development[END_REF][START_REF] Bir | Modal dynamics of large wind turbines with different support structures[END_REF]. Fundamental vibration modes of all the subsystems obtained from the FEM analysis of this wind turbine structure are also reported in these references. Using the natural frequencies and inertia properties reported in [START_REF] Jonkman | Definition of a 5-MW reference wind turbine for offshore system development[END_REF][START_REF] Bir | Modal dynamics of large wind turbines with different support structures[END_REF], unknown equivalent spring stiffness coefficients for vibration DOF α, β, i i i , , ε ϕ θ are calculated.

AERODYNAMIC MODEL

Beam element momentum (BEM) theory [START_REF] Martin | Aerodynamic of wind turbines[END_REF] is used to calculate aerodynamic loads considering wind shear and tower shadow effects. Wind shear is modeled using wind profile power law using mean wind velocity at nacelle height and a parameter of 0.2 for the amount of shear. Tower shadow is modeled using potential theory of flow around a cylinder [START_REF] Martin | Aerodynamic of wind turbines[END_REF]. The tower diameter variation along its height is also taken into account in the calculation of the wind velocity near the tower. Prandtl's tip loss factor and Glauert corrections are considered in the BEM theory used in this study. Aerodynamic loads are calculated at the aerodynamic center of the aerofoil sections of the blade. Velocity triangle without and with considering blade vibrations are shown in Fig. 3. Blade vibration velocities change the relative velocity of wind entering the blade section as shown in Fig. 3(b). As the structural model of the blade is built with respect to pitch axis, due to the offset between aerodynamic center and pitch axis, both the inflow angle and the angle of attack (AOA) depend on torsional vibrations. Expressions for the inflow angle, angle of attack considering blade vibrations are given in Eq. ( 15) & ( 16).

( ) ( )         + + + - = Θ - 2 1 1 1 1 V ' a r V a V tan w if ω (15)
where,

( ) ( ) i i t p AC i i i t p AC i sin a r V cos a r V ε ε θ ε ε ϕ     + Θ - = + Θ + - = + + 2 1 ( ) i t p if aoa ε + Θ - Θ = Θ + (16)
where, aoa if ,Θ Θ are the inflow angle and angle of attack considering blade vibrations;

i i i , , ε ϕ θ   
are the velocities of the i th blade in-plane, out-of-plane and torsional vibrations; i ε is the i th blade torsional vibration; a AC is the distance between the aerodynamic center (A.C.) and pitch axis (P.A.) (shown in Fig. 3) of the blade section. Change in the AOA due to blade vibrations can be approximated according to Eq. ( 17).

( ) ( )

i rel w rel aoa aoa V V a V V V ' a r ε ω - - - + ≈ Θ - Θ = ∆Θ
) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )∆Θ Θ + Θ ≈ Θ ∆Θ Θ + Θ ≈ Θ ∆Θ Θ + Θ ≈ Θ 0 0 0 0 0 0 aoa M aoa M aoa M aoa D aoa D aoa D aoa L aoa L aoa L ' C C C ' C C C ' C C C (18)
) ( ) ( ) ( ) { } i rel w aoa L aoa L rel V V a V V ' a r ' C c C cV L ε ω ρ ρ 2 2 1 2 0 0 0 0 1 1 2 1 2 1 - - - + Θ + Θ ≈ (19) ( ) ( ) ( ) ( ) { } i rel w aoa D aoa D rel V V a V V ' a r ' C c C cV D ε ω ρ ρ 2 2 1 2 0 0 0 0 1 1 2 1 2 1 - - - + Θ + Θ ≈ (20) ( ) ( ) ( ) ( ) { } i rel w aoa M aoa M rel V V a V V ' a r ' C c C V c M ε ω ρ ρ 2 2 1 2 2 2 0 0 0 0 1 1 2 1 2 1 - - - + Θ + Θ ≈ (21)
The aerodynamic loads in Eq. ( 19)-( 21) are introduced in the structural model given in Eq. ( 9) and the EOM considering structural and aerodynamic coupling is given in Eq. ( 22). Aeroelastic natural frequencies can be calculated by transforming { }

) (t q in Eq. (22) to multi-blade coordinates { } ) (t q B
and using the eigenvalue analysis of the system.

{ } [ ] { } [ ] { } { } { } ) t ( f ) t ( f ) t ( q K ) t ( K ) t ( q C ) t ( C ) t ( q ) t ( M Aero Aero Aero + = + + + + ω ω ω    (22)

RESULTS AND DISCUSSION

Data of the model wind turbine [START_REF] Jonkman | Definition of a 5-MW reference wind turbine for offshore system development[END_REF] are inserted into the system of Eq. ( 12) and eigenvalues are obtained at zero speed and by ignoring gravity terms. Natural frequencies of the present 13 DOF model are compared with the values predicted using FAST (aeroelastic simulation tool) [START_REF] Jonkman | Definition of a 5-MW reference wind turbine for offshore system development[END_REF] in Table 1. As FAST doesn't model blade torsional vibrations its natural frequency in the current rigid beam model is compared with the frequency value predicted in [START_REF] Kim | Coupled bending and torsional vibration analysis of flexible wind turbine blades by using assumed modes method[END_REF] where blade is modeled using FEM. Natural frequencies calculated using these two models (FEM vs. 13 DOF model) at zero speed match well and thus this 13 DOF model can be considered as a good representation of dynamic behavior of the wind turbine structure in the low frequency vibration modes.

NREL 5 MW wind turbine operates over a wind velocity range of 3-25 m/s and its rotational speed increases from 7.35 rpm at a wind velocity of 3 m/s to 12.1 rpm at a wind velocity of 11.4 m/s and operates at this rotational speed till 25 m/s wind velocity by changing pitch angle of the blades [START_REF] Jonkman | Definition of a 5-MW reference wind turbine for offshore system development[END_REF]. Natural frequencies of the 13 DOF model considering only the rotational speed variation with wind velocities between 3-25 m/s are shown in Fig. 4. Natural frequencies of the blade DOF change till 11.4 m/s wind velocity due to increase in rotational speed and thereafter remains constant as there is no change in the rotational speed. Aerodynamic loads calculated for mean wind velocities in the range between 3-25 m/s at the nacelle are used to determine aeroelastic natural frequencies of this model which are plotted in Fig. 5. Differences in the natural frequencies in Fig. 4 & 5 are only due to the aerodynamic coupling considered in the latter case.

Aeroelastic natural frequencies calculated using the FEM model of this wind turbine are reported in Ref. [START_REF] Bir | Aeroelastic instabilities of large offshore and onshore wind turbines[END_REF] which are shown in Fig. 6. All vibration modes except tower bending modes change with an increase in wind velocity. Lag modes in Fig. 5 & 6 follow the same trend with increasing wind velocities. Flap modes decrease with an increase in wind velocity in Fig. 6, whereas in Fig. 5, only the collective and regressive flap modes are decreasing with

Figure 1 .

 1 Figure 1. Wind turbine: (a) geometric model, (b) flexible beam model, (c) rigid beam model

Figure 2 .

 2 Figure 2. (a) Geometric model of the wind turbine structure with the considered DOF in the model, (b) Hub and nacelle local coordinate systems in which inertia of masses are defined

  the position vectors of the hub and nacelle C.G.; K γ and C γ are stiffness and damping coefficients of the springs and dampers of the yaw actuation system which control relative position of the rigid hub and nacelle mass assembly with respect to tower; K ψ and C ψ are stiffness and damping coefficients of the equivalent drive shaft.

.

  Total kinetic and potential energies of the system are obtained by summing up individual contributions from all the substructures.

  azimuthal positions of the blades with respect to a fixed reference. The blades are placed symmetrically around the circumference of the hub.

Figure 3 .

 3 Figure 3. Velocity triangle at the blade section (a) without and (b) with considering blade vibrations Expressions for the inflow angle, angle of attack without considering blade vibrations are given in Eq. (13) & (14).

  Θare the inflow angle and angle of attack, ignoring the blade vibrations;t p + Θis the sum of blade pitch and section twist angles; V w and ω are the wind velocity and rotational frequency of the blade; r is the radial distance of the blade section from the hub center; a and a' are the axial and tangential induction factors.

  and pitching moment coefficients at aoa Θ can be expanded using Taylor series expansion about the 0 aoa Θ as given in Eq. (18).

(

  where, C L , C D , C M are the lift, drag and moment coefficients of the aerofoil; C L ', C D ', C M ' are the slopes of the lift, drag and moment coefficient curves. Aerodynamic loads considering these coefficients are expressed in Eq. (19)-(21) where, if written in matrix form the coefficients of the variables into matrices C Aero and K Aero which are known as aerodynamic damping and stiffness matrices.

(

Figure 4 .

 4 Figure 4. Natural frequencies variation considering rotational speed changes with wind velocities using 13 DOF model increase in wind velocity. Aeroelastic flap modes will be accurately predicted by flexible beam models which captures aerodynamic coupling accurately. Differences in the aeroelastic natural frequencies reported in Fig. 5 and Fig. 6 arise due to simplification of flexible structure with rigid beam model and can be attributed to two reasons: centrifugal stiffening and aerodynamic loads. The rigid beam model

Figure 5 .

 5 Figure 5. Natural frequencies variation considering aerodynamic coupling and rotational speed changes at various wind velocities using 13 DOF model

Figure 6 .

 6 Figure 6. Natural frequencies variation considering aerodynamic coupling with FEM model of the NREL 5 MW model wind turbine [10]

Figure 7 .

 7 Figure 7. (a) & (b) Loads exciting blade in-plane vibrations, (c) & (d) loads exciting blade out-of-plane vibrations

Figure 8 .

 8 Figure 8. (a) Tower top fore-aft vibrations, (b) Tower top sideto-side vibrations, (c) Blade tip in-plane vibrations, (d) Blade tip out-of-plane vibrations

  are the position, angular velocity vectors of the hub and nacelle masses m h , m n respectively; T Y (γ) and T Z (ωt+ψ) are the transformation matrices about Y and Z axis (given in Appendix); KE hub , KE nacelle , PE hub , PE nacelle

are the kinetic and potential energies of the hub and nacelle masses; ID h , ID n are transverse mass moment of inertias of the hub and nacelle masses about the local coordinate system defined at their mass centers (shown in Fig.

2(b)

) parallel to OXYZ coordinate system; IP 1 is the sum of hub and generator polar moment of inertias about the low speed shaft

LSS

Gen h IP , IP

; IP n is the nacelle polar moment of inertia;

Table 1 .

 1 Comparison of natural frequencies

		Natural frequency
	Mode	FEM	(Hz)	13 DOF Model
	Tower fore-aft (α)	0.3240 *	0.3237
	Tower side-to-side (β)	0.3120 *	0.3291
	Yaw motion (γ)	-			6.8346
	Drivetrain torsion (ψ)	0.6205 *	0.6099
	Blade regressive lead-lag (θ s )	1.0793 *	1.1170
	Blade progressive lead-lag (θ c )	1.0898 *	1.1399
	Blade collective lead-lag (θ 0 )	-			3.7978
	Blade regressive flap (φ s )	0.6664 *	0.6399
	Blade progressive flap (φ c )	0.6675 *	0.6542
	Blade collective flap (φ 0 )	0.6993 *	0.6879
	Blade regressive torsion (ε s )	-			5.5705
	Blade progressive torsion (ε c )	-			5.5738
	Blade collective torsion (ε 0 )	5.6	~		5.5708

* From Ref.[8],

~ From Ref.

[START_REF] Kim | Coupled bending and torsional vibration analysis of flexible wind turbine blades by using assumed modes method[END_REF]


domain: spi sde

text:
Stefan Gerlach
email: stefan.gerlach@tu-berlin.de

Rita Ugarelli

Paul Uwe Thamsen

Case Study on the Functional Performance of a large Wastewater Pumping Station

Keywords:

In this paper a case study on a large wastewater pumping station is conducted. In a Fault Tree Analysis the critical assets within a wastewater pumping station are identified and the root causes for critical failures are discussed. The relevant Performance Indicators (PIs) for centrifugal pumps are presented and the significance for a PI including wastewater properties and performance under wastewater conditions is discussed. The performance of the critical wastewater pumps in the observed wastewater pumping station is analyzed. An approach for risk assessment for clogging of wastewater pumps is described and a strategy to counter clogging problems is addressed.

INTRODUCTION

Climate change, demographics and water-sensitive behaviour put a great strain on the urban wastewater system by changing amounts, composition and behaviour of wastewater. Especially, combined sewer systems face big changes of loads between dry-weather runoff as the minimum and stormwater as maximum load. Additionally operating conditions, the challenging fluid "Wastewater" and the regulatory background [START_REF] Mudgal | Work on Preparatory studies for implementing measures of the Ecodesign Directive 2009/125/EC, ENER Lot 28-Pumps for private and public wastewater and for fluids with high solids content -Task 1: Definition[END_REF], [START_REF] Hoechel | Analysing sewage handling problems[END_REF], [START_REF] Kroeber | Experiences with speedcontrolled operation of sewage pumps at the Berliner Wasserbetriebe[END_REF] create an operational frame for the main wastewater pumping station that results in a high risk of failures. Especially fibres and textile structures within the wastewater lead to frequent failures in wastewater pumps due to clogging.

In a structured approach based on the frame of Integrated Asset Management (IAM) the case study is oriented along the "Five Core Management Questions" of IAM. The necessary steps are addressed below:

CURRENTS STATE THE WASTEWATER PUMPING STATION

The pumping station analysed in this paper is a large pumping station in an urban area. With a connected catchment area of more than 250,000 inhabitants and about 750km of sewer system this main pumping station delivers about 45,000m³/day of dry-weather run-off to the three connected wastewater treatment plants. The power consumption of the 8 installed pumps ranges from 40 kW up to 455 kW for electrical power and up to 711 kW for diesel-driven peak-load pumps. These machines are used to pump from 100 l/s during nights up to 1450 l/s during rain events. The pump setup is redundant.

Figure 1. Wastewater Pumping station [START_REF] Gerlach | Personally taken photographs[END_REF] The most common way wastewater pumping stations are operated in is on-off-operation. The volume of the suction chamber is used as a buffer. Once a specific water level is reached fixed speed pump operation is started. Once the level falls below a defined minimum pumping is stopped. For high levels parallel operation of multiple pumps can be implemented. The use of Variable Speed Drives (VSD) allows adjusting the pump operation according to the inflow. In the observed pumping station a so-called level-controlled operation is used. Figure 2 shows the schematic layout of a levelcontrolled pump operation. Each box shows the flow rates the specific pump delivers according to the respective level in the pump sump. The left boundary of each box shows the flow rate at minimum pump speed and the right boundary the flow rate at maximum pump speed for the pump or pump combination. Below a specific minimum level in the suction chamber pump operation is stopped. Once the level increases above this minimum the smallest pump is started with minimum speed. If the level in the suction chamber increases further the pumps speed is also increased. If maximum pump speed is reached, the small pump is stopped and the midsized pump starts operation with minimum speed. As long as the level still increases the pump speed is also increased accordingly. The next step is parallel operation of mid-sized pumps. For peak loads and emergencies parallel operation of mid-size and peak-load pumps is implemented.

Figure 2. Scheme for level-controlled operation based on actual operating schemes provided by Berliner Wasserbetriebe

In many cases urban wastewater systems are set up hierarchically with a number of wastewater pumping stations spread over the catchment area that transport the wastewater to one or more main pumping stations. These main pumping stations pump the wastewater through pressure lines to the wastewater treatment facilities. The first step to identify the crucial pumping station is a criticality analysis of the wastewater network. The observed catchment area includes only one main pumping station that is fed by a gravity sewer system. Naturally this main pumping station is identified as most critical pumping station in the system.

The Fault Tree Analysis (FTA) for the wastewater pumping station results in the main events leading to a critical failure of the whole pumping station. Figure 3 shows that out of 8 pumps only 4 pumps contribute to a failure of the whole pumping station under certain conditions. In this context the FTA identifies Pump 1 and Pump 2 as the most critical pumps in the pumping station.

PERFORMANCE OF THE WASTEWATER PUMPING STATION

PERFORMANCE INDICATORS

The evaluation of performance of pumps in pumping systems is generally based on hydraulic characteristics and energy considerations such as energy-efficiency or specific energy E s . [START_REF] Mudgal | Work on Preparatory studies for implementing measures of the Ecodesign Directive 2009/125/EC, ENER Lot 28-Pumps for private and public wastewater and for fluids with high solids content -Task 1: Definition[END_REF] The evaluation of the performance of wastewater pumps however is more difficult. At this moment, the European Commission is implementing energy-efficiency regulation for wastewater pumps in the frame of the European Ecodesign Directive [11]. In a preliminary study it was shown that no standard exists to link the pump performance to anything other than energy-efficiency [START_REF] Mudgal | Work on Preparatory studies for implementing measures of the Ecodesign Directive 2009/125/EC, ENER Lot 28-Pumps for private and public wastewater and for fluids with high solids content -Task 1: Definition[END_REF]. However, the performance of wastewater pumps is evaluated through the following parameters:

• Hydraulic performance • Efficiency • Reliability
Pumps are designed for a specific hydraulic task, i.e. a specific flow rate at a specific head. Any deviation between this design point and the actual operating point results in loss of energy-efficiency.

The operating point itself is predefined by system conditions and characteristic pump curve. Any change of system conditions, i.e. system pressure, or flow conditions within the pump, i.e. clogging, result in a shift of the operating point.

The efficiency of pumps is evaluated in accordance to ISO 9906 [START_REF] Iso | Rotodynamic pumps --Hydraulic performance acceptance tests --Grades 1[END_REF]. This standard defines the testing setup, the testing conditions and the measurement setup to obtain data on hydraulic performance and energy demand. It describes solely the hydraulic performance for clear water operation. A market-study [START_REF] Hoechel | Market study on 4000 wastewater pumps[END_REF] for the most-common impeller types already shows that different impeller types for wastewater pumps result in different expected energyefficiencies for clear water operation. A direct comparison of existing wastewater pumps can only be conducted under the consideration of the real operating conditions, including equal compositions of wastewater. On-going research in cooperation with the pump manufacturers associations EUROPUMP is showing the different performance of different impeller technologies under real wastewater conditions. It is already known that vortex impellers operate with the low energy-efficiencies but also with high reliability for difficult wastewater compositions especially for high textile or fibre load. Other impeller designs however show higher efficiencies but operate with lower reliability. To define a scientific foundation for the comparison of reliability for different impeller types is focus of present research at the TU Berlin.

To increase the significance of wastewater pump testing a new Performance Indicator has to be defined. Within a research project with the pump manufacturers associations the Function-Efficiency-Index FEI is proposed [START_REF] Pöhler | Linking the efficiency to functional performance by a pump test standard for wastewater pumps[END_REF]:

𝐹𝐸𝐼 = (𝛼 • 𝐷 𝐶 + 𝛽 • 𝐷 𝐹 + 𝛾 • 𝐷 𝐿𝑇𝐹 ) • 𝜌 • 𝑔 • 𝐻 • 𝑄 𝑃 (1)
where ρ is the water density, g is gravitational constant, H is lifting height (pressure increase) and Q is the flow delivered by the pump, altogether representing the pump output P Q .

The required power is P. The ration of P Q and P is the pump efficiency for clear water operation. The energy-efficiency is Finally, manufactures have developed many different impeller types and shapes with different behaviors in wastewater operation (Figure 5).

These parameters evidently show that the performance of wastewater pumps is not only defined by the pump type, but additionally by the conditions they are operated in. To this date, there are no existing standards or testing procedures that consider their impact on the actual pump performance.

Actual performance of critical pumps at the observed pumping station 2.2.1 Reliability

The task of the observed pumping station is the safe transport of wastewater from the sewer system to the connected wastewater treatment plants. Apart from energyconsiderations, the hydraulic performance and the reliability of the pumping station are of main concern. The observed pumping station shows a clear distribution of failure events for the last years. While most of the pumps show a low number of clogging events, Pump 1 and Pump 2 show equally high numbers of failures. In the observed period of 12 month, both pumps clogged 78 times, resulting in a high demand for manual cleaning(Figure 6). These number exceeded mechanical or other failures by far. In most cases, the clogging material could be removed through a service opening at the pump. In some cases the pumps had to be disassembled resulting in costs exceeding 1000€ per assignment. The common clogging mechanism is ragging in front of the impeller. Actual observation shows slow accumulations of fibrous material in front of the impellers that subsequently lead to failure of the pumps. It shows that the most critical pumps in the pumping station are also the pumps with the highest rate of failures due to clogging. While for low pump speeds most operating points are not affected by clogging the picture changes for higher pumps speeds. Especially for the highest speeds it is shown that the actual pump performance is affected by clogging to high extend and the observed pump never actually works under design conditions when operated with full pump speed.

Hydraulic performance and Energyefficiency

[Figure 9] This behaviour is supported by level-controlled operation of the pumping station. Once the level in the suction chamber reaches the starting level for Pump 1 or Pump 2, one of the pumps is started at the lowest pump speed of 580rpm. As long as the level in the suction chamber is increasing, the pump speed is gradually increased. In case of beginning accumulation of fibres in the impeller, the actual delivered flow rate decreases resulting in increasing water levels in the suction chamber thus increasing the pump speed. While the accumulation of fibres grows pump speed increases while the actual operating point shift further away from the theoretical Performance Curve. Figure 10 shows that the actual hydraulic capacity of the pumping station is limited. Increased accumulation of fibres affect flow rates to an extent where increasing pump speeds cannot deliver higher flow rates, but flow rates actually decrease. This behaviour is shown for all pump speeds using dimensionless coefficients in Figure 11.

𝜓 = 2 • 𝑔 • 𝐻 𝜋 2 • 𝑛 2 • 𝐷 2 (2)
where g the gravitational acceleration in m/s², H the head in m, n the rotational speed in 1/s and D the impeller diameter in m. The flow coefficient 𝜑 is described by:

𝜑 = 4 • 𝑄 𝜋² • 𝑛 • 𝐷³ (3)
where Q the flow rate in m³/s, n the rotational speed in 1/s and D the impeller diameter in m.

It is shown that the actual pump operation is strongly influenced by clogging for all pump speeds.

The consumed electrical power for the pumps is only measured at the electric feed of the Variable Speed Drives. Since the efficiency of the motors and the variable speed drives are not constant but functions of different operational parameters such as speed, current and torque, a precise determination of the pumps power consumption is not performed at the pumping station. A calculation of energyefficiency for the pumps is not done. The decreased hydraulic power transferred by the pumps due to clogging results in an overall drop of energy-efficiency.

Risks and costs

The risks analysis for the pumping station is done according to an adapted scheme provided by Eiswirth et al. [START_REF] Eiswirth | New methods for water pipeline assessment[END_REF] for the condition assessment for pipes in 2001. Figure 12 was adapted to fit the conditions of the observed pumping station. The risk of failure is a combination of the probability of failures and the consequences connected to the occurring failures.

Figure 12. Risk Assessment for wastewater pumping stations

The probability of pump failure is determined by failures connected to the service life of the pump and external failures due to the difficult fluid to be pumped. Failures that are connected to the service life are countered by detailed maintenance of the pump units. Condition Monitoring is used for predictive maintenance. These procedures are defined by the operator and are not connected to clogging. Failures through externalities for wastewater pumps are hard to predict. The main parameter leading to failures is the wastewater itself in connection to the used pump. At this moment, there are no known deterministic or statistic models that link wastewater composition and pump type to failure probabilities. On the one hand, it is impossible to predict the wastewater composition for any given time; on the other hand, there are no existing performance indicators for wastewater pumps that describe the ability of handling textiles or solids during operation. The proposed FEI will support the decision process for choosing the right wastewater pump for a specific task. In-depth knowledge of wastewater compositions for the planned application are to be provided by the operator. At this time research on mechanical properties of wastewater components is only starting. Anyhow, it is already seen that certain parameters have a negative effect on wastewater pump performance: • Single events: Wastewater transports all sizes of solids that can enter the system at any place upstream the pumping station. This ranges from small debris, to dead animals and bricks from canals, up to car tires or truck tarpaulins. Failure events connected to large fibre accumulations or solids are impossible to predict.

•

Ragging: Ragging is a gradual accumulation of textiles, tissues and fibres in the wastewater pump. This process can be supported by starts and stops of the pump. Gradual ragging of impellers can be the objective of clogging detection and countermeasures at the pump.

•

Wastewater composition and impeller type: At the moment, pump manufacturers describe the capabilities of their pumps for handling different wastewater "classes" very vaguely. Wastewater classes are only described by their origin, e.g. municipal, industrial, and not by their actual load with fibres or textiles. The capability of handling those classes is described by good, neutral or bad.

•

Weather influence: Rain events lead to first-flush effects the lead to increased loads of solids and fibres in the wastewater. The missing link between wastewater composition and actual pump performance does not allow any calculation.

•

Critical operating points and VSD operation: Certain operating conditions aggravate the occurrence of clogging in wastewater pumps. It was observed that low pump speeds and part-load operation increase the risk of clogging [START_REF] Kroeber | Experiences with speedcontrolled operation of sewage pumps at the Berliner Wasserbetriebe[END_REF]. Low pump speeds with low flow rates result in sedimentation of material in the suction chamber. Sediments are remobilised when flow rates increase and can result in higher failure rates.

The consequences of failure at the observed pumping station are divided into external cost driven by service disruption, community disruption and environmental costs and cost directly connected to the clearance of failures. The redundant setup-up of the observed pumping station reduces the risk for external disruption and damage to a minimum. Residents of adjacent buildings however file complaints about vibrations that are connected to operation of clogged pumps. More detailed information was not available to the author.

The costs that are connected to the clearance of failures at the observed pumping station were calculated to an average 310,00€ per assignment in the observed 12 month. Number of failures and connected clearance cost allow a brief calculation of Failure Rates, MTBF, Risk per day and actual Costs per year.

𝐹𝑎𝑖𝑙𝑢𝑟𝑒 𝑅𝑎𝑡𝑒 = 𝐹𝑎𝑖𝑙𝑢𝑟𝑒𝑠 365𝑑 (4)
𝑀𝑇𝐵𝐹 = 1 𝐹𝑎𝑖𝑙𝑢𝑟𝑒 𝑅𝑎𝑡𝑒 (5) 𝑅𝑖𝑠𝑘 𝑝𝑒𝑟 𝑑𝑎𝑦 = 𝐴𝑣𝑒𝑟𝑎𝑔𝑒 𝐶𝑜𝑠𝑡 • 𝐹𝑎𝑖𝑙𝑢𝑟𝑒 𝑅𝑎𝑡𝑒 (6)
𝐶𝑜𝑠𝑡𝑠 𝑝𝑒𝑟 𝑦𝑒𝑎𝑟 = R𝑖𝑠𝑘 𝑝𝑒𝑟 𝑑𝑎𝑦 • 365𝑑

Table 1. Cost calculation for the observed pumping station Table 1 shows that two pump units at the observed pumping station account for two thirds of the overall failure clearance costs during the observed 12 month. In this case, this simplified approach already allows identifying the main objectives for the implementation of risk reduction measures to increase the reliability of the pumping station. The overall significance of the calculated MTBFs is missing. Due to the high influence on the clogging behaviour of the described parameters, additional knowledge to use for future planning cannot be derived. For mechanical failures, low MTBF values create a high pressure of action, resulting in exchange of equipment, change of suppliers or improved maintenance strategies. For wastewater pumping stations the pressure is equally high, while above strategies will most likely not improve the situation.

Plan

High failure risk for the critical pumps at the observed pumping station and a high number of complaints call for action. An exchange of equipment is no option. On the one hand different pumping technology for this specific task is not available on the market and on the other hand a prediction of the actual performance of other impeller technologies for example is at this time not possible. A strategy to reduce the impact of clogging on the actual pump performance should be based on counter-measures. Especially detection of beginning fibre accumulations in front of the impeller connected to automated cleaning procedures show a high potential of reducing complaints due to vibration and reducing cost by decreasing the amount of manual work for cleaning clogged pumps.

In the IMEBA Project several active and passive countermeasures were proposed [START_REF] Thamsen | Diagnose mit aktiver Abwehr von Störungen in Abwasserpumpanlagen[END_REF]. At this point it is planned to implement the Measure Cleaning Sequence at the pumping station. In a two stage approach data is collected to support the set-up of a clogging detection tool, while simultaneously the technical preconditions for implementing VSD-driven "Cleaning Sequence" based on variations of pump speeds are worked out.

Conclusion

A Fault Tree Analysis was used to identify critical pumps for the safe operation of the observed pumping station. The root causes for failures were identified and it was shown that the most critical pumps show the highest rate of failures. An analysis of the hydraulic performance of the wastewater pumping station discovered that clogging not only leads to pump failures but also reduces the attainable flow rates greatly. Additionally clogging results in a loss of energy-efficiency.

A risk assessment scheme is proposed that shows the different influential parameters on pump clogging and the lack of deterministic or statistical models was discussed.

A possible approach to counter pump clogging by implementation of early-detection and an automated cleaning sequence was presented.

ACKNOWLEDGMENTS

This paper is based on work within the KURAS project. The joint research project "KURAS -Concepts for urban

•

  Status: What is the current state of the assets? • Performance: What is the required Level of Service? • Risk: Which of the assets are critical? • Costs: What are my best minimum lifecycle cost CIP and O&M strategies? • Plan: Given above, what are the best long-term funding strategies?

Figure 3 .

 3 Figure 3. Fault Tree Analysis for the observed Wastewater Pumping Station showing clogging as a root cause for critical failure Pump 1 to Pump 4 are the large and peak-load pumps of the pumping station. Two out of these four pumps are needed to deliver the required quantities of wastewater during peak-load or rain operation. While Pump 1 and Pump 2 are driven electrically, Pump 3 and 4 have diesel engines. The allowed operating hours for the diesel engine are limited by clean-air regulations. In this context the FTA identifies Pump 1 and Pump 2 as the most critical pumps in the pumping station.

Figure 4 .

 4 Figure 4. Mean efficiencies of wastewater pumps for different impeller types according to power consumption [4]

Figure 5 .

 5 Figure 5. Different impeller types for wastewater application [9,10]

Figure 6 .

 6 Figure 6. Distribution of clogging events for the installed pumps in a period of 12 month

Figure 7 .

 7 Figure 7. Fibre accumulation in front of the impeller and pump cleaning at the observed pumping station [8]

Figure 8

 8 Figure8shows the theoretical Characteristic Pump Curve for Pump 2 at 780 rpm and the specific System Curve for the connected wastewater pressure line to the connected Wastewater Treatment. Also included are actual operating points for pump speeds between 580rpm and 630rpm in comparison to the scaled pump curves for these speeds. All data was gathered from the operators' datamanagement system over the course of one year. The flow

Figure 8 .

 8 Figure 8. Pump 2: Theoretical Characteristic Pump Curve for Pump 2 and actual operating points for pump speeds 580rpm -630rpm

Figure 9 .

 9 Figure 9. Pump 2: Theoretical Characteristic Pump Curve for Pump 2 and actual operating points for pump speeds 730rpm -780rpm

Figure 10 .

 10 Figure 10. Pump 2: Actual attained flow rates for all pump speeds

Figure 11 :

 11 Figure 11: Pump 2: Pressure coefficient over flow coefficient Pump 2 actual operating points in comparison to the theoretical pump performance curve for all pump speeds

stormwater management and wastewater systems" is funded by the German Federal Ministry of Education and Research and is a research consortium of fourteen cooperation partners. Their contributions to the project's progress and various inputs to the paper are highly appreciated.

NOMENCLATURE


Gaballah

N Kanari

E Allain

I Gaballah

First Version

Keywords:

I. INTRODUCTION

Iron oxides are the major compounds of several valuable metal bearing ores such as ilmenite, chromite, bauxite, etc. They are also the main constituents of the solid wastes and slags generated by metallurgical oxidation operations. The selective extraction of valuable metals from these materials could be achieved using chlorination techniques provided a good understanding of the behavior and reactivity of iron oxides towards different chlorinating agents. For this reason, several investigations concerning the Fe-Cl-O system were carried out, focusing on both the thermodynamic and kinetic aspects.

Details concerning the reaction of iron oxides with different chlorinating agents were described in Reference [START_REF] Kanari | Extraction des métaux de valeur des concentrés de chalcopyrite et de chromite par chloruration[END_REF]. The most usual chlorinating agents were Cl 2 , CO+Cl 2 , COCl 2 , CCl 4 and HCl. The experiments were mostly carried out using thermogravimetric analysis 'TGA'.

The main observations regarding to the chlorination of wüstite were summarized in Table I. It was reported that wüstite starts to react with the chlorinating gas at low temperatures generating several intermediate phases such as FeCl 2 , FeCl 3 , FeOCl, Fe 3 O 4 , and Fe 2 O 3 . On the other hand, the higher oxides of iron (Fe 3 O 4 and Fe 2 O 3 ) were generated even during the chlorination of wüstite in the presence of a reducing atmosphere. However, no studies were devoted to the reaction of wüstite with ferric chloride. HCl 25-850 Formation of FeCl 2 , E a ≈ 111 kJ/mol between 700 °C and 800 °C. [START_REF] Ya | Chlorination of nickel and iron oxides by hydrogen chloride[END_REF] Table II gives some results concerning the chlorination of hematite using various chlorinating agents. The researchers agreed that, whilst using the same chlorinating reagents, the chlorination of hematite starts at higher temperatures than that of wüstite. The final reaction product was found to be FeCl 3 . The values of the apparent activation energy of hematite chlorination varied from about 6 kJ/mol to 190 kJ/mol depending on the chlorinating agent used and on the temperature range explored (Table II).

The present paper focuses on the chlorination and carbochlorination of wüstite and hematite. Both, TG analysis and boat experiments were used. An attempt was made to study the interaction of Fe (1-x) O with FeCl 3 .

Table II 25-900 Oxygen had little from 0 to 40 %. [START_REF] Titi-Manyaka | Thermogravimetric investigation of the chlorination behaviors of some common metals and their oxides[END_REF] HCl 25-900 Formation of FeOCl decomposing from 300 °C to 400 °C. [START_REF] Titi-Manyaka | Thermogravimetric investigation of the chlorination behaviors of some common metals and their oxides[END_REF] HCl 800-1000 E a = 34 kJ/mol, n HCl = 1.

[9]

HCl 25-750 Formation of FeOCl, E a = 5.9 kJ/mol. [START_REF] Fruehan | Rate of chlorination of metals and oxides: Part III. Rate of chlorination of Fe 2 O 3 and NiO in Cl 2 and HCl[END_REF] FeCl 3 25 to > 500

Formation of FeOCl at 225 °C, volatilization of FeCl 3 at 307 °C, and decomposition of FeOCl at 430 °C.

[11]

II. MATERIALS AND EXPERIMENTAL PROCEDURE

Hematite used in this study had a purity higher than 99 % Fe 2 O 3 . Wüstite was prepared from hematite by its reduction using a controlled H 2 -H 2 O atmosphere. The results of scanning electron microscopy 'SEM' and that of X-ray diffraction 'XRD' analyses of these iron oxide samples are given in Figure 1. The phases identified in the wüstite and hematite samples were Fe (1-x) O and Fe 2 O 3 , respectively.

The chlorination tests of samples were performed using the TG and horizontal experimental set-up described earlier in Reference [START_REF] Kanari | Extraction des métaux de valeur des concentrés de chalcopyrite et de chromite par chloruration[END_REF]. Whilst chlorination was investigated using TG technique, a heating rate of furnace of about 25 °C/min was used. The treatment residues were analyzed by XRD.

III. THERMODYNAMIC CONSIDERATIONS

The phase stability diagram of the system Fe-O-Cl [12] is given in Figure 2. This diagram shows that at high chlorine partial pressure and low partial pressure of oxygen the most stable phase is ferric chloride. The ferrous chloride becomes a predominant specie only at low partial pressure of chlorine.

The evolution of free standard energy changes (∆G°) as a function of temperature for the chlorination reactions of wüstite and hematite using different chlorinating agents (Eqs. 1 through 10) is traced in Figure 3. Data listed in Reference [START_REF] Roine | Outokumpu HSC Chemistry for Windows, Version 2.0[END_REF] and [START_REF] Barin | Thermochemical data of pure substances[END_REF] were used for these calculations. Only FeCl 3 was considered as chloride product of the reactions of Fe (1-x) O and Fe 2 O 3 with Cl 2 and Cl 2 +CO. Equations 7 and 8 were envisaged to describe the reactions of Fe (1-x) O with FeCl 3 . The reactions generating FeOCl formation were not considered because FeOCl is expected to be decomposed at temperatures higher than 200 °C [START_REF]Anonymous : Handbook of Chemistry and Physics[END_REF] .

2/3 FeO + Cl 2 → 2/3 FeCl 3 + 1/3 O 2 [1] 8/3 FeO + Cl 2 → 2/3 FeCl 3 + 2/3 Fe 3 O 4 [2] 2 FeO + Cl 2 → 2/3 FeCl 3 + 2/3 Fe 2 O 3 [3] 2/3 FeO + Cl 2 + 2/3 CO → 2/3 FeCl 3 + 2/3 CO 2 [4] 5/3 FeO + Cl 2 + 1/3 CO → 2/3 FeCl 3 + 1/3 Fe 3 O 4 + 1/3 CO 2 [5] 4/3 FeO + Cl 2 + 1/3 CO → 2/3 FeCl 3 + 1/3 Fe 2 O 3 + 1/3 CO 2 [6] 2 FeO + FeCl 3 → 3/2 FeCl 2 + 1/2 Fe 3 O 4 [7] 3/2 FeO + FeCl 3 → 3/2 FeCl 2 + 1/2 Fe 2 O 3 [8] 1/3 Fe 2 O 3 + Cl 2 → 2/3 FeCl 3 + 1/2 O 2 [9] 1/3 Fe 2 O 3 + Cl 2 + CO → 2/3 FeCl 3 + CO 2 [10]
Figure 3 clearly shows that all envisaged reactions, except the one described by Eq. 9, are thermodynamically feasible in the whole temperature range investigated. It is important to note that ferric chloride can react with wüstite, from the thermodynamic point of view, producing higher oxides of iron and ferrous chloride (Eqs. 7 and 8).

IV. RESULTS

A. TG analysis of Fe (1-x) O and Fe 2 O 3 using Cl 2 +N 2 and Cl 2 +CO mixtures

The samples of wüstite and hematite were treated in Cl 2 +N 2 and Cl 2 +CO gas mixtures from room temperature up to 1000 °C. The total gas flow rate was kept constant and equal to 40 L/h (linear gas velocity 'V g ' ≈ 33 cm/minute). The Cl 2 /N 2 and Cl 2 /CO molar ratios were equal to 1.

Figure 4 illustrates the evolution of the percent weight loss 'pct WL' of the wüstite as a function of temperature during its treatment in the two above-mentioned atmospheres. Wüstite started to react with both chlorinating gas mixtures at temperatures lower than 250 °C resulting in a weight gain of about 20 pct at 320 °C. Above this last temperature, only weight losses were observed up to 500 °C. One may point at the curves relative to the wüstite treatments using Cl 2 +N 2 and Cl 2 +CO, they almost overlap up to 550 °C. The XRD results of the residues obtained at this temperature shows the presence of Fe 2 O 3 in both cases. This observation indicates that the global reactions of wüstite chlorination and carbochlorination up to 550 °C are those described by Eq. 3 and 6, respectively.

The weak effect of CO on the chlorination of wüstite at temperatures lower than 550 °C is consistent with observations reported by Reference [START_REF] Pap | TG studies on the reactions of hematite, magnetite and wüstite with COCl 2 and CO+Cl 2[END_REF] and [START_REF] Bertoti | Comparative thermogravimetric study of chlorinations of hematite and wustite[END_REF]

(Table I).

A plateau of pct WL of wüstite sample in Cl 2 +N 2 was observed in the temperature range from about 500 °C to 700 °C. Beyond this temperature range, the pct WL increased with temperature and about 50 pct of sample had reacted at 1000 °C. In the case of carbochlorination, the raise of pct WL of sample was observed from 550 °C and almost total reaction of wüstite with Cl 2 +CO was achieved at temperatures approaching 1000 °C.

The results concerning the treatment of hematite using the same conditions as wüstite are given in Figure 5. Hematite reacted with Cl 2 +N 2 starting from about 700 °C and 30 pct of sample had reacted at 1000 °C. While in the case of its carbochlorination, the reaction had already started at 550 °C and almost total reaction was achieved at 1000 °C. The similar trends of wüstite (see Figure 4) and hematite (see Figure 5) pct WL curves at temperature higher than 500 °C is an indirect confirmation of the formation of Fe 2 O 3 during the treatment of wüstite at low temperature despite the presence of carbon monoxide.

B. Isothermal carbochlorination of Fe (1-x) O and Fe 2 O 3

As mentioned earlier, the TG analyses were performed in non isothermal conditions with a high furnace heating rate (25 °C/min). To avoid this disadvantage, a series of isothermal carbochlorination tests of wüstite and hematite were performed between 350 °C and 550 °C. A carbochlorination gas mixture of Cl 2 +CO (Cl 2 /CO = 1) was used and the reaction time was fixed to be 2 hours. Figure 6 (a) shows the isotherm curves obtained. The similarity of two sample isotherms is clear at temperatures higher than 425 °C. XRD analyses indicated that the residue of wüstite carbochlorination at 350 °C was essentially composed of Fe 2 O 3 and Fe 3 O 4 while only Fe 2 O 3 was identified in the residue obtained at 500 °C. As expected, the carbochlorination residues of hematite at 350 °C and 500 °C contained only Fe 2 O 3 . As the weight losses at 550 °C were about 100 pct for Fe (1-x) O and Fe 2 O 3 , the final gaseous product was obviously FeCl 3 because the ferrous chloride is solid at this temperature [14] .

An attempt was made to estimate the effect of temperature on the carbochlorination of hematite by establishing the Arrhenius diagram. The average reaction rate of Fe 2 O 3 carbochlorination as function of temperature is plotted in Figure 6 (b). An apparent activation energy of about 53 kJ/mol was deduced for the hematite carbochlorination between 350 °C and 550 °C.

C. Reaction of Fe (1-x) O with FeCl 3

As mentioned earlier, the reactions of wüstite with ferric chloride (Eqs. 7 and 8) were thermodynamically favorable. On the other hand, the gaseous product of the reaction of wüstite with chlorine is ferric chloride. To investigate an eventual interaction of Fe (1-x) O with FeCl 3 a series of TG non isothermal tests were performed up to 1000 °C.

The experimental conditions and results of these tests are illustrated in Figure 7. Some treatment residues obtained at 600 °C and 1000 °C were analyzed by XRD to identify the crystallized phases. Results of XRD analysis were displayed in Figure 8. Curve 1 of Figure 7 gives the evolution of pct WL of wüstite versus temperature during its treatment in pure CO. A final weight loss of about 22 pct could be attributed to almost complete reduction of wüstite to iron.

The reaction of wüstite with Cl 2 +N 2 and Cl 2 +CO (curve 2 and 3, Figure 7) started at temperatures lower than 250 °C. The tests were stopped at 600 °C and the XRD analysis results of the obtained residues indicated that hematite was the main phases in both residues. Ferric chloride started to volatilize at temperatures higher than 200 °C in N 2 or CO and about 80 pct of the sample were volatilized at temperatures lower than 400 °C (curve 4 and 5, Figure 7). About 20 pct of the sample were volatilized at temperatures lower than 350 °C when a mixture of Fe (1-x) +FeCl 3 was treated under N 2 (curve 6, Figure 7). Then, the weight loss remained constant up to 700 °C. Beyond this temperature, the sample lost some weight up to 900 °C. XRD analysis of residue obtained at 600 °C (see Figure 8) detected the presence of Fe 3 O 4 (γ-Fe 2 O 3 ?), FeCl 2 . 2H 2 O, and Fe (1-x) O. The ferrous chloride was identified as FeCl 2 . 2H 2 O due to its hygroscopic nature. In the residue of treatment at 1000 °C, Fe 3 O 4 (γ-Fe 2 O 3 ?) was the main phase besides non-reacted Fe (1-x) O. These results compared with the evolution of the vapor pressure of FeCl 2 as function of temperature [14] suggest that ferric chloride reacted with wüstite according to Eqs. 7 and 8.

The same observations were valid for the treatment of a Fe (1-x) +FeCl 3 mixture in CO up to 600 °C (curve 7 of Figure 7 and XRD result of Figure 8). Thus, ferric chloride is able to oxidize wüstite even in presence of carbon monoxide. The sample weight loss was also observed in the temperature range of 700 °C to 900 °C, due to the volatilization of ferrous chloride. It seems that the reduction of iron oxides by CO started when the ferrous chloride was totally volatilized. As shown by curve 7 of Figure 7, another weight loss of sample was observed at temperature higher than 950 °C. The residue obtained at 1000 °C was essentially composed of metallic iron and wüstite.

V. CONCLUSIONS

Chlorination and carbochlorination of wüstite, under non isothermal conditions, started at temperatures lower than 250 °C. Hematite and ferric chloride were the final reaction product. The presence of carbon monoxide had almost not effect during the treatment up to 550 °C.

Hematite seems to be refractory, under non isothermal conditions, to Cl 2 +N 2 gas mixture up to about 700 °C. Its carbochlorination began at about 550 °C and total reaction was achieved at temperatures approaching 1000 °C.

Total reaction of wüstite and hematite producing gaseous ferric chloride was possible during their treatment in Cl 2 +CO for two hours at 550 °C. An apparent activation energy of about 53 kJ/mol was obtained for the carbochlorination of hematite between 350 °C and 550 °C.

Besides chlorine, ferric chloride can react with wüstite at low temperatures generating higher oxides of iron and ferrous chloride despite the presence of carbon monoxide. Figure 3 : Evolution of standard free energy changes as function of temperature for the reactions of wüstite and hematite with different chlorinating agents [12, 13] . Fe Figure 3 : Evolution of standard free energy changes as function of temperature for the reactions of wüstite and hematite with different chlorinating agents [12, 13] .

FIGURE CAPTIONS

(1-x) O * (hkl ) Fe Fe O 1 1) (2 0 0) (2 2 0 * ) 1 2) (1 0 4) (1 1 0) (1 1 3) (0 2 4) (1 1 6) (1 2 2) (2 1 4) (3 0 0) (2 0

Figure 1 :

 1 Figure 1 : Results of SEM and XRD analyses of used samples.

Figure 2 :

 2 Figure 2 : Phase stability diagram of Fe-O-Cl system at 800 °C[12] .

Figure 4 :

 4 Figure 4 : Thermogravimetric analysis of wüstite treatment in Cl 2 +N 2 and Cl 2 +CO gas mixtures.

Figure 5 :

 5 Figure 5 : Thermogravimetric analysis of hematite treatment in Cl 2 +N 2 and Cl 2 +CO gas mixtures.

Figure 6 :

 6 Figure 6 : Isothermal treatment of wüstite and hematite in Cl 2 +CO (a) and Arrhenius plot of hematite carbochlorination between 350 °C and 550 °C (b).

Figure 7 :

 7 Figure 7 : Thermogravimetric analysis of wüstite and FeCl 3 using different conditions.

Figure 8 :

 8 Figure 8 : XRD analyses of residues of wüstite treatment using different conditions.

Figure 1 :

 1 Figure 1 : Results of SEM and XRD analyses of used samples.

Figure 2 :

 2 Figure 2 : Phase stability diagram of Fe-O-Cl system at 800 °C[12] .

Figure 4 :

 4 Figure 4 : Thermogravimetric analysis of wüstite treatment in Cl 2 +N 2 and Cl 2 +CO gas mixtures.

Figure 6 :

 6 Figure 6 : Isothermal treatment of wüstite and hematite in Cl 2 +CO (a) and Arrhenius plot of hematite carbochlorination between 350 °C and 550 °C (b).

Figure 5 :

 5 Figure 5 : Thermogravimetric analysis of hematite treatment in Cl 2 +N 2 and Cl 2 +CO gas mixtures.

Table I .

 I Summary of the Researches Concerning the Chlorination of Wüstite.

	Chlorinating Agent	Temperature Range * , °C	Observations	R **
	CO+Cl 2	27-727	Effect of CO at temperatures higher than 577 °C.	[2]
	CO+Cl 2	27-927	Formation of FeCl 2 , FeCl 3 , FeOCl, Fe 3 O 4 , and Fe 2 O 3 . Carbon	[3]
			monoxide effect observed at temperatures higher than 577 °C.
	COCl 2	27-727	Formation of the intermediate phases : FeCl 2 , Fe 3 O 4 , and Fe 2 O 3 . [2]
	COCl 2	27-927	Formation of FeCl 2 , Fe 3 O 4 , and Fe 2 O 3 .	[3]
	CCl 4	27-927	Formation of FeCl 2 , Fe 3 O 4 , and Fe 2 O 3 .	[3]
	CCl 4	27-927	Formation of FeCl 2 , Fe 3 O 4 , and Fe 2 O 3 .	[4]
	Cl 2	25-800	Reaction occurred at low temperature. Formation of FeCl 3 and	[5]
			Fe 3 O 4 .
	Cl 2	25-800	Rapid chlorination producing FeCl 3 and Fe 3 O 4 .	[6]
	Cl 2	690-780 E a • = 140 kJ/mol.	[7]
	Cl 2	27-927	Formation of FeCl 2 , FeCl 3 , FeOCl, Fe 3 O 4 , and Fe 2 O 3 .	[8]
	Cl 2 +O 2	25-800	Oxygen had little from 0 to 40 %.	[6]
	HCl	25-800	Formation of FeCl 2 .	[6]

  . Summary of the Researches Concerning the Chlorination of Hematite.

	Chlorinating Agent	Temperature Range * , °C	Observations	R **
	C+Cl 2	25-800	Formation of FeCl 2 which volatilized as FeCl 3 .	[5]
	CO+Cl 2	27-727	Reaction started at T > 577 °C with formation of Fe 2 Cl 6 .	[2]
	CO+Cl 2	27-927	E a • = 84 kJ/mol between 327 °C and 527 °C.	[3]
	COCl 2	27-727	Reaction started at T > 427 °C with formation of Fe 2 Cl 6 .	[2]
	COCl 2	27-927	E a = 88 kJ/mol between 237 °C and 297 °C,	[3]
			E a = 60 kJ/mol between 297 °C and 447 °C.
	CCl 4	27-927	E a = 125 kJ/mol between 397 °C and 597 °C,	[4]
			n + CCl 4 = 0.5 at 527 °C.
	CCl 4	27-927	E a = 130 kJ/mol between 452 °C and 552 °C.	[3]
	Cl 2	25-800	Reaction started at 700 °C generating FeCl 3 and O 2 .	[5]
	Cl 2	700-900	E a = 136.8 kJ/mol, n Cl 2 = 0.85.	[9]
	Cl 2	25-900	E a ≈ 96 kJ/mol between 700 °C and 900 °C.	[6]
	Cl 2	597 to	E a = 188 -191 kJ/mol between 597 °C and 777 °C,	[7]
		> 777	E a = 100 kJ/mol at T > 777 °C.
	Cl 2	27-927	E a = 188 kJ/mol between 597 °C and 777 °C,	[3]
			E a = 100 kJ/mol between 777 °C and 807 °C.
	Cl 2	800-1200 E a = 167 kJ/mol, n Cl 2 = 1.	[10]
	Cl 2 +O 2	700-900	n Cl 2 = 2.75.	[9]
	Cl 2 +O 2

ACKNOWLEDGMENTS

This work was performed in the frame of contract N° BRE2-CT92-0173 thanks to the financial support of the European Union (DG-XII). The authors thank Dr. H. L. Schmidt for discussion, suggestion and help. They also would like to thank Dr. J. C. Mugica, Dr. R. Solozabal (INASMET, San Sebastian, Spain) and Dr. M. Coelho (INETI, Lisbon, Portugal) for technical discussions. Authors are indebted to Dr. A. Bonazébi, Dr. M. Djona, Dr. S. Ivanaj, Dr. N. Menad, Dr. N. Mirghaffari for discussions and help on different subjects and to Mrs. C. Richard for technical and administrative support.

Reactions of