Lei Huang

Hua Chen

A NOVEL DESIGN METHOD OF VARIABLE GEOMETRY TURBINE NOZZLES FOR HIGH EXPANSION RATIOS

Keywords: turbine nozzle, aerodynamic design, conformal mapping 1 Inlet 2 Outlet Ref Referential

In variable nozzle geometry turbines (VNT), opening of the nozzles is used to control turbine mass flow and expansion ratio, allowing more turbine power to be generated over wider operating conditions. In turbocharged vehicles, the nozzles are 'closed' to provide high boosts for engine and vehicle acceleration and for engine braking assistance. At the both conditions, high nozzle expansion ratios are created, and shockwaves may generate from the nozzles. These shocks reduce turbine efficiency and they can cause high cycle fatigue (HCF) damage to the downstream rotor blades.

Design of high expansion ratio radial nozzles is difficult for

VNT because transonic flows are very sensitive to small geometry changes, and the large semi-vaneless space created by the nozzles makes the design a tricky business. Shock minimised nozzle designs are therefore often achieved by auto-optimisation technique. While design targets may be achieved, this technique does not offer sufficient insights into how the optimal flow field has been derived, so the same optimisation procedure has to be applied to every new design. In this paper, a new design method that overcomes this problem is proposed. The method first uses a conformal mapping to transfer a radial nozzle from the r- plane into the x-y plane. Mapped nozzle displays amplification of supersonic acceleration and diffusion. This is explained by the curvature changes brought about by the mapping, and a link between the shock strength and the flatness of the suction surface of the mapped nozzle is found. The amplification and the link can be utilised to design nozzles with reduced shock loss in the x-y plane first and then mapped back to the r- plane. Two nozzles for 6:1 expansion ratio were designed in this way and CFD results show a significant reduction of nozzle loss. The nozzles were also checked for fully open condition and no performance penalty was found.

Nomenclature

These shocks generate losses and reduce turbine efficiency.

They will also interact with the downstream rotor and can cause HCF of the rotor blades. The HCF is one of the major concerns in VNT rotor design [START_REF] Chen | Turbine wheel design for Garrett advanced variable geometry turbines for commercial vehicle applications[END_REF]. Efforts were made to design VNT nozzles so that they generate no or weaker shocks and have less losses. Yang et al. [START_REF] Yang | Investigations on the generation and weakening of shock wave in a radial turbine with variable guide vanes[END_REF] proposed an increase of nozzle vane number to provide better guidance to the flow to reduce shock strength. This method however increases nozzle surface friction loss, which is often unacceptable. Zhao et al. [START_REF] Zhan | Numerical Investigation of a novel approach for mitigation of forced response of a variable geometry turbine during engine braking mode[END_REF] suggested adding grooves to the suction surface of the vanes, and their simulation showed that the shocks were weakened and rotor excitation was reduced by 30%, but the nozzle loss was increased. These references dealt with thick vanes used in Honeywell's AVNT, a special type of VNTs that requires thick nozzle vanes. For ordinary VNTs that accept thinner vanes, the optimum vane geometry can and was obtained through auto-optimisation technique [START_REF] Chen | Turbocharger turbine design by autooptimisation[END_REF]. A geometry generating tool coupled with CFD code and driven by an optimisation software was able to generate a geometry that minimises the shock and significantly reduces flow losses. One of such an example is given in Figure 2. It can be seen that the suction side diffusion is minimised. We demonstrate these through an example of improving an existing nozzle designed by the auto-optimisation. The baseline nozzle to be improved is a VNT nozzle for automotive application which was designed using the autooptimisation technique mentioned earlier. Although there is a flow acceleration after geometric throat of the nozzle in the semi-vaneless space along the suction side, the entire flow is subsonic therefore no shocks are produced. There is a small flow diffusion after the acceleration toward vane trailing edge, in order that both the pressure side and the suction side have the same pressure when they meet at the edge. This small diffusion can also be seen in Figure 2. Our design target is a new nozzle that maintains this level of performance at this expansion ratio, while improves upon it at higher expansion ratios when the turbine is in engine braking mode.

The methodology and its application

Conformal mapping

The nozzle has a long semi-vaneless space, and the radial inflow nature of the nozzle makes it difficult to gauge the flow area variation after the nozzle throat and other key geometric features that may affect the flow. So a conformal mapping was first carried out to map the nozzle from the r- plane into the x-y plane, where rref is a reference radius, and rref was taken here as the mean of the maximum and minimum radii of the nozzle. This type of conformal mapping has been used in compressor diffuser vane design [START_REF] Japikse | Centrifugal compressor design and performance[END_REF], and is also used by other industries in radial turbine nozzle design [START_REF][END_REF]. Mapped nozzle is shown in Figure 4a, and CFD predicted Mach no. distribution at the same expansion ratio of 2:1 for this mapped nozzle is given in Figure 4b. While no shock exists before the mapping, a shock is now generated from the suction side of the nozzle vane after the mapping. At the first glance, this result seems to suggest that the conformal mapping may not be useful because it does not reproduce the original flow field.

Since eq. ( 1) is a conformal mapping, the vane angle of the nozzle in Figure 3a is kept in Figure 4a, and this is a useful feature of the conformal mapping. Surface length of the nozzle vanes on the other hand is not maintained after the mapping, and this mainly produces the discrepancy of the flow field. Using eq. ( 1), the relationship between the total differentials of any surface length in the two mapping planes can be obtained, along the suction first accelerates and then diffuses. The only difference seems to be that the one after the mapping displays much stronger acceleration and deceleration due to the reason explained above. This suggests that if a nozzle is optimised in the mapped x-y plan by controlling these accel-eration and deceleration, it may work equally well when it is mapped back to the r- plane using the same conformal mapping. This interpolation of the CFD results offers a possible way to design radial nozzles, because working in the xy plane is always easier than in the r- plane. One may also reason that when nozzle expansion ratio increases from current 2:1 to higher values, flow in the baseline nozzle (in the r- plane) may become similar to one shown in Figure 4b, that is, choking may happen near the throat followed by a supersonic acceleration along the suction surface which then ends with a shock. Figure 6 shows the blade loading of the baseline nozzle when the pressure expansion ratio is 6:1.

           
The similarity of the suction side loading is clear to that of the mapped nozzle in the x-y plane subject to 2:1 expansion ratio (yellow curve in Figure 5). This implies that to design a radial nozzle for high expansion rations, one may use the conformal mapping of eq. ( 1) and design it at a lower expansion ratio in the x-y plane.

To reduce the shock loss and the shock-related excitation, the supersonic acceleration in the semi-vanelss space needs to be reduced. Flow acceleration is controlled by two factors, the area schedule that affects mean flow velocity and surface curvature that influences local acceleration. The flow picture in Figure 4b provides some hints on how the passage area is seen by the flow in this space, and these will be looked into in future. Local flow acceleration is studied first because it provides a direct link between the geometry and the flow locally thus is useful in design optimisation. As will be seen later, the results justify this choice. In the x-y plane, a convex surface will produce local flow acceleration, and in this regard, the surface may be measured by the change of its first derivative. To reduce supersonic acceleration, a surface should be flat so the derivative of the surface should remain constant.

Results and discussion

The two new nozzles were mapped back to the r- plane by conformal mapping of eq. ( 1).

  The turbine of a vehicle turbocharger is subject to a wide flow range and demanding power requirements. At low mass flows, high efficiency and power output are required to improve engine torque and transient response, while a high flow capacity is needed for engine rated power and to reduce engine pumping loss at high speeds. To meet these challenging needs, variable Nozzle Turbines (VNT) employ a nozzle ring upstream of the rotor, by changing the setting angle of the nozzle vanes, different values of nozzle throat area and vane exit angle can be achieved. When an acceleration of the engine or turbine is required, the nozzles are closed to reduce the throat area and make the nozzle exit flow more tangential. So a higher nozzle exit velocity is achieved in a more tangential direction, which enables the rotor to produce more Euler's work. Closing the nozzles can also be used for engine deceleration or braking purpose. This increases the pressure expansion ratio across the nozzles and establishes a high back pressure at engine exhaust manifold, adding to the pumping loss of the engine.When the flow inside the nozzles of a VNT turbine is subsonic or the expansion ratio of the nozzles is well below 2, geometry of the nozzles matters little to the flow loss within the nozzles. When the expansion ≥ 2, the flow in parts of the nozzles may becomes supersonic, and shockwaves could generate from nozzles. Such expansion ratios exist at engine braking or during unconstrained turbine acceleration. Figure1shows calculated Mach no. of a VNT nozzle ring under such a condition, shockwaves generated from the underside or suction side of the nozzles are visible.

Fig. 1 Calculated

 1 Fig. 1 Calculated Mach no. of a VNT nozzle ring under a high expansion ratio [1]

Fig. 2

 2 Fig. 2 Calculated total pressure distribution (left) and static pressure loading (right) of an optimisedVNT nozzleWhile the auto-optimisation method is capable of producing optimal designs, they offer little insights as why a particular geometry is better than others, or how the optimal flow such as the vane loading shown in Figure2is linked to the geometry in the same figure. One of the difficulty in design of radial nozzles for turbocharger VNTs is that the adjacent pair of nozzle vanes usually do not form a proper 'nozzle' in the sense that a large part of the flow region is so called semi-vaneless space (See Figure1). The auto-optimisation technique does not provide a clear design guideline for this region, thereby every new nozzle must be designed using the same black box procedure.

2. 1 Fig. 3

 13 Fig. 3 Baseline VNT nozzle and Mach number distribution at expansion ratio 2:1

  Figure3ashows the nozzle at a closed position and Figure3bgives predict-ed Mach no. distribution at this position and an expansion ratio of 2:1.The computation that produced this result and others followed was carried out with commercial software Fineturbo in which 2D Navier-Stokes equations were solved with the Spalart-Allmaras turbulence closure. Turbine housing and turbine wheel were not included in the calculation, they merely proved the inlet and outlet boundaries. The term of expansion ratio here and afterward refers to the total-tostatic pressure ratio across the computational domain, and not the turbine stage. When the nozzle is closed at low rotor speeds, the mass flow of the turbine stage is small. In this condition the pressure drop across the rotor is small with the most of stage pressure drop happens inside the nozzle. On the other hand, when the nozzle is open at high rotor speeds, the mass flow of the stage is larger and so is the pressure drop of the rotor. In this condition, nozzle expansion ratio is a smaller portion of the stage expansion ratio.

Fig. 4

 4 Fig. 4 Geometry and Mach number distribution of the baseline nozzle at expansion ratio 2:1 after conformal mapping of eq. (1)

Fig. 5 1 2. 3

 513 Fig. 5 Comparison of loading for vanes before and after the conformal mapping, nozzle expansion ratio = 2:1

Fig. 6 1 Figure 7

 617 Fig. 6 Loading of baseline nozzle at 6:1 expansion ratio

Fig. 8

 8 Fig. 8 Geometry of nozzles A, B and C in the x-y plane

Figure 9

 9 compares their geometry with the baseline's. While the suction surface of the new nozzles in the x-y plane is flatter than the baseline, Figure 8, they become curvier in the r- plane. Such a change is less intuitive if the geometry modification was carried out in the r- plane.

Fig. 9

 9 Fig. 9 Geometry of three nozzles in the r- plane

Fig. 10

 10 Fig. 10 Performance of new nozzles at high expansion ratios, nozzles are closed.

Fig. 12 variation at nozzle exit at expansion ratio 6 : 1 ,

 61 Fig. 12 variation at nozzle exit at expansion ratio 6:1, nozzles are closed.

Fig. 13

 13 Fig. 13 Nozzles A, B and C in fully opened position

Table 1 Performance of three nozzles in Figures 7 & 8

 1

Conclusions

A new design method of radial nozzles for high expansion ratios has been developed. It first uses conformal mapping of eq. ( 1) to map the nozzle geometry from the r- plane into the x-y plane then works with the nozzle in the x-y plane to optimise its geometry before mapping the nozzle back to the r- plane.

The total differentials of the correspondent nozzle surface arc length in the two planes of eq. ( 1) are linked through eq.

(2), which can be used to explain the flow field variation before and after the mapping. When the mean radius of vanes is used as the reference radius in the mapping, the curvature of the front part of the suction surface is amplified and the rear part reduced. This creases a strong supersonic acceleration along the surface leading to a shock termination which would not happen in the r- plane under moderate expansion ratios. This suggests that should a nozzle in the xy plane work well under these expansion ratios, it could work equally fine at higher expansion ratios after mapping to the r- plane.

The new method was applied to a VNT nozzle for automotive turbochargers, and two new nozzles were designed. By making the suction surface in the x-y plane flatter than the baseline nozzle after the throat, the new nozzles display less supersonic acceleration and weak shocks than the baseline.

When mapped back to the r- plane, they both show better vane loadings, lower losses than and reduced exit pressure variations to the original nozzle which was previously designed by an auto-optimisation method for a slightly different operating condition.

Current method does not consider area effect after geometric throat. As flow is not fully choked at the throat under moderate expansion ratios, modification of suction surface geometry after the throat will affect aerodynamic throat area, and can change nozzle mass flow under such expansion ratios. This is an area of the new method to be improved. Second thing to be understood is why in the x-y plane, Nozzle B appears more efficient than Nozzle C (0.882 vs. 0.871), but when mapped to the r- plane, Nozzle C is slightly better than Nozzle B (0.9032 vs 0.9027).


Dennis Toebben
email: toebben@ikdg.rwth-aachen.de

Piotr Luczynski

Manfred Wirsum

Wolfgang Mohr

Klaus Helbig

S Y Mpo

O Tatin

Optimized Approach for Determination of the Solid Temperature in a Steam Turbine in Warm-Keeping-Operation

Keywords:

e determination of the temperature distribution in the thick-walled components in steam turbines is increasing in relevance. Due to the growing share of volatile renewable power generation, power plants with a high flexibility and a high integral efficiency are required. e current operational conditions lead to high thermal stresses inside the heavy components and thus to a reduced lifetime. To improve the ability for a fast start-up, the steam turbine can be kept warm during a longer period of stand still. erefore, information about the metal temperature inside the rotor and the casing are crucial. However, the temperature distribution of the inner casing and especially the rotor cannot be measured without high additional effort. us, a calculation model with sufficient accuracy and also manageable calculation effort is required. In the present work, a hybrid -numerical FEM and analytical -approach has been developed to calculate the solid body temperatures of a steam turbine in warm-keeping operation in a most efficient way. e presented model provides accuracy of nearly % compared with Conjugate-Heat-Transfer simulations with a simultaneously reduced calculation effort by a factor of more than .

INTRODUCTION

Today, power plants using steam turbines are responsible for more than 50% of the worldwide power generation [ ]. Most of these power plants have an age of thirty years and more [ ]. Within the last decade, renewable energies are continuously gaining influence in the energy supply sector.

ese renewables lead to a fluctuating power feed-in into the grid which needs to be regulated by the conventional generation system. Most of the installed capacities of steam turbine power plants are designed for basic load with a limited ability for high load ramps and fast start-ups. Consequently, the existing power plants need to be revised and improved to increase the operational flexibility.

Further decreasing operating hours of the conventional power plants and thus an increasing cost pressure in the energy sector are another result of the growing share of renewable energies. For the decision whether a power plant goes on stream, the importance of the start-up costs is still increasing. e base load efficiency is no longer critical, but rather the integral efficiency including the start-up and shutdown procedure as well as base, part and minimum load operation.

One option to improve the flexibility and the integral efficiency is to keep the steam turbine in a warm condition between a shut-down and the next start-up.

e required compensation of the losses can be achieved by the use of heat blankets [ ], gland steam [ ] or even by hot air, which is, for example, passed through the steam turbine [ ]. e last approach is the basis for the present investigations.

For an efficient, fast and lifetime conserving start-up the knowledge of the temperature distribution within the thick walled components of the steam turbine is crucial. Especially the temperature of the rotor is a determining factor which is not measured in commercially operated steam turbines. e determination of the rotor temperature distribution can only achieved by large technical and financial efforts [ ], [ ], [ ].

us, a retrofit of such a measurement system is currently not appropriate for commercial use. Hence, a different solution for the rotor temperature prediction needs to be found.

Considerable improvements of the computing power and numerical simulation tools enable a relatively detailed investigation of the heat transfer and the flow phenomena in steam turbines [ ]. However, detailed fluid-dynamic calculations can only be conducted for a limited number of numerical nodes so far. is leads to either smaller geometries or to coarser numerical meshes and consequently to a reduced accuracy. A highly accurate coupled fluid and solid D-simulation of a whole multistage steam turbine can only be conducted with prohibitive expense. erefore, a simplified calculation model is required.

Several analytical approaches for calculating the heat transfer within a steam turbine are known from literature. Brilliant et al. [ ] developed a separate calculation model for both an impulse and a reaction steam turbine, using differ-ent Nusselt-correlations from literature for various surfaces.

ese Nusselt-correlations are used to calculate the convective heat transfer coefficient (HTC) between steam and the solid structure to serve as an input for further numerical FEM simulations. ey calculate different discrete operating points from start-up till base load. Particularly in the start-up phase the models are not able to reproduce the test data in an appropriate way. Rzadkowski et al. [ ] analyse the transient thermal behavior of a steam turbine during its start-up in discrete time steps. ey use Reynolds-averaged Navier-Stokes equations (RANS) in a CFD simulation for calculating the heat transfer coefficients. Furthermore, Nusselt-correlation taken from literature were used to calculate the heat transfer within the labyrinth seals and due to condensation within the flow channel. ese HTCs were used to couple the flow and the solid body in a FEM model.

A further approach was developed by Pusch et al. [ ]. e aim of the investigation was to simulate the start-up procedure of an intermediate pressure steam turbine to analyze low-cycle-fatique. erefore, a numerical FEM model was developed including the rotor and the inner casing. Various analytical heat transfer correlations were used to calculate the heat transfer. e absorption of heat by the blades and vanes is calculated by a modified correlation which converts the HTC at the blade surfaces (Nusselt-approach developed by Traupel et al. [ ]) to an equivalent HTC at the endwall. us, only the blade and vane root needs to be integrated into the model.

Moroz et al.

[ ] conducted a transient thermo-structural analysis of a combined high pressure / intermediate pressure steam turbine during a start-up. For these investigations a FEM model was developed, which includes the rotor and the inner casing of the steam turbine. e FEM model was separated into several heat transfer zones in which analytical heat transfer correlations from literature were used. In addition to the convective heat transfer, the model is able to calculate condensation within the individual zones. Based on these analytical equations, the temperature distribution of the rotor was analyzed. ese calculations were used to determine the rotor displacement in several discrete time steps during the start-up phase.

Wang et al. [ ] investigate the impact of the steam temperature fluctuation during steady-state operation on the temperature field within the rotor.

e heat transfer between the steam and the rotor is calculated by analytical equations concerning the rotor surface, the seal structure and the blade groove. e heat transfer coefficients serve as an input parameter for a FEM rotor model. It is assumed that all heat fluxes received by the blades are conducted into the rotor through the contact between blade root and blade groove. Similar to the approach of Pusch et al. [ ], the blades are simplified as the blade roots and integrated into the model. In addition, an analytical approach for determining the contact resistance at the blade-rotor-connection depending on the rotational speed is included.

Born et al. [

] developed a numerical model of an intermediate pressure steam turbine for the investigation of steady state operation and natural cooling effects. e model includes the inner and outer casing and the labyrinth seals. e rotor is not modeled. e authors used several analytical heat transfer correlations based on approaches from literature for determining the heat transfer by forced and free convection, radiation and conduction.

All these investigations were made with steam as operational fluid. e present paper deals with an optimized approach for the determination of the solid temperature in a steam turbine during warm-keeping-operation with hot air. Hence, the previous calculation approaches are inappropriate. A simplified Hybrid-FEM model (HFEM) has been developed, which uses analytical heat transfer correlations coupled with a FEM model for time-optimized calculation of the solid body temperature distribution. e model consists out of a repetitive stage including inner casing, rotor, vane, blade and shrouds. In comparison to the approaches known from literature, the present model uses detailed numerical Conjugate-Heat-Transfer (CHT) investigations [ ] to develop suitable heat transfer correlations for the most relevant surfaces concerning the heat transfer. e FEM model was developed on basis of the CHT model grid. e fluid grid was removed and the solid grid was coarsened to optimize the calculation time. us, on the one hand heat transfer and temperature analysis can be made with minimized calculation effort and on the other hand high resolution.

In the first part of the present paper, the model setup and the boundary conditions of the warm-keeping operation with air are introduced. Subsequently, the analytical calculation models for convection and radiation are presented. Finally, a comparison between the results of the numerical CHT model and the Hybrid-FEM model is conducted, which shows the advantages of the simplified model. .

NOMENCLATURE

HYBRID-FEM MODEL

e Hybrid-FEM model uses hybrid calculation approaches: e heat transfer within the solid body is calculated by the use of a Finite Element approach while the heat exchange by forced convection, radiation or contact heat transfer is calculated by a semi empirical approach based on detailed Conjugate-Heat-Transfer (CHT) simulations.

Figure provides detailed information about the model setup of the Hybrid-FEM model. In the case of fast and frequent start-ups and high load ramps, the impact on the lifetime is crucial. erefore, information about the temperature and stress distribution of the turbine in such off-design operations are important. e Hybrid-FEM model can be applied to provide these information. A numerical fluid and solid mesh for one repetitive stage is generated from the basis turbine design. With this mesh, a CHT model is developed which is used for calculating a characteristic map regarding the operating range of interest (e.g. start-up phase, warm-keeping or pre-warming). To secure a appropriate resolution of the thermal boundary layer, the near wall heat transfer is modeled based on the Low-Reynolds-Number method (kω-SST turbulence model), which resolves the details of the boundary layer profile by using very small mesh length scales in the direction normal to the wall. e fluid boundary layer y+ is smaller then . to secure a fine mesh resolution. For transferring these calculation results into the Hybrid-FEM model, modified Nusselt-Correlations are used to fit the characteristic HTC map. In contrast to the requirements of the CHT model, a coarse solid body mesh is sufficient for the Hybrid-FEM model's purposes. Hence, the fluid mesh is removed from the initial mesh and the solid mesh is coarsened.

is repetitive stage mesh can be multiplied to any needed number of turbine stages (neglecting the flow channel expansion). In a final step the adapted analytical heat transfer correlations are integrated into the model. e result is the Hybrid-FEM model which is able to calculate the temperature and stress distribution within large scale turbine geometries with improved accuracy.

is is reasoned by individually modified heat transfer equations, which provide a reduced as well as at least manageable calculation effort for simulating a whole multistage steam turbine.

In comparison to other approaches, the Hybrid-FEM model uses adapted instead of general heat transfer correlations for a specific turbine geometry and specific operation conditions. In the next step of the ongoing investigation the Hybrid-FEM model will be validated against experimental data. Due to the strong basis of the CHT simulations just a few regular operation measurements are sufficient to validate the model. e presented Hybrid-FEM model is used to calculate a warm-keeping operation of a steam turbine. In general, the warm-keeping process is designed to increase the flexibility and reduce the lifetime consumption.

e warm-keeping procedure is characterized by a hot gas (e.g. air) which is passed through the turbine to compensate the heat losses primarily to the ambient and the lubricant oil system. e steam turbine rotates during the warm-keeping process. Besides the conventional rotation direction, it is possible to rotate the steam turbine rotor in reverse direction to increase the warm-keeping efficiency [ ].

e operating case with conventional rotation direction is called quadrant (Q ) and the other one quadrant (Q ).

e heat transfer and flow phenomena in warm-keeping operation have already been investigated [ ]. For this analysis a numerical D single stage CHT-model was used characterized by the following assumptions:

• Steady-state CHT-simulation

• Single blade passage (approx. 2.5 • ) with periodic boundaries in circumferential direction

• Frozen rotor-stator interface

• Fully turbulent ideal gas with constant fluid properties

• No heat radiation

• Low Reynolds kω-SST turbulence model

• Constant material properties for the solid state

• Fluid boundary layer: y+ ≤ .

Toebben et al.

[ ] figured out that in warm-keeping operation most of the heat is transferred by forced convection from the fluid to blades and vanes. ese kind of fins conduct the heat through the blade-rotor respectively vane-casing connection to the main turbine components. erefore, the contact heat resistance between these surfaces is important for calculating the overall heat transfer and the temperature distribution within rotor and casing. Besides convection, heat radiation within the flow channel and between the inner and outer casing occurs. e influence of the heat radiation depends on the surface temperatures. In contrast to the investigations of the flow phenomena and the convective heat transfer in which the radiative heat transfer could be neglected due to similar suface temperatures (steady-state operation), the relevance of radiation increases for transient operations. us, the three heat transfer phenomena convection, radiation and contact heat transfer need to be integrated into the Hybrid-FEM model by using analytical calculation approaches (Fig. ).

e coarsened solid mesh of the single stage Hybrid-FEM model has only 5.1 • 10 4 nodes. In comparison to the solid mesh of the CHT model the reduction factor is . e quality and accuracy of the coarsened mesh was proven in a mesh study.

e reduction of the number of nodes is permi ed because there is no fluid mesh and thus, the interpolation between fluid and solid mesh is no longer crucial. e fluid mesh usually has a much higher number of nodes than the solid mesh to ensure a high resolution of the flow phenomena and the thermal boundary layer. e numerical solver needs to interpolate between the fluid and the solid mesh. Hence, the resolution of the solid mesh is also crucial for the accuracy of the heat transfer calculation. Once the fluid mesh is no longer needed, the solid mesh can be coarsened significantly.

e total mesh reduction factor of the final Hybrid-FEM mesh compared to the fluid (3.6 • 10 6 nodes) and solid (3.2 • 10 6 nodes) CHT mesh is . With this mesh resolution it is possible to simulate a whole steam turbine with acceptable computational effort.

Before going into detail of the analytical calculation approaches, the relevant surfaces of the model should be introduced. e single stage model has eight main surfaces, which are responsible for the main heat exchange between fluid and solid as well as for heat radiation. a simulation with equal temperature on each surface was conducted. To consider the impact of different surface areas, both the absolute and the area related value of the normalized heat flux is given. e blade and the vane (including the tip) receive each about 48% of the heat in warm-keeping operation. is effect is enhanced by the large surfaces of blade and vane compared to the other surfaces, shown by the area related value.

e remaining heat is absorbed by the shrouds of both, blade and vane. Whereas, the surfaces S Housing and R Rotor can be neglected. e separation of the blade surface into blade and blade tip serves to set separate fluid reference temperatures (due to the cooling effect of the shroud) and to include heat radiation between blade tip and blade shroud. For further evaluation only the four main surfaces Vane, Vane Shroud, Blade, Blade Shroud are considered.

For these four main surfaces an analytical approach for the convective heat transfer is developed. e numerically calculated heat transfer coefficients (HTC) hi,n and Nusseltnumber Nu C HT, i, n , which definitions are given in Eq. and Eq. serve as basis. e index i refers to the individual surface and the index n to the inlet plane of the flow channel section in which the regarded blade/vane row (S /R ) is located. Qi is the area averaged heat flux of the surface i with the surface area A i . e reference temperature Tn is the mass flow averaged temperature of the fluid passing the plane n. Ti is defined as the area averaged temperature of the surface i.

hi,n

= Qi A i ( Tn -Ti ) ( ) Nu C HT, i, n = hi,n • l i λ f l uid ( )
e heat transfer within the flow channel of the steam turbine in warm-keeping operation is mainly influenced by the inlet air mass flow ṁ and the rotational speed n. Due to the assumption of atmospheric back pressure and very low mass flow rates, the inter stage pressure drop is very small and thus, the pressure influence is negligible. In order to describe the input parameters, the Reynolds-number in circumferential flow direction (Eq. ) and in axial direction (Eq. ) is defined. e radii r 1 and r 2 are the limits of the inner and outer wall of the flow channel. e characteristic length of the regarded surface is defined as l i .

Re u,i = πn(r 1 + r 2 )l i 60ν ( ) Re a x,i = ṁl i ρπ(r 2 2 -r 2 1 )ν ( )
For the development of the convective heat transfer correlation for the warm-keeping process in two different operations (Q & Q ) a physical correlation approach is used. e investigations have shown that Re a x & Re u are the dominating parameters for the description of the Nusselt-number. erefore, a Nusselt-approach Nu i (Eq. ) was chosen, including empirical factors a, b 1 and b 2 as well as the sum of both Reynolds-numbers and the Prandtl-number, to consider the fluid properties.

Nu i = a i (Re b1 i u, i + Re b2 i a x,i )Pr 1/3 ( )
For each surface i and each operation condition (Q or Q ) one Nusselt-correlation is defined by individual empirical factors a i , b1 i and b2 i . ese coefficients are determined by a minimization of the residual sum of squares based on the CHT results. e Nusselt-correlation of Eq. is plo ed in Fig. for the conventional rotation direction (Q ) and in Fig. for the reversed rotation direction (Q ). ese figures also show the numerically calculated Nusselt-numbers for the convective heat transfer over the axial Reynolds-number with logarithmic scaled axes. In accordance to the findings of the flow analysis [ ], the numerical results are influenced by secondary flow structures. In general, the Nusselt-number rises with increasing axial Reynolds-number and also with increasing circumferential Reynolds-number. Because of the single stage model structure, the vane is only influenced by the rotational speed when the axial Reynolds-number is extremely low. e blade shows a distinct dependency on the rotational speed. However, kinks in the Nusselt-number trend can be observed due to the secondary flow phenomena. In the vane and blade shroud, no distinct distribution of the calculated Nusselt-numbers can be found, especially for the reversed rotational direction (Q ).

e aim of the Nusselt-correlation is to reproduce the trend of the heat transfer based on physical relations. us, the correlation is not able to accurately describe the irregularities regarding changing vortex structures. To evaluate the accuracy of the correlation, the relative root mean square error r RM SE is introduced in Eq. , which represents the relative average deviation of the fi ed response value ŷ(x i ) from the original data value y i . . .

r RM SE = 1 n n i=1 y i -ŷ(x i ) y i 2 ( )
e r RM SE is given in Tab. for each surface and the different operation conditions Q and Q .

e evaluation of the predicted values is made by the same CHT results, which are used for deriving the correlation.

e r RM SE of the S Vane surface is less than ten percent. Due to the nonrotational inlet flow, the Nusselt-numbers of the vane show a distinct trend. e average r RM SE is in a range of to percent. To consider the influence of the different surface areas, the area related average values of the r RM SE are mentioned.

. ANALYTICAL APPROACH: RADIATION

Besides the convective heat transfer, the heat radiation can also significantly influence the general heat exchange. Within the flow channel of a steam turbine, there are several complex geometries.

e radiation heat exchange takes place between all infinitesimal surfaces. Hence, the calculation of these heat exchanges is very complex.

A few commercial numerical FEM so ware tools are able to calculate the radiative heat exchange without a fluid mesh. In the present work, the so ware ANSYS CFX . has been used, which does not provide this opportunity. e numerical CHT simulation (fluid & solid) of simultaneously convection and heat radiation reaches high radiation imbalances. erefore, the heat radiation was not considered in the simulation of the convective heat transfer. However, to use the synergy of the CHT model and the Hybrid-FEM model, ANSYS CFX was chosen as so ware basis.

For the integration of the heat exchange by radiation into the Hybrid-FEM model, an analytical approach has to be found. In accordance with the literature [ ] an approach for the radiative exchange between grey Lambert radiators is chosen. Each of the surfaces, which were defined before, exchanges radiative energy depending on the geometry of the surface and the position within the flow channel. To consider this effect, a view factor F12 has to be found, which represents the share of radiative energy from a surface to the surface . Equation describes the commonly known analytical approach for heat radiation Qrad regarding the view factor and the emissivity ǫ i of the related surface. In conservative estimation the emissivity is assumed as ǫ i = 1 for all surfaces. Due to the oxidation layer, the surfaces within the flow channel have a dark grey color.

Qrad,1→2 = σ(T 4 1 -T 4 2 ) 1-ǫ 1 ǫ 1 A 1 + 1 A 1 F12 + 1-ǫ 2 ǫ 2 A 2 ( )
For the determination of the view factors, the so ware tool ANSYS Fluent was used. erefore, the numerical mesh was imported into Fluent. As calculation method Monte Carlo with the option surface to surface was used. Gas absorption respectively gas radiation is neglected. To consider also the share of radiation, which is exchanged between several stages, the single stage mesh was duplicated three times and connected to a simplified three stage model (without flow channel expansion). e surfaces of each stage were combined to surface groups as shown in Fig.

. e share of radiation energy, which is absorbed by other surfaces such as blades in the second row or parts of the rotor and housing, is less than ten percent. Besides the blade surface, the shroud surfaces (green) also exchanges radiation energy. e blade shroud collects a significant share of the radiation energy from the blade tip. Whereas, the shroud with its several fins absorbs most of the emi ed energy by itself. For the vane (blue) it is similar. us, with two respectively four view factors (incoming & outgoing) for each of these main surfaces approximate % of the irradiation energy can be covered. In Fig. , the continuous arrows demonstrate these main heat exchanges.

e arrows with do ed lines show the secondary exchanges.

For the verification of theses calculated view factors with the solver ANSYS Fluent, the solver ANSYS CFX was used.

e view factors were calculated separately. For this purpose, three simulations for the determination of each view factor were conducted (Monte Carlo, surface to surface). In the first simulation all other surfaces apart from that two for which the view factor should be determined were defined as "openings" with a surface temperature of zero Kelvin. As a result, the emi ed and absorbed heat flux is calculated, including the share of irradiation which is absorbed by the emi ing surface itself. erefore, two further simulations were conducted to each calculate only the irradiation heat flux which is absorbed by the regarding surface itself. e average relative error of these two calculation methods is five percent. is value meets the expected accuracy.

ϕ ϕ ϕ ϕ F 1 F 3 F 2 F 4 F 5 F 6 F 7 F 8
.

MODEL COMPARISON

For the comparison of the HFEM and the CHT model, all operating points were recalculated with the HFEM model. Due to the high numerical imbalances in the energy equations while solving heat radiation, the radiation influence had to be deactivated for no distorting the results. In It can be seen that the results of the area related heat fluxes match in good agreement. For a quantitative analysis, the rRMSE (similar to Eq. ) of each surface is calculated and shown in Tab. . is error estimation takes the general error of the model into account, which is also influenced by the calculation of the fluid reference temperature and the area averaging, in contrast to the rRMSE of the Nusselt-number (Tab. ). It can be noticed that the rRMSE of the heat fluxes is lower than that only of the Nusselt-number. e average area related rRMSE which considers the different surfaces areas and thereby, the impact on the general heat transfer, is in between . and . . . .

_ q i = _ q 0;i (CHT) 0 0.2 0.4 0.6 0.8 1 _ q i = _ q 0;i (HFEM)
_ q i = _ q 0;i (CHT) 0 0.2 0.4 0.6 0.8 1 _ q i = _ q 0;i (HFEM)
In Tab. the calculation effort for a steady-state calculation of the different single stage turbine models is compared. e calculation time is provided in core-h, which represents the number of computing cores multiplied with the calculation time in hours. In contrast to the HFEM, the calculation effort for a CFD (fluid) and a CHT (fluid & solid) is shown.

e calculation effort for a CFD simulation is factor and for the CHT simulation factor higher than the simplified HFEM model (without initial file). us, the HFEM calculation does not need large-capacity computers. . CONCLUSION e operation of turbines under extraordinary conditions becomes more important in the context of flexibility. e numerical simulation of such operating conditions can help to understand the risks as well as advantages.

In this present paper, an optimized approach for the determination of the solid temperature in a steam turbine was introduced. In contrast to other approaches, a detailed Conjugate-Heat-Transfer (CHT) model was used in a first step to calculate a characteristic map of heat transfer coefficients and to adapt empirical factors of a heat transfer correlation. In a next step, a simplified solid body model was created from the CHT model by coarsening the numerical mesh and definition of relevant surface groups.

e analytical calculation approaches for convection and radiation were integrated into this simplified model. is model -so called Hybrid-FEM model (HFEM) -is able to reproduce the results of detailed CHT simulation with a accuracy of above % in less than times of the calculation time. is approach enables the simulation off-design operation -such as the warm-keeping with air -for a whole multistage steam turbine, for which the calculation effort has been exorbitant so far.

In further investigations an analytical approach for the contact heat resistance (see [ ]) will be integrated into the HFEM model. Furthermore, inter stage effects will be analyzed to develop a multistage steam turbine model. e aim is to investigate steady-state warm-keeping as well as transient pre-warming processes for a whole steam turbine with hot air.

Figure .

 . Figure . Process chart -Hybrid-FEM model

  ese main surfaces are shown in Fig. .

.Figure .

 . Figure . Hybrid-FEM model (le ), CHT-Model (right)

Figure .

 . Figure . Heat flux of the main surfaces [ ]

Figure .

 . Figure .Nusselt-number and convective heat transfer correlation -Q

  In Fig. the main surfaces and radiation heat fluxes are illustrated schematically. e blade (red) exchanges heat by radiation mainly with the neighboring blades in circumferential direction and the neighboring vanes in axial direction.

Figure .

 . Figure . Heat exchange by radiation

  Fig. and Fig. the area related heat fluxes of the main

Figure .

 . Figure . Hybrid-FEM model vs. CHT model -Q

Figure illustrates a temperatureFigure .

 . Figureillustrates a temperature contour plot of the CHT model on the le and the HFEM model on the right side. e temperature distribution is nearly identical. e highest deviation is reached at the tip of the blade and vane, due to the highest relative error of the predicted heat transfer coefficients at the shroud surfaces.

Table . r

 . RM SE of Nusselt-correlation (Eq. , )

	Surface	Q	Q
	R Blade	.	.
	R Blade Shroud	.	.
	S Vane	.	.
	S Vane Shroud	.	.
	Average	.	.
	Average (area related)

Table .

 . rRMSE of the area related heat flux

	Surface	Q	Q
	R Blade	.	.
	R Blade Shroud	.	.
	S Vane	.	.
	S Vane Shroud	.	.
	Average	.	.
	Average (area related)

Table .

 . Comparison of calculation effort (steady-state)

	Model Core-h Factor
	HFEM	.
	CHT
	CFD

  estimation of a model for prediction of lcf in steam turbine rotors. ASME Turbo Expo, . [ ] W. Traupel. ermische Turbomaschinen. Bd. . Springer, Berlin, . [ ] L. Moroz, B. Frolov, and R. Kochurov. Steam turbine rotor transient thermo-structural analysis and lifetime prediction. ASME Turbo Expo, . [ ] W.Z. Wang, P. Buhl, A. Klenk, and Y.Z. Liu. e effect of in-service steam temperature transients on the damage behavior of steam turbine rotor. International Journal of Fatique, . [ ] D. Born, P. Stein, G. Marinescu, S. Koch, and D. Schumacher. ermal modelling of an intermediate pressure steam turbine by means of conjugate heat transfer -simulation and validation. ASME Turbo Expo, . [ ] H. D. Baehr and K. Stephan. Heat and Mass Tranfer. Springer-Verlag Berlin Heidelberg, . [ ] D. Toebben, X. E. R. de Graaf, P. Luczynski, M. Wirsum, W. Mohr, and K. Helbig. Test rig for applied experimental investigations of the thermal contact resistance at the blade-rotor-connection in a steam turbine. ISROMAC, .

ACKNOWLEDGMENTS e present work is a part of the joint research program COOREFLEX-Turbo in the frame of AG Turbo. is work is supported by the German ministry of economy and technology under grant number ET B. e authors gratefully acknowledge GE Power for their support and permission to publish this paper.


Dennis Toebben
email: toebben@ikdg.rwth-aachen.de

Xavier E R De Graaf

Piotr Luczynski

Manfred Wirsum

Wolfgang Mohr

Klaus Helbig

S Y Mpo

O Tatin

Test Rig for Applied Experimental Investigations of the Thermal Contact Resistance at the Blade-Rotor-Connection in a Steam Turbine

Studies have shown that in a pre-warming respectively warm-keeping operation of a steam turbine, the blades and vanes transport most of the heat to the thick-walled casing and rotor.

ereby, a bo le-neck arises at the connection between the blade root and the rotor.

e contact heat resistance at these interfaces affects the temperature distribution and thus the thermal stresses in the rotor. e present paper introduces a thermal contact resistance test rig, which is designed to quantify the contact heat transfer at the blade-rotor-connection of a steam turbine. An uncertainty analysis is presented which proves that the average measurement uncertainties are less than one percent. In addition to the test rig, a numerical model of the specimen for the determination of the thermal contact resistance is developed and introduced. Results of several steady-state measurements under atmospheric and evacuated atmosphere using a highly temperature-resistant chromium molybdenum steel are shown. As a main influence parameter the contact pressure is investigated, which is affected by the rotational speed of the turbine. e investigations show a significant contact heat resistance especially at small contact pressures.

INTRODUCTION

Steam turbines in conventional power plants are increasingly strained due to the growing demand of flexible operation. Especially frequent and fast start-up sequences lead to high thermal stress within the thick-walled components like the rotor and the casing. e temperature gradients within the solid structures which are influenced by the heat transfer coefficients at the surfaces are responsible for these high thermal stresses. Before the steam turbine is synchronized with the grid, the heavy components need to be pre-warmed. In a conventional operation superheated steam is used for heating up the steam turbine. Other approaches use hot air, for example to keep the turbine warm or to pre-warm it. In such a operation, the main heat transfer between the main flow and the structure arises at the surfaces of the blades and vanes [ ].

e heat, being conducted to the rotor and the casing, has to pass the contact surfaces between the blade and rotor respectively the vane and the casing. ese contact surfaces constitute a kind of bo le-neck for the heat transport which affects the temperature gradients and consequently the thermal stress.

A recent study [ ] show that the contact conductance has a significant influence on the stress amplitude within the steam turbine rotor. A few experimental and analytical approaches for the determination of the thermal contact re- e experimental values provide the dependence of the TCR on roughness, material structure and the contact pressure.

e results offer an explanation for the large discrepancies that exist between the data of several authors. Burghold et al. [ ] introduce a method to measure contact heat transfer coefficients by means of transient temperature measurements. In these transient measurements an oscillation was observed from which a hysteresis effect arose.

In the present paper a method for applied experimental investigations of the TCR at the blade-rotor-connection in a steam turbine is presented. In the first part of the paper, the test rig and the numerical specimen model is introduced with a detailed explanation of the boundary conditions. Subsequently, the measurement values are discussed in an un-certainty analysis. e third part shows first measurements of the temperature at different measurement positions within the specimen by variation of the contact pressure and the heat flux. In the last part of the present paper the TCR is determined by use of the numerical model and based on the experimental results. analyze the influence of the rotational speed (i.e. contact pressure), the ambient pressure and the contact surface properties (e.g. roughness, oxidation layer, lubricant). Furthermore, the radiation influence via the air pockets, which are located between the contact surfaces can be investigated. e TCR is calculated according to Eq. ( ) in which Q c is the heat flux through the contact surface A c and ∆T c is the temperature difference between the area averaged temperatures of both surfaces which are in contact.

NOMENCLATURE

Symbols

TCR = ∆T c A c Q c = 1 α c ( )
Figures and illustrate the single components of the test rig. e key component is the traction system ( ) in which the specimen is installed. e traction system is located within a pressure chamber ( ). By evacuating the surrounding air within this chamber, the heat losses can be minimized. In addition, the pressure chamber provides the opportunity to investigate the impact of the ambient pressure on the contact heat transfer by pressurizing. e hydraulic system ( ) serves to adjust the traction force and thus the contact pressure. e peristaltic pump ( ), the air cooler ( ) and the tempered water tank ( ) belong to the cooling system. e air cooler dissipates the heat absorbed from the specimen. A tempered water tank and the peristaltic pump ensure the reproducibility of the measurement. Component ( ) is a monitored oil bath in which the reference thermocouples are located.

Figure shows the installed specimens. e experimental set-up consists of two parts. A blade specimen and a rotor specimen with its blade groove are installed in this configuration . During the test, the specimens are surrounded by an additional casing faced by high temperature insulation to reduce radiation heat losses. e blade specimen is heated up by two heating cartridges at the top in a horizontal plane. e rotor specimen is fixed at the bo om. Whereas, above the heating cartridges the traction force is transmi ed. e traction force is provided by a hydraulic stamp and measured by a load cell below the fixed baring. Each measurement signal is compensated by a reference measurement in a oil bath which is monitored by a PT-thermocouple. e traction force is measured by a calibrated load cell as mentioned before. To monitor the boundary conditions within the pressure chamber the temperature and the pressure is measured. e temperature of the water tank as well as the temperature of the heating cartridges are recorded to control the supply system. For a high accuracy the complete measurement chain in the final experimental set-up including the data logging system, ducts and compensation is calibrated. e measurement points are shown in Fig. .  e upper four (pair ) and the lower six measurement points (pair ) are used to measure the heat input and output. Due to balancing the energy flow in this way, a higher accuracy can be achieved compared to the measurement of the power input by the heating cartridges (which is affected by heat losses through the traction rod) and the temperature difference between inlet and outlet of the cooling water. e measurement points of pair to are used to determine the thermal contact resistance. All other temperature measurements are used to investigate the heat transfer through the air pockets and to validate the numerical model which is shown in

[ • C] Ambient pressure - [mbar] Heating input (max.) [W] Conductivity ( - • C) -[W/m/K] Surface roughness (R z )
.

[µm] Used material X CrMoV -Both the blade and the rotor specimen are made of the same material (Table ), which is heat resistant under high load. is material is similar to the material which is used in modern high and intermediate pressure steam turbines. e geometry of the blade has been modified to enable heat input and the tensile strength. For the design of the specimen, numerical Finite Element Methods (FEM) and Conjugate Heat Transfer (CHT) calculations were used to optimize the stability, the homogeneity of temperature and the heating up time.

e blade root design is straightened, however, the contact surfaces represent the original design.

.

NUMERICAL MODEL & ANALYSIS

ere are three reasons for building-up the numerical model of the specimen: Support in the design process, investigation of the heat transfer mechanisms between blade and rotor and in a further step the validation as well as development of a contact heat transfer correlation. During the design process the numerical model was used to design the required heating and cooling conditions and to secure the mechanical stability of the specimen by simulation of the mechanical and thermal induced stresses. In a next step, the geometry parameters were optimized. e width of the specimen conforms with the axial distance between two neighboring repetitive stages. e length of the rotor and the blade specimen was calculated to find an optimum regarding a uniform temperature distribution in the horizontal planes of pair and pair -which guarantees a minimal error in the energy balance -as well as regarding the heating up time or rather the total surface area which affects the heat losses to the ambient. Subsequently, the numerical model was used to find suitable measurement points close to the contact surface but simultaneously with a high temperature difference between opposed measurement points to minimize the relative measurement error. For that purpose, analytical approaches for the determination of the TCR from literature [ ] were integrated into the numerical model. contact areas air pockets can be found. Within theses air pockets heat transfer due to thermal convection, conduction and radiation occurs. For the investigation of the impact of thermal convection the Rayleigh-number (Ra) is calculated using Eq. ( ).

Ra = Gr • Pr = βgρ 2 | ∆T | δ 3 η 2 • ηc p λ ( )
In Eq. ( ) the expansion coefficient is evaluated as β = 1/ TAP and δ is the gap width of the air pocket between blade root and rotor groove. Based on the assumption of an average air temperature of TAP = 250 • C within the air pockets and the temperature difference between the enclosed air and the average contact surface temperature of ∆T = 50K (in a conservative estimation), the Rayleigh-number is calculated to Ra = 3.163. Because of such a low Rayleigh-number it is assumed that thermal convection can be neglected as shown in [ ]. For the calculation of the heat transfer by conduction and radiation through the air pockets, a separate mesh for these cavities is generated as shown in the right chart of Fig.

. For a higher stability of the numerical simulation, the air pockets are defined as solid bodies with the transparency and conductivity of air depending on the average air properties. Figure shows the influence of the thermal conduction and the heat radiation through the air pockets on the total heat transfer from blade to rotor ( Q total ). e calculations are conducted with constant temperature at the heat source (T HC = 500 • C) and at the heat sink (T Cool = 30 • C). With increasing TCR ( Eq. ( )) the impact of heat radiation and conduction through the air exponentially increases. With TCR → 0 the heat transfer through the air pockets ( Q AP ) related to the heat transfer through the contact interfaces ( Q c ) is about % for radiation as well as conduction and about % for only radiation (Fig.

). e effect of heat radiation grows with increasing average temperature level TAP and also with increasing temperature difference of the surfaces in contact.

During a warm-keeping or a pre-warming process, very low rotational speeds of the turbine rotor are possible, depending on the turning engine. In this case, low centrifugal forces as well as low contact pressures occur resulting in high TCRs. us, the heat radiation and the conduction through the AP have to be considered.

. UNCERTAINTY ANALYSIS

In order to quantify uncertainties of the measured data, the full measuring chain has to be analyzed. is analysis is conducted based on the "guide to the expression of uncertainty in measurement" [ ].

e uncertainty analysis can be divided into two parts. e first part includes all Type-N thermocouples mounted on the specimen. e upper thermocouples (pairs -in e lower thermocouples (pair ) are calibrated in a temperature range of 40 • C -90 • C. In accordance to [ ], the calibration is performed through a polynomial fit of th degree. Con-sequently, every thermocouple has individual parameters. After the calibration, the thermocouples are compared through a reference measurement at three temperature points. is reference measurement enables to identify temperature differences between thermocouples a er the calibration. e arrangement of thermocouples in groups with the lowest offset temperature between each other, enables to achieve a very low relative measurement uncertainty (Table ). e approach to examine the relative measurement uncertainty instead of the absolute measurement uncertainty based on [ ] is legitimate, due to the use of temperature differences for determination of the thermal contact resistance. △T R,i describes the maximum temperature difference at all reference 10 -5 10 -4 10 -3 10 -2 10 -1

T CR [m 2 K=W ] 10 -2 10 -1 10 0 _ Q total = _ Q 0;CH T [-] Adiabatic Radiation Radiation + Conduction _ Q AP = _ Q c [-]
T CR [m 2 K=W ] 10 -2 10 -1 10 0 _ Q total = _ Q 0;CH T [-] Adiabatic Radiation Radiation + Conduction _ Q AP = _ Q c [-] Figure .
Heat transfer through the air pockets related to the contact heat transfer (CHT model T HC = 500 . -. . . % e higher measurement uncertainty of the temperature level of the heating cartridge is negligible since it is only used as a security measurement to avoid the overheating of the heating cartridges. In order to achieve reproducible heating temperatures, the thermocouples in the upper level of pair are monitored. Uncertainty in overall measurement results can be induced by the measurement equipment itself, as analyzed above, but also by uncertainties of the measurement procedure.

e first series of measurements consist of operating points (Table ), varying each the ambient pressure (atm, vacuum (∼ mbar)), the traction force (p c ) and the reference temperature at upper level of pair (T H ). All measurements are steady-state. During the measurement data points were collected within a time period of five minutes and finally averaged. To analyze the impact of the ambient pressure on the heat transfer, the pressure in the surrounding atmosphere can be varied. Due to the variation of the contact pressure the effect of increasing rotational speed of the turbine rotor can be simulated and investigated. e variation of the heating temperature enables the simulation of a heating up procedure and to analyze heat radiation influence.

Table . Measurement schedule

p amb atm, vacuum T H [ • C] , , p c [kN]
. , . , . , . , . , . , . , . , .

To prove the steady state measuring conditions, the temperature difference between thermocouple T and T , as a parameter for the entering heat flux is provided in Fig. . is temperature difference fluctuates in a very narrow range up to . K, which is within the band of measurement errors induced by the thermocouples. Hence, the time averaged measurement results can be assumed as steady-state.

. RESULTS & DISCUSSION

e first series of measurements, which is presented in this section, was conducted with a simplified blade specimen. To investigate the influence of the contact pressure isolated from other influences induced by the axial contact interfaces (CI ), these contact areas were removed. Consequently, the investigated blade specimen only has the radial contact interfaces (CI ).

For the evaluation of the measurement results, a dimensionless temperature is defined in Eq. . is temperature ϑ is the average temperature difference of the radial contact surfaces ∆ TCI2 , referenced on the temperature difference at the top and the bo om of the specimen (pair -pair ). pressure p c at the contact interface CI respectively the equivalent rotational speed n (n ∝ p 0.5 c ). e measurements with atmospheric ambient pressure are represented by a continuous line and that ones with evacuated atmosphere by a dot and dash line. e different temperature levels (T H [ • C]) are characterized by separate marker types and colors.

e expected global trend is a decrease in the contact temperature difference ϑ with increasing contact pressure. Both values of ϑ (atm & vac) decrease with increasing temperature level, which indicates the heat radiation influence. At low contact pressure, the impact of the air pockets is significant as shown before in Fig. . is corresponds to the atmospheric results which show a clear temperature level dependency at low contact pressures due to conduction and heat radiation.

e impact of the heat radiation increases with both, the temperature level and the temperature difference between the opposite surfaces. However, with increasing contact pressure, the air pockets lose impact due to the predominant effect of the metal conduction. At the highest measured contact pressure, the measurement results are very close.

In contrast to the very homogeneous trend of the atmospheric measurements, the results under an evacuated atmosphere show an increased sca er. It has been observed that the temperature distribution of the specimen is very sensitive on changes in the pressure level close to vacuum. During the measurements the pressure inside the pressure chamber has fluctuated around %. e used vacuum pump is able to evacuate the pressure chamber on an absolute pressure level of mbar, but without a regulation. is can be a reason for the inhomogeneous measurement results. However, the general trend of both, atmospheric and evacuated ambient pressure, is similar. At low contact pressure the difference of ϑ between the two pressure levels is high because of the high impact of the air pockets and thus, of the ambient pressure. Concerning an evacuated atmosphere it can be assumed that the heat transfer by conduction is negligible. Hence, the remaining heat transfer phenomenon is heat radiation.

Contact pressure [MPa]

Figure shows the overall heat flux through the specimen depending on the contact pressure. e inlet heat flux, which is illustrated by a continuous line, is calculated by the temperature difference at pair and the outlet heat flux (dot and dash line) is calculated by the temperature difference at pair . Both of them are calculated with a temperature depending thermal conductivity. Atmospheric and evacuated ambient pressure measures can be differentiated by different line colors.

As a general trend it can be observed that the heat flux increases with increasing contact pressure and thus, with decreasing TCR. is effect is enhanced by the temperature level. e absolute heat losses -the difference between inlet and outlet -of the measurements with atmospheric pressure also increase with increasing temperature level. is effect can mainly be traced to the increasing heat flux, so that the relative heat losses are primary independent on the temperature level. is underlines the effect of the insulation surrounding the specimen. e heat losses at evacuated ambient pressure confirm this observation. e red lines in Fig. are nearly matching, no ma er the temperature level. At low contact pressure, the outlet heat flux is slightly above the inlet heat flux. is is within the range of measurement error.

. THERMAL CONTACT RESISTANCE

For the determination of the thermal contact resistance the numerical FEM model as well as an one dimensional approach are developed.

A FEM model with integrated thermocouple holes and tetrahedral mesh ( us, the TCR is high as well as the temperature difference between the contact surfaces.

_ Q= _ Q 0 [!]
e FEM model is able to recompute the temperature distribution in good accordance with the measurement data.

e maximal deviation averaged for the measurement series is 2 • C.

For the examination of the TCR, also a simplified onedimensional ( D) approach is used and compared with the results of the FEM model (Fig . ). e D model combined with a TCR correlation, which will be developed, enables the fast and simple calculation of the heat which is conducted from the blades to the rotor of a steam turbine. For the D approach it is assumed that the heat is transferred only in one dimension, orthogonally to the contact surface between . CONCLUSION e present paper deals with an experimental set-up for applied investigations of the thermal contact resistance (TCR) at the blade-rotor-connection in a steam turbine. erefore, a test rig was designed, which is able to measure the temperature at different measurement positions within a rotor and blade specimen. e ambient pressure can be reduced for minimal heat losses and increased measurement accuracy. For the quantification of the measurement uncertainties, an uncertainty analysis has been conducted. e relative error of the relevant temperature measurements (T 1 -T 26 ) is less than . %.

First measurements with a specimen with radial contact interfaces were conducted and introduced. e results show a distinct trend with increasing contact pressure and increasing temperature level. For the determination of the TCR an D approach as well as a numerical FEM model were developed.

e results and a comparison of both models were presented.

e TCR reaches high values with decreasing contact pressure. us, the assumption of full metal contact between the blade root and the rotor especial during operations with low rotational speed, such as start-up / cool-down or warm-keeping / pre-warming, is not justifiable.

In future investigations an analytical approach for the TCR will be developed based on the conducted and further measurements. erefore, the influence of the contact surface properties (roughness and oxidation layer) will be analyzed.

e final analytical correlation of the TCR will be integrated into a turbine model [cf. [ ]] for simulating steam turbine warm-keeping operation.

  sistance (TCR) are known from literature. Fenech, Rohsenow, Mikic, and Yovanovich are pioneers in research on TCRtopics ([ ], [ ], [ ]). e first investigations were conducted based on research on spacecra and other space vehicles. A basis for further analytical approaches provide Cooper et al. [ ]. ey developed a calculation model for conforming (nominally flat) rough surfaces. A detailed overview on the research done on this field from to is made by Madhusudana [ ]. e latest summary regarding TCR was released in in the second edition of the book thermal contact conductance wri en by [ ]. Ustinov et al. ([ ] & [ ]) developed a new model for calculating the real contact area.

Figure

  Figure shows the test rig, which is designed for the measurement of the thermal contact resistance (TCR) at the contact areas of a turbine blade. Different designs of blade or vane roots can be investigated. e test rig offers the possibility to

Figure .

 . Figure . CAD model

Figure .

 . Figure . Test rig

  Fig. . e relevant test rig and material properties are summarized in Table .

  Figure . Supply scheme

Figure .

 . Figure . Measurement positions

Figure .

 . Figure .CHT simulation with hexahedral mesh (le ), hexahedral mesh (center), contact interfaces (CI) and air pockets (AP) (right)

  Fig. ) are calibrated in a temperature range of 100 • C -500 • C.

Figure .

 . Figure . Total heat transfer from blade to rotor regarding different heat transfer phenomena through the air pockets (CHT model T HC = 500 • C, T Cool = 30 • C)

Figure .

 . Figure illustrates the measurement results as comparative temperature difference ϑ on the Y-Axis over the contact

Figure .

 . Figure .Average dimensionless temperature differences ϑ at contact interface CI

  Fig. ) is used for the numerical investi-

Figure .

 . Figure .Incoming and outgoing heat fluxes for different measurement series

Figure . Figure . Figure .

 ... Figure . Temperature distribution, calculated with FEM model (p c = 4.8 MPa, vac, T H = 150 • C)

Table .

 . Test rig and material properties

	Contact pressure (max.)	[MPa]
	Average temperature at CS	-

  Tmin,i indicates the minimum temperature difference within pair i throughout the operating range of the test rig. So, the stated relative error (△T R,i /△ Tmin,i ) in Table represents the maximum relative error of the thermocouple pairs.

	Table . Absolute and relative error of thermocouples
	mounted on the specimen
				△T R,i [K] △ Tmin,i [K] rel. err. [%]
	Pair ( -) *	.	.	.
	Pair		.	.	.
	Pair		.	.	.
	Pair		.	.	.
	Pair		.	.	.
	Pair		.	.	.
	* Difference between pair and pair .
		e second part of the uncertainty analysis includes every
	measurement equipment used to monitor the experimental
	test rig.	is equipment is also used to reproduce certain
	measurement points during series of measurements. Table
	shows that the temperature level (T HC ) can be reproduced
	in a range of . K in maximum. e reproducibility of the
	tensile force has a maximal relative error of . % and the
	PT	has an uncertainty of .	K.
	Table . Uncertainty analysis of remaining measurement
	equipment
				Measurement range U .	rel. err.
	PT	[ • C]	-	.	. %
	T HC [ • C]	-	.	. %
	p [mbar]		-	.	. %
	p c [MPa]

• C, T Cool = 30 • C) measurement points of pair i. e term △

ACKNOWLEDGMENTS

e authors gratefully acknowledge GE Power for their support and permission to publish this paper.


Bo Tan

L Scott Stephens
email: l.stephens@uky.edu

Analysis of End Face Rotary Valve Seals Including Viscoelastic Effect

Keywords:

Certain machines utilize end face rotary valves to distribute liquids from port to port. In these designs the end faces not only possess ports for distribution but also serve as end face mechanical seals to seal in the coolant. One such application is in the regenerator section of a magneto-caloric refrigeration cycle, where the cost, power usage and leakage are critical factors. Cost measures encourage the use of non-traditional seal materials such as elastomers and plastics, creating a potential for dynamic induced phenomena due to viscoelasticity. In this paper, a novel end face rotary valve seal design is introduced. The analytical model of viscoelastic effects for PTFE seals is formulated and the dynamic oscillation of the seal is measured experimentally. Finally, a brief fluid film thickness discussion is presented considering the seal ring dynamics and viscoelastic effects.

INTRODUCTION

A novel end face mechanical seal design is proposed for use in the regenerator section of a magneto-caloric refrigeration cycle in reference [START_REF] Benedict | Magneto caloric heat pump with variable magnetization[END_REF]. The whole system is shown in Figure 1 where the heat pump is 100. The exploded view of heat pump is in Figure 2. The seals are installed at 136 and 138 in Figure 2. This novel seal design combines the functions of sealing and fluid switching (valving), which simplifies the design and reduces the cost of the overall heat pump making it a potential solution to replace conventional compressor based cooling cycles in refrigeration. The potential benefits of this include less pollution, lower power usage, higher machine reliability and lower cost. However, to realize these benefits, the seal design must be optimized with respect to cost, power loss, leakage control, and seal life. One motivation of this research is to understand the operation of the seal such that an optimal design can be found. This will enable the new refrigeration system to move further towards commercialization in a product.

The key difference between this novel seal design and a conventional end face mechanical seal is that holes (ports) are machined into the face of the stationary seal ring and slots are machined into the face of the rotary seal, allowing the sealed fluid to transfer, axially, between the rotating and stationary seal rings. The primary sealing function is therefore to keep the fluid in the high pressure ports from leaking to adjacent low pressure ports and to the environment on both the outer and inner diameters of the seal rings. A second key difference is that the materials for the proposed seal faces are plastics, polymers and soft metals instead of the traditional very hard materials such as tungsten carbide and carbon. Indeed, reference [START_REF] Stieha | Analysis of an end face mechanical seal with internal valving for low power applications[END_REF] has conducted initial experimental work demonstrating that the material pairs of PTFE vs. Delrin and PTFE vs. Bronze show particular promise with respect to sealing ability and power loss. Figure 3 shows the structure of the sealing system. The upper ring, a), is the rotary face and the lower ring, b), is the stationary ring. Part c) is the sealing ring insert which provides the face material of the stationary ring. The seal insert is fixed with the stator and allows liquid to flow through the stationary slots and to be distributed to each port as the upper ring rotates. Four of the ports are supply ports operating at a pressure of about 50 psi and four of the ports are return ports operating at about 3 psi. There are 8 ports in the rotating ring and 4 distribution slots in the stationary ring. Therefore, at any one time, 2 slots are supply, 2 are return and 4 are blocked (closed). Some experiments and trials have been obtained in [START_REF] Stieha | Analysis of an end face mechanical seal with internal valving for low power applications[END_REF], which indicates PTFE vs Delrin result in low power loss and leakage rate. Though much analysis exists for conventional rotary seals, rotary valve seals are only mentioned and found in a few patents. For example, Nichols et al claim a type of fluid switching valve and its sealing faces are fluorocarboncontaining polymer and Tungsten Carbide/Carbon [START_REF] Nichols | Long lifetime fluid switching valve[END_REF]. Also, there are similar seals which are developed by Nichols [START_REF] Nichols | Dual random access, threeway rotary valve apparatus[END_REF][START_REF] Nichols | Segmenting valve[END_REF]. In the early 2010s, Moeller et al and Wan developed seal faces that allow a particular fluid to flow, but is applied as a static seal that can be rotated to a different setting [START_REF] Moeller | Fluid Switching Valve[END_REF][START_REF] Wan | Fluid selection valve[END_REF]. In the early 1990s, a seal surface on a shaft is developed by Stich, where fluid flows through the shaft and then is directed circumferentially to the outer diameter of the shaft and expelled into the static sealing area [START_REF] Stich | Rotary coupling for fluid between stationary and rotating machine parts[END_REF]. In the 1960s there was the development of a rotary distributing valve claimed by Carson et al that utilized valves to transfer the fluid stream from one conduit to any other conduits [START_REF] Frank | Rotary valve[END_REF].

In the areas of analysis and experimentation there is little or no information for this seal type. This paper attempts to fill this gap by generating an analytical model of the proposed seal including viscoelasticity and some supporting experimental work. It is essential to understand the viscoelastic behavior of PTFE used as an end face rotary valve seal, as this may affect the sealing performance. A dynamic excitation of the seal always occurs while it rotates due to the surface roughness, waviness, misalignment and other factors. Beyond a certain frequency, the seal will not be able to follow the contact motions of the seal surface and form errors, due to the viscoelastic effects of the seal material [START_REF] Stakenborg | Visco-Elastohydrodynamic (VEHD) Lubrication in Radial Lip Seals: Part 1-Steady-State Dynamic Viscoelastic Seal Behavior[END_REF][START_REF] Van Leeuwen | Visco-Elastohydrodynamic (VEHD) Lubrication in Radial Lip Seals: Part 2-Fluid Film Formation[END_REF]. Under these circumstances, the oscillating clearances may occur between the seal stationary and rotary rings. For a long term purpose, it is meaningful to find a lubrication model to describe the behavior of this type of seals so that simulation can be accomplished and advice for choosing sealing materials.

Generalized Viscoelastic Model Viscoelastic Effects

The viscoelastic effects include two basic concepts, i.e. creep and relaxation, which are shown in Figure 6. The strain response shows that it is time dependent for a viscoelastic material under a constant stress, while relaxation behavior is shown as a stress response under constant strain [START_REF] Findley | Creep and relaxation of nonlinear viscoelastic materials with an introduction to linear viscoelasticity[END_REF].

Viscoelastic Model

A Prony series representation is widely used for linear viscoelastic material functions. There are many efficient and accurate numerical methods of inter-conversion for Prony series using experimental data from selected polymer materials [START_REF] Park | Methods of interconversion between linear viscoelastic material functions. Part I-A numerical method based on Prony series[END_REF]. Here, the seal is PTFE and its mating face is Delrin. The creep compliance D is characterized using the generalized Kelvin model. It is given as [START_REF] Benedict | Magneto caloric heat pump with variable magnetization[END_REF] where Dg the glassy compliance, η 0 the zero-shear or longtime viscosity, D j retardation strengths, and τ j retardation times which are positive constants. It should be noted that η 0 ∞ for a viscoelastic solids, i.e. PTFE here. where E e the equilibrium modulus, E i relaxation strengths, and ρ i relaxation times are all positive constants. This expression is often referred as a Prony Series. Note that Equation 1 and Equation 2 can be converted.

Presmoothing Experimental Data

The constants in Equation 1 are obtained by fitting the expression from experimental data. Various methods of curve fitting have been proposed. For example, the collocation method was proposed by [START_REF] Schapery | A simple collocation method for fitting viscoelastic models to experimental data[END_REF] and the multidata method by [START_REF] Cost | A multidata method of approximate Laplace transform inversion[END_REF]. More recently, a simple and efficient Power-law presmoothing method was proposed by [START_REF] Park | Fitting Prony-series viscoelastic models with power-law presmoothing[END_REF] and numerical method of inter-conversion was also proposed by [START_REF] Park | Methods of interconversion between linear viscoelastic material functions. Part I-A numerical method based on Prony series[END_REF].

It is more efficient to fit Prony Series models by using presmoothing than direct fitting because of the computational efficiency associated with the exponential basis functions of a Prony Series. A modified power law" (MPL) is implemented, ˆ can be specified as priori. i D ˆ must be determined [START_REF] Park | Fitting Prony-series viscoelastic models with power-law presmoothing[END_REF]. The least squares method and nonlinear regression are applied for obtaining the parameters in Equation 3.

     M i n i i g t D D t D 1 ) 1 ( ) (  (3)

Prony-series Fit to Presmoothing data

Similarly, the coefficients in Equation 1 are obtained based on presmoothing curve. The curve fitting results are shown in Figure 7. The experimental data is for the PTFE ring under constant compressive stress of 6.9 MPa [START_REF] Fluoroproducts | Teflon® PTFE fluoropolymer resin: Properties Handbook[END_REF]. We use 7 terms for MPL (M=7) and n=1.5 in Equation 3. After normalization, the compliance coefficients of Prony-series fit to presmoothing data are shown in Table A numerical inter-conversion is obtained for E e and E i by use of the relationship between complex functions. D g and τ j , D j (j = 1,…,7), η 0 are known and the target constants are obtained later by a root-finding method [START_REF] Park | Methods of interconversion between linear viscoelastic material functions. Part I-A numerical method based on Prony series[END_REF], thus, other unknown sets of parameters are determined by using the complex functions relationship , ρ i is relaxation time, ω k denotes the discrete value of angular frequency. Using the rootfinding method, ρ i are shown in Table 2.

The coefficients of creep compliance are known and hence we can solve ρ i and {D} for E i . Therefore, we have obtained the coefficients of modulus E i . E(t) from Equation 2. The result is shown in Figure 8 and it is theoretically consistent. In the next section, the relaxation response in frequency domain will be discussed.

Time response of linear viscoelastic material under oscillating loading

The oscillating load is considered and it is typically caused by the surface properties of the rotary and seal faces. Then, the strain is given as follows if the stress is an input given as E'(ω) is storage modulus and E"(ω) is loss modulus. We can easily draw Relaxation modulus vs frequency in Figure 9. The loss angle δ = arctan(E"/E') = 0.001 (rad). It can also be interpreted as a measurement of the materials damping level and it is a relatively low value here. So the elastic property dominates this viscoelastic material. If other mating materials are considered in more general cases, the loss angle δ would be a critical parameter. Since the total stress σ t is written as [START_REF] Stakenborg | Visco-Elastohydrodynamic (VEHD) Lubrication in Radial Lip Seals: Part 1-Steady-State Dynamic Viscoelastic Seal Behavior[END_REF] where σ s is static stress and σ d is dynamic stress, we can determine whether a separation would happen or not due to the viscoelasticity. More specifically, assume that compressive stress is greater than 0. If total stress is less than 0, it means that separation happens. In this application, the rotating upper mount is 60±1 rpm, which corresponds to

) ( 4 . 714 ) ( * MPa E  
. Figure 11 shows experimental results for the dynamic oscillation of the seal during operation. According to this figure, the amplitude is nearly 10 2 μm at r m =25.4 mm. Thus, the oscillating load will cause 6% strain (Dynamic strain d  =d/seal thickness=99.9 μm /1562.1 μm =0.064) at r m = 25.4 (mm), then the amplitude of dynamic stress is 45.5 MPa (dynamic stress

d d E    ) ( * 
), which is greater than the static pressure 0.19 MPa (i.e. 40 lb, 1.448 in 2 , seal No.1 in Table 3). These conditions cause separation. The measurement and calculation of dynamic stress amplitude will be illustrated next.

Experimental Amplitude of Strain Oscillation

It is critical to discuss the amplitude of strain caused by oscillating load. Assume that the upper mount rotor is rigid and the lower mount stator face (PTFE) is deformed. For example, in our experiment, the rotor is Delrin and it is stiffer than PTFE [START_REF] Fluoroproducts | Teflon® PTFE fluoropolymer resin: Properties Handbook[END_REF][START_REF] Guide | Module III[END_REF]. Our test rig schematic illustration is shown in Figure 10. Torque sensor, flow-rate sensor, rpm sensor and etc. are installed and monitored by dSPACE. Two proximity sensors are located under the stator and rotor respectively. The subtraction of the amplitude of dm1 and d m2 is to determine the amplitude of dynamic displacement of the seal, defined as d, and the average surface roughness Ra corresponding to rotor upper face and contact face and PTFE contact face are checked. From the rigidity assumption for Delrin, the strain is ε=|d|/t 1 , where t 1 is the thickness of soft seal. Here the seal thickness is 1.562±0.012 mm. Figure 15 shows 6 points are measured on contact face and top face for Delrin (upper mount) and 3 points are measured on PTFE (lower mount) for surface roughness by using Zygo NV 5000. The parallelism is used as the input to the material dynamic response analysis. The parallelism, e, of rotors is measured. The surface properties are shown in Figure 13, 14 and these indicate that the oscillation is not caused by Ra because Ra << 10 2 μm in pre-run conditions. Two sets of data are provided in Table 3 Thus, the oscillation amplitude is considered as a tilted upper mount from assembling the upper mount even though a ball joint is designed between rotor and motor.

   tan )) ( cos( t r h    (11)
Thus, h(r, θ, t) without considering viscoelasticity is given as

0 h h h   
. Here, tanα is the ratio of oscillation amplitude over the distance between proximity sensor and 3 â .

t t   2 ) ( 
, because the angular velocity of the system is 2π (rad/sec). If viscoelastic effect is considered, then film thickness could be given as follow.

)] ) ( cos( )) ( [cos(            t t d h ( 12
)
where d is the amplitude of displacement of seals, δ is loss angle.

The coefficients of friction are shown in

Discussion and Conclusions

The experimental CoF and the slight wear from post-run samples show that it should be treated as mixed lubrication and contact model needs to be considered [START_REF] Minet | A deterministic mixed lubrication model for mechanical seals[END_REF][START_REF] Hamrock | Fundamentals of fluid film lubrication[END_REF]. As already discussed that viscoelastic effects would cause separations at 60 rpm, the clearance is assumed to be fully flooded. This could be considered in mixed lubrication model for future study. We can conclude the following.

Figure 1 .Figure 2 . 1 ]Figure 3 .Figure 5 .

 12135 Figure 1. The schematic illustration of an exemplary heatpump positioned in an exemplary refrigerator[START_REF] Benedict | Magneto caloric heat pump with variable magnetization[END_REF]

Figure 6 .

 6 Figure 6. Creep and Relaxation A generalized Maxwell model is given by

  M are all constants. n can be estimated as 1.5 by the glassy compliance. i

Figure 7 .Figure 8 .

 78 Figure 7. Prony Series Fit to Presmoothed Strain

  system of linear algebraic equations [A]{D} ={B}. The derivation in details are found in[START_REF] Park | Methods of interconversion between linear viscoelastic material functions. Part I-A numerical method based on Prony series[END_REF]

  phase lag. One can obtain the following familiar relationship:

Figure 9 . 1 

 91 Figure 9. Relaxation Modulus vs ω (rad/sec) The strain amplitudes are shown in Figure 11 and Figure 12 measured by SKF Ryton CMSS 665, the uncertainty of this sensor is 27 . 1  μm. The sampling frequency is 500 Hz. A low pass filter is used. The amplitude in Figure 11 and the geometric relationship in Figure 10 (b) show |d m1 | is 100 μm level at r m =25.4 mm. Similarly, the amplitude of |d m2 | at r m =25.4 mm are obtained. The period is 1 second and this is consistent with 60 rpm. |d| is close to 10 2 μm.From the rigidity assumption for Delrin, the strain is ε=|d|/t 1 , where t 1 is the thickness of soft seal. Here the seal thickness is 1.562±0.012 mm. Figure15shows 6 points are measured on contact face and top face for Delrin (upper mount) and 3 points are measured on PTFE (lower mount) for surface roughness by using Zygo NV 5000. The parallelism is used as the input to the material dynamic response analysis. The parallelism, e, of rotors is measured. The surface properties are shown in Figure13, 14 and these indicate that the oscillation is not caused by Ra because Ra

1 Figure 12 . 1 Figure 13 . 1 Figure 14 . 1 Figure 15 . 4 .

 1121131141154 Figure 11. Dynamic Oscillation Amplitude of d 1 SealNo.1

Figure 17 .

 17 Figure 17. Top View of Tilted Upper Mount Coordinate System

1 . 4 .

 14 The procedures of modeling the viscoelastic materials for this end face seal are developed, though PTFE is dominated by elastic properties. The viscoelastic effect cause separations under certain circumstances at 60 rpm. 2. Dynamic oscillation of the end face rotary seal has been found and experimental setup is illustrated for measuring the oscillation. 3. Film thickness with and without viscoelasticity are given according to the experimental setup and amplitude measurement of oscillation. A tilted upper mount has been explored and model. However, it is worth to try to avoid the tilted upper mount for further study. Mixed lubrication model is justified. The CoF and the variation of wear height indicate that Elastohydrodynamic lubrication and boundary lubrication conditions should be considered for further research.

Table 1 .

 1 . Since there is an integral relationship, Equation4, between relaxation modulus E(t) and creep compliance D(t), Creep Compliance Coefficients for Prony Fits

	(4)

Table 2 .

 2 Relaxation modulus and relaxation times

	i	τ i	ρ i	E i (GPa)
	e	-	-	1.45E-01
	1	1e+01	1.00E+01	4.83E-05
	2	1e+02	4.36E+01	4.07E-01
	3	1e+03	8.22E+02	5.15E-02
	4	1e+04	9.48E+03	1.30E-02
	5	1e+05	1.00E+05	2.71E-12
	6	1e+06	7.17E+05	6.97E-02
	7	1e+07	8.35E+06	2.77E-02

Table 3 .

 3 . Data sets for oscillating load.

	Seal No.	1	2
	Weight (kg)	18.14	36.28
	d1	max	1.631	1.791
	(mm)	min	1.200	1.400
	dm1 (μm)	113	103
	d2	max	2.342	2.173
	(mm)	min	2.337	2.168
	dm2 (μm)	1.27	1.54
	e (mm)	0.0508	0.0635
	d (μm)	99.90	86.40

Table 4

 4 T is the toque, R is radius of a point on seal face, W load is the normal load on seal face, T is measured by Interface Transducer SM-100 .The seal has slight wear (about 1 μm) both under 18.14 kg (40 lbs.) and 36.28 kg (80 lbs.). The CoF indicates that the lubrication regime is boundary and thin film lubrication. Hence, the lubrication model is justified as mixed lubrication and is viewed as a mixture of Elastohydrodynamic lubrication and boundary lubrication[START_REF] Minet | A deterministic mixed lubrication model for mechanical seals[END_REF][START_REF] Szeri | Fluid film lubrication[END_REF]. If Equation 12 is less than 0, then contact models should be than 0, then Equation12is applied in Elasto-hydrodynamic lubrication.

	by using

where

Table 4 .

 4 Coefficients of FrictionFigure 16. Tilted Upper Mount Coordinate System

Acknowledgments

Thanks to the sponsorship from GE appliances (Haier) and my colleague Mr. Joseph K. Stieha. His kindness and patience help enlighten me.