Prediction of acoustic response on tube-rows with bias-flow using linearized
Navier-Stokes equations in frequency domain
Wei NA(1) , Susann BOIJ(2) , Aswathy SURENDRAN(3) , Dong YANG(4) , Aimee MORGANS(5)(1,2) KTH Royal Institute of Technology, Department of Aeronautical and Vehicle Engineering, Stockholm, Sweden, wein@kth.se(1,2) KTH Royal Institute of Technology, The Marcus Wallenberg Laboratory for Sound and Vibration Research (MWL), Sweden(3,4,5) Imperial College London, Department of Mechanical Engineering, UK, a.surendran@imperial.ac.uk

Abstract
Heat exchangers are widely used in industries where heat needs to be transferred from one fluid to another fluid.
For example, there are plate-fin heat exchangers in gas turbine engines, shell and tube heat exchangers in oil
refineries and tube bundle heat exchangers in domestic boilers. Among the different types of exchangers, tube
bundle exchangers are the most commonly used heat exchange equipments. In this paper, heat exchanger withthe structure of tube rows with bias flow and possible acoustic interaction at the tube row is studied for the
acoustic reflection and transmission. The numerical methodology solving the linearized Navier-Stokes equationin the frequency domain is used. It has the advantage of taking into account the flow effects, viscous losses aswell as thermal losses in the acoustic propagation. The simplified geometry for the heat exchanger investigatedin this paper is a two-dimensional rectangular duct with two half cylinders with a bias flow going through thegap between the cylinders. In the current study, the acoustic response is predicted numerically with the cold
flow only and compared to the experimental data, as a preparation for the next study with the hot flow.
Keywords: Heat Exchangers, Bias-flow, LNSE, Frequency domain

1

INTRODUCTION

Tube-rows exist widely in engineering applications, for example, in tube bundle heat exchangers, in domestic
boilers, in air coolers. Multiple physical phenomenon can happen when flow passing through the tube rows.
As addressed in [1], from the acoustic point of view, acoustic waves can scatter on the tube rows, leading to
the sound transmission and reflection; flow separation due to the cross-flow, leading to the generation of thevorticity waves, etc.

Figure 1. Schematic of three-dimensional tube-rows
with bias-flow.

Figure 2. The simplified two-dimensional numerical geometry of the tube-row heat exchangers investigated in the paper, where R = 10 mm,
H p =25 mm, Lu = 40R and Ld = 40R, respectively.

An array of cylinders in cross-flow can be modelled as a single cylinder in a hard-walled flow-duct through2134

using the image sources method [1]. Therefore, the three-dimensional heat exchanger with tube-rows as Fig.
(1) shown is simplified to the two-dimensional rectangular duct with two half cylinders as Fig. (2) illustrates. In
the paper, the acoustic response is investigated numerically with the geometry as Fig. (2) shown with bias-flow
at different Mach numbers.
When the flow is passing through the two half cylinders within the duct, a jet flow is formed. Recirculation
zones are formed as well by shear layers along the jet flow. Such shear layers are unstable to low frequency
perturbations, energy is transferred from the acoustic waves to the vorticity waves. It is of interest to investigate
the flow-acoustic interaction for the tube-row heat exchangers.

2

NUMERICAL METHODOLOGY

The linearized Navier-Stokes equations are derived from the full compressible Navier–Stokes equations. A frequency domain approach has been taken by prescribing harmonic time-dependence of the perturbed quantities.
In this way, any perturbed quantity can be represented as q0 (x, ω,t) = Re{q̂(x)eiωt }, where q̂ is a complex quantity and ω is the angular frequency. The frequency domain linearized Navier-Stokes equations used in the paper
are presented as:iω ρ̂ + ∇ · (ρ0 ûu + ρ̂uu0 ) = 0iωρ0 ûu + (ûu · ∇)uu0 + (uu0 · ∇)ûu + ρ̂(uu0 · ∇)uu0 = ∇ · σ(1)T̂ρ0C p (iω T̂ + ûu · ∇T0 + u0 · ∇T̂ ) + ρ̂C p (uu0 · ∇T0 ) − (iω p̂ + ûu · ∇p0 + u0 · ∇ p̂) − (uu0 · ∇p0 ) = ∇ · (k∇T̂ ) + Φ
T0
with

2
σ = − p̂pI + µ(∇ûu + (∇ûu)T ) + (µB − µ)(∇ · ûu)I3andρ̂ = ρ0 (

p̂
T̂
− )
p0 T0

where a hat ˆ indicates a perturbed quantity, a subscript zero indicates mean flow quantities, ρ is the density,
p is the pressure, u is the velocity vector, T is the temperature, ω is the angular frequency, µ is the dynamic
viscosity. µB is the bulk viscosity, C p is the heat capacities for constant pressure, k is the thermal conductivity,
I is the unit matrix and Φ is the viscous dissipation of mechanical energy due to the shear viscosity and thefluid relaxation losses, which is a function of velocity, dynamic viscosity, and bulk viscosity, defined as:2
Φ = [µ(∇ûu + (∇ûu)T ) + (µB − µ)(∇ · ûu)I] · ∇ûu
3
Boundary conditions are needed for each dependent variable. Different boundary conditions are set up for the
top and the bottom duct walls as well as the walls of the heat exchangers. The following boundary conditions
are used in simulations.
At the top and bottom duct walls, a rigid, symmetrical, adiabatic wall boundary condition is imposed as:n̂n · ûu = 0, n̂n · ∇ρ̂ = 0 and − n̂n · (k∇T̂ ) = 0(2)

In the vicinity of the heat exchangers, at the wall of two half cylinders, rigid, no-slip, isothermal wall boundary
conditions are applied as:ûu = 0, n̂n · ∇ρ̂ = 0 and T̂ = 0(3)
where n̂n is the unit vector normal to the duct wall. It should be noted that as the velocity components are
restricted to zero on all edges.

2135

3

BACKGROUND MEANFLOW

Linearization of the Navier–Stokes equations is performed by splitting the quantities of the total field into
the mean part and the perturbed part, where the mean flow field is supposed to be obtained first for thefollowing acoustic simulations. In this paper, the steady state mean flow is simulated with a compressible
Reynolds-Averaged Navier–Stokes (RANS) solver, with the SST turbulence model. Figure (3) presents thesteady state mean flow field of a simplified two-dimensional heat exchanger geometry when inlet flow speed isuin = 11.5 m/s, yielding a Mach number of Ma = 0.033. In the computational dynamic fluid (CFD) simulations,the velocity is prescribed at the inlet of the duct; the top and bottom of duct wall, the symmetric boundary
condition is applied; on the two half cylinders, the non-slip wall boundary condition is applied together withthe wall functions at y+ ≈ 1; at the outlet of the duct, an ambient pressure boundary condition is applied.

Figure 3. Mean flow field with the inlet flow speed uin = 11.5 m/s, yielding Ma = 0.033: (a) shows the flow
velocity component u0x and streamlines of the steady-state meanflow in [m/s]; (b) shows the flow velocity
component u0y in [m/s]; (c) shows the pressure p0 in [Pa] and (d) shows the turbulent eddy viscosity µT in[Pa·s].
It is observed in Fig. (3) that a jet with a maximum flow speed of about 60 m/s, i.e. Ma ≈ 0.175 is formed
downstream of the two half cylinders due to the geometrical contraction. The jet has a length of about 2–4
duct diameters, expanding into a channel flow, a recirculation zone thus exhibits. In the recirculation zone, theflow reattaches at about 3-5 duct diameters downstream of the two half cylinders. Further downstream, the flow
is again the unperturbed channel flow. The sound–vortex interaction is assumed to happen at the shear layer ofthe jet flow.

2136

Figure 4. Schematic of the CFD mesh and the acoustic mesh in the vicinity of the half cylinder: (a) CFD
mesh; (b) acoustic mesh.

4

MAPPING BETWEEN CFD AND ACOUSTICS MESH

In most numerical cases, the CFD mesh and the acoustic mesh are two independent meshes. The CFD meshneeds to be generated with respect of the characteristics of the flow field, while, the acoustic mesh can begenerated with the rule of thumb that, 7 numerical cells per wave length are enough. In addition, the boundary
layer mesh in the CFD mesh and in the acoustic mesh does not need the same resolution. For example, in the
CFD simulations, in the near wall region, the flow boundary layer mesh is generated under the considerationof the dimensionless wall distance y+ with 7-8 cells, however, the acoustic boundary
layer mesh is generatedp
within 4-5 cells, the thickness is frequency dependent and calculated by δA = (2µ)/(ωρ), where ρ is densityof the fluid, µ is the dynamic viscosity and ω is the angular frequency.
When the CFD and the acoustic models are not solved on the same computational mesh, numerical noise canhappen as the terms containing gradients of the background mean flow variables can become very noisy. Therefore, careful mapping of the CFD solution from the CFD mesh onto the acoustics mesh should be done in order
not to introduce non-physical numerical noise into the acoustics solution. In the heat exchanger simulations, twodifferent mesh are used as Fig. (4) shows. The acoustic mesh used for the solutions of the perturbed quantities
is triangular mesh with about 54 000 elements.
A mapping technique with mapping equations that smooths the background flow variables is used. In the
mapping, the mean background flow pressure p0 , velocity field u0 and turbulent viscosity µT variables are
mapped onto the corresponding variables on the acoustics mesh p0(aco) , u0(aco) , and µT (aco) . The mapping andsmoothing is achieved by solving the Eq. (4):p0(aco) − p0 = δ h2 ∇ · ∇p0(aco)u0(aco) − u0 = δ h2 ∇ · ∇u0(aco)(4)

2µT (aco) − µT = δ h ∇ · ∇µT (aco)where the term on the right hand side adds smoothing using isotropic diffusion. The amount of diffusion iscontrolled by the parameter δ (a constant that can be tuned) and the mesh size squared, h2 . The term in somesense corresponds to source term stabilization as known from CFD.
The comparison of the velocity profile between the CFD mesh and the acoustic mesh with mapping techniqueis plotted in Fig. (5). The mapping technique generally works well, there are only discrepancies at the sharpedge of the velocity profile, as Fig. (5) shown.

2137

CFD mesh
acoustic mesh

0.01

0.005

y-axisy-axis

0.005

0

0

-0.005-0.005

-0.01

-0.01

-10

CFD mesh
acoustic mesh

0.01

0

10

20

30

40

50

-10

0

10

x-Velocity Component(a) at x = 0.05 m, located in the jet region.

20

30

40

50

60

x-Velocity Component(b) at x = 0.15 m, located in the donwstream channel flow.

Figure 5. Comparison of the velocity profile between the CFD mesh and the acoustic mesh.

5

ACOUSTIC SIMULATIONSAs the mean flow field is calculated already, the perturbed field including both acoustic as well as vorticity
perturbations is solved through the Eqs. (1). The simulations are carried out in COMSOL MultiphysicsV5.3,
which is a commercial finite element method (FEM) solver. The Mach number and the frequencies are chosencorresponding to those measured in the experiments of [1]. The frequency of the incident wave is 150 Hz,
250 Hz and 350 Hz, respectively. The seven mean flow velocities are: uin = 5.0, 7.5, 9.5, 11.5, 13.5, 14.5 and 15.5
m/s, corresponds to the Mach number: Ma = 0.015, 0.021, 0.027, 0.033, 0.039, 0.042 and 0.045.
As shown in Fig. (6) and Fig. (7), at the cylinders, acoustic waves are scattered, part of the acoustic waves
are transmitted downstream and part of the acoustic waves are reflected upstream. Besides that, some of theacoustic energy is transferred to the vorticity waves. The vorticities are convected by the mean flow, thosebecome even stronger in the shear layer of the flow.
In order to obtain the acoustic scattering properties for the heat exchangers with tube-row geometries, and tocompare with experimental results, the perturbed field is decomposed into upstream and downstream propagating
acoustic waves. The acoustic properties are then formulated as a scattering matrix by using the Two Microphone
technique (TMM), more details about the TMM can be found in [2]. The magnitudes of the enthalpy scatteringmatrix calculated from simulations as well as experimental results are shown in Fig. (8). In the experiments,
usually the acoustic pressure is measured, however, the enthalpy scattering matrix is present in this paper byusing the convention below [1]:



 

 P

1−M2
1
T1→2 1+MR
2
T1→2 RP2
1+M
1+M
2


 2 
=
P
1−M2
1+M1
RP1
T2→1
R
T
1
2→1
1−M1
1−M1
|
{z
}
|
{z
}
Measurements(5)

Scattering elements plotted in Fig. (8)M1 refers to the Mach number at the inlet of the duct, while, M2 refers to the Mach number further downstreamof the heat exchanger.
As can be seen from Fig. (8) the results from the simulations are in good agreement with the experimental
results from [1]. The numerical results predict the scattering of enthalpy precisely when the incident wave is
imposed from the downstream side, corresponding to the scattering matrix elements |T2−1 | and |R2 |. There are2138

Figure 6. Perturbed field at frequency f = 2000 Hz for Ma = 0.033 (a) particle velocity perturbations in [m/s](b) density perturbations in [kg/m3 ] (c) temperature perturbations in [K] (d) pressure perturbations in [Pa].

Figure 7. Perturbed field at frequency f = 350 Hz for Ma = 0.033 (a) particle velocity perturbations in [m/s](b) density perturbations in [kg/m3 ] (c) temperature perturbations in [K] (d) pressure perturbations in [Pa].

2139more discrepancies for the scattering matrix element |R1 | compared with the experimental results, and this ispossibly due to inaccuracies in the simulations of the meanflow, i.e., the thickness of the shear layer or thereattachment point of the flow. It has been noticed that |R1 | seems to be the scattering matrix element thatis most sensitive to mean flow errors [3]. In the acoustic numerical simulations, it is always a difficulty toreproduce the same flow profile as the measurements.

Figure 8. Magnitudes of the enthalpy scattering matrix, plotted as a function of Ma number at frequencies
f = 150, 250, 350 Hz, respectively. Dashed lines with markers: simulations, solid lines: experimental resultsfrom [1].

6

CONCLUSIONSIn this paper, a numerical study of the acoustic scattering from tube-row heat exchangers with bias-flow atseven different Mach numbers is carried out. The simulations are performed by means of solving the linearized
Navier-Stokes equations (LNSE) in frequency domain, taking into account the hydrodynamic effect as well as
the viscous effect. The numerical methodology has been improved through solving the energy equation, insteadof using the isentropic relation for the acoustic field compared with the previous work at KTH [4], [5], [6].
The numerical results presented in this paper clearly show that the sound-vortex interaction is resolved by thecurrent numerical methodology, the sound vortex is generated along the shear layer of the jet flow and convectedby the mean flow downstream of the heat exchangers, the energy is thus transferred from the hydrodynamic
mode to the vorticity mode. In addition, the acoustic scattering matrix in terms of enthalpy is compared tothe experimental results when bias-flow is present, and good agreement is found. It indicates that the flow2140

effects with the Mach number up to Ma = 0.045 are not negligible on the acoustic response in tube-row heatexchangers. For future work, the acoustic response with incident waves at higher frequencies will be examined
and the heat transfer will be included in the heat exchanger simulations.

ACKNOWLEDGEMENTS
This research was financially supported by the Marcus Wallenberg Laboratory for Sound and Vibration Research(MWL) at KTH, Sweden.

REFERENCES[1] Aswathy Surendran, Maria A Heckl, Luck Peerlings, Susann Boij, Hans Bodén, and Avraham Hirschberg.
Aeroacoustic response of an array of tubes with and without bias-flow. Journal of Sound and Vibration,
434:1–16, 2018.
[2] Hans Bodén and Mats Åbom. Influence of errors on the two-microphone method for measuring acoustic
properties in ducts. The Journal of the Acoustical Society of America, 79:541–549, 1986.
[3] Emma Alenius. CFD of duct acoustics for turbocharger applications. 2010.
[4] W. Na, G. Efraimsson, and S.Boij. Simulations of the scattering of sound waves at a sudden area expansion
in a 3D duct. In ICSV 21: The 21st International Congress on Sound and Vibration, 2014.
[5] Wei Na. Frequency domain linearized navier-stokes equations methodology for aero-acoustic and thermoacoustic simulations. 2015.
[6] A. Kierkegaard, S. Allam, G. Efraimsson, and M. Åbom. Simulations of whistling and the whistling potentiality of an in-duct orifice with linear aeroacoustics. Journal of Sound and Vibration, 331:1084–1096,2012.

2141