Abstract
Increasing demands on the capabilities of engine thermo-fluid dynamic simulation and the ability to accurately predict both performance and acoustics have led to the development of several approaches, ranging from fully 3D to simplified 1D models. The quasi-3D approach is proposed as a compromise between the time-demanding 3D CFD analysis and the fast 1D approach; it allows to model the acoustics of intake and exhaust system components, used in internal combustion engines, resorting to a 3D network of 0D cells. Due to its 3D nature, the model predicts high-order modes, improving the accuracy at high frequencies with respect to conventional plane-wave approaches. The conservation equations of mass and energy are solved at cell centers, whereas the momentum equation is applied to cell connections including specific source term to account for the of sound-absorbing materials and perforated elements. The quasi-3D approach has been validated by comparing the predicted transmission loss to measured data for a number of standard configurations typical of internal combustion engine exhaust systems: a reverse-flow chamber and series chambers with perforates and resistive material.
1. Introduction
Internal combustion engines are the source of mechanical, combustion, and gas dynamic noise. In particular, gas dynamic or pulse noise is related to the unsteady flows in the intake and exhaust systems, induced by the cylinder gas exchange process [1]. The level and quality of noise radiated from the open ends can be controlled by different arrangements of pipe systems and silencers to achieve the required vehicle exterior and interior sound characteristics by attenuating or enhancing certain spectral components. Today numerical simulation codes are very useful during the design and optimization process of both intake and exhaust manifolds and mufflers to quickly define the best geometries which can be finally refined experimentally. The attenuation features of simple and complex acoustic filters can be predicted by 1D and 2D-3D fluid dynamic/acoustic simulation codes with different levels of complexity, both in the time and frequency domain. On one hand, 1D linear/nonlinear codes are nowadays widely used to calculate the acoustic performances of mufflers [2, 3]. Certainly 1D linear acoustic codes are mainly applied for this purpose [4–6] due to their simplicity, accuracy, and very short computational times. However, as the flow speed becomes not negligible and the sound pressure level of the signal reaches values of 160 dB, as it can happen inside a muffler for internal combustion engines, the magnitude of the incident pressure perturbation makes the linear approximation not valid anymore [7, 8]. Therefore, when the dynamics of the system is nonlinear, time domain simulation tools become attractive. For instance, in the case of perforated tube silencers, the fluid dynamics of the holes exhibit nonlinear features when exposed to sound pressure levels typical of internal combustion engines. For this reason, time-domain simulation tools are preferably applied to better reproduce complex pulsating flows associated with mean velocity, since they directly account for these features in the fundamental equations they are based on [7, 9]. Therefore, in the last decade 1D nonlinear, time domain models have become preferential tools for the simulation of the high amplitude wave motion in the silencer ducts with mean flow, in order to predict the acoustic response of the system and the noise radiated when it is coupled to the engine source itself. Many researchers have pointed out the successful development and validation of advanced silencer modeling by a time-domain, nonlinear, one-dimensional approach, including complex elements like perforates [10], flow reversals [11], sound absorptive materials [12], with both engine or white noise/impulse excitation [13, 14]. The available approaches are based mainly on the segmentation of the perforated pipe and the related cavity [15, 16], each one characterized by having specific advantages and constrains, such as accuracy, predictivity, and required computational time. With this type of approach, it is some times mandatory to know how different arrangements act acoustically, in order to model them correctly using the appropriate modeling elements. This means that they require an high level of knowledge by the user, in order to build the correct 1D equivalent schematic of the device. On the other hand, linear and nonlinear 1D approaches are not adequate to capture the nonplanar wave motion in complex geometry systems and the related higher order modes, due to the intrinsic limitation of the plane-wave assumption. To this purpose, a simple or complex shape device (air-box or muffler), characterized by a significant multidimensional wave motion in a certain frequency range, can be successfully simulated resorting to a multiD linear/nonlinear model. MultiD linear codes, based on the finite element method (FEM) or boundary element method (BEM) [17–19], are commonly used for acoustic simulations. Nonlinear 3D approaches have also been exploited by some authors showing good potentialities for the prediction of silencer attenuation characteristics with and without mean flow [20–24]. However, the adoption of 3D nonlinear models in not so widespread as 3D linear models, one reason is certainly related to the significant computational effort, typically two-three orders of magnitude greater for CFD models. In particular, in cases when perforated baffles and pipes are present, the combination of a fine mesh resolution, required to capture correctly the porosity, and the reduction of time step, due to stability issues, may lead to a significant increase of computation runtime. Conversely, the higher computational cost of time-domain nonlinear multiD codes can be balanced by the wider applicability of simulation and validity of the results, ranging from the muffler Transmission Loss (TL), including high amplitude wave effects, to acoustic and fluid dynamic performances of silencing systems under engine pressure pulsations. Nowadays nonlinear 1D and CFD fluid dynamic codes are becoming rather common for acoustic/fluid dynamic applications [13, 22, 25, 26], since they allow the prediction of silencer performances and radiated pulse noise when the duct system is connected to the engine source itself, submitted to high amplitude pressure waves and significant mean-flow. These conditions are directly taken into account in the fundamental conservation equations of nonlinear models. The counterpart is represented by linear acoustic models (1D or multiD), applicable in the case of small (acoustic) pressure perturbations, which frequently is not acceptable in I.C. engines pipes. Linear acoustic codes are characterized by reduced computational times, generally 2 or 3 orders of magnitudes with respect to nonlinear CFD codes. Furthermore, in recent years the great potential of hybrid 1D-multiD fluid dynamic (nonlinear) models has been recognized. Integrated 1D and multiD computational domains can be built to carry out the simulation of the whole engine intake and/or exhaust duct systems including the silencers, in order to dedicate the 3D approach to model the complex shapes and reserve the 1D approach to simulate simple pipe systems [26, 27]. As an alternative to the fully multiD CFD simulation, a nonlinear quasi-3D approach (also referenced as “3Dcell” method) has been investigated recently [25, 28, 29], in order to develop an intermediate tool between 1D and 3D time-domain nonlinear models, representing a good compromise in terms of computational effort and prediction accuracy of nonplanar wave motion. The quasi-3D method relies on a network of 0-D elements, such as volumes, connected to the others by means of specific connecting elements, which allow to define the characteristic lengths and shape of the device. Due to its 3D nature, this method is able to predict higher order modes, improving the accuracy at high frequencies, with respect to conventional 1D plane-wave approaches. The quasi-3D approach has been extended to the treatment of perforated elements and of sound absorbing materials, modeling their impact on the propagation of acoustic waves with specific source term included in the governing equations. This work describes the development of a quasi-3D method focusing on the modeling of perforates and sound absorbing material as well as its validation against experimental data for a number of standard configurations. Both the simulation and measured data for the described components have been carried out without mean flow.
2. Conservation Equations
The numerical method is based on the formulation of the conservation equations of mass, momentum, and energy for unsteady flows. Since the viscosity of the gas both in the intake and exhaust system is very low, it can be neglected without introducing an excessive approximation, leading to the Euler formulation [30]: This system of equations is then closed with the perfect gas equation of state.
3. Numerical Method
The method adopted for the solution of the Euler set of equations operates on a staggered computational grid (Figure 1), in which the intensive properties are evaluated in the cell, while velocity and momentum are defined in the ports.
The conservation equations (1), (2), (3) are integrated according to a pseudostaggered arrangement [31, 32]. According to this definition, the fluxes at the port elements are used for the solution of the scalar quantities defined over the cell element. The solution procedure adopted is an explicit time marching staggered leapfrog method in which the momentum equation is used for determining the port flux. This assumption leads to the definition of three different time levels (, , and ) and the integration of the three conservation equations (1), (2), (3) is suitably distributed over them in order to achieve a low diffusive method with second order accuracy both in space and time. In particular, considering the mass equation, the integration is performed over the cell volume and the time-step : By means of the Gauss theorem, the transport term is reduced to a surface integral, giving the following: which, on the basis of the finite volume discretization, becomes a summation of the mass flows through the ports (). Taking into account the time marching stencil, the mass conservation becomes In a similar way, the discretized energy equation has the following expression: where is the enthalpy evaluated in the port element. The integration of the momentum equation is performed referring to a control volume defined by the port surface and the distance between the two neighboring cell centers (Figure 1). Since the velocity vector in the port is directed normally to the surface, according to the particular grid adopted, the momentum equation reduces to a 1D problem arbitrarily oriented in space. By naming the coordinate along the normal to the port, the integration of the momentum equation over the control volume leads to the following: where the subscripts and indicate, respectively, the quantities in the left and right cells, whereas is a generic source term which may assume different expressions, depending on the type of the port that has to be modeled. The solution of the momentum equation allows to determine the port mass flow This quantity is then used in the mass and energy conservation equations to march ahead in time. In particular, referring to the solution of the energy equation, the total enthalpy flow is expressed as the product between the mass flow and the enthalpy value at the port location: In case the achievement of second order accuracy in the description of contact surfaces is addressed, the value of can be estimated resorting to linear interpolation between the enthalpy in the neighboring cells. Otherwise, when the method stability is mandatory and the low accuracy resolution of contact surfaces does not affect the parameter of interest, the adoption of an upwind approach, by assigning the enthalpy of the left or right cell depending on the flow direction, may be exploited.
According to the definition of pseudostaggered arrangement, to allow the construction of the port momentum equation (9), the velocity must be defined also in the centroid of the cell element. The cell momentum is a vectorial quantity and its direction is determined by the momentums of all the ports linked to the cell weighted by the their distance from the cell center: The cell velocity is then determined from the cell momentum as follows: In Figure 2 it is shown how the staggered leapfrog method is applied to a staggered grid arrangement. The computational control volumes for the momentum equation are represented as small squares, while the cell control volumes are represented by circles. They are centered at different locations both in the scale of space and time.
In this way the fluxes calculated in the ports at a generic time step can be used to update the solution in the cells from time step to time step . The adoption of time-centered differences for each conservation equation, along with space-centered differences, allows the achievement of second order accuracy both in space and time.
4. Diffusion Term on Momentum Equation
The method discussed so far does not satisfy the stability requirement, especially in cases where pressure gradients are significant. For instance, simulations with mean flow or with high amplitude pressure perturbations traveling along the calculation domain (pulsating flow generated by an engine) will be affected by nonphysical oscillations. In order to apply the method to the solution of unsteady flows in the engine intake and exhaust system, a stabilization technique has been applied to reduce this spurious behavior. Due to the specific structure of the solution procedure, in which the solution of the momentum equation on the port is used to determine the fluxes for mass and energy equations on the cell, it is possible to stabilize the method acting only on the momentum equation. Therefore, the idea is to add a diffusion source term to the momentum equation (hence the name DTM) which will limit the computed mass flux at the port. In the DTM approach the momentum equation (2) is modified by adding a source term as follows: The tensor is a function of the local velocity fields and of the spatial and temporal discretization, and is defined as: where is a scalar quantity and represents the diffusion coefficient.
If the velocity field is smooth the contribution of the diffusion term is negligible. However, when large gradients in the solution occur, this term increases in order to damp the numerical spurious oscillations. Projecting (12) along the port direction (), it reduces to a 1D problem arbitrarily oriented in space. In this case only the component of velocity in the port direction () is considered, therefore the term assumes the following expression: The discretization of (12) along the port direction will lead to the following simplified 1D equation: This definition requires the evaluation in the cells of the diffusivity coefficient and of the terms and (generically referred to as hereafter). The coefficient is evaluated on the basis of the local velocity, the time-step and the mesh size: The term is a vector defined in the cell as the product between the diffusion coefficient and the gradient of the mass flow in each direction: The scalar product between and the port direction gives the term used in (15). The gradient of mass flow in the cell is calculated from the mass flow computed in the attached ports. In particular, it has been computed considering the projections of the port mass flows along each direction , , : The cell dimensions (, , ) have been determined starting from the lengths of the attached ports, since, according to the staggered arrangement, these are properties of the connecting elements and are not initially defined on the cells.
5. Perforates Modeling
Perforated elements have been modeled by extending the properties of the connector element. The perforated connector is characterized by the number of holes , hole cross section , hole length , and by the total connector surface . The approach developed for the perforated elements solves the momentum equation for only one single hole, in order to calculate its mass flow, and calculates the connector mass flux by knowing the number of holes afterwards. As shown in Figure 3, the momentum equation for one single hole is solved over a control volume defined by the hole cross section and the length . This length is determined on the basis of end corrections as a function of the hole length and hole diameter , according to the following relation: Usually is considered to be 0.85 in the case of one orifice in an infinite plate. To account for the mutual interference of neighboring holes the value of porosity can be considered [33, 34]:
This length correction is necessary to take into account the inflow-outflow phenomena and the fact that the flow field is extended outside the hole. In this particular case the value of used for perforated elements is 0.8 [35]. The momentum equation (9) applied to the hole control volume will become as follows: where , , are, respectively, the density, the pressure and the normal velocity at the control volume boundaries. In particular and have been assumed equal to the same quantities in the neighboring cell, while the velocity has been interpolated between the connector and cell centers as follows: The term in the momentum conservation equation (21) accounts for the gas-wall friction and is determined as follows [1]: where the friction factor is defined as [1]:
6. Absorptive Material Modeling
The 3D cell model has been extended to account for the effect of sound absorbing material. This material usually fills up cavities and gaps created inside the muffler, as a consequence the presence of porous material is considered as a cell property. The flow interacts with the porous media transferring to the fibers its momentum. The acoustic energy associated to pressure waves is then dissipated as heat. To take into account this dissipation, the momentum equation has been modified by including an additional source term, which reproduces the effects of the flow resistance operated by the absorptive material. The momentum equation, which in the staggered stencil is solved on the connector, has been modified by adding two terms accounting for the resistivity properties of the two neighboring cells: where expresses the material resistivity. This property is the most important characteristic of a porous material and can be defined by the following equation: The material resistivity can be calculated by the following correlation [36]: where is the dynamic viscosity of the gas, is the mean square average fiber radius, the ratio between the material packing density and the material density, and is the gas velocity. The , , , coefficients depend on the type of sound absorptive material used. For example, for rock wool the following values are proposed: , , , [37].
7. Results
The validation of the developed submodels for the 3Dcell approach has been carried out by comparing the predicted transmission loss to the measured values of different configurations: a reverse-flow chamber with pipe extension, an expansion chamber with perforated elements, and an expansion chamber with extended outlet duct and absorptive material.
7.1. Reverse-Flow Chamber
The first validation case is a reverse-flow chamber with extended inlet and outlet [26]. This chamber has a length of 494 mm and a diameter of 197 mm. The inlet and the outlet pipes are extended into the chamber, respectively, by 257 mm and 17 mm and their diameter is 50 mm. The 3Dcell model was built using a network of 5782 cells, resulting in an average cell size of 17 mm in each direction. As shown in Figure 4, the adopted mesh reproduces the circular section of the chamber in such a way that the maximum transversal dimension of the model does not exceed the dimension of the chamber diameter (197 mm). This allows the model to capture the occurrence of eventual higher order modes [38, 39]. In order to correctly capture the length of the duct extensions, a non uniform mesh has been adopted in the longitudinal direction and the chamber has been divided into three parts with slightly different mesh spacing: 17 mm the fist one, 17.14 mm the second one and 16.92 mm the third one (Figure 5). Simulations has been performed imposing a white noise perturbation at the inlet of the system and an anechoic termination at the outlet. Four oscillation periods, referring to the lowest frequency (20 Hz) used in the calculation, have been simulated, resulting in approximately 5 minutes of simulation runtime.
Figure 6 shows the comparison between the predicted transmission loss and the measured values for the reverse-flow chamber. Generally the agreement can be regarded as good, since the abatement peaks and pass bands are correctly captured over a wide range of frequency, without showing any frequency shift. In particular, the first resonant peak at 414 Hz corresponds to the quarter wave resonance of the reverse-flow length minus the inlet extension, while the second peak at 974 Hz is related to the inlet extension. The transparency frequencies (354 Hz, 720 Hz, 1062 Hz, 1460 Hz) are associated to the main longitudinal resonance. Higher order modes are excited at frequencies higher than 1000 Hz.
7.2. Chamber with Perforated Elements
The chamber with perforated elements has a length of 305 mm and a diameter of 114 mm. The perforated duct is characterized by a diameter of 50 mm, a thickness of 1.6 mm and a porosity of 4.6% [10]. The diameter of the holes is 3.2 mm and the corrected length adopted in the simulation resulting from the formula proposed by Torregrosa is 4.8 mm [35].
As shown in Figures 7 and 8, the cluster of calculation cells is a 7 by 7 by 20 structure and the circular section is approximated by a square shape one. In particular, the reconstruction of the chamber section by means of a square shape is due to the fact that the higher order modes are triggered at a frequency out of the range of interest. For this kind of geometry (cylindrical and symmetric), the cut off frequency of cylindrical higher order modes is given by the following formula [38]: Hence, resorting to a square reconstruction the mesh size can be adjusted (coarsened in this case) in order to achieve the best compromise between accuracy and computation time. The number of periods adopted for the simulation is once again 4, resulting in approximately 0.5 minute of simulation runtime. The comparison between predicted TL and measured data is shown in Figure 9. The 3Dcell simulation code correctly predicts the transparency frequencies associated to the main longitudinal resonance (600 Hz) and the resonances peaks associated to the perforated element.
7.3. Chamber with Absorbing Material
The validation of the sound absorptive material was carried out resorting to measured data provided by the ACC in Graz [40]. The silencer considered is an expansion chamber with an extended outlet pipe, whose characteristic dimensions are summarized in Figure 11. The space between the extended pipe and the muffler canning has been filled with sound absorbing material, without recurring to the use of a perforated plate. The 3Dcell cell matrix is built using a network of 1522 cells with an average mesh spacing of 2 cm.
As shown in Figures 10 and 12, the mesh adopted reconstructs the circular section of the chamber in order to better capture eventual higher order modes. A nonuniform mesh in longitudinal direction has been adopted to correctly locate the outlet section connector. The same number of periods adopted in the previous cases have been simulated, resulting in approximately 1 minute of simulation runtime. Comparisons between predicted TL and measured data are reported in Figures 13 and 14, respectively, associated to the silencer configurations without and with the dissipative material. In the first case it is shown a good agreement, since the 3Dcell correctly predicts the first transparency frequency at 556 Hz, due to the main longitudinal resonance, and the peak at 940 Hz due to the quarter wave resonance associated to the outlet extension. The configuration with the dissipative material, characterized by a resistivity of 8000 N s/m4, has been correctly predicted as well, accurately locating the resonance and transparency frequencies. Furthermore, a sensitivity test has been performed on the resistivity of the material. In particular, it has been increased by one and two orders of magnitude and the results have been compared, showing that the interface between the sound absorbing material and the empty region of the series chamber behaves as a partially reflecting surface. Figure 15 shows the effect of varying the resistivity of the sound absorptive material. It can be noted that by increasing the resistivity up to a value of 300000 N s/m4, the surface behaves more like a rigid wall and the TL resembles the shape of a simple expansion chamber whose main resonance frequency is given by the length of the empty chamber (shorter than the baseline).
8. Conclusions
In this paper a general quasi-3D method for the acoustic modeling of silencers has been presented. The aim of this work was to show the potentialities of a quasi-3D approach for the prediction of simple and complex muffler performances with sound absorbing material and perforated elements. This model has been validated on three different muffler configurations, namely a reverse-flow chamber, an expansion chamber with one perforated pipe and an expansion chamber with sound absorptive material. The comparison between the calculated TL and the measured data has pointed out the capabilities of the model in predicting of the acoustic phenomena occurring in these mufflers. In particular, the model is able to predict the frequencies of the quarter wave resonances associated to inlet and outlet extensions, the peaks associated to perforated elements, and the dissipative effect produced by the absorptive material. The presented model is based on the pure geometrical reconstruction the mufflers and does not require the creation of an acoustically equivalent model. Moreover, due to the coarse level of mesh adopted, the computational burden of the simulations is very low if compared to other 3D CFD methods, allowing the adoption of DOE for fast muffler optimization.
Acknowledgment
The research of the Acoustic Competence Center is carried out with the support and funding of AVL List GmbH in Graz, Austria.