Estimation of exhaust gas aerodynamic force on the variable geometry turbocharger actuator: 1D flow model approach

Energy Conversion and Management 84 (2014) 436–447

Contents lists available at ScienceDirect

Energy Conversion and Management journal homepage: www.elsevier.com/locate/enconman

Estimation of exhaust gas aerodynamic force on the variable geometry turbocharger actuator: 1D flow model approach Fayez Shakil Ahmed a, Salah Laghrouche a,⇑, Adeel Mehmood b, Mohammed El Bagdouri a a b

Laboratoire IRTES-SET, Université de Technologie de Belfort-Montbèliard (UTBM), Belfort, France Laboratoire MIPS, Université de Haute Alsace (UHA), Mulhouse, France

a r t i c l e

i n f o

Article history: Received 11 November 2013 Accepted 29 March 2014

Keywords: Diesel engine Exhaust aerodynamic force Variable turbine geometry 1D flow modeling

a b s t r a c t This paper provides a reliable tool for simulating the effects of exhaust gas flow through the variable turbine geometry section of a variable geometry turbocharger (VGT), on flow control mechanism. The main objective is to estimate the resistive aerodynamic force exerted by the flow upon the variable geometry vanes and the controlling actuator, in order to improve the control of vane angles. To achieve this, a 1D model of the exhaust flow is developed using Navier–Stokes equations. As the flow characteristics depend upon the volute geometry, impeller blade force and the existing viscous friction, the related source terms (losses) are also included in the model. In order to guarantee stability, an implicit numerical solver has been developed for the resolution of the Navier–Stokes problem. The resulting simulation tool has been validated through comparison with experimentally obtained values of turbine inlet pressure and the aerodynamic force as measured at the actuator shaft. The simulator shows good compliance with experimental results. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction Fuel economy and reduction in exhaust emission levels are the major challenges of today’s automobile industry [1,2]. These challenges translate directly onto technological improvements in vehicles. Coupled between the intake and exhaust, turbochargers have a direct influence on the engine’s performance in terms of both, the fuel efficiency and emission quality [3]. Yet the performance of the turbochargers themselves is usually calculated through empiric methods such as interpolation of characteristic maps. These maps differ from each other for different class of turbochargers and they must be adapted, requiring expensive tests on specialized rigs [4–6]. The information in these maps is limited to static flow characteristics at the outlets, making it impossible to determine the effects of flow inside the turbocharger. This is specially critical for the turbine section, of Variable Geometry Turbochargers (VGTs), as static characteristic maps provide no information about the effect of flow on the variable geometry vanes. In reality, the geometry of a VGT turbine inlet is determined by the angle of vanes at the inlet, which are controlled by an actuator. As these vanes regulate the quantity of exhaust flow and its angle of incidence on the impeller, they are subjected to a ⇑ Corresponding author. Tel.: +33 612832335. E-mail addresses: [email protected] (F.S. Ahmed), salah. [email protected] (S. Laghrouche), [email protected] (A. Mehmood). http://dx.doi.org/10.1016/j.enconman.2014.03.080 0196-8904/Ó 2014 Elsevier Ltd. All rights reserved.

significant resistive force by the exhaust gas. Estimation of this aerodynamic force on the vanes and its effect on the VGT actuator requires exhaust flow modeling inside the turbine as well. While the task appears daunting, it has now become absolutely necessary for VGT control at levels of precision that are compatible with the energy efficiency requirement of today. In the literature review, such models can be found for multistage compressors and turbines used in the heavy industry [7,8], based on wave energy conversion [9]. Analytical and lumped volume (0D) models are simple but not accurate [10]. In the automotive industry, such models are discussed only to estimate the exhaust pressure by considering engine dynamics, i.e. engine speed and load or turbocharger speed [11–13]. Other methods involve the lumped volume approach where the turbocharger is considered as an object with single constant volume [14,15]. The drawback of this approach is that we require turbocharger maps for further analysis and the quantities like air density and flow rate are not available. On the contrary 2D or 3D models are too complex to be solved by numerical analysis in real-time for control objectives [16–18]. Therefore, 1D models of the air flow through a turbocharger have received more attention of researchers. The existing 1D models deal mostly with fixed geometry compressors and turbochargers [19]. The modeling of aerodynamic forces is also limited to impellers [20,21]. As the role of aerodynamic forces on variable geometry vanes cannot be neglected in VGTs [22], it is important to obtain accurate yet computationally simple

F.S. Ahmed et al. / Energy Conversion and Management 84 (2014) 436–447

437

Nomenclature

a te b1 bmetal s db di

q fs f fp ft ! ! r ! !

sf

a An As Ao CL CD Cp e ep e0 Ec ! F a

Fluid inlet and outlet angle Vane width at exit Angle of attack Camber line angle Maximum inter-vane distance Stagger angle Incidence Fluid density Pressure loses due to secondary effects Pressure loses Pressure loses due to profile Pressure loses due to play Tensor of constraints Tensor of viscous constraints Velocity of sound Channel cross section Channel cross section at start Channel cross section at exit Lift coefficient Drag coefficient Specific heat at constant pressure Energy Maximum thicknesses of the blade Total energy Kinetic energy External force on the fluid Tip to maximum camber distance

exhaust-flow model of the VGT turbine geometry section, in order to calculate and compensate these forces in future control applications [23,24]. This paper presents the development of a simulation tool aimed at estimation and prediction of the aerodynamic force on the vanes of the variable turbine geometry of VGTs, using 1D exhaust gas flow model. This tool has been designed primarily for automobile manufacturers, enabling them to preempt the aerodynamic force at basic design level, in order to improve the overall efficiency of the final engine system as well as emissions. It also aids in scientific endeavors of developing a better understanding of the dynamics of the modern engine and its air-path. The major development in this paper is the 1D exhaust gas flow model through the variable turbine geometry section, designed specifically for turbochargers. The flow model is derived from Navier Stokes equations for 1D flow, including losses due to geometry, friction and impeller blade, commonly known as source terms. The differential equations of this model are resolved with an implicit numerical solve, therefore its stability is guaranteed. The exhaust flow model provides the pressure and flow rate throughout the variable geometry section, which permits to calculate aerodynamic force on the turbine geometry vanes and on the actuator. The impact of exhaust gases on the VGT vanes is transformed to actuator’s stem via calculation of its components through the mechanical linkage between the vanes and the actuator. This simulation tool can easily be integrated in conventional automotive simulation platforms like AMESim and control design softwares like Matlab. Its performance is validated through comparison with experimentally obtained values of the turbine inlet pressure and aerodynamic force measured on the actuator shaft. The rest of the paper is organized as follows: In Section 2, the mechanism of the VGT is introduced. The one dimensional CFD model of air through the channel between the exhaust manifold

h 0 h k b2 m t M o p r p0 Q Q ext R Re Re ref S c Ds t0 Dt ! u ! V Vm Vu ! W v b

Enthalpy Total enthalpy Blade clearance Angle of attack at exit Arc length Maximum width Mach number Opening inter-vane distance Pressure Radius Total pressure Source terms External forces Perfect gas constant Reynolds number Reference Reynolds number Cross sectional area. Chord length Variation in entropy Total temperature Time interval ! Normal to m Velocity vector Meridional velocity Tangential velocity Velocity vector in reference frame Volume Maximum camber

and the turbine is presented in Section 2 and the resulting force on the actuator is calculated in Section 3. In Section 4, simulation results of 1D CFD model are compared with experimental results for validation. These results were obtained at different speeds and loads with and without pressure losses. Section 5 concludes the results. 2. Aerodynamic force modeling The power or torque of a diesel (compression-ignition) engine can be improved by forcing more air into the cylinder, achieving high Air to Fuel ratio (AFR) [25]. The air mass-flow is increased by boosting the pressure in the intake manifold [26] through turbochargers. Turbochargers consist of a turbine and a centrifugal compressor, both connected by a rigid shaft. The turbine runs on the exhaust gas flow, utilizing energy that would otherwise have been wasted. The compressor, running through the turbine forces air into the intake manifold, providing the pressure boost. 2.1. VGT mechanism In modern turbochargers, the compressor and turbine speed is controlled by varying the geometry of the turbine inlet, thereby changing the exhaust flow characteristics. A variable geometry turbocharger (see Fig. 1) has moveable vanes located around the turbine, which change the turbine inlet geometry in accordance with the engine operating point. The position of the VGT vanes is controlled by a dedicated electrical or electro-pneumatic linear VGT actuator that is connected to the vanes through a rack and pinion crankshaft and a unison ring. The latter ensures that all the vanes move together. The aerodynamic force exerted by the exhaust gas flow on the turbine geometry vanes acts as a resistive force to the actuator

438

F.S. Ahmed et al. / Energy Conversion and Management 84 (2014) 436–447

fluid in the vertical direction is zero V l ¼ 0. Henceforth, we get the following equation

ZZZ

Fig. 1. VGT with an electro-pneumatic actuator system.

and inhibits proper control. The internal flow characteristics of the turbine geometry need to be known in order to estimate this force. In practice we only know the gas properties before it enters the turbine. Therefore, in order to determine the pressure, density and temperature inside the turbine, the gas flow through the turbine inlet vanes needs to be modeled. Once these variables are found, the torque acting on the turbine vanes can be calculated and the net force acting on the actuator can be derived from force transfer analysis through the mechanical linkages. The torque on the vanes is transferred to the crank and actuator via two levers. The primary lever is connected to the actuator crank and the secondary lever is connected to all the vanes through a unison ring. In Fig. 2 an isometric view of the turbine and linkage is shown. In this work our focus will be on the prediction of exhaust air pressure and the aero forces acting on the pneumatic actuator, caused by exhaust air.

2.2. Hypothesis and problem formulation with Navier–Stokes equations The Navier–Stokes equations are used to model the aerodynamic force exerted by the exhaust gases coming from the exhaust manifold and passing through the turbocharger. These equations represent the conservation laws for mass, energy and momentum of the fluid in all directions [27]. The Navier–Stokes equations can be reformed to obtain 1D model of the fluid with some assumptions on the fluid’s direction of flow. The velocity vector ! for fluid is V ¼ ½V m V u V l T with condition that the velocity of the

Fig. 2. Variable turbine geometry: vane actuation mechanism.

dv ¼

Z

S dm:

ð1Þ

The term S is the cross section of the channel through which the gas passes, and dv ¼ dm du dl is an elementary volume. m; u and l represent respectively, the meridional, tangential and vertical axes of the fluid flow, as explained in Fig. 3. For the perfect gas the tensor for the viscous forces can be ! ! ! ! consider as equal to zero i.e. sf ¼ 0 . Additionally, for initial development of Navier–Stokes equations, it is consider that external source terms are zero. They will be introduced in the model later on. 2.3. Model derivation 1D flow models are equivalent to the generalization of fluid flow in a thin tube or a narrow channel. This interpretation permits us to analyze the properties of the fluid; velocity, density and energy change only along the tube or the channel. The law of conservation of mass for a fluid can be derived from the Navier–Stokes equation. The equation is also known as the continuity equation and its differential form is as follows [28]

@Sq @SqV m þ ¼ 0: @t @m

ð2Þ

The terms q and V m are the density and meridional velocity of the fluid, respectively. The momentum equation is based on the principle of conservation of momentum. In the derivation, we assume that V l ¼ 0, i.e. the fluid has only tangential V u and meridional V m velocity components. Therefore, the momentum equation obtained by the equivalent meridian velocity (velocity of fluid in the direction of fluid flow) is given as

  2 @SqV m @Sq p=q þ V m V 2 @r @S þ ¼ Sq u þp : @m @m @t r @m

ð3Þ

The momentum equation for the tangential component of the velocity is given as

@SqV u @SqV m V u V m V u @r þ ¼ Sq : @t @m r @m

ð4Þ

Fig. 4 shows the direction of fluid through a variable geometry turbocharger. Cartesian components of the velocity are also shown in this figure. The direction of air flow at the turbine outlet is perpendicular to the flow at inlet. In this figure, m is the meridian coordinate and r is the radius of the fluid channel. The conservation

Fig. 3. Velocity components of air-flow through a vane.

F.S. Ahmed et al. / Energy Conversion and Management 84 (2014) 436–447

439

2.4. Source terms The source terms in the Euler equations are used to represent the influence of the blade rows and the flow-path on the fluid. The source terms vector Q is the vector of all contributions which were initially neglected in the derivation. It may be further divided into four distinct categories by their physical meanings

Q ¼ Qg þ Q b þ Q f þ Q c;

ð10Þ

where Q c is the cooling source term and it is neglected in this derivation i.e. Q c ¼ 0. The other terms, for example Q g ; Q b and Q f are the source terms due to the geometry of the volute, impeller blade force and viscus friction, respectively. These terms are described below. Fig. 4. Exhaust gas flow through VGT.

of energy, like conservation of mass and momentum, represents that total energy in a fluid remain constant. In differential form, it is expressed as 0

@Se0 q @SqV m h V m V u @r þ ¼ SqX ; @m @t @m

ð5Þ

0

where e0 and h are the total energy and enthalpy of the fluid. Eqs. (2)–(5) can be combined together and can be written in the nondimensional vector form

@U @E þ ¼ Q; @t @m

ð6Þ

where U; E and Q are conservative variable vector, flux vector and the source term vector respectively. In the vector from these terms are given as

3 2 2 3 qV m S 3 0 qS   7 6 6 V u W u @r @S 7 6 qV S 7 7 6 p þ qV 2m S 7 6 Sq r @m þ p @m m 7 6 7; Q ¼ 6 7: U¼6 7; E ¼ 6 7 7 6 6 V V @r m u 4 qV u S 5 S q 5 4 qV m V u S 5 4 @m r @r qe0 S SqV m V u X @m q V m h0 S 2

ð7Þ T

Let U ¼ ½U 0 U 1 U 2 U 3  , then the conservative flux vector E can be written as a function of components of the vector U. By rewriting the flux vector, we get

E0 ¼ U 1 ; 1 U 21 U 22 þ E1 ¼ ðc  1Þ U 3  2 U0 U0 E2 ¼ E2 ¼

þ

U 21 ; U0 ð8Þ

U1 U2 ; U0

cU 3 

!!

c  1 U 21 2

U2 þ 2 U0 U0

!!

2.4.1. Geometry of the flow-path The source term Q g , represents the variation in fluid velocities and pressure drop across the exhaust air flow-path. The air flow direction through the vein between the exhaust manifold and the turbocharger is shown in Fig. 4, where r represents the mean-line radius. This term can be obtained by considering the right hand side of Eqs. (2)–(5). The final form of geometrical source term is given by

2

3

0

6 Sq V u V u @r þ p @S 7 6 @m @m 7 r Qg ¼ 6 7: 4 Sq V m V u @r 5 r

0 2.4.2. Evaluation of the output angle The turbine impeller and the VGT vanes, both modify the fluid angle by directing it from their leading edge to the trailing edge, as seen in Fig. 6. Fluids angles a and b are defined as follows: The rotating blades of the impeller change the angle b, while the VGT vanes change the angle a. The following model has been proposed in [7] for evaluating the output angle of the fluid:  For M 2 < 0:5: the output angle is determined by the nomogram 7.  For 0:5 < M 2 < 1, a linear transition   is carried out.  For M2 ¼ 1 we get: a2 ¼ cos1 AAos . 2.4.3. Blade source terms This term is implemented to represent the isentropic influence of blade rows. The force applied by the blades onto the fluid can be deduced from the knowledge of the velocity triangles, see Fig. 5. An angular momentum balance over the blade row provides the tangential component of the blade force. Here we have used the tangential velocities along the blade deflection induced by the e u i (see Fig. 7). given component for V

e u i ¼ Vm i ; V tanðai Þ

U1 : U0

ð12Þ

In order to find a numerical solution of Eq. (6), the Jacobian flux vector A needs to be calculated. For this, the flux vector E needs to be written as a function of U. The Jacobian vector is obtained by taking @E the partial derivatives of E with respect to U [18], i.e A ¼ @U , we get 3 2 0 1 0 0 ðc3ÞV 2m þðc1ÞV 2u 6 ðc  3ÞV m ðc  1ÞV u c  1 7 7 6 2 7: A¼6 7 6 V V 0 V mV u u m 5 4    V 2u 2 2 c2 2c3 2 a2 a2 V þ V  V  ð c  1ÞV V c V  V m m u m m u m 2 2 c1 c1 2 ð9Þ

ð11Þ

@m

Fig. 5. Velocity triangle.

440

F.S. Ahmed et al. / Energy Conversion and Management 84 (2014) 436–447

Fig. 9. Forces on the vane.

where X is the angular velocity of the blades. The blade force in the tangential direction is given by

Fig. 6. Blade geometry.

e u i  r i1 V u i1 Þ F bu i ¼ ðr i1 V

qm : Si r i ðmiþ1  mi1 Þ

ð14Þ

The forces acting on the blade are perpendicular to the referential ! velocity W . Henceforth, it can be assumed that the axial and tangential forces do not generate losses

F bm i W m i þ F bu i W u i ¼ 0;

ð15Þ

where W m is the relative velocity vector. By analyzing the velocity triangles for the fluid through the single stage turbocharger, it can be deduced that in the meridional direction relative meridional velocity is equal to the absolute velocity, thence W m i ¼ V m i . Therefore, the blade force in the meridional direction is given by

F bm i ¼ F bu i

Wu i : Vm i

ð16Þ

The source term related to the action of impeller blade on the turbocharger vanes is given below

2

3

0

7 7 7: 5

6 F S 6 bm i i Qbi ¼ 6 4 F bu i Si

ð17Þ

r i XF bu i Si

Fig. 7. Gas outlet angles vs. cos1 ðo=sÞ (straight back blades operating at low mach numbers).

2.4.4. Pressure losses As described previously, the source term obtained from friction depends upon the pressure losses in the turbocharger. These losses are mainly due to the turbine geometry, play in its components and some secondary effects. In [7], two different methods are presented to determine these losses. These methods are described below. 2.4.5. Interpolation based method The pressure losses in the vanes increase the losses due to tip clearance in a blade row of a turbine increases, which results in a decrease of turbine efficiency. The gas mass flow and the power output corresponding to a fixed turbine speed and pressure ratio also change. Typical profile losses on conventional blades intermediate between nozzle and impulse blades may be interpolated in the following manner

fp ¼

Fig. 8. Finite element division of VGT vane.

where ai is the angle of the fluid with respect to blade comber line. The index i comes from the discretization of the blade in finite number of elements, as presented in Fig. 8. By writing the Euler equation and the equation of variation in energy, we get (see Fig. 9)

e u i  r i1 XV u i1 Þ r i F bu i X ¼ ðr i X V

qm ; Si ðmi  mi1 Þ

Y pðb1 ¼0Þ þ



b1

a2

2



!   ep  Y pðb1 ¼a2Þ  Y pðb1 ¼0Þ ; 0:2c

ð18Þ

where ep is the pitch to chord ratio. The term Y p represents the profile losses at different inlet angles for blades. These losses are determined through the mapping between inlet and outlet angles for both the blade and fluid: and the profile. Stodola quoted a simplified method to determine the losses associated with the tip clearance. The formula to determine clearance losses is given by 1:4

ð13Þ

ft ¼ 6:26

k ; hc

ð19Þ

441

F.S. Ahmed et al. / Energy Conversion and Management 84 (2014) 436–447

where hc is the blade height and k is the radial clearance. In this approach the secondary losses are neglected, therefore we considered that

fs ¼ 0:

ð20Þ

2.4.6. Angle ratio based method In this method, the profile losses are determined with more sophisticated and more precise model. It depends upon the ratio of angles a2 and b1 presented in Section (2.4.2), as explained in [19,29].

fp ðdi¼0Þ ¼

!  2 b   ep a12 b1  Y pðb1 ¼0Þ þ Y pðb1 ¼a2Þ  Y pðb1 ¼0Þ 0:2c a2

ð21Þ

with Y pðb1 ¼0Þ and Y pðb1 ¼a2 Þ are obtained from the maps between pitch to chord ratio and profile loss coefficient. The pressure losses due to blade tip clearance ft is given by the following equation

 2 k C L cosða2 Þ2 ft ¼ 0:5 ; hc cs cosðam Þ3

To obtain the source term that contains the momentum and energy equations can be determined from the axial and tangential components of the friction force.

F fm ¼ qt

R lnð1  fÞ Vm qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; Dm 2 V m þ ðV u  r XÞ2

ð29Þ

where Dm is the small step size in the meridional direction. The tangential component of the force is given by

R lnð1  fÞ V u  rX qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : Dm 2 V m þ ðV u  r XÞ2

F fu ¼ qt

ð30Þ

The source term related to friction forces is given by the following expression

2

0

6 F S 6 fm i i Qf i ¼ 6 4 F fu i Si

3 7 7 7: 5

ð31Þ

r i X F fu i Si

ð22Þ

where C L is the lift coefficient. The term am is the mean gas angle. Its expression for this angle is given in [28]

2.5. Numeric resolution of 1D Navier Stokes equations

 2 C L cosða2 Þ2 fs ¼ k s cosðam Þ3 c

The objective of this solver is to calculate the point force exerted by the gas on the vanes, and then to find out the total force experienced by the actuator. Let us recall Eq. (6), which represents fluid flow in matrix form

ð23Þ

0 2 1 A2

B A C with k ¼ f @ 1þ1 I:D A is the ratio of drag and lift coefficient. It can also O:D

be mapped by the ratio of inlet and outlet flow area

which can be written as

Here I:D ¼ r  h2c and O:D ¼ r þ h2c . The inlet and outlet flow area depend upon the blade inlet and the fluid outlet angle

A1 ¼ An 0 cosðb1 Þ; A2 ¼ An 2 cosða2 Þ:

ð24Þ

The ratio CsL c can also be represented as a function of fluid inlet, outlet and mean angle

CL s c

¼ 2 ðtanða1 Þ  tanða2 ÞÞcosðam Þ:

ð25Þ

2.4.7. Friction source terms In the Euler equations, the entropy production is due to a separate distributed friction force F f in the momentum equations. This force causes the loss in fluid pressure when it moves through the channel. The pressure drop due to this friction may cause the variation in momentum and energy loss in the fluid. This friction force is aligned with the relative velocity by the following expression. The negative sign indicate that it has the opposite direction

! ! ! W F f ¼ jF f j ! : jW j

ð26Þ

The change in entropy of the fluid can be represented by the following formula

Ds ¼ R lnð1  fÞ;

@U @U þA ¼ Q; @t @m

ð32Þ

@E where A ¼ @U is the Jacobian of E. This equation will be discretized in order to solve it numerically. As explained in [28,19], in order to guarantee the convergence of the solver, the positive and negative components of the vector flow A need to be separated. This is achieved by separating the jacobian of E in two matrices, one composed of positive eigenvalues Aþ and the other composed of negative eigenvalues A . Starting first by the diagonalization of A, we suppose that there exists a matrix X and a diagonal matrix D such that AX ¼ XD, with

2

3

Vm

6 6 D¼6 4

7 7 7; 5

Vm a þ Vm

where a is the velocity of sound. Let Dþ and D be the matrices composed respectively of positive and negative eigenvalues of A. In our case, the flow is subsonic, i.e. V m < a. Therefore we get

2

Vm

6 6 Dþ ¼ 6 4

Vm a þ Vm 0

3

2

7 7 7; 5

6 6 D ¼ 6 4

3

0

7 7 7: 5

0 0 a þ V m

ð34Þ Let Aþ ; A ; Eþ ; E given as

Aþ ¼ XDþ X 1 ; A ¼ XD X 1 ; Eþ ¼ Aþ U;

ð28Þ

ð33Þ

a þ V m

ð27Þ

where R is the gas constant and f represents the pressure losses due to profile, play in the turbine geometry and secondary effects. The force related to these losses can be evaluated from the Gibbs relation [30], which is given below

! Ds jF f j ¼ qt : Dm

@U @E þ ¼ Q; @t @m

E ¼ A U:

ð35Þ

442

F.S. Ahmed et al. / Energy Conversion and Management 84 (2014) 436–447

It is clear that in the supersonic case, i.e. V m > a, we have Dþ ¼ D et D ¼ 0. This implies that

Aþ ¼ A;

A ¼ 0;

Eþ ¼ E;

E ¼ 0:

The index ðnÞ is no longer explicit and DU will be expressed in terms of elements calculated in previous iteration. Rewriting this equation in form of matrix sums,

       @  þ @  þ A þ A   B DU ¼ Dt E þ E  Q : I þ Dt @m @m

2.6. Implicit numerical resolution Once the matrices have been calculated, we need to determine an appropriate numeric solver for the partial differential equation. This requires discretization of time as well as space. In each iteration, the vector U will be calculated at all discrete nodes of each vane. The nodes together form a finite element mesh representing a vane. The operator @U is estimated by @t

ð41Þ

Using first order backward difference for Aþ and Eþ , and first order forward difference for A et E , we obtain for i ¼ 1;    ; Imax

  þ  þ þ A  Aþi1 Aiþ1  Ai I þ Dt i þ  Bi DU mi  mi1 miþ1  mi  þ  Eþ  Eþi E  Eþi1 þ iþ1  Qi : ¼ Dt i mi  mi1 miþ1  mi

ð42Þ

@U DU U ðnþ1Þ  U ðnÞ  ¼ ; @t Dt Dt

We have considered: Dmþ ¼ mi  mi1 and Dm ¼ miþ1  mi . This implies that

where Dt is the discretization time, and U ðnÞ is the value of the vector U calculated at the time nDt. As mentioned before, E is expressed in the form E ¼ Eþ þ E , in order to ensure solver convergence. This implies that



DU @ ¼ ðEþ þ E Þ þ Q ; @m Dt

ð36Þ

where the approximation of @E backward difference and the @m  approximation of @E requires forward difference, explained in @m Table 1. These methods can be formulated by explicit and implicit formulations [31]. In explicit methods, the vector U ðnþ1Þ is estimated from precedent values only, i.e. U n ; Q n ; En , as shown in the following equation

 ðnÞ U ðnþ1Þ  U ðnÞ @E ¼ þ Q ðnÞ : @m Dt

ð37Þ

On the other hand, implicit methods estimate the vectors U; Q ; E together, for the instance n þ 1, using the value of U calculated in the preceding iteration. These methods, as shown in the following equation, are more complicated, yet more precise and stable than explicit methods. Therefore, in our work we have developed the solver using implicit method.

ð38Þ

Using Taylor expansion, Eq. (38) is simplified as

 Q

ðnÞ

@Q ðnÞ þ DU @U

! ¼ 0;

ð39Þ

where DU ¼ U ðnþ1Þ  U ðnÞ . Considering the Jacobian of Q ; B ¼ @Q , we obtain the following @U implicit form: ðnÞ

I þ Dt

@A  BðnÞ @m

!!

!

ðnÞ

DU ¼ Dt

ð43Þ

Dt þ A ; Dmþ i1  þ  Ai Ai ;   B AAi ¼ I þ Dt i Dmþ Dm Dt  AP i ¼ A ; Dm iþ1 þ Eþ  Eþi E  Eþi1 RHSi ¼ i þ iþ1  Q i: mi  mi1 miþ1  mi

AM i ¼ 

@E  Q ðnÞ : @m

ð40Þ

AMi DU i1 þ AAi DU i þ AP i DU iþ1 ¼ RHSi : From here, the problem is reduced to calculating Imax linear equations with Imax unknowns, i.e.

AA0

6 6 AM 1 6 6 6 6 6 6 6 6 6 4

3

AP 0 AA1 .. .

..

.

..

.

..

.

.

..

.

..

.. 2 6 6 6 6 ¼6 6 6 4

Notation

First order forward difference

FOFDðX i Þ

First order backward difference Second order forward difference

FOBDðX i Þ

Second order backward difference

..

.

. AAImax 1

7 7 7 7 7 7 7  DU 7 7 7 7 AP Imax 1 5

AMImax AAImax RSH0  AM0 RSH1 .. . .. .

ð45Þ

3 7 7 7 7 7: 7 7 5

RSHImax 1 RSHImax  AP Imax

Table 1 Numerical differentiation algorithms. Method

ð44Þ

The Eq. (42) can be rewritten in the following compact form:

2

 ðnþ1Þ U ðnþ1Þ  U ðnÞ @E ¼ þ Q ðnþ1Þ : @m Dt !

   þ  Dt þ Ai Ai A D U þ I þ D t   B DU i1 i1 i Dmþ i1 Dmþ Dm    þ  Eþiþ1  Eþi Dt  Ei  Eþi1 þ A D U ¼ D t þ  Q iþ1 i : Dm iþ1 mi  mi1 miþ1  mi

Let us define

þ

DU @ @EðnÞ EðnÞ þ þ DU Dt @m @U



This system can be solved easily for DU, and then for U ðnþ1Þ , we get:

U ðnþ1Þ ¼ U ðnÞ þ DU:

@X i @m @X i @m

X i ¼ mX iþ1 iþ1 mi

SOFDðX i Þ

@X i @m

aþ1Þ X iþ1 aðaþ1ÞX iþ2 ¼ aðaþ2ÞXiaþð ; ðaþ1Þðmiþ1 mi Þ iþ2 miþ1 a ¼ mmiþ1m i

SOBDðX i Þ

@X i @m

a1Þ X iþ1 þaðaþ1ÞX i2 ¼ aðaþ2ÞXi ð ; aðaþ1Þðmiþ1 mi Þ i1 mi2 a ¼ mmim i1

i1 ¼ mX ii X mi1 2

2

It can be seen that AM 0 and AP n of Eq. (45) are the boundary conditions, hence known. According to the stability theory presented in [19], for subsonic flow, static density needs to be imposed on the upstream, static pressure on the downstream and static fluid speed in both sides. The other variables can be interpolated:

F.S. Ahmed et al. / Energy Conversion and Management 84 (2014) 436–447

8 q0 imposed > > > V u 0 imposed > > : p0 ¼ 2p1  p2

transfer as a function of the actuator position, let us consider a unit net torque (1 N m) to obtain a generalized expression. 3.1. Relation between crank handle angle and actuator position Let A be the position of the rotation axis of the crank handle (in 2D, with the direction Y oriented along the crank handle axis towards the center of the turbine). We measure:

8 qi ¼ 2qi1  qi2 > > > V u i ¼ V u i1 V u i2 > > : pi imposed

xðAÞ ¼ 1:433 and yðAÞ ¼ 27:308:

The calculation ends when the residual of @U becomes inferior to a @t certain threshold. As their exist many definitions for this residual, we consider the maximum and average residuals:

!

ðnþ1Þ ðnÞ

U

i;j  U i;j

resmax ¼ max

; ðnÞ



U i;j

!

ðnþ1Þ ðnÞ

Im ax X

U 3 1 X

i;j  U i;j

resav g ¼

:

ðnÞ

4Imax i¼0 j¼0

U

ð46Þ

Once the calculation is terminated; the resultant vector U allows us to calculate the physical quantities, like V u ; V m ; q, on the blades. The meridian and tangential forces on the finite elements of the vanes are then calculated as follows

dF u ðmÞ ¼ qðmÞV m ðmÞhsðmÞ

Also, let C be the point of connection between the crank handle and actuator shaft. We take

~ L1 ¼ kACk;

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x0 ðCÞ2 þ y0 ðCÞ2 ¼ L1 ¼ 15:895; qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi L ¼ x00 ðCÞ2 þ y00 ðCÞ2 ; tanðdÞ ¼

x0 ðCÞ ; y0 ðCÞ

with

ð47Þ

x0 ðCÞ ¼ xðCÞ  xðAÞ; 00

x ðCÞ ¼ xðCÞ  xðBÞ;

where h et s are the height and distance between two vanes. The total torque acting on the vanes is then obtained by integrating the force on each finite element. The resultant net torque is the sum of the torques on each vane, i.e. the cross product between the displacement ~ r and vector sum of the meridian and F tangential forces ~ n X ~ Fi: r i K~

xðBÞ ¼ 100:212 and yðBÞ ¼ 4:602;

In Fig. 10 geometrical position of the actuator and the crank are shown to determine the relationship between crank handle angle and the actuator position.

3. Calculation of aerodynamic force on the actuator

@V u dm; @m @V m dF m ðmÞ ¼ qðmÞV m ðmÞhsðmÞ dm; @m

Let B be the position of the point where the actuator shaft enters the actuator, and where the actuator allows small rotations in the shaft to accommodate the crank handle rotation.

~ L ¼ kBCk;

i;j

~ s¼

443

ð48Þ

i¼1

This torque is transferred to the actuator shaft through the mechanical linkages. In order to find a mathematical relationship for this

y0 ðCÞ ¼ yðCÞ  yðAÞ; y00 ðCÞ ¼ yðCÞ  yðBÞ;

and

d ¼ c þ arctan



 1:433 ; 27:308

where c is the angle of the crank handle measured when the vanes are closed. 3.2. Relation between crank handle angle and lever direction Let h be the crank handle angle, and b be the angle of vane levers. These are shifted, and vary equally:

Fig. 10. Physical positions of primary lever pivot point (A), crank rotation axis (C) actuator position (B).

444

F.S. Ahmed et al. / Energy Conversion and Management 84 (2014) 436–447

Fig. 11. Turbocharger geometry with all dimensions.

Fig. 12. Engine test bench.

b ¼ h  b3 ; where b3 is the direction of levers, measured when the vanes are closed. In our case, the aerodynamic force tends to open the vanes:

b ¼ h  7:456: 3.3. Relation between the levers and unison ring Let us define r2 as the radial position of the point of contact between the unison ring and the ‘‘secondary lever’’ of the crank handle, and r 3 as the radial position of the point of contact between the unison ring and any vane lever. The maximum values of r 2 and r3 are r 20 and r 30 respectively, both obtained when the direction of

the secondary lever is purely radial to the unison ring. We obtain the following relationships:

    b b sin ; 2 2     b b sin : r 3 ¼ r 30  2r 30 tan 2 2

r 2 ¼ r 20  2r 20 tan

3.3.1. Total torque on the rim The total torque on the rim of the unison ring can be found by

T rim ¼ n

r3 cosðu3  a þ hÞ; L3

F.S. Ahmed et al. / Energy Conversion and Management 84 (2014) 436–447

445

Fig. 13. A simulator based on C compiler code to solve 1D flow equations.

with n is number of vanes, a is the displacement angle of the crank handle, L3 is length of vane levers and u3 is the angle of contact between levers and unison ring. In Fig. 11, different dimensions of the turbocharger are shown to determine the angle of each part associated with actuator. 3.3.2. Net torque on the crank handle The net torque on the crank handle can be found by

T crank

Let w be the angle between the shaft and the ‘‘primary lever’’ of the crank handle, i.e. between segments BC and AC:

 w ¼ d  n; d ¼ arctan

Force on actuator: Finally, the force acting on the actuator is obtained as

Lrod ¼

L2 ¼ T rim cosðu2  a þ hÞ; L3

where L2 is the length of the ‘‘secondary lever’’ of the crank handle and u2 is the angle between the secondary lever and the unison ring.

 x0 ðCÞ : 0 y ðCÞ



1000T crank

cosða0  a þ hÞ



: sinðwÞ L1

where a0 is the angle between the axis of the actuator shaft and the normal of the primary lever, when the vanes are open. 4. Model validation

3.4. Angle between the crank handle and actuator shaft Let n be the direction of the shaft:

 00  x ðCÞ : n ¼ arctan 00 y ðCÞ

The engine testing was carried at the ‘‘R&D moteurs’’ facilities with the DV6TED4 engine (4 cylinder, 16 valve, turbocharged diesel). The test bench (Fig. 12) is equipped with GT1544V turbocharger, manufactured by Honeywell GARRETT. The extensive instrumentation allows us to measure many quantities that would

Fig. 14. Comparison of turbine inlet pressure, engine: 2000RPM, 100 Nm.

446

F.S. Ahmed et al. / Energy Conversion and Management 84 (2014) 436–447

Fig. 15. Comparison of turbine inlet pressure, engine: 2000RPM, 150 Nm.

Fig. 16. Comparison of aerodynamic force on actuator, engine: 2000RPM, 100 Nm.

Fig. 17. Comparison of aerodynamic force on actuator, engine: 2000RPM, 150 Nm.

otherwise not be measured in automobiles. For our interest, the force sensor installed on the actuator shaft and the pressure sensor at the turbine inlet is used for validation of the simulator.

Data acquisition was done in two modes: slow acquisition (sampling frequency of 10 Hz) and medium acquisition (sampling frequency of 1000 Hz) on an AVL system. The aerodynamic force

F.S. Ahmed et al. / Energy Conversion and Management 84 (2014) 436–447

estimate was obtained by implementing the calculations discussed about 1D Solver in this article. In this section, comparative results between the estimated and measured force are discussed, along with the turbine inlet pressure measured in the engine and estimated through the 1D Model. A simulation tool was developed in QtCreator, in which, the 1D flow model was implemented in form of C++ class libraries. The graphical user interface of this tool is shown in Fig. 13. It is used to configure the initial conditions of the model, i.e. the initial temperature and pressure at the turbine inlet, its speed and vane opening position. The comparison results between estimated and measured turbine inlet pressure are shown in Figs. 14 and 15. The error is around 20% when losses are not included in the calculations and less than 10% when they are included. It can be seen that the results deviate at larger actuator displacements. This is because when the vane is nearly closed, the small throat area causes high turbulence in the actual airflow, which is not taken into account in the calculations. Figs. 16 and 17 show the comparison of the actuator force at different engine loads. It can be seen that the pressure losses play an important part in the estimation of aerodynamic force at high torque. The addition of losses improves the correlation between the estimated result and the experimental result. At lower torques, the effect of losses is minimal. The minor deviation in these results can be attributed to the distribution of mechanical tolerances between the mechanisms. 5. Conclusion In this article, a 1D turbine airflow model was developed in order to estimate the aerodynamic force exerted by the exhaust gas flow on the VGT vanes. The model was discretized in order to facilitate its implementation, and numerical approximation of its partial derivatives was also discussed. Then, the aerodynamic force was estimated on the VGT vanes and transferred to the actuator shaft, in order to obtain its influence on the actuator. Pressure losses in the turbine caused by its geometry are also considered. The model was used to develop a simulation tool for prediction of the aerodynamic force as well as the internal mass-flow states of the turbine. Experimental validation with pressure losses shows that the tool predicts the aerodynamic force accurately. However, the estimation of turbine inlet pressure requires detailed modeling of turbulent air-flow at narrow apertures. This simulation tool will be used to predict the effect of the aerodynamic force on the air-path control and to develop robust controllers that can successfully counteract it. Future research works will be focused on developing real-time versions of this model, for possible integration in the control directly, leading towards enhancement the overall turbocharger performance and engine efficiency. References [1] Jacobs T, Assanis D, Filipi Z. The impact of exhaust gas recirculation on performance and emissions of a heavy duty diesel engines. Journal Society of Automotive Engineers Inc. SAE paper 2003-01-1068; 2003. [2] Rabih O, Rafic Y, Jean-Claude C, Rachid O. New indicated mean effective pressure (imep) model for predicting crankshaft movement. Energy Convers Manage 2013;52:3376–82.

447

[3] Osaka K, Ibaraki S. Development of the high performance and high reliability vg turbocharger for automotive applications. Mitsubishi Heavy Indust Tech Rev 2006;43(3). [4] Guzzella L, Amstutz A. Control of diesel engine. IEEE Control Syst Mag 1998;18(5):53–71. [5] Rajoo S, Martinez-Botas R. Variable geometry mixed flow turbine for turbocharger: an experimental study. Int J Fluid Machine Syst 2008;1(1):155–68. [6] Tancrez M, Galindo J, Guardiola C, Fajardo P, Varnier O. Turbine adapted maps for turbocharger engine matching. Exp Therm Fluid Sci 2011;35(1):146–53. [7] Ainley DG, Mathieson GCR. An examination of the flow and pressure losses in blade rows of axial flow turbines. Ministry of supply technical memorandum no. 2891. Aeronautical Research Council, London; 1951. [8] Kacker S, Okapuu U. A mean line prediction method for axial flow turbine efficiency. J Eng Power 1982;104:111–9. [9] Setoguchi T, Santhankumar S, Takao M, Kim TH, Kaneko K. A performance study of a radial turbine for wave energy conversion. J Power Energy 2001;216:15–22. [10] Jung M. Mean-value modelling and robust control of the airpath of a turbocharged diesel engine, Phd thesis. University of cambridge, United kingdom; 2003. [11] Hoge G. Exhast pressure modeling for feedforward active control of internal combustion engine exhaust noise. J Automobile Eng Part D 1998;213:109–17. [12] Macian V, Lujan JM, Bermudez V, Guardiola C. Exhaust pressure pulsation observation from turbocharger instantaneous speed measurement. J Measure Sci Technol 2004;15:1185–94. [13] Desantes J, Galindo J, Guardiola C, Dolz V. Air mass flow estimation in turbocharged diesel engines from in-cylinder pressure measurement. Exp Therm Fluid Sci 2010;34(1):37–47. [14] Kolmanovsky I, Moraal P, Nieuwstadt M, Stefanopoulou A. Issues in modelling and control of intake flow in variable geometry turbocharged engine. In: 18th IFIP conf. on sys. modeling and optimization; 1997. p. 436–45. [15] Serrano JR, Arnau FJ, Dolz V, Tiseira A, Cervello C. A model of turbocharger radial turbines appropriate to be used in zero- and one-dimensional gas dynamics codes for internal combustion engines modelling. Energy Convers Manage 2008;49:3729–45. [16] Dunham J, Came PM. Improvements to the Ainley–Mathieson method of turbine performance prediction. ASME J Eng Power 1970;92(3):252–6. [17] Du W, Leonard O. Numerical simulation of surge in axial compressor. Int J Rotat Machine 2012;12:1–10. [18] Leonard O, Adam O. A quasi one dimensional cfd model for multistage turbomachines. J Therm Sci 2007;16(2):1–15. [19] Yee HC. Numerical approximation of boundary conditions with applications to inviscid equations of gas dynamics. NASA technical memorandum 81265, Ames Research Center, Moffett field, California; 1981. [20] Raetz H, Kammeyer J, Natkaniec CK, Seume JR. Numerical investigation of aerodynamic radial and axial impeller forces in a turbocharger. ASME Turbo Expo, Vancouver, Canada; 2011. [21] Jinnai Y, Arimizu H, Tashiro N, Tojo M, Yokoyama T, Hayashi N. A variable geometry turbocharger for passenger cars to meet European union emission regulations. Mitsubishi Heavy Indust Tech Rev 2012;49(2):17–26. [22] Mehmood A, Laghrouche S, El Bagdouri M. Modeling identification and simulation of pneumatic actuator for vgt system. Sens Actuat: Phys 2011;165:367–78. [23] Aungier RH. Centrifugal compressors: a strategy for aerodynamic design and analysis. ASME Press, ISBN-13: 978-0791800935; 2000. [24] Love AC, Agnew G, Szczyrba B. Adaptive variable geometry turbocharger strategy, patent number: US 2009/0123272 A1; 2009. [25] Hiereth H, Prenninger P. Charging the internal combustion engine. SpringerVerlag, ISBN-13: 978-3211330333; 2007. [26] Heywood J. Internal combustion engine fundamentals. McGraw-Hill, ISBN-13: 978-0070286375; 1988. [27] Temam R. Navier–Stokes equations: theory and numerical analysis. American Mathematical Society, ISBN-13: 978-0821827376; 2000. [28] Hoffman KA, Chiang ST. Computational fluid dynamics, vol. 1. Engineering Education System, ISBN-13: 978-0962373107; 2000. [29] Ainley DG, Mathieson GCR. A method of performance estimation for axial flow turbines. Reports and memorandum no. 2974, Aeronautical Research Council, London; 1951. [30] Jennions IK, Stow P. A quasi three dimensional turbomachinery blade design system: part 1 throughflow analysis. ASME J Eng Gas Turb Power 1983;107(2):301–7. [31] Morton KW, Mayers DF. Numerical solution of partial differential equations: an introduction. Cambridge University Press, ISBN-13: 978-0521607933; 2005.