The application of automatic differentiation techniques for the prediction of rotor-dynamic coefficients of labyrinth seals

The article describes the development of bulk-flow code for the prediction of rotor-dynamic coefficients of labyrinth seals. The code is based on the so-called single control volume approach by Childs and Scharrer [1] and the the forces are evaluated using the automatic differentiation technique. The resulting code is very simple and provides reasonable predictions of stiffness and damping coefficients at short computational time.


Introduction
Labyrinth seals are commonly used in turbines and compressors to dissipate energy and reduce the amount of leakage flow.Unfortunately the rotor vibrations cause a non-uniformity in the circumferential pressure distribution which in turn produces forces acting on the rotor and can cause the rotor to become unstable.The prediction of rotor-dynamic coefficients and leakage flow for a straightthrough teeth on stator labyrinth seal is the main subject of this work.
There are several methods for prediction of forces generated by the leakage flows based on different models ranging from single control volume bulk-flow code (see e.g.[1][2][3]), through more complicated two or three control volumes codes [4], to full 3D simulations using finite differences or finite volumes [5][6][7].The later method based on the simulation of 3D turbulent flows gives the most complete description of the flow field and allows arbitrary geometric configurations of the seals but this comes at the price of an excessive amount of computational resources.On the other hand the single control volume bulkflow method based on the solution of cavity-averaged flow quantities predicts the leakage mass flow rate as well as the rotor-dynamic coefficients in much simpler way at much lower price.
The original bulk-flow codes [1,2] calculate the forces acting on disturbed rotor using perturbation technique, i.e. each relevant quantity is expressed as φ(t, θ) = φ 0 + φ 1 (t, θ) and a system of equation for zeroth order components φ 0 and first order components φ 1 (t, θ) is built and solved.In our code we simplify this method using automatic differentiation technique using complex extension of Clifford's dual numbers.a e-mail: Jiri.Furst@fs.cvut.cz

Single control volume model
The swirling gas passing through the labyrinth seal enters the seal at high pressure through the clearance between the first tooth and the opposite wall and enters to first cavity.Its circumferential momentum is altered due to friction with the cavity walls and then it continues through second clearance to next cavity.Once the gas passes several cavities, it leaves the seal by the last clearance at much lower pressure.The flow in each labyrinth cavity is almost uniform in radial and axial directions with exceptions of boundary layers and "jets" downstream clearances.The non-uniformity of the flow in circumferential direction due to small rotor vibrations is however very important and causes the rotor-dynamic forces.
The single control volume code assumes that the gas pressure as well as the circumferential velocity in each cavity are independent of the radial and axial coordinates, the temperature is constant across the seal and the eccentricity of the rotor is small relative to radial clearance.Using this assumptions the following system of equations can be developed: and Here t is the time, θ is the angular coordinate denoting the position in the seal, ρ i is the density in i-th cavity, A i = L i (B +CR) is the transverse surface area, V i is the circumferential velocity, ṁi+1/2 is the mass flow rate per circumferential length between i-th and i+1st cavity, p i is the pressure, τ ri and τ si are the wall shear stresses at the rotor and stator part of the cavity, ar i and as i are the dimensionless rotor and stator lengths (ar i = 1 and as i = (2B+ L i )/L i for teeth on stator configuration, and ar i = (2B + L i )/L i and as i = 1 for teeth on rotor).The geometry of the seal is given by the shaft radius R s , teeth height B, teeth pitch L i , clearance CR and by the number of teeth NT .The wall shear stresses are modeled as turbulent flow in a smooth pipe using Blasius correlation where Dh i = 2L i (B + CR)/(L i + B + CR) is the hydraulic diameter.The constants m 0 and n 0 are given for turbulent flow between smooth surfaces as m 0 = −0.25 and n 0 = 0.079.The density ρ i is related to the pressure p i using the ideal gas law p i /ρ i = rT .The boundary conditions are as follows: the pressure p 0 and the swirl velocity V 0 are prescribed at the inlet (i = 0) and the outlet pressure p NT is prescribed downstream the last tooth (i = NT ).

Axisymmetric case
In the case of central position of the rotor and stationary uniform inlet and outlet parameters the equations (1,2) have a steady axisymmetric solution (i.e.∂/∂θ = 0.).The solution is calculated with time marching using explicit Euler method: In order to analyze the stability of the numerical method the equation ( 8) is rewritten using equation of state and the model for mass flow rate as The analysis will be further simplified by taking constant value of C 0 = 0.716 as proposed by Eser an Kazakia in [2].The stability of the numerical method as well as the maximum principle for the pressure follows from the monotonicity of H c .Hence if then the method will be stable and produces monotone distribution of the pressure p i ≤ p i−1 .First two inequalities are satisfied always whereas the third one leads to stability condition for time step magnitude

Case with eccentric rotor
In order to evaluate the rotor-dynamic coefficients one assumes rotor with small eccentricity orbiting with precession speed Ω around its nominal position.The flow field is then no more axisymmetric and one has to include all terms in the eq.( 1), (2).Codes developed by Childs and Scharrer [1], Eser and Kazakia [2], or by Malvano, Vatta, and Vigliani [3] use perturbation technique splitting all relevant quantities to its mean value corresponding to axisymmetrical case and to a perturbation due to rotor eccentricity φ(t, θ) = φ 0 + φ 1 (t, θ) and developing the system of equation for first order variables φ 1 .In contrary to their approach we employ an automatic differentiation technique which greatly simplifies the code.

Automatic differentiation via dual numbers
The "classical" dual numbers extend the real numbers by adding a new element δ (similar to imaginary unit ι for complex numbers) with the property δ 2 = 0.It means that each dual number x is represented by two real numbers: the main value x v and the perturbation x p , hence x = D(x v , x p ) := x v + δx p .Assume x = x v + δx p and y = y v + δy p are two dual numbers and α a real number.
Then one can define the dual number arithmetics as follows x/y := Next, one can extend all differentiable real functions using the following rule: where f and f at the right hand side are the real function and its derivative.Let's define two operators V(x) := x v and P(x) := x p giving the main value and perturbation part of a dual number.Then one can calculate function value and first derivative of any differentiable function as: Note that the formula holds even for compound functions and thus it can be used for calculation of derivatives of any (smooth) algorithms.

Dual numbers with complex rotating perturbations
However one can use the above mentioned technique for computing first derivatives, we propose an extension of dual numbers which is more appropriate for the calculation of rotor-dynamic coefficients.We assume real main value corresponding to centered position of the rotor and complex perturbation including rotation term e ι(θ−Ωt) : Here φ v ∈ R is the main value and φ p ∈ C is the complex perturbation quotient.
Note that this extension uses the same arithmetic operations as the "classical" dual numbers.Let x = D * (x v , x p ), y = D * (y v , y p ), and α ∈ R. Then x + y = D * (x v + y v , x p + y p ), (24) αx = D * (αx v , αx p ) (25) and any smooth real function f can be extended for duals as Let's define the operator V(x) := x v in the same way as for the "classical" dual numbers and the operator P(x) := x p returns only the perturbation quotient.

Linearization using extended dual numbers
The above mentioned dual numbers with complex rotating perturbations can be easily used for automatic linearization of the code for the case of eccentric rotor making precession along a circular orbit with radius at the precession speed Ω.The (perturbed) clearance is then and it can be represented by an extended dual number hence the following identity holds Let's denote by tilde dual number representations of each relevant quantities, e.g.Ãi = L i (B + CR).The extended dual number representation of unknown flow parameters ρi , pi , Ṽi can be obtained with time marching method similar to (8), ( 9) including corrections for ∂/∂t and ∂/∂θ according to following rules: The final time marching method is then EFM 2015 The last two terms in eq. ( 34) and three terms in eq. ( 35) are the contributions from partial derivatives ∂/∂θ and ∂/∂t.The tilde denotes here extended dual numbers.Note that the mass flux rates, wall shear stresses, and isothermal gas law is also evaluated using extended dual arithmetics.
The steady solution is calculated using timeindependent boundary conditions for p0 = D * (p 0 , 0), Ṽ0 = D * (V 0 , 0), and pNT = D * (p NT , 0).The forces acting on the rotor are then evaluated in rotating coordinate frame neglecting wall shear stresses by integrating the pressure perturbation: The above mentioned method has been implemented using operator overloading and polymorphic functions in Julia language [8].

Rotor-dynamic coefficients
The forces acting on the rotor are evaluated for several precession speeds and finally the rotor-dynamic coefficients (i.e.stiffness coefficients K, k, and damping coefficients D, d) are calculated using least-squares method using linear model (see [9]) Here F r and F t are the radial and tangential components of the force, is the rotor displacement, and Ω is the precession angular speed.

Validation
In order to validate our code a simulation of flow through a straight seal with stator teeth are performed.The geometrical parameters as well as the flow parameters are given in the table 1.The forces acting on the rotor are calculated for 10 precession speeds Ω in the range 0 rpm to 8000 rpm and the rotor-dynamic coefficients are calculated using least-square method.The figures 2 and 3 show the cross-coupled stiffness and direct damping coefficients (hence the components of circumferential force F θ ) for several inlet swirl velocity ratios V 0 /(Ω r R s ).The results obtained with current code labeled by "Fuerst, one volume" are compared to experimental results of Childs and Scharrer [10], "classical" one   control volume code [1], and with the two control volumes code [11].One can see that the results for cross coupled stiffness and direct damping corresponds very well to experimental data as well as to results obtained with other codes.
On the other hand there is quite large discrepancy in direct stiffness K obtained with current code and experimental data (see figure 4).This discrepancy is caused probably by the fact that the radial force component shows non-linear dependency on the precession speed (see fig-    coupled damping c depends strongly on the precession assumed during computation.Note that this is not the case for the circumferential component of the force which is almost linear, see figure 6.
Unfortunately there is no or very little informations about precession speeds for the experimental data as well as for the results of Forte and Latini [11].Nevertheless the description of the experimental apparatus states that  the rotor makes linear side-to-side motion [10].Similarly the "classical" perturbation based codes including those of [11] assume elliptical motion.In both cases the elliptical as well as the linear motion can be decomposed onto two circular orbits with opposite precession speeds.
Therefore the calculation is done here also for two opposite precession speeds Ω =−8000 rpm and 8000 rpm.The rotor-dynamic coefficients are then evaluated directly from two forces without using least-square method.
One can see that there is almost no difference in circumferential components k and C but the agreement with experimental data as well as with the computations of other authors is much better in this case.

Conclusion
The original automatic differentiation technique based on dual numbers with complex rotating perturbations allows rapid development of a code for dynamic analysis of flows in labyrinth seals without the need of tedious and errorprone linearization.The code is further simplified using modern programming techniques (operator overloading, polymorphic functions, ...).
The code is validated against experimental data [10] and numerical results of other authors using two variants EFM 2015 02025-p.5 of "classical" bulk-flow method for the case of a lookthrough seal with 16 stator teeth.The results show that the cross-coupled stiffness and direct damping coefficients are predicted in accordance to available data.Due to nonlinear dependence of radial force on the precession speed the agreement of direct stiffness coefficient with available data depends highly on the precession speeds used during the computation.

DOI: 10 Figure 1 .
Figure 1.Geometry of straight through seal with stator teeth.

ure 5 )
and therefore both the direct stiffness K and cross-EPJ Web of Conferences 02025-p.4

Table 1 .
Parameters of test seal