Analysis of non-linear losses in a parallel plate thermoacoustic stack

Abstract

Purpose

This study aims to analyse the non-linear losses of a porous media (stack) composed by parallel plates and inserted in a resonator tube in oscillatory flows by proposing numerical correlations between pressure gradient and velocity.

Design/methodology/approach

The numerical correlations origin from computational fluid dynamics simulations, conducted at the microscopic scale, in which three fluid channels representing the porous media are taken into account. More specifically, for a specific frequency and stack porosity, the oscillating pressure input is varied, and the velocity and the pressure-drop are post-processed in the frequency domain (Fast Fourier Transform analysis).

Findings

It emerges that the viscous component of pressure drop follows a quadratic trend with respect to velocity inside the stack, while the inertial component is linear also at high-velocity regimes. Furthermore, the non-linear coefficient b of the correlation ax + bx² (related to the Forchheimer coefficient) is discovered to be dependent on frequency. The largest value of the b is found at low frequencies as the fluid particle displacement is comparable to the stack length. Furthermore, the lower the porosity the higher the Forchheimer term because the velocity gradients at the stack geometrical discontinuities are more pronounced.

Originality/value

The main novelty of this work is that, for the first time, non-linear losses of a parallel plate stack are investigated from a macroscopic point of view and summarised into a non-linear correlation, similar to the steady-state and well-known Darcy–Forchheimer law. The main difference is that it considers the frequency dependence of both Darcy and Forchheimer terms. The results can be used to enhance the analysis and design of thermoacoustic devices, which use the kind of stacks studied in the present work.

Citation

Di Meglio, A., Massarotti, N., Rolland, S. and Nithiarasu, P. (2024), "Analysis of non-linear losses in a parallel plate thermoacoustic stack", International Journal of Numerical Methods for Heat & Fluid Flow, Vol. 34 No. 1, pp. 353-377. https://doi.org/10.1108/HFF-02-2023-0071

NOMENCLATURE

Symbol

Meaning Unit

a

Linear coefficient of the correlation av + bv² Pa · s/m²

b

Non-linear coefficient of the correlation av + bv² kg/m⁴

c_F

Forchheimer coefficient −

C_VM

Virtual Mass coefficient −

d

Solid half-plate thickness or total derivative m

e

Internal energy J/m³

f₀

Fundamental frequency of pressure input Hz

f_v

Average Viscous function −

h

Viscous function −

i

Imaginary unity −

k_f

Fluid thermal conductivity W/m/K

L_stack

Stack length m

L_tot

Total length of the computational domain m

p

Pressure Pa

t

Time s

T

Temperature K

v

Velocity vector m/s

x_s

Stack position m

y₀

Half of the fluid channel width m

Greek letters

α

Thermal diffusion m²/s

δ_v

Viscous penetration depth m

∂

Partial derivative

θ

Phase shift between pressure gradient and velocity −

κ

Permeability m²

κ ˜c

Complex permeability m²

µ

Dynamic viscosity Pa · s

v

Kinematic viscosity m²/s

ξ

Fluid particle displacement m

π

Dimensionless numbers group −

ρ

Density kg/m³

τ

Time scale coefficient s

ϕ

Porosity−

ω₀

Angular frequency Hz

∇(·)

Gradient (Divergence)

Subscript

i

Intrinsic

s

Superficial

VM

Virtual Mass

KE

Kinetic Energy

1

First order variable

0

Mean-time variable

Abbreviation

CFD

Computational Fluid Dynamics

DF

Darcy-Forchheimer

REV

Representative Elementary Volume

SPL

Sound Pressure Level

1. Introduction

Oscillating flows in porous media are encountered in a broad range of real-life applications including medicine (Keith Sharp et al., 2019), geophysics (Zhang et al., 2019), acoustics (Horoshenkov, 2017) and thermoacoustics (Chen et al., 2021). In general, heat and fluid flow in porous media can be studied either at the microscopic level, by directly solving the Navier Stokes Equations (NSE) with an appropriate numerical model, by using analytical solutions, or at the macroscopic scale (Massarotti et al., 2003). In the latter case, to reduce the prohibitive computational costs, the NSE are averaged over a Representative Elementary Volume (REV) and the additional variables require for a closure model. The Darcy–Forchheimer model (Whitaker, 1986; Whitaker, 1996) is one of the most used closure models of the momentum equation in the available literature. It consists of two source terms: a viscous linear term (Darcy term, in which the permeability of the medium has to be specified) and a non-linear form drag (Forchheimer term). Examples, in which the porous media approach has been adopted to derive analytical or powerful numerical solutions, can be found in recent refs (Bhatti and Powell, 2023; Bhatti and Öztop, 2022; Wakif et al., 2022; Rasool et al., 2023; Wakif, 2023; Elboughdiri et al., 2023a; Sharma et al., 2023; Zhang et al., 2023; Elboughdiri et al., 2023b), thus demonstrating the importance of porous media in multiphysics problems such as nano fluids, biodynamics and heat/mass transfer applications.

Such a model can only be applied when frequency tends to zero, as the oscillating inertial contributions to the flow can be neglected. As frequency increases, the inertia-induced phase shift between pressure gradient and velocity has to be taken into account. Such a phase shift can be intrinsically considered in the frequency domain with a complex permeability, introduced by Johnson et al. (1987). In the time domain, several attempts have been made to adapt the Darcy steady-state model to oscillating flows. Within this framework, three approaches can be found in the literature. The first approach is the so-called “VANS” (volume average Navier Stokes), valid only for low frequencies and uses a time scale coefficient in the unsteady term obtained by scaling the steady-state Darcy model (Zhu and Manhart, 2016). The second, named virtual mass approach, extends the previous model by considering an additional term dependent on porosity to consider the phase shift between pressure gradient and velocity (Lowe et al., 2008). However, this coefficient is required to be determined experimentally to close the model. Zhu et al. (2014) proposed a numerical calculation, via direct numerical simulations, to determine the virtual mass coefficient by scaling a kinetic energy equation. The third method is the one proposed by Di Meglio et al. (2022) to analytically close this problem, without performing prior simulations at the microscopic level, by linking the real and imaginary part of Johnson’s complex permeability to the classical steady-state source term and the virtual mass correction, respectively.

Similar systematic studies of the Forchheimer term’s frequency dependence have not been conducted, mainly due to the lack of analytical or semi-analytical solutions. However, some attempts were made in the past, in acoustics, to consider non-linear propagation losses in porous media. For example, Melling (Melling, 1973) theoretically characterised the non-linear trend of the impedance (as the ratio between pressure gradient and velocity) of a perforated plate stack by using stationary data. An analogous approach was used to evaluate the effect of sudden contraction-expansions in oscillating flows (Clark et al., 2007). Both cases, based on a steady-state parametrization, did not consider the effect of frequency. Wilson et al. (1988) instead introduced a Forchheimer-like correction factor for the complex density and compared the numerical solutions of the non-linear wave equations with available experimental data for porous media (foams) in terms of the acoustic impendence against the Sound Pressure Level (SPL). In this acoustic framework, Wang et al. (2009) successfully adopted a semi-empirical Forchheimer coefficient to better predict the sound absorption coefficient against the SPL. Ge et al. (2012) carried out experiments on the non-linear impedance of a mesh screen stack discovering that the resistance part of the impedance (directly related to permeability and Forchheimer term) follows a non-linear trend while the reactance (imaginary part of the impedance, directly related to the above-mentioned virtual mass coefficient) stays linear also at high SPLs. Therefore, three parameters were used to fully describe porous media flow in this case: static permeability, virtual mass correction and Forchheimer term. Furthermore, they were found to be independent of frequency, but dependent on porosity and stack length (due to tortuosity). In the works of Liu (2005), the viscous and inertial pressure gradients were obtained experimentally for a uniform cross-section stacks (square channels) and a mesh screen regenerator. In both cases, the viscous pressure gradient is well-fitted by the Darcy–Forchheimer law. In particular, the mesh-screen regenerator experienced higher non-linear losses than the squared channels stack due to tortuous geometry. For a simple tortuous geometry, such as a cylinder bank (or transversal pin array stack) in a previous computational fluid dynamics (CFD) work (Di Meglio and Massarotti, 2022a), it has been confirmed that non-linear losses are not negligible when Reynolds number (based on the pin diameter and maximum velocity) is higher than 1.

A significant knowledge gap exists in the available literature on an in-depth characterisation of the Darcy–Forchheimer correlation in oscillatory flows, even for simple porous media geometry such as the parallel plate stacks. Oscillatory flows through parallel plate stacks have been only studied at the microscopic level. For instance, Jaworski et al. (2009) carried out an experimental campaign to study the entrance effects, as well as the vortex-shedding phenomena coming from the discontinuity between the free fluid region and the stack zone (Shi et al., 2010). In the latter article, a classification of the fluid flow microstructures using three dimensionless numbers was proposed (Jaworski et al., 2009). Moreover, Shi et al. (2011) numerically calculated a correlation between the Strouhal and Reynolds numbers. However, none of the above works express the non-linear Forchheimer losses synthetically. The need of a DF-based macroscopic model is important for rapidly designing systems and reducing the computational costs. An important practical application where these flows occur is thermoacoustics. Here, the porous media is generally called stack, and it is responsible for energy conversion between heat and acoustic power (and vice-versa). The parallel plate stack is among the most efficient stacks, as testified by the numerous articles in which this stack is used (Marx and Blanc-Benon, 2004; Liu et al., 2019; Ke et al., 2010). To date, there are two main numerical tools available for studying thermoacoustics: DeltaEC and the finite volume-based CFD. The first one is one-dimensional tool useful for analysis/design in the frequency domain. Here, the Forchheimer term is not included but only a correction factor is implemented through stationary data fitting (Clark et al., 2007). Among CFD-based models, as highlighted in the results of a recent review article (Di Meglio and Massarotti, 2022b), no case study is available that uses the Forchheimer term for this kind of stacks. In this framework, the difference between the microscopic and macroscopic models (with the only Darcy/linear term) in terms of pressure amplitude is not negligible (Di Meglio and Massarotti, 2023).

Based on the above observations, it can be concluded that a detailed study of the oscillatory flows using a macroscopic analysis is not available. Thus, the present study aims to:

• Develop an appropriate CFD model, including the stack and the surrounding thermoacoustic device in which it is operated. This is an essential requirement to fully characterise the non-linear behaviour of the porous media.

• Develop an appropriate post-processing procedure to analyse the viscous losses (where velocity is phased in time with the pressure gradient across the stack) separately from the oscillating inertial effects. This is done by analysing the results in the frequency domain.

• Develop a numerical correlation between the pressure gradient (both components, that one in phase and the other one delayed by 90° with velocity) and velocity to obtain the Darcy–Forchheimer relationship in oscillatory flows. This is done by varying the pressure amplitude. (The correlation are obtained using the CFTool library in MATLAB.)

• Study the effects of varying the frequency of the boundary condition on the Forchheimer term. This is done to have a clearer picture of the Forchheimer term similar to what was done in the past for the Darcy term.

• Study the effect of the stack porosity on the Forchheimer coefficient. This is done to identify the major sources of non-linearities in such stacks.

The results of this work show that the component of pressure drop in phase with velocity, follows a quadratic trend with respect to velocity, as expected for a similar steady-state non-linear porous media model. On the other hand, the inertial pressure gradient, the component of the pressure drop 90° out of phase with velocity, continues to be linear also for higher velocities in the stack. Furthermore, the dependence of the Forchheimer coefficient on the frequency and porosity is highlighted for the first time. Finally, the analysis of local correlations (inside and outside the stack) between pressure gradient and velocity confirms that the main sources of nonlinearities are the geometrical discontinuities of the section area at the ends of stacks and the finite length of the stack itself.

The subsequent sections of the article are organised as follows. In Section 2, the governing equations both at the microscopic and macroscopic scales are presented by reviewing the microscopic solution of the parallel plate stack in the frequency domain as well as the concept of complex permeability. In Section 3, the numerical setup and the postprocessing procedure of the data are explained, whereas in Section 4, the results are presented and discussed. Finally, conclusions are drawn in Section 5.

2. Governing equations

2.1 Microscopic scale

At the microscopic scale, oscillating flows are governed by the unsteady and compressible NSEs. They are reported here below, for a Newtonian fluid, in the conservative form, suitable to be solved with a finite volume approach implemented in OpenFOAM (Greenshields, 2015). The scope of application considered in this work is in the field of thermoacoustic stack/heat exchangers. Their Mach number based on acoustic peak velocity in the simulations is always less than 0.1 and compressibility effects could be neglected in the momentum equation (2) as well as the viscous dissipation in energy equation [37]. However, in such kinds of problems, the propagation of disturbances at a finite speed of sound requires considering compressibility effects both in the mass balance equation (1) (∂ρ/∂t) as well as the time derivative of pressure in the energy equation (3) (∂ρ/∂t). The energy equation is written in terms of sensible enthalpy h (Greenshields, 2015). Pressure p, density ρ and temperature T are related by the ideal gas law:

(1) ∂ρ∂t+∇⋅(ρv)=0

(2) ∂(ρv)∂t+∇⋅(ρvv)=−∇p+μ∇2v+∇⋅[μ(∇v+∇Tv)]−23∇(μ(∇⋅v))

(3) ∂(ρh)∂t+∇⋅(ρvh)+∂(ρK)∂t+∇⋅(ρvK)−∂p∂t=∇⋅(α∇h)

where v, µ, α are the velocity vector, dynamic viscosity and thermal diffusivity of the medium, respectively.

2.2 Macroscopic scale

To reduce computational costs, a general domain composed of both free fluid and porous media zones is introduced. This requires an additional pressure gradient source term in the momentum equation to simulate the forces exchanged microscopically between the fluid and solid matrix. Once a homogenous and periodic REV is identified in steady-state flows, the Darcy–Forchheimer law can be derived by the volumetric average of the momentum equation (VANS approach), as shown in (Vafai, 2005):

(4) ∇〈pi〉=−(μκ+cFκρϕ|〈v〉s|)〈v〉s

where κ is the permeability [m²], c_F the dimensionless Forchheimer coefficient, ϕ the porosity of the porous medium. The notation 〈〉_i,s refers to intrinsic (i) and superficial (s) averaged quantities, respectively.

In oscillating flows, the vast majority of theoretical studies focus on the linear term [first term in equation (4)]. Most of them aimed at extending the Darcy approximation of equation (4), by introducing a time scale coefficient τ in the unsteady term, as follows (Zhu and Manhart, 2016):

(5) −κμ∇〈p〉i=(1+κρϕcF|〈v〉s|)〈v〉s+τ∂〈v〉s∂t

According to the VANS approach (Zhu and Manhart, 2016), such a time scale coefficient can be derived by a scaling analysis, and it is equal to:

(6) τVANS=κνϕ

This approach is consistent with the steady-state solution but fails at high-frequency. In the virtual mass approach (Lowe et al., 2008), a correction factor to the previous formula is added to take into account the unsteady effects:

(7) τVM=κνϕ(1+CVM1−ϕϕ)

However, the virtual mass constant, C_VM, has to be defined experimentally depending on the frequency. Zhu et al. (2014) and Zhu and Manhart (2016) proposed a numerical calculation of the time coefficient τ based on the scaling of the kinetic energy equation:

(8) τKE=κνϕ〈v⋅v〉i〈v〉i〈v〉i

This approach has demonstrated the best agreement with analytical solutions, but it requires prior numerical simulations at the pore scale to be practically implemented (Bhatti and Öztop, 2022). According to the latest extension to the method, presented by Di Meglio et al. (2022), the concept of the complex permeability κ˜c, in terms of its real and imaginary part, can be useful to separate the viscous effect, in which the pressure gradient is time-phased with velocity, and the oscillating inertial effects in which the pressure gradient is 90° out of phase with velocity by definition. The time-domain linear macroscopic governing equation, limited to the porous media zone, is the following, as demonstrated in a previous work (Di Meglio et al., 2022):

(9) 1ϕρ ∂〈v〉s∂t{ 1+12ϕδv2⋅J[1κ˜c]} =−∇〈p〉i−μ⋅ℜ[1κ˜c]〈v〉s

where δv=2ν/ω is the viscous penetration depth (v is the kinematic viscosity) and ℜ(ℑ) represents the real (imaginary) part operator. The complex permeability for a parallel plate stack is (Di Meglio et al., 2022):

(10) κ˜c=ϕδv21−fv2i fv,fv=tanh⁡[(1+i)y0δv](1+i)y0δv

where y₀ is the half-channel height. The formula expressed in the equation (10), in terms of both real and imaginary parts, will be used as an analytical reference to verify the linear trend between pressure gradient and velocity in the stack.

3. Proposed method

3.1 Base reference: non-linear macroscopic model

Inspired by the non-linear Darcy–Forchheimer model for steady state [equation (4)], where the non-linear term is proportional to |〈v〉_s|〈v〉_s, and by the theoretical approach used in by Wilson et al. (1988), the general relation between macroscopic pressure gradient and superficial velocity in oscillating flows can be written as follows (Wilson et al., 1988):

(11) −∇〈p〉i=(aR+bRLstack|〈v〉s|)〈v〉s+1ϕρ ∂〈v〉s∂t{ 1+aV+bVLstack|〈v〉s|}

Equation (11) constitutes the analytical reference for the analysis carried out in this work. In the cases of negligible non-linear contributions, the link between equation (11) and the linear macroscopic momentum equation (9) is uniquely determined. The coefficients (a_R, a_V) are confirmed to be consistent with the formulation and analytical value presented in equations (9) and (10), whereas the other coefficients (b_R, b_V) are investigated more deeply for the first time (subscripts R and V stand for real and virtual, in analogy with the above-explained concept of virtual mass). The choice of using the stack length L_stack, that will be discussed more deeply in Sections 3 and 4, allows one to distinguish the linear viscous losses, independent of the stack length, from the non-linear losses. As the latter are determined by the stack discontinuities, they become more significant for smaller stacks. In the frequency domain, equation (11) for each frequency can be split into two different equations: the viscous pressure gradient in phase velocity [the first term of equation. (11)], and the inertial pressure gradient, 90° out of phase with velocity [second part of equation (11)].

3.2 Problem setup

3.2.1 Geometry and mesh.

The computational domain consists of a 2D pipe of length L_tot = 0.5 m. At a horizontal distance x_s = 0.1 m the stack (L_stack = 30 mm), composed of three fluid channels whose width is 2 y₀ = 1 mm and equally spaced by a gap depending on the porosity, (ϕ = 0.5, 0.6, 0.7, 0.8) is placed, as shown in Figure 1 (Jaworski et al., 2009; Chen et al., 2020). The computational domain and mesh are pictured in Figure 1. The mesh is made of three sub-meshes: a quadrilateral mesh inside the stack, a triangular mesh in the “transition zones” just outside the stack ends and a mapped mesh far away from these regions. The size of the smallest elements in the stack and transitions zones has been chosen according to the thinnest viscous penetration depth δ_v at 350 Hz.

3.2.2 Boundary and initial conditions.

As a boundary condition, a sinusoidal input pressure has been implemented on the left-hand side of the computational domain as follows:

(12) p(t)=p1sin⁡(ω0t)+p0

where p₁ and ω₀ = 2πf₀ are the pressure input and angular frequency and whose values are specified in Table 1 for all parametric simulations carried out in this work. The reference pressure p₀ has been fixed to 10⁵ Pa. The opposite side, on the right-hand side, adiabatic wall boundary conditions are applied. A no slip and isothermal at 300 K boundary condition is considered on the stack walls. Cyclic boundary conditions (periodic) are constrained on the upper and lower surfaces of the domain to mimic the transversal periodicity of the system, as only three solid plates and fluid channels are simulated in this work. A previous literature work demonstrated that three plates are the best compromise between accuracy and computational efforts of the simulations (Jaworski et al., 2009; Chen et al., 2020). Regarding the initial conditions, zero velocity, pressure and uniform temperature are used. A coarser mesh is used initially until the volumetric average velocity over the stack and the pressure drop calculated become periodic. Then, these results are used as initial conditions for the simulations using the finest mesh.

3.2.3 Numerical solver.

The governing equations presented in the previous section are solved within the open-source environment OpenFOAM, by using the so-called application “buoyantPimpleFoam” an unsteady and compressible thermal and fluid flow solver. A second-order discretization scheme has been used for the convection and diffusion terms, whereas a first-order implicit scheme has been used for time discretization (Table 2).

The right-hand side boundary has been assumed to be a solid wall. The stack walls are assumed to be isothermal so that the energy equation in the solid stack plate can be neglected (Di Meglio et al., 2021; Swift, 1988). Finally, periodic boundary conditions have been applied at the top and bottom boundary surfaces.

No turbulence model has been implemented because the highest Reynolds number based on the viscous penetration depth is always less than 500. This is a well-known and accepted critical Reynolds number value considered in the available literature for oscillating flows (Ohmi and Iguchi, 1982) when y₀ > δ_v (high-frequency limit). For y₀ < δ_v, the reference critical Reynolds number is about 2,300 at steady state (Swift, 2003). On the other hand, other experimental results show that some transitional turbulent effect arises even at lower Reynolds (Jaworski et al., 2009). In this work, a laminar model has been considered.

3.3 Post-processing procedure

Velocity and pressure are sampled at every time step of Δt = 10⁻⁵ s and for at least 10 acoustic periods to guarantee that the maximum frequency capturable is higher than the input frequency f₀ and a good frequency resolution. The post-processing sampling intervals do not correspond to the time steps used in the simulation, which are lower and adaptive such that the Courant number is less than 0.1. The numerical data are post-processed as follows. At the global level, the surface average of the pressure at the left and right ends of the stack are extracted (Figure 1), as well as the volumetric average of the velocity over the stack. The Fast Fourier Transform of these quantities are performed to isolate the fundamental frequency ω₀, used in the inlet boundary condition. The pressure gradient is calculated as the difference between the right and left surface average pressure divided by the total length of the stack. The velocity can be assumed to be at zero-phase without harming the generality of the procedure, because only the relative phase shift between pressure gradient and velocity, indicated as θ in Figure 2, matters. Therefore, the pressure drop can be decomposed along the velocity but opposite direction, that represents viscous pressure drop and in a vertical component, as illustrated in Figure 2.

From a mathematical perspective:

(13) {Δpviscous=|Δp|⋅cos⁡(θ)Δpinertial=|Δp|⋅sin⁡(θ)

For each frequency and porosity, both viscous and inertial pressure drops are correlated to the velocity at the fundamental frequency in the “CFTool” package provided by MATLAB to search the coefficients a and b of the following analytical expression:

(14) ∇pviscous,inertial|x=aR,Vv1+bR,VLstackv12

3.4 Dimensionless analysis

Performing a dimensionless analysis of the problem discussed in this study becomes essential in order to minimise the variables influencing the pressure drop. This, in turn, reduces the number of numerical simulations needed to accurately replicate the phenomena.

The pressure gradient can be written as a function of the following dimensional variables:

(15) ∇p=∇p (v,ν,y0, d, Lstack,ω)

where d is the half-thickness of the plate between stack passages. According to the Buckingham theorem (Misic et al., 2010), the dimensionless version of the equation (15) can be expressed as a function of the number of the independent variables (6) minus the number of the fundamental units involved in the specific problem (length, mass, time):

(16) ∇p*=∇p* (π1,π2,π3)

where π_i are the three dimensionless numbers.

The first dimensionless quantity is the Reynolds number, based on the superficial velocity as velocity scale and y₀ as length scale. The second dimensionless number must consider the frequency ω and several choices have been adopted previously in the available literature such as the ratio y₀/δ_v, or the fluid particle displacement |ξ| = |v_s|/(ω₀L_stack). However, the first one does not allow considering the effect of the stack length. The second instead does not fit the post-processing procedure due to the simultaneous dependence on the velocity and frequency. To overcome these problems, the ratio between ξ and Re has been chosen. The last dimensionless parameter is the porosity, as a combination of y₀ and d:

(17) {π1:Re=vy0νπ2:ξRe=2(δvy0)2y0Lstackπ3:ϕ=y0y0+d

The pressure gradient is made non-dimensional using the factor 12ρ0u2/y0  .

4. Results and discussion

4.1 Verification and validation

The reliability of the numerical model is assessed via mesh sensitivity analysis for verification, then via comparison with experimental data for the validation phase. The mesh sensitivity analysis is performed using relevant physical quantities for the present work, the volumetric velocity over the stack and the viscous pressure drop Δp_viscous, derived as explained in the previous Section 3.3, for the highest frequency (350 Hz) and Reynolds number in the stack (corresponding to p₁ = 10,000 Pa) investigated in this work. Six different meshes have been used. Results illustrate both parameters are independent of the mesh size when the number of finite volumes is approximately 150,000, as shown in Figure 3.

Furthermore, a validation step has been carried out to ensure that numerical models can reproduce the reality by comparing the CFD results with experimental data taken from the available literature for oscillating flows. More specifically, the experimental setup described by Berson et al. (2008) has been reproduced numerically, being very similar the numerical setup built in this article. For two different input pressures (500 Pa, 2,000 Pa), the instantaneous numerical velocity profiles along the section at x = 0.2227 m are compared to the experimental data for eight frames within an acoustic period. Figure 4 shows the comparison of velocity profiles along with experimental error bars. As seen, a good agreement between experimental and simulation data is observed, especially at the middle of the channel. Minor deviations close to the walls could be attributed to asymmetry induced by experimental settings.

4.2 Viscous and inertial pressure drop against velocity

The viscous and inertial pressure drop, calculated according to the equation (13) and evaluated at the fundamental frequency (350 Hz), are plotted against the volumetric average velocity in Figure 5. The solid line shown in Figure 5 represents the numerical fit characterised by an R² value above 0.9992. In addition to this, the prediction bounds (with a confidence interval of 95%) are also included to assess the quality of the curve fitting.

Overall, the viscous pressure drop follows a parabolic trend. However, the non-linear contribution is negligible as the velocity approaches zero. In Figure 5, the slope of the solid line represents the real part of the inverse of the complex permeability [equation (10)]. As seen in Figure 5(b), the inertial pressure drop is linear also at the highest velocities. Here, the slope of the curve is the imaginary part of the inverse of the complex permeability. This first finding is in agreement with what has been found in previous work focusing on a transversal pin array stack (Di Meglio and Massarotti, 2022a) or in other experimental works (Liu, 2005). As a result, in the equation (11) the virtual mass coefficient comes from only the linear contribution (b_V ≈ 0), while the non-linearities contributes to the viscous component of pressure drop (b_R ≠ 0).

(18) vs→0, aR→μ⋅ℜ[1κ˜c];aV→μ⋅ℑ[1κ˜c]

This result is valid for all frequencies investigated. The curves fitted by the numerical data in terms of dimensionless pressure gradient (named friction factor) against the Reynolds number are shown in Figure 6. At low Reynolds numbers, the higher the frequency, the larger the friction factor. This is valid up to a Reynolds number approximately equal to 300, where the trend is inverted. The lowest frequency curve is characterised by the highest friction factor value as the non-linear contribution b starts to have an influence.

4.3 Permeability and Forchheimer coefficients versus frequency

The linear viscous losses have a well-known behaviour with frequency, the a-coefficient increases with frequency, as analytically predictable from equation (10). Once the a-coefficient is constrained to its asymptotic linear value, an effective analysis of the behaviour of the b-coefficient with frequency can be performed. Both trends are pictured in Figure 7.

The b-coefficient can be made dimensionless to obtain a Forchheimer term c_F resembling that of the steady-state correlations [equation (4)], using the density and the square root of permeability as length scale. Although this has been done using the analytical permeability for each frequency, the curves confirm the dependence of the b-coefficient on the frequency (Figure 8).

As mentioned in Section 3.4, π₂ can better explain the physics behind the trend pictured in Figure 9. For π₂ → 0, the Forchheimer coefficient tends ideally to zero. This is due to the fluid particle displacement being much smaller than the stack length for a finite Reynolds number. Furthermore, non-linearities caused by geometrical discontinuities are negligible compared to the viscous linear losses due to solid and impermeable walls of the fluid channel. This condition can also be obtained at high frequencies (small particle displacement) or large stack lengths. On the other hand, at low frequency (or stack length, π₂ → ∞), the fluid particle displacement can be comparable or higher than the stack length and non-linearities are more pronounced. CFD data can be fitted by an exponential law passing through the origin and reaching a saturation value, as shown in Figure 9.

4.4 Forchheimer coefficients versus porosity

At 50 Hz, four different porosities, from 0.5 to 0.8, are investigated. As for dimensionless frequency, the permeability has a linear and analytical closure with respect to porosity [equation (10)]. The trend of the Forchheimer coefficient against frequency, pictured in Figure 10, depicts that the lower the porosity the higher the non-linear losses but with a non-linear trend. This is expected because the velocity gradients at low porosities are more significant than those experienced at high porosity. Furthermore, such a finding is in qualitative agreement with the dependence of the concentrated pressure drop caused by sudden contraction/expansion for steady-state flows on the porosity. The CFD-data are fitted with accuracy by a rational law passing through the point of coordinate (ϕ = 1, c_F = 0) with a single constant c as shown in Figure 10.

A summary of the correlations found in this work is provided in Table 3, where the structure of the fitting law and the numerical coefficients are reported for both the viscous and oscillating inertial components of the pressure drop.

4.5 Analysis of the proposed approach

All the above results and the consequent correlations have been obtained according to the procedure described in Section 3.3 by calculating the pressure gradient as the pressure difference (decomposed along the average velocity and evaluated at the fundamental frequency) between the ends of the stack divided by the stack length. The velocity has been calculated as a volumetric average over the stack. This method is straightforward, and it is generally adopted not only in numerical simulations (at least in steady-state flows) but also in experiments, and it allows capturing non-linear losses. In porous and porous-like media with no tortuosity, the non-linearity arises from the transition between regions, which is captured in the proposed method. This is clarified in Figure 11, where three types of pressure drop line graphs are plotted against the superficial velocity. The first one in blue (Procedure i, Figure 11) represents the above-described procedure, i.e. by calculating the pressure drop as the difference between the pressure at the stack ends, considering all non-linear losses including the end effects:

(19) Δpi=|pleft−pright|

The second type is presented by evaluating the local pressure gradient at the middle of the stack multiplied by the length of the stack (red line in Figure 11, Procedure ii) as follows:

(20) Δpii=|dpdx|middle|⋅Lstack

Even though the trend of the pressure drop against velocity is still clearly non-linear, the numerical values are significantly smaller than those in the previous curve as the end effects are not considered. Finally, the linear Darcy pressure drop is shown as a reference by the yellow line in Figure 11 (Procedure iii) and calculated by using the analytical real part of the inverse of permeability:

(21) Δpiii=μ⋅ℜ[1κc]vs⋅Lstack

This is also pictured in Figure 11 to understand how much a simple linear Darcy model can underestimate the total pressure drop in oscillating porous media flows. For the highest velocity, the percentage difference between the adopted approach (i) and the linear reference (iii), that represents the current approach to model this kind of stacks, is approximately 300%.

Furthermore, looking at the distribution of the viscous pressure along the stack (f = 50 Hz, ϕ = 0.6) pictured in Figure 12, for the highest pressure amplitude, the pressure drop at the stack/free fluid discontinuities is comparable to that obtained inside the stack. This would be underestimated not only with a purely Darcy model but also by evaluating the pressure gradient at the middle of the stack. However, if the microstructure of the stack is modelled with a single porous box, the local peaks at the discontinuities between the porous and free fluid regions are not locally reproduced at the macroscopic level but are equally distributed along the stack length.

While the component of pressure in phase of velocity is characterised by an abrupt pressure drop at the stack region, for velocity is not the same, as shown in Figure 13. The “effective” velocity (Darcy velocity dived by the porosity) is plotted for three different pressure amplitudes. From the velocity in Figure 13 and pressure in Figure 12, the acoustic power density (i.e. the mechanical power flux) can be plotted. In the frequency domain, this is equal to the product of pressure, velocity amplitudes and their phase shift θ (Swift, 1988).

(22) W=∫Tp(t)u(t)dt=12|p1||u1|cos⁡(θ)=12p1,viscu1 

From an energy balance perspective, the acoustic power at the left surface (x = 0) shown in Figure 14 represents the net energy entering into the system, even though both pressure and velocity are oscillating with a zero-time average value. This energy budget is converted by the stack into heat that in turn is released to the environment at the temperature level imposed by the boundary condition. The temperature contours (Figure 15) as well as the velocity field (Figures 15 and 16) in the stack region illustrate the phenomena. Both fields are averaged over last five acoustic periods of the simulations. It is apparent that temperature is approximately constant, while some local entrance/exit effects are visible in the time average of velocity vector field. Pairs of counter-rotating vortices are formed both before the stack (in the purely free fluid zone) and inside the stack, as shown in Figure 16.

5. Conclusions

While modelling steady-state flows through porous media has a solid understanding in the literature, oscillatory flow works focus mostly on purely linear Darcy losses, and the equivalent non-linear Forchheimer effects in general porous media model have not been studied with a similar and in-depth systematic approach. To address this gap, this work investigates the presence of non-linear Forchheimer effects in a porous media made of parallel plate channels inserted in a resonator tube. The specific objective of this article is to provide numerical correlations between pressure gradient and velocity and characterise the Forchheimer term against frequency and porosity. Data are obtained by performing parametric CFD simulations at the microscopic level for various frequency or porosity values.

The instantaneous velocity distribution along a fluid channel between the parallel solid plates is validated against experimental data for two different pressure amplitudes. The results of the parametric simulations against the Reynolds number show that the viscous of pressure drop experiences a parabolic trend as in stationary flows. On the other hand, the out-of-phase component of pressure drop remains linear even at higher velocity regimes. For velocity tending to zero, the relationship between pressure-drop and velocity is in perfect agreement with the Darcy law extended to oscillatory flows. While the linear coefficient of the parabolic law has a well-known frequency-dependent pattern, the present work highlights, for the first time, the dependence of the Forchheimer term on the frequency. More specifically, CFD data show that the b-coefficient increases for lower frequencies, tending to an asymptotic value. This can be interpreted as the consequence of increased fluid displacement thus increasing the viscous contribution, while oscillating inertial contributions are dominant at high frequencies. Furthermore, sensitivity analysis by varying the porosity quantitatively demonstrates that the origin of non-linearities for this type of stack having no tortuosity is the geometrical discontinuity at the stack ends: the lower the porosity, the larger the Forchheimer coefficient. Moreover, a local analysis of the correlation between pressure gradient and velocity has shown that the non-linear Forchheimer effects are also visible locally at the middle of the stack.

The study highlights the requirement to account for turbulent onset at the transition into and out of the porous region. While the Reynolds number inside the channel has been always lower than the most accepted definition of the critical Reynolds number in the available literature for oscillating flows, the study could be extended by the use of an appropriate turbulence model to provide further insight into the three-dimensional turbulence structures that arise.

The findings of this study, in addition to a better understanding of oscillatory flows in porous media theory, can be exploited to improve the design tools of thermoacoustic stacks and heat exchangers, especially at high amplitude regimes where non-linear losses must be considered. More specifically, future applications maybe:

• Generalise the results for uniform cross-section stacks, having a uniform cross section in the oscillating flow direction, such as the circular pore and longitudinal pin array stacks.

• Implement the correlations obtained in this work in linear-based proprietary codes/software such as DeltaEC for higher fidelity results (in the frequency domain).

• Implement the same (macroscopic) correlations in CFD-based simulations to reduce the computational costs of a full microscopic approach.

Figures

Figure 1.

Domain and computational mesh of the CFD simulations

Figure 2.

Phasor diagram: pressure gradient (viscous and inertial) and velocity in the frequency domain

Figure 3.

Mesh sensitivity analysis for case p₁ = 10,000 Pa, f0 = 350 Hz, ϕ = 0.5

Figure 4.

Validation against experimental data by Berson et al. (2008)

Figure 5.

Viscous (a) and inertial (b) pressure gradient at f₀ = 350 Hz, ϕ = 0.5

Figure 6.

Dimensionless viscous pressure gradient at different frequencies plotted against the Reynolds number

Figure 7.

a-coefficient and b-coefficient for viscous pressure gradient against frequency

Figure 8.

Forchheimer coefficient (expected value in black and upper/lower bounds with a confidence interval of 95% in blue)

Figure 9.

Forchheimer coefficient vs the second dimensionless number π₂ involving both stack length and frequency

Figure 10.

Forchheimer coefficient against porosity at f₀ = 50 Hz

Figure 11.

Pressure drop against velocity for f₀ = 50 Hz and ϕ = 0.6 according to the three procedures described above

Figure 12.

1D distribution of the viscous pressure at the fundamental frequency 50 Hz, ϕ = 0.6

Figure 13.

1D distribution of the velocity at the fundamental frequency 50 Hz, ϕ = 0.6

Figure 14.

1D distribution of acoustic power the fundamental frequency 50 Hz, ϕ = 0.6

Figure 15.

Time average temperature field

Figure 16.

Time average velocity field

Table 1.

Boundary conditions

Boundary	Momentum equation	Energy equation
Left surface	Oscillating pressure (equation (12))	Adiabatic
Right surface	Wall	Adiabatic
Stack walls	No-slip	Isothermal
Lateral surfaces	Periodicity	Periodicity

Source: Table by authors

Table 2.

Summary of the numerical cases simulated with CFD

CASE	f0 [Hz]	ϕ	p1[Pa]
1–7	10	0.5	100, 500, 1000, 2000, 5000, 10,000, 20,000
8–35	50	0.5, 0.6, 0.7, 0.8	280, 560, 1,125, 2250, 4500, 9,000, 18,000
36–42	100	0.5	100, 500, 1,000, 2000, 3,000, 4,000, 5,000
43–54	200	0.5	10, 100, 500, 1,000, 1,500, 1,800, 2,000
55–63	350	0.5	100, 200, 400, 1,000, 2,000, 4,000, 8,000

Source: Table by authors

Table 3.

Summary of the main numerical correlations

Pressure gradient	Velocity (π1)	Frequency (π2	Porosity (π3)
(dpdx)viscous	aRx + bRx2 aR[Pasm2], bR[kgm4]	aR=μ ℜ[1y02k˜c]	aR=μ ℜ[1y02k˜c]
		bR⋅y0ρ0=0.028(1−1268π2)	bR⋅y0ρ0=0.025(1−π3)π3
(dpdx)inertial	aVx	aV=μ ℑ[1y02k˜c]	aV=μ ℑ[1y02k˜c]
		bV = 0	bV = 0

Source: Table by authors