Development of Heat Exchanger Modeling Capability for a Finite-Volume Aeroelasticity Solver

Mitchell, Sam; Tartaruga, Irene; Stapelfeldt, Sina

doi:10.1115/1.4067380

Abstract

Heat exchangers are frequently used in aero-engines and are known to significantly affect the surrounding steady and unsteady flow. In certain applications, they may thereby also influence the aeroelastic stability of upstream or downstream components, but there is limited research on this in the public domain. This article aims to demonstrate the influence of heat exchangers on unsteady flows relevant to aeroelastic problems. This is achieved by developing heat exchanger modeling capability for an in-house finite volume aeroelasticity solver, for which heat exchanger is represented as a porous medium, as this is the established approach in existing aerodynamic studies using commercial computational fluid dynamics (CFD) software. The governing equations for a Darcy–Forchheimer porous media model suitable for unsteady and compressible flows are presented, which are derived by the application of volume-averaging theory to the Navier–Stokes equations. The implementation of this model within the time integration method used for the solver is then described and verified by comparison of results for steady flows against an established commercial CFD solver, where close agreement between both in-house and commercial solvers has been observed. Lastly, a preliminary demonstration of the capability to model unsteady heat exchanger flows is presented by application to an aeroacoustic problem, where the interaction of the pressure waves and the heat exchanger is investigated.

1 Introduction

Computational aeroelasticity is a crucial aspect of modern aircraft engine design. It involves coupling aerodynamic and structural simulations to predict complex and inherently unsteady aeromechanical phenomena such as compressor flutter and forced response. It is known that upstream positioned heat exchangers shall significantly affect the flow experienced by a downstream component; for example, in a heat pump where the fan is installed downstream of the heat exchanger, the interaction between the nonuniform turbulent heat exchanger exit flow field and the fan has been found to increase the emitted sound [1,2]. Furthermore, it is known that intake flow nonuniformity is likely to increase aerodynamic forcing on an aero-engine fan blade and may lead to large vibration responses and failure due to high-cycle fatigue, for example, as investigated in Refs. [3,4]. The presence of heat exchangers may also influence fan or compressor aerodynamic damping. From studies on fan flutter stability [5,6], it is known that changes in acoustic impedance upstream of the fan blade, caused for example by an intake opening, decrease the stability of the fan by reflecting acoustic waves generated by the blade vibration. This is shown in Fig. 1 from Ref. [7], where the aerodynamic damping of a fan blade without intake, with a short and with a long intake is compared for a range of fan speeds. A similar effect is created with porous liners, used for sound control and to influence aeroelastic stability [7–9]. Since heat exchangers also significantly alter the upstream flow, it can be expected that they similarly affect aerodynamic damping. To the authors’ knowledge, the influence of heat exchangers on the aeroelastic behavior of downstream aero-engine turbomachinery components has however not yet been investigated in existing literature. This highlights an important emerging area of research, especially since some modern aero-engines may feature this configuration; for example, precooled engines in which the turbo core is positioned downstream of a precooler heat exchanger [10,11].

Fig. 1

View large Download slide

Influence of intake length on aerodynamic damping for a range of fan speeds, from Ref. [7]

One approach to including heat exchangers in a computational fluid dynamics (CFD) simulation would be to directly simulate the entire heat exchanger and the flow within, although the mesh required to capture the intricate geometry of the heat exchanger core would be impractically large and complex to generate. It is therefore more common to instead model the heat exchanger as a porous medium, enabling the capture of its macroscopic effects with a practical mesh size. This was the method chosen by Missirlis et al. [12], who investigated a U-shape heat exchanger within the exhaust nozzle of a turbofan engine. They determined a pressure drop law for their model by experiments and then applied it to two-dimensional computational simulations. More recently, Missirlis et al. [13] developed a porous media model based on a modified form of the Darcy–Forchheimer equation and applied it to steady flow through an aero-engine heat exchanger. In both cases, close agreement between experimental and numerical results for pressure drop and flow-field properties, demonstrating the effectiveness of the porous media-based method for modeling steady flows through heat exchangers.

Other researchers who have used this method include Alshare et al. [14], who studied a problem representative of a shell-and-tube heat exchanger. Their results suggested that the use of the porous media model is an adequate approximation fully resolving the heat exchanger geometry. Wang et al. [15] used porous media modeling to study a plate-fin heat exchanger, where they found their simulations valuable toward suggesting improvements to its design. Yang et al. [16] showed good agreement between porous media modeling and experimental data for predicting pressure drop and heat transfer of a rod-baffle shell-and-tube heat exchanger. Musto et al. [17] investigated a heat exchanger within an oil cooler of a propeller aircraft. They showed close reproduction of experimental pressure drop and heat rejection trends using simulations with porous media modeling. An and Kim [18] applied a porous media model to efficiently analyze a fin-and-tube heat exchanger relevant to air-conditioning systems. Ding et al. [11] applied a porous media model to study a compact tube heat exchanger for a hypersonic precooled aero-engine, for which they analyzed the effects of varying inlet pressure distortions. In each of the discussed studies, a commercial CFD solver was chosen, most often fluent. This emphasizes the existing ability of such solvers to include porous media models. These studies, however, do not involve any aeroelastic applications; they are limited to aerodynamic and thermal analyses of steady flows through the heat exchanger, with little to no consideration of unsteady flows.

There are of course many other applications of porous media modeling apart from heat exchangers, and the resultant advancements in modeling theory are also valuable to our application. In particular, models suitable for unsteady and compressible flows are of particular relevance to our work, since turbomachinery aeroelasticity problems regularly involve such flows. A particularly relevant example is the study of Jarauta et al. [19], who developed a Darcy–Forchheimer-based porous media model for steady compressible flows by extending volume-averaging theory to the compressible Navier–Stokes equation. They validated their model using experimental data and applied it to accurately predict two-dimensional in-plane and through-plane channel flows relevant to fuel cells. Prokein et al. [20] developed a porous media model suitable for unsteady and compressible flows with heat transfer and applied it to study transpiration cooling in supersonic flows.

The aim of this work is to develop heat exchanger modeling capability for an aeroelastic solver and apply it to both steady and unsteady flow cases. This is toward the longer-term goal of simulating industrial turbomachinery aeroelasticity problems involving heat exchangers. To do this, we shall utilize the existing foundations of porous media modeling for heat exchanger aerodynamic analysis, along with recent advancements of porous media modeling theory for unsteady compressible flows. First, we shall present a porous media model suitable for unsteady compressible flows, based on the volume averaging approach [21]. To derive it, we extend Jarauta et al.’s [19] recent porous media model for steady compressible flows, by adding the unsteady temporal derivative to the continuity and momentum equations, and including the energy equation. This initial study will focus mainly on the mechanical effects of the heat exchanger, and not heat transfer; thus, a local thermal equilibrium (LTE) shall be assumed for this model.

We shall then describe the implementation of this model within the numerical scheme of an established aeroelastic solver, au3d [22]. This is an in-house code for the solution of turbomachinery aeroelasticity problems such as flutter and forced response [5,23]. Its aerodynamic solver is capable of solving steady and unsteady viscous compressible flows on unstructured hybrid grids. To do this, the compressible Favre-averaged Navier–Stokes equations together with the Spalart–Allmaras turbulence model are solved with an implicit dual time-stepping scheme and a node-centered finite-volume scheme for spatial discretization [22]. The implementation of the presented compressible porous media model within au3d shall next be verified using steady flows around two simplified yet industrially relevant geometries representing a ducted heat exchanger. The obtained results shall be compared against results from fluent, as it is known to perform well for steady situations [11,13,15–17], and we thus can expect direct agreement with its porous media model for these cases.

We shall then present an unsteady application, where we shall investigate the interaction of acoustic waves and the modeled porous media. This is already an important field of study, such as for application to sound absorption [24] and geotechnics [25], including the use of volume-average-based porous media modeling [26]. This study shall be limited to single-harmonic, planar acoustic sources, which provides a controlled setting for initial study. We shall apply an established eigenmode analysis method [27], enabling the detection and quantification of the prevalent acoustic modes, and so phenomena such as attenuation and reflection caused by the presence of the modeled porous medium.

2 Flow Model

As mentioned in the previous section, the established approach to including heat exchangers in CFD studies is to model them as a porous medium. In this section, we use the volume-averaging theory (VAT) to present a Darcy–Forchheimer porous media model suitable for unsteady, compressible flows.

2.1 Volume-Averaging Theory.

To begin the volume-averaging process [21,28,29], we first consider the porous medium as a two-phase flow: one solid phase (⁠ $s$ ⁠) representing the structure and one fluid phase ( $f$ ⁠) representing the pores. As shown in Fig. 2, the pore size $d$ is assumed to be much smaller than the macroscopic system length $L$ ⁠.

Fig. 2

View large Download slide

Illustration of porous medium from different length scales

At any point in space within the porous medium (⁠ $x$ ⁠), a small averaging volume $V$ with radius $r_{0}$ may be defined, also known as the representative elementary volume (REV). It is comprised of a fluid component $V_{f} (x)$ and a solid component $V_{s} (x)$ ⁠, such that $V = V_{f} + V_{s}$ ⁠.

There are two ways to obtain a volume-average value of property existing within medium, both involving a volume integral across the REV. These are shown for a generic fluid property

ϕ_{f}

in Eqs. and , respectively, and are known as the phase average

⟨ ϕ_{f} ⟩

and intrinsic phase average

⟨ ϕ_{f} ⟩^{f}

⁠. These two averages are related by a material property known as the fluid porosity

ε_{f} = ⟨ ϕ_{f} ⟩ / ⟨ ϕ_{f} ⟩^{f} = V_{f} / V

[28,29].

⟨ ϕ_{f} ⟩ = \frac{1}{V} \int_{V_{f}} ϕ_{f} d V

(1)

⟨ ϕ_{f} ⟩^{f} = \frac{1}{V_{f}} \int_{V_{f}} ϕ_{f} d V

(2)

The spatial averaging theorem [28] shown in Eq. relates the average of a spatial gradient to the spatial gradient of an average. Here,

A_{f s}

is the total surface area comprising the boundaries between fluid and solid phases within the averaging volume, and

n_{f s}

is the unit normal vector pointing from the fluid to solid phase.

⟨ \nabla ϕ_{f} ⟩ = \nabla ⟨ ϕ_{f} ⟩ + \frac{1}{V} \int_{A_{f s}} ϕ_{f} n_{f s} d A

(3)

Also of importance shall be the relations between products of average quantities [19]. With the reasonable assumption that the pore size is much smaller than the size of the averaging volume (⁠

d << r_{0}

⁠), Eqs. and are valid. Here,

ψ

represents another generic quantity within the porous medium.

⟨ ϕ_{f} ψ_{f} ⟩ \approx \frac{1}{ε_{f}} ⟨ ϕ_{f} ⟩ ⟨ ψ_{f} ⟩

(4)

⟨ ϕ_{f} ψ_{f} ⟩^{f} \approx ⟨ ϕ_{f} ⟩^{f} ⟨ ψ_{f} ⟩^{f}

(5)

2.2 Application to the Conservation Equations.

We now apply the volume-averaging theory summarised above to the compressible Navier–Stokes system of equations, to obtain a set of governing equations that model the volume-averaged behavior of the fluid within the medium.

2.2.1 Mass.

First, the phase average of the continuity equation is taken, as shown in Eq.

\frac{1}{V} \int_{V_{f}} \frac{\partial ρ_{f}}{\partial t} + \nabla \cdot (ρ_{f} u_{f}) d V = 0

(6)

Since the medium is assumed rigid, its volume is independent of time, which means the temporal derivative of

ρ

may be removed from the integration. This, along with the application of the spatial averaging theorem (Eq. ) to the spatial derivative, leads to Eq.

\frac{\partial ⟨ ρ_{f} ⟩}{\partial t} + \nabla \cdot ⟨ ρ_{f} u_{f} ⟩ + \int_{A_{f s}} n_{f s} \cdot u_{f} d A = 0

(7)

As the solid (⁠

σ

⁠) phase is impermeable, the surface integral within the above equation is zero. This along with applying Eq. results in Eq.

\frac{\partial ⟨ ρ_{f} ⟩}{\partial t} + \nabla \cdot (\frac{1}{ε_{f}} ⟨ ρ_{f} ⟩ ⟨ u_{f} ⟩) = 0

(8)

As will be discussed later, it is beneficial to express average velocity in terms of the intrinsic phase average, and average density in terms of the phase average [19]. This results in the volume-averaged continuity equation that shall be used within the porous media model.

\frac{\partial ⟨ ρ_{f} ⟩}{\partial t} + \nabla \cdot (⟨ ρ_{f} ⟩ ⟨ u_{f} ⟩^{f}) = 0

(9)

2.2.2 Momentum.

Next, the conservation of average momentum within the porous medium is derived, by applying VAT to the Navier–Stokes equation. This leads to Eq. , with average pressure and viscous stress given by Eqs. and , respectively.

\begin{aligned} \frac{1}{ε_{f}} \frac{\partial ⟨ ρ_{f} ⟩ ⟨ u_{f} ⟩}{\partial t} + \nabla \cdot (\frac{1}{ε_{f}^{2}} ⟨ ρ_{f} ⟩ ⟨ u_{f} ⟩ \otimes ⟨ u_{f} ⟩) \\ + \frac{1}{V} \int_{A_{f s}} (ρ_{f} u_{f} \otimes u_{f}) n_{f s} d A = \nabla \cdot (- ⟨ p_{f} ⟩ I + ⟨ σ_{f} ⟩) \\ + \frac{1}{V} \int_{A_{f s}} (- p_{f} I + σ_{f}) n_{f s} d A \end{aligned}

(10)

⟨ p_{f} ⟩ = ⟨ ρ_{f} R T_{f} ⟩

(11)

σ_{f} = μ_{f} (\nabla ⟨ u_{f} ⟩ + (\nabla ⟨ u_{f} ⟩)^{T}) + λ_{f} (\nabla \cdot ⟨ u_{f} ⟩) I

(12)

The surface integral on the left hand side (LHS) must be equal to zero due to the impermeability of the solid phase. The surface integral on the right hand side (RHS) is however unknown and must be approximated. Following Whitaker [29] and Jarauta et al. [19], the Darcy–Forchheimer equation that governs flow through a porous medium is now used as a closure. Here,

K_{D}

and

K_{F}

are tensors containing the viscous and inertial permeability coefficients of the porous medium, respectively.

\begin{aligned} \frac{1}{V} \int_{A_{f s}} (- p_{f} I + σ_{f}) n_{f s} d A \\ \approx - μ_{f} K_{D}^{- 1} ⟨ u_{f} ⟩ - K_{F}^{- 1} ⟨ ρ_{f} ⟩ | ⟨ u_{f} ⟩ | ⟨ u_{f} ⟩ \end{aligned}

(13)

Substitution of the above equation into Eq. and expression of average velocity in terms of the intrinsic phase average give the volume-averaged Navier–Stokes equation shown in Eq. .

\begin{aligned} \frac{\partial ⟨ ρ_{f} ⟩ ⟨ u_{f} ⟩^{f}}{\partial t} + \nabla \cdot (⟨ ρ_{f} ⟩ ⟨ u_{f} ⟩^{f} \otimes ⟨ u_{f} ⟩^{f}) \\ = \nabla \cdot (- ⟨ p_{f} ⟩ I + ⟨ σ_{f} ⟩^{f}) - μ_{f} K_{D}^{- 1} ε_{f} ⟨ u_{f} ⟩^{f} \\ - K_{F}^{- 1} ⟨ ρ_{f} ⟩ | ε_{f} ⟨ u_{f} ⟩^{f} | ε_{f} ⟨ u_{f} ⟩^{f} \end{aligned}

(14)

It should be noted also that due to the choice of averaging variables $, ⟨ ρ_{f} ⟩, ⟨ u_{f} ⟩$ ⁠, experimental permeability data should also be corrected as in the study of Jarauta et al. [19], for example, $K_{D_{model}} = ε_{f} K_{D_{experiment}}$ and $K_{F_{model}} = ε_{f}^{2} K_{F_{experiment}}$ ⁠.

2.2.3 Energy.

Third, the volume-averaging procedure is applied to the conservation equation of total energy in the fluid phase. This is shown below, with the dependency of total energy on temperature exposed. Several intermediate steps have been omitted for brevity. Heat transfer between the fluid and solid phases is dependent on the interphase temperature difference and is represented using the volumetric heat transfer coefficient

h_{v}

[20,30].

\begin{aligned} \frac{\partial}{\partial t} (⟨ ρ_{f} ⟩ (c_{v_{f}} ⟨ T_{f} ⟩^{f} + \frac{1}{2} ⟨ u_{f} ⟩^{2^{f}})) + \nabla \cdot \\ (⟨ ρ_{f} ⟩ ⟨ u_{f} ⟩^{f} (c_{v_{f}} ⟨ T_{f} ⟩^{f} + \frac{⟨ p_{f} ⟩}{⟨ ρ_{f} ⟩} + \frac{⟨ u_{f} ⟩^{2^{f}}}{2})) \\ - \nabla \cdot (ε_{f} k_{f} \nabla ⟨ T_{f} ⟩^{f} - ⟨ σ_{f} ⟩ \cdot ⟨ u_{f} ⟩^{f}) = h_{v} (⟨ T_{s} ⟩^{s} - ⟨ T_{f} ⟩^{f}) \end{aligned}

(15)

To derive the conservation of energy within the entire spatially averaged medium, the conservation of energy within the solid phase must be taken into account also. The energy equation for this phase is as follows [20]:

ρ_{s} c_{p_{s}} \frac{\partial T_{s}}{\partial t} = \nabla \cdot (k_{s} \nabla T_{s})

(16)

Taking the spatial average of this leads to the following. Note here that since this phase is assumed to be a homogeneous solid, we consider the density of the porous medium

ρ_{s}

as a constant.

ε_{s} ρ_{s} c_{p_{s}} \frac{⟨ T_{s} ⟩^{s}}{\partial t} = \nabla \cdot (ε_{s} k_{s} \nabla ⟨ T_{s} ⟩^{s}) + h_{v} (⟨ T_{f} ⟩^{f} - ⟨ T_{s} ⟩^{s})

(17)

In this initial study, we shall next assume the simplified situation of LTE, so that the temperature of the fluid and solid phases are equal (⁠

⟨ T_{s} ⟩ = ⟨ T_{f} ⟩

⁠). This assumption allows the combination of Eqs. and into a single conservation equation for the volume-average energy in both phases, as shown in Eq. [30]. Here, the volume-average total energy within the fluid and solid phases is shown in Eqs. and , respectively.

⟨ E_{f} ⟩^{f} = c_{v_{f}} ⟨ T_{f} ⟩^{f} + \frac{⟨ u_{f} ⟩^{2^{f}}}{2}

(18)

⟨ E_{s} ⟩^{s} = c_{p_{s}} ⟨ T_{s} ⟩^{s}

(19)

\begin{aligned} \frac{\partial}{\partial t} (⟨ ρ_{f} ⟩ ⟨ E_{f} ⟩^{f} + ε_{s} ρ_{s} ⟨ E_{s} ⟩^{s}) \\ + \nabla \cdot (⟨ u_{f} ⟩^{f} (⟨ ρ_{f} ⟩ ⟨ E_{f} ⟩^{f} + ⟨ p_{f} ⟩)) \\ - \nabla \cdot (k_{s} ε_{s} \nabla ⟨ T_{s} ⟩^{s} + k_{f} ε_{f} \nabla ⟨ T_{f} ⟩^{f} + ⟨ σ_{f} ⟩^{f} \cdot ⟨ u_{f} ⟩^{f}) = 0 \end{aligned}

(20)

The conduction terms in Eq. may also be grouped into a single term with an effective conductivity

k_{e f f}

as shown below.

\begin{aligned} ε_{s} k_{s} ⟨ T_{s} ⟩^{s} + ε_{f} k_{f} ⟨ T_{f} ⟩^{f} = (ε_{s} k_{s} + ε_{f} k_{f}) ⟨ T_{f} ⟩^{f} = k_{eff} ⟨ T_{f} ⟩^{f} \end{aligned}

(21)

To summarize, it is the VANS equations (9), (14), and (20) that govern the conservation of spatially averaged mass, momentum, and energy within the porous medium. It is therefore these equations that shall be solved numerically within porous medium regions of the CFD simulation domain.

2.3 Boundary Conditions

2.3.1 Porous–Fluid Interfaces.

To model interfaces between regions of fluid and porous medium, such as inlets and outlets to a heat exchanger, we extend the single domain method [19] to include the unsteady LTE energy equation also. The first step to this is to consider that the fluid fluxes through the interface must be conserved, which provides the relationships between the standard fluid variables on the fluid side of the interface to the spatially averaged variables on the porous side. Consider for example the mass and convective momentum fluxes at an interface

Σ

as shown below in Eqs. and .

ρ u \cdot n |_{Σ} = ⟨ ρ_{f} ⟩ ⟨ u ⟩^{f} \cdot n |_{Σ}

(22)

(ρ u \otimes u) |_{Σ} = (⟨ ρ_{f} ⟩ ⟨ u_{f} ⟩^{f} \otimes ⟨ u_{f} ⟩^{f}) |_{Σ}

(23)

Solving the above two equations results in the continuity of density and velocity normal component at the boundary [19]. In our case, we also assume a continuity of the shear stress tensor and thus the tangential component of the velocity also. We have found this assumption suitable so far for our applications to primarily axial flows, though in other situations, it may be worthwhile to consider additional constraints on the stress tensor [19]. This leads to the following boundary conditions for density and velocity:

ρ |_{Σ} = ⟨ ρ_{f} ⟩ |_{Σ}

(24)

u |_{Σ} = ⟨ u_{f} ⟩^{f} |_{Σ}

(25)

By applying the same technique to the convective energy flux in Eq. , we can determine continuity of energy and thus temperature at the boundary as shown below. This approach to the thermal boundary condition is suitable for our current applications to LTE problems. For future studies involving non-LTE problems, the use of explicit and more detailed thermal boundary conditions should be investigated, such as in existing studies [20].

ρ u E \cdot n |_{Σ} = ⟨ ρ_{f} ⟩ ⟨ u_{f} ⟩^{f} ⟨ E_{f} ⟩^{f} \cdot n |_{Σ}

(26)

T |_{Σ} = ⟨ T_{f} ⟩^{f} |_{Σ}

(27)

It is now clear that the particular choice of spatially averaged variables $(⟨ ρ_{f} ⟩, ⟨ u_{f} ⟩^{f}, ⟨ T_{f} ⟩^{f})$ in which to cast the governing equations of the model means that we may achieve a continuous variation of flow variables across the porous–fluid interface with the single domain method, as there is no need to apply a porous jump condition at the interfaces. This is in fact the main reason behind this choice of variables; options requiring a porous jump condition shall invoke discontinuities at the interface that may be detrimental to numerical stability [19]. The conservative finite-volume spatial discretization used by au3d respects the continuity of fluxes such as those shown in Eqs. (22), (23), and (26) [22]; thus, its utilization at the porous–fluid interface implicitly imposes these boundary conditions in a straightforward way.

2.3.2 Solid Walls.

Currently, either inviscid or viscous wall boundary conditions may be enforced at solid surfaces of the heat exchanger exterior, for which au3d’s existing functionality is utilized. For the latter, the option of employing a wall model to reduce resolution demand at the surfaces is also available. Inlet and outlet surfaces to porous regions of the domain may feature either conformal or nonconformal mesh nodes, thanks to au3d’s interpolation plane functionality.

3 Numerical Implementation

In this section, we shall discuss the implementation of the above porous media model within au3d. This involves primarily the discretization of the volume-averaged Navier–Stokes (VANS) equations: Eqs. (9), (14), and (20) in space and time. These VANS equations are similar to the Favre-averaged Navier–Stokes equations already solved by au3d [22]; therefore, much of the existing code features may be utilized to solve the VANS equations also. The following section therefore closely follows the original presentation of au3d’s numerical scheme [22], though the required modifications to solve the VANS equations instead shall be highlighted.

First, Eqs. , , and are reexpressed in integral form and gathered into a matrix system as shown below. Here, the vector

W = {[ρ, ρ u, ρ E]}^{T}

contains the conservative flow variables, and the functions

F

and

S

contain the flux and source terms respectively. Notice that in Eq. , the LTE assumption has resulted in the inclusion of the solid phase energy in the time derivative. To include this in the following implementation, we define a function

\tilde{W}

that provides a set of alternative conservative variables to reflect this, which are then easily calculable from the standard conservative variables, i.e.,

\tilde{W} = {[ρ, ρ u, ρ E + ε_{s} ρ_{s} E_{s}]}^{T}

⁠.

\frac{\partial}{\partial t} \int_{Ω} \tilde{W} (W) d Ω + \oint_{Γ} F (W) \cdot n d Γ = S (W)

(28)

A spatial discretization of the above system is then performed using au3d’s existing finite-volume method; details of which can be found in the original publication [22]. Importantly, it follows an edge-based approach, is suitable for unstructured grids, and includes the option of first- or second-order spatial accuracy. The above system is thus transformed to an ordinary differential equation (ODE) in time, as shown below for an internal mesh node

_{I}

⁠, where

Ω_{I}

is the control volume and

R_{I}

is the residual operator.

\frac{d}{d t} (Ω_{I} {\tilde{W}}_{I} (W)) + R_{I} (W) = 0

(29)

3.1 Temporal Discretization.

We make use also of the existing implicit dual-time temporal discretization method used in au3d. Some of the main steps are covered here, though more details can be found in the original au3d publication [22] or other prior applications of this method [31,32]. The first step is to consider a fully implicit evaluation of the above system and approximate the time derivative with a second-order backward difference discretization with step

Δ t

⁠, where

n

denotes the discrete-time level [31,32].

\begin{aligned} \frac{3 Ω_{I}^{n + 1} {\tilde{W}}_{I} (W^{n + 1}) - 4 Ω_{I}^{n} {\tilde{W}}_{I} (W^{n}) + Ω_{I}^{n - 1} {\tilde{W}}_{I} (W^{n - 1})}{2 Δ t} + R_{I} (W^{n + 1}) = 0 \end{aligned}

(30)

Next, a derivative of

W^{n + 1}

with respect to a pseudo-time

t^{*}

is introduced as shown below. The term

R^{*}

⁠, known as the unsteady residual, is shown in Eq. . The steady-state solution of Eq. is also the solution of Eq. ; thus, one can obtain

W^{n + 1}

by solving Eq. , for example, with a time-marching approach [31].

\frac{\partial W_{I}^{n + 1}}{\partial t^{*}} + \frac{1}{Ω_{I}^{n + 1}} R_{I}^{*} (W^{n + 1}) = 0

(31)

\begin{aligned} R_{I}^{*} (W^{n + 1}) & = \frac{3 Ω_{I}^{n + 1} {\tilde{W}}_{I} (W^{n + 1}) - 4 Ω_{I}^{n} {\tilde{W}}_{I} (W^{n}) + Ω_{I}^{n - 1} {\tilde{W}}_{I} (W^{n - 1})}{2 Δ t} \\ + R_{I} (W^{n + 1}) = 0 \end{aligned}

(32)

It is now important to emphasize the first modification from au3d’s existing method, which is the inclusion of the additional source term

S^{PMM}

present in the VANS equations. This is now included in the residual in addition to centrifugal and coriolis forces

S^{c}

⁠, as shown in Eq. . Currently, this comprises of only a momentum component, as can be seen in Eqs. and , though future work will include heat sources also. It should be noted here that the terms within the flux operator

F_{I}

that depend on thermal conductivity must now be computed using

k_{e f f}

as shown in Eq. .

R_{I} (W^{n}) = F_{I} (W^{n}) + S_{I}^{c} (W^{n}) + S_{I}^{PMM} (W^{n})

(33)

S_{I}^{PMM} = [\begin{matrix} 0 \\ - μ K_{D}^{- 1} ε_{f} u - K_{F}^{- 1} ρ | ε_{f} u | ε_{f} u \\ 0 \end{matrix}]

(34)

Now, an implicit evaluation of Eq. is performed, including a first-order backward difference approximation of the pseudo-time derivative, where

m

represents the discrete pseudo-time level. As is conventional, we now use

W_{I}

to refer to the approximation of

W_{I}^{n + 1}

[31].

\frac{Δ W_{I}}{Δ t^{*}} = \frac{W_{I}^{m + 1} - W_{I}^{m}}{Δ t^{*}} = - \frac{1}{Ω_{I}} R_{I}^{*} (W^{m + 1})

(35)

An approximation of the unsteady residual at the next pseudo-time level

m + 1

can be obtained via linearization around the current level as shown below [22,31].

R^{*} (W^{m + 1}) \approx R^{*} (W^{m}) + \frac{\partial R^{*}}{\partial W} Δ W

(36)

The Jacobian of the unsteady residual may be calculated as the derivative of Eq. with respect to the conservative variables

W

⁠, leading to the following equation. The next modification is now apparent; the presence of the Jacobian of the alternative variable function

\tilde{W}

with respect to

W

⁠. This matrix is calculable analytically and thus does not introduce significant computational demand to obtain.

\frac{\partial R^{*}}{\partial W} = \frac{3 Ω}{2 Δ t} \frac{\partial \tilde{W}}{\partial W} + \frac{\partial R}{\partial W}

(37)

The next modification is also visible in the Jacobian of the residual as shown below; there is now an additional component including the derivative of the VANS source term. This matrix is also calculable analytically, by considering Eq.

\frac{\partial R}{\partial W} = \frac{\partial F}{\partial W} + \frac{\partial S^{c}}{\partial W} + \frac{\partial S^{PMM}}{\partial W}

(38)

In the same way as shown in the original publications [22,31], substitution of Eqs. and into Eq. leads to the linear matrix system shown below. The solution of this system provides

W^{m + 1}

⁠, and thus, a converged solution will be equal to

W^{n + 1}

⁠. We employ au3d’s existing methodology to solve this system, consisting of a Newton iteration scheme including internal point-Jacobi iterations [22].

[(\frac{Ω}{Δ t^{*}} + \frac{3 Ω}{2 Δ t} \frac{\partial \tilde{W}}{\partial W}) + \frac{\partial R}{\partial W}] Δ W = - R^{*} (W^{m})

(39)

4 Results for Steady Flows

Now we shall present some recent results obtained from the application of the above Darcy–Forchheimer LTE porous media model in au3d. Three simplified geometries shall be used, which are representative of a general flow through a duct containing a heat exchanger. The first two cases involve steady flows, for which we shall compare results from both au3d and fluent. For these two cases, we can expect a near-exact agreement between both solvers, thus obtaining close results will serve as a strong verification of the implementation of porous media modeling capability in au3d.

4.1 Case 1: Duct With Porous Plug.

The porous-plug case provides a highly simplified setting, where solely the effect of the porous medium may be examined and any uncertainty due to complex geometry and flow features minimized. Figure 3 shows the geometry of this case; a cuboid duct of length $L = 1.6$ m is divided into four volume regions. The second region is filled with a porous medium that spans the entire $y - z$ cross section, representing a simplified heat exchanger. The freestream flow direction is aligned with the $x$ -axis. The pressure-driven flow is defined by an inlet at atmospheric conditions (total pressure of $1.013 \times 10^{5}$ Pa and total temperature of 288 K), and an outlet static pressure of $6 \times 10^{4}$ Pa. The Reynolds and Mach numbers based on the freestream flow upstream of the heat exchanger and the domain length $L$ are approximately $Re = 6 \times 10^{6}$ and $M = 0.16$ ⁠. These conditions are representative of those experienced in industrial heat exchanger flows, for example, the Mach number is within proximity to existing aero-engine studies [11,17]. This means also that these conditions are most suitable for our verification using results from fluent, since it is at these conditions for which the accuracy of the porous media model of fluent is most validated in existing studies.

Fig. 3

View large Download slide

Domain geometry of case 1: (⁠ $a$ ⁠) $x - y$ mid-plane and (⁠ $b$ ⁠) $y - z$ plane at $x = 1.0$

The porous medium in this case is defined by the viscous and inertial resistance (inverse permeability) coefficients shown in Table 1, and a porosity value of $ε = 0.7$ ⁠. While in practice these values must typically be evaluated with an experimental curve-fitting procedure [13,17], for our current verification purposes, we have chosen values that are representative of an industrial aero-engine heat exchanger, for example, they are the same order as in studies of industrial cases [13,17]. Nevertheless, these coefficients will vary between the specific heat exchanger core configuration and conditions, and so the robustness of the current model to accommodate variation in these coefficients is demonstrated in later sections. It is also important that the inverse permeability in the $y$ and $z$ directions is much greater than in the $x$ direction, representing an effectively solid medium in $y$ and $z$ ⁠.

Table 1

Inverse permeability coefficients of applied porous medium

	Inverse permeability
	Viscous (⁠ $K_{D}^{- 1}$ ⁠) $(1 / m^{2})$	Inertial (⁠ $K_{F}^{- 1}$ ⁠) (1/m)
xx	1.11 × 10⁷	70
yy	2.111 × 10¹¹	0
zz	2.111 × 10¹¹	0

4.1.1 Numerical Setup.

As shown in Fig. 4, a three-dimensional structured mesh of 761117 hexahedral cells was generated using the icem meshing software. The mesh for each volume region was generated independently, and the whole-domain mesh was subsequently obtained by joining these four meshes. The mesh is nonconformal at the interfaces of each volume region, as intended to maintain generality of application. Interpolation planes are employed to transfer the solution at the three interfaces between volume regions. The same mesh was used for both au3d and fluent simulations.

Fig. 4

View large Download slide

Case 1 mesh

4.1.2 Analysis of Steady Flow-Field.

Figure 5 shows the variation in primitive fluid variables along the $x$ -axis centerline of the cuboid channel. It is clear to see that the au3d solution agrees very closely with the fluent solution. The large pressure drop through the porous media (⁠ $0.5 < x < 0.7$ ⁠) is visible in both Figs. 5 and 6. A corresponding increase in velocity and a decrease in density and temperature are also visible in Fig. 5. It should be mentioned that while the compressible porous media model implemented in au3d solves for the intrinsic phase-averaged velocity $⟨ u_{f} ⟩^{f}$ ⁠, fluent’s porous media model here solves for the phase-averaged velocity $⟨ u_{f} ⟩$ ⁠. The agreement in velocities $⟨ u_{f} ⟩^{f}$ and $⟨ u_{f} ⟩$ of each solver, respectively, between $0.5 < x < 0.7$ shown in Fig. 5 and therefore means that the intrinsic phase-averaged velocity predicted in the porous medium by au3d and fluent differs from each other by a factor of porosity $ε$ (since $⟨ u_{f} ⟩ = ε ⟨ u_{f} ⟩^{f}$ ⁠). This should be expected considering the additional factor of $ε$ introduced by considering the density as a spatially variable function (rather than a constant) in the volume-averaging process presented in Sec. 2, as is appropriate for a compressible model [19,20]. Such close agreement in results shows that the losses due to the porous medium are modeled very similarly in au3d and fluent, thus serving as a good initial verification of the porous media modeling implementation in au3d.

Fig. 5

View large Download slide

Variation of case 1 flow variables along the $x$ -axis centerline

Fig. 6

View large Download slide

Case 1 static pressure in $x - y$ center plane: (⁠ $a$ ⁠) au3d and (⁠ $b$ ⁠) fluent

4.2 Case 2: Duct With Porous Cube.

This case features a channel of the same dimensions as the previous case, although now the porous medium does not comprise the entire cross section. As can be seen in Fig. 7, instead there is a cube of porous medium located within the center of the channel, with a fluid region around it acting as a bypass. The $x - y$ and $x - z$ walls of the porous cube are impermeable and viscous, representing a solid casing around the exterior of the cube. For this case, there are three volume regions, the second of which is the porous medium.

Fig. 7

View large Download slide

Domain geometry of case 2: (⁠ $a$ ⁠) $x - y$ mid-plane and (⁠ $b$ ⁠) $y - z$ plane at $x = 0.6$

Again a pressure-driven flow is enforced by maintaining constant values of total pressure and total temperature at the inlet and static pressure at the outlet. The same values of $p_{0} = 101300$ Pa and $T_{0} = 288$ K were maintained at the inlet, although this time a higher static pressure of $p = 80, 000$ was applied at the outlet. The freestream Reynolds and Mach numbers based on the inlet conditions and domain length are approximately $Re = 10^{7}$ and $M = 0.32$ ⁠.

4.2.1 Numerical Setup.

A hexahedral mesh was used for this case as well, for which the number of cells is 2,309,660. The mesh is again nonconformal at volume region interfaces. Figure 8 shows a cut of the mesh at the mid-plane, so that the cube of porous medium within the interior is also visible. The mesh regions surrounding the porous cube are significantly refined, especially at the surfaces and corners as there will be boundary layers here. Fluid regions feature wall-modeled boundary conditions and turbulence modeling with the Spalart–Allmaras model. Wall-modeled boundary conditions are also applied to the exterior side of the previously mentioned $x - y$ and $x - z$ walls of the porous cube in this case. Turbulence was not modeled within the porous region, neither was the near-wall boundary layer model active.

Fig. 8

View large Download slide

Case 2 mesh: cut at $x - y$ mid-plane

4.2.2 Analysis of Steady Flow Field.

As for the previous case, the variation of primitive flow variables along the $x$ -axis centerline of the channel is shown in Fig. 9 for this case. The results of au3d and fluent again agree well in general, further verifying the implementation of the porous media model within au3d. In particular, the pressure drop through the porous cube is clear. This similarity is also visible by comparing the velocity contours shown in Fig. 10. The differences in flow features in comparison to the previous case are also apparent in these figures. Figures 9 and 10 show that the blockage effect and localized pressure rise caused by the cubic porous medium and the duct walls result in diversion and acceleration of the flow around the cube. This has caused the separation of the boundary layer at the leading edges of the cube walls, with the separation bubble extending downstream over the walls toward the trailing edges, and subsequently, a wake is formed aft of the cube. The flow properties in the wake show minor differences between each solver; the most probable cause of which are theoretical differences in solver-specific formulation, such as the specific wall model used for the viscous boundary condition, or differences in the amount of numerical dissipation resulting from differences in the numerical schemes. This is thus not of concern in the verification of the heat exchanger modeling capability.

Fig. 9

View large Download slide

Variation of case 2 flow variables along the $x$ -axis centerline

Fig. 10

View large Download slide

Case 2 $x$ velocity in $x - y$ center plane: (⁠ $a$ ⁠) au3d and (⁠ $b$ ⁠) fluent

5 Results for Unsteady Flows

Now that the implementation of the porous media model in au3d has been sufficiently verified; next, we shall apply it to an unsteady flow case representing a simplified heat exchanger, the unsteady porous-plug case. The purpose of this case is to investigate the effects of a heat exchanger modeled as a porous medium on the unsteady flow field. It is similar to the steady porous-plug case, although the outlet pressure boundary condition now includes a sinusoidally fluctuating component, representing upstream propagating acoustic waves, for example, as a result of a downstream flutter event. The geometry is similar to the porous-plug case as shown in Fig. 3; the porous medium is located again between $0.5 \leq x \leq 0.7$ and covers the entirety of the duct cross section. The main difference is that the duct is now annular, with a radial extent of $0.3$ m, as shown in Fig. 11. This simplifies the application of the unsteady outlet boundary condition and the eigenmode analysis presented in later sections, since for the code used here these are implemented in cylindrical coordinates. The porous model coefficients are the same as shown in Table 1, since as mentioned previously these are representative values of industrial aero-engine heat exchangers.

Fig. 11

View large Download slide

Case 3 mesh and initial condition

5.1 Case 3: Unsteady Porous Plug

5.1.1 Numerical Setup.

A part-annulus structured mesh was generated for this case using an in-house meshing tool, as shown in Fig. 11. An axial cell size of $Δ x = 0.005$ m was chosen to give sufficient spatial resolution of the imposed acoustic waves. Similarly, a time-step of $Δ t = 1 \times 10^{- 5}$ s was used to guarantee sufficient temporal resolution of the waves.

The inlet and outlet boundaries are enforced using the Riemann-invariants method, as to minimize acoustic reflections. Furthermore, though not shown in Fig. 11, an additional coarse nozzle region was included upstream of the duct (before $x = 0$ ⁠), to further minimize any artificial reflections of the upstream-traveling acoustic wave. The outlet features a sinusoidally varying unsteady pressure component of amplitude 300 Pa, at a frequency of 1000 Hz, representing an acoustic source. Periodic boundaries are applied in the circumferential direction, since in this study the applied acoustic source is circumferentially uniform. Inviscid wall boundary conditions are applied at the hub and casing surfaces.

An initial condition was obtained from a prior steady simulation, the pressure contours of which are also shown in Fig. 11. The unsteady simulation was then performed across 50 acoustic periods, for which the solution was sampled for analysis at every time-step across the latter 25 periods.

5.1.2 Spatial Variation of Unsteady Pressure.

Figures 12 and 13 show the unsteady pressure component at the casing surface and through the centerline of the domain respectively. Both the standard configuration featuring the porous region and an equivalent duct flow without the porous region are shown. In the latter, the acoustic waves propagate upstream throughout the entire domain with low attenuation, while for the former, the acoustic waves are significantly attenuated as soon as the porous medium is encountered. Figure 14 shows the amplitude frequency spectrum of the unsteady pressure field at two axial positions, further showing that the frequency of the acoustic wave remains unchanged at 1000 Hz throughout the domain; yet, the amplitude is heavily diminished upstream of the porous medium.

Fig. 12

View large Download slide

Case 3 unsteady pressure at casing surface: (a) duct without porous region and (b) duct with porous region

Fig. 13

View large Download slide

Case 3 unsteady pressure through $x$ -axis centerline

Fig. 14

View large Download slide

Amplitude frequency spectrum of case 3 unsteady pressure

It is also apparent in Fig. 13 that the amplitude is consistently very near to 300 Pa throughout the entire domain for the case without the porous medium, whereas when the porous medium is included the amplitude varies slightly in the axial direction, downstream of the porous region. There is a variation of approximately 40 Pa in amplitude relative to the case without the porous medium, at twice the spatial wavelength of the source. In the next section, it shall be shown that this is due to an acoustic reflection caused by the modeled porous medium.

5.1.3 Eigenmode Analysis.

To explore the interaction of the porous medium and the acoustic source in further depth, we shall use the method of eigenmode analysis, which is an established and important method within turbomachinery aeroacoustics [5–7,27]. This method assumes a flow solution comprised of small perturbations in the primitive flow variables

Q_{P} = [ρ^{'}, u^{'}, p^{'}]

over a steady mean flow. The perturbation flow field is assumed to consist of radial eigenmodes of the form shown in Eq. . Each mode is characterized by a radial function

Q_{P} (r)

⁠, a temporal frequency

ω

⁠, a circumferential mode number

m

⁠, and an axial wave number

k

[27].

Q_{P} = e x p (i ω t + i m θ + i k x) Q_{P} (r)

(40)

A generalized eigenvalue problem (GEP) is then constructed using the linearised conservation equations in cylindrical coordinates, with radial flow gradients calculated numerically by a finite-difference method. The numerical solution of the constructed GEP then provides the complex axial wavenumbers $k = k_{r} + i k_{i}$ and the corresponding eigenvectors of modes with a specified frequency and mode number. For subsonic flow, the modes may be grouped into four categories: upstream propagating acoustic, downstream propagating acoustic, downstream propagating entropy, and downstream propagating vorticity. Most important to this study shall be the acoustic modes, which are characterized by large pressure content in their eigenvectors. The imaginary component of the axial wavenumber can then further classify the mode; a positive value (⁠ $k_{i} > 0$ ⁠) characterizes a downstream propagating evanescent (cut-off) mode, and a negative value (⁠ $k_{i} < 0$ ⁠) characterizes an upstream propagating evanescent mode. A zero value (⁠ $k_{i} = 0$ ⁠) characterizes a propagating (cut-on) mode. In practice, a small imaginary component may be added to the frequency input of the GEP (⁠ $ω_{i} \approx - 10^{- 5} ω_{r}$ ⁠), resulting in a small imaginary component to the wave number of the mode, thus revealing its direction of propagation.

We shall next apply this method to the numerical results of the unsteady porous-plug application case, using an in-house code that has been validated and applied extensively to industrial turbomachinery aeroelasticity problems [5–7]. This enables the detection and quantification of acoustic phenomena such as attenuation and reflection induced by the modeled porous medium, through analysis of the modal content as described above.

First, the modes were determined by inspection of the axial wave numbers. As expected, the upstream propagating acoustic mode introduced by the sinusoidal outlet boundary condition is clearly detected. Furthermore, there is a downstream propagating acoustic mode with significant pressure content, implying an acoustic reflection caused by the porous medium. For conciseness, we shall refer to the upstream and downstream propagating modes as upstream and downstream, respectively, for the remainder of this discussion. As also expected, the entropy and vorticity modes were also detected, along with evanescent acoustic modes of low pressure amplitude. Figure 15 shows the wave numbers of the upstream and downstream modes at an axial position near the outlet, along with the entropy and vorticity modes, and some of the evanescent acoustic modes.

Fig. 15

View large Download slide

Case 3 axial wave numbers (k) calculated at outlet

The upstream and downstream modes are present also in the equivalent duct without the porous medium, shown in Fig. 12(a), although the pressure magnitude of the downstream mode in this case is negligibly small. To confirm these observations, Table 2 quantitatively compares the wave numbers and the radial average of the pressure amplitude (⁠

| p |

⁠) of the upstream and downstream modes, at an axial position downstream of the porous medium. To calculate

| p |

⁠, consider that for a radial mode

n

with circumferential mode number

m

⁠, the pressure content may be calculated at axial position

x

and radial position

r

as in Eq. [27]. Here,

a_{m n}

is the modal amplitude, and the modal pressure value is given by the fifth element of the eigenvector (⁠

u_{5_{m n}}

⁠).

p_{m n} (x, r) = a_{m n} (x) u_{5_{m n}} (r)

(41)

Table 2

Axial wave numbers and radial average pressure amplitude of upstream and downstream propagating modes, calculated at $x = 1.5$ m

	Upstream	Downstream
k	Porous	25.7 − 0.00026i	−14.6 + 0.00015i
	Empty	25.7 − 0.00026i	−14.6 + 0.00015i
$\| p \|$ (Pa)	Porous	300.5	37.2
	Empty	300.5	0.00329

Considering that the mode (and so

p_{m n}

⁠) is complex-valued, the amplitude of

p_{m n}

is then averaged over the set of discrete radial positions

r = 1, 2, \dots, N_{r}

to give the radially averaged pressure amplitude of the mode. Since in this study the imposed sinusoidal pressure boundary condition is radially uniform, the expected

| p |

value is the amplitude of the planar sinusoidal wave.

| p | = \frac{1}{N_{r}} \sum_{r = 1}^{N_{r}} | p_{m n_{r}} |

(42)

Next, the eigenmode analysis was performed at varying axial positions throughout the domain, to investigate the axial variation of the pressure amplitude. Figure 16 shows that the average pressure amplitudes of both modes reduce significantly between the fluid regions upstream and downstream of the porous medium. The amplitude of the downstream mode reduces to a negligibly small value upstream of the porous medium, confirming that this mode represents the reflected acoustic wave. It should be noted that the analysis was not performed within the region of porous medium (⁠ $0.5 \leq x \leq 0.7$ ⁠), since this would violate the assumption of an axially uniform base flow required for the theoretical validity of the eigenmode analysis.

Fig. 16

View large Download slide

Radially averaged pressure amplitudes of acoustic upstream and downstream propagating modes for case 3

For future aeroelastic application, it is important to further investigate the influence of the porous media model on the amplitude of the reflected wave, since as explained in Sec. 1, many aeroelastic problems are known to be sensitive to acoustic reflections. Since the pressure amplitudes of the upstream and downstream modes may now be precisely quantified using the eigenmode analysis, it is possible to define the reflection coefficient

R

⁠, as shown in Eq. . For this, the established acoustic theory for the propagation of a planar wave across an interface between two media of different acoustic impedance is used; in this case, the incident wave propagates from the fluid region downstream of the porous medium, into the porous region, through the porous–fluid interface. The pressure amplitude of incident wave (⁠

P_{i}

⁠) is given by the upstream mode and that of the reflected wave (⁠

P_{r}

⁠) by the downstream mode.

R = \frac{P_{r}}{P_{i}} = \frac{| p |_{down}}{| p |_{up}}

(43)

First, the variation of the reflection coefficient with the porous media permeability coefficients is investigated, since as mentioned previously these are in practice highly dependent on factors such as the configuration of the heat exchanger core and flight conditions. A series of cases were studied in which the axial components of the viscous (⁠ $K_{D_{x x}}$ ⁠) and inertial (⁠ $K_{F_{x x}}$ ⁠) permeability coefficients were varied independently. The $y$ and $z$ components were fixed for this study, since for the considered configuration they are defined to be large enough to represent an impermeable boundary in these directions. Furthermore, this was a good opportunity to demonstrate the robustness of the presented porous media model and its implementation, since stable solutions were obtained across the range of variation in both coefficients $K_{D_{x x}}$ and $K_{F_{x x}}$ ⁠. Figure 17 shows that the reflection coefficient varies nonlinearly with either viscous or inertial permeability coefficient, away from the value of $R = 12.5 %$ for the standard permeability coefficient values in Table 1. As expected, a higher value of either permeability coefficient invokes a stronger acoustic reflection. By comparing Figs. 17(a) and 17(b), it is seen that the reflection coefficient is significantly more sensitive to the inertial resistance than to the viscous resistance. The reflection coefficient increases only slightly over an increase of several orders of magnitude in $K_{D_{x x}}$ ⁠, whereas it decreases to near-zero within the limit of $K_{F_{x x}} \to 0$ ⁠, showing that the reflected wave is predominantly caused by the inertial resistance of the porous medium.

Fig. 17

View large Download slide

Variation of reflection coefficient with permeability coefficients: (⁠ $a$ ⁠) viscous and (⁠ $b$ ⁠) inertial

Next, the variation in reflection coefficient with the frequency of the incident acoustic wave is considered, across a range of frequencies commonly of importance to aeroelastic problems. For the following cases, the permeability coefficients were maintained at the standard values shown in Table 1. Figure 18 shows that the reflection coefficient was found to decrease significantly with an increase in frequency. Furthermore, as the frequency is decreased from the standard value of 1000 Hz, the increase in reflection becomes highly nonlinear; it is expected that this is due to the acoustic wavelength increasing toward the scale of the domain length, whereas for the standard frequency and above, the wavelength is significantly shorter than the domain length.

Fig. 18

View large Download slide

Variation of reflection coefficient with frequency of incident wave

Lastly, the sensitivity of the reflection coefficient to variation in the amplitude of the incident wave was studied, in which cases where the amplitude was increased above the standard value of 300 Pa were considered. Again the permeability coefficients were maintained at the values in Table 1. Figure 19 shows that a decrease in the reflection coefficient was observed with an increase in amplitude; however, this decrease is so small (⁠ $< 0.5 %$ ⁠) that it should not be considered significant, and the reflection coefficient was found to be effectively constant across the range of amplitudes tested. It should be noted that since the reflection coefficient is the ratio of reflected to incident amplitudes, the amplitude of the reflected wave therefore increases linearly with the amplitude of the incident wave.

Fig. 19

View large Download slide

Variation of reflection coefficient with amplitude of incident wave

6 Conclusions

This study presents the current progress in the development of heat exchanger modeling capability in the aeroelastic solver au3d. The method of choice was to model the heat exchanger as porous medium, which is an established approach in CFD. A porous media model suitable for compressible and unsteady flows at local thermal equilibrium was presented, along with its implementation within the solver. The model was then verified using steady flow situations and representative cases, showing very close agreement with the commercial solver fluent. We then presented an application to a representative unsteady flow problem featuring a heat exchanger modeled as a porous medium, involving the interaction of acoustic waves and the modeled porous medium. An eigenmode analysis was performed to analyze the resultant acoustic field in terms of its modal content. The effects of varying permeability of the porous medium, as well as the frequency and amplitude of the acoustic source, were explored. Future work will involve the extension of the current capability to nonlocal thermal equilibrium models, thus progressing to models allowing for the inclusion of heat addition or rejection by the modeled heat exchanger. The pursuit of experimental reference data for the validation of unsteady heat exchanger flows should also be considered in the future.

Acknowledgment

The authors thank Rolls-Royce plc for supporting this research as part of the ACRE program.

Conflict of Interest

There are no conflicts of interest.

Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

Nomenclature

$k$ =: Thermal conductivity (W/m/K)
$p$ =: Static pressure (Pa)
$u$ =: Velocity (m/s)
$E$ =: Specific energy (J/kg)
$T$ =: Static temperature (K)
$K_{D}$ =: Darcy (viscous) permeability coefficient (m²)
$K_{F}$ =: Forchheimer (inertial) permeability coefficient (m)

Greek Symbols

$ε$ =: Porosity
$μ$ =: Dynamic viscosity (kg/m/s)
$ρ$ =: Density (kg/m³)
$σ$ =: Viscous shear stress tensor (N/m²)

Superscripts and Subscripts

s =: Solid phase
f =: Fluid phase
n =: Discrete time level
m =: Discrete pseudo-time level

Dimensionless Groups

M =: Mach number
Re =: Reynolds number

References

1.

Dietrich

,

P.

, and

Schneider

,

M.

,

2023

, “

Aeroacoustic Simulations of an Axial Fan With Modelled Turbulent Inflow Conditions

,”

Inter. J. Turbomach. Propulsion Power

,

8

(

2

), p.

13

.

Google Scholar

Crossref

2.

Lucius

,

A.

,

Schneider

,

M.

,

Schweitzer-De Bortoli

,

S.

,

Gerhard

,

T.

, and

Geyer

,

T.

,

2019

, “

Aeroacoustic Simulation and Experimental Validation of Sound Emission of an Axial Fan Applied in a Heat Pump

,”

Proceedings of the 23rd International Congress on Acoustics

,

Aachen, Germany

,

Sept. 9–13

, pp.

1851

–

1858

.

3.

Rabe

,

D. C.

,

Williams

,

C.

, and

Hah

,

C.

,

1999

, “

Inlet Flow Distortion and Unsteady Blade Response in a Transonic Axial-Compressor Rotor

,”

International Symposium on Air Breathing Engines

,

Florence, Italy

,

Sept. 5–10

,

ISABE 99-7287

.

4.

Pan

,

T.

,

Yan

,

Z.

,

Lu

,

H.

, and

Li

,

Q.

,

2021

, “

Numerical Investigation on the Forced Vibration Induced by the Low Engine Order Under Boundary Layer Ingestion Condition

,”

Aeros. Sci. Technol.

,

115

, p.

106802

.

5.

Zhao

,

F.

,

Nipkau

,

J.

, and

Vahdati

,

M.

,

2016

, “

Influence of Acoustic Reflections on Flutter Stability of An Embedded Blade Row

,”

J. Power Energy

,

230

(

1

), pp.

29

–

43

.

Google Scholar

Crossref

6.

Zhao

,

F.

,

Smith

,

N.

, and

Vahdati

,

M.

,

2017

, “

A Simple Model for Identifying the Flutter Bite of Fan Blades

,”

ASME J. Turbomach.

,

139

(

7

), p.

071003

.

Google Scholar

Crossref

7.

Lee

,

K.-B.

,

Wilson

,

M.

, and

Vahdati

,

M.

,

2017

, “

Numerical Study on Aeroelastic Instability for a Low-Speed Fan

,”

ASME J. Turbomach.

,

139

(

7

), p.

071004

.

Google Scholar

Crossref

8.

Liu

,

H.

,

Azarpeyvand

,

M.

,

Wei

,

J.

, and

Qu

,

Z.

,

2015

, “

Tandem Cylinder Aerodynamic Sound Control Using Porous Coating

,”

J. Sound Vib.

,

334

, pp.

190

–

201

.

Google Scholar

Crossref

9.

Sun

,

Y.

,

Wang

,

X.

,

Du

,

L.

, and

Sun

,

X.

,

2020

, “

Effect of Acoustic Treatment on Fan Flutter Stability

,”

J. Fluids Struct.

,

93

, p.

102877

.

10.

Varvill

,

R.

,

Duran

,

I.

,

Kirk

,

A.

,

Langridge

,

S.

,

Nailard

,

O.

,

Payne

,

R.

, and

Webber

,

H.

,

2019

, “

Sabre Technology Development: Status and Update

,”

8th European Conference for Aeronautics and Space Sciences (EUCASS)

,

Madrid, Spain

,

July 1–4

,

EUCASS 2019-307

.

11.

Ding

,

W.

,

Eri

,

Q.

,

Kong

,

B.

, and

Wang

,

C.

,

2022

, “

The Effect of Inlet Pressure Distortion on the Performance of an Axisymmetric Compact Tube Heat Exchanger With Radial Counter Flow Type for Hypersonic Pre-Cooled Aero-Engine

,”

J. Mech. Sci. Technol.

,

36

(

6

), pp.

3181

–

3191

.

Google Scholar

Crossref

12.

Missirlis

,

D.

,

Yakinthos

,

K.

,

Palikaras

,

A.

,

Katheder

,

K.

, and

Goulas

,

A.

,

2005

, “

Experimental and Numerical Investigation of the Flow Field Through a Heat Exchanger for Aero-Engine Applications

,” pp.

440

–

458

.

13.

Missirlis

,

D.

,

Donnerhack

,

S.

,

Seite

,

O.

,

Albanakis

,

C.

,

Sideridis

,

A.

,

Yakinthos

,

K.

, and

Goulas

,

A.

,

2010

, “

Numerical Development of a Heat Transfer and Pressure Drop Porosity Model for a Heat Exchanger for Aero Engine Applications

,”

Appl. Therm. Eng.

,

30

(

11–12

), pp.

1341

–

1350

.

Google Scholar

Crossref

14.

Alshare

,

A. A.

,

Simon

,

T. W.

, and

Strykowski

,

P. J.

,

2010

, “

Simulations of Flow and Heat Transfer in a Serpentine Heat Exchanger Having Dispersed Resistance With Porous-Continuum and Continuum Models

,”

Int. J. Heat Mass Transfer

,

53

(

5-6

), pp.

1088

–

1099

.

Google Scholar

Crossref

15.

Wang

,

W.

,

Guo

,

J.

,

Zhang

,

S.

, and

Yang

,

J.

,

2014

, “

Numerical Study on Hydrodynamic Characteristics of Plate-Fin Heat Exchanger Using Porous Media Approach

,”

Comput. Chem. Eng.

,

61

, pp.

30

–

37

.

Google Scholar

Crossref

16.

Yang

,

J.

,

Ma

,

L.

,

Bock

,

J.

,

Jacobi

,

A. M.

, and

Liu

,

W.

,

2014

, “

A Comparison of Four Numerical Modeling Approaches for Enhanced Shell-and-Tube Heat Exchangers with Experimental Validation

,”

Appl. Therm. Eng.

,

65

(

1–2

), pp.

369

–

383

.

Google Scholar

Crossref

17.

Musto

,

M.

,

Bianco

,

N.

,

Rotondo

,

G.

,

Toscano

,

F.

, and

Pezzella

,

G.

,

2016

, “

A Simplified Methodology to Simulate a Heat Exchanger in An Aircraft’s Oil Cooler by Means of a Porous Media Model

,”

Appl. Therm. Eng.

,

94

, pp.

836

–

845

.

Google Scholar

Crossref

18.

An

,

C. S.

, and

Kim

,

M. H.

,

2018

, “

Thermo-Hydraulic Analysis of Multi-Row Cross-Flow Heat Exchangers

,”

Int. J. Heat Mass Transfer

,

120

, pp.

534

–

539

.

Google Scholar

Crossref

19.

Jarauta

,

A.

,

Zingan

,

V.

,

Minev

,

P.

, and

Secanell

,

M.

,

2020

, “

A Compressible Fluid Flow Model Coupling Channel and Porous Media Flows and Its Application to Fuel Cell Materials

,”

Transp Porous Med

,

134

, pp.

351

–

386

.

Google Scholar

Crossref

20.

Prokein

,

D.

,

Dittert

,

C.

,

Böhrk

,

H.

, and

von Wolfersdorf

,

J.

,

2020

, “

Numerical Simulation of Transpiration Cooling Experiments in Supersonic Flow Using OpenFOAM

,”

CEAS Space J.

,

12

, pp.

247

–

265

.

Google Scholar

Crossref

21.

Whitaker

,

S.

,

1966

, “

The Equations of Motion in Porous Media

,”

Chem. Eng. Sci.

,

21

(

3

), pp.

291

–

300

.

Google Scholar

Crossref

22.

Sayma

,

A. I.

,

Vahdati

,

M.

,

Sbardella

,

L.

, and

Imregun

,

M.

,

2000

, “

Modeling of Three-Dimensional Viscous Compressible Turbomachinery Flows Using Unstructured Hybrid Grids

,”

AIAA J.

,

38

(

6

), pp.

945

–

954

.

Google Scholar

Crossref

23.

Stapelfeldt

,

S.

, and

Vahdati

,

M.

,

2019

, “

Improving the Flutter Margin of an Unstable Fan Blade

,”

ASME J. Turbomach.

,

141

(

7

), p.

071006

.

Google Scholar

Crossref

24.

Yang

,

M.

, and

Sheng

,

P.

,

2017

, “

Sound Absorption Structures: From Porous Media to Acoustic Metamaterials

,”

Annu. Rev. Mater. Res.

,

47

, pp.

83

–

114

.

Google Scholar

Crossref

25.

Horoshenkov

,

K. V.

,

Hurrell

,

A.

, and

Groby

,

J. -P.

,

2019

, “

A Three-Parameter Analytical Model for the Acoustical Properties of Porous Media

,”

J. Acous. Soc. Am.

,

145

(

4

), pp.

2512

–

2517

.

Google Scholar

Crossref

26.

Satcunanathan

,

S.

,

Zamponi

,

R.

,

Meinke

,

M.

,

Van de Wyer

,

N.

,

Schram

,

C.

, and

Schröder

,

W.

,

2019

, “

Validation for a Model for Acoustic Absorption in Porous Media

,”

Proceedings of the 48th International Congress and Exhibition on Noise Control Engineering

,

Madrid, Spain

,

June 16–19

.

27.

Moinier

,

P.

, and

Giles

,

M. B.

,

2005

, “

Eigenmode Analysis for Turbomachinery Applications

,”

J. Propul. Power

,

21

(

6

), pp.

973

–

978

.

Google Scholar

Crossref

28.

Whitaker

,

S.

,

1986

, “

Flow in Porous Media I : A Theoretical Derivation of Darcy’s Law

,”

Trans. Porous Media

,

1

, pp.

3

–

25

.

Google Scholar

Crossref

29.

Whitaker

,

S.

,

1996

, “

The Forchheimer Equation: A Theoretical Development

,”

Trans. Porous Media

,

25

, pp.

27

–

61

.

Google Scholar

Crossref

30.

Whitaker

,

S.

,

1986

, “

Local Thermal Equilibrium: An Application to Packed Bed Catalytic Reactor Design

,”

Chem. Eng. Sci.

,

41

(

8

), pp.

2029

–

2039

.

Google Scholar

Crossref

31.

Dubuc

,

L.

,

Cantariti

,

F.

,

Woodgate

,

M.

,

Gribben

,

B.

,

Badcock

,

K. J.

, and

Richards

,

B. E.

,

1998

, “

Solution of the Unsteady Euler Equations Using an Implicit Dual-Time Method

,”

AIAA J.

,

36

(

8

), pp.

1417

–

1424

.

Google Scholar

Crossref

32.

Jameson

,

A.

,

1991

, “

Time Dependent Calculations Using Multigrid, With Applications to Unsteady Flows Past Airfoils and Wings

,”

10th Computational Fluid Dynamics Conference

,

Honolulu, HI

,

June 24–26

.

Google Scholar

Crossref

Development of Heat Exchanger Modeling Capability for a Finite-Volume Aeroelasticity Solver

Abstract

1 Introduction

2 Flow Model

2.1 Volume-Averaging Theory.

2.2 Application to the Conservation Equations.

2.2.1 Mass.

2.2.2 Momentum.

2.2.3 Energy.

2.3 Boundary Conditions

2.3.1 Porous–Fluid Interfaces.

2.3.2 Solid Walls.

3 Numerical Implementation

3.1 Temporal Discretization.

4 Results for Steady Flows

4.1 Case 1: Duct With Porous Plug.

4.1.1 Numerical Setup.

4.1.2 Analysis of Steady Flow-Field.

4.2 Case 2: Duct With Porous Cube.

4.2.1 Numerical Setup.

4.2.2 Analysis of Steady Flow Field.

5 Results for Unsteady Flows

5.1 Case 3: Unsteady Porous Plug

5.1.1 Numerical Setup.

5.1.2 Spatial Variation of Unsteady Pressure.

5.1.3 Eigenmode Analysis.

6 Conclusions

Acknowledgment

Conflict of Interest

Data Availability Statement

Nomenclature

Greek Symbols

Superscripts and Subscripts

Dimensionless Groups

References

Get Email Alerts

Cited By

Related Articles

Related Proceedings Papers

Related Chapters

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

This Feature Is Available To Subscribers Only