Y Rae

A Benaarbia

J Hughes

W Sun

Experimental characterisation and computational modelling of cyclic viscoplastic behaviour of turbine steel

Keywords: Unified Viscoplasticity, Hysteresis Behaviour, High-Temperature Steel, Turbine

Fully reversed strain controlled low cycle fatigue and creep-fatigue interaction tests have been performed at ±0.7% strain amplitude and at three different temperatures (400°C, 500°C and 600°C) to investigate the cyclic behaviour of a FV566 martensitic turbine steel. From a material point of view, the hysteresis mechanical responses have demonstrated cyclic hardening at the running-in stage and subsequent, hysteresis cyclic softening during the rest of the material life. The relaxation and energy behaviours have shown a rapid decrease at the very beginning of loading followed by quasi-stabilisation throughout the test. A unified, temperature-and rate-dependent visco-plastic model was then developed and implemented into the Abaqus finite element (FE) code through a user defined subroutine (UMAT). The material parameters in the model were determined via an optimisation procedure based on a genetic solver. The multi-axial form of the constitutive model developed was demonstrated by analysing the thermomechanical responses of an industrial gas turbine rotor subjected to in-service conditions. A sub-modelling technique was used to optimise the FEA. A 2D global model of the rotor with a 3D sub-model of the second stage of the low pressure turbine were then analysed in turn. The complex transient stress and accumulated plastic strain fields were investigated under realistic thermo-mechanical fatigue loading (start-up and shut-down power plant loads). The sub-model was then used for local analysis leading to identification of potential crack initiation sites for the presented types of blade roots.

Nomenclature

C C C

Fourth Rotor, Finite Element Modelling

Introduction

The worldwide energy consumption is constantly increasing and to satisfy this demand a range of power plants is employed. Peaking power plants often use industrial gas turbine engines for electricity production. The turbine rotor operates under extreme conditions that give rise to a number of damage mechanisms due to creep-fatigue, fatigue-relaxation and thermo-mechanical fatigue solicitations. The investigation of such mechanisms is therefore crucial for understanding their effects which can lead to better designs of the turbine and improvement of its operating conditions. A typical rotor consist of two basic blocks: a compressor and a turbine. The later is split into two separate sections: high pressure (HP) and low pressure (LP) (see Fig. 1a). A section of the second stage of the low pressure turbine is shown in Fig. 1b. In HP and LP sections fir tree shaped root blades are often used. The way that the blade root is secured in the rotor groove is shown in Fig. 1c. FV566 martensitic steel (X12CrNiMoV12-3) is a common material for turbine components (e.g.

blades, disks, rotors etc.) that was specifically designed for operation within high temperature environments and elevated pressures [START_REF] Smithells | Metals Reference Book[END_REF]). A number of authors have reviewed the high temperature material performance of FV566 by a wide range of microstructural surveys, investigating cyclic softening behaviour and highlighting the importance of interaction between creep and fatigue [START_REF] Perkins | The influence of inclusions on the fatigue performance of a low pressure turbine blade steel[END_REF]; [START_REF] Turnbull | Comparative evaluation of environment induced cracking of conventional and advanced steam turbine blade steels. part 2: Corrosion fatigue[END_REF]; [START_REF] Seumangal | Influence of the heat treatment procedure on the stress corrosion cracking behaviour of low pressure turbine blade material FV566[END_REF]). Some other turbine rotor materials have also been investigated [START_REF] Holdsworth | Development of thermal fatigue damage in 1CrMoV rotor steel[END_REF]; [START_REF] Colombo | Thermo-mechanical fatigue tests on uniaxial and component-like 1CrMoV rotor steel specimens[END_REF]; [START_REF] Kupkovits | Thermomechanical fatigue of a directionally-solidified Ni-base superalloy: Smooth and cylindrically-notched specimens[END_REF]; [START_REF] He | Effect of loading rate on low-cycle fatigue properties of turbine rotor steel[END_REF]).

At low temperatures, the deformation of the majority of rotor materials can be regarded as rate independent [START_REF] Hormozi | Experimental study of type 316 stainless steel failure under LCF/TMF loading conditions[END_REF]), and their behaviour can be described by the classic plasticity theory [START_REF] Hill | A variational principle of maximum plastic work in classical plasticity[END_REF][START_REF] Hill | The Mathematical Theory of Plasticity[END_REF]; [START_REF] Prager | The theory of plasticity -a survey of recent achievements[END_REF]). At high temperatures, creep mechanisms are activated and deformation becomes rate dependent [START_REF] Marahleh | Creep life prediction of service-exposed turbine blades[END_REF]). The modelling of creep and plasticity can be performed based on unified and non-unified viscoplasticity approaches. In the unified theories, the inelastic strain represents both the rate-dependent creep (viscous part) and rate-independent plasticity [START_REF] Chaboche | Modelization of the strain memory effect on the cyclic hardening of 316 stainless steel[END_REF]). In non-unified viscoplasticity theories, creep and plasticity are considered separately and resulting deformations are treated individually, as a result the interaction between them is missing [START_REF] Contesti | Description of creep-plasticity interaction with non-unified constitutive equations: application to an austenitic stainless steel[END_REF]). Numerous viscoplastic models for metals have been developed, which can be split into two major categories:

physics-based taking microstructure evolution into account [START_REF] Nemat-Nasser | A physically-based constitutive model for bcc crystals with application to polycrystalline tantalum[END_REF]; [START_REF] Voyiadjis | A physically based constitutive model for fcc metals with applications to dynamic hardness[END_REF]; [START_REF] Abed | Constitutive modeling of nitrogen-alloyed austenitic stainless steel at low and high strain rates and temperatures[END_REF]; [START_REF] Zhang | Deformation behavior and microstructure evolution of the Cu-2Ni-0.5Si-0.15Ag alloy during hot compression[END_REF]) and phenomenological that are based on empirical observations [START_REF] Ramberg | Description of stress-strain curves by three parameters[END_REF]; [START_REF] Anand | Constitutive equations for hot-working of metals[END_REF]; [START_REF] Rusinek | Constitutive relations in 3-D for a wide range of strain rates and temperatures -application to mild steels[END_REF]; [START_REF] Arghavani | A 3-D phenomenological constitutive model for shape memory alloys under multiaxial loadings[END_REF]; [START_REF] Shutov | A viscoplasticity model with an enhanced control of the yield surface distortion[END_REF]; [START_REF] Lin | A phenomenological constitutive model for high temperature flow stress prediction of Al-Cu-Mg alloy[END_REF]). Such approaches have been used to model the mechanical behaviour of many engineering materials [START_REF] Frederick | A mathematical representation of the multiaxial Bauschinger effect[END_REF]; [START_REF] Chaboche | A review of some plasticity and viscoplasticity constitutive theories[END_REF]; [START_REF] Shutov | Finite strain viscoplasticity with nonlinear kinematic hardening: Phenomenological modeling and time integration[END_REF]; [START_REF] Besson | Continuum models of ductile fracture: A review[END_REF]; [START_REF] Drozdov | Viscoelasticity and viscoplasticity of polypropylene/polyethylene blends[END_REF]; [START_REF] Yan | A viscoelastic-viscoplastic model of titanium structures subjected to thermo-chemo-mechanical environment[END_REF]; [START_REF] Ahmed | Unified viscoplasticity modeling for isothermal lowcycle fatigue and fatigue-creep stress-strain responses of Haynes 230[END_REF]; [START_REF] He | A dislocation density based viscoplastic constitutive model for lead free solder under drop impact[END_REF]; [START_REF] Zhou | A damage-coupled unified viscoplastic constitutive model for prediction of forming limits of 22MnB5 at high temperatures[END_REF]; [START_REF] Rodas | Crystal viscoplasticity model for the creep-fatigue interactions in single-crystal ni-base superalloy CMSX-8[END_REF]). Some of them have been applied to anisothermal fatigue situations [START_REF] Zhang | Anisothermal cyclic plasticity modelling of martensitic steels[END_REF]; [START_REF] Saad | Cyclic plasticity and creep of power plant materials[END_REF]; [START_REF] Barrett | An improved unified viscoplastic constitutive model for strain-rate sensitivity in high temperature fatigue[END_REF]). In the field of high temperature materials, the Lemaitre-Chaboche viscoplastic model is commonly employed [START_REF] Chaboche | Modelization of the strain memory effect on the cyclic hardening of 316 stainless steel[END_REF]). Several authors have applied this model to study the responses of some turbine rotors under fatiguecreep conditions [START_REF] Nayebi | Analysis of unified continuum damage mechanics model of gas turbine rotor steel: Life assessment[END_REF]; Benaarbia et al. (2018a)). Moreover, multiple viscosity functions are available in literature to model the viscoplastic flow, a number of which have been reviewed by [START_REF] Chaboche | A review of some plasticity and viscoplasticity constitutive theories[END_REF]). These can include: hyperbolic sine type functions [START_REF] Miller | An inelastic constitutive model for monotonic, cyclic, and creep deformation: Part i -equations development and analytical procedures[END_REF]), power functions [START_REF] Robinson | A unified creep-plasticity model for structural metals at high temperature[END_REF]), exponential functions [START_REF] Bodner | Constitutive equations for elastic-viscoplastic strainhardening materials[END_REF]), Onera exponential function [START_REF] Chaboche | A review of some plasticity and viscoplasticity constitutive theories[END_REF]), etc.

The main objective of this study was to develop a modelling methodology for industrial gas turbine rotors under realistic loading conditions using a viscoplasticity based framework. The constitutive model used in this study has been derived from the original Chaboche-Lemaitre model [START_REF] Chaboche | On the plastic and viscoplastic constitutive equationspart i: Rules developed with internal variable concept[END_REF]). It was implemented into the Abaqus FE software by developing a corresponding UMAT subroutine. The overall stress distributions, temperature contours and viscoplastic deformations within the gas turbine have been examined by using a simplified 2D axisymmetric model of the full structure. Several locations within the rotor have been considered for detailed analysis, where high temperature and mechanical loadings are applied. Investigation of local structures, which may govern the rotor performance, has been conducted using a 3D sub-model of a second stage LP turbine, where the fir-tree blade attachment and adjacent geometry were modelled.

This paper is organised as follows. Section 2 presents the experimental observations followed by a mechanical analysis of the FV566 cyclic behaviour. Section 3 discusses the theoretical basis of the constitutive model used to conduct the FEA. It also presents the numerical implementation and the uniaxial validation of the model. Section 4 provides details of the material used, geometry of the model, thermo-mechanical loadings, boundary conditions for the FEA, followed by a discussion of the numerical results from the FEA. Concluding remarks are made in Section 5.

Experimental Observations

Material and experimental procedure

Cylindrical threaded ISO 12106 specimens of FV566 martensitic stainless steel with 5mm gauge diameter and 10mm gauge length were used for experiments (see Fig. 2a). The chemical composition of the material is given in Table 1. The material has been cut from an area close to the centreline of a service-aged gas turbine rotor which had previously been subjected to 90,000 operational hours at a maximum speed of 3000rpm and temperature of approximately 420˚C.

Fully reversed uniaxial strain-controlled low cycle fatigue tests with (dwell-type trapezoidal waveform, denoted hereafter by (DWT), as shown in Fig. 2c) and without (saw-tooth triangular waveform, referred by (SWT), as shown in Fig. 2b) tensile strain hold period were conducted at four different temperatures (300˚C, 400˚C, 500˚C and 600˚C) with a strain rate of 0.01%.s -1 , a strain amplitude of ±0.7% and a loading ratio of R = -1. The hold time was applied at the peak tension for a period of 150 minutes for the first cycle and 5 minutes for every consecutive cycle. The experiments were performed on the Tinius Olsen H25KS electromechanical testing machine (see Fig. 3) fitted with a Tinius Olsen DSCCTOL 25 kN pancake style load cell and a Severn Thermal Solutions SF2107 furnace. The temperature in the furnace was controlled using K-type thermocouples, while the tests were strain controlled using a pair of LVDT-0582 linear voltage differential transformer extensometers.

C

Si Mn Cr Mo Ni V S Fe 0.16 0.038 0.668 11.9 1.68 2.52 0.298 0.006 Remainder Microstructure of the as received material as well as fracture surfaces of failed specimens were imaged with a field-emission gun scanning electron microscope (FEG-SEM) Philips XL30 using secondary electron (SE) imaging and acceleration voltage of 20 kV. For microstructure imaging the material surface was prepared using conventional grinding and polishing steps followed by final polishing using colloidal silica, it was then chemically etched using Villella reagent.

It can be seen in the scanning electron micrograph of the as received material, Fig. 4, that the microstructure is composed of martensite. No obvious damage in terms of voids or cavities have been identified on the surface as a result of ageing.

µm

Tensile tests

Some mechanical properties of the FV566 steel were derived from monotonous tensile tests. The stress-strain data obtained from the experiments are shown in Fig. 5 for the selected temperatures of 20˚C, 350˚C, 400˚C and 500˚C. The basic mechanical properties are summarised in Table 2.

According to the mechanical responses, it is observed that the material exhibits a change in deformation behaviour when changing temperature: i) the maximum reachable load increases when the applied temperature decreases, ii) the tensile strength of 600˚C was notably lower than that of 350˚C indicating softening induced by the viscous induced-thermal effect, iii) the ultimate engineering strain increases by about 16% when going from 350˚C to 600˚C highlighting an increase in ductility. All these experimental findings will be consolidated with cyclic results to highlight most of the leading mechanisms commonly observed during the inelastic deformation of rotor steels.

[MPa]

[%] The topography of the rupture surfaces within the cross-sections were characterised using SE imaging. Overall morphology of failure for both SWT and DWT tests is shown in Fig. 6 for cross-sections of specimens tested at 600˚C. High magnification observation of the fracture surface reveals that during the crack propagation stage, SWT loading is causing beachmarks (see Fig. 6c), while DWT loading leads to noticeable voids spaced out in a similar manner (see

Mechanical behaviour

Cyclic strain-stress hysteresis loops for the three selected temperatures are plotted in Fig. 7.

Results obtained from saw-tooth loading are presented in the left hand column and dwell-type loading results are given in the right column. In each case, we investigated the hysteresis responses by introducing three mechanical parameters to characterise variations in the hysteresis loops: the change in stress, ∆σ/2, which reflects the cyclic softening (Fig. 8a), the hysteresis area, A h , which represents the deformation energy lost during the fatigue process (Fig. 8b) and the viscous stress, σ v , which expresses the stress relaxation response (Fig. 8c). The main features extracted from the experimental findings are presented below.

Cyclic softening. The cyclic softening was noticeably observed in all hysteresis responses; the change in stress amplitude decreased continuously with the increase of cycles. Three different stages were distinguished in the stress amplitude curve. These stages are characterised by their own mechanical behaviour. The first stage shows a relatively quick non-linear decrease in both tensile and compressive peak stresses (approximately 14-18% of the fatigue life) due to an increase in inelastic strain. The second stage corresponds to a steady state where the rate of the cyclic softening is roughly constant without real stabilisation of the stress amplitude. The last stage shows a sharp decrease in the stress amplitude due to viscoplastic strain and damage localisation effects, rapidly leading to the macroscopic crack as the curves at their right extremity show in Fig. 8a. Another interesting aspect shown in Fig. 8a is that the rate and duration of cyclic softening were strongly affected by the applied temperature. It is clear that the lower the applied temperature, the faster the softening rate and the smaller the number of cycles needed to reach complete softening. Both the type of loading and the temperature have a strong effect on the total number of cycles to failure. The specimens under saw-tooth loading last longer than those under dwell-type loading; the difference is increasing as the temperature increases.

Deformation energy. For the saw-tooth tests, the deformation energy (area of hysteresis loops,

A h ) increased at the very beginning of loading (the first ten cycles), then slightly decreased during the whole test and rapidly dropped off during the last fatigue stages (see Fig. 8b). For the 400˚C test, the hysteresis area seems to continuously increase until a certain number of cycles (around 55 cycles), and then decrease during the last stage of fatigue life. For dwell-type tests, all findings indicate a gradual decrease of the deformation energy during the fatigue process, except the 400˚C test where the hysteresis area seems to slightly increase during the first hundred cycles. Another observation can be seen in the final fatigue stages; it is noted that the hysteresis loops became crescent shaped with a prominent pointed tail in the compressive region. This observation has been reported for many high-temperature steels (see for example [START_REF] Saad | Cyclic plasticity and creep of power plant materials[END_REF]; Benaarbia et al. (2018b)).

Stress relaxation.

The stress relaxation response during the tensile hold time for the three selected temperatures was compared and shown in Fig. 8c. The reduction in stress during the 2.5 hour dwell period was 262.2 MPa, 273.5 MPa, 155.1 MPa and 66.2 MPa for 600˚C, 500˚C, 400˚C and 300˚C, respectively. This leads to conclude that the viscous stress relaxation magnitude during the application of load is strongly affected by the temperature; it increases with the increase of temperature. When stress relaxation values at 300°C and 400°C are compared, it can be established that influence of creep becomes increasingly significant at temperatures above 400°C. It was also noted that the stress relaxed to approximately 162.9MPa for 600˚C; a value close to what would usually be considered as the yield stress of the investigated material at 600˚C (circa 168.42MPa). This observation can still be accommodated in the standard unified elastoviscoplastic strain decomposition (i.e. = e + vp ). A final observation that was extracted from the relaxation responses is that the rapid decrease in the viscous stress was mostly happening within the first few minutes, then the stress reached a quasi-stabilised value (the quasi-saturated state of stress relaxation seems to be achieved within 5 minutes). The subsequent 5 minutes long hold period responses support this point further. It has been seen also that the magnitude of relaxed stress decreased with the increase of cycle number (this means that the relaxation response was affected by the cyclic softening due to repeated loading).

Viscoplasticity-based constitutive model

Model development

Assuming small perturbations, the total material behaviour is given by two contributions:

the elastic behaviour subjected to Hooke's law and viscoplastic behaviour that complies with classical viscoplasticity theory. In order to describe the viscoplastic deformation, it is assumed there is a convex in the space of conjugate variables inside which the behaviour is elastic and the boundary is described by a yield function. The main constitutive equations of the proposed model, restricted to a von Mises yield criterion, are given below:

= e + vp , (1a)
σ σ σ = C C C : : : -vp , ( 1b
) f = J 2 (σ σ σ, X X X) -R (T, λ(T)) -k(T), (1c)
˙ vp = λ ∂ f ∂σ σ σ = 3 2 (σ σ σ -X X X) J 2 (σ σ σ, X X X) g vp λ, X X X, R; T , (1d)
where , e and vp are the total strain tensor, the elastic strain tensor and the viscoplastic strain tensor, respectively. σ σ σ stands for the stress tensor, C C C represents the fourth-order elastic stiffness tensor while T is the absolute temperature. The function f denotes the yield surface, k represents the initial size of the yield surface while X X X and R stand for the overall back stress tensor (for kinematic hardening) and drag stress (for isotropic hardening), respectively. Following the normality assumption of viscoplasticity, the rate of viscoplastic strain tensor, ˙ vp , is normal to the tangent to yield surface while its magnitude is identical to the viscoplastic multiplier rate, λ. The second invariant, J 2 , is defined such that J 2 (σ σ σ, X X X) = 3 2 (σ σ σ -X X X) : : : (σ σ σ -X X X) while the symbol σ σ σ represents the deviatoric part of σ σ σ and ": : :" the product contracted twice. The viscous function, g vp , chosen hereafter obeys the power law form and is given such that:

g vp := z -1 (T) f m(T) if f > 0 ; 0 if f ≤ 0, (2)
where z(T) and m(T) are the viscoplastic resistance and hardening exponent functions which are both considered temperature dependent.

The back-stress tensor is employed to describe the motion of the center of the yield surface.

It enables modelling of the strain hardening rate of the mechanical response. It is considered hereafter to obey the Armstrong-Frederic form, including linear hardening and dynamic recovery coupled with an additional temperature rate term. The Chaboche-type kinematic hardening rule with two backstresses is adopted where the rate expression of each component is given below:

X X X = ∑ i X X X (i) , (3a) Ẋ X X (i) (T) = α (i) (T) 2 3 β (i) (T) ˙ vp -X X X (i) (T) ṗ + α (i) -1 (T)α (i) ,T (T) + β (i) -1 (T)β (i) ,T (T) X X X (i) (T) Ṫ, (3b)
where ṗ stands for the accumulated inelastic strain rate (assumed to be identical to the viscoplastic multiplier rate in this case), while α(T) and β(T) are kinematic hardening temperature dependent functions. The constitutive model is derived to simulate the strain hardening rate observed within the hysteresis loops, but also to account for the cyclic softening (decrease of stress range with increasing number of cycles, see Fig. 8a). Thus, the isotropic hardening drag stress, R, is considered such that:

Ṙ(T, p(T)) = (η(T) -R(T, p(T))) ζ(T) ṗ + R(T, p(T)) η(T) -1 η ,T (T) + ζ(T) -1 ζ ,T (T) Ṫ, (4)
where η(T) and ζ(T) are isotropic hardening temperature dependent functions to be determined from the change in stress amplitude during the cyclic loading at isothermal conditions.

UMAT implementation

The numerical implementation of the constitutive model is performed based on the return mapping algorithm [START_REF] Simo | Computational Inelasticity[END_REF]). The procedure is split into two steps: the first step is called elastic prediction where the viscoplastic strains and other internal variables do not evolve (the total strain is assumed entirely elastic), while the second step is named viscoplastic correction where the total strain is assumed fixed and internal variables evolve. The later step requires the numerical integration of the rate form equation of the viscoplastic strain. This is commonly done by using the Newton Raphson algorithm.

In the following, ∆t stands for the time increment and δt the Newton Raphson increment. Using Hooke's law and the unified viscoplasticity framework (only the deviatoric part of the stress tensor affects the plastic behaviour due to the fact that the hydrostatic pressure remains constant), the following expressions are derived:

σ σ σ = σ σ σ tr -C : ∆ C : ∆ C : ∆ vp , (5a)
J 2 = J tr 2 -3G ṗ∆t, n n n = n n n tr , (5b)
Φ := ṗ -( f /z) m , with f = J 2 -R -k. ( 5c
)
where G stands for the shear modulus, while σ σ σ tr , J tr 2 , n n n tr and Φ are the trial stress, trial second stress invariant, trial flow direction and viscoplastic residue, respectively.

The nullity of the residual function in Eq. ( 5c) can be used to solve the accumulated plastic strain when f > 0. This is done iteratively using the Newton-Raphson's scheme. After rearrangement, the Newton increment of the accumulated plastic strain rate can be written as follows:

δ ṗ = - m Rω + ∑ i Ω (i) n : X n : X n : X (i) Ṫ + Φ ∆t -1 z m f 1-m z m f 1-m + m∆t 3G + ζ (η -R) + ∑ i α (i) β (i) -n : X n : X n : X (i) , (6)
where

ω = η -1 η ,T + ζ -1 ζ ,T and
Ω (i) = β (i) -1 α (i) ,T + α (i) -1 β (i)
,T . The Newton Raphson increment of the accumulated plastic strain, δp, can be determined and thus updated by δp = ṗ∆t. When the norm of Φ is less than a tolerance, the iterative scheme ends and the time increment of the viscoplastic strain, ∆ ∆ ∆ vp , is updated by ∆ ∆ ∆ vp = n n n ṗ∆t while the total stress, σ, is corrected using Eq. (5a).

Model parameters identification

Material parameters were calibrated from the aforementioned isothermal tests through the following procedure:

1. The initial values of the viscous parameters were determined from the stress relaxation data in the first cycle from dwell tests. As a first estimate, it is assumed that the cyclic hardening does not significantly affect the viscous constants (i.e. m and z). Thus, the viscous stress, σ v , satisfies ln (σ v ) = m -1 ln ( ṗ) + ln (z), where the accumulated plastic strain rate, ṗ, is approximated to equal to ṗ = ˙ vp ≈ -E -1 σ (the total strain rate is held equal to zero in the stress relaxation period) with E standing for the Young's modulus. By assuming σ > X, the viscous stress, σ v , in the first long dwell period can be rearranged

as σ v = σ -X -R -k.
Plotting ln (σ v ) vs. ln ( ṗ) allows identification of both viscous resistance, z, and exponent parameter, m. The initial viscous constants at the three selected temperatures are presented in Fig. 9 (a) and Table 3.

2. The initial estimate of the isotropic hardening parameters was derived from the change in stress amplitude in the saw-tooth data. After integration of the drag stress expression (Eq. ( 4)) under isothermal conditions, the drag stress can be linked to the accumulated plastic strain through R = η (1exp (-ζ p)), where η is taken as the maximum stress difference between the first and (quasi)-stabilised cycles. In such a case, the value of the drag stress was determined as the difference between the maximum stress at the considered cycle and the maximum stress of the first cycle, while the accumulated plastic strain was considered approximately twice the viscoplastic strain range for each cycle. The initial results of the isotropic hardening parameters extracted from the identification process were presented in Fig. 9 (b) and Table 3. 3. The initial values of the kinematic hardening parameters were obtained from the hysteresis loops for different temperatures. Combining the analytical integration of Eq. (3b) and the expression of the yield function in Eq. (1c) allows to write:

σ = σ v + ∑ i α (i) 1 -exp -β (i) vp + R + k. ( 7
)
The partial derivative of σ -R in Eq. ( 7) with respect to the viscoplastic strain vp (here denoted by VPKD1), by assuming as a first approximation that the first backstress has a negligible effect on the hardening of the second backstress and that the viscous and yield stresses are independent of the viscoplastic strain, allows the determination of the second backstress parameters (α (2) and β (2) ). The differentiation of σ -R -X 2 with regard to the viscoplastic strain vp (here denoted by VPKD2) leads to the identification of the first backstress parameters (α (1) and β (1) ). By construction, β (1) is considered to be much higher than β (2) to account for both slow and rapid hardening. The results of kinematic parameters are presented in Fig. 10 and Table 3. All these values were considered as starting parameters in the optimisation process.

- Figure 10: The identification of hardening parameters at the three selected temperatures. VPKD1 and VPKD2 stand for the partial derivatives of σ -R and σ -R -X 2 , respectively, in Eq. ( 7) with respect to the viscoplastic strain, vp .

Model validation

The proposed unified viscoplastic model was thus used to simulate the temperature dependent cyclic deformations of the rotor steel under the strain-controlled saw-tooth conditions. The model uses a total of 10 material parameters to describe the hysteresis responses under uniaxial loading conditions. These parameters were determined through the optimisation procedure by considering a set of initial estimates of the isothermal tests as a staring point (see Table 3). The optimisation procedure was based on a Matlab routine that uses a genetic algorithm coupled with the Levenberg-Marquart algorithm. The main advantage of using this algorithm is to determine a global minimum for the optimisation process and avoid the local minima. The considered cost function for optimisation is expressed as follows

F c (P) = ∑ ne i=1 ∑ ned j=1 Ξ i σ num i,j -σ exp i,j 2 , ( 8
)
where σ num and σ exp are the numerical and experimental stresses, respectively. P is the parameter to be identified, ne and ned are, respectively, the number of experimental data and the number of data points sampled from the particular data set, while Ξ i is a weighting function that normalises contributions from the i th test.

Readers should note that linear relationships have been fitted to each material parameter in order to establish a set of temperature dependent functions to make easily the numerical implementation of the constitutive equations. The following expressions were then chosen for the numerical calculations E = -0.3726 [-],

α 1 = -0.06525 T + 52.352 [-]
,

β 1 = -6.1245 T + 3989.3 [MPa], α 2 = -0.4393 T -150.24 [-]
,

β 2 = -0.2914 T + 273.73 [MPa], η = -0.1678 T -9.0067 [MPa], ζ = -0.00505 T + 4.945 [-].
Fig. 11 shows the experimental and simulated stress-strain curves for the three investigated temperatures (400˚C, 500˚C and 600˚C) with a strain range of ±0.7%. Generally speaking, the proposed model represents very well the experimental cyclic behaviour. The hysteresis responses are well captured for the 600˚C test, while slightly shifted at the 100th cycle for both 400˚C and 500˚C. The hardening behaviour is well reproduced by the proposed model for most of the selected cycles. It is observed from Fig. 12 that the stress range evolution at each temperature can also be described very well by the proposed model for most of the cycles. Since damage is not considered in this work, the cyclic softening is a bit overestimated after 20 cycles for both 400˚C and 500˚C tests.

Finite Element Methodology

During the FEA of the turbine rotor, the decision on whether to use the whole 3D model or a simplified 2D model depends on many factors, such as the complexity of the geometry (e.g.

seal grooves, cooling holes, etc.), the required accuracy and the applicability of the results. The 2D modelling approach has been used by many authors for the structural analysis of turbine rotors and has shown a good effectiveness [START_REF] Banaszkiewicz | Stress corrosion cracking of a 60MW steam turbine rotor[END_REF]; [START_REF] Wang | The effect of in-service steam temperature transients on the damage behavior of a steam turbine rotor[END_REF]; [START_REF] Liu | Transient analysis and design improvement of a gas turbine rotor based on thermal-mechanical method[END_REF]; Benaarbia et al. (2018a)). Although 2D models require less time and are relatively easier to implement, they do not accurately represent the complexity of the real problem and have a number of limitations. For example, distribution of loads within the blade groove and locations of cooling holes cannot be accurately modelled in 2D space where axisymmetric conditions do not apply [START_REF] Zboinski | Numerical research on 3D contact problems of turbomachinery blade attachments in the elastic range[END_REF]). On the other hand, 3D modelling can be used to achieve accurate results in the areas of concern where the local structure and loading can be represented. A sub-modelling process is often applied to reduce the computational time, while maintaining accuracy of results. The boundary conditions used in the sub-model are thus extracted from the global model where a coarse mesh is applied. In the case of the rotor structure under investigation, a blade root is fitted in the disk groove (see Fig. 1c) and exerts a pressure to its surface during the turbine operation, leading to elevated contact pressures. Since local stress concentrations are expected to occur in these areas, a sub-model is needed to highlight the evolving of stress and strain heterogeneities.

The viscoplasticity-based constitutive model was then used to investigate the behaviour of the rotor with a sequentially coupled thermal-stress analysis using a global and sub-modelling process. During the simulation the rotor was subjected to creep-fatigue during a number of cycles of its typical loading sequence. The process of the FEA is as follows:

1. Set-up the rotor FE global model (full rotor geometry, coarse mesh, boundary conditions).

2. Run a heat transfer step to determine the temperature distribution within the global model (thermal properties and loads).

3. Perform a thermo-mechanical stress analysis of the global model using the thermal and mechanical loads, as well as temperature dependent material properties.

4. Set-up the second stage turbine FE sub-model (local geometry, fine mesh boundary conditions).

5. Conduct a thermo-mechanical stress analysis of the sub-model using temperature and displacement fields from the global model.

Global model parameters

The detailed geometry of the rotor shown in Fig. 1a was adapted for a more efficient analysis.

This was done by removing small local features such as blade root grooves and sealing structure protrusions. The extent of rotor geometry used as a global model is shown in Fig. 13a. Cyclic symmetry boundary conditions were applied to the cut sides of the slice to make a virtual 360 degree model.

Thermal loading

During its operation the turbine experiences numerous start-ups, transients due to the sudden load variations and shut-downs. At the compressor inlet, outside air is drawn in, keeping the rotor at relatively low temperature. As the air is compressed its temperature rises, while at the compressor outlet the rotor surface typically experiences temperatures of around 550˚C. A realistic typical turbine loading pattern with representative spatial surface temperature variation and rotational speed magnitude was obtained from the power plant instruments and is shown in Fig. 14. The loading pattern consists of a number of repeated cycles representing a typical daily routine of a peaking power plant with 14 operational hours and 10 standby hours. A single loading cycle that typically lasts one day comprises of a start-up, a dwell period at full load and shut-down followed by a standby period (see Fig. 14). Start-up takes about 20 minutes; at this point the turbine reaches its full speed of 3000 rpm and approximately half of the operational temperature. It takes another 5 hours for the temperature to reach its operational value. Before the shut-down, the temperature is halved approximately 10 minutes prior, then the rotational speed is gradually reduced for about 40 minutes until the rotor stops. The outer areas of the rotor in the turbine section are exposed to hot gases from the combustion chamber. This in combination with the forced cooling through the cooling holes produces a varying temperature distribution along the rotor. To determine the temperature distribution within the rotor global model, a thermo-dynamic assessment was performed by applying the spatial temperature distribution in Fig. 13a to the surface of the rotor. The compressor intake air temperature was conservatively taken to be 11 °C. The mesh set up for the rotor global model is shown in Fig. 13a.

This comprises of 16718 elements. For heat transfer analysis, the model was partitioned into 30 sections to aid the meshing process and a DC3D20 20-node quadratic brick element was mostly used. For some areas with more complex shapes, DC3D10 10-node quadratic tetrahedron elements were also used. The heat transfer step in the analysis was then used to determine the temperature field within the global model.

Mechanical loading

Fig. 13a shows a section of the rotor identifying applied mechanical loads and boundary conditions. The mechanical load comprised of the rotor centrifugal force and equivalent pressure loads at the turbine blade locations, both of which vary together with the rotational speed. For the definition of the rotor centrifugal force, the material density of 7750 kg/m 3 and a variable angular velocity up to 3000 rpm were used as inputs. The blades were not modelled explicitly as this would make the model unnecessarily complex, which are not the main focus of this work.

Instead an equivalent pressure load was applied to the blade root locations, it was calculated using the blade mass and its centre of gravity. An axial displacement constraint was applied at the far edge of the rotor. The temperature field previously found in the heat transfer step, together with mechanical loads and temperature dependent properties were then used to conduct a thermo-mechanical stress analysis of the global model. The stress analysis was performed using C3D20 20-node quadratic brick elements and C3D10 10-node tetrahedron elements.

Sub-model parameters

The second stage of the LP turbine shown in Fig. 13b was chosen for further analysis as a sub-model to get a more accurate response in the area of interest. The additional sub-model displacement boundary conditions are set at the cut areas, where the deformation from the global model would otherwise act on. Only the root part of the blade was modelled and a centrifugal force equivalent to that of the blade aerofoil was applied using a pressure load on top of the root. The blade root is fitted in the turbine groove and engaged by means of three teeth.

It exerts a pressure on the metal surface during the turbine operation. This creates a contact problem, therefore a contact interface was specified between the root and groove. A close up of the mesh near the contact interface is shown in Fig. 13c. The contact areas due to blade root teeth can be divided into six sets which are indicated as locations A1-A3 and A1'-A3'. The blade groove is skewed with respect to the turbine axis as shown in Fig. 13d, where the complexity of the blade root and groove requires 3D modelling. The low pressure turbine sub-model mesh was developed using 32246 elements for the turbine section and 20110 elements for the blade root. A penalty method contact algorithm with 0.3 friction coefficient was used for the contact simulation. Adjacent to the interface region, where high stress and strain gradients are to be expected, the mesh was refined using smaller elements. The temperature distribution for the sub-model was interpolated from the global model. Temperature and displacement fields from the global model were applied as boundary conditions at the cut boundaries of the sub-model.

The centrifugal load was applied to the whole volume of material. A thermo-mechanical stress analysis of the sub-model was then conducted.

Results and discussion

Temperature and stress field analysis

Based on the analysis described in the previous section, temperature and von Mises stress fields were obtained for the global model. Fig. 15 shows the response of the rotor after 10 loading cycles. Two different time points in the loading cycle (P1 and P2 ) were considered for the FEA (see Fig. 14). P1 is located after a period of cooling, at a point just before the start-up when the rotor is not spinning, and P2 is located at the end of the hold period at full load and maximum temperature. Although limited to just a few degrees over the surrounding temperature, some residual heat was observed at the point P1 within the inner area near the centreline and at the turbine outlet area (see Fig. 15a), with outer areas equalising with the surrounding temperature.

The von Mises stress contour at P1 shows strong presence of spatial heterogeneities in central areas of the inner and outer surfaces of the rotor (see Fig. 15c). The highest values of von Mises residual stress are mostly located near the grooves and bends due to the non-uniformity of the geometry. The main area of high stress occurs at the transition between the LP and HP turbines, as this end of the rotor has an open ended bore, it is likely caused by the bending forces during cooling due to the presence of the hollow inner bore of the low pressure turbine. The temperature contour for P2 is shown in Fig. 15b, at his point the temperature within the rotor is at its operational equilibrium. The stress contour for P2 shows that the distribution of stress field has changed (see Fig. 15d). The magnitudes of the stresses on the outer surface have reduced while on the inner section of the rotor they have intensified.

Conclusions

In the present study, fully reversed uniaxial strain controlled low cycle fatigue tests at high temperatures with and without dwell period have been conducted on a FV566 service-aged turbine steel. Changes in stress ranges, hysteresis area and viscous stress have been investigated

to characterise the cyclic mechanical behaviour of the FV566 steel. The heterogeneities start to appear as soon as the mechanical loading is applied and intensify as the cycles progress. The rate of increase was found to vary for the three locations observed, with a highest rate of increase for both the von Mises stress and plastic strain located within the fir-tree blade groove bottom contact area. The regions inspected would see effects from both creep and fatigue in service and therefore a good indication of the possible critical regions could be obtained. This could aid in identification of component areas in need of modification or redesign to help in elimination of the problematic areas.

Figure 1 :

 1 Figure 1: Gas turbine rotor components showing a) compressor, low and high pressure turbine sections, b) 2nd stage of the low pressure turbine, c) a fir-tree joint through which mechanical loads are transferred from the blade to the rotor.

Figure 2 :Figure 3 :

 23 Figure 2: (a) Geometry and dimensions of the specimen (all measures given in mm). Strain-time waveform employed for (b) saw-tooth and (c) dwell-type loading patterns.

Figure 4 :

 4 Figure 4: An FEG-SEM image of the chemically etched FV566 material in as received aged state showing typical martensitic microstructure.

Figure 5 :

 5 Figure 5: Engineering stress versus engineering strain plots showing the effect of temperature on the tensile properties of the FV566 steel.

FigFigure 6 :

 6 Fig.6d).

Figure 7 :

 7 Figure 7: Stress-strain loops for saw-tooth waveform at temperatures of (a) 400˚C, (b) 500˚C and (c) 600˚C. Stress-strain loops for dwell-type waveform at temperatures of (d) 400˚C, (e) 500˚C and (f) 600˚C.

Figure 8 :

 8 Figure 8: (a) Stress range evolution, (b) deformation energy evolution, (c) stress relaxation during 1st and 2nd cycles of the selected dwell-type waveforms.

Figure 11 :

 11 Figure11: Comparison of the experimental and simulated stress strain curves for the ∆ = 1.4% strain-controlled saw-tooth tests at the three selected temperatures. Different life stages were considered: the first cycle (left), the 20th cycle (middle) and the 100th cycle (right).

Figure 12 :

 12 Figure 12: Comparison of the numerical and experimental stress amplitude evolution for the three selected temperatures: 400˚C (left), 500˚C (middle) and 600˚C (right).

Figure 13 :Figure 14 :

 1314 Figure 13: Loading, boundary conditions and mesh for a) a rotor slice global model b) a 1/86 segment of the 2nd stage LP turbine sub-model and c) fir-tree blade root geometry with six sets of contact areas A1-A3 and A1'-A3' (pressure loads from adjacent dovetail style blade roots can be seen on either side of the fir-tree root). d) Blade root skew angle.

Figure 15 :Figure 16 :

 1516 Figure 15: Temperature and von Mises stress results of the rotor for the 10th cycle.

Figure 17 :

 17 Figure 17: Evolutions of the a) von Mises stress and b) accumulated plastic strain for the 3 locations (Loc1, Loc2 and Loc3) during the first 10 cycles.

  From test findings several observations and conclusions have been drawn: (i) The hysteresis behaviour of the FV566 steel shows cyclic hardening at the running-in stage and subsequent accelerated cyclic softening during the rest of the material life. (ii) A noticeable effect of the applied temperature and the type of loading on the lifespan of the FV566 steel has been observed (when compared to saw-tooth tests, an increase in sample life has been noticed in dwell-type tests). (iii) Relaxation and energy behaviours have demonstrated fast decreases at the very beginning of loading followed by quasi-stabilisation throughout the rest of the test. (iv) A significant stress reduction during dwell tests have been observed, demonstrating the presence of viscous effects due to high temperature conditions.A unified temperature-dependent viscoplastic constitutive model was then presented for which 10 material parameters have been derived from the experimental data. These have been identified via an optimisation procedure based on a genetic solver. The proposed model seems to account for the strain hardening, cyclic softening and relaxation behaviour. The validation of the model has shown a good correlation between the hysteresis responses from the model and experiments for the first hundred cycles. A 2D FE simulation performed on an industrial gas turbine rotor was then conducted in Abaqus using the corresponding UMAT subroutine. The results of the 2D analysis have investigated the spatial distributions of the von Mises stress before and after 10th operational loading cycles. These have indicated that after the unload the main stress concentration areas were on the rotor surface while after a period at a full load the distribution changed towards the rotor centre. A 3D sub-model at the second stage of the LP turbine was then designed to perform an accurate local study of creep-fatigue effects by looking at the spatial distributions of the von Mises stress, temperature and accumulated plastic strain contours.

Table 1 :

 1 The chemical composition of the investigated FV566 martensitic stainless steel (wt.%) determined using a spark emission spectrometer.

Table 2 :

 2 Basic mechanical properties of the investigated martensitic stainless steel FV566 (yield strength, σ y0.2 , ultimate tensile strength, σ UTS , uniform elongation, and area reduction at rupture, A r ).

	T	E	σ y0.2	σ UTS Elongation A r
	[˚C] [GPa] [MPa] [MPa]	[%]	[%]
	20	214	886	963	16	47
	350	182	705	793	11	50
	500	167	582	650	17	63
	600	140	420	456	27	90
	2.3. Cyclic tests
	2.3.1. Rupture behaviour

Table 3 :

 3 A summary of the viscoplastic material constants (determined by applying the optimisation procedure to the saw-tooth and dwell-type data). OV, IE, UB and LB are the optimised values, initial estimates, upper and lower bounds, respectively.

	Mechanical feature	Symbol Unit Temperature IE	UB	LB	OV
	Elasticity	E	GPa 400˚C	180.51	IE×50 IE×10 -2 170.63
				500˚C	166.76	IE×50 IE×10 -2 141.45
				600˚C	120.85	IE×50 IE×10 -2 96.11
		k	MPa 400˚C	380.53	IE×1.5 IE×0.5	382.32
				500˚C	269.25	IE×1.5 IE×0.5	265.31
				600˚C	168.41	IE×1.5 IE×0.5	168.32
		ν	-	400˚C	-	-	-	0.30
				500˚C	-	-	-	0.30
				600˚C	-	-	-	0.30
	Non linear viscosity	z	MPa 400˚C	818.93	IE×50 IE×10 -2 1100.05
				500˚C	932.62	IE×50 IE×10 -2 23663.01
				600˚C	963.91	IE×50 IE×10 -2 48195.51
		m	-	400˚C	66.84	IE×50 IE×10 -2 6.48
				500˚C	21.85	IE×50 IE×10 -2 2.13
				600˚C	10.31	IE×50 IE×10 -2 1.73
	Kinematic hardening α 1	-	400˚C 500˚C	256.69 124.09	IE×50 IE×10 -2 29.91 IE×50 IE×10 -2 12.41
				600˚C	168.61	IE×50 IE×10 -2 16.86
		β 1	MPa 400˚C 500˚C	13348.24 IE×50 IE×10 -2 1374.85 12682.31 IE×50 IE×10 -2 1256.32
				600˚C	1590.02	IE×50 IE×10 -2 149.96
		α 2	-	400˚C 500˚C	185.29 219.94	IE×50 IE×10 -2 17.96 IE×50 IE×10 -2 84.46
				600˚C	226.39	IE×50 IE×10 -2 105.82
		β 2	MPa 400˚C 500˚C	1426.61 1547.91	IE×50 IE×10 -2 143.51 IE×50 IE×10 -2 155.34
				600˚C	843.63	IE×50 IE×10 -2 85.23
	Isotropic hardening	η	MPa 400˚C	-60.68	IE×0	IE×10 2	-72.31
				500˚C	-107.81 IE×0	IE×10 2	-100.54
				600˚C	-112.43 IE×0	IE×10 2	-105.87
		ζ	-	400˚C	0.31	IE×0	IE×10 2	2.53
				500˚C	0.38	IE×0	IE×10 2	3.21
				600˚C	0.47	IE×0	IE×10 2	1.52

Acknowledgements

This work was supported by the Engineering and Physical Sciences Research Council (Grant numbers: EP/L016206/1, EP/N509991/1 and EP/K021095/1). The funding is provided through the EPSRC Centre for Doctoral Training in Innovative Metal Processing (IMPaCT impact.ac.uk), EPSRC/TSB IMPULSE project and EPSRC Flex-E-Plant project. The work was also partly sponsored by the RWE Generation (UK). This support is gratefully acknowledged. The authors would also like to acknowledge RWE for permission to publish this paper.

We also thank the following partners for their the valuable contributions: GE Power, Doosan Babcock Limited, Centrica plc., EDF Energy (West Burton Power) Limited, Uniper Technologies Limited, Goodwin Steel Castings Limited, NPL Management Limited, R-MC Power Recov-ery Limited, RWE Generation UK plc., Scottish and Southern Energy (SSE) plc., Siemens Industrial Turbomachinery and TWI Limited.


Study of MIMO Channel Matrices Correlation to

Optimize Resource Allocation Algorithms in Multi-Users 5G

Malo Manini, Cédric Gueguen, Rodolphe Legouable, Xavier Lagrange b-com, University of Rennes1/IRISA, IMT Atlantique/IRISA France Email: {malo.manini, rodolphe.legouable,}@b-com.com; {cedric.gueguen}@irisa.fr;{xavier.lagrange}@imt-atlantique.fr;

Abstract-Resource allocation is a major research field in wireless networks. The main challenge is to find the most suitable user for each time and frequency resource to reach the quality of service requirements. The necessity to always increase the capacity of wireless networks leads to systems using more and more antennas. Namely, the use of Massive-MIMO (Multiple Input Multiple Output) allows to simultaneously schedule several users on the same frequency resource. This is a new paradigm for resource scheduling. In this context, resource allocation algorithms have to include in their decision mechanism the choice of users association on the same time and frequency resource. A possible strategy, seen in the state of the art, is to base this user grouping on their channel-matrix correlation to minimize the interference level. In this paper we show the limited use of this approach in a system using block diagonal precoder.

Index Terms-5G, Massive-MU-MIMO, Channel correlation, Resources allocation.

I. INTRODUCTION

One of the key technologies for 5G radio transmission is Massive-MIMO (Multiple Input Multiple Output) [START_REF] Wang | Cellular architecture and key technologies for 5g wireless communication networks[END_REF]. Having a larger number of antennas at the base station can significantly improve the link budget by focusing the energy in the chosen direction [START_REF] Telatar | Capacity of Multi-antenna Gaussian Channels: Capacity of Multi-antenna Gaussian Channels[END_REF]. Focusing the energy can also be used to serve multiple users (MU) on the same time and frequency resource. This opens the way to conceiving new scheduling algorithms in order to take a real advantage from Massive-MIMO. The objective of this paper is to provide a better understanding of the MU mode in a massive MIMO context by comparing its total system rate with the one given by a traditional single-user (SU) approach.

More precisely, we compare two strategies. In the SU mode, only one user is served on a given time-frequency resource unit. There is no intracell interference. Several users are served by a classical round-robin scheduler. In the MU mode, several users are served on the same time-frequency unit. The total system rate is expected to be higher. However, there is intracell interference, which reduces the bit rate. A standard approach to estimate the interference level is to consider the correlation between the user-channel matrices [START_REF] Ji | Low Complex User Selection Strategies for Multi-User MIMO Downlink Scenario[END_REF], [START_REF] Li | MU-MIMO downlink scheduling based on users' correlation and fairness[END_REF], [START_REF] Shimbo | Control Overhead Reduction Method Employing Frequency Correlation for MU-MIMO-OFDM THP with User Scheduling[END_REF]. Our secondary objective is to check whether the correlation is a good indicator of the total system rate or not.

I II. SYSTEM DESCRIPTION

We consider one base station (BS) and several user equipment (UE) and we study the downlink transmission of one resource block (RB). The channel is assumed to be constant for all sub-carriers within a resource block. We define B as the bandwidth of the resource block, U as the set of all UEs in the cell, I as the set of all UE using simultaneously the same RB, and n t , n r the number of antennas at the transmitter (BS) and at the receiver (UE), respectively. The channel between the set of antennas at the BS and the set of antennas at UE i is modelled by complex matrix H i . The precoding matrix at the BS is represented by matrix F i and the digital combining matrix at UE i is represented by W * i . The white noise average power is denoted by σ 2 N . Note that the interference from neighbour base stations can be integrated in σ 2 N . The channel matrix of a UE i is given by H i [t] ∈ C nr×nt and is perfectly known at the eNB. Different models can be used to compute the H i [t] matrix, such as the "one-ring" model [START_REF] Shiu | Fading correlation and its effect on the capacity of multielement antenna systems[END_REF] or extended Saleh-Valenzuela model [START_REF] Gustafson | On mm-Wave Multipath Clustering and Channel Modeling[END_REF]. The total number of parallel transmissions that can be made on the same RB by one eNB is given by n t . The maximum number of transmissions on the same RB to a given UE is min(n t , n r ). The time is considered discrete and ∆t is equivalent to a timeslot. We suppose the coherence time of the channel greater than ∆t and therefore the scheduling time.

1) Transmission model: At a given time, scheduling can be viewed as an indicator function: δ i (t) where i is the UE index, k the RB index and δ i (t) ∈ [0, min(n t , n r )] gives the number of streams on resource k allocated to i.

We have the following constraint:

i∈U δ i (t) ≤ n t . (1)
The potential number of usable resources is given by n t ×K, where K is number of resource blocks. We assume that there is no limitation regarding the number of RF chains and the processing capacity. To illustrate the previous equation, we can consider a number of UEs equal to n t × K, for example n t × K = 100. In this case UEs will only have 1 RB among 100 on average per timeslot. It is then safe to say that in most cases, δ i (t) = 0.

For the sake of clarity we omit t for this part. We denote the precoding matrix associated to channel H i as F i ∈ C nt×δ i,k . The dimension of F i depends on the number of antennas at the eNB and the selected streams. The signal transmitted to UE i is :

x i = F i s i (2)
where i ∈ U and s i [k] is the δ i × 1 transmitted vector of symbols at subcarrier k. Note that s i [k] is null if δ i = 0. The global signal transmitted by the base station is:

X = j∈U ,δj >0 F j s j (3)
At the receiver, the signal is affected by the channel matrix H i , noise corrupting the received signal n of power σ 2 and I external cell interferences:

y i = H i X signal + I + n external interf erences+noise (4)
Due to the fact that we are considering a Multi-User transmission (MU-MIMO) the transmitted signal for other UEs scheduled on the same resource block has to be considered as internal interference:

y i = H i F i s i signal + j =i,δj >0 H i F j s j + (I + n) interf erences+noise (5)
Those internal interferences depend on the precoders F j orthogonality to the channel matrix H i . The signal at the receiver after combining is finally given by:

z i = W * i H i F i s i + W * i j =i,δj >0 H i F j s j + W * i (n + I) (6)
where

W * i [k]
is the digital combining matrix [START_REF] Alkhateeb | Frequency Selective Hybrid Precoding for Limited Feedback Millimeter Wave Systems[END_REF] and (.) * the conjugate transpose of a complex matrix.

2) Precoding techniques: Two precoding mechanisms are considered: a) SVD: A MIMO channel H i of a user i can be decomposed using the singular-value decomposition (SVD) as follow:

H i = U i Σ i V * i (7)
The optimal precoder in SU-MIMO for i is F i [START_REF] Telatar | Capacity of Multi-antenna Gaussian Channels: Capacity of Multi-antenna Gaussian Channels[END_REF]:

F i = V i Λ 1/2 (8)
where Λ 1/2 is a diagonal matrix of which the element λ i scales the power transmitted.

b) Block Diagonal precoding [START_REF] Boccardi | A Near-Optimum Technique using Linear Precoding for the MIMO Broadcast Channel[END_REF][10]: In MU-MIMO several users can be scheduled simultaneously on the same resource block. In such condition, users can experience interferences from other users. The objective of the Block Diagonal precoder is to eliminate UEs interferences when they are scheduled in a MU-MIMO context. Therefore, all users that are simultaneously scheduled have to be consider in the precoding process. We define HT i the concatenation of channel matrices of all users in I except i:

HT i = [H * 1 ...H * i-1 H * i+1 ...H * I ] * (9)
where I = card(I). Note that HT i is a n r × (I -1) rows and n t columns complex matrix. We are using the singular-value decomposition on HT i .

HT i = U T i Σ i [V (1) i V (0) i ] * (10)
where

[V (1) i V (0)
i ] is a n t × n t matrix. V

(1) i contains vectors corresponding to nonzero singular values and V (0) i contains vectors corresponding to zero singular values.

HT i Ṽ (0) i = Ũi Σi [ Ṽ (1)
i

Ṽ (0) i ] * (11)
The total precoding matrix is given by:

T BD = [ Ṽ (0) 1 V (1) 1
Ṽ (0) 2 V (1) 2 .
..

Ṽ (0) I V (1)
I ]Λ 1/2 (12)
where Λ 1/2 is a diagonal matrix of which the element λ i scales the power transmitted.

3) User Equipment throughput: For the MU mode, we consider the Block Diagonal (BD) precoding technique [START_REF] Khan | A simple block diagonal precoding for multi-user MIMO broadcast channels[END_REF] because the BD precoder is focused on interference management and therefore limits the reduction of the bit rate due to the intracell interference. The bit rate for UE i is given by [START_REF] Abdo | MU-MIMO Downlink Capacity Analysis and Optimum Code Weight Vector Design for 5g Big Data Massive Antenna Millimeter Wave Communication[END_REF] :

R i = B log 2 1 nr + W * i H i F i (W * i H i F i ) * W * i (σ 2 1 nr + j =i,δj >0 H i F j ( H i F j ) * )W i (13)
where 1 nr is an identity matrix of size n r .

In the SU mode, δ i is non-zero for only one value of i. Thus, there is no intracell interference and the rate C i is maximized. Equation ( 13) is simplified as:

C i = B log 2 1 nr + W * i H i F i (W * i H i F i ) * W * i (σ 2 1 nr )W i (14)
Our objective is to compare the total system rate for MU and SU strategies. With MU, n t /n r are served at the same time. With SU, only one user is served. Serving n t /n r users is done by considering n t /n r successive slots. We define n s as the number of simultaneously scheduled users, where n s ≤ n t /n r . The spectral efficiency gain χ is:

χ = E ns i=1 R i ns i=1 C i /n s , ( 15
)
where E is the mathematical expectation. We consider a large set of random deployments of terminals. Thus, H i is a random matrix. In the following, for different configurations, we study χ, which is thus our main performance indicator.

In MU-MIMO a good interference management is crucial to fully profit from the technique. As seen in [START_REF] Li | MU-MIMO downlink scheduling based on users' correlation and fairness[END_REF] and [START_REF] Ji | Low Complex User Selection Strategies for Multi-User MIMO Downlink Scenario[END_REF] the channels matrices correlation is commonly used as an indicator for interferences level between UEs. The correlation of two UEs channel matrices indicates the interference level if they are using simultaneously the same resource block. The correlation of two UEs is defined by [START_REF] Czink | Spatial separation of multi-user MIMO channels[END_REF]:

ξ(i, j) = |tr( H i H * j )| || H i || F || H j || F (16)
where (.) * denotes the conjugate transpose operation and ||.|| F the Forbenius norm of a matrix.

III. SIMULATION SET UP

Two user equipment (UE) are considered. The channels of the UEs are mutually independent. A channel for a user i is generated using the extended Saleh-Valenzuela model [START_REF] Gustafson | On mm-Wave Multipath Clustering and Channel Modeling[END_REF] where the obstructed-line-of-sight (OLOS) parameters are given in table I We consider one sector in a typical 3-sector configuration: the angle between the terminal and the base station is between 0 • and 120 • horizontally and between -45 • and 45 • vertically. We assume that the transmitter has a full knowledge of the channel for each UE. Each antenna settings is tested over an average of 1 million configurations, where each configuration is a new independent set of random variables.

IV. IMPACT OF THE NUMBER OF ANTENNAS ON THE CORRELATION DISTRIBUTION

The number of antennas is a main strength of Massive-MIMO. The more antennas there are, the more precise the spatial separation will be. However, due to technical limitations, such as precoding complexity or room space available for mobile devices, the number of antennas is limited. 3GPP fixes the number of antennas at the transmitter and at the receiver in the system evaluation framework. In our performance evaluation, we study different configurations to understand the impact on the correlation. In the following example we show the impact of the number of antennas, at the transmitter n t and at the receiver n r , on UEs correlation repartition. For each sample we compute the channel correlations of two UEs.

Figure 1 shows the empirical PDF of the correlation. The CDF can be easily deduced from the PDF but is not shown due to space constraints. Several observations can be highlighted configurations n r = 2&n t = 64 and n r = 8&n t = 16 have a nearly similar correlation distribution. Having a greater number of antennas allows the system to have more narrow beams, resulting in a lower average correlation. With n r = 8&n t = 64 there is a larger number of low correlation than with n r = 2&n t = 16.

V. IMPACT OF THE CORRELATION ON THROUGHPUT

The main objective when designing a resource allocation algorithm is to increase the capacity of the system. To understand if the correlation is a good indicator for a scheduling process, we study the correlation impact on throughput. For each sample, we calculate the throughput gain from the MU-MIMO allocation, that is given by 15. Figure 2 shows χ values, computed in (15), for different antenna configurations. Each antenna configuration is represented tree boxes with 0.1, 0.2 and 0.3 correlations. Concerning the throughput gain, for example, with n r = 2&n t = 16, between 0.1 and 0.3 correlations, less than 3% of capacity is lost on average. The Block Diagonal precoding technique is built to withstand the increase of interferences. When the correlation increases, even if the spectral efficiency decreases, the precoder is able to contain the degradation. With an antenna configuration n r = 2&n t = 16 at a correlation of 0.2 they may be more than 10% between the highest and lowest value of throughput. This difference depends on the degree of freedom available in the precoding process. To reduce the deviation the ratio nt nr should be larger. The main issue in this example, is the ability to predict the gain. In all configurations there is an intersection between throughput values, meaning that the same throughput gain can be experienced with different correlation values. Note that in the case of n r = 8&n t = 16 and n r = 4&n t = 8, it is not even profitable to schedule two UEs at the same time. They will both experience a smaller throughput than if they were scheduled in Single-User mode.

VI. MOBILITY IMPACT ON THE CORRELATION BETWEEN UES

In order to be able to use the correlation as an indicator, the correlation should be stable in frequency and in time when users are moving. The more variable the channel, the more frequent the correlation should be updated. There is thus a feasibility issue: requiring more Channel State Information feedbacks and more CPU resources. As seen in the simulation setup, we are considering two users with independent channel matrices. The results are computed on one snapshot in each configuration. One user is static during the simulation, while the other user is moving and its channel adjusts accordingly. The purpose of this experiment is to determine if the correlation is stable enough to be calculated as often as possible for two UEs during a period of time and over a given bandwidth.

Figure 3 shows the correlation evolution between two users depending on the frequency. Over 100 MHz bandwidth, the correlation significantly varies in amplitude but slowly in time. This shows that the variation of phase offset is not negligible over a large bandwidth. However, the rate variation is slow enough and only a few correlation calculations are needed by timeslot. Another variable factor is the displacement over time. Figure 4 shows the correlation evolution when a user moves from its original position. Distance is represented in λ, which is the wavelength, and 10λ at 26 GHz is travelled in roughly 84 ms. We can observe a high variability of the correlation with the distance changes, with a similar amplitude gap compares to figure 3. When the variation in frequency and distance are combined, the correlation stability might not be sufficient for an accurate calculation between scheduling decisions.

VII. CONCLUSION

Correlation is often used in the user selection process as an indicator to maximize the global throughput. In this study, we show the correlation repartition, the relation between correlation and throughput and the mobility impact on the correlation when using a Block Diagonal precoder. The repartition of the correlation depends on the product of the number of antennas at the transmitter and at the receiver. A higher product reduces the probability of high correlation values and diversity. The ratio between the number of receive and transmit antennas determines the variability of the throughput gain. A large number of transmitting antennas compared to number of receiver antennas makes the throughput gain more predictable. Finally, the correlation variability quickly changes over frequency and time, when a user is mobile. There is a substantial number of correlation calculations needed over time. As a summary, a user selection process, using Block Diagonal precoding, can barely benefit from a correlation indicator in different configurations, due to high variability with minimum effect on the throughput. MU-MIMO allows more than two users sharing the same resource block. In future works, we will study how a higher number of simultaneously scheduled users impacts the correlation and the corresponding throughput.

Fig. 1 :Fig. 2 :

 12 Fig. 1: Correlation repartition depending on antenna numbers

Fig. 3 :Fig. 4 :

 34 Fig. 3: Correlation over 100 MHz bandwidth at 26GHz

  .

	Λ	5	ns	Cluster arrival rate
	λ	1	ns	Ray arrival rate
	Γ	8.7	ns	Cluster decay rate
	γ	4.7	ns	Ray decay rate
	σ	0.1	rad	Intra-cluster angles standard deviation

TABLE I :

 I Channel model parameters, from[START_REF] Gustafson | On mm-Wave Multipath Clustering and Channel Modeling[END_REF]