1D Modeling Approach for Heat Transfer in Packed Beds with Embedded Heat Sources

: To improve heat transfer in packed bed temperature swing adsorption processes for direct air capture, reactors with embedded cylindrical heating pipes have been developed. Optimization of these inherently dynamic systems currently requires a transient two-dimensional (2D) model, which is computationally very expensive. In the present work, a method is therefore developed to translate a 2D fixed bed geometry into an equivalent one-dimensional (1D) model by placing line heat sources along the direction modeled in 1D. To determine these sources’ strength, a Nusselt correlation is required. It is found that for staggered cylinder configurations, the single cylinder correlation works well. For in-line configurations, an analytical correction factor is successfully developed to account for the effect of the thermal wake of the upstream cylinders on the heat transfer around a cylinder. The 2D to 1D translation approach is then tested with three different cylinder-packing geometries at varying Pe ́ clet numbers and for steady-state and dynamic simulations. For the steady-state simulations, the 1D model has a maximum deviation of 10% in the bed mean temperature and for the outlet temperature from the 2D results (scaled to the maximum temperature difference), thus showing good agreement. For the dynamic simulations, the deviation is below 20% for most conditions, showing reasonably good agreement. The merit of the translation approach becomes apparent when looking at the computational time: the 1D model calculations are a factor of 500 faster than the 2D calculations.


INTRODUCTION
Packed beds are encountered in many chemical engineering applications and have been used for a long time in the field.One application of packed beds that has become very relevant recently because of climate change mitigation is their use in temperature swing adsorption processes for direct air capture (DAC).During the regeneration of the sorbents used in such processes, heat must be supplied to the sorbent particles, and as such, heat transfer is important in these systems.To increase the heat transfer rate during desorption in fixed bed reactors, designs with embedded cylindrical heating pipes have been developed, for example, by Bajamundi et al. 1 and Schellevis et al. 2 The heating pipes in these designs are evenly spaced in the direction perpendicular to the flow.In these systems, coupled transport of momentum, heat, and mass takes place.Because of the geometry, these transport processes are clearly two-dimensional (2D) phenomena.Furthermore, fixed bed adsorption systems are inherently transient.These factors combined make modeling fixed bed adsorption systems numerically challenging and computationally intensive.
For design, scaling up, and optimization purposes under varying process conditions, it is desired to have a reasonably fast model of such adsorption systems so that design calculations and sensitivity analyses can be done quickly and without too much computational effort.A one-dimensional (1D) model would fit these purposes.However, in a 1D model, it is not trivial to describe the heat transfer phenomena induced by flow around an arrangement of pipes. 3In this article, a method for translating 2D geometries with evenly spaced heating pipes in a packed bed to a representative 1D model is therefore developed.Such a 1D model has the benefits of much faster calculation times and easier implementation (especially of coupled transport processes), while still accounting for the 2D nature of the heat transfer.A similar dimension reduction approach is proposed for heat transfer inside particles in a fluidized bed by Luo et al. 4 In this work, 1D intraparticle heat transfer limitations are captured in a zero-dimensional (0D) particle model by use of a modified heat transfer coefficient correlation.This approach, however, is not directly applicable to the problem in the present work as the system is very different.
Before we go into the approach of describing 2D heat transfer in a 1D model, some definitions and previous work regarding heat transfer in fixed beds are discussed.There are two relevant phases in a fixed bed system: the solid phase (sorbent particles) and the fluid phase flowing between the particles.In the solid phase, heat can be transferred via conduction, while convection is the most significant transport mechanism in the fluid phase.The temperatures of the fluid and solid phases are not necessarily the same: if the heat transfer resistance between the two phases is significant, a temperature difference between them will arise.Frequently, however, heat transfer between the fluid and particles is fast enough, so the temperatures of the two phases can be assumed equal.This assumption is often made and is called the local thermal equilibrium (LTE) assumption.In this work, the LTE assumption is also made, so the fluid and particles in the fixed bed can be described as a single phase. 5nother important aspect of the packed bed flow is the flow regime. 5A key dimensionless parameter governing this is the Darcy number.For low Darcy numbers (Da < 10 −3 ), the flow can be accurately described by Darcy's law, neglecting any inertial effects.For higher Darcy numbers, these inertial effects must be included by including, for example, the Forcheimer extension or the Ergun model. 5In this work, only the Darcy flow is considered, but the analysis presented in this article is likely to hold in non-Darcy flow regimes as well, if the appropriate Nusselt correlations are used and the additional dispersion effects are included.
Initially, the literature is reviewed to assess if an appropriate Nusselt correlation for heat transfer around banks of cylinders in a packed bed is available.Heat transfer around single cylinders and cylinder banks in packed beds has been investigated experimentally, analytically, and numerically by a number of studies in the literature.Cheng 6,7 derived an analytical relation (eq 1) for the Nusselt number around a single cylinder in porous medium under steady Darcy flow conditions at high Pećlet numbers Pe ( 1 ) D using the boundary layer approximation and with the LTE assumption.The Pećlet number is, in this case, defined as the rate of heat transfer by convection over the heat transfer rate by conduction. 5,8Care has to be taken that the effective thermal diffusivity (averaged between the fluid and particles) is used. 5,6Alternative steady-state analytical correlations using different assumptions were later derived by Shigeo, 9 Romero, 10 Pop and Yan, 11 and Magyari and Keller. 12Nasr et al. 13 experimentally investigated heat transfer around a single cylinder in a packed bed and found a reasonable agreement between experimental results and the analytical correlation of Cheng. 7 Pe 1.0157 A significant difference between the above-mentioned steadystate situation and the adsorption-based processes targeted here is the changing temperature field in the bed as well as the flow variation in the different phases of the adsorption/desorption cycle.Analytical relations for transient forced convection over single cylinders in porous medium are therefore also investigated.Shigeo 14 developed a relatively simple expression by dividing the time domain in a conduction-dominated and convection-dominated part.Sano 15 Thevenin 16 has numerically analyzed transient forced convection heat transfer from a cylinder in a packed bed using the Darcy−Forcheimer model and found that the time to reach steady-state slightly depends on the Darcy number.Al-Sumaily et al. 17−19 conducted a series of analyses on steady-state and transient heat transfer from a single cylinder in a packed bed, modeling the fluid and solid phase separately with a local thermal nonequilibrium (LTNE) model.They found that their LTNE model described the experimental results of Nasr et al. 13 better than the analytical relation of Cheng. 7In conclusion, the above-cited literature results indicate that a fair description of heat transfer around a single cylinder in a fixed bed is possible.For banks of cylinders in a fixed bed, however, the situation is different.
Banks of cylindrical tubes in free flow (no particles around them) have been investigated in a number of studies.Z ̌ukauskas 20 gives an overview of experimental studies of heat transfer from cylinders and cylinder banks in the absence of porous media.They found that for both in-line and staggered banks of cylinders, the overall Nusselt number increases with the number of rows toward an asymptotic maximum value.Heat transfer in the first row is the poorest and similar to a single cylinder at low Reynolds numbers.For subsequent rows, the heat transfer increases due to the turbulence induced by the upstream rows.Gnielinski found similar results. 3,21Both Z ̌ukauskas and Gnielinski developed correlations to calculate the Nusselt number for a bank of cylinders in free flow.Fowler and Bejan 22 found that for banks of cylinders in the absence of porous media at very low Reynolds numbers (1−30), the opposite trend of what Z ̌ukauskas and Gnielinski observed occurs.In this case, the Nusselt number of the first row is highest, and it decays to an asymptotic value.This effect is attributed to flow development in the first rows of the bank.
For small in-line and staggered configurations of cylinders embedded in porous medium, only a few studies have been done.Layeghi 23 numerically analyzed forced convection from a staggered cylinder bank embedded in wooden porous medium.They found that the Nusselt number for the first row is highest and then decreases for subsequent rows, similar to the results of Fowler and Bejan. 22Al-Sumaily 24 numerically analyzed four cylinders in a staggered and in-line configuration in a fixed bed, using a local thermal nonequilibrium model.It was found that for in-line cylinders, the Nusselt number for the downstream cylinders is significantly decreased compared to the first row of cylinders, while for staggered cylinders, this effect is less pronounced.The difference in the Nusselt number decreases with increasing spacing between the cylinders.Sayehvand et al. 25 analyzed two in-line cylinders in a packed bed and found that the

Industrial & Engineering Chemistry Research
overall Nusselt number for the two cylinders is significantly lower than the single cylinder Nusselt number from Cheng's correlation. 7n summary, many studies have investigated heat transfer around single and multiple cylinders in a fixed bed.For a single cylinder, easy-to-use correlations to calculate the Nusselt number (both steady-state and dynamic) exist.For banks of cylinders in fixed beds, no universally applicable correlation for the Nusselt number seems to be available due to the interactions of the cylinders in such a configuration.Especially for in-line configurations, heat transfer seems to decrease for downstream rows.Also, for low Pećlet numbers (where conduction becomes significant), no satisfactory correlation was found.Finally, dynamic heat transfer around the banks of cylinders in a fixed bed is hardly investigated.For the 1D modeling approach in the present work, having a sufficiently accurate Nusselt number correlation for a cylinder bank configuration is a necessity, and therefore, such a correlation will be developed in Section 2.2.
Based on the literature review, we conclude that there is no satisfactory way of translating the 2D heat transfer problem at hand into a 1D model yet.In this work, such an approach (which Figure 1 schematically shows) is therefore developed.Figure 1 shows that from a given 2D geometry, a representative 1D model will be constructed with heat sources distributed over the axial direction.Whereas in the 2D model at the tube walls, the temperature can be set to a constant value, in the 1D model, this does not work.Therefore, in the 1D model, a heat transfer coefficient or Nusselt number is required that accounts for the heat transfer resistance from the tube walls to the fluid in the bed.Finding a suitable correlation for this from 2D (or 3D) simulations as reliable input for the 1D modeling is a main objective.In this work, first, a literature review about heat transfer from cylinders in a fixed bed is therefore done.Subsequently, a 2D model is constructed.Third, the translation approach from a 2D model to a representative 1D model is discussed.It will be shown that for in-line cylinders, no satisfactory Nusselt correlation is available in the literature.Therefore, such a correlation based on an analytical approach is developed as the next step.Finally, the results of the Nusselt number calculation and the overall 2D to 1D translation approach are discussed.

MODELING METHODS
In the present work, two heat transfer models are created: a 2D model using commercial finite element software (Section 2.3) and a 1D model in Python (Section 2.4).Both models are able to describe the steady-state as well as the transient behavior of the system.The 2D model is used to model heat transfer around banks of tubes in a packed bed under varying conditions and in different geometries (Section 2.1).A geometry-based procedure is then used to translate the 2D geometry into an equivalent 1D model.Then, the 1D model is run under the same conditions as the 2D model, and the heat transfer results are compared using several key performance indicators.As mentioned in the Introduction, for the 1D model, a correlation for the Nusselt number as a function of the Pećlet number is required.Such a correlation is developed in Section 2.2 and tested against a 2D model describing heat transfer around cylinders in series in a channel.
2.1.Geometries.Three different geometric cases are considered in this work, which are displayed in Figures 2−4.
In case 1, the tubes are placed on a rectangular grid, in line with the flow direction.In case 2, the pipes are placed in a staggered configuration.Case 3 places the pipes in a part staggered, part inline configuration, where the free flow path between the pipes is minimized.In a follow-up study, more different geometries (including combinations of the three cases in this work) could be investigated.

Nusselt Number Calculation.
The Nusselt number that is required in the 1D model can, of course, be calculated from the 2D simulation results.It is, however, desirable to be able to calculate the Nusselt number via correlations, so that no computationally expensive 2D simulations are required at all.It will be shown later that under a large range of conditions, the single-cylinder correlation by Cheng 7 (see eq 1) gives acceptable    Industrial & Engineering Chemistry Research results.However, for low Pećlet numbers and especially for inline cylinder configurations, the correlation is considerably off, and for these regions, alternative correlations are developed in this section.

Low Pećlet Number.
For Pećlet numbers smaller than or close to unity, heat transfer around the cylinder is dominated by conduction in the direction perpendicular to the cylinder wall.This contradicts the assumption of a thin boundary layer by Cheng, 7 making their correlation inapplicable to this situation.A conduction-corrected correlation for a cylinder, as part of a bank of cylinders in a fixed bed, is therefore derived here.In the limiting case of a Pećlet number of zero, with the regular assumptions (T → T ∞ when r → ∞), there exists no steady-state solution for this system.However, assuming that at a certain distance s from the cylinder wall (which depends on the tube bank geometry), the temperature equal to T ∞ allows solving the steady-state conduction heat balance by integrating Fourier's law.This leads to the following expression for the Nusselt number for a distance s from a cylinder with a diameter D ( ) Defining the dimensionless distance s̃= s/D, this expression can be rewritten as Analogous to the well-known Nusselt correlations for heat transfer around a sphere, eqs 1 and 4 can be combined to yield the following overall expression for conduction and convection around a cylinder Equation 5 is later found to give acceptable results for the Nusselt number over the entire range of Pećlet numbers in staggered banks of cylinders.
2.2.2.In-Line Cylinder Configurations.It will turn out later (see Section 3.1) that for geometric cases 2 and 3 (staggered cylinder configurations), eq 5 gives acceptable results for the Nusselt number over the whole range of Pećlet numbers.For inline cylinders (such as in geometric case 1), the correlation is less accurate, due to the influence of the thermal wake of the upstream cylinders.In order to have a more accurate description of heat transfer around cylinders placed in series in a packed bed, an analytical approach is worked out in this section.The main premise for the approach in this work is that the decrease in heat transfer from subsequent cylinders is not a hydrodynamic effect but a thermal driving force one.Hence, it is assumed that the Nusselt number is essentially unchanged for subsequent cylinders, but the driving force for heat transfer is distorted due to the thermal wake of the upstream cylinders.Note that the present situation in a fixed bed is different from the free flow around in-line cylinders: in the latter case, the flow profile after a cylinder is distorted (e.g., a vortex sheet will form), influencing the heat transfer around downstream cylinders.In a fixed bed, however, the hydrodynamics are governed mostly by the interaction of the fluid with the particles of the fixed bed, yielding a very similar flow profile around all cylinders.Therefore, the flow profile hardly varies for subsequent downstream cylinders.A sketch of the system is shown in Figure 5.
At reasonably high Pećlet numbers (Pe D ≫ 1), the driving force for heat transfer from a single cylinder in a packed bed is only influenced by its wall temperature and the upstream fluid temperature.Thus, in this case, the driving force is ΔT = T wall − T̅ mix.cup,before , where T̅ mix.cup,before is the mixing-cup averaged fluid temperature just upstream of the cylinder.However, for multiple cylinders in series, the actual driving force of cylinder n will be influenced by the wake of the upstream n 1 cylinders.The actual driving force for heat transfer from cylinder n will therefore be ΔT n = T wall − T̅ wake,n−1 , where T̅ wake.,n−1 is the average temperature of the fluid in the wake of the cylinders upstream of n.Based on these observations, two driving force definitions are introduced Here, ΔT obs.,n is the observed driving force for heat transfer from cylinder n based on the fluid mixing-cup average temperature over the entire cross-section of the packed bed and hence dependent on W (see Figure 5), while ΔT act.,n is the actual driving force cylinder n experiences.Note that when the distance between subsequent tubes tends to infinity, both driving forces become equal as the wake disappears.When doing a one-dimensional heat transfer simulation, ΔT obs.,n is always known, while ΔT act.,n is unknown.Thus, the goal at this point is to find an expression that relates the actual driving force to the observed driving force.Then, with this driving force, the heat flux from a cylinder can be calculated with To derive the aforementioned relation between the actual and observed driving forces, an equation for the wake temperature after a series of cylinders is required.If it is assumed that the heating pipe can be considered a line source buried in the fixed bed, the following expression is obtained for the thermal wake after the heat source (assuming Pe D ≫ 1) 5 Here, T ∞ is the fluid temperature upstream of the heat source, q′ is the strength of the heat source (the amount of heat transferred from the first cylinder to the fluid) in W m −1 , x is the distance from the heat source, parallel to the flow, and y is the distance from the heat source perpendicular to the flow.To ease the analysis, the following dimensionless quantities are introduced Using these dimensionless variables, eq 10 is cast in a dimensionless form i k j j j j j y This temperature distribution that is formed after the first cylinder affects the driving force for heat transfer around the second cylinder.Since the hydrodynamic conditions around both cylinders are the same as those argued above, the Nusselt number is the same as well.However, because of the wake of the first cylinder, the driving force for heat transfer around the second cylinder is different.To quantify this effect, eq 15 is integrated over a certain dimensionless vertical length Y ̃to get the average wake temperature i k j j j j j y Two variables in eq 16 still require a value: the horizontal distance x̃and the vertical distance to average over Y ̃.For x, the most logical choice is the dimensionless pitch pt (heart-to-heart distance) of two subsequent cylinders.It makes sense that Y̅ must be proportional to the thickness of the thermal boundary layer around the cylinder of interest.For flow around a cylinder in a packed bed, the boundary layer thickness scales according to Taking C = π/2 (so that the boundary layer thickness becomes proportional to the swept arc length D

2
) was found to give the best results. 3,7,26Substituting the relation for Y ̃in eq 16 gives i k j j j j j j j y This relation for the wake temperature holds for a system with two cylinders.To describe a system with more cylinders in series, it is assumed that the wakes of subsequent cylinders are superimposed on each other so that their contributions add.In the equation form, this means that for the wake temperature after the nth cylinder, the following relation holds i k j j j j j j j j y { z z z z z z z z Pe q n i p To relate the average wake temperature and the mixing-cup average temperature, we derived an expression for the mixingcup average fluid temperature after n cylinders.For the system at hand (a fluid flowing through a channel while being heated by several point sources in a row), the following expression is trivially derived (assuming again that Pe D ≫ 1) In this expression, W ̃is the dimensionless width of the channel (W/D) and q̃i is the dimensionless heat flux coming from the ith cylinder.
Now that expressions for the wake temperature and the mixing-cup average temperature have been derived, the goal of finding a relation between the observed temperature difference and the actual (wake-induced) temperature difference is revisited.In the dimensionless form, a relation between ΔΘ act.,nand ΔΘ obs.,n is sought, where 1 Since ΔΘ obs.,n is known when doing calculations in a 1D model, but ΔΘ act.,n is required to calculate the heat flux from the nth cylinder, a correction factor is derived to relate the two quantities.For this correction factor, two approaches are investigated.The most straightforward approach is to relate the temperature differences directly, so that Another, an equally valid approach to define a correction factor is The temperature difference correction factor ϕ and the temperature correction factor ψ are related as follows For both correction factors, an analytical expression can now be derived, based on their definitions, eqs 17 and 19.For the temperature difference correction factor, this yields i k j j j j y 1

Industrial & Engineering Chemistry Research
This equation is inconvenient in its present form, as it requires the heat fluxes of all sources up to cylinder n − 1 to calculate the correction factor for cylinder n.Replacing the fluxes with their dimensionless expressions qñ Equation 28 still requires a lot of upstream data (the upstream correction factors and mixing-cup average temperatures) to calculate the temperature difference correction factor for cylinder n and is therefore still inconvenient for practical use.It is not straightforward to simplify eq 28 even further; therefore, the attention is turned to the temperature correction factor ψ instead.Taking the definition of ψ in eq 24 and substituting eqs 17 and 19 gives i k j j j j y This equation also contains the fluxes around the upstream cylinders, but a simple assumption about the flux allows a large simplification of the relation.The flux for subsequent cylinders will decay because the driving force decreases.This is due to both the increasing influence of the wake of the upstream cylinders and the mixing-cup average fluid temperature increases.Therefore, it is assumed that the flux around cylinder i follows the proportionality relation in eq 31 In eq 31, q i is the mean flux over cylinders 1 to n.This proportionality approximation is compared to 2D simulation data of ten cylinders in a row in Figure 6 and works reasonably well.Substituting qĩ in eq 30 with the approximation and canceling the mean flux terms yield i k j j j j y Equation 32 is a relation that does not depend on any upstream fluxes or temperatures and therefore allows relatively easy calculation of the temperature correction factor.Additionally, the equation is based on an analytical derivation and, in addition to the flux approximation in eq 31, does not contain any additional fit parameters.It is good to note that at very low Pećlet numbers, eq 32 can give ψ values lower than unity, suggesting enhancement of heat transfer due to the thermal wake.This is obviously unrealistic, and therefore, when ψ is calculated lower than unity, a value of unity should be assumed.Equation 32 is expected to perform well for medium to high Pećlet numbers (10 or larger) based on the assumptions made during its derivation.For low Pećlet numbers, the agreement is likely less good, but for steady-state 1D simulations, the correction factor is less important under these conditions as heat transfer is dominated by conduction then.It is also expected that at very low pitches, the correlation performs worse, as the assumption of treating the cylinders as line sources is not valid in that case due to overlapping boundary layers.

2D Model.
For the 2D model, COMSOL Multiphysics version 5.5 is used to solve the mass, momentum, and energy conservation equations.The 2D energy balance is given by eq 33.The momentum and continuity equations are also solved but not shown here.For the fluid flow, the "free and porous media flow" model is used.The energy balance is solved using the 'Heat Transfer in Porous Media' module.All physical properties are made independent of the temperature and pressure.During the solution of the model, first, the flow field is solved for a steadystate situation.Then, the energy balance is solved using a transient solver using the steady-state flow field obtained in the first step.
As boundary conditions for the 2D problem, the following expressions are used T T at cylinder walls: In these expressions, n⃗ is the unit of the outward normal vector to the boundary.
To compare the 1D model performance to the 2D results, we required to calculate the Nusselt number from the 2D simulations.This number is calculated via eq 37, where the logarithmic mean temperature difference between the fluid in the bed and the cylinder walls is used as a driving force.In eq 37, n cyl is the number of heating tubes in the 2D domain, ΔT lm is the logarithmic mean driving force (see eq 38), q h ″ is the heat flux normal to the cylinder walls, and s is the coordinate along the circumference of a pipe.This Nusselt number is then inserted as input in the 1D model.2.4.Translation to 1D.In the 2D geometry, every point inside the sorbent bed is located at a certain distance from the nearest heating pipe wall.Heat from the pipe walls, which is mainly supplied via conduction in the direction perpendicular to the walls, has to be transported over this distance.Thus, a key parameter in heat transfer from the pipe wall to the bed is the mean heat transfer distance.If the geometry is simplified to a 1D model, this heat transfer distance must be accounted for.However, in a 1D model, there are no pipes and no gradients perpendicular to the net flow direction.The only way, therefore, to include the heat transfer distance is to place point heat sources distributed over the length of the 1D domain.These sources are then separated by twice the mean heat transfer distance calculated from the 2D geometry.
To calculate the mean heat transfer distance, the 2D geometry is taken, and a large number of points is evenly distributed over the domain.The points that lie within the heating pipes are excluded from the analysis.Then, for all points, the distance to the nearest cylinder wall is calculated and the average is taken to be d̅ .To assess the distribution of the heat transfer distance throughout the domain, a distribution plot is also made (see Figure S1 in the Supporting Information).Finally, the mean heat transfer distance is calculated separately for the points on the inlet and outlet boundaries to yield the values d̅ inlet and d̅ outlet , respectively.
Another aspect of the 2D to 1D translation is the determination of the 1D domain length.It would be easiest to select the length of the 2D packed bed in the flow direction as the length of the 1D domain.However, as shown in Figure 7, this will not take into account the effect of the volume occupied by heating pipes on the fluid residence time (which should obviously be the same in 2D and 1D).Furthermore, the total masses of gas and particles in the two situations must be equal.To achieve the latter, first, the heating tubes are removed from the 2D geometry by defining the tube fraction ε tube as in eq 39 and then calculating the modified domain height and width via eq 40 and eq 41.This process is schematically shown in Figure 7.To achieve an equal residence time in the 1D and 2D geometries, the 1D superficial velocity is defined as in eq 42, where u in, 2D is the superficial inlet velocity in the 2D system.
Next, the heat sources must be placed in the 1D domain.Ideally, the first heat source is placed d̅ inlet from the inlet, the last heat source is placed d̅ outlet from the outlet, and the other heat sources are separated by 2 d̅ .However, in most cases, this arrangement will not fit exactly in the domain length L 1D .Thus, a geometric factor (f geom , see eq 44) is introduced, by which the heat transfer distances are multiplied to fit the sources in the domain.To calculate this factor, the optimum number of sources is required and is calculated with eq 43.The final placement of the heat sources is shown schematically in Figure 8. Figure 9 schematically shows the entire translation approach.

Industrial & Engineering Chemistry Research
balance equation for this case is shown in eq 45.In this equation, besides the transient, diffusion, and convection terms, a term related to the heat sources is included.In this term, δ s (x) is a Dirac-delta function, which has the value of one at the location of a heat source and zero otherwise.a tube is the specific surface area of the tubes, defined in eq 46.The heat transfer coefficient h t,tube can be calculated from the 2D simulations or predicted  For ease of analysis, the energy balance is nondimensionalized, as shown in eq 47.As the characteristic length, the cylinder diameter is chosen.Furthermore, Θ is the dimensionless temperature (defined in Table 1), τ D is the dimensionless time (defined in Table 1), x̃is the dimensionless axial coordinate, Pe D is the Pećlet number based on the effective thermal conductivity (defined in Table 1), at ube is the dimensionless specific surface area (defined in eq 51), and Nu D is the Nusselt number.The dimensionless heat balance is solved numerically, using a finite-volume discretization in space and an implicit numerical method in time.The spatial grid is shown schematically in Figure 10.The domain is divided into N = 1000 cells with a constant width of Δx.The discretized flux equation is shown in eq 48, with the semidiscretized heat balance shown in eq 49.δ n (x) in eq 49 equals zero in the cells without a heat source and unity in cells with a heat source.In these cells, the discretized heat source as shown in eq 50 is added, taking care that the total heat source is not affected by the cell size.
The system of ODEs in time resulting from the semidiscretized heat balance (eq 49) is solved with the IDA solver from the SUNDIALS suite, 27 using the Assimulo 28 and NumPy 29 packages with Python.The IDA solver uses an implicit, variable-order backward differentiation formula (BDF) method.The results of both the 2D and 1D simulations were processed and visualized using the NumPy and Matplotlib Python packages.

RESULTS AND DISCUSSION
The results are discussed in four parts.First, the calculation of the Nusselt number using correlations for staggered and in-line cylinder arrangements is discussed.Second, the heat transfer distance distribution in the 2D geometries is discussed.Third, steady-state simulations are done for the three geometries, and the 2D simulation results are compared to the 1D model results.Finally, this is repeated for transient simulations.For both the steady-state and transient simulations, the Nusselt number used in the 1D model was calculated from the 2D simulation results and the correlations derived in this work.

Nusselt Number Calculation.
First, for the three geometric cases, the Nusselt numbers calculated from the 2D simulations are compared to the analytical Nusselt correlation given by eq 5 as a function of the Pećlet number.The Nusselt numbers from the 2D simulations are calculated in eq 37. Figure 11 shows the comparison.It is clear that for cases 2 and 3 (staggered geometries), the single cylinder correlation gives very acceptable results (a deviation of 15% or less).For case 1 (in-line geometry), however, the correlation is considerably off.
To analyze why the correlation is off for cylinders in series, a simple geometry of five cylinders in series in a fixed bed is modeled (see Figure 5).Then, the dimensionless heat flux distribution around the cylinders as a function of the angle ϑ is calculated from the simulation results and plotted in Figure 12.The calculated distributions are compared with the analytical  results of Cheng, 6 shown by the solid, black line in Figure 12.It is immediately clear that the flux profile the first cylinder differs from the profiles around the other cylinders.Furthermore, the profile around the first cylinder closely resembles that of the analytical solution for a single cylinder in a packed bed.This shows that heat transfer around the first cylinder is hardly affected by the downstream cylinders.However, it is also clear that the heat transfer around the second to the fifth cylinders is affected significantly by the upstream cylinders.This is because the temperature profile of the fluid flowing toward these cylinders is very different from that of the first cylinder.This explains the different Nusselt numbers encountered for geometric case 1.In fact, heat transfer around downstream cylinders decreases with the position in the array, just as Fowler et al. 22 found for banks of cylinders at very low Reynolds numbers.
Next, the analytical correction factor approach for the cylinders in series is validated.To test eq 32, 2D simulations are done with ten cylinders in series in a channel with a dimensionless width of 20 (see Figure 5 for a schematic view of the system).As boundary conditions, eqs 34−36 were used.The dimensionless cylinder pitch was varied between 1.1 and 16, and the Pećlet number was varied between 1 and 1000.From the 2D simulations, the temperature correction factor ψ was calculated and compared to the analytical results from eq 32.
Figure 13 shows a comparison of the analytical temperature correction factor (presented in eq 32) and the factor calculated from the 2D simulation results as a function of the Pećlet number and the dimensionless pitch for the second cylinder in the row.First of all, the correction factor increases with increasing Pećlet number, which is expected as at constant pitch, the thermal wake will be narrower at a higher Pećlet number.This impedes heat transfer, thus explaining the higher temperature correction factor.Second, the correction factor decays as the pitch increases.This is also in line with expectations and the work of Al-Sumaily 24 as at a higher pitch, the thermal wake of a cylinder has more time to dissipate, and thus, the influence on downstream cylinders will be less severe.Overall, the analytically derived temperature correction factor matches the numerical correction factor calculated from the 2D simulations very well under all conditions.This is remarkable given the assumptions made during the analytical derivation.
To further assess the performance of eq 32, a parity plot of the analytical ψ values versus the values calculated from 2D simulations for different Pećlet numbers, dimensionless pitches, and cylinders in the row is shown in Figure 14.The color in the plot shows the Pećlet number corresponding to the correction factor.It is clear that under most conditions, the analytical temperature correction factor ψ (being the difference in the observed driving force versus the actual driving force for heat transfer) is well within 10% of the correction factor calculated from the 2D simulations.Only at very low Pećlet numbers (Pe D < 5), the agreement is less good, which is expected based on the assumptions done during the derivation of eq 32.Since the correction factor also performs reasonably well at low Pećlet numbers, it will be used in the 1D model under these conditions as well.As a future improvement, a more accurate correlation for the correction factor could be developed for the low Pećlet (conduction-dominated) situation.
The success of this strategy is also clearly illustrated when looking at the results for the 1D versus 2D simulation comparison in Figures 17 and S-5 in the Supporting Information, where the lines marked with "corr." in the legend show the results from the 1D model, where the Nusselt number from eq 5 is used combined with the analytical correction factor from eq 32.The agreement of the 1D model with the 2D model is even slightly better with the purely analytical Nusselt/ correction factor approach as compared to using the Nusselt number calculated from the 2D simulations in the 1D model.This is because when using the correlation, a different Nusselt number is used for the different rows (due to the correction   Industrial & Engineering Chemistry Research factor), while the average Nusselt number from the 2D simulations is taken as constant.This shows that the Nusselt number for use in the 1D model can be calculated using analytical correlations, without having to do 2D simulations under a wide range of conditions.

Steady-State Simulations.
Steady-state simulations are done with the 2D model for the three geometric cases, all with varying Pećlet numbers in the range of 1−1000.1D simulations are done with both the Nusselt number calculated from the 2D simulations (eq 37) and the analytical Nusselt numbers calculated from eq 5 (cases 2 and 3) and eq 5 in combination with eq 32 (case 1).Figures 2, 3, and 4 show the 2D temperature profiles for a Pećlet number of 50 for the three geometric cases.The temperature profiles are very similar to those for Stokes flow around a bank of cylinders without a packed bed, as is expected from the literature. 5This is because the packed bed suppresses the formation of a vortex street behind the cylinders that would be expected for the laminar free flow around a cylinder.
To compare the temperature profiles from the 2D and 1D simulations, the 2D results are collapsed along the x direction to yield an average axial profile.This is done with eq 52, which calculates the mixing-cup average temperature as a function of the axial coordinate.In the equation, z is the axial coordinate.Figure 15 shows the collapsed 2D temperature profiles (solid lines) and the 1D temperature profiles (dashed lines) over the dimensionless axial coordinate for case 1 at three different Pećlet numbers.In the 1D lines, the location of the heat sources is clearly visible, and the 1D temperature profiles are slightly below the 2D profiles, but otherwise, the agreement between the 1D and 2D results is good for all three Pećlet numbers.The same plot is made for case 2 in Figure 16, with very similar results.For case 3, the results are very similar to those for case 2, so the plot is not shown here.The underprediction of the temperature profiles in the 1D model could be caused by the difference in "heated length" between the 1D and 2D models in the sense that in the 1D model, the cylinder is reduced to a point source (the "kinks" in Figure 15), whereas in the 2D model, there is a certain length in the x direction that is heated by the cylinders.Another cause could be the distance of the inlet to the first row of cylinders, which is different in the 2D and 1D models due to geometrical The accuracy of the 1D model could possibly be improved by improving the estimation of the inlet distance and by replacing the point sources with line segments in the axial direction to smooth the temperature profiles.
To make a more complete comparison between the 2D and 1D steady-state results, Figures 17−19 show the mean and outlet dimensionless temperatures for cases 1−3, respectively, and Figures S-5−S-7 in the Supporting Information show the relative error in the 1D model results.In the 1D simulations denoted by "corr", the analytical expression for the Nusselt number was used; for the simulations marked with "sim", the Nusselt number calculated from 2D simulations was used.In all cases, the 1D model predictions are slightly below the 2D results, but the error is within 10% for the mean and outlet temperatures.The calculation time of the 2D results was 110 s on a 24-core AMD Ryzen Threadripper 3960X CPU using all cores, whereas the 1D results were calculated in only 5 s on a single core of an AMD Ryzen 5 3600 CPU.Accounting for the number of cores used, this means that the 1D model is about 500 times faster than the 2D model, while still having a maximum deviation of only 10%, demonstrating the benefit of the translation approach.Comparison of steady-state cup-mixed average temperatures of the 2D and 1D simulations as a function of the Pećlet number for case 1, using the Nusselt number calculated from 2D simulations (denoted by "sim") and from eq 5 (denoted by "corr").

Transient Results.
Besides the steady-state performance, it is also relevant to assess the transient performance of the 1D translation approach, especially in the context of temperature swing adsorption applications, such as direct air capture, which are inherently transient.The transient simulations in 1D are done in two different ways: one with a constant Nusselt number and one with a time-dependent Nusselt number.In the first case, the steady-state Nusselt number calculated from the corresponding 2D simulation is used, whereas in the second case, eq 2 is used with Pe 1.0157 D replaced by the steady-state Nusselt number from the 2D simulation.The first case is denoted by an "s" in the figures in this section with the second case by a "d".
Figure 20 shows the dimensionless outlet temperatures for case 1 over time for three different Pećlet numbers.As expected based on the results in Section 3.2, the steady-state temperatures do not exactly match.The transient behavior, however, is very similar, and no large deviations are visible from these figures.At small times, the 1D temperature rises faster than the 2D temperature.This is because in the 1D model, the distance from the last heat source to the system outlet is smaller than in the 2D model due to the translation approach.The 1D temperature profiles with a dynamic Nusselt number ("d") rise faster than those with a steady Nusselt number ("s").This is expected based on eq 2. For very small Pećlet numbers, it seems that the steady Nusselt number simulations agree better with the 2D results, but for larger Pećlet numbers, it is not immediately obvious from these figures which approach is more accurate.Figures S-8 and S-10 in the Supporting Information show the same graphs for cases 2 and 3, respectively.The results are similar to those for case 1, showing that the transient model works independent of the geometry.
More or less, the same observations hold for Figure 21, which shows the dimensionless mean temperatures from the 2D and 1D simulations of case 1 over time for three different Pećlet numbers.The profiles have the same shape, even though the steady-state temperatures differ between the 2D and 1D models.The mean temperature profile is slightly more smooth than the outlet temperature profile at high Pećlet numbers (comparing Figure 20 and Figure 21).This is because the outlet temperature profile reflects the start-up effect of the heat sources (e.g., rapid heat transfer at the start of the simulation because the domain is still cold), which gives temperature peaks that are convected to the outlet.Again, in the 1D dynamic Nusselt number simulations, the steady state is reached faster than the steady Nusselt number simulations, but it is not clear from these figures which approach is best.Similar results are obtained for cases 2 and 3 in Figures S9 and S11 in the Supporting Information.
For temperature swing adsorption, the mean bed temperature is of interest.Therefore, in Figures 22−24, the dimensionless time to reach 50 and 90% of the steady-state mean temperature    is plotted for the geometric cases 1−3 for varying Pećlet number.For the 1D curves, the 1D steady-state temperature is considered, while for the 2D curves, the 2D steady-state temperature is taken, even though both steady-state temperatures differ somewhat.This is done to assess the dynamic performance independent of the steady-state deviations of the 1D model.Again, the "s" for the 1D curves shows that a steady Nusselt number was used, while the "d" shows that a dynamic Nusselt was used.From the figures, it is clear that the time to reach 90% of the steady-state temperature predicted by the 1D model is very close compared to the one calculated by the 2D model for all geometric cases and especially for the steady Nusselt simulations.A slightly larger deviation of the 1D model is apparent in the time to reach 50% of the steady-state time; however, agreement (especially using the steady-state Nusselt simulations) is still acceptable.
From Figures 22−24, it becomes also clear that in most cases, using the steady-state Nusselt number in the 1D model gives more accurate results, especially considering 90% times.Most probably, the cause for this is the distribution of the bed around the heating pipes in the 2D model (see also Figure S1 in the Supporting Information), which causes the bed very close to the tubes to heat faster in 2D, causing the lower t 50% in 2D.However, the material much further away from the tubes takes a lot longer to heat up in 2D, which causes the t 90% s to be similar in 1D and 2D.This effect obviously cannot be taken fully into account in the 1D model, hence the deviation.At small Pećlet numbers, in the dynamic Nusselt 1D simulations, the 50 and 90% steadystate temperatures are reached too fast in all three geometric cases, meaning that the dynamic interplay between the Nusselt number and heat transfer distance d̅ used in the 1D simulations leads to too fast heat transfer predictions.At high Pećlet numbers, both the steady and dynamic Nusselt 1D simulations predict a too short 90% steady-state time, although the steady Nusselt simulations are closer to the 2D results.The trends in the 1D and 2D results also differ here: while for the 2D results, the 90% time increases at high Pećlet numbers, it flattens or even decreases for the 1D results.This is because at high Pećlet numbers, it takes comparatively more dimensionless time before the furthest corners in the 2D domain are fully heated, an effect which cannot be fully captured in the 1D model.For the 50% steady-state time, at large Pećlet numbers, the steady Nusselt number simulations also have better agreement.For most conditions, the error of the 1D dynamic model is below 20%.
Similarly, Figures 25, 26, and S12 in the Supporting Information show the time to reach 50 and 90% of the steadystate outlet temperature for the three geometric cases as a function of the Pećlet number.At low Pećlet numbers, the agreement between the 1D and 2D simulations both for the 50 and 90% times is good, especially for the steady-Nusselt number ("s") simulations.At larger Pećlet numbers, the agreement is less good with the 1D model generally underpredicting the 50 and 90% steady-state times.Only for case 1, the 1D model    overpredicts the 50% steady-state time.This is probably caused by the flow profile in case 1 (where the tubes are fully in line).In this case, in the 2D model, the flow can easily go in between the tube rows, leading to a fast development of the outlet temperature.In the other cases, the flow and temperature profiles have to develop in between the staggered tubes, which takes longer.In the 1D model, this difference cannot be included properly, leading to the different behavior.Regarding the steady Nusselt versus the dynamic Nusselt simulations, the same general observation as for the mean temperature curves holds: generally, the steady-state Nusselt number model is more accurate as the dynamic Nusselt number model overpredicts the Nusselt number.
The results show that for regular tube arrangements, the 2D to 1D translation approach gives acceptable results with a steadystate error of under 10% and a dynamic error of below 20%.The computational costs of the 1D model are a factor 500 less than those for the 2D model, both for the steady-state and dynamic simulations showing the benefit of this approach.Furthermore, in the 1D model, it is easily possible to add chemical reactions or adsorption kinetics, species transport, and nonconstant physical properties without too much additional computational effort.In the 2D model, including such features (especially chemical reactions) leads to a much more computationally difficult model and therefore much longer calculation times.

CONCLUSIONS
In the present work, a method is developed to capture 2D heat transfer effects in regular cylinder banks in a 1D model.The benefit of this approach is a significant decrease in computational effort when modeling such systems while retaining a reasonable accuracy.One area where this approach can be used successfully is when optimizing fixed bed direct air capture systems, where 2D heat transfer is important, but including chemical reactions and mass transfer in a 2D model leads to prohibitively large computational times.
Correlations for calculating the Nusselt number for a cylinder bank embedded in a fixed bed were investigated.For a staggered cylinder arrangement, the analytical single cylinder Nusselt correlation was sufficiently accurate; however, for in-line arrangements, the thermal wake significantly decreases heat transfer around downstream rows.An analytical correction factor for the single cylinder correlation was therefore derived, tested, and proven to be successful.
The 2D to 1D translation approach was tested by modeling three geometries in 2D and translating them to a 1D model.Then, the mean and outlet mixing-cup averaged temperatures as predicted by the 2D and 1D models were compared for both steady-state and dynamic cases.For the steady-state results, the 1D model results (mean bed temperature and outlet temperature) are within 10% of the 2D results for all Pećlet numbers between 1 and 1000 and geometries.
The dynamic performance was tested by comparing the times to reach 50 and 90% of the steady-state temperature in the 2D and 1D models.Especially, the 90% time predicted by the 1D model is in good agreement with the 2D results with an error generally less than 20%.The deviation of the 50% time is still acceptable but somewhat larger, probably because here, the irregular parts of the domain (corners, etc.) play a larger role.Looking at the computational costs, the 1D model is approximately 500 times faster than the 2D model.This, together with the good agreement with the 2D results, shows the merit of the translation approach.■ ACKNOWLEDGMENTS This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Figure 1 .
Figure 1.Schematic view of the translation approach.

Figure 2 .
Figure 2. Case 1 geometry with a temperature profile for Pe D = 50.

Figure 3 .
Figure 3. Case 2 geometry with a temperature profile for Pe D = 50.

Figure 4 .
Figure 4. Case 3 geometry with a temperature profile for Pe D = 50.

Figure 5 .
Figure 5. Schematic view of the cylinders in series in a packed bed in a channel.

Figure 6 .
Figure 6.Flux over cylinder i divided by the average flux for ten cylinders in a row with a pitch of 2 at different Pećlet numbers (color scale) as calculated from 2D results compared with the relation i 1 1 + .

2 . 5 .
1D Model.In the 1D model, the reactor is modeled in the net flow direction, which is the axial direction.The energy

Figure 7 .
Figure 7. Correction of the bed length and diameter when the heating pipes are not included in the bed.Note that in both domains, the area of the white part is equal.

Figure 8 .
Figure 8. Placement of the heat sources in the 1D domain.

Figure 9 .
Figure 9. Simplified schematic depiction of the translation from the 2D model to the 1D model.

Figure 10 .
Figure 10.Schematic view of the spatial grid used in the discretization of the one-dimensional transport equations.

Figure 11 .
Figure 11.Comparison of the Nusselt number calculated from the correlations and calculated from the 2D simulation results.

Figure 12 .
Figure 12.Dimensionless heat flux distribution around five in-line cylinders in a row.

Figure 13 .
Figure 13.Analytical temperature correction factor (lines) compared to the factor calculated from 2D simulations (markers) as a function of dimensionless pitch and the Pećlet number for two cylinders in series.

Figure 15 .
Figure 15.Comparison of steady-state cup-mixed average axial temperature profiles of the 2D and 1D simulations for different Pećlet numbers for case 1.

Figure 16 .
Figure 16.Comparison of steady-state cup-mixed average axial temperature profiles of the 2D and 1D simulations for different Pećlet numbers for case 2.

Figure 17 .
Figure 17.Comparison of steady-state cup-mixed average temperatures of the 2D and 1D simulations as a function of the Pećlet number for case 1, using the Nusselt number calculated from 2D simulations (denoted by "sim") and from eq 5 (denoted by "corr").

Figure 18 .
Figure18.Comparison of steady-state cup-mixed average temperatures of the 2D and 1D simulations as a function of the Pećlet number for case 2, using the Nusselt number calculated from 2D simulations (denoted by "sim") and from eq 5 (denoted by "corr").

Figure 19 .
Figure19.Comparison of steady-state cup-mixed average temperatures of the 2D and 1D simulations as a function of the Pećlet number for case 3, using the Nusselt number calculated from 2D simulations (denoted by "sim") and from eq 5 (denoted by "corr").

Figure 20 .
Figure 20.Comparison of transient outlet temperature profiles for different Pećlet numbers for case 1.

Figure 21 .
Figure 21.Comparison of transient mean temperature profiles for different Pećlet numbers for case 1.

Figure 22 .
Figure 22.Comparison of the time to reach 50 and 90% of the steadystate mean temperature with varying Pećlet numbers for case 1.

Figure 23 .
Figure 23.Comparison of the time to reach 50 and 90% of the steadystate mean temperature with varying Pećlet numbers for case 2.

Figure 24 .
Figure 24.Comparison of the time to reach 50 and 90% of the steadystate mean temperature with varying Pećlet numbers for case 3.Figure 25.Comparison of the time to reach 50 and 90% of the steadystate outlet temperature with varying Pećlet numbers for case 1.

Figure 25 .
Figure 24.Comparison of the time to reach 50 and 90% of the steadystate mean temperature with varying Pećlet numbers for case 3.Figure 25.Comparison of the time to reach 50 and 90% of the steadystate outlet temperature with varying Pećlet numbers for case 1.
15veloped several expressions for the transient Nusselt number around a single cylinder, such as an infinite series expression for the transient Nusselt number valid for all time values and a wide range of Pećlet numbers and a simpler relation for large Pećlet numbers and all times.Equation 2 shows the latter expression, where τ D is the dimensionless time, defined as It is clear that this dynamic expression equals the steady-state expression in eq 1, multiplied by a time-dependent expression that decays to one when τ D → ∞.Sano15found that this expression has an error of less than 3% for Pe D ≥ 400.For small Pećlet numbers, the accuracy of this expression decreases due to the increased contribution of conduction.In these conditions, Sano15advises to use the eulerized small time solution from their work.
== and E is the complete elliptic integral of the second kind.

Table 1 .
List of Symbols Used in This Paper in Alphabetical Order