Adaptive Mesh Refinement Strategies for Cost-Effective Eddy-Resolving Transient Simulations of Spray Dryers

Abstract

The use of adaptive meshing strategies to perform cost-effective transient simulations of spray drying processes is evaluated. These simulations are often computationally expensive, given the large differences between the characteristic times of the central jet and those of the unsteady flow developed at its periphery. Managing the computational cost through the control of the grid resolution by regions is inadequate in many of these applications since the grid resolution requirements change dynamically within the domain. These conditions are related to the unsteady nature of the flow in both the central jet and the flow recirculation zones. Therefore, the application of adaptive mesh refinement (AMR) strategies is recommended. In this paper, general AMR criteria based on relative errors are evaluated by testing three mesh adaptation criteria: velocity gradient, pressure gradient, and vorticity. This evaluation is performed using a low-cost turbulence model with eddy resolution (DDES) in two different types of drying chambers, in which experimental measurements are available. The use of AMR exerts appreciable effects on decreasing computational costs and contributes to the capture of large eddies in critical regions. The present approach provides an appropriate balance between solution accuracy and computational cost. By using a correct AMR configuration, it is possible to obtain results similar to those obtained on a fixed grid but reducing the computational costs by 3 to 5 times.

1. Introduction

Spray drying represents a well-known example of particle-laden flow where a hot-gas jet interacts with a spray cloud to produce high-quality powders. General aspects of spray drying processes are well documented in the literature, e.g., Santos et al. [1]. The interaction between the gas and liquid phases is complex [2] and involves several thermo-physical processes affecting the properties of the powder [3]. These processes are strongly related to the gas-flow dynamics inside the drying chamber [4], and its modeling is essential to predict the performance of a spray dryer [5]. A large number of previous studies using computational fluid dynamics (CFD) have assumed a steady-state approach to model the airflow inside the drying chamber, such as Woo et al. [6], but the current trend towards using a transient approach is now widely accepted [4]. The assumption of a steady-flow field has been found to be inaccurate, as demonstrated by various experimental and CFD studies.

In a pilot-scale spray chamber, Southwell and Langrish [7] used visual techniques, which revealed a swirl-dependent chaotic behavior of the flow in significant portions of the drying chamber. They also visualized a central vortex that exhibited temporal and spatial unsteadiness. Using laser-Doppler anemometry, Southwell and Langrish [8] further investigated the effect of swirl on flow stability. Their findings showed that, for all studied swirl settings (ranging from 0° to 45°) at the inlet, no swirl angle resulted in steady flow behavior within the drying chamber. They also reported the existence of slow, large, and unsteady recirculation regions that altered the flow direction. These flow instabilities are more pronounced in short-form spray dryers [9] with lower swirl values at the inlet jet. Lebarbier et al. [10] employed image analysis techniques (light sheet) to study the flow patterns in a 1/31 scale model of an industrial-sized spray chamber. These observations also revealed a strong time-dependent precessing behavior in the central jet.

Transient CFD studies have also confirmed the time-dependent nature of the airflow. In a CFD study of gas flow inside a pilot-scale spray dryer with and without swirl, Langrish et al. [11] explored the time-dependent nature of flow patterns. In this study, they found that even under zero inlet swirl conditions, the central jet exhibited slight precession around the axis, challenging the validity of the steady-state assumption. In a three-dimension transient study of gas-particle interactions in an industrial-scale spray dryer, Jin and Chen [12] reported low-frequency oscillations of the central jet and three-dimensionality in the flow field, leading to a more realistic dispersion of particles within the chamber. Using a large-eddy simulation approach (LES) to model airflow in a semi-industrial-sized spray chamber, Jongsma et al. [13] reported the existence of large-scale flow structures that exhibited a transient nature, particularly when the central jet was deflected towards the walls. This mechanism is particularly significant for an accurate prediction of the wall deposition of particles.

A further drawback of the steady-state approach is its apparent inability to reach a converged state, a matter debated by Fletcher and Langrish [14] and verified by Lebarbier et al. [10]. This limitation also emerged in preliminary simulations conducted in this study. While some authors have reported residual convergence in steady-state simulations, coarse grids may attenuate these instabilities [9], hence failing to replicate critical flow phenomena within the drying chamber and potentially resulting in misleading results. If sufficient computational resources are available, a transient approach within a three-dimensional domain should be used for CFD simulations of spray dryers [2].

The main difficulty of applying the transient approach lies in its numerical cost [15]. This is due to the ranges of flow time-scales inside the drying chamber [16], which vary from microseconds in the inlet diffuser to seconds in the external recirculation zones and in the precessing motion of the jet. As a consequence, the number of resolved time steps required by the solution must be long enough to develop instabilities, settling the flow to its quasi-steady limits [9] and capturing sufficient information. In previous studies, simulation times vary widely, making it difficult to determine whether conditions for developed flows were achieved or whether run times were selected under computational constraints. Examples of such simulation times reported in pilot- and semi-industrial-size chamber simulations vary between $10\mathrm{s}$ in Fletcher and Langrish [14], Gimbun et al. [17], $20\mathrm{s}$ in Benavides-Morán et al. [18], $90\mathrm{s}$ to $180\mathrm{s}$ in Jongsma et al. [13], $190\mathrm{s}$ in Jin and Chen [12], $197\mathrm{s}$ in Saleh and Nahi Saleh [19], and up to $280\mathrm{s}$ in Langrish et al. [11]. Computational costs become even more relevant for industrial-scale simulations with atomization of liquids, where parametrization and process optimization require running numerous simulation cases. These limitations are a relevant challenge for future numerical development to facilitate the use of the CFD technique on a more routine basis [15].

An alternative to reduce computational costs is the use of adaptive mesh refinement (AMR) strategies. Mesh adaptation is accomplished by dividing or coarsening groups of cells following a refinement criterion [20]. The criterion must be related to the physical features of the flow problem and the turbulence approach used. Examples of criteria discussed in the literature include the ratio of modeled to resolved turbulent kinetic energies [21], the residual velocity magnitude and the vorticity [20], and the gradient-based relative error [22]. Although AMR has not yet been used in spray dryer simulations, it can provide relevant advantages due to the unstable nature of the flow, which causes the position of the main flow structures to change with time. For eddy-resolving approaches, such as detached eddy simulation (DES) and large eddy simulation (LES), the capture of large eddies requires grid adaptation along their path. For this reason, it would be impractical to attempt a mesh optimization using a static-grid approach.

For related cases such as jets and separated flows, the use of AMR has given very good results. Muralidharan and Menon [23] studied a transverse jet into a cross-flow configuration with and without reaction using mixed AMR criteria based on the pressure gradient and the vorticity. The results agreed very well with direct numerical simulation (DNS) data. Mozaffari et al. [24] successfully used fast- and slow-evolving AMR criteria along with a hybrid LES/RANS scheme to model the flow over a backward-facing step. The backward-facing step is a classic flow problem that greatly resembles the discharge of annular jets, typically found in spray drying chambers. Still, for modeling of full-size spray drying processes, there is uncertainty as to which AMR criterion might work best. For example, Wackers et al. [25] mentions that a pressure gradient-based refinement criterion has advantages over a velocity gradient-based refinement criterion when the simulation interest is far from the boundary layer, as is the case of spray dryer modeling. However, Gutiérrez Suárez et al. [26] found that refinement in the high-gradient zones downstream of the discharge, where the high-velocity jet interacts with its surroundings, has a direct effect on the accuracy of the results. As for AMR parameters, such as refinement level and threshold, they are not sufficiently discussed in all studies, and it is assumed that their adjustment varies for each specific application.

This paper proposes the application of AMR methods to reduce computational costs in transient simulations of spray dryers. To the authors’s knowledge, the application of AMR methods in spray dryer simulations is a novel approach. Consequently, there is uncertainty about whether using specific AMR configuration parameters that have demonstrated effectiveness in other applications will be suitable in this context. This study evaluates three refinement criteria and two AMR configuration parameters to model airflow in two distinct spray drying chambers operating with a co-current flow configuration. The criteria for evaluation are based on the relative errors associated with velocity gradient, pressure gradient, and vorticity, while the AMR configuration parameters are the refinement level and threshold. The first drying chamber is of pilot size, has an annular inlet with Re = 19,600, and has been studied by Gutiérrez Suárez [27] and by Benavides-Morán et al. [18]. The second drying chamber is of semi-industrial size, has an annular inlet with Re = 41,200, and was studied by Kieviet [3]. For both drying chambers, the cited studies provide experimental airflow data obtained by hotwire anemometry (HWA), serving as a base for validation. To account for the unsteady nature of the internal airflow, a transient modeling approach in a fully three-dimensional domain is used. Turbulence resolution capability is provided by a delayed-DES model, which is suitable for the simulation of spray dryers [26]. Special attention is paid to the full development of flow instabilities. The primary objective is to determine an AMR configuration that reduces the computational costs of transient CFD simulations of spray dryers while maintaining solution accuracy.

This paper is organized into five sections. Section 2 provides a description of the numerical methods employed in the study. In Section 3, the computational configuration is discussed, including details about the domain, base grid, and boundary conditions. Section 4 outlines the numerical experiments and the simulation methodology employed. The results and discussion are presented in Section 5, where the evaluation of each AMR criterion, threshold values, and refinement limits are discussed. Finally, Section 6 presents the conclusions derived from the study.

2. Numerical Methods and Turbulence Modeling

2.1. Governing Equations and Turbulence Modeling

The governing equations for this problem are the incompressible-flow unsteady Navier–Stokes equations for continuity (Equation (1)) and momentum (Equation (2)):

$$\frac{\partial u_i^{\ast }}{\partial x_i}=0$$

(1)

$$\frac{\partial u_i^{\ast }}{\partial t}+\frac{\partial \left(u_i^{\ast }{u}_j^{\ast }\right)}{\partial x_j}=-\frac{1}{\rho }\frac{\partial p^{\ast }}{\partial x_i}+\left(\left[\nu +{\nu }^{\ast }\right]\frac{\partial u_i^{\ast }}{\partial x_j}\right).$$

(2)

Due to the use of a hybrid LES/RANS turbulence model, these equations may adopt unsteady RANS and filtered (LES) formulations through the domain (see [26,28]). In URANS, an averaging operation is performed on pressure and velocity terms ${\varphi }_n^{\ast }=〈{\varphi }_n〉$, and the modeled viscosity ${\nu }^{\ast }$ is replaced by the turbulent viscosity ${\nu }_T$, which is calculated by the RANS turbulence model. In LES mode, spatial filtering is used, with ${\varphi }_n^{\ast }=\overline{{\varphi }_n}$, and the modeled viscosity ${\nu }^{\ast }$ is replaced by the sub-grid scale (SGS) viscosity, with ${\nu }^{\ast }={\nu }_{\mathrm{SGS}}$, calculated by the delayed-DES model.

The delayed-DES (DDES) model is classified as an interfacing LES/RANS turbulence formulation [29]. This model is characterized by an improved behavior in near-wall regions, avoiding an early switch from RANS to LES mode [30]. The early switching may happen in conditions of low-flow instabilities or thick boundary layers. This near-wall functionality under very coarse grids has been verified in a prior work by Gutiérrez Suárez et al. [26]. In the SST-DDES model, the transport equation for the turbulent kinetic energy (k) is defined as:

$$\frac{\partial \left(\rho k\right)}{\partial t}+\frac{\partial \left(\rho u_{j}k\right)}{\partial x_j}=\rho p-D_{\mathrm{DDES}}^k+\frac{\partial \left[\left(\mu +{\mu }_t{\sigma }_k\right)\frac{\partial k}{\partial x_j}\right]}{\partial x_j}.$$

(3)

Here, the term $D_{\mathrm{DDES}}^k$ represents a relation between RANS and LES turbulence length scales:

$$D_{\mathrm{DDES}}^k=\rho k^{3/2}/l_{\mathrm{DDES},}$$

(4)

with

$$l_{\mathrm{DDES}}=l_{\mathrm{RANS}}-f_{\mathrm{d}\max (0,l_{\mathrm{RANS}}-l_{\mathrm{LES}})}.$$

(5)

In this equation $l_{\mathrm{RANS}}$ is the length scale defined for the SST $k-\omega$ model as $l_{\mathrm{RANS}}=\sqrt{k}/\left(C_{\mu }\omega \right)$. $l_{\mathrm{LES}}$ is defined in terms of a model constant $C_{\mathrm{DES}}$ and the local maximum grid spacing $\Delta$ as $l_{\mathrm{LES}}=C_{\mathrm{DES}}\Delta .$ The empirical delay function $f_{\mathrm{d}}$ protects the boundary layers from switching into LES mode. This function includes a correction ratio, which includes the molecular and turbulent viscosities, the strain-rate and vorticity vectors, and the distance to the nearest wall. Additional details of the DDES model, including the values of its constants, can be found in Gritskevich et al. [31].

2.2. Adaptive Mesh Refinement

The quantitative error analysis performed by Gutiérrez Suárez et al. [26] when numerically studying an annular swirling jet gives insight about key zones, where refinement would have a larger impact on improving the solution accuracy. These zones were found to be high-shear regions, where the high-velocity jet interacts with its lower-velocity surroundings and regions with high velocity and pressure fluctuations, such as those found beneath the bluff-body. In order to prioritize the refinement in these regions, three sensors are used for adaptation, of which two are based on gradient (velocity and pressure fields) and the third on vorticity.

A general refinement strategy based on Bathe and Zhang [22] is implemented in this work. This criterion determines whether or not an element should be refined/unrefined according to its local error. For any grid element, the local error $c_l$ is defined in terms of its mean edge length ${\Delta }_l$ and the norm of a flow-solution variable in the element $F_e$ as:

$$c_l={\Delta }_l\parallel F_e\parallel .$$

(6)

The flow-solution variable $F_e$ will be defined in this paper in terms of the velocity gradient $\nabla U$, the pressure gradient $\nabla p$, or the vorticity $\mathsf{\Omega }$. For the whole domain, an averaged error $c_{\mathrm{avg}}$ can be determined from the local error $c_l$ of each of the grid elements as:

$$c_{\mathrm{avg}}=\frac{1}{N_e}\sum _e{\Delta }_l\parallel F_e\parallel ;$$

(7)

here, $N_e$ represents the actual number of grid elements. Defining the error relation for any grid element as:

$$c_{\mathrm{rel}}=c_l/c_{\mathrm{avg}},$$

(8)

the decision whether to refine a grid element or not is established in terms of a given refining criterion $c_{\mathrm{crit}}$. If $c_{\mathrm{rel}}>c_{\mathrm{crit}}$, the local element is refined. The selection of this refinement criterion is dependent on the case and physical phenomenon represented. Its value should be weighted according to the maximum level of refinement and the computational resources available.

Cell unrefinement occurs when $c_{\mathrm{rel}}<c_{\mathrm{un}}$, where $c_{\mathrm{un}}$ represents an unrefinement criterion, assumed in this study as $c_{\mathrm{un}}=0.75c_{\mathrm{crit}}$. This selection is experience-based, enabling a balance between the number of refined elements and the desired refinement coverage. Increasing this value may amplify instabilities and lead to additional numerical fluctuations, whereas decreasing this value would keep cells refined where the relative error is not significant, resulting in inefficiency in the grid adaptation process. The selected value presented a good tradeoff between these discussed effects across all the studied cases.

Because the value of $c_{\mathrm{rel}}$ changes for different flow velocities and annular jet configurations, in this paper, a better definition is given using the Reynolds number as:

$$\begin{array}{c}{C}_{\mathrm{crit}}=c_{\mathrm{crit}}Re,\\ C_{\mathrm{rel}}=c_{\mathrm{rel}}Re.\end{array}$$

(9)

2.3. Numerical Schemes and Methods

The open-source CFD software OpenFOAM V1706 was used to solve the governing equations (Equations (1) and (2)). The viscous terms were resolved using DES-based turbulence models. The discretization methods for different terms are presented in

Table 1. Numerical schemes and coupling methods used in the present study.

Schemes	Description
Temporal terms	SOUE (Adams Moulton), second-order implicit, described in Moukalled et al. [33].
Gradient—$\nabla \varphi$	CDS (central differencing scheme).
Gradient—$\nabla \mathit{u}$	Limited-CDS.
Advective terms—div$\left(\varphi \right)$	Hybrid scheme, described in Gutiérrez Suárez et al. [26], Travin et al. [34].
RANS mode	LUD (linear upwind differencing).
LES mode	CDS.
Pressure–velocity coupling	PIMPLE algorithm, described in Holzmann [32].

. All temporal and spacial terms used are second-order accurate. The pressure–velocity coupling was achieved using the PIMPLE algorithm. The PIMPLE algorithm, described in Holzmann [32], is a combination of SIMPLE and PISO algorithms. This algorithm allows running transient problems under large Courant numbers. For pressure and velocity, the settings used for residual control include a solver tolerance of $1.0\times {10}^{-05}$, a relative tolerance of $0.001$, and a maximum number of iterations of 10 (per outer PIMPLE loop). Adaptive mesh refinement (AMR) is implemented through relative error criteria (see Equations (6)–(9)). These equations were incorporated into the OpenFOAM solver pimpleDyMFoam and compiled as a new solver for the purposes of the present study.

3. Computational Configuration

3.1. Flow Configuration and Domain

The modeling of the airflow (cold flow) is performed in two spray drying chambers, whose geometries and domains are presented in Figure 1. Both drying chambers are of co-current type and have an annular inlet. A significant portion of the annular duct is included in order to develop the flow before it is discharged into the chamber. The flow characteristics and annular jet geometries for both dryers are presented in

Table 2. Flow and geometry characteristics of the annular entrances for the pilot-size and semi-industrial-size spray drying chambers.

Description	Pilot-Size Spray Drying Chamber	Semi-Industrial Spray Drying Chamber
Volumetric flow	$0.0342$ ${\mathrm{m}}^3{\mathrm{s}}^{-1}$	$0.377$ ${\mathrm{m}}^3{\mathrm{s}}^{-1}$
RPMs at the discharge	238	0
Swirl number	0.052	0
$R_o$	0.0355 m	0.244 m
$R_i$	0.0260 m	0.2024 m
$D_h=4A_{wet}/2\pi D_{avg}$	0.0190 m	0.084 m
Re	19,600	41,200
Reference air density	$1.1$$\mathrm{kg}{\mathrm{m}}^{-3}$	$1.23$$\mathrm{kg}{\mathrm{m}}^{-3}$
Kinematic air viscosity	$1.42\times {10}^{-5}$${\mathrm{m}}^2{\mathrm{s}}^{-1}$	$1.48\times {10}^{-5}$${\mathrm{m}}^2{\mathrm{s}}^{-1}$

. For the semi-industrial-size dryer, the flow velocity measured at the discharge [35] was used as a reference. With this value and assuming a developed turbulent velocity profile [36], the volumetric flow rate upstream $Q=0.377$ ${\mathrm{m}}^3{\mathrm{s}}^{-1}$ was calculated. This value is approximately 10 times higher than that of the pilot-size chamber. Because Kieviet et al. [35] removed the upper section of the air disperser, where the swirl vanes are located, this case was modeled without swirl.

For the pilot-size dryer, the values reported in Benavides-Morán et al. [18] were used, with the exception of the swirl number at the jet discharge, which was obtained from the experimental measurements of Gutiérrez Suárez [27]. The swirl number $S_n$ at the discharge was calculated using the reported value of 238 RPM through the following equation:

$$S_n=\frac{{\int }_0^R\rho U_x{U}_{\theta }{r}^{2}dr}{R{\int }_0^R\rho U_x^{2}rdr}.$$

(10)

Table 2 presents the flow and geometry characteristics of the annular entrances for the pilot-size and semi-industrial-size spray drying chambers.

The same boundary condition (BC) types were used for both drying chambers. For the air inlet, a constant volumetric flow condition was used. As was done in Gutiérrez Suárez et al. [26], turbulent variables k and $\omega$ were defined through a fixed-value BC using very low initial values ($1\times {10}^{-6}).$ These values are expected to develop throughout the inlet duct, all the way to discharge. In the walls, the velocity was defined through a no-slip BC, which enforces a Dirichlet-type BC $(U_{x,\theta ,r}=0)$. Turbulence variables k and $\omega$ were defined in this zone through standard wall functions.For the outlet, a zero-gradient boundary condition was applied to the velocity and turbulence variables, while the pressure was defined through a fixed-value BC, with $p_{\mathrm{outlet}}=0.$

3.2. Grids

The base grids for both drying chambers are presented in Figure 2. Both grids are structured and include a prerefined (to level 1) region covering the inlet duct. (The level of refinement refers to the number of refinement operations performed on the base grid element and subsequently to the new elements created after each refinement operation. Therefore, a level 2 refined cell is created from a base cell (level 0) that was refined 2 times. During this operation, 64 cells are created from the original cell.) The operation was carried out to improve the accuracy and stability of the solution in the near-field of the jet during the initialization of the cases. Upon activation of the AMR, the refinement of this region was directly controlled by the mesh adaptation algorithm.

For both grids, the elements are concentrated in the jet discharge and in the outlet duct. A positive growth factor from this point towards the radial boundaries and the outlet was implemented. The growth factor was adjusted to maintain values of $y^{+}$ between 100 and 300 on the chamber walls following activation of the AMR, assuming a flow recirculation velocity of 0.5 m/s. Due to the use of adaptive meshing, the base grids (level 0, without prerefinements) are very coarse: the mesh of the pilot-size chamber has 27,550 cells and the grid of the semi-industrial-size chamber has 39,540 cells. The distribution of elements in axial, tangential, and radial directions is presented in

Table 3. Total and direction-wise number of elements for the level 0 base grids of the pilot-size and semi-industrial-size drying chambers (without prerefinement). The elements of the air disperser are not included in the directional count.

Spray Drying Chamber	Elements	x	$\mathit{\theta }$	r
Pilot-size chamber	27,505	88	19	16
Semi-industrial-size chamber	39,540	52	19	28

. For the pilot-size chamber dryer, the mesh is completely structured and is composed of hexahedral cells. For the semi-industrial-size chamber, the outlet pipe crosses the conical section of the chamber and makes the construction of the mesh more difficult. In a previous study, Anandharamakrishnan et al. [38] used a hexahedral mesh for the cylindrical section and a tetrahybrid mesh for the bottom conical section. In the present study, a structured mesh composed of hexahedral cells is used as the background mesh, and the outlet pipe is incorporated using split-hexahedral cells.

Due to AMR, the number of grid elements varies dynamically depending on the value of the non-dimensional refining criteria $C_{\mathrm{crit}}$ (see Equation (9)) and on the maximum refinement level allowed ${\mathrm{ref}}_{\max }$. The aim of grid adaptation is to produce a numerical solution equivalent in accuracy to that which could be obtained if the entire grid was refined but using fewer elements. For this reason, grid independence is not evaluated by consecutively increasing the overall element level but by testing different AMR configurations and their effect on the accuracy of the results.

For reference, a global refinement of the base grids (up to level 2) would result in a number of elements comparable to that reported in other studies (see

Table 4. Grid sizes, turbulence model, and time treatment in previous CFD studies of the pilot-size and semi-industrial-size drying chambers. For comparison purposes, the values used in the present study include a globally refined level 2 base grid. See Table 5 for the effective grid sizes when using AMR.

Drying Chamber	Author	Year	Grid Size	Turbulence	Time
Pilot size	Gutiérrez Suárez [27]	2015	12,672 (2D-axi-symmetric)	$k-ϵ$	steady
	Benavides-Morán et al. [18]	2020	2,200,000	SST $k-\omega$ and SAS	transient
	Present study (level 2 refined base grid)	2023	1,750,000	SST-DDES	transient
Semi-industrial size	Huang et al. [40]	2006	not reported	$k-ϵ$	steady
	Anandharamakrishnan et al. [38]	2010	294,237	$k-ϵ$	steady
	Saleh [39]	2010	840,000	$k-ϵ$	transient
	Noor and Jolius [41]	2011	503,000	SA-DES	transient
	Saleh and Nahi Saleh [19]	2011	516,000	RNG $k-ϵ$	transient
	Gimbun et al. [17]	2015	420,000	SA-DES	transient
	Keshani et al. [42]	2015	1,223,675 (Case A)	RNG $k-ϵ$	transient
	Present study (level 2 refined base grid)	2023	736,000	SST-DDES	transient

). For the pilot-size chamber, the resulting number of elements is similar to that reported by Benavides-Morán et al. [18]. For the semi-industrial-size chamber, the number of elements is higher than that reported by most previous studies and is similar to that used by Saleh [39]. However, the effective grid sizes when using AMR are significantly smaller (see

Table 5. Study cases for the pilot-size and semi-industrial-size spray chambers.

Pilot-size drying chamber
Case	Number of cells	AMR flow variable	$C_{\mathrm{crit}}$	Max. refinement level	Turbulence
C (base)	1,750,000	fixed grid	-	Level 2 globally refined base grid	DDES
C1	300,000	$\nabla U$	370 (ref)	2	DDES
C2	300,000	$\nabla p$	370	2	DDES
C2-1	78,000	$\nabla p$	3700	2	DDES
C2-2	180,000	$\nabla p$	3700	3	DDES
C2-3	380,000	$\nabla p$	370	2 using a level 1 globally refined base grid	DDES
C3	300,000	$\mathsf{\Omega }$	370 (ref)	2	DDES
Semi-industrial-size drying chamber
S1	170,000	$\nabla U$	12,360	2	DDES
S2	220,000	$\nabla U$	4120	2	DDES
S3	378,000	$\nabla U$	1236	2	DDES

).

4. Numerical Experiments

4.1. Study Cases

Table 5 presents the main characteristics of the studied cases. The initial seven cases are specifically centered around the pilot-size drying chamber. Case C, serving as a base case, employs a fixed mesh (without AMR) consisting of 1.75 million elements, achieved by a global refinement of all elements of the base grid to level 2. This reference case enables the evaluation of results variations upon the activation of AMR.

Cases C1, C2, and C3 evaluate the effect of using different flow-solution variables as the basis for AMR. In these cases, the maximum refinement level is 2, and the refinement threshold is adjusted to obtain around 300,000 grid elements. Cases C2-1 and C2-2 evaluate the effects of the threshold values and maximum refinement level over the solution. Case C2-1 maintains the variable flow solution of case C2 ($\nabla p$) but uses a higher refinement threshold, significantly reducing the number of typical grid elements. Case C2-2 also uses the same variable flow-solution as C2 and C2-1, but the maximum level of refinement is set to level 3. In this case, the refinement threshold is significantly increased to reduce the number of refined elements and obtain a computational cost similar to cases C1, C2, and C3. For Case C2-3, a fixed prerefinement is performed on the entire base grid (to level 1) and a further level of refinement is allowed (up to level 2). For the semi-industrial-size chamber, three cases were studied using the same variable flow-solution ($\nabla U$). These cases were selected to determine the effect of the threshold value on the numerical results in a different drying chamber.

For all cases with AMR, the mesh adaptation was triggered every five time steps ($t_{\mathrm{AMR}}=5\Delta t$), with a two-cell refinement buffer around refined cells. This choice was experience-based and aimed to balance computational efficiency while preventing information loss caused by a slow response of the adaptation method to rapidly changing flow features. The maximum Courant number was set to ${\mathrm{Co}}_{\max }=0.85$, primarily reached at the jet discharge, where the refinement behavior tended to be uniform. In other regions within the drying chamber, the maximum Courant number was significantly lower than the inverse of the adaptation time (${\mathrm{Co}}_{\max }<1/t_{\mathrm{AMR}}$), ensuring adequate time for the AMR method to react to dynamic flow changes. The choice of ${\mathrm{Co}}_{\max }$ was based on findings from a previous study [26]. While the PIMPLE algorithm permits the use of higher ${\mathrm{Co}}_{\max }$ values, configuring it in this manner is not recommended due to the generation of artificial fluctuations at the jet discharge and an excessive dissipation downstream [26].

4.2. Simulation Times and Sampling of Data

The simulation times in this study are described in terms of a non-dimensional characteristic time $t_{\mathrm{c}}$ defined as:

$$t_c=\frac{t_s}{t_f}=\frac{t_s}{Q/V_c}.$$

(11)

Here, $t_s$ represents the time elapsed since the start of the simulation, $t_f$ is the time to fill the chamber, Q is the volumetric flow rate at the inlet, and $V_c$ is the volume of the chamber. For the pilot-size chamber, $t_f=12.81$ s, and for the semi-industrial-size chamber, $t_f=23.85$ s. The characteristic time $t_{\mathrm{c}}$ represents the non-dimensional time required to fill the chamber once. All cases are initialized from $t_c=0$ to $t_c=0.5$ without AMR and using the basic grid presented in Figure 2. Then, from $t_c=0.5$ to $t_c=2.0$, all the elements are refined to level 1. After $t_c=2.0$, the AMR method is activated and the cases are run up to $t_c=5.0$ using the specific configurations described in Table 5.

Data collection by the probes is active during the entire simulation. However, only the data obtained after $3.5t_c$ are used for the analysis of the internal flow pattern. The criteria to determine these sampling times are presented in the results Section 5.1. The data collected by the probes are the instantaneous flow velocity and the turbulent kinetic energy calculated by the RANS component of the SST-DDES model. These data are sampled at 100 Hz using probes distributed in different axial and radial positions. Figure 3 presents the axial sampling positions and radial paths for both drying chambers. For the pilot-size chamber, the measured axial positions are $x=(0.1,0.2,0.3,0.6)$m. A total of 15 probes were set per axial position. Because the flow pattern in this drying chamber is asymmetric (as the annular duct is not concentric), the same measurement positions reported in previous studies [18,27] were used, varying the tangential position for each radial position. For the semi-industrial-size chamber, the axial positions are $x=(0.3,0.6,1.0,2.0)$m, with the same radial arrangement used by Kieviet et al. [35]. The mean-velocity profiles were obtained by averaging the data sampled by each probe. This sampling frequency was ensured by a dynamic adjustment of the simulation time step size $\Delta t$.

5. Results

5.1. Computational Costs

The simulations were conducted on the JFF2 cluster at the CTTC (Heat and Mass Transfer Technological Center) of the Polytechnic University of Catalonia. This cluster comprises a total of 128 nodes, each equipped with 2 AMD Opteron Quad Core processors and 16 GB of RAM. These nodes are interconnected via an Infiniband DDR 4X network with a latency of 2.6 microseconds and a bandwidth of 20 Gbits/s. Cases C1 to C3 and S1 to S3 were executed using six processors, while cases C2-1 and C2-2 used four processors. For case C (base case), 16 processors were employed. This selection of cores per case was made to maintain an optimal range of cells per processor across all computational cases. Consequently, it helped reduce bottleneck effects associated with parallel computing, enabling a consistent comparison of computational costs. Further information on HPC scalability and bottleneck effects in OpenFOAM can be consulted in Culpo [43].

Figure 4 presents the computational costs in terms of effective processor hours ($PH$) per characteristic time $t_{\mathrm{c}}$. For similar grid sizes, costs are slightly higher in the pilot-size dryer than in the semi-industrial dryer. This is probably because its annular section is comparatively smaller and thinner, thus requiring a smaller time step. Compared to the base case (C), the reduction in computational cost when using AMR is between 3.2 and 11 times. Case C2-2 (up to three refinement levels) is an exception. This indicates that solving a larger fraction of the turbulence brings with it a substantial increase in computational cost, even if fewer grid elements are used.

Compared to the base case (C), the reduction in computational costs primarily results from a decrease in the number of grid elements, significantly lowering the cost of each time step. On the other hand, the cost is not affected by the time step size, which remains unchanged. This occurs because the regions with the largest influence on the time step size (where the Courant number is the highest) tend to remain refined throughout the simulation. There is a minor additional cost associated with running the AMR algorithm, which counts as running an additional time step in terms of cost.

5.2. Flow Development and Sampling Times

This section presents the results of the time-dependent flow development in the pilot- and semi-industrial-size drying chambers. These results support the criteria to sample the data presented in Section 4.2. The development times are represented by a characteristic time $t_{\mathrm{c}}$ and are presented for cases C and S2 in Figure 5. The mean axial velocity $U_x$ and its rms fluctuations are used to evaluate the evolution of the flow development. These values are presented in intervals using a moving window average.

Figure 5a presents the mean $U_x$ results for the pilot-size drying chamber. At the centerline ($r/R=0$), a stabilization in the sampled velocity values is observed after $t_{\mathrm{c}}=2.5$. At the recirculation region ($r/R=0.96$), there is an increasing trend in the velocity magnitude, which means that the recirculation zone is becoming stronger. For the same drying chamber, rms fluctuations of the axial velocity $U_x$ are presented in Figure 5b. It is observed that after $t_{\mathrm{c}}=2$, the values of the rms fluctuations in the centerline and the recirculation region achieve a semi-stable behavior.

Figure 5c presents the mean $U_x$ results for the semi-industrial-size drying chamber. In this case, a clear trend is not observed in the time evolution of the results of the measurement points at $r/R=0$ and $r/R=0.866$. A slight decrease in the axial velocity values is observed in the central jet around $t_{\mathrm{c}}=3$, which is related to an increase in velocity magnitude in the recirculation zone. However, after this point, the velocity shows an increasing trend with a wavelike pattern. Because of this, the measurement was extended from $t_{\mathrm{c}}=5\mathrm{to}{t}_{\mathrm{c}}=8$. The results of this extended measurement (average velocity over the range) are shown as pointed horizontal lines in Figure 5c. It is observed that the recirculation flow becomes stronger and that the velocity in the central zone ($r/R=0$) takes an intermediate value. Figure 5d presents the sampled values of rms fluctuations. A clear increasing trend is observed at the center point ($r/R=0$) and recirculation zone ($r/R=0.866$). The increasing trend is more marked at the recirculation point, although the slope decreases slightly after $t_{\mathrm{c}}=1.2$.

These results show that the flow structure inside the drying chamber requires considerable time to fully develop. They also present evidence that the fluctuations are self-sustained by finding a correlation between the fluctuations in the recirculation zone and at the central point. Based on a qualitative analysis of the results, for both cases, after $3.5t_c$, a decrease in the variation of axial velocities and rms fluctuations per segment is observed. Therefore, this time is selected to start the sampling of data. Given the results presented here, where the most significant variations in velocities and rms fluctuations occur before $t_{\mathrm{c}}=1$, the choice of measurement times is considered sufficient for the purpose of this study. For simulations of other spray dryers under different operating conditions, a trade-off between computational costs, grid refinement, and measurement time should be made.

5.3. Typical Refinement Regions

5.3.1. Pilot-Size Spray Chamber

Figure 6 presents the transient behavior of the airflow inside the pilot-size chamber obtained from the base case. The ability of the model to resolve turbulence is observed, in which detached eddies are captured just downstream of the discharge (see Figure 6a). These eddies grow in size downstream as the jet decays. Regarding the velocity along the radial direction at $x/H=0.367$ (see Figure 6b), strong flow asymmetry is observed from the centerline. This asymmetry is a usual behavior in actual spray dryers, but in this particular case, it is most likely amplified by the defect in the concentricity of the annular entrance (see Figure 1c). This leads to a strong deflection of the jet, which sometimes can result in velocities close to 0 in the centerline.

Figure 7a–c presents the effect of the refinement criteria ($\nabla U,\nabla p,\mathsf{\Omega }$) on the relative error fields $C_{\mathrm{rel}}$. This visualization is important because it gives insight into the nature of each of the criteria and the regions that would be selected for refinement ($C_{\mathrm{rel}}>C_{\mathrm{crit}}$) in case of using adaptive meshing.

The first evaluation of these relative error fields is performed on a computational solution developed on a fixed grid (case C), which is presented in Figure 7a,b. In the near-field of the jet (Figure 7a), the differences are mainly observed in two regions: in the inlet diffuser, including the jet discharge, in the air disperser, and in the cells near the wall. The regions marked for refinement by the vorticity and velocity gradient criteria are quite similar. Both criteria refine cells along a large section of the inlet duct, but this refinement is only performed in near-wall cells. These criteria also exhibit high values of $C_{\mathrm{rel}}$ in the high shear stress zone downstream of the discharge, where the high-velocity flow interacts with the surroundings. In some cases, near-wall cells on the upper walls of the chamber are marked for refinement. This behavior is undesirable because such cells do not interact with high-flow velocities. It occurs mainly because of the high aspect ratio of the elements in this zone (see Figure 2a). Since the mean element length ${\Delta }_l$ in Equation (6) is calculated from the cubic root of its volume ${\Delta }_l=\sqrt[3]{V_c}$, in cells with a high aspect ratio, the magnitude of the local error is overpredicted. Regarding the pressure gradient criterion, high values of $C_{\mathrm{rel}}$ are presented in the entire lower section of the air disperser, near the jet discharge. This behavior differs from that of the other two criteria. Additionally, it does not present high values of $C_{\mathrm{rel}}$ in the upper wall due to the condition imposed $\nabla p=0$ on the walls, which hinders the refinement of these cells.

Downstream of the jet discharge (Figure 7b), the $C_{\mathrm{rel}}$ fields generated by the $\nabla U$ and $\mathsf{\Omega }$ criteria remain similar. In both cases, some elements near the lower walls are refined, which is desirable if high velocities are present in those areas. In detail, the $\nabla U$ criterion presents a slightly higher coverage of the far zone of the jet. To reach the same coverage, the $\mathsf{\Omega }$ criterion would require a $C_{\mathrm{crit}}$ threshold significantly lower than 240.

Figure 7c presents the effective $C_{\mathrm{rel}}$ fields obtained using AMR for cases C1, C2, C3, C2-1, and C2-2 taken at $t_c=3.0.$ In general, for cases C1 to C3, the behavior in the inlet duct zone and at the jet discharge is very similar to that obtained for the base case (Figure 7a). Downstream, however, the zone selected for refinement ($C_{\mathrm{rel}}>C_{\mathrm{crit}}$) exhibits a smaller extent, covering less than half the height of the chamber and terminating abruptly. This limitation arises because the base grid in this instance is too coarse to raise the local error above the refinement threshold. Regarding C2-1 and C2-2 cases, they presented a much lower coverage, especially for case C2-2, limiting itself to the discharge and the very-near region of the jet.

Figure 7d presents the effective refined regions (in terms of the cell level) for all the cases related to the pilot-size dryer. A two-cell refinement buffer is implemented around each refined cell to ensure smooth grid transitions. For cases C1 to C3, this approach effectively enhances the refinement area around the jet discharge, containing more cells than those initially selected for refinement ($C_{\mathrm{rel}}>320$). However, this strategy proves less effective for the other cases (C2-1 and C2-3), particularly for case C2-2, where a refined region is not retained at the jet discharge. Because the level of refinement in the discharge is crucial for modeling the development and decay downstream of the annular jet, this behavior is deemed unsatisfactory. Downstream of the discharge, the extent of the refined zone is similar to that selected for refinement ($C_{\mathrm{rel}}>C_{\mathrm{crit}}$).

Overall, case C2-3 presents a superior behavior to all the other AMR cases in this region by avoiding an abrupt transition between the refinement zones. This behavior aligns with the $C_{\mathrm{rel}}$ fields observed for the base case (Figure 7a,b), further affirming the previously mentioned influence of the background grid on grid adaptation behavior. Near the output duct, the pressure gradient criterion generates a homogeneous refinement, which contrasts with the velocity gradient and vorticity criteria exhibiting similar operating characteristics to those observed in the inlet duct.

5.3.2. Semi-Industrial-Size Spray Chamber

Figure 8 presents the temporal evolution ($t_c=4.5,5.5,6.5$) of velocity contours and refinement regions for case S2. Figure 8a shows the instantaneous contours of absolute velocity. In the high-velocity jet, the formation of detached eddies is observed just downstream of the annular blockage. The size of these eddies increases downstream as the jet widens and decays. Some recirculatory phenomena in the periphery of the jet are also observed. These flow phenomena may be related to the zones chosen for refinement, as shown in Figure 8b. These zones are presented as contours of the relative error $C_{\mathrm{rel}}$. As shown in Equations (8) and (9), the values of $C_{\mathrm{rel}}$ and $C_{\mathrm{crit}}$ determine the refining operations through the domain. As expected, the higher relative errors are present in high-shear zones, which are related to the jet, eddies, and the walls of the domain. Compared to the pilot-size dryer, in this case, the entire jet discharge region is not selected for refinement. This is mainly due to two factors: the low production of turbulence upstream (a turbulence intensity at the discharge of 2%), which is influenced by the absence of inlet swirl, and a slow reattachment of the jet. This results in low-velocity gradients in the central section and hence low levels of refinement. The delay in jet reattachment is further exacerbated by the open bluff-body cover (it is not completely clear in the article of Kieviet et al. [35] if the cover of the bluff-body was left open; it was assumed that it was from the interpretation of the drying chamber figures and the computational grid presented in that study), generating large pressure fluctuations that prevent the two sections of the annular jet from merging.

Downstream of the discharge, some regions marked for refinement follow the large-detached eddies resolved by the DDES model. These refinement regions travel with the eddy until it loses energy and is dissipated. Regions marked for refinement are also observed around the outlet duct, including the pipe walls and below the outlet duct, where the flow experiences a sudden change in direction. Recirculation regions are not selected for refinement due to the low energy of recirculating eddies. However, some near-wall elements are marked for refinement in regions where flow impingement occurs. This occurs especially in the section where the shape of the chamber changes from cylindrical to conical and in the areas opposite the outlet duct. A correct capture of the flow behavior in these regions appears to have a significant effect on flow stability, as mentioned by Keshani et al. [42], making their refinement desirable.

5.4. Turbulence Resolution at the Jet Discharge

Table 6. Turbulence-related parameters measured at the discharge of the jet for all study cases of the pilot-size spray drying chamber.

Reference	rms ${\mathit{U}}_{\mathit{x}}$	${\mathit{k}}_{\mathbf{LES}}$	${\mathit{k}}_{\mathbf{RANS}}$	k	$\mathit{\Gamma }$
C	0.0052	0.0186	0.0789	0.097	0.191
C1	0.0073	0.0271	0.2443	0.271	0.100
C2	0.0141	0.0233	0.3289	0.352	0.066
C3	0.0110	0.0265	0.0950	0.122	0.218
C2-1	0.0146	0.0156	0.0247	0.040	0.387
C2-2	0.0416	0.1239	0.0492	0.173	0.716
C2-3	0.0711	0.0225	0.07508	0.095	0.237

presents some parameters related to turbulence resolution measured at the jet discharge for the pilot-size dryer. These parameters are the rms fluctuations of the axial velocity (rms $U_x$), the resolved turbulent kinetic energy ($k_{\mathrm{LES}})$, the modeled turbulent kinetic energy $k_{\mathrm{RANS}}$, the total turbulent kinetic energy k, and the resolved turbulent kinetic energy fraction $\Gamma$. The resolved fraction is calculated from the resolved turbulent kinetic energy $k_{\mathrm{LES}}$ and the modeled turbulent kinetic energy $k_{\mathrm{RANS}}$ as:

$$\Gamma =\frac{k_{\mathrm{LES}}}{k_{\mathrm{LES}}+k_{\mathrm{RANS}}}.$$

(12)

According to Pope [21], for a well-resolved LES simulation, $\Gamma >0.8$. For the base case (C), it is observed that the resolved fraction of turbulence is very low compared to a LES study ($\Gamma$ = 0.191). This means that only a small number of very large turbulent structures are being resolved and that most of the turbulent kinetic energy corresponds to the modeled fraction ($k_{\mathrm{RANS}}$). The low fraction of resolved turbulence is a consequence of the coarseness of the grids and the very low inlet swirl values. Increasing this resolved fraction would imply radically increasing the number of elements in the inlet duct region, which results in very large computational costs.

The effect of using the different AMR criteria (C1, C2, C3) is presented next. For rms $U_x$, $k_{\mathrm{LES}}$, and $k_{\mathrm{RANS}}$, a substantial increase is observed when using any of the criteria, especially for $\nabla p$ (case C2), where the generated rms $U_x$ fluctuations are much higher than those predicted with the fixed grid, and results in a generalized increase in the resolved turbulent kinetic energy $k_{\mathrm{LES}}$. Since the resolution of the grid is much the same as that of the reference case (case C), it is evident that a problem of artificial fluctuations is occurring, most likely due to drastic changes in the level of refinement in this critical zone. This leads to interpolation errors, which, in turn, cause artificial fluctuations. The increase is more marked in case C2 $(\nabla p$ criterion), possibly because, upstream of the discharge, the grid level changes more abruptly than for the other cases, as shown in Figure 7. For cases C1 and C3, the increase in rms $U_x$ does not produce an increase in $\Gamma$ because the increase in $k_{\mathrm{RANS}}$ is proportionally higher, making the $\Gamma$ value lower than that of the base case (C). In case C3 ($\mathsf{\Omega }$ criterion), the increase in resolved turbulent kinetic energy $k_{\mathrm{LES}}$ is greater than the increase in modeled turbulent kinetic energy $k_{\mathrm{RANS}}$, resulting in a total turbulent kinetic energy value similar to that of the base case.

Regarding cases C2-1 and C2-2, using a pressure gradient refinement criterion, and different threshold and maximum refinement level values, the following results are obtained. For case C2-1, using a larger threshold and consequently a smaller number of elements than in case 2, the rms $U_x$ values are increased. This occurs due to the lower coverage of refined regions in the input diffuser, limiting themselves to the jet discharge. Therefore, near the jet discharge, sudden changes in the refinement level and the mode of operation (RANS or LES) occur constantly, generating strong artificial oscillations. On the other hand, since most of the diffuser is unrefined to an extremely coarse, level 0 grid, the production of modeled turbulent kinetic energy $k_{\mathrm{RANS}}$ is reduced. This is due to the incorrect calculation of wall shear stresses and low shear stresses between cell faces. The resulting overprediction of rms $U_x$ and underprediction of $k_{\mathrm{RANS}}$ artificially increase the resolved fraction $\Gamma$. For case C2-2, where a higher level of refinement (level 3 refinement) is allowed, something similar happens. The resolved fraction of turbulence reaches a value of $\Gamma =0.71$, which is close to the required value for LES. Compared to the base case, the resolved turbulent kinetic energy $k_{\mathrm{LES}}$ increases 12 times, and the modeled turbulent kinetic energy $k_{\mathrm{RANS}}$ decreases by half. This increase in the resolved fraction is related to the finer grid but also to the same artificial fluctuations presented in case C2-1. For case C2-3, running under a more refined background grid, the values of $k_{\mathrm{RANS}}$ and k are similar to those obtained in the base case, while $k_{\mathrm{DES}}$ and $\Gamma$ are slightly overpredicted. This indicates that the effect of the artificial fluctuations is less marked than in the previous cases, affecting only the resolved portion of the turbulence.

The various turbulence-resolving behaviors presented in the discharge have effects downstream. Visualization of coherent structures in the flow is used to observe these effects. Figure 9 presents these coherent structures downstream of the discharge. The coherent structures are colored by the instantaneous pressure value and are generated as iso-surfaces of the so-called Q-criterion. For the base case (without AMR) presented in Figure 9a, some vortical structures can be visualized, which break into smaller structures at the end. Comparatively, for the $\nabla U$ criterion presented in Figure 9b, no such larger structures are observed, although the amount of smaller structures is qualitatively similar. For the pressure gradient criterion (Figure 9c), much larger structures are observed than in the base case, which is related to the stronger fluctuations (due to artificial effects already explained) and to higher rms values. Compared to the $\nabla U$ case, the $\nabla p$ case seems to present a smaller amount of smaller structures, which would explain why the $k_{\mathrm{LES}}$ values are not so different. Regarding the $\mathsf{\Omega }$ criterion, qualitatively, the behavior is similar to that of the base case, but the structures seem to be larger. Hence, the increase in the resolved fraction of turbulence would be due to the increase in the size of these structures due to artificial effects. From these results, it is concluded that abrupt transitions in refinement levels at the discharge are not desirable due to the negative effect downstream.

5.5. Effects of the AMR Approach on the Mean Axial Velocity Profiles and Turbulence Resolution in the Pilot-Size Spray Chamber

5.5.1. Mean Axial Velocity Profiles

Figure 10 presents the effects of the different AMR criteria on the radial profiles of the axial velocity $U_x$ in the pilot-size chamber. For comparison purposes, the experimental HWA data reported by Gutiérrez Suárez [27] are included. Overall, no significant differences are observed between the results predicted by the reference case without AMR (case C) and the cases using AMR (C1, C2, C3). It is observed that at lower axial positions (starting at $x/H=0.22$), case C1 ($\nabla U$ criterion) underpredicts the jet decay rate, while case 2 ($\nabla p$ criterion) predicts a higher decay rate than the reference case. This is an effect of the overprediction on the turbulent kinetic energy in the discharge. Compared to experimental measurements, the CFD results underpredict the jet decay and widening, which is evident from $x/R=0.147$. This behavior has been observed in computational studies of average velocities for this and other spray chambers and is explained in part by the mesh resolution requirements in the near-field of the annular jets.

Still, for the pilot-size spray chamber, these results are an improvement over those obtained in a previous study [18], as presented in Figure 11.

5.5.2. Turbulence Statistics

Figure 12 presents the radial profiles of the fraction of the resolved turbulent kinetic energy $\Gamma$ for the base case (case C) and AMR cases (C1, C2, and C3) obtained in different axial positions downstream of the jet discharge. In general, it is observed that the base grid (case C) presents the best overall results in terms of the portion of turbulence resolved. This result is expected since it has all elements refined to level 2 in a fixed grid. For all cases, close to the discharge ($x/H=0.073$ and $r/R$ = 0 to 0.15), the resolved fraction $\Gamma$ varies between 0.8 and 0.9. This stands in contrast to the results presented in Table 6, where much lower values were obtained. For all cases, this value becomes larger at lower axial positions, being close to 1.0. This is explained by the transition to larger turbulent structures, which are better captured in lower axial positions by maintaining grid resolution. In this zone, case C2 (criterion $\nabla p$) predicts a higher fraction of turbulence than the reference case. This is most likely due to the abrupt change in the level of refinement in the outflow duct, which generates some artificial oscillations. Outside the central jet ($r/R>0.2$), all criteria would be expected to have a lower resolved fraction than the base case (case C) since they do not explicitly refine the recirculation zone. However, the $\nabla p$ criterion predicts at various times higher resolved fraction values than the base case, indicating that the artificial fluctuations generated in the discharge also affect the recirculation zone. For this zone, the vorticity criterion offers the best results.

5.6. Effects of the Refinement Threshold and Cell Level

5.6.1. Pilot-Size Spray Chamber

Axial Velocity Profiles

Figure 13 presents the effects of using a lower threshold (case C2-1) and different refinement settings (cases C2-2 and C2-3) on the mean axial velocity $U_x$. For comparison purposes, data from cases C and C2 are included. Overall, cases C2-1 and C2-2 offer a similar behavior at all axial positions, predicting lower axial velocities than the base case (C). On the other hand, the C2-3 case predicts very similar results to those obtained by the base case. The apparent improvement in the results predicted by cases C2-1 and C2-2 is most likely due to incorrect turbulent kinetic energy production in the jet discharge (see Table 6), along with low refinement coverage in the jet (see Figure 7d), as discussed in the next section.

Turbulence Statistics

Figure 14 presents the radial profiles of the fraction of the resolved turbulent kinetic energy $\Gamma$ for cases C, C2, C2-1, C2-2, and C2-3. In the high jet velocity zone, cases C2-1 and C2-2 present a similar behavior, predicting lower resolved fraction values than the reference cases (C and C2). This is especially noticeable in the high shear stress zone that marks the periphery of the jet with its surroundings ($x/H=0.147,0.22$; $r/R\approx 0.25$), indicating a significative coarseness of the grid. On the other hand, in these regions, case 2-3 presents similar results to those of the base case. In the recirculation zone, which starts at $r/R\approx 0.25$ in the upper lines and $r/R\approx 0.45$ in the lower lines, different behaviors are presented. In this region, cases C2-2 and C2-3 tend to predict higher resolved fractions than that of the base case. This behavior is explained by Figure 15, showing that the increase of the resolved fraction in the recirculation region is related to a decrease in the turbulent kinetic energy $k_{\mathrm{RANS}}$ values. On the other hand, case C2-1 underpredicts the resolved fraction in all regions and axial positions. It is concluded that the discrepancies obtained are due to insufficient mesh coverage and abrupt transitions in the mesh refinements in the initial region of the jet, which are more evident for cases C2-1 and C2-2 than for C2-3. Moreover, turbulence resolution quality indicators should be analyzed comprehensively. For example, for C2-1 and C2-2, the maximum discrepancy with case C results in the resolved fraction of the near-jet region ($x/H=0.147,r/R<0.2$) being less than 10%, while for the modeled kinetic energy ($k_{\mathrm{RANS}}$), it is more than 250% (Figure 15).

5.6.2. Semi-Industrial-Size Spray Chamber

Axial Velocity Profiles

Figure 16 presents the radial profiles of the mean axial velocity at five different axial positions. Overall, the CFD results follow well the trend of experimental data. Compared with the experimental data, in the central jet, in the near-field region ($x=0.3\mathrm{m}),$ all cases predict well the characteristic shape of the annular jet, with higher velocities below the annular entrance and lower below the bulk-body. From $x/H=1.0$, all cases overpredict the experimentally measured mean velocity values. Although experimental measurements of mean velocity are unavailable at $x/H=1.4$ and $x/H=2.0$, the CFD predicted values appear similar to the maximum measured values, thus are assumed to be higher than the mean velocity. The differences in the results among different AMR thresholds are very small. It is observed that the S2 and S3 cases behave almost identically, while the S1 case predicts slightly lower decay rates in the central region. This is most evident at $x/H=0.536$ and $r/R=0$, where case S3 predicts a velocity 10% larger than cases S1 and S2. In the recirculation zone, the velocity predicted by case S3 tends to be lower than for cases S1 and S2. In summary, using a higher refinement threshold (case S3) reduces the decay rate of the main jet and generates lower velocities in the recirculation zone.

Figure 17 compares the results obtained from case S3 with the results obtained by other studies. Details of these studies can be consulted in Table 4. At $x=1.0$ m, the velocity profile predicted by case S3 is very similar to the one reported by the transient simulation of Saleh [39], which used three times more grid elements. At $x=2.0$ m, case S3 has a much better prediction of the radial velocity profiles than those reported by Anandharamakrishnan et al. [38] and by Huang et al. [40], which use a steady-state approach. In these studies, the asymmetry of the jet appears exaggerated, which might be related to a non-converged solution. A more extensive comparison with results from other studies was not performed since the vast majority of those do not use a three-dimensional transient approach. In addition, many papers do not report results at lower axial positions, limiting the discussion to upper lines ($x=0.3;0.6$), where discrepancies are less significant.

Velocity Intervals: Comparison with Experimental Data

Figure 18 presents the radial profiles of axial velocity intervals at different axial positions. In general, using lower $C_{\mathrm{crit}}$ threshold values increases the maximum axial velocity at most measurement positions. This is especially noticeable in the high-shear zone, where the jet interacts with its surroundings, and in the recirculation zone. In high-shear zones ($r/R\approx 0.15)$, the maximum velocity values predicted with case S3 are up to 80% higher than those of case S1. Nonetheless, around the jet and towards the walls, all cases follow well the pattern of the measured velocity intervals.

6. Conclusions

In this study, adaptive mesh refinement (AMR) is used in conjunction with a hybrid LES/RANS turbulence model to simulate the airflow in two different spray drying chambers. The use of AMR allowed a significant computational cost reduction. Compared to the base case (fixed grid with two refinement levels), the cases evaluated showed a 70% to 91% reduction in computational cost. This computational cost reduction becomes critical due to the long simulation times required to develop the internal flow ($t_c>2.5$). The AMR approach presented in this study can be used in parametric studies of spray drying processes and other confined jet applications where computational costs become a major limiting factor. In this way, it is possible to run numerous simulation cases using a high-end workstation. However, cost reduction strategies must be accompanied by a careful application of the AMR methods to avoid an impact on the accuracy of the results.

In general, the effect of different factors on mesh adaptation behavior must be assessed for all the criteria evaluated. In the absence of swirl at the discharge (as is the case in the semi-industrial-size dryer), the radial coverage of the refining zone is very different from that obtained in a low-swirl jet (pilot-size chamber). The size of the background grid (level 0) is also influential, since if it is excessively coarse, a sufficiently high relative error will not be generated to trigger the mesh adaptation. Finally, the refinement threshold must also be balanced against the maximum refinement level allowed. Based on the results of this study, it does not always seem appropriate to reduce refinement coverage by a higher level of refinement. It is estimated that the ideal refinement threshold should generate a grid containing between 25% and 40% of the elements of a fixed grid. The use of lower threshold values, such as that of case S3, which generates a grid with a proportion greater than 40%, does not seem to bring a noticeable improvement in the results.

In detail, the differences in results when using various AMR criteria ($\nabla U$, $\nabla p$, $\mathsf{\Omega }$) are small but significant. At first glance, the pressure gradient-based criterion allows for a better modeling of the jet discharge by ensuring its complete refinement. However, the transition between different refinement levels is abrupt, generating artificial oscillations. Moreover, as it focuses solely on a narrow region, at the jet discharge, it inaccurately models the kinetic energy $k_{\mathrm{RANS}}$ upstream. Both the $\nabla U$ and $\mathsf{\Omega }$ criteria generate a smoother transition between refinement levels in the inlet duct but are not consistent in refining to the highest level of the section with the highest axial velocity in the discharge. Due to these sporadic changes in the level of refinement, some artificial fluctuations are generated but of lesser intensity than in the case of the $\nabla p$ criterion. To eliminate or reduce these fluctuations, there are several alternatives: The first one is to carry out a fixed prerefinement in the cells around the inlet diffuser. In this way, AMR-controlled unrefinement would not be allowed in such cells, and these artificial fluctuations would be eliminated. According to the results of Mozaffari et al. [24], it is recommended to extend the prerefining to the recirculation zone under the bluff-body to remove truncation errors due to the high flow instability; the second is to increase the number of buffer cells between cells with different refinement levels. In this way, the transition between cell levels would be smoothed, and the results of the $\nabla p$ criterion would be enhanced; the third alternative is to combine the $\nabla p$ criterion with one of the other two to generate a more homogeneous refinement in the discharge and in the near-region of the jet; the fourth is incorporating nonlinear terms in the relative error calculation, such as $C_{\mathrm{rel}}=\sqrt[3]{c_{\mathrm{rel}}Re}$. This would smooth the transition between refinement levels and increase the coverage of transition cells.

Downstream of the discharge, the cases with higher jet refinement coverage (C1–C3 and S1–S3) manage to capture by dynamic grid adaptation the transient motion of the jet. The operational differences between the AMR methods are mainly observed in their interaction with the walls. The criterion $\nabla p$ does not perform refinements on the chamber walls (except for the inlet and outlet ducts). The $\nabla U$ and $\mathsf{\Omega }$ criteria do perform refinements on the walls. This can be an advantage in improving the accuracy of the results where there is impingement between the gas or particles with the chamber walls or internal ducts, as in the case of the semi-industrial-size dryer. Regarding the axial velocity profiles, all cases predict good results compared to experimental results and those obtained in previous studies. However, it is necessary to properly adjust the AMR method to reduce the generation of artificial oscillations in the jet discharge, which can potentially increase the turbulent kinetic energy dissipation rate and artificially increase the jet decay rate. To better assess this impact and improve the applicability of AMR methods along turbulence-resolving models in spray drying simulations, it is desirable that future work include experimental measurements of turbulence inside the chamber.