A model of impeller whirl for baffled mixing vessels

Journal of Fluids and Structures ] (]]]]) ]]]–]]]

Contents lists available at ScienceDirect

Journal of Fluids and Structures journal homepage: www.elsevier.com/locate/jfs

A model of impeller whirl for baffled mixing vessels G.C. Cudmore, A.G.L. Holloway n, A.G. Gerber University of New Brunswick, Fredericton, New Brunswick, Canada

a r t i c l e i n f o

abstract

Article history: Received 20 December 2012 Accepted 19 January 2015

A model of impeller whirl that provides stability limits and mean square radial deflections for industrial mixing vessels is introduced. The model has a linear form and includes contributions to the mass, damping, and stiffness that arise from fluid forces acting on the impeller. The fluid forces were derived from Large Eddy Simulation (LES) of a mixer with a pitched blade impeller rotating at speed N, and undergoing prescribed circular orbits of frequency Ω. Simulations were performed for Reynolds numbers, 10 r Re r 10 000, and subsynchronous whirl frequency ratios,
1 r Ω=2πN r 1. Model predictions of impeller whirl instability and mean square vibration amplitude were compared to experimental data from a laboratory scale mixer. The simplicity and economy of the new model allows its use in on-line diagnostics and operational planning. & 2015 Elsevier Ltd. All rights reserved.

Keywords: Mixing vessel Pitched blade impeller Whirl Computational fluid dynamics Linear stochastic model Fluid added coefficients

1. Introduction The objective of the present study was the development of a mathematical model for small amplitude vibration of rotating impellers that could be used as a diagnostic and operational planning tool for mixing equipment. One requirement of the model was that it could run faster than real time and could be applied to potentially thousands of operating conditions in an economical manner. A second requirement of the model was that it provides stability limits and the mean square amplitude of the random impeller orbits directly without the need for a time domain solution. These requirements effectively rule out direct, fully coupled computational fluid–structure interaction models. All of the above requirements can be met by a linear lumped-parameter approach of minimal dimensionality provided that the important effects of the fluid forces on the orbiting impeller can be accurately modelled using the concept of fluid added mass, damping and stiffness. Such an approach has been applied to the modelling of pump and turbine vibration as described by Jery et al. (1985) but it has not been applied to mixing equipment where the flow is highly turbulent and the impeller orbits are predominantly random. The method derives the fluid forces for a prescribed circular orbit of the impeller and projects them onto the impeller displacement, velocity and acceleration which are harmonic functions. In the present study Large Eddy Simulation (LES) was used to obtain the fluid forces on the orbiting impeller rather than experiment as used by Jery et al. (1985). An advantage of this modelling approach is that once the fluid added coefficients are determined for a given mixer geometry they may be used over the entire range of operating conditions without recourse to further time consuming CFD simulations. Using CFD to extract the fluid added coefficients requires that the level of CFD resolution be adequate to extract the pertinent features of the mixing system as well as efficient so that it can support simulations for each new impeller geometry. The present study provides initial guidance on these matters.

n

Corresponding author.

http://dx.doi.org/10.1016/j.jfluidstructs.2015.01.010 0889-9746/& 2015 Elsevier Ltd. All rights reserved.

Please cite this article as: Cudmore, G.C., et al., A model of impeller whirl for baffled mixing vessels. Journal of Fluids and Structures (2015), http://dx.doi.org/10.1016/j.jfluidstructs.2015.01.010i

G.C. Cudmore et al. / Journal of Fluids and Structures ] (]]]]) ]]]–]]]

2

Nomenclature Uppercase C ¼ C s þ C f total diagonal damping, N s m
1 ½C damping matrix, N s m
1 Cf fluid added diagonal damping, N s m
1 CP pressure coefficient, Cs structural diagonal damping, N s m
1 CT thrust coefficient D impeller diameter, m E modulus of elasticity, N m
2 E½r 2  expected value of mean square impeller displacement, m2 FN; FT components of fluid forces in the orbital frame, N FN ; FT average forces over one orbit period, N FN ; FT average of the orbital average fluid forces over a number of orbits, N Fp buckling load, N HðωÞ transfer function I second moment of shaft cross sectional area, m4 K ¼ K s þ K f total diagonal stiffness, N m
1 ½K stiffness matrix, N m
1 Kf fluid added diagonal stiffness, N m
1 Ks structural diagonal stiffness, N m
1 M ¼ M s þ Mf total diagonal mass, kg ½M mass matrix, kg Mf fluid added diagonal mass, kg MI impeller mass, kg Ms structural diagonal mass, kg Msh shaft mass, kg N rotational speed of impeller, rev s
1 Np power number P pressure, N m
2 P time-averaged pressure, N m
2 Re Reynolds number SðωÞ ¼ Sf ðωÞ þ Si ðωÞ power spectral density of excitation forces Sf ðωÞ power spectral density of fluid excitation forces Si ðωÞ power spectral density of rotating imbalance force S0 low frequency approximation of Sf ðωÞ impeller volume, m3 Y; Z components of impeller displacement in inertial frame, m Y_ ; Z_ components of impeller velocity in inertial frame, m s
1 Y€ ; Z€ components of impeller acceleration in inertial frame, m s
2

Lowercase c ¼ cs þcf total off-diagonal damping cf fluid added off-diagonal damping cs structural off-diagonal damping e mass eccentricity, m f frequency, Hz fX axial thrust on impeller, N f Y; f Z components of the fluid forces in the horizontal plane of the inertial frame, N g ¼9.8 gravitational acceleration, m s
2 k ¼ ks þkf total off-diagonal stiffness, N m
1 kf fluid added off-diagonal stiffness, N m
1 ks structural off-diagonal stiffness, N m
1 l shaft length, m m ¼ ms þ mf total off-diagonal mass, kg t time, s Greek symbols Ω Γ δj ϵ θ λ μ μeff μsgs ρ ω ωd ωnb

orbit rotational speed, rad s
1 mesh stiffness phase angle of the jth eigenmode, rad impeller orbit eccentricity impeller angular position, rad damping exponent dynamic viscosity, Pa s effective dynamic viscosity, Pa s dynamic viscosity (subgrid scale stresses), Pa s fluid density, kg m
3 angular frequency, rad s
1 whirl frequency of the jth eigenvalue, rad s
1 bending natural frequency of impeller and shaft in air, rad s
1

Vectors in complex plane A ¼ A þ ia eigenvector of the linear model, m C ¼ C þ ic complex damping i2 ¼
1 K ¼ K þik complex stiffness M ¼ M þim complex mass s ¼ λ þ iω eigenvalues of the linear model, s
1 r ¼ X þiY impeller displacement vector Superscript n

þ

quantity normalized using shaft rotation rate N quantity normalized using the natural bending frequency of shaft in air ωnb

A schematic of a baffled mixing vessel with pitched blade impeller is shown in Fig. 1. The role of the impeller is to circulate the fluid while the radial baffle plates enhance mixing. Fluid forces acting on the impeller arise from instantaneous asymmetric blade loading that results from fluid turbulence or impeller deflection.

Please cite this article as: Cudmore, G.C., et al., A model of impeller whirl for baffled mixing vessels. Journal of Fluids and Structures (2015), http://dx.doi.org/10.1016/j.jfluidstructs.2015.01.010i

G.C. Cudmore et al. / Journal of Fluids and Structures ] (]]]]) ]]]–]]]

3

Fig. 1. Schematic of a standard industrial mixing vessel.

2. Literature review The study of Jery et al. (1985) reports measured fluid forces and moments acting on a pump impeller undergoing a circular orbit of fixed frequency, Ω. A small orbit radius was selected to not invalidate the small perturbation assumptions that underlie a linear model, while still allowing the hydrodynamic forces resulting from whirl to be distinct and separable from the neutral position forces. The experiment was conducted with a variety of operating conditions including 3 different impellers, 7 volutes, pump speeds in a range, 500 rpm oN o 2000 rpm, and flow coefficients, 0 o C Q o 0:11. The authors expressed their force data in the orbit reference frame, FN and FT, from which added mass, damping and stiffness matrices could be extracted using a least squares quadratic fit over a range of orbit frequency ratios, Ω=2πN. A similar study to the above by Guinzberg et al. (1994) examined the effect of leakage flows in centrifugal pumps and found that increasing the flow rate in the pump reduced the destabilizing tangential force acting on the impeller. It was also observed that the tangential forces are inversely related to the clearance in the pump. Brennen and Acosta (2006) reported on similar experiments conducted on axial flow inducers and noted that the force data becomes erratic at low flows, while still resembling the results of centrifugal pumps at higher flows. This case of low flow has some similarity to the present mixer applications because of the flow turbulence and large clearances. The present work relies on CFD simulations to obtain the fluid forces acting on the orbiting impeller and it is therefore important that CFD has a demonstrated capacity to predict mixer flow fields to a reasonable accuracy. It must also have the flexibility to accommodate complex geometry and operating conditions involving stationary and rotating frames of reference, and turbulent flow fields. There have been numerous computational studies of stirred mixing vessels but here we restrict our attention to those that are most closely related to the objectives of the present study. An early study by Gotz et al. (1997) of a flat bottomed, 4 baffled mixing vessel with a six-bladed pitched blade impeller used a rotating impeller modelled using transient momentum source terms in the inertial frame. The results of these simulations were judged to be quite promising, showing a close match to experimentally determined velocity fields and the authors considered the method to be a practical engineering tool in need of further refinements. The calculations however require experimentally determined lift and drag coefficients for the impeller and this represents a serious challenge to the method in turbulent mixer flow. In another study, Brucato et al. (1998) performed a thorough examination of three different methods of modelling a mixer for CFD calculation purposes: (i) impeller boundary condition method, where the system is run as a steady state Please cite this article as: Cudmore, G.C., et al., A model of impeller whirl for baffled mixing vessels. Journal of Fluids and Structures (2015), http://dx.doi.org/10.1016/j.jfluidstructs.2015.01.010i

G.C. Cudmore et al. / Journal of Fluids and Structures ] (]]]]) ]]]–]]]

4

simulation and effects of impeller motion are added by manipulating CFD equation coefficients, (ii) inner–outer iterative procedure where a portion of the domain near the impeller is set to rotate, while the outer area remains stationary, and an overlap region is used to unify the two domains, and (iii) sliding grid method (similar to Berger et al., 2003), where a full section of the domain near the impeller is set to rotate while the baffles and outer tank remain stationary and grid techniques are applied to the interface to ensure flow continuity between the two sections. It was shown that, while the sliding grid method had the largest computation time requirements, it also provided the most reliable results. A simulation using the sliding grid method and large eddy simulation was performed by Zadghaffari et al. (2009) as well, and demonstrated very similar flow patterns to experimental results obtained with a Laser Doppler Anemometry (LDA) study. Considering the requirements of a full mixer analysis, Berger et al. (2003) performed work on a fully coupled CFD mixer simulation using an analytical model for impeller deflection and inserting this into a conventional CFD simulation with a deforming mesh and rotating frame for the impeller and a stationary frame for the baffles, separated by a sliding grid interface. The model developed for impeller deflection showed close agreement with analytical results and it was suggested that this approach significantly reduces computational time over fully coupled finite element methods. Simulations of a six-blade Rushton impeller in low Reynolds number conditions were performed by Rice et al. (2006) for Reynolds numbers of Re ¼ 1; 10; 28. The computed flow field was validated using experimental velocity data obtained using Laser Anemometry. It was found that for Reynolds numbers beyond Re ¼10 the stagnant, slug flow found at lower Reynolds numbers was replaced by a more conventional radial flow at higher Reynolds numbers. At present the application of CFD to mixers can resolve time averaged 3D flow fields with reasonable accuracy, including interactions between baffles and impellers. The sliding mesh approach with fluid structure coupling accomplished with techniques involving moving meshes within a rotating frame of reference being the most common. The present work will utilize a number of these proven approaches but will replace the two-way coupling of fluid and impeller motion with prescribed impeller orbit motion as used by Jery et al. (1985). The expressed purpose of this calculation being to generate the fluid added mass, damping and stiffness representing the forces that act on the impeller in response to the orbiting motion. This contribution to the literature serves to align the investigation of mixers to previous studies in centrifugal pumps and turbines. The only available experimental study of impeller vibration in baffled mixing vessels that could be used to evaluate the model results is that of Kippers and Holloway (2014) which describes the vibration of a pitched, four blade impeller in a baffled vessel operating over a range of rotation rates, fluid viscosities, and impeller diameters. Impeller vibration was inferred from shaft strain and high speed images of the impeller's motion. One of the findings was that the impeller displacement in the horizontal plane, in turbulent flow conditions, could be described by a rotationally symmetric normally distributed random variable. The radial variance of this distribution was found to increase exponentially with impeller rotation rate and to depend on the direction of impeller rotation. For some operating conditions there was evidence of whirl instability beyond a threshold impeller rotation rate. Results for the largest impeller were chosen for comparison in the present study because it had the smallest blade-baffle clearance and potentially the strongest blade-baffle interaction. 3. Model of mixer flow field The present simulations focus on the fluid forces acting on the impeller and this allows for some model simplifications to speed computation. For example, species transport, chemical reactions and heat transfer are not important considerations and simulations only need to model a single component fluid phase with a constant temperature. Modelling fluid flow under these assumptions requires that single-phase mass and momentum conservation equations in the form of the Navier–Stokes equations be employed, which will be discretized for solution using a finite volume approach. Additionally, a turbulence model will be needed that can accurately represent the energy containing scales responsible for the fluid forces. Note that the fine scales of turbulence responsible for mixing and chemical reaction prediction are only of indirect importance in this study. Mass and momentum conservation statements used in the finite volume method are given here for incompressible flow (ANSYS CFX, 2009b):

and

 ∂

ρ uj
wj ¼ 0; ∂xj

ð1Þ

     ∂ρui ∂

∂P ∂ ∂ui ∂uj þ ρ uj
wj ui ¼
þ μeff þ þSui ; ∂xj ∂xi ∂xj ∂t ∂xj ∂xi

ð2Þ

where ui is the velocity field in tensor form, P is the fluid pressure, μeff is an effective dynamic viscosity which can include turbulent mixing influences and Sui is a source term. At higher Reynolds numbers where turbulence is active, μeff will be equated to a sub-grid scale turbulent viscosity and the velocity and pressure field become spatially filtered quantities consistent with a Large Eddy Simulation (LES) formulation. The computational model will employ a stationary and rotating frame for baffles and impeller respectively. In the rotating subdomain, Eq. (2) above for momentum requires centrifugal and Coriolis body forces be applied and these are added to the source term Sui accordingly. Furthermore, in order to account for the moving mesh surrounding the moving impeller and shaft the formulation must solve the equations for mass and momentum conservation independent of mesh movement within each time step. This is accomplished using an Arbitrary Eulerian–Lagrangian (ALE) approach (Hirt et al., 1974) Please cite this article as: Cudmore, G.C., et al., A model of impeller whirl for baffled mixing vessels. Journal of Fluids and Structures (2015), http://dx.doi.org/10.1016/j.jfluidstructs.2015.01.010i

G.C. Cudmore et al. / Journal of Fluids and Structures ] (]]]]) ]]]–]]]

5

where all convective transport is relative to the mesh velocity, wj. In the ALE approach a conservation equation for space is needed to ensure the change in control volume size over time, when a moving mesh is active, is balanced with the volume change resulting from control-volume surface motion, Z Z d dV ¼ wj dnj ; ð3Þ dt VðtÞ S where the mesh velocity field wj is derived by solving the Laplacian equations for displacement of the nodal coordinate Δxi ,   ∂ ∂Δxi Γ ¼ 0: ð4Þ ∂xj ∂xj The mesh stiffness, Γ, is chosen to ensure mesh quality and the impeller orbit provides boundary conditions to Eq. (4). The LES method takes advantage of the universal, isotropic properties of smaller scale eddies in turbulent flows and models them while directly resolving the larger energy containing eddies. The constant mesh density chosen for the mixer model (as discussed subsequently) is useful for supplying a consistent filter function for the LES method. Thus the LES method will resolve the larger scale eddies in the bulk of the mixing vessel and use subgrid-scale model for μsgs to estimate the subgrid scale stresses. In boundary layer regions an LES WALE (wall-adapted local eddy-viscosity) model was used due to its improved performance (Nicoud and Ducros, 1999). The WALE model is also computationally efficient in complex geometries requiring no secondary filtering or application of damping functions near walls. It is expected that the flow separation from the impeller and baffles will be determined more by the geometric features than the surface shear stress distributions as the boundary layer modelling used is considered to be quite coarse.

4. Model of impeller whirl The equations that describe small dynamic motions of the impeller in the horizontal plane can be written according to Adams (2001) as )  ( € )  ( _ )   ( f Y ðtÞ M
m C
c Y K
k Y Y ; ð5Þ þ þ ¼ f Z ðtÞ m M c C Z k K Z_ Z€ T

where ½M, ½C and ½K are the mass, damping and stiffness matrices respectively, and ff Y ðtÞ; f Z ðtÞg is the forcing vector. Eq. (5) for the impeller motion displays a unique anti-symmetric character that allows reduction to a single ordinary differential equation for the complex vector displacement, r ¼ Y þ iZ. Following Vance (1988) we rewrite Eq. (5) in complex form as Mr€ þ Cr_ þ Kr ¼ fðtÞ;

ð6Þ

with the complex quantities: M ¼ M þ im, C ¼ C þ ic, K ¼ K þik. A limitation of this method is that the excitation force, fðtÞ ¼ f Y ðtÞ þ if Z ðtÞ, must be a rotating phasor with a phase difference of either 7 901 or a random variable which is statistically axisymmetric in the (Y, Z) plane as described by Fuhrmann (1999). The diagonal mass term M is given by M ¼ Mf þ Ms ;

ð7Þ

where Mf is the fluid added mass and Ms is the structural mass which includes the mass of the impeller, MI, and a representative portion of the shaft mass, Msh. For a cantilever shaft we use M s ¼ M I þ 12 M sh :

ð8Þ

The off-diagonal component of mass, m, will be taken as zero, a standard assumption in rotordynamic analysis according to Adams (2001), because an off-diagonal mass causes unrealistic instability at a high rate of rotation relative to the natural frequency. Consistent with this assumption the results of the CFD simulations presented in a later section found that the offdiagonal mass, mf, is essentially zero. The diagonal damping term C is given by C ¼ Cf þ Cs;

ð9Þ

where Cf is the diagonal fluid added damping. The structural damping, Cs, was determined from the logarithmic decrement of a bump test of the shaft in air. It depends on shaft material, size and bearing support. The off-diagonal damping is entirely due to fluid forces so that c ¼ cf :

ð10Þ

The diagonal stiffness, K, combines fluid and structural components K ¼ Kf þ Ks;

ð11Þ

Please cite this article as: Cudmore, G.C., et al., A model of impeller whirl for baffled mixing vessels. Journal of Fluids and Structures (2015), http://dx.doi.org/10.1016/j.jfluidstructs.2015.01.010i

G.C. Cudmore et al. / Journal of Fluids and Structures ] (]]]]) ]]]–]]]

6

where Kf is the diagonal fluid added stiffness. The structural stiffness of the shaft, Ks, is calculated according to Berger et al. (2003) as   3EI f Ks ¼ 2 1 þ x ; ð12Þ FP l where l is the shaft length, I is the second moment of area of the shaft cross-section, and E is the modulus of elasticity of the shaft. The buckling load, FP, can be calculated as FP ¼

π 2 EI 4l

2

:

ð13Þ

The net thrust on the impeller, fx, is given by f X ¼ C T ρD4 ð2πNÞ2 þ M I g þ M sh g
ρ g;

ð14Þ

where CT is the thrust coefficient for a rotating pitched blade impeller, ρ is the fluid density, D is the impeller diameter, g is the acceleration due to gravity, and is the volume of the shaft and impeller. The sign of CT changes with the direction of shaft rotation such that axial load is compressive for N 4 0 and tensile for N o 0. The off-diagonal stiffness, k, includes both fluid and structural contributions k ¼ kf þks ;

ð15Þ

where ks describes the follower force caused by internal shaft friction (Vance, 1988) ks ¼ 2πNC s :

ð16Þ

The fluid added mass, damping and stiffness, Mf, Cf, cf, Kf, and kf, will be determined from LES simulation of the fluid flow in the mixing vessel with an orbiting impeller. The solution of Eq. (6) has the form r ¼ Aest where A ¼ A þia. The eigenvalues have the form sj ¼ λj þiωdj where λj are the damping exponents and ωdj are the whirl frequencies (j¼1, 2, 3, 4). Specifically pffiffiffi pffiffiffi
C 7 R cos ðγ=2Þ
c 7 R sin ðγ=2Þ 7i ; ð17Þ s¼ 2M 2M where tan γ ¼

2Cc
4Mk C 2
c2
4MK

;

ð18Þ

and qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð19Þ ðC 2
c2
4MKÞ2 þ ð2Cc
4MkÞ2 : pffiffiffi The same sign for 7 R is selected for each term. It should be noted that M, C, c, K, and k are functions of shaft rotation rate, impeller shape and diameter, fluid viscosity and density. In addition, the off-diagonal stiffness and damping, k and c, reverse sign with a reversal of shaft rotation rate. In the absence of off-diagonal damping, c, Eq. (17) yields a result equivalent to analytical solutions presented by Vance (1988). The corresponding orbit modes can be written in terms of the eigenvector magnitude, jAj j, and phase, δj, as R¼

Y j ¼ eλj t jAj j cos ðωdj t þ δj Þ;

ð20Þ

Z j ¼ eλj t jaj j cos ðωdj tÞ:

ð21Þ

and

The magnitude of the transfer function, jHðωÞj, for the present linear system is given by jHðωÞj2 ¼

1 ð
Mω2
cω þ KÞ2 þ ðCω þkÞ2

:

ð22Þ

It is asymmetric with respect to the direction of whirl, ω, because of the imaginary stiffness and damping, k and c. The mean square radial displacement of the impeller can be calculated according to Newland (1975) using the transfer function and the power spectral density of the excitation force Z 1 E½r 2  ¼ jHðωÞj2 SðωÞ dω; ð23Þ
1

where SðωÞ ¼ FðωÞFn ðωÞ and Fn ðωÞ is the conjugate of the Fourier transform of the excitation forces fðtÞ arising from rotating imbalance and fluid excitation forces. The power spectral density of the rotating impeller imbalance can be expressed as Si ¼ ðM s eω2 Þ2 δðω þ2πNÞ;

ð24Þ

Please cite this article as: Cudmore, G.C., et al., A model of impeller whirl for baffled mixing vessels. Journal of Fluids and Structures (2015), http://dx.doi.org/10.1016/j.jfluidstructs.2015.01.010i

G.C. Cudmore et al. / Journal of Fluids and Structures ] (]]]]) ]]]–]]]

7

where e is the mass eccentricity and δðωþ 2πNÞ ¼ 1 when ω ¼ 2πN and zero otherwise. The nature of the imbalance force is such that it has the same frequency as the impeller rotation and it rotates in the same direction as the impeller. The mass eccentricity is determined using the mass of the impeller and the mean geometric shaft eccentricity. The power spectral density for the fluid excitation forces in turbulent flow has a broad band random character that was derived from LES simulation of the fluid flow in the mixing tank with a rotating impeller. The low frequency portion of SðωÞ is particularly important because the natural frequency of the impeller vibration is approximately 4 Hz. In this range, SðωÞ is approximately constant but is subject to large statistical uncertainty which arises from the limited number of long period events that occur within a simulated time record. For both these reasons, Sf ðωÞ will be approximated for the purpose of estimating the mean square impeller motion as Sf ðωÞ ¼ S0 ;

ð25Þ

where S0 depends on the fluid properties, impeller rotation rate and impeller and tank geometry. An important characteristic of the fluid excitation forces is that the phase relationship between fY and fZ spans a wide range, unlike the imbalance force, which is radially directed with a phase angle between fY and fZ of þ 901. Substituting the above power spectral density estimates into Eq. (23) the mean square radial displacement of the impeller can be calculated accordingly as Z 1 jHðωÞj2 dω: ð26Þ E½r 2  ¼ ðmeÞ2 ð2πNÞ4 jHð2πNÞj2 þSo
1

5. CFD results 5.1. Computational model The details of the mixing vessel geometry and operation required the creation of two separate domains to be joined: a rotating inner domain where a moving ALE mesh relative to the impeller is active, and an outer domain with a stationary mesh. Flow between the two domains is handled using a sliding grid interface. This allows flow between the rotating and stationary sections without requiring an exact node to node match between the two domains, and this is crucial for the transient runs of the simulation. This configuration is consistent with a transient rotor–stator sliding grid method often employed in turbomachinery applications. The outer domain consists of the outer tank walls, including a ring portion of the bottom of the tank and a ring portion of the lid, as well as the four baffles. All of these walls are treated as no-slip walls in the simulation. The inner domain includes the impeller shaft, the impeller, the remainder of the top and bottom walls of the tank and the surrounding fluid. The impeller and shaft as well as the lid and bottom of the tank are all treated as no-slip walls. The shaft and impeller are assigned to rotate at a set rate and orbit (which requires the ALE model be activated) depending on the particular simulation. The inner domain was configured to rotate 1/256th of a rotation every time step to ensure that Fourier analysis of complete rotations are computationally efficient. As a result, the physical time covered by a time step is a function of the shaft rotation frequency. Each simulation was begun as a steady state simulation using a frozen rotor model (baffles and impellers fixed relative to one another but still solved in different frames of reference) to obtain initial conditions for the fully transient runs with impeller orbit. These steady state simulations are run until all momentum and mass residuals are below an RMS value of 1E
4. Subsequently transient simulations are started using the transient rotor–stator model. This model fully accounts for the interaction between the two domains at every time step as the inner domain is clocked relative to the outer. A maximum of 10 coefficient loops are performed within each time step, once again targeting RMS residual values below 1E
4. At each time step data pertaining to pressures and forces were calculated and written out in the form of normal and tangential forces acting on the orbiting impeller. The CFD solution was obtained using an unstructured isotropic tetrahedral mesh to provide a consistent filter for the LES model. The meshing can be visualized via the major surface elements in Figs. 2 and 3. ANSYS CFX (2009a) was used with a nominally second order accurate finite volume discretization on the mesh which is subsequently solved in a coupled manner using an algebraic multi-grid to maintain linear solution time with increasing mesh sizes. Using CFD to generate the normal and tangential forces under small amplitude impeller displacements and then extracting fluid added mass, stiffness and damping for use in a low-order rotordynamic model suggests some lee-way in the resolution of the CFD model. What is sought is essentially aggregated information on the fluid dynamics. It should be noted that obtaining statistically converged force data from an unsteady mixing process can be computationally very expensive if the mesh resolution is unnecessarily large. To determine an appropriate mesh resolution the integrated forces on the impeller were evaluated for mesh sizes of 270k, 870k and 1.7M nodes at Reynolds numbers of 10 and 1000, and orbit frequencies in the range,
1 o Ω=2πN o1. Note that these node counts are equivalent to 1.35M, 4.35M and 8.5M elements respectively in other popular CFD software because CFX forms its control volumes around element vertices (or nodes) as opposed to using elements for the controlvolumes. Results of this study suggested that the percent change in the impeller forces (from the smallest mesh to largest) was found to be in the range of 10%. In addition, it was found that the use of the 270k node mesh, over a 1.7M mesh, allowed for unsteady simulations to be statistically converged in 2 orders of magnitude less CPU time and therefore provided the Please cite this article as: Cudmore, G.C., et al., A model of impeller whirl for baffled mixing vessels. Journal of Fluids and Structures (2015), http://dx.doi.org/10.1016/j.jfluidstructs.2015.01.010i

G.C. Cudmore et al. / Journal of Fluids and Structures ] (]]]]) ]]]–]]]

8

Fig. 2. Surface meshes for the impeller and shaft.

opportunity to reduce the statistical error of the fluid added coefficients which is relatively large for time average quantities in turbulent mixer flow. Considering both mesh resolution errors and statistical errors a 270k mesh was considered a practical compromise to generate the fluid added coefficients over a wide range of frequency ratios and Reynolds number that is required in this work. The details of the mesh independence study are provided in Section 5.7 following the introduction of the fluid forces, model coefficients, and statistical treatment of the data. The subsequent success in applying the fluid added coefficients calculated from this approach to experimental data, as described in the next sections, suggests that high resolution simulations may not be needed to deduce the forces acting on the impeller. This will be a matter that will be studied in more detail in future research. 5.2. Pressure and velocity fields All of the CFD simulations were performed on the model mixer geometry used in the experiments of Kippers and Holloway (2014). The mixing tank had 4 equally spaced baffles running the full length of the tank and the impeller had a 254 mm diameter and 451 blade pitch. In the computations, the impeller rotation rate was fixed at N¼ 269 rpm, the fluid density ρ was 998 kg/m3, and the fluid viscosity μ was varied to generate a range of Reynolds number, 10 rρND2 =μ r10 000. The CFD simulations used a prescribed circular impeller orbit of frequency, Ω, so the impeller mass and shaft stiffness do not have any role in the simulations or interpretation of the results. A rendering of the instantaneous velocity field on a horizontal plane near the middle of the impeller for Re ¼10 000 is shown in Fig. 4. Overall, the velocity field is disordered with significant velocity vectors in the baffled outer region of the vessel and small vortices forming near the baffles. There is a trend toward outward radial flow which would need to be balanced by inflow from planes above and below it. The impeller is eccentric with the blade-baffle clearance being smaller near the bottom of the image. The velocity field on a meridional plane is shown in Fig. 5 at an instant when the blades are midway between the baffles. The velocities are greatest in the vicinity of the blades and a tip vortex is evident on left hand side where the blade is moving away from the reader. Please cite this article as: Cudmore, G.C., et al., A model of impeller whirl for baffled mixing vessels. Journal of Fluids and Structures (2015), http://dx.doi.org/10.1016/j.jfluidstructs.2015.01.010i

G.C. Cudmore et al. / Journal of Fluids and Structures ] (]]]]) ]]]–]]]

9

Fig. 3. Surface meshes for full tank.

Fig. 4. Computed instantaneous velocity vector field for a plane located near the center of the impeller, as viewed from the top of the vessel for Re¼10 000. The impeller is rotating clockwise.

Fig. 6 shows a view of the pressure field from the same perspective as Fig. 4. Areas of high pressure can clearly be seen on the front of the baffles as the blades approach, in CW rotation, and areas of low pressure on the back of the baffles. The lowest pressures occur at the bottom of the figure where the blade-baffle clearance is least and a vortex core trails from the Please cite this article as: Cudmore, G.C., et al., A model of impeller whirl for baffled mixing vessels. Journal of Fluids and Structures (2015), http://dx.doi.org/10.1016/j.jfluidstructs.2015.01.010i

G.C. Cudmore et al. / Journal of Fluids and Structures ] (]]]]) ]]]–]]]

10

Fig. 5. Computed instantaneous velocity vector field for a meridional plane passing through the center of the vessel, as viewed from the side of the vessel for Re¼ 10 000.

Fig. 6. Computed instantaneous pressure field for a plane located near the center of the impeller, as viewed from the top of the vessel for Re ¼10 000. The impeller is rotating clockwise.

blade tip. High pressures can be seen on the blade surfaces with the highest being on the right side. This asymmetry of the pressure field is responsible for the net force on the impeller in the horizontal plane and is to be expected because the flow currents within the tank generate different angles of flow incidence to each of the blades. The pressure field in a meridional plane for this case is shown in Fig. 7. Overall the pressure below the impeller is highest, which is consistent with the flow being pushed down in the center of the tank and returning upwards between the baffles. This pressure difference occurs across the impeller but is nonuniform with the greatest pressure difference across the right hand blade suggesting that the lift on this blade is greatest. 5.3. Experimental validation of CFD results The power number, _ W Np ¼ 5 3 ; ρD N

ð27Þ

and coefficient of thrust, CT ¼

fX ρð2πNÞ2 D4

;

ð28Þ

Please cite this article as: Cudmore, G.C., et al., A model of impeller whirl for baffled mixing vessels. Journal of Fluids and Structures (2015), http://dx.doi.org/10.1016/j.jfluidstructs.2015.01.010i

G.C. Cudmore et al. / Journal of Fluids and Structures ] (]]]]) ]]]–]]]

11

Fig. 7. Computed instantaneous pressure field for an XZ plane passing through the center of the vessel, as viewed from the side of the vessel for Re¼10 000.

Table 1 Computed power numbers and thrust coefficient. Re

NP

CT

10 100 1000 10 000

5.9 1.6 1.6 1.4

0.020 0.012 0.015 0.016

were derived from CFD simulations for different Reynolds numbers and the results are summarized in Table 1. The computed power numbers are practically independent of Reynolds number above Re ¼100, as indicated in the literature reviewed by Oldshue (1983), and the value at high Reynolds number compares well with NP ¼1.5, measured by Kippers and Holloway (2014). Pressure fluctuations were made by Kippers and Holloway (2014) at the midpoint of the mixing tank wall 0.14 m from the bottom of the mixing vessel and in 3 angular positions: 2π=45 rad ahead of the baffle, 2π=45 rad behind the baffle, and at π=4 rad from the baffle. Data was sampled over 135 revolutions of the impeller at a rate of 250 Hz. Results for the pressure coefficient, P
P CP ¼ 1 2 ; 2 2ρD ð2πNÞ

ð29Þ

where P is the time averaged pressure, are provided in Figs. 8 and 9. The values of CP, indicated in these polar plots correspond to the measured pressure at θ ¼ 0 when one of the 4 impeller blades are at the angular position, θ ¼ 2πNt. The impeller is rotating CCW in this diagram. For example, in Fig. 8, when the impeller blade is located at 11π=6 ahead of the baffle the pressure at the transducer location, θ ¼ 0, is a maximum. When the impeller blade is directly across from the baffle the pressure at θ ¼ 0 is a minimum. Figs. 8 and 9 all show 4 fold symmetry because there are 4 equally spaced blades. The pressure values at the transducer position were derived from CFD simulations run at 224 rpm for 40 revolutions without impeller eccentricity using a 270k node mesh. Figs. 8 and 9 show qualitative agreement between the CFD and experimental results with differences in the maxima of C P o 20%. However, Fig. 10 shows the CFD results significantly overpredict the maximum values of CP and underpredict the minimum values that occur immediately behind the baffle as the blade passes. In evaluating this comparison one should consider that the absence of impeller eccentricity in this simulation would affect the blade baffle clearance. 5.4. Calculating the fluid forces The fluid forces acting on the orbiting impeller can be described as ( ) " # f Y ðtÞ BYY BYZ Y ¼ ; f Z ðtÞ BZY BZZ Z

ð30Þ

Please cite this article as: Cudmore, G.C., et al., A model of impeller whirl for baffled mixing vessels. Journal of Fluids and Structures (2015), http://dx.doi.org/10.1016/j.jfluidstructs.2015.01.010i

G.C. Cudmore et al. / Journal of Fluids and Structures ] (]]]]) ]]]–]]]

12

Fig. 8. Comparison of tank wall pressures at θ ¼ 0 from experiment and CFD simulation as a function of impeller angular position. The baffles are located 2π=45 rad after the pressure transducer and Re ¼ 3200. Solid line represents experimentally derived CP, dashed lines are the 95% confidence interval for that data, and solid circles are CFD results for CP. The impeller is rotating counter-clockwise.

Fig. 9. Comparison of tank wall pressures at θ ¼ 0 from experiment and CFD simulation as a function of impeller angular position. The baffles are located π=4 rad before the pressure transducer and Re ¼ 3200. Symbols as in Fig. 8.

where the matrix ½B may be decomposed into components proportional to the impeller displacement, velocity and acceleration (Jery et al., 1985) (

BYY

BYZ

BZY

BZZ

)

Y Z

"

 ¼


Mf


mf

mf

Mf

#(

Y€ Z€

)

Please cite this article as: Cudmore, G.C., et al., A model of impeller whirl for baffled mixing vessels. Journal of Fluids and Structures (2015), http://dx.doi.org/10.1016/j.jfluidstructs.2015.01.010i

G.C. Cudmore et al. / Journal of Fluids and Structures ] (]]]]) ]]]–]]]

13

Fig. 10. Comparison of tank wall pressures at θ ¼ 0 from experiment and CFD simulation as a function of impeller angular position. The baffles are located 2π=45 rad before the pressure transducer and Re¼3200. Symbols as in Fig. 8.

Fig. 11. Definition of the circular impeller orbit and orbital (N, T) reference frame. The impeller center lies on the circle at the origin of the orbital frame. Y and Z are coordinate directions in the inertial reference frame shown here at the center of the mixing vessel.

"


Cf


cf

cf

Cf

#(

Y_ Z_

)

"


Kf kf

#
kf  Y : Kf Z

ð31Þ

The prescribed impeller rotation and orbital motion are described using three reference frames: (1) the inertial reference frame, (2) the rotating reference frame attached to the impeller, and (3) the orbital reference frame associated with the impeller orbit. The inertial reference frame remains fixed in its orientation to the baffles. The rotating reference frame rotates relative to the inertial frame at the shaft rotation rate, N. The impeller center is offset by eccentricity ϵ from the geometric center of the tank and for the CFD simulation it is forced to travel in a circular path with an orbit frequency Ω, while the impeller rotates at rate 2πN. The orbital reference frame has coordinate directions normal and tangential to the circle of the orbit at the instantaneous position of the impeller, as shown in Fig. 11. The inertial coordinates of the impeller Please cite this article as: Cudmore, G.C., et al., A model of impeller whirl for baffled mixing vessels. Journal of Fluids and Structures (2015), http://dx.doi.org/10.1016/j.jfluidstructs.2015.01.010i

G.C. Cudmore et al. / Journal of Fluids and Structures ] (]]]]) ]]]–]]]

14

center are defined as Y ¼ ϵ sin Ωt;

ð32Þ

Z ¼ ϵ cos Ωt:

ð33Þ

and

The transformation used to convert the force data from the inertial reference frame (Y–Z) to the orbital reference frame (N–T) has the form F N ¼ f Y sin Ωt þ f Z cos Ωt;

ð34Þ

F T ¼ f Z sin Ωt
f Y cos Ωt:

ð35Þ

Inverting the above transform and substituting with Eq. (31) into Eq. (30) give " # " #( ) M f
mf F N
F T sin Ωt Y€ ¼
mf Mf FT FN cos Ωt Z€ " #( ) " # C f
cf K f
kf Y Y_ :

cf Cf kf Kf Z Z_

ð36Þ

Finally, substituting Eqs. (32) and (33) and their derivatives into Eq. (36), the forces in the normal and tangential reference frame can be expressed as F N ¼ ϵðM f Ω2
cf Ω
K f Þ;

ð37Þ

F T ¼ ϵð
mf Ω2 þC f Ω
kf Þ:

ð38Þ

In this study the fluid added mass, damping and stiffness, Mf, Cf, cf, Kf, and kf, will be obtained as the least squares quadratic fit of the forces acting on the impeller, in the normal–tangential orbital reference frame, over a range of orbit frequencies, Ω. The fact that the fluid added stiffness, damping, and mass are proportional to the 0th, 1st, and 2nd powers of Ω is the basis of this composition method. 5.5. Treatment of the unsteady force data The process of determining the fluid added mass, damping and stiffness from Eqs. (37) and (38) requires the determination of average values for FN and FT over an orbit period. These average force components were derived from CFD simulations conducted over the range of Reynolds numbers and whirl frequencies shown in Table 2. Before force records were considered for analysis a sufficient number of shaft rotations were simulated to ensure the decay of any startup transients in the solution. Following this each simulation was extended until the running average of the force components became independent of the number of orbits included. For the Re ¼10 cases this required only a single orbit period but for the high Reynolds number turbulent flows extended simulations and statistical analysis were required. In either case, the average normal force component over a single orbit was calculated as ! n X FN ¼ F Ni =n; ð39Þ i¼1

where n is the number of data points per orbit for the case being considered, which may range from 256 points to 1024 points. For the random cases the average force value F N was then averaged over a number of orbit periods, that is 0 FN ¼ @

m X

1 F Nj A=m;

ð40Þ

j¼1

Table 2 Parameters used for numerical simulations. Ω=ð2πNÞ Re 10 100 1000 10 000


1.25


1


0.5

✓

✓ ✓ ✓ ✓

✓ ✓ ✓ ✓


0.25

✓

0.25

0.5

1

✓

✓ ✓ ✓ ✓

✓ ✓ ✓ ✓

Please cite this article as: Cudmore, G.C., et al., A model of impeller whirl for baffled mixing vessels. Journal of Fluids and Structures (2015), http://dx.doi.org/10.1016/j.jfluidstructs.2015.01.010i

G.C. Cudmore et al. / Journal of Fluids and Structures ] (]]]]) ]]]–]]]

15

Fig. 12. Normal and tangential forces acting on the impeller in the orbital frame of reference with Ω=ð2πNÞ ¼
0:5 and Re¼ 10. The results are plotted for one full orbit, equivalent to 2 full rotations of the shaft. Mean value indicated by bold line, standard deviation by dashed lines.

Fig. 13. Normal and tangential force components acting on the impeller in the orbital frame of reference with Ω=ð2πNÞ ¼
0:5 and Re ¼ 10 000. The data record includes 12 full orbits, equivalent to 24 rotations of the shaft. Mean value indicated by bold line, standard deviation by dashed lines.

where m is the number of orbit periods being analyzed. Convergence of the averages is tested using a difference between subsequent averages of the orbit forces, as ΔF ¼ jF Nj
F Nðj
1Þ j:

ð41Þ

The average is considered converged when ΔF o 5:0  10
4 N. This value was chosen to ensure relatively small statistical effects on the calculated coefficients without necessitating excessive computation time. To obtain statistical convergence in the case with Re ¼10 000 and Ω=ð2πNÞ ¼
0:5 required 12 full orbits and 24 rotations of the shaft. In the high Reynolds number turbulent cases more frequency ratios were simulated to increase the confidence in the curve fitting of FN and FT versus Ω and hence in the fluid added mass, damping and stiffness. A sample of the fluid forces in the orbit frame for the Re ¼10 and Ω=ð2πNÞ ¼
0:5 case expressed in non-dimensional form, F nN ¼ F N =ρD2 ð2πNÞ2 , is shown in Fig. 12. The fluctuations of the force signal correspond to the interactions between the blades and baffles which occur 4 times per shaft revolution. Fig. 13 shows a sample of data from the Re¼10 000 and Ω=ð2πNÞ ¼
0:5 case in the orbit frame for comparison. The flow is fully turbulent and the periodicity observed in the Re¼10 case is much less pronounced. The frequency content of the force for the Ω=ð2πNÞ ¼ 7 0:5 case at Reynolds numbers of 10 and 10 000 is shown in Fig. 14. At each Reynolds number there is a strong peak at a frequency of 14.08 rad/s which corresponds to the orbit frequency jΩj ¼ 0:5ð2πNÞ. In order to generalize the power spectral density of the fluid excitation forces we express it in the

Please cite this article as: Cudmore, G.C., et al., A model of impeller whirl for baffled mixing vessels. Journal of Fluids and Structures (2015), http://dx.doi.org/10.1016/j.jfluidstructs.2015.01.010i

G.C. Cudmore et al. / Journal of Fluids and Structures ] (]]]]) ]]]–]]]

16

Fig. 14. Spectral analysis of the computed forces acting on the impeller for Re ¼ 10 (A, B) and 10 000 (C, D), Ω=ð2πNÞ ¼
0:5 ðA; CÞ and Ω=ð2πNÞ ¼ þ 0:5 ðB; DÞ.

non-dimensional form
 Snf ωn ¼

n

Sf ðω Þ ; ρ2 ð2πNÞ3 D8

ð42Þ

where the angular frequency is ωn ¼ ω=2πN. In the absence of impeller orbit the dimensionless power spectral density Snf ¼ f ðReÞ, and one would further expect that at high Reynolds numbers, Snf was independent of Reynolds number. Using Fig. 14 we infer that the low frequency estimate, Sno ¼ 5  10
7 , and so the dimensional power spectral density in the low frequency range would be S0 ¼ ρ2 ð2πNÞ3 D8 Sno ;

ð43Þ

where the shaft rotation rate 2πN o 100 rads=s (N o 16 Hz).

5.6. Calculating the fluid added coefficients The fluid added mass, damping and stiffness coefficients are defined as M nf ¼

Mf ρD3

C nf ; cnf ¼

; C f ; cf

ρD3 2πN

ð44Þ

;

ð45Þ

Please cite this article as: Cudmore, G.C., et al., A model of impeller whirl for baffled mixing vessels. Journal of Fluids and Structures (2015), http://dx.doi.org/10.1016/j.jfluidstructs.2015.01.010i

G.C. Cudmore et al. / Journal of Fluids and Structures ] (]]]]) ]]]–]]]

17

Fig. 15. Plots of average normal force acting on the impeller as a function of whirl frequency and Reynolds number:    Re ¼1000; – – Re¼10 000.

Re¼ 10; - - - Re¼100;

Fig. 16. Plots of average tangential force acting on the impeller as a function of whirl frequency and Reynolds number:    Re¼1000; – – Re¼10 000.

Re¼10; - - - Re¼ 100;

and n

K nf ; kf ¼

K f ; kf ρD3 ð2πNÞ2

:

ð46Þ

The normal and tangential forces on the impeller, Eqs. (37) and (38), are re-written in terms of these coefficients as F nN ¼ ϵn ½M nf ðΩ=2πNÞ2
cnf ðΩ=2πNÞ
K nf ;

ð47Þ

and n

F nT ¼ ϵn ½
mnf ðΩ=2πNÞ2 þ C nf ðΩ=2πNÞ
kf ;

ð48Þ

where the normalized impeller eccentricity is ϵn ¼ ϵ=D. Plots of the average normal and tangential force versus orbit frequency ratio, Ω=ð2πNÞ, created using the CFD simulation data at Reynolds numbers of Re¼10, 100, 1000 and 10 000 are shown in Figs. 15 and 16. A quadratic fit in terms of Ωn was applied as described in Section 5.4 to obtain the coefficients representing the added mass, damping and stiffness coefficients. The normal force trends displayed in Fig. 15 for increasing Re include a decreasing curvature of the quadratic fit (proportional to the added mass, M nf ), a decrease in the slope at zero whirl (proportional to the off-diagonal damping coefficient cnf ), and a nearly constant intercept at zero whirl resulting in a fixed stiffness coefficient, K nf 
0:05. The tangential force trends in Fig. 16 for increasing Re include zero curvature for all Re (implying that the off-diagonal mass mf ¼0, consistent with the maxim that off-diagonal mass must be zero), a nearly constant slope and intercept for Re 4100 n implying that the diagonal damping coefficient, C nf  0:1, and the off-diagonal stiffness coefficient, kf  0:08, is nearly independent of Reynolds number. Table 3 provides a summary of the fluid force coefficients deduced from Figs. 15 and 16. These coefficients can be used to determine the fluid added mass, damping and stiffness for geometrically similar impellers rotating at different speeds and in Newtonian fluids having different values of ρ and μ. However, such a generalization is limited to tank geometries with similar baffle arrangements and impeller to tank diameter ratio. Please cite this article as: Cudmore, G.C., et al., A model of impeller whirl for baffled mixing vessels. Journal of Fluids and Structures (2015), http://dx.doi.org/10.1016/j.jfluidstructs.2015.01.010i

G.C. Cudmore et al. / Journal of Fluids and Structures ] (]]]]) ]]]–]]]

18

Table 3 Summary of fluid added coefficients determined from CFD simulations. Re

M nf

C nf

cnf

K nf

kf

10 10a 100 1000 10 000

0.0615 0.0598 0.0216
0.0072
0.0368

0.2980 0.2988 0.1003 0.0907 0.0987

0.1325 0.1419 0.0717 0.0554 0.0308


0.0535
0.0523
0.0490
0.0543
0.0524

0.2775 0.2764 0.0847 0.0733 0.0740

a

n

Case run with double value of eccentricity.

Table 4 The effect of mesh refinement on normalized force magnitude, jF n j. Re

10 1000

Num. of Nodes 270k

870k

1700k

7.77  10
4 3.08  10
4

8.49  10
4 2.86  10
4

8.81  10
4 2.72  10
4

Table 5 The effect of mesh refinement on the fluid added coefficients. Fluid Added Mass Coefficients

Re ¼ 10 270k

Re¼ 1000 870k

Δ%

270k

870k

Δ% 30

M nf

0.062

0.067

8.5


0.0061


0.0087

C nf

0.31

0.30

3.0

0.089

0.091

K nf


0.054


0.056

4.1


0.056


0.056

0.00

cnf


0.13


0.14

4.6


0.054


0.056

4.1

n


0.28


0.29

5.7


0.071


0.074

3.6

kf

2.0

5.7. Mesh independence of the fluid added coefficients In support of developing a CFD methodology with a model resolution sufficient to extract the fluid added coefficients a mesh sensitivity study was undertaken. Since the goal was not to determine detailed flow physics of the mixer, it was expected that a coarser mesh system could be employed. In order to test the mesh dependence of the fluid forces derived from CFD, 3 different meshes were created with 270k, 870k, and 1.7M nodes. Simulations with Re ¼10 and Re¼1000 were qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n2 run on all 3 meshes for a frequency ratio of Ω=2πN ¼ 0:5 and the resulting orbit average fluid force, jF n j ¼ F nT2 þ F N , is compared in Table 4. The effect of mesh refinement on the fluid added coefficients was evaluated by solving cases with Ω=2πN ¼ 0; 7 0:5, and 71.0 and Re¼10 and 1000 using the 270k and 870k node meshes. The results are summarized in Table 5 with the percent change between successive meshes defined as Δ¼

jK nf ;270k
K nf ;870k j ; K nf ;870k

ð49Þ

for example, where K nf ;270k is the fluid added stiffness coefficient computed on the 270k node mesh and K nf ;870k is the fluid added stiffness coefficient computed on the 870k node mesh. From Tables 4 and 5 one can see that the solutions are converging for all Reynolds number as the number of nodes in the computational model increases. An additional level of mesh refinement would help quantify the residual mesh dependence, but this was impractical given the computational requirements for a larger mesh. For example, a single orbit period of data on the 1.7M node mesh requires 4 days of continuous calculations on 16 processors and approximately 10 orbits of data are needed to obtain an average force for a single value of Ω=2πN in the Re ¼1000 case. This study provided relative measures of the errors associated with calculations performed on each mesh and for laminar flow with Re¼10, it can be considered conclusive. However for the turbulent case with Re ¼1000 the LES turbulence model used in this study may also impact the accuracy of the results. Knowledge of the turbulence energy generation and dissipation length scales in the system would have been valuable for comparison purposes (Pope, 2000). Please cite this article as: Cudmore, G.C., et al., A model of impeller whirl for baffled mixing vessels. Journal of Fluids and Structures (2015), http://dx.doi.org/10.1016/j.jfluidstructs.2015.01.010i

G.C. Cudmore et al. / Journal of Fluids and Structures ] (]]]]) ]]]–]]]

19

In the interest of balancing calculation accuracy with computation time it was determined that the 270k node mesh will be used for the complete study. Most of the results given by this mesh differ from the 1.7M node mesh by less than 10%, and the solution can be completed in 1/100th the time required for the 1.7M node mesh. The changes in added mass at the higher Re number are 30% but its value is small and not significant relative to the structural mass. Should more accurate results be required a larger mesh study could be conducted at a later date, however the lowest mesh refinement should provide an acceptable level of accuracy for the present engineering model. 6. Linear model results In this section the fluid added and structural mass, damping and stiffness are combined with the linear model results of impeller whirl to make predictions of the stability and mean square amplitude of vibration that will be compared to the experiments of Kippers and Holloway (2014). The structural mass, damping, and stiffness were Ks ¼1485 N/m, Cs ¼ 5.27 N s/m, M I ¼ 2:22 kg and M sh ¼ 0:277 kg. The dependence of the fluid added mass, damping and stiffness on shaft turning speed were deduced from the Reynolds number dependent coefficients presented in Table 3 as follows: M f ¼ ρD3 M nf ;

ð50Þ

C f ¼ ρD3 2πjNjC nf ;

ð51Þ

cf ¼ ρD3 2πjNjcnf ;

ð52Þ

K f ¼ ρD3 ð2πNÞ2 K nf ;

ð53Þ

n

kf ¼ ρD3 ð2πNÞ2 kf :

ð54Þ

The resulting real components of the eigenvalues, λ=ωnb , given by Eq. (17) as a function of shaft rotation speed, N þ ¼ 2πN=ωnb , are negative, as shown in Fig. 17, for the range
1:28 oN þ o 1:31 indicating stability. Outside of this range one of the eigenvalues is positive, indicating that the response to an input at the corresponding frequency would grow exponentially within the amplitude limits of a linear model. There is a slight asymmetry in the linear model between values for positive and negative N þ owing to the effects of thrust from the pitched blade impeller; an effect that will be neglected in the following. The imaginary portion of the eigenvalues, shown in Fig. 18, indicates that there are two distinct whirls, as indicated by the solid and dashed lines, within the stable range. The higher frequency whirl (both forward and backward) is less damped than the lower frequency whirl. This difference is entirely due to the off-diagonal fluid damping cf which increases in magnitude with N þ and acts to reduce stability of the higher frequency whirl. The modes of the impeller orbit for impeller turning rates of jN þ j ¼ 0.5, 1.0, and 1.3 are shown in Figs. 19 through 21, respectively. Each figure shows two orbits corresponding to the two eigenvalue pairs. All of the orbits are elliptical spirals that tend towards the origin as t increases. The orbit associated with eigenvalue s1 is less damped than that of eigenvalue s2 consistent with Fig. 17. The difference in damping and the ratio of the major and minor axes of the ellipse increases with increasing jN þ j until at the maximum stable turning rate, jN þ j ¼ 1:3, one is an undamped elongated ellipse and the other rapidly decaying toward the origin.

Fig. 17. Plot of the real and imaginary components of the eigenvalues obtained using fluid coefficients of Table 3 as a function of shaft rotation. Vertical lines indicate the limits of Re 4 100.

Please cite this article as: Cudmore, G.C., et al., A model of impeller whirl for baffled mixing vessels. Journal of Fluids and Structures (2015), http://dx.doi.org/10.1016/j.jfluidstructs.2015.01.010i

20

G.C. Cudmore et al. / Journal of Fluids and Structures ] (]]]]) ]]]–]]]

Fig. 18. Plot of the real and imaginary components of the eigenvalues obtained using fluid coefficients of Table 3 as a function of shaft rotation. Vertical lines indicate the limits of Re 4 100. Synchronous whirl indicated by

.

Fig. 19. Impeller whirl mode shapes associated with a turning rate of jN þ j ¼ 0:5. Arbitrary scale and arbitrary orientation. Arrows indicate direction of increasing time.

Fig. 20. Impeller whirl mode shapes associated with a turning rate of jN þ j ¼ 1. Arbitrary scale and orientation. Arrows indicate direction of increasing time.

Please cite this article as: Cudmore, G.C., et al., A model of impeller whirl for baffled mixing vessels. Journal of Fluids and Structures (2015), http://dx.doi.org/10.1016/j.jfluidstructs.2015.01.010i

G.C. Cudmore et al. / Journal of Fluids and Structures ] (]]]]) ]]]–]]]

21

Fig. 21. Mode shapes associated with s1 (left) and s2 (right) for a turning rate of jN þ j ¼ 1:3. Arbitrary scale and orientation. Arrows indicate direction of increasing time.

Fig. 22. Magnitude of the transfer function, H þ ðω þ Þ for various shaft turning rates: Forward whirl corresponds to ω þ 40 and backward whirl to ω þ o 0.

N þ ¼ 0:0, – – N þ ¼ 0:5,    N þ ¼ 1:0, and 
N þ ¼ 1:3.

The magnitude of the transfer function, H þ ðω þ Þ ¼ HðωÞM s ω2nb , as a function of the whirl speed and direction, ω þ , is shown in Fig. 22 for N þ ¼ 0:0, 0.5, 1.0, and 1.3. For N þ ¼ 0 the magnitude of the transfer function is symmetric about ω þ ¼ 0 with maxima at jω þ j ¼ 1 as would be expected for a vibrating shaft in the absence of turning. Increasing the shaft rotation rate, however, produces asymmetry in the magnitude of the transfer function. For example, the magnitude of the maxima occurring at a frequency around ω þ ¼ 0:8 decreases, suggesting that the damping of forward whirl grows with N þ . Whereas for ω þ o 0 the maximum of the transfer function shows an initial drop in amplitude from N þ ¼ 0:0 to N þ ¼ 0:5, and then increases with N þ as the threshold of instability is approached. For N þ o 0 the transfer function would be the same as that plotted in Fig. 22 except with the H þ ðω þ Þ reflected across the ω þ ¼ 0 line. The mean square radial deflection of the impeller for an arbitrary impeller mass eccentricity can be determined using Eq. (23) and is shown in Fig. 23 for clockwise and counterclockwise rotation. The amplitude increases rapidly with shaft rotation rate to a maximum at N þ ¼ 0:75 before gradually declining. No instability is predicted because the rotating imbalance force rotates at the same rate as the impeller and in the same direction and therefore only produces forward whirl (ω þ 40) which, as described in Fig. 22, is diminished at all ω þ 4 0 as the impeller rotation rate increases. The mean square radial deflection for an arbitrary uniform spectral density of fluid forces, Sno , has been calculated using Eq. (23). The results shown in Fig. 23 indicate that as the shaft rotation rate increases in either direction the mean square response increases exponentially up to the points of instability at N þ ¼ 7 1:3. According to the linear model the fluid excitation forces generate whirl instability (whereas the imbalance force cannot) because the fluid excitation has both forward and backward driving components due to the random phase between the force components fY and fZ. Please cite this article as: Cudmore, G.C., et al., A model of impeller whirl for baffled mixing vessels. Journal of Fluids and Structures (2015), http://dx.doi.org/10.1016/j.jfluidstructs.2015.01.010i

G.C. Cudmore et al. / Journal of Fluids and Structures ] (]]]]) ]]]–]]]

22

Fig. 23. Mean square amplitude of impeller vibration resulting from rotating imbalance forces and random fluid excitation forces according to the linear model. Amplitudes are proportional to the mass eccentricity and the power spectral density of fluid forces respectively. Table 6 Fluid added and structural mass, damping and stiffness at the onset of instability. Mf Cf cf Kf kf


0.625 51.6 16.9
876 1234

Ms Cs cs Ks ks

2.35 4.9 0.0 1484
4.9

kg kg/s kg/s kg/s2 kg/s2

7. Discussion Whirl instability of a pitched blade impeller and its avoidance are practical concerns in the design and operation of baffled mixing vessel. The linear model presented in this study predicts that whirl stability will occur where the real part of Eq. (17) is negative (Cudmore et al., 2013) 2

Mk
cCk
C 2 K o0:

ð55Þ þ

At the onset of instability, N ¼
1:3 or N ¼
5:1 Hz, the corresponding values for the fluid added and structural mass, damping and stiffness have the values listed in Table 6. At this point, the total mass, M, is dominated by the structural contribution, the damping terms, C and c, and the off-diagonal stiffness, k, are dominated by the fluid contribution, and the diagonal structural stiffness, Ks, has been significantly offset by a negative fluid stiffness, Kf. Using these approximations Eq. (55) can be written as 2

Mkf
cf C f kf
C 2f ðK s þ K f Þ o0:

ð56Þ

The first term of Eq. (56) is positive, and therefore destabilizing, and is proportional to the off-diagonal stiffness and impeller mass. The second term is stabilizing but smaller than the first and both terms increase in proportion to N4 with a fixed ratio. The third term, proportional to the square of rotational speed and the diagonal damping, is then responsible for maintaining stability of the whirling motion. However, K f p N 2 , and has a negative sign, so the total stiffness, K ¼ K f þ K s , decreases with increasing N until the stability limit is reached. One can therefore conclude that whirl instability in this case results from (1) the presence of off-diagonal stiffness and (2) a combination of negative fluid stiffness and inadequate structural stiffness. A result that has practical implications for shaft selection. It is interesting to note that when Kf ¼0 the third term of (56) is proportional to N2 while the destabilizing terms are proportion to N4 so that instability must arise at high rotational speeds. Conversely, stability at high N, where fluid forces are dominant, can only occur in the presence of positive values of Kf. The stability criterion, Eq. (56), also predicts instability at N þ ¼ 1:3 because each term of Eq. (56) increases with even powers of N. The above description of impeller whirl instability is similar to that reported by Vance (1988) for a rotating disk with cylindrical containment with the primary difference being the presence of the baffles in the present case. The present problem also has some analogous features to fluid elastic instability of tube arrays as described in Paidoussis et al. (2011). This includes off-diagonal stiffness elements of opposite sign coupling the degrees of freedom and negative fluid added stiffness that can overwhelm the structural stiffness to produce static divergence. Kippers and Holloway (2014) reported the mean square radial displacement, E½r þ  where r þ ¼ r=D, for the 254 mm impeller used in the present CFD study provided in Fig. 24. For the range N þ o 0, the amplitude rises exponentially with decreasing N þ . The highest speed that could be tested without excessive contact between blades and baffles was in the Please cite this article as: Cudmore, G.C., et al., A model of impeller whirl for baffled mixing vessels. Journal of Fluids and Structures (2015), http://dx.doi.org/10.1016/j.jfluidstructs.2015.01.010i

G.C. Cudmore et al. / Journal of Fluids and Structures ] (]]]]) ]]]–]]]

Fig. 24. Estimated mean square vibration amplitude for counter-clockwise and clockwise whirl calculated using Sno ¼ 5:4  10
7 . response for forward whirl, – – mean square response for backward whirl,  experimental results due to Kippers and Holloway (2014).

23

mean square

vicinity of N þ ¼
1:2. For N þ 4 0 the amplitude increases up to N þ ¼ 0:75 after which it declined before rising with a reduced slope. The measured mean square amplitude is consistent with the existence of instability and its approximate onset predicted by the model for N þ o 0. However, the symmetry of the instability, with respect to shaft turning direction, N, predicted by the model was not present in the experiment. The reason for the experimentally observed asymmetry is still being studied but one hypothesis is that the unequal support of the impeller by the cantilevered shaft is responsible. This effect was not included in the linear model because the CFD simulation, on which the fluid added mass, damping and stiffness were based, only considered an impeller that followed a prescribed circular orbit in the horizontal plane. The mean square radial deflection due to random fluid excitation forces can be determined from the linear model, using Eq. (26). The results obtained using the low frequency approximation to the power spectral density, Sno ¼ 5  10
7 , are compared in Fig. 24 to the experimental data of Kippers and Holloway (2014). Note that the imbalance forces in this case are negligible due to the high grade of balance used in the experiment (average geometric eccentricity 0.53 mm 70.19 mm). For the range N þ o 0, the predictions are qualitatively correct with both measured and predicted amplitudes increasing with impeller speed towards the point of instability at N þ 
1:3. The comparison is less favorable for N þ 4 0 with the model prediction in agreement with the data for N þ o 0:75 beyond which the experimental mean square displacement falls away from the predictions. Kippers and Holloway (2014) also found that the mean square radial deflection could be described by a rotationally symmetric normal distribution and that the impeller orbits near the threshold of instability (Nn ¼
1:2) contained significant structure including elongated ellipses that bear some similarity to those shown in Fig. 21. In order to further explore the issue of asymmetry in the mean square amplitude with respect to shaft turning direction consider the transfer function shown in Fig. 22. It has a strong asymmetry about ω þ ¼ 0 so that the backward whirl cases (ω þ o0) become less damped with increasing shaft speed while the forward whirl cases (ω þ 4 0) become more damped. These contributions to the radial mean square displacement due to negative whirl and positive whirl can be separated as Z 0 E½r þ 2 ω þ neg ¼ S0 jHðωÞj2 dω; ð57Þ
1

E½r þ 2 ω þ pos ¼ S0

Z

1

jHðωÞj2 dω:

ð58Þ

0

Both Eðr þ 2 Þω þ neg and Eðr þ 2 Þω þ pos are plotted in Fig. 24 for the full range of N þ . Note that when ω þ and N þ have the same sign it is forward whirl and when they have opposite signs it is backward whirl. The backward whirl matches the measured mean square amplitude quantitatively in the interval
1:2 oN þ o 0 while the sum of forward and backward whirl contributions matches in the interval 0 o N þ o0:75. In the interval N þ 40:75 the forward whirl gives the better representation of the data. 8. Summary and conclusions A linear model of impeller whirl which used fluid added mass, damping and stiffness to model fluid forces on the orbiting impeller was introduced as a potential method to predict whirl instability and mean square radial impeller deflection. The fluid forces acting on a rotating pitched blade impeller were derived from CFD simulations in which the Please cite this article as: Cudmore, G.C., et al., A model of impeller whirl for baffled mixing vessels. Journal of Fluids and Structures (2015), http://dx.doi.org/10.1016/j.jfluidstructs.2015.01.010i

24

G.C. Cudmore et al. / Journal of Fluids and Structures ] (]]]]) ]]]–]]]

impeller followed a prescribed circular orbit over a range of orbit frequencies. The fluid added coefficients of mass, damping, and stiffness for whirling motion of the impeller depend on Reynolds number only and were found to have nearly constant values in the Reynolds number range 100 oRe o10 000 for the geometry studied. It was observed that for rates of rotation greater than half the natural bending frequency, jN þ j ¼ j2πN=ωnb j4 0:5, the fluid added mass is relatively small, the damping and off-diagonal stiffness are dominated by the fluid forces, and the fluid added diagonal stiffness is the same order as the structural stiffness and of opposite sign. The linear model predicts stability of the impeller motion within the range
1:3 o N þ o 1:3 which is consistent with experimental data. The model predicts whirl instability at the boundaries of this interval which is observed experimentally for N þ ¼
1:2. An examination of the model eigenvalues suggests that the instability is the result of negative fluid added stiffness and that the safe operating speed of the impeller can be increased by increasing the structural stiffness of the shaft. The predicted mean square radial impeller deflection was calculated using the linear model and the power spectral density of the fluid forces as determined from LES simulation. The predictions were found to be within 15% of the experimental results for shaft speed range
1:2 o N þ o0:75. However, for N þ 4 0:75 the predicted mean square response was qualitatively different from the measured values which showed a decrease in radial deflection in the vicinity of N þ ¼ 0:75. The method for determining the fluid added coefficients that has been presented can be extended to different mixer configurations. An advantage of this method over a fully coupled fluid–structure interaction solution of the problem is that the fluid added coefficients, once determined, can be used to make predictions for impeller vibration for a variety of operating conditions of the same mixer configuration. Acknowledgments The authors acknowledge financial support of NSERC and NOVA Chemicals through a Collaborative Research and Development Grant (No. CRDPJ 320426). References ANSYS CFX, 2009a. Discretization of the Governing Equations, Discretization and Solution Theory. ANSYS CFX-Solver Theory Guide, Canonsburg, PA. ANSYS CFX, 2009b. Governing Equations, Basic Solver Capability Theory. ANSYS CFX-Solver Theory Guide, Canonsburg, PA. Adams, M.L., 2001. Rotating Machinery Vibration. Marcel Dekker Inc., NY. Berger, T., Fischer, M., Strohmeier, K., 2003. Fluid-structure interactions of stirrers in mixing vessels. Journal of Pressure Vessel Technology 125, 440–445. Brennen, C.E., Acosta, A.J., 2006. Fluid-induced rotordynamic forces and instabilities. Structural Control and Health Monitoring 13, 10–26. Brucato, A., Ciofalo, M., Grisafi, F., Micale, G., 1998. Numerical prediction of flow fields in baffled stirred vessels: a comparison of alternative modeling approaches. Chemical Engineering Science 53 (21), 3653–3684. Cudmore, G.C., Holloway, A.G.L., Gerber, A.G., 2013. Whirl instability of a rotating impeller in a baffled mixing vessel. In: Proceedings of PVP2013. Fuhrmann, D.R., 1999. Complex random variables and stochastic processes. In: Vijay, E., Madisetti, K., Williams, D.B. (Eds.), Digital Signal Processing Handbook, CRC Press, Boca Raton. Gotz, S., Sperling, R., Liepe, F., Jembere, S., 1997. Numerical determination of the three-dimensional velocity distribution in a baffled pitched blade impeller stirred vessel. Chemical Engineering & Technology 20, 596–605. Guinzberg, A., Brennen, C.E., Acosta, A.J., Caughey, T.K., 1994. Experimental results for the rotordynamic characteristics of leakage flows in centrifugal pumps. Journal of Fluids Engineering 116, 110–115. Hirt, C.W., Amsden, A.A., Cook, J.L., 1974. An arbitrary Lagrangian–Eulerian computing method for all flow speeds. Journal of Computational Physics 14, 227–253. Jery, B., Brennen, C.E., Caughey, T.K., Acosta, A., 1985. Forces on centrifugal pump impellers. In: Proceedings of the Second International Pump Symposium. Kippers, N., Holloway, A.G.L., 2014. Experiments on the whirling of pitched blade impellers in baffled mixing vessels. Journal of Fluids and Structures 49, 29–52. Newland, D.E., 1975. An Introduction to Random Vibrations and Spectral Analysis. Longman Group Limited, London. Nicoud, F., Ducros, F., 1999. Subgrid-scale stress modelling based on the square of the velocity gradient tensor. Flow, Turbulence and Combustion 62, 183–200. Oldshue, J., 1983. Fluid Mixing Technology. McGraw-Hill Publications Co, New York. Paidoussis, M.P., Price, S.J., deLangre, E., 2011. Fluid-Structure Interactions: Cross-Flow-Induced Instabilities. Cambridge University Press, Cambridge, UK. Pope, S.B., 2000. Turbulent Flows. Cambridge University Press, Cambridge, UK. Rice, M., Hall, J., Papadakis, G., Yianneskis, M., 2006. Investigation of laminar flow in a stirred vessel at low Reynolds numbers. Chemical Engineering Science 61, 2762–2770. Vance, J.M., 1988. Rotordynamics of Turbomachines. John Wiley & Sons, Inc., NY. Zadghaffari, R., Moghaddas, J.S., Revstedt, J., 2009. A mixing study in a double-Rushton stirred tank. Computers & Chemical Engineering 33, 1240–1246.

Please cite this article as: Cudmore, G.C., et al., A model of impeller whirl for baffled mixing vessels. Journal of Fluids and Structures (2015), http://dx.doi.org/10.1016/j.jfluidstructs.2015.01.010i