Yoann Robert

Alban Leroyer
email: alban.leroyer@ec-nantes.fr

Sophie Barré

• Patrick Queutey

Michel Visonneau

Validation of CFD simulations of the flow around a full-scale rowing blade with realistic kinematics

Keywords: Computational fluid dynamics, Fluid-structure interaction, Rowing, Validation

This article deals with the validation of the modelling and numerical simulation of a rowing stroke, by means of CFD. Simplified but realistic strokes were performed in a towing tank with a rotating arm and a real flexible oar. Those laboratory conditions are better controlled than those of in situ trials. An FSI procedure is developed to take into account the oar bending, which is essential in the physics of this flow. The results show that this numerical framework is able to reproduce qualitatively the real flow including the breaking of the free surface around the blade and the transport of the air cavity behind it. The profiles of forces are well reproduced, with propulsive forces overestimated by 5-12% for their maxima. The study also focuses on the computation of the uncertainties. It is highlighted that, even for this well-controlled experimental equipment, the uncertainties on the quantities of interest are of about 11%. In other words, the experimental uncertainty covers the numerical errors. So, this numerical modelling is validated and can be used for design and optimisation of blades and oars, or to contribute to the better understanding of the boat-oar-rower system and its dynamics.

Introduction

Over the last three decades, CFD simulations have been increasingly used in the industry and the academic research. This trend is explained by the ongoing development of the physical models and the numerical methods, allowing the results to be accurate enough to compete with experimental means such as wind tunnels and towing tanks. Together with these scientific advances, IT resources have been growing to the pace of Moore's law, which has eased CFD usage. This tool has also been used intensively in sports engineering [START_REF] Hanna | CFD in sport-a retrospective; 1992-2012[END_REF], especially in Formula One or sailing, where budgets allocated to the development are very high. In sports like rowing, where funding is less substantial, progress goes slower. However, some promising work has been produced.

For the non-specialist reader, the following definitions are given. A rowing stroke is a cycle composed of two parts: the propulsion, where the oar blade is immersed, and the recovery, where it is in the air. The propulsion can also be divided into several phases: the catch, where the blade rapidly goes through the free surface (i.e. the air-water interface), the drive, where the rower transmits the power in a very short time, and the catch or finish, where the blade is removed from the water. We refer to a real rowing stroke when the cycle given above is respected and the elements of the following list of physically relevant considerations are at least taken into account: the blade kinematics is imposed by a human rower; it is modified due to the oar flexibility and the FSI; the free surface is well described; the dynamics of the boat-oars-rower system is well represented.

In rowing, the boat and the oar blades are crossing the free surface, and the flow around those two bodies are complex and atypical in naval hydrodynamics. The flow around the blade is particularly challenging to model. In fact, the blade is moving with six degrees of freedom (DOFs). During the catch, the blade breaks the free surface and possibly projects drops at the impact. Then, during the drive, a ventilation phenomenon appears, leading to some air trapped behind the blade and at its extremity. This air cavity follows the blade during almost all the rowing stroke. The oar bending is also a key element of the flow and must be handled with care. Indeed, the deflection modifies the incident angle defined as the angle between the oar axis and the direction of the velocity of the water relative to the blade. Another characteristics is the violent dynamics of the stroke, where the incident angle quickly varies.

The scope of the following state of the art is limited to the research works related to numerical simulations. However, it is not so restricted, since most of the quoted articles refer to experimental works chosen as configurations for the necessary step of validation of the numerical model, as the present article does.

The first reference showing the usage of CFD for the study of the flow around the blade dates back to 1993. Videv and Doi [22] performed simulations with a two-dimensional monofluid configuration. As the blade is fully immersed and the plane of study is horizontal, this results in a profile comparable to a hydrofoil. With this very simple approach, the authors conclude on the importance of capturing the very unsteady nature of the flow. Despite this highly important consideration, a lot of authors consider that the study of the flow can be done with a sequence of discrete stationary orientations of the oar. For instance, Coppel et al. [START_REF] Coppel | Simulating the fluid dynamic behaviour of oar blades in competition rowing[END_REF] use the RANS solver ANSYS Fluent to run simulations essentially of stationary nature, without free surface, using the threedimensional geometries of a flat blade and a Big Blade. The rotation of the blade is done with a cylindrical sliding mesh. As the configuration radically simplifies the physical problem, the numerical results for the drag and lift coefficients get far from the experimental ones of Caplan and Gardner [START_REF] Caplan | Optimization of oar blade design for improved performance in rowing[END_REF]. Sliasas and Tullis [START_REF] Sliasas | The dynamic flow behaviour of an oar blade in motion using a hydrodynamics-based shell-velocity-coupled model of a rowing stroke[END_REF] also replicate the same case using the same solver but taking account of the free surface and therefore obtain better results.

Kinoshita et al. [START_REF] Kinoshita | Rowing velocity prediction program with estimating hydrodynamic load acting on an oar blade[END_REF] were the first ones to carry out threedimensional unsteady computations for this flow. However, the configuration was simple and unrepresentative of a real rowing stroke since the blade rotates around a vertical axis passing at mid-chord, with a uniform flow velocity and without free surface. Sliasas and Tullis [START_REF] Sliasas | Numerical modelling of rowing blade hydrodynamics[END_REF] also performed unsteady computations using ANSYS CFX and demonstrated that the maxima of unsteady forces are significantly higher than the steady ones. Yet the results are unphysical at some point, especially for the propulsive force which is momentarily and unexpectedly negative. This default was corrected in another paper published later [START_REF] Sliasas | A hydrodynamics-based model of a rowing stroke simulating effects of drag and lift on oar blade efficiency for various cant angles[END_REF]. Leroyer et al. [START_REF] Leroyer | Experimental and numerical investigations of the flow around an oar blade[END_REF] performed simulations reproducing a simplified rowing stroke but with a dynamic similar to a real one, especially for the free surface deformation. In this case, the oar is rigid and the blade is flat and at scale 0.7. The results are in good agreement with experimental values. The authors continue their analysis of the flow [START_REF] Leroyer | Influence of free surface, unsteadiness and viscous effects on oar blade hydrodynamic loads[END_REF] and conclude that the free surface and the unsteadiness are of paramount importance but viscous contributions are negligible. They also set up a modelling of a real rowing stroke [START_REF] Leroyer | Validation of CFD simulations for unsteady violent free surface flows : the case of the hydrodynamics around rowing blades[END_REF], also used by Robert et al. [15]. The kinematics of the blade comes from in situ recordings, i.e. where an instrumented boat was used by an elite rower on a river. Those measurements serve as input data for the simulations. The oar flexibility is also taken into account, with linear and angular deflections proportional to the main bending moment. In this strong coupling-here for the first time in the rowing literature-the position and the orientation of the blade are updated at each non-linear iteration of the flow solver ISIS-CFD, in function of the hydrodynamic forces and moments acting on the blade. The obtained results show that the dynamic of the rowing stroke is well captured, despite a difference with the experimental results for the propulsive moment which grows quickly at a given angular interval and reaches up to 15-20%. This deviation might come from an inaccurate capture of physical phenomena (transport of the air cavity attached to the blade, flow separation, etc.). However, since such a default does not appear in the present results, it is more likely due to higher uncertainties on the measurements aboard the instrumented boat and because some DOFs are not measured.

The aim of the article is to show that the numerical modelling that we initiated [15] and improved is able to finely reproduce a realistic rowing stroke obtained in laboratory conditions, where uncertainties on measurements are smaller than in situ trials.

Methods

Experimental configuration

In 1999 an experimental campaign was carried out in the towing tank of the LHEEA Laboratory. The experimenters have designed accurate, controlled and repeatable test means to generate a simplified rowing stroke with a realistic kinematics. An arm allows the rotation of the oar and the carriage of the towing tank is used to produce a uniform translation (see Fig. 1). This configuration enables to decouple the study of the propulsive system (i.e. the oar blade) from the motor (i.e. the rower). Therefore one can observe, characterise and model the hydrodynamic forces and moments in function of several parameters.

The towing tank had a length of 70 m in 19991 , a width of 5m and a depth of 3m . The length is sufficient to perform a rowing stroke, considering the acceleration and deceleration phases of the carriage. Given the length of the oar, the axis of rotation is shifted from the midplane to increase spacing between the blade and the vertical wall. The blade chord is small enough relative to the water depth to neglect the vertical confinement.

A six-component dynamometric platform is used to measure all the forces and moments and is installed on the carriage. Strain gauges are located on the sensitive part of the platform. The internal forces and moments of the transmission are not taken into account. In this way the hydrodynamic, aerodynamic and inertial forces are measured. Specific tests are carried out with the oar in the air to model aerodynamic and inertial forces. Those force models are verified by comparing with a trial performed in the air using a kinematics identical to one of those chosen for the trials (oar in water). Using those models, the hydrodynamic forces are retrieved by subtracting the aerodynamic and inertial forces from the measured forces.

A rotating arm, whose frame is lightweight and made of a carbon-fibre composite material, was specifically designed for those trials. Like in the real case, the oar is linked to this device at two points: one at the handle location and the other one at the oarlock location. The oar orientation is defined using the Cardan angles with the configuration given in Fig. 2. Three successive rotations enables to define the oar orientation. The first rotation of an angle 0 called sweep angle (classically called yaw angle) and around the vector changes the original basis ( , , ) to the basis ( , , ) , with = . The second rotation of an angle 1 called immersion angle (classically called pitch angle) and around the vector changes the basis ( , , ) to the basis ( , , ) , with = . Finally, the third rotation of an angle 2 called roll angle and around the vector changes the basis ( , , ) to the basis ( , , ) , with = . The roll angle 2 and the immersion angle 1 are constant during the trial and the sweep angle 0 varies from about 0 to 180

• .
The geared motor unit consists of a high-specific-power brushless motor. A digitally controlled drive imposes a 200-Hz-sampled speed command to the motor. An example of the laws imposed in the trials, composed of the five phases explained hereafter, is given in Fig. 3. As the carriage cannot be vertically moved during the test, the blade is immersed from the beginning to the end of the stroke. Furthermore, the trial cannot begin with the oar already oriented at the attack angle because it would cause unsuitable perturbations. So the phase I starts with the oar almost aligned to the carriage direction, at an angle earlier found that does not generate any transversal force, and the angular velocity is such that the least possible flow perturbations are created. The phase II models the catch, corresponding to a rapid increase of the velocity at 0 = 35

• . The phase III is performed at a constant kinematic efficiency (detailed in the Sect. 4). The phase IV is a transition to the finish, in which the blade should get out of water in a real stroke and where the velocity rapidly decreases at 0 = 120

• . The phase V is done at a moderate velocity to avoid needless stress on the structure and the sensors.

The fact that this rowing stroke is simplified is because the blade is always immersed and because the kinematics is not imposed as a human would do: it is not the same kinematics and the motor imposes the oar orientation regardless of the hydrodynamic forces applied on the blade (i.e. there is no dynamic feedback reaction, as a rower is naturally submitted to).

Different values of the immersion and roll angles are tested, as well as for the carriage velocity V c . The list of the trials used in this article is presented in Table 1. The variable c 2 corresponds to the opposite of 2 , to match the usual convention in rowing. The vertical reference is such that, when the oar is oriented with 2 = 0

• and 1 = 9.5

• , the top edge of the blade is awash.

A 12-bit Keithley DAS1602 computer board is used to record the signals of the three forces and three moments at a given point, the angular position and velocity, with a sampling rate of 500 Hz. As the recording is subject to noise of different natures coming from the experimental environment (e.g. the platform and rolling noise), a low-pass filtering with a cutoff frequency of 11 Hz is realised.

Numerical configuration

The fluid solver used is ISIS-CFD. It is developed by the METHRIC2 group of the LHEEA Lab and available as a part of the FINE™/Marine computing suite, which is dedicated to marine applications. This is an incompressible unsteady Reynolds-averaged Navier-Stokes (RANS) solver, based on a fully unstructured finite-volume method to build the spatial discretisation of the conservation equations. The pressure-velocity coupling is obtained through a Rhie & Chow SIMPLE-like method: at each time step, the velocity is updated from the momentum equations and the pressure is given by the law of mass conservation, transformed into a pressure equation. An arbitrary Lagrangian Eulerian (ALE) formulation is used to take into account the modification of the fluid spatial domain [START_REF] Leroyer | Experimental and numerical investigations of the flow around an oar blade[END_REF]. It is associated with robust and fast grid deformation techniques [START_REF] Leroyer | Numerical methods for RANSE simulations of a self-propelled fish-like body[END_REF]. The temporal discretisation scheme is the backward difference formula of order 2 (BDF2) when dealing with unsteady configurations. For each time step, an inner loop, also called non-linear loop, associated to a Picard linearisation is used to solve the non-linearities of the system and to converge all the sequential coupled equations. Free-surface flows are addressed with interface-capturing methods, by solving a convection equation for the volume fraction of water, which is discretised with specific compressive discretisation schemes [START_REF] Queutey | An interface capturing method for free-surface hydrodynamic flows[END_REF]. The code is fully parallel using the message passing interface (MPI) protocol. An automatic adaptive grid refinement technique [START_REF] Wackers | Adaptive grid refinement for hydrodynamic flows[END_REF] and a sliding grid method are also included.

The computational domain is a rectangular cuboid of dimensions 5m × 5m × 3m whose centre initially coincides with the blade root, i.e. the intersection of the shaft and the oar blade, as illustrated in Fig. 4a. This domain follows rigidly the body of the oar. It is large enough since the use of a larger domain has given identical results. This is explained by the fact that a rowing stroke takes place in a relatively confined space, as it can be seen in Fig. 6 or Fig. 2 of [15]. Indeed, the blade trajectory draws a flattened loop. The free surface location is by definition the z = 0 m plane. All the meshes are created with the software Hexpress™. As it can be seen in Fig. 4b,c, the mesh is refined vertically around the free surface and in a more isotropic manner inside a box surrounding the blade. To perform a grid convergence study three levels of mesh fineness are tested. The approximate numbers of cells for each grid are given in Table 2.

The rigid body motion is planar because the axis of rotation of the arm is perpendicular to the translation vector of the carriage. The reference point chosen for the rotation is the blade root. The temporal evolution of the sweep angle comes from the integration of the angular velocity given in Fig. 3. As the blade root is not located on the axis of rotation, it means that two translations have to be imposed to the blade: one colinear to and another one along . Before the stroke generation, a transient phase has to be added. It contains a 1-s acceleration of the carriage with an adjustment of the oar orientation3 and it is followed by a 1.5-s phase to make the flow converge.

A Dirichlet condition is applied on the four vertical boundaries of the domain, where a null velocity is imposed (i.e. a far field condition). For the two faces at the top and the bottom, another kind of Dirichlet condition is used: the pressure is set to its atmospheric value. A wall function is applied on the blade and the shaft.

The temporal discretisation scheme is the backward differencing scheme of second order (BDF2) to deal with this unsteady configuration. To ensure the stability of the FSI coupling, the convergence of the non-linear loop is achieved when a reduction of 4 orders of the residuals is obtained or, by default, when 15 non-linear iterations are performed. The pressure is solved using the multigrid BoomerAMG solver [START_REF] Henson | BoomerAMG: a parallel algebraic multigrid solver and preconditioner[END_REF], also with 4 orders, or 50 iterations to the maximum. For the simulation of a real rowing stroke, the use of a Euler model is sufficient because the viscous effect is negligible, as shown by Leroyer et al. in [START_REF] Leroyer | Influence of free surface, unsteadiness and viscous effects on oar blade hydrodynamic loads[END_REF]. However, for the particular rowing stroke treated in this article, the viscous effect is not negligible anymore for the first and the last phases. Indeed, for the phases I ( 0 close to 0

• ) and V ( 0 close to 180 • ), the oar is approximately aligned with the direction of the carriage and sees a large magnitude of incident velocity. The k -SST turbulence model of Menter [START_REF] Menter | Zonal two-equation kturbulence models for aerodynamic flows[END_REF] is then used.

In addition to having a sufficiently fine mesh to capture the free surface around the blade, the Courant number also has to remain below 0.3 so that the interface-capturing BRICS scheme used for the resolution of the volume fraction keeps full compressivity properties. This constraint implies that the time step is limited and the CPU time increases consequently. A specific time-integration method [START_REF] Robert | Numerical simulation and modelling of tridimensional free-surface flows[END_REF][START_REF] Robert | Time-splitting procedure for the resolution of the volume fraction for unsteady flows[END_REF], which is called sub-cycling algorithm, was developed in order to speed up unsteady computations. It is based on the resolution of the convection equation for the volume fraction in several successive steps in which the associated time step (and thus the Courant number) is smaller. The other equations remain unchanged and are solved using a larger global time step. Since the resolution of the volume fraction has a very limited CPU cost, this technique enables to decrease the total computation time, without loss of accuracy. Validations were carried out on two-and three-dimensional test cases. It has lead to the conclusion that, with this procedure, the CPU time is reduced by a factor usually between 3 and 4.

Structural model of the oar and FSI coupling

Structural model of the oar

The flexibility of the oar shaft is taken into account. Thus the imposed kinematics of the blade is modified and depends on the fluid loads which apply on this one, according to a law given hereafter. In return the flow is also altered, so are the hydrodynamic forces. The FSI coupling which has been developed can be qualified as strong. Indeed there is a convergence loop in the FSI coupling, coinciding with the resolution of the fluid solver. The variables are exchanged and updated at each non-linear iteration [START_REF] Leroyer | Numerical methods for RANSE simulations of a self-propelled fish-like body[END_REF][START_REF] Yvin | Added mass evaluation with a finite-volume solver for applications in fluidstructure interaction problems solved with co-simulation[END_REF].

The inertial effects of the oar are small compared to the fluid forces acting on the blade, that is why a resolution of the static equilibrium of the structure is suitable. The linear (f) and angular ( ) deflections are only considered in the (D, , ) plane associated to the frame linked to the oar. Deflections outside this plane are far weaker therefore neglected. Two models have been tested for the modelling of the oar bending, which is described in Fig. 5. The first one, denoted as simplified model, uses a proportionality between the bending moment at the sensors location (point C, merged with point B in this configuration) and the linear and angular deflections, as given by the system (1). The second model, designated by complete model, comes from laws of the beam theory and links directly the hydrodynamic force F and moment M (at the point D) to the deflections, as shown by the system (2). The oar is fixed to the rotating arm at two points A and B representing, respectively, the handle and the oarlock locations, separated by a distance L i = 0.800 m , also known as inner lever. The distance between the axis of rotation and the blade root (point D) along the oar axis is L e = 1.558 m . The point E is a rough approxi- mation of the point where the hydrodynamic force is applied and it is located at a distance l = 0.258 m of the point D. In fact it is mostly the point where loads were applied on a test bench for the determination of the coefficients K lin = 3.4241 ⋅ 10 3 N and K ang = 2.948 ⋅ 10 3 N ⋅ m.

The coefficients of the complete model depend on geometrical and mechanical parameters, according to equations (3)

(1)

f = 1 K lin M fz (C) = 1 K ang M fz (C) (2) f = aM(D)+bF = cM(D)+dF (3) a = 1 (EI G z ) lin L e 2 + L i 3 L e , b = 1 (EI G z ) lin (L e + L i ) L 2 e 3 , c = 1 (EI G z ) ang L e + L i 3 and d = 1 (EI G z ) ang L e 2 + L i 3 L e ,
where EI G z is the product of the Young's modulus by the second moment of area along the z and at the point D. It is also called flexural rigidity and its value along the oar axis was studied by Laschowski et al. [START_REF] Laschowski | Modelling the deflection of rowing oar shafts[END_REF]. Different values are chosen for the linear or angular deflections for a better fitting:

(EI G z ) lin = 4.3872 ⋅ 10 3 N ⋅ m 2 and (EI G z ) ang = 3.3778 ⋅ 10 3 N ⋅ m 2 .

FSI coupling

This FSI coupling cannot be used as is. Because of the added mass effect, a part of the fluid forces on the blade depends on the acceleration and may lead to the divergence of the numerical coupling [START_REF] Söding | How to integrate free motions of solids in fluids[END_REF]. To keep it stable, the position of the blade is updated at each non-linear iteration including a physical relaxation depending on an accurate computation of the added mass, whose procedure comes from Yvin et al. 25] but adapted for the quasi-static approach used here. The computation of the 6 × 6 added-mass matrix MA is done through the resolution of a Laplace operator to obtain the real added mass.

Stabilisation using the added mass for the FSI coupling with the simplified model

For the simplified model, defined by (1), we modify two DOFs of the blade, and , but both are proportional. The stabilisation of one of them leads to the stabilisation of the other one. So, we can work only on the stabilisation of one DOF.

According to the principle of angular momentum on the blade, written at the sensor location (point C) and in the basis associated to the blade, we have the following relation:

where I is the moment of inertia along the third axis of the basis, M z the moment of the hydrodynamic forces at the point D along the same axis and M in is the moment of iner- tia due to the non-Galilean nature of the considered frame. Furthermore, the fluid forces and moments are composed of two parts, as it appears in (5): the first one is independent of (4)

I ̈𝛼 = M z -K lin f + M in ,
Fig. 5 Modelling of the bending of the oar shaft the acceleration (written with the overline), the second one depends on it through the added-mass matrix MA.

Here we consider that the component dependent on the added mass mainly responds to f , which is here equal to ÿ , i.e. the acceleration normal to the blade. In other words, all the terms of MA except MA 22 are neglected. This hypothesis is linked to the model, based on a deformation of the oar due to an orthogonal force applied on the blade. Besides, the model is expressed using the moment M z . So, we choose a value MA 62 = MA 22 L e and the system (5) becomes As previously mentioned, inertial terms are negligible. As a result, I ̈𝛼 and M in can be set to zero in [START_REF] Hanna | CFD in sport-a retrospective; 1992-2012[END_REF]. Using the quasi- static approach, this equation can be written as follows, with n 1 being the temporal iteration, k the previous iteration of the non-linear loop and k + 1 the actual one:

To make the writing clearer, n + 1 is omitted all quantities are expressed at the corresponding temporal iteration. We add to both sides of the last equation the added-mass term MA 22 L e f , at two different values of non-linear itera- tion to decrease the dependency of the right-hand side on the acceleration: This manipulation stabilizes the FSI coupling, by acting as physical relaxation of the phenomenon (see equation [START_REF] Leroyer | Numerical methods for RANSE simulations of a self-propelled fish-like body[END_REF]). It recovers the initial equation ( 7) at convergence. Since the BDF2 is used for the time integration, the increment of position f k+1 = f k+1 -f k is linked to the increment of accelera- tion 𝛿 fk+1 = fk+1 -fk . The link is a first-order approximation on the time step, expressed as follows:

where e c is a coefficient of BDF2 and its value is e c = 3 2 t , where t = t n+1t n . Then, equation ( 8) can be rewritten as follows:

(5)

⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ F x F y F z M x M y M z ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ F x F y F z M x M y M z ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ - ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝
⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ẍ ÿ z rx ry rz ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ , (6) M z = M z -MA 22 L e f . (7) K lin f n+1 k+1 = M z n+1 k . ( 8
) K lin f k+1 + MA 22 L e f k+1 = M z k + MA 22 L e f k . (9) 𝛿 fk+1 = e c 2 𝛿f k+1 , (10)
K lin f k+1 + e c 2 MA 22 L e f k+1 = M z k + e c 2 MA 22 L e f k .
Using the temporary value of f, denoted f tmp and obtained with the expression of the model (1) such as f tmp = M z k ∕K lin , we can rewrite [START_REF] Leroyer | Validation of CFD simulations for unsteady violent free surface flows : the case of the hydrodynamics around rowing blades[END_REF], after some manipu- lations, in a shorter way:

with the relaxation coefficient r given by Using the model [START_REF] Leroyer | Validation of CFD simulations for unsteady violent free surface flows : the case of the hydrodynamics around rowing blades[END_REF], the new value of the angular deflections is given by the relation For the first non-linear iteration (i.e. for k = 1 ), of each tem- poral iteration, the values f 0 and 0 come from an extrapola- tion at the second order of the converged values at the two previous time steps.

Stabilisation using the added mass for the FSI coupling with the complete model

In the complete model, f and are independent. The system (2) that defines them in function of F y and M z can be written under a matrix form and inversed, such as with We consider that the force F y and the moment M z depend on the diagonal term of the added-mass matrix and its corresponding DOF (respectively, f and ̈𝛼 ) and also on an extra- diagonal term associated to the other DOF (respectively, ̈𝛼 and f ). Therefore the expression (5) can be simplified under the form: with (11)

f k+1 = rf tmp +(1 -r) f k , (12) r = 1 1 + e c 2 MA 22 L e K lin . ( 13
) k+1 = K lin K ang f k+1 . ( 14
) f k+1 = B M z F y k with B = ab cd , (15) B -1 f k+1 = M z F y k , (16)
B -1 = 1 𝛥 B d -b -ca and 𝛥 B = det (B)=ad-bc< 0. ( 17
) M z F y = M z F y -C f ̈𝛼 ,
In a similar way than the simplified model, we add the added-mass terms on both sides of the matrix system (15) and we use the BDF2 to transform the increment of acceleration into an increment of position, through the e c 2 coefficient. This gives the following expression: with Finally, we can retrieve a more interesting expression with a left multiplication both sides by D -1 :

(18) C = MA 62 MA 66 MA 22 MA 26 . ( 19
) D f k+1 = M z F y k + e c 2 C f k , ( 20
) D = B -1 + e c 2 C . ( 21
) f k+1 = D -1 M z F y k + e c 2 D -1 C f k .
This expression is not developed because it would unnecessarily take space. It should be noted that, even with developments, a simple expression in form of ( 11) cannot be found, where a relaxation coefficient would appear. Nevertheless the procedure acts similarly as a relaxation operator.

Results

Figure 6 shows six snapshots of the flow the rowing stroke 154, obtained by a numerical simulation performed with a coarse grid of 1.2. ⋅ 10 6 cells. The characteristics of a real rowing strokes are well reproduced. In particular, one can observe the breaking of the free surface and the creation of an air cavity behind the blade. Two grid convergence studies were carried out with the three levels of mesh fineness of Table 2. They are based on two different experimental trials whose parameters are voluntarily chosen as different as possible in the available database: the smallest velocity against the largest simultaneously Fig. 6 Visualisation, from a stationary location, looking in the hull-to-bank direction, of a simulation of a simplified rowing stroke. The free surface is coloured according to its vertical coordinate z (in m) with the largest immersion angle against the smallest (trials 67 and 154). As the results are similar, only the study of the trial 154 is presented. Figure 7 shows the evolutions of the hydrodynamic forces F x , F y and F z and the moment M z rela- tive to the axis of rotation, in function of the sweep angle 0 . The global shapes of the force curves of the experimental trials are well captured by the simulations. The propulsive force and moment have similar profiles and the same observations can be made. Numerical simulations overestimate from 5 to 12% the maxima. The error between experiments and simulations grows in the first half of the rowing stroke, till 0 ≈ 90 • , then it decreases up to 0 ≈ 130

• and it strongly increases until the end. The maxima of the force F y are also overestimated, by the same order of magnitude. However, the profile is disturbed in the interval [70

• , 90 • ]
. Even if it is the least well fitted, the shape of the F z profile is globally captured.

All the results obtained with the coarse, medium and fine meshes are really close, as visible in Fig. 7. The discretisation error of the coarse mesh is sufficiently limited and its maximal value occurs for 0 ≈ 90

• for the propulsive force F x and moment M z . For this particular angle, the decrease of the discretisation error between the coarse mesh and the fine mesh is only of about 2% (relatively to the experimental value). So the coarse mesh is dense enough and has been chosen for the other computations later described. No study of time step convergence has been done because all simulations were carried out with a Courant number C o < 0.3 for the area surrounding the free surface, allowing the use the compressive property of the BRICS scheme in order to capture the free surface. The total CPU time needed for those computations is between about 1000 h (on 56 cores) for the coarse grids and 10000 h (on 224 cores) for the fine grids.

Using the quite rich database of results from the experimental campaign, two sets of rowing strokes were selected where only one parameter varies in each of them: the immersion and roll angles. The objective is to analyse the influence of those parameters and to check if the trends obtained experimentally are similar to those of the numerical simulations.

Figure 8 shows the influence of the immersion angle on the force F z and moment M z . F z has been chosen rather than the propulsive force F x or the transversal force F y because these curves show similar trends compared to the one of M z , so the comments would be identical. Curves of the same marker shape represent the same immersion angle. Those having filled markers associated with continuous lines are for the numerical signals and those with empty markers and Fig. 7 Study of the grid convergence based on the angular evolution of the hydrodynamic forces F x , F y and F z and the moment M z relative to the axis of rotation dashed lines refer to the experiment data. Again, differences between experimental and numerical data are observed, especially for the low magnitude component F z where the relative error is larger. However, the influence of 1 follows the same law in both cases: the maxima of the forces and moments are increasing functions of the immersion angle, for the tested values. With Fig. 9, the same remark can be drawn regarding the roll angle but only for the force F z . The propulsive moment M z has an opposite behaviour. Again the experimental trends are well captured by the simulations, even for the low magnitude vertical component of the force.

The influence of the models for the oar bending can be seen in Fig. 10, where the propulsive force F x is plotted along the sweep angle for the simplified and the complete models. As expected the complete model gives a smaller error. At 0 = 98

• , the error indeed decreases from 14 to 12% . Whatever the model, it clearly appears that the oar flexibility has a considerable effect on the propulsive force. When the oar is rigid, the force responds earlier in time Fig. 8 Influence of the immersion angle on the angular evolution of the hydrodynamic force F z and moment M z Fig. 9 Influence of the roll angle on the angular evolution of the hydrodynamic force F z and moment M z Fig. 10 Influence the model of oar bending on the angular evolution of the hydrodynamic force F x (or in sweep angle) and undergoes more variations. The maximal values for the linear and angular deflections are, respectively, of about 12 cm and 8

• for the simulation with the complete model. They are high enough to explain those different behaviours.

Discussion

The results for this configuration in laboratory are less error prone than those obtained in real conditions [15], because there are more uncertainties on the latter. In the first configuration the propulsive force and moment are overestimated by 5-12% while they are underestimated by 15-20% in the real situation. This margin of error is satisfactory given that the flow is highly unsteady with a violent free surface deformation. To the authors' knowledge no other computations with this level of complexity have been published, making the comparison impossible. That is why we analyse the error in this discussion with a focus on its sources and its quantification.

Inventory of the possible sources of error

We propose a list of the main sources of error for both configurations: trials in laboratory and on site.

A first class of error is the experimental one. Unlike the in situ configuration where neither the roll of the boat nor the altitude of the blade were measured, here all the kinematics is known and available to be numerically reproduced. In addition to that, trials in laboratory benefits from a stable environment, which is not necessary the case of in situ measurements. However, even if uncertainties are more controlled in laboratory, they can affect the accuracy of derived quantities. Although the intrinsic error margins of each sensor is sufficiently small (lower than one percent) with a good enough repeatability, this can lead by propagation to higher uncertainties on some quantities of interest, as proved in 4.2. The accuracy and the repeatability of the measurements are here critical points due to the specific kinematic chain leading to the blade motion with respect to water. Regarding the in situ configuration, this issue is far more difficult to tackle since the velocity of the hull has large and abrupt variations. The velocity sensor itself and/or the associated low-pass filtering of its output signal may lead to a loss of frequential contents and a slight but still excessive smoothing. As a consequence, the reproduction of the real kinematics which is used to drive the simulation, especially the resulting acceleration may be altered and cause significant differences when comparing the fluid force, due to the fact that a large part of it depends on its acceleration with respect to water.

The numerical errors constitute a second class of error. First there is the modelling error, or the difference between the real flow and the mathematical model chosen to describe it. In this case the model for the resolution of the flow is made of the RANS equations, or the Euler model (null viscosity and slip condition on walls) for the real rowing stroke. Those models allow an accurate description of the physics, at least for the angular interval of interest (i.e. till the finish phase). However, those models fail at the end of this simplified stroke, i.e. for sweep angles such as 𝜓 0 > 130

• . Some additional investigations using more advanced turbulent models, such as DES (detached eddy simulation) or LES (large eddy simulation), could be useful, even if the cost would be high. Another modelling error should be examined, i.e. the boundary conditions. The effect of the vertical walls has been neglected but, as the blade comes near the lateral wall (about two chord lengths), it might be interesting to model the confinement effect. Then there is the discretisation error, that is to say the difference between the exact solution of the mathematical model and the one of the discretised equations implemented in the solver. Here we checked that the influence of the grid size is weaker than the observed difference with the experimental measurements. Finally, there is the convergence error, or the error between the exact solution of the discretised equations and the obtained solution after a finite number of iterations or when a given order of convergence is reached. In our computations, the convergence of the forces within each time step is checked a posteriori.

Study of the uncertainties on the forces

The hydrodynamic forces (or moments) applied on the blade depend on the incident velocity of the flow, which should be expressed at a point where the moment is null. As this characteristics cannot be always verified, one often chooses to express the velocity at an arbitrary point, approximately at the middle of the blade. In this study, this point is I (identical to E, see Fig. 11) and distant of L I = 1.815 m of the oar- lock. For the sake of simplicity only one DOF is considered for the boat and one for the oar. Furthermore the flexibility of the oar is neglected. The boat moves at a constant velocity V b = 4.202m ⋅ s -1 along the (O, ) axis. The oar rotates relative to the boat at the point B along the (B, 0 ) axis, at an angular velocity ̇𝜓 0 given by the law expressed in Fig. 3, and is inclined on the water of an angle 1 = 9.5

• . So this configuration is similar to the trial 154.

The incident velocity i of the flow relative to the blade is the opposite of the velocity of the point I belonging to the blade relatively to the water ∈ ∕ and is obtained by where ∈ ∕ is simply the velocity of the boat expressed in ℜ 0 , = V b , and ∈ ∕ = L I cos 𝜃 1 ̇𝜓 0 is expressed in the frame ℜ 1 associated to the oar (see Fig. 11). The incident velocity, expressed in ℜ 1 is then equal to With a quantity called kinematic efficiency (see Barré and Kobus [START_REF] Barré | Comparison between common models of forces on blades and forces measured by towing tank tests[END_REF]) and defined by norm of the incident velocity can be expressed by and the incident angle, defined by the angle between (equivalent to and in this planar motion) and (see Fig. 11), is equal to

(22) i =-∈ ∕ =- ∈ ∕ + ∈ ∕ , (23)
V i = ⎛ ⎜ ⎜ ⎝ -V b cos 𝜓 0 -L I cos 𝜃 1 ̇𝜓 0 + V b sin 𝜓 0 0 ⎞ ⎟ ⎟ ⎠ ℜ 1 . (24) 𝜂 = V b sin 𝜓 0 L I cos 𝜃 1 ̇𝜓 0 , (
) V i = V b √ 1 + 1 1 -2 sin 2 0 , 25
) i = arctan tan 0 1 -1 . (26
Let be a quantity depending on N parameters k , k ∈ [ [ 1, N]] , each of them estimated with an absolute error k , 𝛥 k > 0 . The absolute error on can be estimated by the formula The parameters of V i and i are 0 , ̇𝜓 0 , 1 , V b and L I . To sim- plify the developments, we consider 0 and ̇𝜓 0 as independ- ent. Moreover, we only select the parameters V b , ̇𝜓 0 and L I for this study because a posteriori results show minor influences of 1 and 0 and it makes the formulation (28) shorter. We apply (27) for V i and i. After some developments, the corresponding relative (for V i ) and absolute (for i) errors can be expressed as and For 0 = 90 • , we have i = 90

• and i = 0

• (see equations ( 26) and ( 29)). This means that the drag and propulsive forces are equal for this particular orientation. The quantification of the uncertainty on the propulsive force will only be given for this case because it is the more interesting one and a general approach would be more intricate. Using dimensional analysis, V i (considered as constant for this specific case) and F D are linked by where is the water density, A the projected area of the blade and C D the drag coefficient. These three quantities are con- sidered constant. As a result, the relative uncertainty on the propulsive force F x follows the following inequality As the goal is to quantify the uncertainties, values have to be given for the intrinsic uncertainties of the parameters. Since the real values are not known, we have chosen them in a quite fair manner. Indeed they are all three of about half a percent, which is supposed to give confidence for an

(27) ≤ N ∑ k=1 | | | | k | | | | k . ( 28
) 𝛥V i V i ≤ V b V i 2 | | | | | 1 - sin 2 𝜓 0 𝜂 | | | | | 𝛥V b V b + 1 𝜂 | | | | 1 𝜂 -1 | | | | sin 2 𝜓 0 𝛥L I L I + 𝛥 ̇𝜓 0 ̇𝜓 0 (29) 𝛥i ≤ | | tan 𝜓 0 | | 𝜂 1 + tan 2 i 𝛥V b V b + 𝛥L I L I + 𝛥 ̇𝜓 0 ̇𝜓 0 . ( 30
) C D = F D 1 2 AV i 2 , ( 31
) for 0 = 90 • , F x F x ≤ 2 V i V i .
evaluation of the quantity of interest. The values are documented on Table 3. Using them we can plot the evolution of V i and

V i V i
relative to the sweep angle 0 , as done in Fig. 12.

V i is maximal and equal to V b at both extremities of the stroke. Between them, V i monotonically decreases, reaches its minimal value around 0 ≈ 90 • (i.e. when the propulsive force and moment are maximal) and increases till the end. The relative error has an inverted behaviour and it is therefore maximal when V i is minimal. The corresponding value is of about 5.5%, indicating a relative error on the forces of 11%. This value is close to the gap obtained between the experimental and numerical results (i.e. 5-12%). In fact the uncertainty is higher because there should be an additional term on the right-hand side of (31), coming from the unsteady nature of C D in this kind of flow which is not taken into account in equation (30). It demonstrates that our numerical simulations are robust enough and that the step of validation of our frame of modelling has been completed.

Figure 13 allows to understand the proportion of each parameter in the evolution of the relative error on V i and likewise for the forces. The curves are plotted using the right-hand side of the inequality (28). The parameter V b has a dominant influence at both extremities of the stroke while L I and ̇𝜓 0 prevail more in the central part and behave similarly. The three parameters have equal influences for 0 ≈ 45 • and 0 ≈ 115

• . At 0 = 90 • where the propulsive and moment are maximal, the parameters L I and ̇𝜓 0 account each for 38% of V i , while V b is responsible for 24% of this uncer- tainty. At 0 = 59

• and 0 = 113 • , V b has no influence on V i . It is explained by the fact that i is perpendicular to b for those two orientations and, at the first order, a variation of V b does not change V i but only the angle i.

Conclusions

The frame of modelling and numerical simulation of a rowing stroke was tested and validated using an experimental database coming from trials performed in a towing-tank, whose measurements are less prone to uncertainties than those of experiments carried out in situ [15]. The counterpart of this configuration lies in the simplification of the rowing stroke, where two DOFs are lost for the propulsive phase: the rotation allowing the immersion of the oar on the water and the vertical translation. In addition the carriage of the towing tank moves at a constant velocity while a rowing boat has an unsteady kinematics.

In this more controlled configuration better results were obtained. The simulations are able to capture the profiles of the forces and to reproduce the trend induced by the variations of two parameters: the immersion and roll angles of the oar. They overestimate the maximum of the propulsive force F x and moment M z of 5-12% with the three levels of mesh fineness tested, meaning that the 1.2-million-cell mesh is sufficient for those simulations. The errors between experimental and numerical results are close to the uncertainty of 11% obtained in a study with only three parameters whose intrinsic uncertainty were only of about 0.5% for each. On the one hand, this shows that the simulation framework is a robust tool for the simulation of a rowing stroke. On the other hand, the study of the uncertainties indicates that this kind of experiments is a very delicate task, even for a priori very controlled and stable environments.

A new model for the oar bending taking separately into account the hydrodynamic force and moment acting on the blade has been added to the existing procedure for the FSI coupling, allowing a more accurate evaluation of the linear and angular deflections. As these two DOFs are independent in this model, it may lead to oscillations followed by a divergence of the simulations. To avoid this situation, a method to stabilise the coupling has been developed, whose principle is to couple the relaxations of those variables including the off-diagonal terms of the added mass matrix.

The next step of this ongoing work is a co-simulation including the resolutions of both flows, the one around the blade and the one around the hull, and the resolution of the dynamics of the multibody boat-oars-rower system. The latter will be inspired of the work of Rongère et al. [START_REF] Rongère | Dynamic modeling and simulation of rowing with a robotics formalism[END_REF]. This frame will constitute the ambitious but necessary basis of a fine modelling for the analysis and the optimisation of the performance, like the influence of parameters.

Fig. 1 Fig. 2 Fig. 3

 123 Fig.1Sketch of the experimental configuration used to generate a simplified rowing stroke in a towing tank with a real oar blade (nonreal dimensions)

Fig. 4

 4 Fig. 4 Views of the computational domain and a 1.2 ⋅ 10 6 -cell grid used for the numerical simulations of a simplified rowing stroke. a Computational domain. b Sectional view at X 0 = 0.3 m. c Surface mesh of the Hatchet-kind blade (external face)

Fig. 11

 11 Fig. 11 Definition of the incident velocity on the oar blade

Fig. 13

 13 Fig. 13 Evolution of the contributions of the parameters ̇𝜓 0 , L I and V b on the absolute error of the incident velocity in function of the sweep angle

Table 1

 1 List of the experimental trials used and their parameters

	# Trial	V c (ms -1 )	c 2 ( • )	1 ( • )
	59	4.0	0	10.5
	61	4.0	4	10.5
	63	4.0	10	10.5
	67	4.0	8	11.5
	99	4.1	4	11.5
	109	4.1	4	9.5
	122	4.1	4	10.5
	154	4.2	6	9.5

Table 2

 2 Number of cells (expressed in millions) of the grid used for the simulations

	Grids		Coarse	Medium	Fine
	Trials	67	1.2	2.7	12.3
		154	1.2	2.7	11.4
		Others	0.8-1.3	-	-

  MA 11 MA 12 MA 13 MA 14 MA 15 MA 16 MA 21 MA 22 MA 23 MA 24 MA 25 MA 26 MA 31 MA 32 MA 33 MA 34 MA 35 MA 36 MA 41 MA 42 MA 43 MA 44 MA 45 MA 46 MA 51 MA 52 MA 53 MA 54 MA 55 MA 56 MA 61 MA 62 MA 63 MA 64 MA 65 MA 66

Table 3

 3 Parameters chosen for the uncertainty study, with their nominal values and their absolute and relative uncertainties † For the maximum value of ̇𝜓 0 Evolution of the norm of the incident velocity V i and its relative error

	Parameter	Value	Absolute uncertainty	Relative
				uncer-
				tainty
	V b	4.202m ⋅ s -1	2cm ⋅ s -1	0.48%
	̇𝜓 0	180.5 deg/s †	1 deg/s	0.55% †
	L I	1.815 m	1 cm	0.55%

The towing tank has a length of 140 m today.

The acronym METHRIC comes from the French for Modelling of High-Reynolds Incompressible Turbulent Flows and Couplings

If the trial starts with an oar set to an angle 0 = 0 • , it may generate large oscillations and disturb the flow. To avoid this behaviour, the oar is set to a slightly different angle, such that the lift force is null.

Acknowledgements This work was granted access to the HPC resources of GENCI (Grand Equipement National de Calcul Intensif) under the allocation A0022A00129 and was supported by a grant from the Région Pays de la Loire through the project ANOPACy.


Christian Wagner
email: c.wagner@tum.de

Thomas Thümmel

Daniel Rixen

S Y Mpo

O Tatin

Experimental Prediciton of Instability in Rotor Seal Systems using Output Only Data

Keywords: Rotordynamics, stability diagnosis, rotor seal system, active magnetic bearing

In common high-speed turbomachinery or pumps, usually contactless seals are used to separate different fluids and pressure levels. e seal gap causes a fluid flow, which generates rotational speed dependent forces to the rotor in the deflectional and tangential directions. ey act like cross-coupled stiffness and can excite the rotor system to large vibrations, the instability occurs similar to the 'oil whip' phenomenon in journal bearings. us, a diagnosis method is necessary to ensure safe operation. To predict the speed limit of the machinery, two methods are presented in this research: an analysis using active magnetic bearing excitation and an output only damping ratio estimation method. e methods are experimentally evaluated on a test rig and using numerical simulations.

INTRODUCTION

Seals in turbines, centrifugal pumps or compressors are commonly used to separate different fluids and pressure levels. Because of the high rotational speeds of common turbomachines, contactless seals, such as floating ring, labyrinth or seal gaps are inserted between the rotating and the stationary parts. e clearance which always exists around these contactless seals permits fluid flow through the gap. For an eccentric rotor position, the fluid-velocity distribution inside the seal becomes asymmetric, which generates forces on the rotor.

ese can result in stiffening, restoring, and damping effects, and in particular tangential forces, caused by the whirling fluid flow inside the seal gap. ese are responsible for cross-coupling terms in the stiffness matrix.

ey cause a rotor vibration at the onset speed of instability, leading to large displacements and to a breakdown of the system, similar to the 'oil-whip' phenomenon in journal bearings. Furthermore, the rotor-system or the seal condition can change during lifetime. us, an experimental diagnosis method is necessary to avoid rotor instability and to ensure safe and stable long-term operation.

A modern approach is to predict the stability limit of the entire system using output only run up measurements. Observing the parameters during normal operation, a change of the stability limit can be detected and used as an indicator for monitoring.

e first focus of this contribution is an introduction to rotor-seal system instability using a rotor model and rotordynamic seal coefficients, based on [ ]. e second focus is the comparison of several experimental methods on the test rig with and without additional excitation using an active magnetic bearing (AMB). is offers possibilities for actual turbo-machinery applications.

In the recent decades, a lot of effort has gone into modeling and identifying the dynamics of rotor seal systems. Amongst the first researchers to model the effects of seals using rotordynamic seal coefficients were . MODELING -

h s = m s 0 0 m s ẍ ÿ + c s 2m s γΩ -2m s γΩ c s ẋ ẏ + k s -m s γ 2 Ω 2 γΩc s -γΩc s k s -m s γ 2 Ω 2 x y ( )
where the rotational speed of the rotor is Ω. .

Dynamic Behavior

To analyze the dynamic behavior of the rotor seal system, alongside a time integration simulation, the eigenvalues are also calculated. erefore we write for the coupled system (rotor + two seals), as in [ ]:

M 0 0 M ẍ ÿ + Cc -cC ẋ ẏ + Kk -kK x y = h ( )
Where for improved readability the summarized coefficients are:

M = m r + 2m s C = 2c s ; c = 4m s γΩ K = k r + 2(k s -m s γ 2 Ω 2 ); k = 2γΩc s
It can be seen that the cross coupled terms k, c lead to tangential forces, which can destabilize the system. ese forces transmit energy from the rotor rotation to the bending motion. Assuming a symmetrical rotor system, we substitute x and y with the complex coordinates z = x + j y and F z = h x + jh y , where j is the imaginary unit. An eigenvalue analysis, for determining the system's stability and natural frequencies is common; F z = 0 and z = ẑe λt are used:

M λ 2 + Cλ + K -j (cλ + k) = 0

. Simulation results

To present the influence of the seals on the rotor system, the simplified J model is used with the parameters corresponding to our test rig: ) shows the natural frequency, the critical speed, of the 'dry' (seal-less) rotor. It can be seen that the first natural frequency is constant, which leads to a critical speed at 2316 rpm. e Campbell diagram (see fig. ) shows the behavior of the natural frequency influenced by the seals. It follows that the first natural frequency in forward mode ω 1 f is lower than half the rotational speed, the backward mode ω 1b is less and can also be seen here. Figure shows the damping ratio D = δ/ω o of the rotor seal system, ω 0 is the undamped natural frequency. e damping ratio D has a zero crossing at 8222 rpm, above which the system becomes unstable. At the onset speed of instability, the rotor starts to vibrate at its first natural frequency, which dominates the displacement response (see waterfall plot of time simulation with unbalance and random noise in fig. ). is presented, then two experimental methods for predicting the onset speed of instability are shown. e first method uses an active magnetic bearing (AMB) for excitation (as a reference method), the second one uses displacement sensor output only data to determine the system's damping ratio.

. EXPERIMENTAL INVESTIGATIONS

. Seals Test Rig e stability analysis method applied to rotor seal systems is clearly described in [ ]. Here it is used as a reference method to compare the output only approach.

. . Analysis Principle

To apply this analysis method, an AMB to excite the rotor seal system with a force F z = h x + jh y in the forward whirl direction (same direction as the rotor rotation Ω) is used. e measured response is the rotor displacement signal z = x+ j y.

Transforming the excitation F z and the response z into the frequency domain (using F z = Fz e jωt and z = ẑe jωt ), the transfer function G(ω) can be calculated as follows:

G(ω) = ẑ Fz = 1 -Mω 2 + cω + K + j (Cω -k) ( )
e system becomes unstable ( ẑ arises), when the denominator of eq. ( ) nullifies. is is equivalent to the zero crossing of the damping ratio D, shown above. erefore the transfer function G is separated into real and imaginary part as follows:

Re{G(ω)} = -Mω 2 + cω + K (-Mω 2 + cω + K ) 2 -(Cω -k) 2 ( ) Im{G(ω)} = -Cω + k (-Mω 2 + cω + K ) 2 -(Cω -k) 2 ( )
In this form, it is possible to determine the zero crossing frequencies ω re0 and ω im0 :

0 = -Mω 2 re0 + cω re0 + K ( ) 0 = -Cω im0 + k ( )
Solved for ω re0 and ω im0 , inserting M 's model:

ω re0 = γ 1 + m r /2m s Ω ± k r + 2k s m r + 2m s + ( 4m 2 s γ 2 (m r + 2m s ) 2 - γ 2 2 + m r /m s )Ω 2 ( ) ω im0 = γΩ ( )
e structure of the equations shows that for low rotational speeds, ω re0 is almost equal to the undamped natural frequency ω 0 and exhibits linear behavior at high rotational speeds Ω. ω im0 is linearly proportional to the rotational speed.

. . Experimental Results

e analysis on the test rig is performed at several constant rotational speeds Ω (in the save operating range). en the rotor is excited by the AMB with F z in forward whirl direction using a sweep (cosine in x, sine in y) about the frequency range of interest (from 10 Hz up to 100 Hz) for calculating the transfer function G(ω) using fast Fourier transformation ( ).

e zero crossing frequencies of the real and imaginary parts of G(ω), ω re0 and ω im0 must be found (see fig.

). ese zero crossing frequencies show almost linear be- havior dependent on the rotational speed Ω, as described above. When both curves cross, the instability occurs. Using a simple linear fit (least squares method) to the points and an extrapolation, the onset speed of instability can be predicted based on measurements in the safe operating range (see fig.

).

. Damping Ratio Estimation e damping ratio estimation using output only data is based on the investigations of [ ] and [ ]. To remove the dominating unbalance response of the rotor system, an interpolation method, presented by [ ] for an operational modal analysis case, is applied.

. . Theoretical Methodology

If the damping ratio becomes zero, the rotor seal system becomes unstable and the amplitudes increase rapidly. e idea of this output only method (without AMB excitation) is to fit the shape of the natural frequency amplitude (see the waterfall plot in fig. ) and to calculate the damping ratio of the system. erefore a stationary rotor run-up is performed. e displacement signals x and y, transformed to z = x + j y are used for the analysis. e advantage of this transformation is its suppression of the backward whirl motion in the positive spectrum to focus on the forward one (the system's instability is a vibration in the forward whirl motion). To perform curve fi ing, the dominating unbalance response at the frequency of the rotational speed Ω, must be removed. For this an interpolation method proposed by J [ ] is used. It is possible to fit a PT-model (MCK model) curve to the measured amplitude data:

|G(ω)| = a (1 -ω ω 1 2 ) 2 + (2D ω ω 1 ) 2 ( )
using the first natural frequency ω 1 and the scaling factor a. .

DISCUSSION AND COMPARISON OF THE APPLIED METHODOLOGIES

In this article, two experimental methods for predicting the onset speed of instability are presented. e predicted rotational speeds for both experimental methods and the bulk flow simulation results are compared in the table:

Name

Co-ad Output-Only Simulation Predicted speed in rpm 8609 8015 8223

Comparing the waterfall plots of the simulation (just 222 rpm before the instability occurs) and the measurement, fig.

and fig. , the trend in the increase in the amplitude at the system's first natural frequency is similar. is implies that at 7000 rpm the test rig is very close to the instability (the peak of the first natural frequency is actually higher than the unbalance response). is behavior shows that the estimated onset speeds of instability are too high compared with the real system. To this extent the presented Output-Only method gives a be er prediction of the stability limit than the Co-ad method.

Additional comments by the authors: e early stage of the presented methodologies and measurements imply that we cannot make a comment on the uncertainties of the results. Future work will be a sensitivity analysis using different fluid parameters, seal geometries etc. To improve the prediction methodology, the amplitude of mass unbalance should taken into account for the run-up measurement. A higher rotor amplitude leads (due to nonlinear effects in the seal gap, see [ ]) to a more stable system and can influence the prediction. e ability of the MCK model to provide an accurate fit for calculating the damping ratio is well in our simple rotor test rig case, where the first bending mode is the dominating motion. In a real turbomachinery case with multiple resonance peaks the user has to pick out the 'critical' natural frequencies (a coupled structure-seal simulation can be used to identify the 'unstable' natural frequencies to observe). To include the output only estimation in a condition monitoring system, only rotor displacement signals have to be observed.

e damping ratio and its rotational speed dependency can be estimated during run-up and operation. Observing the available stability margin can be used as a characteristic monitoring value.

. CONCLUSION e presented experimental methods for characterizing the stability limit of a rotor seal system are in the stable and safe operating range.

e first method uses an AMB to excite the system at several discrete rotational speeds for measuring the transfer function of rotor displacement to AMB force. An analysis of real and imaginary part according to the rotational speed leads to an estimation of the onset speed of instability.

e second presented method uses Output-Only data (rotor displacement). A PT-model fit is used to identify the system's damping ratio. e rotational speed dependency of the damping ratio leads to an estimation of the stability limit, which can be used as a monitoring value.

e methodologies are examined to a numerical model and to the test rig.

ACKNOWLEDGMENTS

is project is supported by the Ludwig Bölkow Campus and the Bavarian State. e friendly and effective cooperation between the partners of the research project is much appreciated.

  B [ ] and C [ ], who solve the bulk-flow equations and create analytical short seal solutions. M [ ],[ ] improves the seal models and includes the eccentricity influence. Other experimental investigations using AMB's in rotorseal systems are [ ], [ ] and [ ]. Using a levitating rotor to examine the seal is performed by [ ] and [ ]. e idea of using output only data to identify the system's damping ratio is also used in different ways by [ ], [ ] and [ ]. A lot of investigations take place in the field of operational modal analysis. An interesting method to remove the unbalance response is proposed by [ ]. e fundamentals of rotordynamics, modeling and simulation are clearly described in [ ], [ ] and [ ].

  Figure . Rotor model

  e relation of m s , c s , k s ,γto the seal geometry, fluid properties and friction factors is well described in [ ]. Note also that c s , k s ,γ are nonlinear functions of x and y, see [ ].

  Figure . 'Dry' rotor Campbell diagram

  Figure . Rotor seal Campbell diagram

Figure . 1 Figure .

 .1. Figure .Damping ratio of the rotor seal system

  Figure . Seals test rig

Figure .

 . Figure .Co-ad analysis of simulation data

Figure .

 . Figure . Co-ad analysis of measured data

Figure .

 . Figure . Waterfall plot of measured amplitude data of the test rig

Figure . Figure .

 .. Figure . Amplitude spectrum of the measured rotor seal system; Raw data, interpolation to remove the unbalance response and PT-model fit

  Identification of Rotordynamic Seal Coefficients by Means of Impedance Matrix and an Optimization Strategy. In Proceedings of the XVII International Symposium on Dynamic Problems of Mechanics, São Sebastião, SP, Brazil, . ABCM. [ ] Zachary Sco Zutavern. Identification of Rotordynamic Forces in a Flexible Rotor System Using Magnetic Bearings Identification of Rotordynamic Forces in a Flexible Rotor System Using Magnetic Bearings. PhD thesis, Texas AM, . [ ] I. F. Santos and Peter Kjaer Svendsen. Noninvasive Parameter Identification in Rotordynamics Via Fluid Film Bearings-Linking Active Lubrication and Operational Modal Analysis. Journal of Engineering for Gas Turbines and Power, ( ): , . [ ] W. Tsunoda, C. Wagner, T. Berninger, T. ümmel, and D. Rixen. Measurement method of damping ratio as stability diagnosis for rotor-seal system with fluidic random excitation. In e th International Education Forum on Environment and Energy Science, San Diego, California, USA, . [ ] C. Hunter Cloud, Eric H. Maslen, and Lloyd E. Barre . Robert Gasch, Rainer Nordmann, and Herbert Pfützner. Rotordynamik. Springer Berlin Heidelberg, Berlin, Heidelberg, . [ ] Agnes Muszynska. Rotordynamics. CRC Press ( May ), .

					[ ] A Muszynska. Improvements in Lightly Loaded Rotor /
					Bearing and Rotor / Seal Models. Journal of Vibration,
					Acoustics, Stress, and Reliability in Design,	: -,
					.
					[ ] Klaus Kwanka. Dynamic Coefficients of Stepped
					Labyrinth Gas Seals. Journal of Engineering for Gas
					Turbines and Power,	: -,	.
					[ ] Manuel Gaszner. Rotordynamische Charakterisierung von
					Dichtungssystemen zur Anwendung in Kra werksdampf-
					turbinen. PhD thesis, Technical University of Munich,
					.
					[ ] J .G. M. David, G. R.	iago, and A. N. C. P. Fer-
					nando. Damping Ratio Estimation Techniques for Rotordynamic
					Stability Measurements. Journal of Engineering for Gas
					Turbines and Power,	:	--	-,	.
					[ ] Leonardo Baldassarre, Alberto Guglielmo, Michael
					Catanzaro, Raphael Timbó Silva, and Marcelo Accorsi
					Miranda. Operational Modal Analysis Application for
					the Measure of Logarithmic Decrement in Centrifugal
					Compressor. In th Turbomachinery Symposium (	),
	[ ] Agnieszka Muszynska. Whirl and whip-rotor/bearing	Houston, Texas,	.
	stability problems. Journal of Sound and vibration,	[ ] Joel Johansson, Henrik Samuelsson, Niels Jacobsen, and
	( ): -,	.			Anders Angantyr. Operational Modal Analysis of Large
	[ ] H.F. Black and D.N. Jenssen. Black, Jenssen Dynamic Hybrid Bearing Characteristics.pdf. In Proceed---ings of the Institution of Mechanical Engineers, Conference	-pole Rotating Machinery. In rd International Op-erational Modal Analysis Conference, pages -, Portonovo, AN, Italy, .
	Proceedings, Vol -, pages -, London, England,	[ ]
	. SAGE Publications Sage UK.
	[ ] Dara W. Childs. Dynamic Analysis of Turbulent Annu-
	lar Seals Based On Hirs' Lubrication Equation. ASME,
	Lubrication Technology,	: -,	.


Bernd Wagner
email: bernd.wagner@dlr.de

Julian D Pauw

Lucrezia Veggi

S Y Mpo

O Tatin

Assessment of Performance Variation on the Axial and Radial Forces in Turbopump Configurations for Liquid Rocket Engines

Keywords:

In the frame of the KonRAT project, selected topics of turbopump design and operation are considered. In the work package of secondary the following topics are considered systems gap flow, bearing, and seals. Especially the bearings of liquid rocket engines turbopumps are highly loaded due to harsh environment in which they are operating.

us, accurate knowledge of the acting radial and axial forces in the nominal operation is necessary for the design of the bearing unit and is subsequently important for safe operation of the turbomachinery.

is requires, in turn, an understanding of the flow in the rotor stator cavities, since a pressure distribution and thus forces can be derived therefrom. A program was developed for the calculation of the acting forces. It is generally applicable for different turbopump configurations. Random numbers are used to vary the input values. e results are presented.

INTRODUCTION

In modern liquid propellant rocket engines (LRE) turbopumps (TP) are inevitable to provide the combustion chamber with a high propellant massflow at high pressure. Both are necessary for the performance of the engine and an efficient, lightweight stage e.g. due to a low tank pressure, thus a reduced structure weight. In order to meet the engine requirements in terms of performance the turbopump assembly (TPA) must be able to operate at high speed, work over a wide operational range, and has a high power density. Depending on the configuration and the possible engine cycle, a TPA consists of inducers, pumps, turbines, a seals, bearings, and sha system, and a housing with flanges. Each part is designed to operate under harsh boundary conditions and with a low system weight. In order to minimize the stresses on heavily loaded parts like the bearings, the forces acting in radial and axial direction should be estimated. Furthermore, during the design process of the TP parts, the balance of the axial forces can be achieved.

erefore, the aerodynamic and hydraulic forces acting on each part of the TPA have to be calculated and brought into relation. e results are an important starting point for further design details and first rotordynamic analysis.

. PROJECT BACKGROUND e project KonRAT has been initialized 2015 as part of the national effort to build up further competencies in turbomachinery technology for LRE. KonRAT is a German abbreviation for application of rocket propulsion components for aerospace transport systems. Key objectives of this project are establishing design procedures for TP for cryogenic oxygen, investigations on selected problems in turbopumps and, finally, additive manufacturing methods for aerospace applications like turbopump parts and rocket engine valves. e turbopump branch of the project is mainly performed at the Technical University of Munich, for more information see [START_REF] Pauw | Untersuchungen zu Turbopumpen an der Technischen Universität München im Rahmen des KonRAT Projekts[END_REF].

e work on TP includes state of the art design procedures [START_REF] Veggi | A study on the design of LOX turbopump inducers[END_REF][START_REF] Souverein | On the effect of axial turbine rotor blade design on efficiency: a parametric study of the Baljé-diagram[END_REF]. In a newly erected test bench measurements with heated water are performed for TP components [START_REF] Pauw | Design procedure of a turbopump test bench[END_REF]. Experimental and numerical observations on seals and bearing and their influence on the rotordynamic are also ongoing [START_REF] Ch | Rotordynamic effects in turbopumps for space propulsion systems -first minimal models and experimental validation[END_REF][START_REF] Ch | Rotordynamik: Modellierung und einfluss von schrägkugellagern für hochdrehzahlanwendungen[END_REF]. Detailed analysis of seal and bearing elastohydrodynamic are performed with sophisticated fluid structure interaction solvers.

e DLR contributes to the project with observations on secondary flows, including the axial and radial forces resulting from these flows.

Two baseline configurations are considered in the project. First, a future upper stage in the 120 kN class for launchers is investigated, second, a main stage engine for sub-orbital flight vehicle. e configurations are summarized in table 1. .

CALCULATION OF ACTING FORCES

e major part of the forces on the bearings is caused by the pressure forces acting on the outer walls of flow surrounded rotating parts. In addition, forces and moments come from momentum and weight, as well as vibrations, unbalance and disturbances in the flow. e used methods to calculate these forces are briefly presented here.

. Axial Thrust

First, reference is made to the axial loads, also called axial thrust or force, which is produced by the inducers, the impellers and the turbines. Simple ways to estimate these forces are shown here for no flow through configurations.

e axial forces have two causes. First, forces are caused by the deflection of the flow as an impulse force on the side walls of the rotating components, which can generally be formulated

dF = d dt (dm • c). (1)
Second, pressure forces act on the side walls, which are described as follows for axis-symmetric cases F = 2π prdr.

(

) 2
From this, the axial forces can be derived from the pressure distribution of surfaces of all the individual components.

. . Inducer Typically, inducers are installed in liquid rocket turbopumps, which have a pure axial construction without a shroud. ey are therefore to be treated as axial rotors, for which the axial force can be determined by

F ax,in = π 4 (d 2 2 -d 2 n )ρgH in (3)
In Gülich [START_REF] Gülich | Centrifugal Pumps[END_REF] this equation is additionally multiplied by a correction factor between 1 and 1.1. In case the inducers are equipped with a variable diameter an impulse change of the flow and the resulting force on the inducer can be considered. However, this additional force is very low compared to the other forces and is ignored here.

. . Impeller

e axial forces acting on the impeller are again composed by impulse forces and pressure forces acting on the outer walls of the impeller. Here, only the case of a wheel side wall which is not flowed through is to be considered. is is an acceptable simplification for the design considerations discussed here. According to Gülich [START_REF] Gülich | Centrifugal Pumps[END_REF] from equation 2 the forces on the hub and shroud of the impellers can be determined in the axial direction with

F = πr 2 2 1 -x 2 ∆p imp - ρ 4 u 2 2 k2 1 -x 2 2 (4)
Where x is the respective ratio between the inner and outer diameters of the impeller side wall, and k is the average rotational factor. e geometrical influences on the pressure profile acting on the side wall are mainly introduced by this factor. Since in the present case exclusively no flow through is considered in the rotor stator cavities, the rotation factor can be expressed by

k 0 = 1 1 + A B * c f,w c f,imp with A = 1 cos δ w r w r 2 5 - r i r 2 5 + 5 t ax r 2 r w r 2 4 B = 1 cos δ R 1 - r i r 2 5 + 5 t ax r 2 r i r 2 4 * 1 + r w -r i t ax tan δ w - r 2 -r i t ax tan δ R (5)
If the angles α between the side walls of the impeller and the housing side walls are neglected and r i is small compared to the impeller radius r 2 , equation 5 is simplified to

k 0 = 1 1 + rw r2 2 rw r2 + 5 tax r2 c f,w c f,imp (6)
e impulse force is also very small for the impeller and is calculated by

F imp = ρQc 1m . (7)
Depending on the installation situation, the static pressure at the impeller inlet and the unbalanced sha thrust over the entire turbopump can also be considered. Axial forces resulting from small pressure affected surfaces and seals are neglected in this preliminary study.

. . Turbine

For the calculation of the forces acting on the turbine a general approach is applied in order to meet the possible configurations of the turbines. is includes the two causes of the acting forces from equations ( 1) and ( 2) and is formulated as

F ax,t = πα part 360 • (p 1 -p 2 ) + (c 2 1,ax ρ 1 -c 2 2,ax ρ 2 ) . (8)
A partial admission inflow of the turbine can be taken into account by the angle α part .

is would cause an additional moment on the sha . Equation 8, however, represents only an estimation of the forces, a more precise determination of axial forces acting on the turbine disks should be carried out with numerical flow solvers like for example in Dorney [START_REF] Dorney | Calculation of turbine axial thrust by coupled CFD simulations of the main flow path and secondary cavity flow in an SLI LOx turbine[END_REF].

. Radial Forces

Since the flow in turbopumps is axis-symmetrical, most of the forces acting in the radial direction are mutually compensated. Radial forces are only caused by a disturbance of rotational symmetry. In the following, a difference is made between stationary forces, whose effective direction is determined by the geometry of the turbopump, and dynamic forces.

e major part of the radial loads acting on an impeller is usually produced by the spiral-shaped outlet at the last stage. It can be calculated as follows:

F R = k R ρHd 2 b 2,ges (9)
Where b 2,ges is the width of the impeller at the outlet, including the hub and shroud wall thickness. e radial thrust coefficient k R can be estimated according to Stepanoff [START_REF] Stepanoff | Radial-und Axialpumpen. eorie, Entwurf, Anwendung[END_REF] with

k R = 0.36 1 - Q Q * . ( 10
)
In this case the radial force acts in the direction of a point which lies approximately 60 • downstream of the housing tongue when partial load is reached.

e maximum value of k R was used for the calculation in the program. More detailed methods for the calculation of the radial forces in spiral housings like e.g. from Domm and Hergt [START_REF] Domm | Radial forces on impeller of volute casing pumps[END_REF] are not yet taken into account and will be implemented in the future.

. Weight Forces -Gravitation

Even if the proportion of the weight forces at the acting axial and radial forces is quite small, it is calculated in the following form.

F g,x = -m(g + a) sin α w F g,y = -m(g + a) cos α w . ( 11
)
In addition to the gravitational acceleration g, the lioff acceleration a of the launcher or of the stage can also be added. e mass m is determined individually on the basis of the given dimensions and material for each component, but it can also be input directly if precise values are known from CAD models. Typically, turbopumps are installed at an angle of 90 • or 270 • to allow the pumps to be fed in the direction of acceleration. More rarely is a horizontal arrangement which means in the convention in here: 0 • or 180 • .

. Dynamic Forces -Unbalance e equations above consider only steady-state forces. During the operation of rocket engine turbopump strong dynamic loads also affect the rotor. Vibration and acoustics loads emerging from the engine and the rocket structure act from outside on the turbopump. e turbopump itself is exposed to hydraulic and aerodynamic vibrations which couples with rotordynamics and unbalance vibrations. e direction of action and/or the magnitude may change considerably over time. Most of these vibrations and forces are difficult to assess analytically in the preamble in a precise way. However, for forces caused by unbalance (mass asymmetry), an experimentally proven first formula for well-balanced, fast-rotating runners can be found in [START_REF] Sigloch | Strömungsmaschinen. Grundlagen und Anwendungen[END_REF] as

F un = (0.02...0.04)d 2 mn 2 . ( 12
)
0.04 was used in the program as prefactor, whereby the calculated unbalance forces are strongly overpredicted. e general application of this equation for the fast-rotating components in space propulsion drives is not entirely clear and variation of this factor has a linear impact on the prediction of dynamic forces due to unbalance.

. APPLIED METHOD . Program for Axial and Radial Forces

A tool has been developed to estimate the forces acting on inducers, impellers, and turbines in both axial and radial directions. e calculation of forces acting on the hydraulic parts in axial direction are based on Gülich [START_REF] Gülich | Centrifugal Pumps[END_REF], where the impulse, pressure forces acting on hub and shroud, and unbalanced axial force on sha are considered for no flow through at a first step. Radial forces are determined only for the volute with a simple method proposed by Stepanoff [START_REF] Stepanoff | Radial-und Axialpumpen. eorie, Entwurf, Anwendung[END_REF]. Aerodynamic forces on the turbine are estimated considering impulse forces and if present pressure difference. Additional moments due to partial admission of the turbine can be taken into account, too. Forces due to unbalance of the parts are roughly evaluated. Finally, the weight of each part is calculated depending on the orientation of the TP in the engine and the acceleration of the rocket stage. Based on the dimension and position of the parts, the resulting bearing forces are calculated.

Figure . Example of turbopump configuration

e different configurations of existing TPAs have been evaluated in a small literature review and the tool is able cover most types of modern TPA featuring inducer, up to three stage impellers for both oxidizer and fuel pumps together with a maximum of a three stage turbine on one sha . In figure 1 an example of a possible TP configuration can be seen. e convention, which is used for the acceleration and orientation, can also be found there.

. Variation of Input Values

A main feature of the tool is the generation of random numbers to take the uncertainties of performance values into account.

is Monte Carlo like approach is necessary in an early stage of design. Each input value for the calculation of the forces can be provided with a certain range of uncertainty.

e main difference to a classical Monte Carlo method is, that the set of random number is created with the following Gaussian distribution.

f (x) = 1 σ √ 2π e -1 2 ( x-µ σ ) 2 (13)
e cycle of random number generation is shown in fig- ure 2b. e uncertainties can be set globally, but also individually to accommodate the detailing and reliability of the input values. is function can also be used to assess the influence of input values on the entire design. A detail description of the tool can be found here [START_REF] Wagner | Untersuchungen zu Sekundärsystemen in Turbopumpen für Flüssigkeitsraketenantrieben[END_REF]. An enhancement of the tool is the interface to a simple bearings stiffness tool based on a general theory for elastically constrained bearings [START_REF] Jones | A general theory for elastically constrained ball and radial roller bearings under arbitrary load and speed conditions[END_REF] and the generation of input values for a dry rotordynamic analysis, which can be started directly from the tool.

e statistical procedure has been implemented instead of an analytical evaluation for the following reasons:

• Keep the tool very flexible and allow individual adjustments during the design process

• Provide probability curves of the loads for the bearing design .

RESULTS AND DISCUSSION

e method has been applied during the design of a liquid rocket expander cycle upper stage engine LOx turbopump with roughly 120 kN thrust. e design methods and configuration of this TP are described in Veggi [START_REF] Veggi | Numerical and experimantal activities on liquid oxygen turbopumps[END_REF], it is an alternative design approach to the work of Maeding and Souverein [START_REF] Maeding | A preliminary design study for an expander LOX turbopump[END_REF][START_REF] Souverein | Design and tool anchoring for a 120kN expander cycle rocket engine LOX turbopump[END_REF]. 2.5 bar Temperature at pump inlet T 1 90 K Total pressure at pump outlet p t,3 70 bar e 120 kN turbopump consists of an inducer, an impeller, a fixed bearing, the sealing unit, a loose or floating bearing and a turbine, all together on one sha . For the calculations the design values are available and a direct interface to design tools has been implemented. e dimensions for bearing, seal, and sha were estimated based on experience.

. Forces Acting on the Bearings

According to previous investigations [START_REF] Wagner | Untersuchungen zu Sekundärsystemen in Turbopumpen für Flüssigkeitsraketenantrieben[END_REF], a se ing of one hundred thousand random numbers, which leads to converged solution, and a global deviation from the input values of 20 % was used.

is strong deviation is certainly very conservative.

e two bearings are viewed separately from each other. In the present case, the fixed bearing is assigned to the pump, the floating bearing to the turbine. An inverted configuration could also be possible.

Figure 3 shows the histogram of values of the bearing forces acting on the fixed and loose bearing. All results are In the present case of 20 % variation values of axial forces with negative signs can occur. us, the direction of loads on the bearing would change. Frequent changes would increase the load on the bearing and limit the lifetime. erefore, these cases should be avoided at any time.

Radial forces are also shown in figure 3. In the fixed bearing (top figure), the static radial force are mainly caused by the spiral housing. e forces are in the same order of the axial force. is results in a balanced discharge of force into the bearing when, for example, angular contact ball bearings are used. e dynamic radials loads are plo ed in red and they are half of the radial once.

ey result from the expected unbalance of the rotating mass of the inducer, impeller, and partly the sha according to equation 12. However, the total estimated radial force, thus, the sum of the proportions of the static and dynamic forces are shown in yellow.

e loose bearing, in contrast, shows a different behav-

ior. e fraction of dynamic forces is very high, as it can be seen in figure 3(b). e loose bearing is situated close to the turbine and the assumed unbalance of the turbine disk according to equation 12. e high rotational speed, the large diameter and the necessary mass lead to these high values.

Up to now, the turbine disk design is not fixed yet. erefore, the mass is estimated within the program. A constant shear stress profile is assumed and defines the width of the disk depending on the radius. From this derived geometry, the mass can be calculated for high strength steel. In the present case of 120 kN engine, the turbine is designed as full admission turbine. us, the turbine with its axis-symmetric inflow does not create static radial forces or moments on the rotor.

e acting forces on the floating bearing results from the volute forces and are transferred by bending. Hence, their portion is very small.

. Variation of Parameters

e results in figure 3 show already the capability of the program to use random numbers for the calculation of axial and radial forces. Due to uncertainties in the prediction of the input values for the estimation of those forces, this procedure was chosen. During the design and development process, more precise values will be available and the variation of random number can be individually adapted and/or generally lowed. In the following sections the influence of both procedures will be discussed.

. . Global Variation 20 input values in the equations 3 to 12 can currently be varied in the program. is includes also the mass of the rotating parts. In this section the global se ing for the deviations in standard distribution in equation 13 is varied for the axial forces. Results for four different global se ings of uncertainties are presented in figure 4. e impact of the different values for the variation is clearly visible. A variation of 10 % already lead to a three σ deviation of more than 50 % for the axial thrust. In a comparison performed by Beck [17] between simulations and the results of equation 4 for a LOx TP with comparable dimensions, the offset ranged between 10 % to 17 %, depending on the geometrical conditions of the rotor stator cavity without flow through. erefore, 10 % variation in the values are not unlikely and the 20 % variation can be treated as an extreme case.

. . Individual Variation A er the global variation of the 20 values the influence of variation of each value is studied. Following procedure is applied for this assessment. Each values is individually varied while keeping the other values constant as the input value and then the next value is varied. For each calculation the three σ value of the forces on both fixed and loose bearing is recorded and plo ed. Now, the forces are normalised by the mean value in order to visualise clearly the deviations. In figure 5 the results of this study are shown for the fixed bearing. It should be noted that a different scale has been used for figure 5a, because the impact is higher there and greater values occur. Considering the axial thrust, the head of inducer and impeller have a strong impact. e other values for equation 4 also influence the axial force strongly. e variation of the pressure ratio, here as two pressure values (p int and p outt ) through the turbine contributes to the axial thrust. Some minor influence is given by the volume flow Q. In case of static radial forces, there is a strong dependence on the head rise H which can be seen from equation 9.

e geometrical values are kept constant for the estimation of the radial forces due to the volute around the impeller. e mass of the parts has some contribution to the radial forces, too. eir influence on the dynamic forces can been seen in the red bars, where up to now only unbalance is considered. e results for the loose bearing are shown in the two graphics of figure 6. Of course, there is no contribution to axial forces. e static radial force mainly results from bending and are transfered via the sha from the volute. In the 120 kN turbopump, the loose bearing is close to the turbine. erefore, the unbalance forces resulting from the variation of m t are very high. In figure 6, the values of exceed the scale of the y-axis.

is visualisation was chosen on purpose so that the values of the other results are still in a reasonable size.

. CONCLUSION

A program was developed to calculate the bearing forces in turbopumps of liquid rocket engines. It is applicable to most turbopump configurations. is was confirmed by intensive literature research on turbopump designs. e operation of the program and results were discussed with a LOx turbopump for a 120 kN upper stage engine. e results are influenced by two areas. On the one hand, the input values for the program depend on the quality of the design of the hydraulic and aerodynamic components, on the other hand, the validity and precision of the equation used for estimating the acting forces. Both uncertainties were taken into account by the generation of random numbers around the input values with a given variation.

e variance was adapted individually in order to investigate the influence on the resulting forces. e input values become more precise through each design loop. Numerical simulations and experiments can also be used to prove both the empirical equations used and the input values and, if necessary, further adapt them.

e presented considerations represent only one area in the investigations on secondary systems and are only one aspect in the selection of the turbopump configuration. With regard to the bearing forces, several points have to be mentioned which have a strong influence on the results and the configurations, but have not been taken into account here:

• Calculation of impeller side walls with regard to superimposed flow

• Adjustments by axial thrust compensation

• Observations outside the nominal operating point

• Stiffness of ball bearings due to the distribution of forces

• Consideration of rotordynamic effects

Work on these points is already taking place. For example, a module to deal with superimposed flow in the rotor stator cavity is available and has been validated against numerical simualtions. It has been applied on active axial thrust balancing in LRE TP [START_REF] Maier | Analytical and numerical assessment of axial thrust balancing systems in liquid rocket engine LOx turbopumps[END_REF]. An integration and loosely coupling of programs for the calculation of bearing stiffness and / or rotordynamic natural frequencies of the turbopump have also already been implemented and tested. e secondary systems remain an interface between the disciplines working on the turbopump design and contribute to the improvement of the entire system.

Figure .

 . Figure . Flow chart for random number generation

Figure .

 . Figure . Results of bearing loads normalised by axial forces (blue), radial forces statically (green), dynamic radial forces (red), combined radial forces (yellow)

Figure .

 . Figure . Influence of variation on axial forces

Figure .

 . Figure . Results of normalised bearing loads on the fixed bearing for individual variation of input values, axial forces (blue), static radial forces (green), dynamic radial forces (red)

Figure .

 . Figure . Results of normalised bearing loads on the loose bearing for individual variation of input values, radial forces statically (green), dynamic radial forces (red).

α

  AbbreviationsDLR German Aerospace Center -Deutsches Zentrum für Lu -und Raumfahrt LOx Liquid oxygen LRE Liquid Rocket Engine TP Turbopumps TPA Turbopumps Assembly TUM Technical University of Munich Symbols α Angle between impeller side walls o and the housing α w Orientation of turbopump µ Mean of distribution ρ Density σ Standard deviation b 2,ges width of impeller at outlet c 1m Meridian velocity inflow c f Skin friction coefficients a Launcher or stage acceleration c Absolute velocity d Diamter F Force g Gravitational acceleration. g = 9.81m/s 2 H Head rise k Average rotational factor m Mass n Rotational speed p Pressure Q Volume rate r Radius t Axial clearance u Circumferential speed x between the inner and outer diameters Subscripts 2 At impeller outlet

Table .

 . Applications of the KonRAT project

	Application	Launcher	Sub-orbital vehicle
	Type	Upper stage Main stage (boost)
	Cycle	Expander	Gas generator
	Fuel combination	LH2/LOx	LCH4/LOx
	rust /t	11-15	40-50
	Specific impulse /s >450	>340

Table .

 . LOx turbopump nominal operating conditions

	Property	Value Unit
	Rotational speed N	20000 rpm
	Nominal mass flow rate ṁ	25	kg/s
	Total pressure at pump inlet p t,1

ACKNOWLEDGMENTS

e KonRAT project within the Ludwig Bölkow Campus is funded by the Bavarian Ministry of Economic Affairs and Media, Energy and Technology under grant number: LABAY83E.

is financial support and the good cooperation between the project partners is gratefully acknowledged.