- Email: [email protected]

Contents lists available at SciVerse ScienceDirect

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

Analysis of a directional hydraulic valve by a Direct Numerical Simulation using an immersed-boundary method Antonio Posa a,⇑, Paolo Oresta b, Antonio Lippolis b a b

Department of Mechanical and Aerospace Engineering, The George Washington University, 801 22nd Street, N.W., Washington, DC 20052, USA Dipartimento di Meccanica, Matematica e Management, Politecnico di Bari, Viale Japigia 182, 70126 Bari, Italy

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 16 December 2011 Received in revised form 26 June 2012 Accepted 20 July 2012 Available online 17 October 2012

The improvement of the hydraulic valves depends on the careful analysis of the coherent structures driving the motion of the working ﬂuid. In the past those devices have been studied by experimental tests; during the last 15 years also several numerical works have been presented, solving the ﬂow on body-ﬁtted computational grids by RANS methods. In this study a different approach is proposed for the axisymmetric analysis of a directional valve (4/3, closed center): whereas the RANS techniques are based on the time-averaged equations of the ﬂow, in the present work the unsteady Navier–Stokes equations have been solved using the Direct Numerical Simulation (DNS), which provides important details on the instantaneous structures of the ﬂow, affecting the valve performance. Furthermore, while in the previous numerical studies the computational domain has been discretized by conformal grids, in this case the ﬂuid-body interaction has been represented by an immersed-boundary (IB) method on a Cartesian grid. The analysis of the discharge coefﬁcient and the ﬂow forces for different openings s and pressure drops Dp is presented in this paper. The behavior of those global parameters is justiﬁed also considering the time-averaged and the instantaneous ﬁelds. For small openings and pressure drops the ﬂow is steady and attached to the wall of the discharge chamber on the side of the restricted section. When s and Dp are increased the jet separates at the restricted section and it re-attaches downstream (Coanda effect), keeping the steady state. Finally, for large openings and pressure drops the ﬂow becomes strongly unsteady: it is organized like a free jet and is dominated by large vortices. Ó 2012 Elsevier Ltd. All rights reserved.

Keywords: Direct Numerical Simulation Directional hydraulic valves Finite-difference methods Immersed-boundary methods

1. Introduction In mechanics there are several devices where a ﬂow through a restricted section is produced. Also every hydraulic component is based on the losses generated on a restricted section. It is critical to evaluate accurately those losses and the corresponding discharge coefﬁcient, summarizing their effects: the aim is a careful analysis of the behavior of the device, in order to improve its design. Therefore the literature presents many studies dealing with the ﬂows through restricted sections of hydraulic valves and the global parameters describing them. Originally those studies were based on the experimental work due to Merrit [15]. However those global parameters depend on the time evolution of the coherent structures driving the motion of the working ﬂuid. The topology of these structures is challenging to investigate using experiments, since the measurements are affected by the highly unsteady dynamics ⇑ Corresponding author. Tel.: +1 3018252066. E-mail addresses: [email protected] (P. Oresta), [email protected] (A. Lippolis).

(A.

Posa),

[email protected]

0196-8904/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.enconman.2012.07.012

of the working ﬂuid inside a chamber with a restricted section. Then, thanks to the remarkable developments achieved by the Computational Fluid Dynamics (CFD) during the last years, several analyses on the hydraulic valves have been performed numerically [1–3,5–8,10–13,23]; they allowed to gain a better insight into many phenomena, with signiﬁcant savings in comparison with the experimental approach, above all for the optimization of the geometries, requiring several prototypes when an experimental method is used. However those numerical works have been based on the RANS formulation of the Navier–Stokes equations. This method is not fully satisfactory, since it utilizes a set of parameters to tune a turbulence model: whereas the values of those parameters have been widely analyzed in the literature for the external ﬂows (above all in aeronautics), for the internal ﬂows they are strongly dependent on the considered ﬂow problem, hence their right deﬁnition is much more problematic. This issue implies that an accurate numerical modeling of turbulence in the ﬁeld treated in this study is not trivial. According to the authors of the present paper the improvements of the computational resources and of the numerical

498

A. Posa et al. / Energy Conversion and Management 65 (2013) 497–506

Nomenclature

al, cl, ql coefﬁcients of the time discretization scheme #2 mean angle between the velocity vector and the axial direction on the restricted section of the valve ^ u provisional velocity at the intermediate time level f forcing term in the momentum equation RHS sum of the discretized viscous, convective and pressure gradient terms of the momentum equation uF velocity in the ﬂuid domain u velocity vector VC velocity on the surface of the immersed-boundary velocity boundary condition at the interface points VI Dp pressure drop Dpmax maximum pressure drop for each opening Dt time step m kinematic viscosity of the working ﬂuid / pseudo-pressure used in the projection operation q density of the working ﬂuid discharge coefﬁcient Ce

methods allow a different approach to the simulation of the considered ﬂows: they can be analyzed by the direct solution of the Navier–Stokes equations (DNS), which avoids the time-averaging of the RANS approach and simulates the unsteady evolution of the coherent structures of the ﬂow. Hence no turbulence model has been required for the simulations treated in this paper: this is useful to improve the accuracy of the numerical results, since there are no approximations due to the turbulence modeling. The whole range of spatial and temporal scales has been solved from the smallest dissipative scale up to the largest one, because the grid and time steps used in the present study are small enough to represent the dynamics of the smallest and fastest vortices. Usually the CFD methods discretize the computational domain by body-ﬁtted grids, whose geometry is conformal to the body over which the ﬂow develops: the boundary conditions on the solid object are deﬁned at the nodes of the computational mesh on the surface of the same object. On the contrary, in this work the interaction between the ﬂuid and the body is represented by an immersed-boundary method [14,16]. This allows to use regular grids, Cartesian or cylindrical: the body is ‘‘immersed’’ inside the computational grid and the no-slip boundary conditions are enforced at some interface nodes by a discrete source term in the momentum equation, without the requirement that those points are on the surface of the body. A more detailed description of the IB technique will be provided in Section 3. The IB method decreases the computational cost associated with the grid generation, which becomes absolutely negligible. Moreover, the cells distortion is avoided: the discretization of the differential equations of the ﬂow is simpliﬁed and the numerical accuracy is improved. This is very beneﬁcial, above all if the IB method is utilized with an eddy-resolving technique, as the Direct Numerical Simulation, requiring high accuracy and good conservation properties. Another important feature of the IB approach is the possibility to simulate the ﬂow around moving boundaries without the need to modify the computational grid, since there is no requirement that the grid is conformal to the surface of the bodies. On the contrary, a body-ﬁtted method needs to generate a new mesh at every time instant. This causes not only higher computational costs, as a new grid is created and the solution has to be interpolated; the same interpolation process is also responsible for a reduced accuracy of the solution.

D d F G Gmax K L l Nr Nz p Q Qmax Re s t U

external diameter of the spool internal diameter of the spool axial force on the spool mass ﬂow maximum mass ﬂow for each opening non-dimensional axial force on the spool characteristic length of the ﬂow problem index of the last time level number of grid nodes along the radial direction number of grid nodes along the axial direction ratio pressure/density volumetric ﬂow rate maximum volumetric ﬂow rate for each opening Reynolds number opening of the valve time characteristic velocity of the ﬂow problem

Therefore the body-ﬁtted methods are less accurate and more time consuming, compared with the immersed-boundary ones, above all in the case of ﬂows involving moving bodies. Then in this paper a series of DNS simulations of a directional hydraulic valve (4/3, closed center) is presented; to deﬁne the conditions on the interface ﬂuid-body an immersed-boundary method has been used. The choice of that kind of valve is justiﬁed by the experience in this ﬁeld by the authors [1,2,10–12], by its geometry, which is almost axisymmetric, and by the low Reynolds numbers, due to the small dimensions and the high viscosity values of the mineral oils. Considering that the axial symmetry of the directional valves is almost exact, as it is broken only by the adduction and the discharge chambers, an axisymmetric study has been performed, using an approximation which has been assumed acceptable. This 2D approach has been chosen to decrease the computational time of the simulations and to carry out a wide series of tests. The low Reynolds numbers made it possible to satisfy the stability conditions for the direct solution of the Navier–Stokes equations with a computational effort small enough. A code originally developed by Verzicco et al. [22] has been used. This paper is organized as follows: further details about the numerical method are provided in Section 2, while the immersed-boundary approach is discussed more carefully in Section 3; the valve geometry, the computational grid and the simulated openings and ﬂow rates are described in Section 4; about the results of the simulations, in Section 5 the dependence of the global parameters on the opening and the pressure drop is analyzed, whereas in Section 6 and in Section 7 the effects of the time-averaged and the instantaneous ﬂow topology on the global parameters are treated. Finally, in Section 8 the conclusions of the present study are reported.

2. Numerical method The analysis of the valve has been performed by the numerical solution of the incompressible Navier–Stokes equations. In this case the continuity and the momentum equations are:

ru¼0

ð1Þ

@u 1 þ r uu ¼ rp þ r2 u þ f @t Re

ð2Þ

A. Posa et al. / Energy Conversion and Management 65 (2013) 497–506

In Eqs. (1) and (2) u is the velocity vector, t the time variable, p the ratio between pressure and density, f the forcing term due to the immersed-boundary and Re = UL/m the Reynolds number, based on the kinematic viscosity of the ﬂuid m, a reference length L (the opening s of the valve) and a reference velocity U (the mean radial velocity on the restricted section). The velocity U has been evaluated as the ratio between the ﬂow rate and the cylindrical area of the restricted section. For large openings of the valve the axial component of the velocity is not negligible, but it does not change the order of magnitude of the Reynolds number and its physical meaning. In this work the openings are in the range 0.1–2.0 mm and the kinematic viscosity is equal to 50 cSt, then the Reynolds number is between 50 (smallest opening and ﬂow rate) and 3000 (largest opening and ﬂow rate). These values are quite small, allowing a direct numerical simulation of the ﬂow inside the valve. The Navier– Stokes equations have been solved without time-averaging, used by the RANS methods. Therefore the time evolution of the ﬂow has been simulated and its statistics have not been evaluated by averaged equations, but using the instantaneous ﬁelds of the DNS solution. The differential equations have been discretized in time by a fractional-step method: the viscous terms have been treated implicitly by a Crank–Nicolson scheme, while for the convective terms an explicit Adams–Bashfort scheme has been used. The accuracy is second order in time. All the spatial derivatives have been approximated by second-order central ﬁnite-differences on a staggered grid. This criterion of discretization avoids the numerical dissipation: diffusion is only due to physical terms, the viscous ones. This conservation property is crucial for unsteady methods, as DNS, which simulate the dynamics of the ﬂow, improving the accuracy of the solution in comparison with RANS. The discretized momentum equation at the intermediate time level l + 1/2 is provided in the following equation:

499

3. The immersed-boundary method A brief description of the immersed-boundary method has been already presented in the Introduction and more particulars will be provided in this section. In this method the body is immersed into a regular grid, Cartesian or cylindrical, whose nodes are tagged as interior (inside the immersed-boundary), exterior (in the ﬂuid domain) and interface points. In the present implementation of the IB method the interface nodes are the ones in the ﬂuid domain adjacent to one or more interior points. The different types of grid nodes have been represented in Fig. 1 for a body with a circular shape. The no-slip boundary condition is enforced at the interface points by a discrete source term f in the momentum equation (Eq. (2)). This forcing term is evaluated according to Eq. (7), as proposed by Mohd-Yusof [17].

f

lþ12

¼

1 V I ul RHS lþ2 Dt

ð7Þ

Then f at the intermediate time level is established considering the velocity u at the previous time level, the velocity boundary 1 condition VI at the interface points and the term RHS lþ2 . This is the right hand side of Eq. (3), excluding f; then it is composed of the discretized viscous, convective and pressure gradient terms. Using the method proposed by Balaras [4], the velocity VI at the interface point I (Fig. 2) has been evaluated by a linear interpolation,

^ ul al 2 ^ u r ðu ¼ al rpl ½cl r ðuuÞl þ ql r ðuuÞl1 þ Dt 2Re þ ul Þ þ f

lþ12

ð3Þ

^ is a provisional non-solenoidal velocity, l the index In Eq. (3) u of the previous time level, Dt the time step; al, ql and cl are the coefﬁcients of the numerical scheme. ^ does not satisfy the The velocity at the intermediate time level u continuity. Then, it is projected onto a divergence-free vector ﬁeld by Eq. (4).

^ al Dt r/ ulþ1 ¼ u

ð4Þ

Fig. 1. Tagging of the grid points as: : interior; h: interface; : exterior.

The continuity condition implies that the divergence of the velocity ﬁeld at the time level l + 1 must be equal to zero. Then, this leads to Eq. (5).

r2 / ¼

^ ru al Dt

ð5Þ

In the present study the discretized form of the elliptic Eq. (5) corresponds to a penta-diagonal problem, since the simulations have been performed on a 2D domain. This problem has been solved by a generalized cyclic reduction algorithm [19], using the FISHPACK package. Finally, substituting Eq. (4) into Eq. (3), it is possible to update the pressure ﬁeld for the time level l + 1, according to Eq. (6).

plþ1 ¼ pl þ /

a l Dt 2Re

r2 /

ð6Þ

More details on the numerical solution of the Navier–Stokes equations are reported in [21].

Fig. 2. Evaluation of the velocity VI at the interface points: linear interpolation between the ﬂuid domain and the immersed-boundary.

500

A. Posa et al. / Energy Conversion and Management 65 (2013) 497–506

along the normal to the body, between the ﬂuid point F and the surface of the immersed-boundary, where the velocities are respectively uF and VC. The point F has been placed one grid cell from the interface node I, along the outward direction normal to the surface of the body. As the grid geometry is not conformal to the shape of the immersed-boundary, in general F is not a grid point; therefore the velocity uF has been estimated by a bilinear interpolation, using the closest surrounding nodes of the computational grid (P1, P2, P3, P4 in Fig. 2). It is interesting to notice that the evaluation of the forcing term according to Eq. (7) is equivalent to enforce in Eq. (3) the velocity boundary condition VI at the interface points. This is evident ^ ¼ VI. substituting Eq. (7) into Eq. (3), giving u Further details on the immersed-boundary technique can be found in [14,16]. Recent works using this method are discussed in [9,18,20].

Fig. 3. Meridian section of a 4/3 closed center directional valve. The body valve is represented in dark gray, the spool in light gray, the simulated domain in blue. The axis of symmetry is shown by ——. The black arrows indicate the direction of the ﬂow through the metering sections, the red one the shift of the spool during the opening of the valve.

4. Computational set up In Fig. 3 a meridian section of a 4/3 closed center directional valve is shown during its opening. The shift of the spool (light gray) inside the body valve (dark gray) opens the connection ports between P and A and between B and T respectively. The high pressure oil entering P goes to the actuator through the metering section PA, deﬁned by the edges of the body valve and the spool. The ﬂuid coming from the actuator ﬂows to the reservoir through the metering section B-T. In Fig. 3 the blue1 area represents the simulated domain. The geometry of the valve is axisymmetric, except for the adduction and the discharge connections, that are not distributed along the whole azimuthal direction: they are located at a particular angular position, while the rest of the chambers is bounded radially by the walls of the body valve, as shown in Fig. 3. Then the numerical simulations have been performed using axisymmetric conditions on a meridian section. Thus the circumferential ﬂows, related to the speciﬁc positions of the inﬂow and the outﬂow connections, have not been taken into account in the present study. In Fig. 4 the detail of the simulated domain is provided, with its dimensions and the locations of the inﬂow/outﬂow sections. This ﬁgure shows also the axial opening s, deﬁned by the position of the spool respect to the body valve. It is useful to notice that the immersed-boundary approach is especially suitable for the present geometry, which is aligned with the grid lines and makes it easy to cluster points near the walls and on the restricted section. The computational domain has been discretized using a grid with Nr = 359 and Nz = 419 points respectively along the radial (vertical) and the axial (horizontal) directions. For clarity in Fig. 5 the grid is plotted representing only one point every four. Fig. 5 shows that the regular grid is ﬁner in the area where larger gradients occur, i.e., on the restricted section and in the discharge chamber. In fact, the stretching is strong upstream the restricted section, but it is much more limited downstream. The minimum step, along both the radial and the axial directions, is equal to 0.02 mm. Several tests have been carried out using lower and higher resolution levels. Some comparisons with the ﬁelds on a ﬁner computational grid for the maximum pressure drop and the mean opening have been reported in Section 6 to prove the accuracy of the results presented in this paper. In this numerical study four different ﬂow rates for ﬁve openings of the valve have been treated. The maximum volumetric ﬂow 1 For interpretation of colors in Fig. 3, and in the following ones the reader is referred to the web version of this article.

Fig. 4. Simulated geometry. Its dimensions are provided in mm.

rate Q for s = 1.0 mm has been chosen on the basis of the value reported in [12]. The other three have been estabpﬃﬃﬃ levels of ﬂowprate ﬃﬃﬃ lished using the factors 1= 2, 1/2 and 1= 8. Then, the reduced values correspond approximately to 70%, 50%, and 35% of the original one. At the other openings of the valve the ﬂow rates have been set considering a linear dependence on s. A summary of the simulated openings s and mass ﬂows G is given in Table 1. It is important to highlight that in Table 1 the maximum mass ﬂow Gmax is a function of the opening. Furthermore, the density q of the working ﬂuid, relating G and Q, has been established equal to 850 kg/m3. The computational time is strongly different for steady and unsteady ﬂows. As it will be shown later, the steady case pertains to the smallest opening for every ﬂow rate and the opening s = 0.5 mm for G = 50%Gmax and G = 35%Gmax. All the other simulations present an unsteady solution. For them the computational time has been increased to have accurate statistics of the ﬂow.

Fig. 5. Computational grid. Only 1 point every 4 is shown along both directions.

A. Posa et al. / Energy Conversion and Management 65 (2013) 497–506

501

Table 1 Mass ﬂows in kg/s for the simulated openings. s (mm)

100%Gmax

70%Gmax

50%Gmax

35%Gmax

0.1 0.5 1.0 1.5 2.0

0.320 1.600 3.200 4.800 6.400

0.226 1.131 2.263 3.394 4.525

0.160 0.800 1.600 2.400 3.200

0.113 0.566 1.131 1.697 2.263

The computations have been performed by a code parallelized in OpenMP: four processors have been used for each simulation. In the case of steady solution the mean computational time has been 11 h 46 min, while for the other simulations the mean has been 21 h 59 min. 5. Global parameters In this section the behavior of the main global parameters for various openings and ﬂow rates is described. These parameters, as the discharge coefﬁcient Ce and the axial force on the spool F, or its corresponding non-dimensional parameter K, are critical in the valve design. Here for the simulations with an unsteady solution their time-averaged values are reported. Such parameters are able to provide an overall information about the performance of the valve. Therefore their role is central in the valve design. The discharge coefﬁcient is deﬁned as the ratio between the actual ﬂow rate and the theoretical one, as follows:

Ce ¼

Q qﬃﬃﬃﬃﬃﬃ pDs 2Dqp

ð8Þ

In Eq. (8) D is the external diameter of the spool. The other quantities have been already deﬁned above. The theoretical ﬂow rate is deﬁned by the Borda hypothesis, according to which the ﬂow is isentropic up to the restricted section and isobaric downstream. Then, the discharge coefﬁcient is different from 1 because of four main reasons: (1) the section of the vena contracta is smaller than the restricted section of the valve; (2) the losses upstream the restricted section; the ﬂow is not effectively isentropic; (3) the effect of the boundary layer; (4) the pressure recovery downstream the restricted section. The force parameter K represents the non-dimensional axial force on the spool. It is deﬁned as follows:

K¼

F 2

2

p D d Dp 4

ð9Þ

where d is the internal diameter of the spool. Fig. 6 shows the trend of Dp as a function of the mass ﬂow for each opening s. The function Dp = Dp(G) is almost parabolic for all

Fig. 6. Evolution of the pressure drop for variable mass ﬂows. Also the parabolic regression lines are plotted. : s = 0.1 mm; : s = 0.5 mm; : s = 1.0 mm; : s = 1.5 mm; : s = 2.0 mm.

Fig. 7. Discharge coefﬁcient Ce as a function of the valve opening for different ﬂow rates. : Q = 100%Qmax (Dp = 100%Dpmax); : Q = 70%Qmax (Dp = 50%Dpmax); : Q = 50%Qmax (Dp = 25%Dpmax); : Q = 35%Qmax (Dp = 12.5%Dpmax).

the openings of the valve, because of the small variations of the discharge coefﬁcient, as it will be shown later. The performance of the valve is described better by the discharge coefﬁcient: Fig. 7 provides Ce as a function of the valve opening for different pressure drops. For large openings it is substantially not dependent on the pressure drop. On the contrary, for small openings the pressure drop affects more the discharge coefﬁcient. When s = 0.1 mm the role of the boundary layer is signiﬁcant and variable for different pressure drops; larger values of Dp cause higher velocities on the restricted section and decrease the inﬂuence of the boundary layer on the ﬂow rate. Therefore at the smallest simulated opening the highest Ce is associated with the largest Dp. For s = 0.5 mm the behavior is modiﬁed. This is due to the Coanda effect. When the pressure drop and the ﬂow rate are small the Coanda effect occurs: the jet re-attaches downstream the restricted section on the left wall of the discharge chamber. When Dp and Q are increased the Coanda effect tends to decay: in the discharge chamber the ﬂow becomes strongly unsteady and large vortices are generated, causing signiﬁcant losses and smaller values of Ce. This justiﬁes the inversion of the inﬂuence of the pressure drop on the discharge coefﬁcient from s = 0.1 mm to s = 0.5 mm. Fig. 7 shows also that the discharge coefﬁcient is decreasing with increasing openings of the valve; this trend is due to the variation of the angle of the jet on the restricted section. In fact, the angle #2 between the mean velocity vector on the restricted section and the axial direction is smaller for larger openings. This implies an increasing difference between the section of the vena contracta and the restricted section, reducing the value of Ce. However it is interesting to highlight that the evolution shown in Fig. 7 is not monotonic: a maximum for s = 0.5 mm and small pressure drops is produced. It can be explained considering the Coanda effect. It is not able to change strongly the discharge angle, but it decreases the losses downstream the restricted section in comparison with the ones generated by larger ﬂow rates and valve openings. Furthermore from s = 0.1 mm to s = 0.5 mm the inﬂuence of the boundary layer on the restricted section is reduced: the effect on the discharge coefﬁcient is beneﬁcial. Those phenomena imply the production of a local maximum of Ce at about s = 0.5 mm. The above considerations about the discharge angle are veriﬁed by Fig. 8. The dependence of #2 on the valve opening is evident. The discharge angle decreases with increasing openings. Fig. 8 proves also the negligible inﬂuence of the pressure drop on the discharge angle. It is possible to notice a slight variation of #2 only when s = 0.5 mm; for the smallest pressure drop the Coanda effect occurs, but it is almost vanished for the largest one. It is able to increase the discharge angle, but this phenomenon is not very evident, although a change is found in comparison with the

502

A. Posa et al. / Energy Conversion and Management 65 (2013) 497–506

Fig. 8. Discharge angle #2 as a function of the pressure drop for different openings of the valve.

behavior for the other openings of the valve. For larger values of s the Coanda effect is not veriﬁed, while for the smallest one the jet is always attached to the left wall of the outﬂow chamber, without separation at the restricted section and reattachment downstream. Fig. 9 represents also the dependence of the discharge coefﬁcient on the pressure drop. That ﬁgure allows to complete the description of the phenomena which affect Ce. For the largest simulated openings (s = 1.0 mm, s = 1.5 mm and s = 2.0 mm) the trend is almost constant. In those conditions the Coanda effect does not occur, the inﬂuence of the boundary layer is negligible, since those openings are large, and for variable values of Dp the discharge angle does not change, as proved by Fig. 8. The phenomena at the smallest openings are more interesting. For s = 0.5 mm Ce decreases with increasing pressure drops. This depends on the decay of the Coanda effect: the discharge angle slightly decreases and, above all, the losses in the outﬂow chamber are enhanced. Differently the slope of Ce for s = 0.1 mm is positive. This is due to the reduced inﬂuence of the boundary layer when the ﬂow rate is increased. In this case for every pressure drop the jet is attached, without the production of dissipative structures in the discharge chamber, and the angle #2 is constant. Fig. 10 highlights the linear evolution of the axial force on the spool as a function of the pressure drop. As proved in [12], F is only dependent on the valve opening and the pressure drop, assuming that the viscous force on the body valve is negligible, the inﬂow is radial and Ce and #2 are constant. Then the results of Fig. 10 are consistent, showing the accuracy of the simulations: Figs. 8 and 9 proved that the hypothesis of Ce and #2 independent on Dp is substantially correct, above all for large values of s; therefore for each opening of the valve F is roughly a linear function of Dp. In Fig. 11 the inﬂuence of the opening of the valve on the axial force on the spool is presented, using the non-dimensional parameter K. This allows to eliminate the dependence of F on the pressure

Fig. 10. Axial force on the spool as a function of the pressure drop for different openings of the valve. Also the linear regression lines are plotted. : s = 0.1 mm; : s = 0.5 mm; : s = 1.0 mm; : s = 1.5 mm; : s = 2.0 mm.

Fig. 11. Non-dimensional axial force on the spool as a function of the valve opening for different pressure drops. Also the linear regression line is plotted. : Q = 100%Qmax (Dp = 100%Dpmax); : Q = 70%Qmax (Dp = 50%Dpmax); : Q = 50%Qmax (Dp = 25%Dpmax); : Q = 35%Qmax (Dp = 12.5%Dpmax).

drop. The evolution of K is almost proportional to s, but there are some deviations from the linear regression line shown in Fig. 11. Those occur above all at the opening s = 0.5 mm for small pressure drops, since in those conditions Ce and #2 present more signiﬁcant variations. In fact, the Coanda effect is responsible for an increase of the discharge angle and the discharge coefﬁcient with decreasing Dp (Figs. 8 and 9). The ﬁrst phenomenon tends to reduce the values of F and K; the effect of the second one is opposite. However the increase of Ce is more substantial than the one of #2 . Hence for s = 0.5 mm and small pressure drops the trend of K is more than linear. In Table 2 the maximum values of the time-averaged pressure drop between the inﬂow and the restricted section have been provided, compared with the ones between the inﬂow and the outﬂow, reported in parentheses. They are useful to prove that cavitation cannot occur for large openings, since the pressure values on the restricted section are higher than the ones at the outﬂow, because of signiﬁcant dissipative phenomena in the discharge chamber. For small openings the pressure on the re-

Table 2 Maximum time-averaged pressure drops between the inﬂow and the restricted section for each simulated opening and ﬂow rate. The pressure drops between the inﬂow and the outﬂow are reported in parentheses. The values are provided in MPa.

Fig. 9. Discharge coefﬁcient as a function of the pressure drop for different openings of the valve.

s (mm)

100%Gmax

70%Gmax

50%Gmax

35%Gmax

0.1 0.5 1.0 1.5 2.0

5.726(4.689) 5.200(5.594) 5.518(6.235) 6.028(6.961) 6.454(7.992)

2.909(2.416) 2.615(2.375) 2.806(3.080) 3.025(3.559) 3.178(4.016)

1.491(1.261) 1.328(1.097) 1.412(1.546) 1.539(1.784) 1.648(1.955)

0.774(0.671) 0.697(0.502) 0.711(0.735) 0.783(0.916) 0.833(0.966)

A. Posa et al. / Energy Conversion and Management 65 (2013) 497–506

stricted section is lower than the one at the outlet, but, considering the typical values of the outﬂow pressure for industrial valves, it is possible to assume that no cavitation phenomena occur in the working conditions treated in the present paper. 6. Time-averaged ﬂow ﬁelds The time-averaged velocity ﬁelds are plotted for Dp = 12.5%Dpmax and Dp = 100%Dpmax, considering the ﬁve simulated openings of the valve (Fig. 12). In the ﬁrst case (s = 0.1 mm) the ﬂow is steady and attached to the left wall of the discharge chamber (a), with high outﬂow angles #2 (about 75°) and a small axial force on the spool. The small ﬂow rate implies low values of velocity over the whole computational domain, except for the restricted section, where a peak is generated. The wall jet boundary layer makes the ﬂow organized in a coherent ligament without mixing phenomena responsible for energy losses. This behavior increases the discharge coefﬁcient (Ce 0.7) in comparison with its value at the other openings (for instance, Ce 0.55 for s = 2.0 mm). Those ﬂow conditions occur with both the smallest and the largest simulated ﬂow rates (a1 and a2). With increasing openings s the jet separates at the restricted section and it re-attaches downstream (Coanda effect), keeping the steady state for small values of ﬂow rate and pressure drop (b1); when Dp and Q are increased the Coanda effect tends to

Fig. 12. Time-averaged velocity ﬁelds. a: s = 0.1 mm; b: s = 0.5 mm; c: s = 1.0 mm; d: s = 1.5 mm; e: s = 2.0 mm. 1: Q = 35%Qmax (Dp = 12.5%Dpmax); 2: Q = 100%Qmax (Dp = 100%Dpmax). The scales range from 0 (blue) to 38 (red) m/s on the left side and from 0 (blue) to 110 (red) m/s on the right side. (For interpretation of the references to colour in this ﬁgure legend, the reader is referred to the web version of this article.)

503

decay and the ﬂow becomes unsteady, spreading the momentum of the jet in the discharge chamber (b2). Then the Coanda effect occurs with low values of valve opening s and pressure drop Dp. For increasing Dp the tail of the jet ligament detaches from the side wall, producing a mixing layer along the downstream region. Some energy is lost in the mixing region and the discharge coefﬁcient decreases; for instance, when s = 0.5 mm, Ce = 0.77 for the minimum Dp and Ce = 0.65 for the maximum Dp. In the second condition the ﬂow is unsteady and the Coanda effect is not veriﬁed at every time instant. Fig. 12 conﬁrms also that the inﬂuence of the pressure drop on the angle #2 at the restricted section is negligible: even the Coanda effect is not able to affect substantially that angle; in b1 #2 is only slightly larger than in b2. For s = 0.5 mm the discharge angle is strongly different from p/2; this causes a signiﬁcant axial component of the force on the spool. In the case of large s and Dp the ﬂow is organized like a free jet: it becomes statistically more symmetric, reducing its interaction with the left wall of the outﬂow chamber. This is already evident when s = 1.0 mm (c). However at the largest opening (e) the inﬂuence of the right wall is not negligible. For the largest pressure drop (e2) it is also possible to notice the production of some structures upstream the restricted section, close to the edges of the body valve and, above all, of the spool. The angle #2 and the discharge coefﬁcient Ce are decreased in comparison with the ones at the small openings. The trend of Ce with increasing s is due to the one of the discharge angle, but also to the wide mixing regions where the energy is dissipated and whose generation is enhanced by high values of s and Q, as proved by the time-averaged ﬁelds in Fig. 12c–e. In Fig. 13 the time-averaged velocity ﬁeld provided by a simulation on a ﬁner grid is shown, to prove the accuracy of the results presented in this paper using the grid described in Section 4. The enhanced computational mesh of Fig. 13 is twice ﬁner than the one of Fig. 12, both along the radial and the axial directions. The ﬁeld of Fig. 13 refers to s = 1.0 mm and Q = 100%Qmax, then it has to be compared with the one of Fig. 12c2. The agreement is satisfactory: upstream the restricted section, in the area of steady ﬂow, which is the one computationally less demanding, no difference can be noticed; downstream the restricted section only minor variations are shown at the center of the discharge chamber from the coarser to the ﬁner grid. In Fig. 14 also the comparison between the ﬁelds of the turbulent kinetic energy is proposed for the same opening and ﬂow rate of Fig. 13. Only the detail of the discharge chamber is represented, since the turbulent kinetic energy upstream the restricted section is negligible. Even in Fig. 14 the agreement between the results on the coarser and the ﬁner grids is good. In both cases the highest levels of turbulence are generated in the areas where the jet interacts with the ﬂuid inside the outﬂow chamber, causing the unsteady production of vortices. Only very small differences can be observed in Fig. 14 between a and b: on the coarser grid the values of turbulent

Fig. 13. Time-averaged velocity ﬁeld for s = 1.0 mm and Q = 100%Qmax (Dp = 100%Dpmax) on a computational grid composed of 717 837 nodes respectively along the radial (vertical) and the axial (horizontal) directions. The scale ranges from 0 (blue) to 110 (red) m/s. (For interpretation of the references to colour in this ﬁgure legend, the reader is referred to the web version of this article.)

504

A. Posa et al. / Energy Conversion and Management 65 (2013) 497–506

Fig. 14. Turbulent kinetic energy for s = 1.0 mm and Q = 100%Qmax (Dp = 100%Dpmax) in the discharge chamber. The scale ranges from 0 (blue) to 2500 (red) m2/s2. a: grid 359 419; b: grid 717 837. The local origin of the axes is at the left edge of the restricted section. The scale of the tics is in mm. (For interpretation of the references to colour in this ﬁgure legend, the reader is referred to the web version of this article.)

kinetic energy are slightly lower near the outﬂow, because of the stretching of the computational mesh. 7. Instantaneous ﬂow ﬁelds The results treated in this paper refer to eddy-resolving simulations (DNS); the instantaneous ﬂow dynamics has been simulated, differently from the usual RANS approaches, which solve the timeaveraged Navier–Stokes equations. No model has been used to represent the turbulence of the ﬂow: all the scales have been solved. This method is useful not only to provide several details on the instantaneous dynamics, giving a better insight into the ﬂow properties. A DNS is also able to improve the estimate of the statistics, since they can be evaluated directly from the instantaneous results, whose accuracy is based on computational codes with optimal conservation properties, as in the present study. In this section some instantaneous ﬁelds will be presented, to describe the unsteady nature of the ﬂow for large values of valve opening and pressure drop. The time evolutions of the axial force on the spool for both the critical opening (s = 0.5 mm, for which the ﬂow becomes unsteady) and the maximum one are shown in Fig. 15, where the results with the maximum Dp are considered. In this study the unsteadiness of

Fig. 15. Time-history of the axial force acting on the spool for Q = 100%Qmax (Dp = 100%Dpmax) and two different openings: s = 0.5 mm ( ) and s = 2.0 mm ( ). The dashed lines refer to the time-averaged values.

the ﬂow occurs starting from s = 0.5 mm, when the ﬂow rate and the pressure drop are large enough; for the smallest simulated opening the ﬂow is always steady. The vorticity ﬁelds are used to highlight the different ﬂow regimes at the opening s = 0.5 mm for the minimum and the maximum values of the pressure drop (Fig. 16). The physics is almost steady upstream the restricted section, then the most interesting area of the instantaneous ﬁelds is downstream, the one represented in Fig. 16. In a, for the minimum Dp, corresponding to the smallest ﬂow rate, the vorticity ﬁeld is characterized by a steady coherent structure on the left wall of the outﬂow chamber (Coanda effect), while in b the ﬂow is unsteady and it is governed by large vortices populating the domain downstream the restricted section; the Coanda effect is not completely vanished, but the jet is not able to stay attached along the whole radial extent of the left wall. The unsteadiness of the ﬂow is much more intense for s = 2.0 mm than for s = 0.5 mm. For the largest opening the production of unsteady structures downstream the restricted section is stronger, as proved in Fig. 17, where two instantaneous vorticity ﬁelds are presented. As in Fig. 16 they refer respectively to the simulations performed with the minimum and the maximum ﬂow rates. Therefore the standard deviation of the force evolution in time is much higher for larger openings, as shown in Fig. 15. It is useful to point out that an accurate estimate of the peak values of the force on the spool is crucial for the proper design of the driving system of the valve. Only an unsteady simulation can provide this kind of information, since the maximum force depends on the time evolution of the ﬂow. In order to explain the strong time-dependence of F, also the instantaneous velocity ﬁelds of Fig. 18 are represented for both s = 0.5 mm and s = 2.0 mm. They are associated with minima (on the left) and maxima (on the right) of F. Those ﬁelds pertain to the maximum pressure drop, as Fig. 15, and prove that the time evolution of the axial force on the spool is due to signiﬁcant variations in time of the discharge angle. The minima of F occur for the largest angles #2 , for which the axial component of the momentum ﬂow on the restricted section is decreased, whereas the maxima are caused by the smallest discharge angles. The increased unsteadiness of the ﬂow, related to higher values of ﬂow rate and pressure drop, is shown better in Figs. 19 and 20. They provide some instantaneous velocity ﬁelds with s = 1.0 mm respectively for the maximum and the minimum simulated pressure drops. Fig. 19 emphasizes the strong time-dependence of the ﬂow. The unsteadiness of the discharge angle is evident; #2 decreases from a

Fig. 16. Instantaneous vorticity ﬁelds in the outﬂow chamber for a valve opening s = 0.5 mm. a: Q = 35%Qmax (Dp = 12.5%Dpmax); b: Q = 100%Qmax (Dp = 100%Dpmax). The scale ranges from 105 (blue) to +105 (red) 1/s. (For interpretation of the references to colour in this ﬁgure legend, the reader is referred to the web version of this article.)

A. Posa et al. / Energy Conversion and Management 65 (2013) 497–506

505

Fig. 17. Instantaneous vorticity ﬁelds in the outﬂow chamber for a valve opening s = 2.0 mm. a: Q = 35%Qmax (Dp = 12.5%Dpmax); b: Q = 100%Qmax (Dp = 100%Dpmax). The scale ranges from 105 (blue) to +10 5 (red) 1/s. (For interpretation of the references to colour in this ﬁgure legend, the reader is referred to the web version of this article.)

Fig. 19. Instantaneous velocity ﬁelds in the outﬂow chamber for s = 1.0 mm and Q = 100%Qmax (Dp = 100%Dpmax). The scale ranges from 0 (blue) to 110 (red) m/s. (For interpretation of the references to colour in this ﬁgure legend, the reader is referred to the web version of this article.)

Fig. 18. Instantaneous velocity ﬁelds in the outﬂow chamber for Q = 100%Qmax (Dp = 100%Dpmax). a: s = 0.5 mm; b: s = 2.0 mm. The scales range from 0 (blue) to 110 (red) m/s in a and from 0 (blue) to 140 (red) m/s in b. The ﬁelds on the left refer to time instants of force minima, the ones on the right to time instants of force maxima. (For interpretation of the references to colour in this ﬁgure legend, the reader is referred to the web version of this article.)

to d: the jet from the restricted section becomes more axial. Furthermore, the local velocity minima, corresponding to blue dots in the velocity ﬁelds of Fig. 19, are the centers of vortical structures: it is shown that those vortices affect the whole discharge chamber for the considered valve opening and pressure drop. Also Fig. 20 refers to the case s = 1.0 mm, but Dp = 12.5%Dpmax. The same time step of Fig. 19 has been considered to plot the instantaneous ﬁelds. It is shown that the unsteadiness is substantially decreased in comparison with Fig. 19. The time-dependence

Fig. 20. Instantaneous velocity ﬁelds in the outﬂow chamber for s = 1.0 mm and Q = 35%Qmax (Dp = 12.5%Dpmax). The scale ranges from 0 (blue) to 38 (red) m/s. (For interpretation of the references to colour in this ﬁgure legend, the reader is referred to the web version of this article.)

506

A. Posa et al. / Energy Conversion and Management 65 (2013) 497–506

of the discharge angle is reduced and the production of coherent structures is much less intense: in Fig. 20 there are only two vortices in the discharge chamber, respectively on the left and on the right sides of the jet; they occur because of the viscous interaction of the same jet with the stagnating ﬂuid. The stagnation areas are larger than in Fig. 19 and are located above all on the right side of the chamber; the convection of the vortices in Fig. 20 is slower than in Fig. 19, because of the lower velocity values.

8. Conclusions In this paper the ﬂow through a directional valve (4/3, closed center) has been studied. The analysis of the discharge coefﬁcient, the discharge angle and the ﬂow forces for different openings s and pressure drops Dp has been provided. The ﬂow has been simulated in a 2D computational domain using the Direct Numerical Simulation: all the scales of the ﬂow have been solved. The interaction between the ﬂuid and the valve has been modeled by an immersed-boundary technique. In the case of strongly unsteady ﬂows, as the one analyzed in this study, the Direct Numerical Simulation is able to describe carefully the time-dependent coherent structures, avoiding the numerical dissipation and preserving the energy budget. This is a crucial feature, since by a DNS the statistics are more accurate in comparison with other techniques, based on turbulence models and dissipative numerical methods, as the Reynolds Averaged Navier– Stokes (RANS) approaches. The results of the simulations showed that the axial force on the spool is roughly proportional to the pressure drop and to the valve opening. The discharge coefﬁcient depends on the opening s, whereas the pressure drop affects its value only in a narrow range close to s = 0.5 mm. The inﬂuence of the pressure drop on the discharge angle is almost negligible, while increasing openings of the valve are responsible for a signiﬁcant decrease of #2 . The behavior of these global parameters has been explained on the basis of the instantaneous and the time-averaged ﬂow ﬁelds. For small openings and pressure drops the ﬂow is steady and attached to the wall of the outﬂow chamber on the side of the restricted section. For larger values of s and Dp the jet separates at the restricted section and it re-attaches downstream (Coanda effect), keeping the steady state. Finally, in the case of large openings and pressure drops the ﬂow is strongly unsteady and organized like a free jet; it is dominated by large coherent structures in the whole discharge chamber. For large values of s and Dp the time evolution of the force on the spool is characterized by a considerable standard deviation: for instance, when s = 2.0 mm and Dp = 100%Dpmax the maximum axial force is 40% higher than the time-averaged one. That peak is critical for the proper design of the driving system of the valve and is dependent on particular instantaneous conditions, for which the discharge angle is decreased in comparison with its average. Therefore the standard deviation and the highest value of F can be reduced controlling the time evolution of the ﬂow in the discharge chamber.

Acknowledgment The authors are grateful to CASPUR (Consorzio interuniversitario per le Applicazioni di Supercalcolo Per Università e Ricerca) for providing computational resources. References [1] Amirante R, Vescovo Del G, Lippolis A. A ﬂow forces analysis on an open center hydraulic directional control valve sliding spool. Energy Convers Manag 2006;47(1):114–31. [2] Amirante R, Del Vescovo G, Lippolis A. Evaluation of the ﬂow forces on an open-centre directional control valve by means of a computational ﬂuid dynamic analysis. Energy Convers Manag 2006;47(13–14):1748–60. [3] Amirante R, Moscatelli PG, Catalano LA. Evaluation of the ﬂow forces on a direct (single stage) proportional valve by means of a computational ﬂuid dynamic analysis. Energy Convers Manag 2007;48(3):942–53. [4] Balaras E. Modeling complex boundaries using an external force ﬁeld on ﬁxed Cartesian grids in large-eddy simulations. Comput Fluids 2004;33(3):375–404. [5] Borghi M, Milani M, Paoluzzi R. Inﬂuence of notch shape and number of notches on the metering characteristics of hydraulic spool valves. Int J Fluid Power 2005;6(2):5–18. [6] Borghi M, Milani M, Paltrinieri F. The inﬂuence of the notch shape and number on proportional directional control valve metering characteristics. SAE Trans – J Comm Veh 2005;113:147–56. [7] Bottazzi D, Franzoni F, Milani M, Montorsi L. Metering characteristics of a closed center load-sensing proportional control valve. SAE Paper – Int J Comm Veh 2010;2(2):66–74. [8] Bottazzi D, Franzoni F, Milani M, Montorsi L. Analysis of a hydraulic valve by means of a transient multidimensional CFD approach. In Proceedings of the 7th international ﬂuid power conference. Aachen (Germany), 22–24 March 2010. p. 409–21. [9] de Tullio MD, Nam J, Pascazio G, Balaras E, Verzicco R. Computational prediction of mechanical hemolysis in aortic valved prostheses. Euro J Mech – B/Fluids 2012;35:47–53. [10] Del Vescovo G, Lippolis A. Three-dimensional analysis of ﬂow forces on directional control valves. Int J Fluid Power 2003;4(2). [11] Del Vescovo G, Lippolis A. CFD analysis and optimization of a compensated spool valve geometry. In: Proceedings of the 2nd international conference on computational methods in ﬂuid power technology. Aalborg (Denmark), 2–3 August 2006. [12] Del Vescovo G, Lippolis A. A review analysis of unsteady forces in hydraulic valves. Int J Fluid Power 2006;7(3):29–39. [13] Franzoni F, Milani M, Montorsi L. A CFD multidimensional approach to hydraulic component design. SAE Trans – J Comm Veh 2007;116:246–59. [14] Iaccarino G, Verzicco R. Immersed boundary technique for turbulent ﬂow simulations. Appl Mech Rev 2003;56(3):331–47. [15] Merrit HE. Hydraulic control systems. New York: John Wiley & sons; 1967. [16] Mittal R, Iaccarino G. Immersed boundary methods. Ann Rev Fluid Mech 2005;37:239–61. [17] Mohd-Yusof J. Combined immersed-boundary/B-spline methods for simulations of ﬂow in complex geometries. Annual Research Briefs, Center for Turbulence Research, University of Stanford; 1997. p. 317–28. [18] Posa A, Lippolis A, Verzicco R, Balaras E. Large-eddy simulations in mixed-ﬂow pumps using an immersed-boundary method. Comput Fluids 2011 ;47(1):33–43. [19] Swarztrauber PN. A direct method for the discrete solution of separable elliptic equations. SIAM J Numer Anal 1974;11(6):1136–50. [20] Vanella M, Rabenold P, Balaras E. A direct-forcing embedded-boundary method with adaptive mesh reﬁnement for ﬂuid-structure interaction problems. J Comput Phys 2010;229(18):6427–49. [21] Verzicco R, Orlandi P. A ﬁnite-difference scheme for three-dimensional incompressible ﬂows in cylindrical coordinates. J Comput Phys 1996 ;123(2):402–14. [22] Verzicco R, Mohd-Yusof J, Orlandi P, Haworth D. Large eddy simulation in complex geometric conﬁgurations using boundary body forces. AIAA J 2000;38(3):427–33. [23] Yang R. Hydraulic oil ﬂow wall shear effect on valve actuator ﬂow forces. In: Proceedings of the 2nd international conference on computational methods in ﬂuid power technology. Aalborg (Denmark), 2–3 August 2006.