A low–storage Runge–Kutta OpenFOAM solver for compressible low–Mach number ﬂows: aeroacoustic and thermo–ﬂuid dynamic applications

. A solver for compressible Navier–Stokes equations is presented in this paper. Low-storage Runge-Kutta schemes were adopted for time integration; on the other hand the ﬁnite volume approach available within OpenFOAM library has been adopted for space discretization. Kurganov-Noelle-Petrova approach was used for convective terms, while central schemes for di ↵ usive ones. The aforementioned techniques were selected and tested in order to allow the possibility of solving a broad range of physical phenomena with particular emphasis to aeroacoustic and thermo-ﬂuid dynamic problems. Indeed, that standard OpenFOAM solution techniques produce an unacceptable dissipation for acoustic phenomena computations. Non–reﬂective boundary treatment was also considered to avoid spurious numerical reﬂections. The reliability and the robustness of the solver is proved by computing several benchmarks. Lastly, the impact of the thermal boundary conditions on the sound propagation was analyzed.


Introduction
In recent years computational fluid-dynamics (CFD) was widely employed in several scientific and industrial research fields: the constant improvement and availability of computational resources allowed to face increasingly complex real-life problems. In this scenario, it is worth noting that compressible Navier-Stokes equations (NSE) can describe a broad range of scientific interest physical phenomena. The aim of this work is to apply a finite volume solver for compressible NSE, named caafoam, under development by the authors of this paper. caafoam is able to fully resolve the aeroacoustic waves propagation, heat transfer and their interaction with a direct numerical simulation (DNS) approach. DNS was selected in order to make public the possibility to performs this kind of computations to general CFD users since, in this moment, similar techniques can be used only by a limited number of academic developers. The direct simulation of the aeroacoustic sound is a challenging task because the sound pressure is usually much smaller than the ambient one, [1], thus numerical dissipation and dispersion have to be controlled to avoid sound production depletion. Moreover standard boundary conditions produce acoustic waves spurious reflections; to fix these issues low-dissipative discretization schemes and non-reflecting boundary treatment are strictly required, [2]. For this reason in literature we can ⇤ e-mail: v.dalessandro@univpm.it ⇤⇤ e-mail: m.falone@pm.univpm.it ⇤⇤⇤ e-mail: l.giammichele@pm.univpm.it ⇤⇤⇤⇤ e-mail: s.montelpare@unich.it find several papers about the adoption of high-order finite di↵erence (FD), [3], finite volume, [4], and, more recently, Discontinuous Galerkin (DG), [5] methods for aeroacoustics. In similar approaches Runge-Kutta (RK) methods were employed for time-integration. Di↵erently, in this work we adopt a density-based solver relying on the collocated finite volume method for space discretization available in OpenFOAM library, [6]. Low-storage high-order RK schemes are adopted for time-integration, while an artificial sponge-layer was used as non-reflective boundary approach. Our solver was already tested with several benchmarks to prove its reliability in aeroacoustic applications, [7]: in particular aeolian tones radiated from rigid bodies in an uniform velocity inlet flow were investigated. In this work the reliability of caafoam for forced convection heat transfer problems was also tested. Lastly, we investigate the impact of the thermal boundary conditions on the sound propagation; it will be shown that the wall temperature increment can reduces the lift pulsation and increase the drag generated the by Karman vortex street that is shed behind blu↵ bodies. Similar results were already found for the lift and drag modification on low Reynolds number operating airfoils [8][9][10]. In this context lift reduction is a very important element since it can lead to an aeroacoustic perturbation decay and in a modification of the sound directivity. This paper is organized as follows: in Sec. 2 we present our flow model; in Sec. 3 the key numerical ingredients of caafoam are descried. Sec. 4 is devoted to the numerical results; and Sec. 5 contains the conclusions.

Flow model
The flow model equations, based on the set of conservative variables = (⇢, ⇢u i , ⇢E) T read: where F c ( ) and F v ( , r ) are convective and viscous fluxes respectively given by: The viscous stress tensor is computed using the constitutive relation for Newtonian fluids, while for heat flux vector components the Fourier postulate is used: The pressure is computed assuming the ideal gas equation of state as thermodynamic model, while the fluid temperature is evaluated starting from the total internal energy according to: The nonphysical term, in the right-hand side of eq. 1 represent the expression of the artificial sponge layer; it produces a damping of the flow variables near the external boundaries to a known reference solution, re f . The term σ is a spatially varying coefficient defined as follows: where L sp is the layer width, d the minimum distance from the inflow/outflow boundaries, σ 0 is a constant value and n an integer parameter controlling the shape of the sponge profile. The optimal sponge layer design is a complex problem; larger sponges are better than small ones with the same strength. Indeed, they damp flow features more quietly, [11]. However, larger sponge layers require larger computational domains, thus they are more expensive. Mani, [12], faced a theoretical and numerical analysis of artificial sponge-layers. He provided several practical guidelines to CFD/CAA practitioners in order to avoid the sponge failure. In particular the non reflective boundary implementation is based on the following parameter: where -⌘ target defines the sponge strength and it is expressed in dB. This kind of sponge has been proved suffice to perform direct simulation of the aerodynamic sound. Moreover the sponge width must be set with the following constraint: where f is sound disturbance frequency and a 1 is the speed of sound, [12]. For all the computations presented in this paper we have fixed n = 2 in eq. 5, while the dimensionless parameter ⇣ L sp · f ⌘ /a 1 is 0.5 in order to limit the computational load.

Numerical approximation
A collocated cell-centered finite volume method, avaiable within OpenFOAM library, is adopted to solve the governing equations implemented in caafoam. The centralupwind scheme developed by Kurganov-Noelle-Petrova (KNP), [13], was used for convective terms approximation, while standard central schemes for di↵usive contributions are adopted. It is important to point out that KNP method is a second-order accurate scheme and it does not employ a Riemann solver technique for numerical fluxes approximation. Thanks to these features it allows to limit numerical dissipation without increase the computational load, [14]. For time integration a 2N storage explicit Runge-Kutta (ERK) scheme was implemented in order to exploit its good performance for large scale computations. In particular we have considered a five stage, fourth order accurate, ERK scheme available in the open literature, [15]. In our computational experience this scheme is really appealing since it allows to use a maximum Courant number equal to 1. Thus, it allows to exploit good parallel behaviour of the explicit approaches without a severe time-step size restriction.

Results
In this paper we have considered several flow/acoustic problems. A first aim of the presented cases is to ascertain the reliability of the numerical solver discussed in this paper for areoacoustics and forced convective heat transfer. Secondly we also present several results concerning the impact of thermal boundary conditions on sound generation and propagation. The flow past a square cylinders at di↵erent Reynolds numbers, Re = 90, 120, 150 as well as past two square cylinders in tandem configuration at Re = 150 were addressed. For the isolated square cylinder three wall boundary conditions for energy equation were considered: (@T/@n| w = 0, T w = 2T 1 , T w = 3T 1 ), while only the first two ones for the tandem configuration. In all the aforementioned cases the Mach number of the undisturbed flow is M 1 = 0.2, the ratio of specific heats γ = 1.4 and the Prandtl number, Pr, is fixed equal to 0.75. In the following we present the results in terms of the dimensionless parameters related to fluid dynamic, heat transfer and acoustic fields: the drag and lift coefficients, the Strouhal number, Sr, the Nusselt number, Nu, the fluctuating pressure, p 0 , and its root mean square, p rms . The drag and lift coefficients are defined by eq. 8: Being the cases treated in the paper markedly in unsteady regime the solutions are compared on the basis of the statistics of the force coefficients: lift coefficient root mean square C L,rms and time averaged drag coefficient, hC D i. The amplitude of oscillations of the force coefficients over their mean value are also used: The local Nusselt number Nu is calculated as: where T w is the wall temperature, and T 1 is the freestream temperature. The average Nusselt number, Nu avg , can be obtained by integrating the Nu along the body surface. The dimensionless fluctuating pressure is defined as: Polar plots containing the root mean square of p 0 , p rms , will be showed to investigate the nature of the sound in the far field. For this purpose, acoustic statistics were sampled over a dimensionless time of about u 1 T/D = 100.

Validation tests
For the sake of validation the flow past a square cylinder in a uniform velocity inlet was analysed. The sound radiated by the Karman vortex street that is shed behind the cylinder at Re = 150, under an adiabatic wall boundary condition, was already discussed in the open literature by Inoue et al., [16]. The acoustic field for this configuration can be observed in Fig. 1. On the other hand the same geometry was widely adopted in order to characterize the forced convection heat transfer for several Reynolds numbers, [17,18], so this problem is an appropriate benchmark.
In this study, a fully structured computational grid was used: the grid cells were clustered near the cylinder walls using a dimensionless first cell height, y c /D, equal to 5 · 10 −3 and the computational domain was extended up to 200 D from the center of the cylinder. The obtained computational mesh has a total cells number, N c , of 4.4 · 10 6 . For validation purposes, we have also considered the flow and the sound generation around two square cylinders in tandem configuration with L/D = 2, where L is distance between the cylinder centers and D their side length. The computational domain was generated in order to place the far field at 200 D from the cylinders center midpoint. The grid consists of about 4.2 · 10 6 quadrilateral orthogonal cells; near the cylinders walls a grid refinement was performed adopting y c /D = 10 −2 . In both the configurations, acoustic waves reflections at the far-field boundaries are removed using a the sponge layer strength of 40 dB, while its dimensionless thickness is 0.5 to limit the computational resource requirements. The computations presented in this work were run using 256 CPU-cores on MARCONI-A2 HPC system hosted by CINECA. MARCONI is a NeXtScale cluster consisting of 3600 compute nodes with 1 Intel Knights Landing (KNL) processor having 68 cores with 1.40 GHz. Each node is equipped with 96 GB of RAM and 16GB of Multi-channel Dynamic Random Access-Memory (MCDRAM).

Aeroacoustics
The polar plot containing the root mean square of p 0 is showed in Fig. 2. The data are collected over a circumference built around the square cylinder and having a radius r = 75 D, according with Inoue et al., [16]. It is clear that our approach ensures a good reconstruction of the acoustic far field whereas rhoCentralFoam, also equipped with the sponge layer, produces an un-physical jagged behavior of the p rms . Moreover caafoam shows a directivity of ✓ = 79 • that is in good agreement with the Inoue et al. prediction, i.e. ✓ = 80 • . Di↵erently rhoCentralFoam overestimates the directivity with a response of ✓ = 89 • . Regarding the aerodynamic parameters, our solver produced good results as showed in Tab. 1. From this point of view, also rhoCentralFoam is in accordance with the reference values.
As regards the cylinders in tandem configuration the dipolar nature of the sound field generated by the interaction of the flow and the two cylinders is clear from the directivity plot showed in Fig. 3. Note that the plot is referred to a circle having r/D equal to 80. caafoam results are in good agreement with literature reference data, [16], and a directivity equal to 68 • is highlighted. Lastly we want to remark that, also for this test case, the non-reflective version of rhoCentralFoam is inappropriate for the far field sound field computation. Indeed, it is evident that, in the rear zone of the acoustic field, rhoCentralFoam is not in agreement with the data reported in [16].

Heat Transfer
A uniform temperature boundary condition was employed on the square cylinder surface to investigate the forced   [17,18], for Re = 90, 120, 150. Despite Fig. 4 and Fig. 5 may suggest that our solver underestimates the Nu, it is clear that Nu-Re trend is correctly reproduced. It is worth noting that also the Strouhal number and the aerodynamic coefficients show the same trend observed in literature, as noticeable in Fig. 6. In our opinion the slight deviation between our results and the reference data is due to the compressibility e↵ects introduced by our solver which are not take into account in the incompressible approach used in [17,18].

Thermal effects on aeolian tones
The acoustic field radiated from a square cylinder in a cross laminar flow is generated by the periodical vortical structures shedding. This phenomenon causes a pressure fluctuation on the cylinder surfaces leading to a unsteady lift/downforce generation. Also the drag is influenced by the Karman vortex street, showing a downstream/upstream pulsation. These perturbations have a sound quadrupole nature; however, the lift dominance yields a typical bipolar acoustic field, [19]. Moreover, it is important to remark that an acoustic dipole generates a directional perturbation: it depends on the cosine of the angle between the observer and the dipole axis, [20]. Because of the alternation of vortices detachment and their location downstream to the cylinder, the dipole axis is not aligned with the undisturbed flow direction. As already noted in [8][9][10], the wall temperature increase produces a drag force growth and a instantenuos lift force reduction. It is recognized that this evidence is related to a thermal e↵ect on the pressure and to a local dynamic viscosity modification. In this work we separate the pressure and viscosity contributions in order to evalute their reciprocal influence. For this reason, in all the computations presented in this section, µ is fixed to a constant value. Furthermore, to estimate the acoustic pressure modification we consider a circle having r = 40 D to avoid an excessive sound wave decay in the far field. This condition enabled us to fix the computational domain extension at r = 150 D, reducing the cells number to: N c ' 3 · 10 6 . Looking at Fig. 10 is evident, for the isolated square cylinder, that the wall temperature raise produces, according to the literature, an increase of hC D i and a reduction of ∆C L and ∆C D as reported in Fig. 11. Thus, it is obviuous that on the square cylinder surface the pressure fluctuations are dumped for higher wall temperature values resulting in a far field sound abatament, as also depitected in Fig. 1, 7 and 8. For the isolated square cylinder, this behaviour was observed for all the Re numbers investigated. It is im- The last point is related to a greater impact of the adverse pressure gradient (due to thermal e↵ects) compared to convective terms. According to the sound generation process two di↵erent e↵ects act to modify the dipolar axis orientation: the shedding frequency and the vorticity di↵usion. For all the three Re numbers a monotone directivity diminishing should be expected because of the shedding frequency reduction by the wall temperature e↵ect, see Fig. 6. Indeed the vortices distance grows in the main flow direction. Dif-

Conclusions
This paper addresses the application of a finite-volume solver for compressible Navier-Stokes equations, named caafoam, which is able to fully resolve aeroacoustic waves propagation, heat transfer and their interaction with a direct numerical simulation approach. Our solver adopts low-storage explicit Runge-Kutta schemes for time-integration. KNP scheme was used for convective terms, while standard central schemes were considered for di↵usive ones. The non-reflective boundary treatment was achieved using an artificial sponge-layer approach. It was selected because it is simple to be coded, robust, non-sti↵, In this paper the reliability of caafoam was proved for both heat transfer and aeroacoustic problems; in particular aeolian tones radiated from rigid bodies in an uniform velocity inlet flow were treated. In this context is important to put in evidence that standard OpenFOAM solvers (equipped with non-reflective boundary treatment) are not able to correctly resolve the acoustic waves. Furthermore we have investigated the impact of wall boundary conditions for energy equation on sound propagation in fluids. We have showed that the wall temperature increment reduces the aeroacoustic sound due to lift oscillations amplitude drop. This e↵ect is clearly correlated to the pressure field modification since viscosity e↵ects were deliberately removed. On the other hand the wall temperature a↵ects the sound directivty in a contrasting way: when the convective e↵ects are not significative the directivity increases with the temperature, while at higher convective regimes we obtain the opposite phenomenon.
Ongoing work is devoted to the evaluation of the influence of the dynamic viscosity to the aeroacoustic sound decay and to enlarge similar investigations to other flow regimes.

Acknowledgements
We acknowledge the CINECA Award N. HP10CQ7QS7 YEAR 2019 under the ISCRA initiative, for the avail-