Thermo-Hydraulic Characteristics of Micro-Scale Surface Roughness Topology of Additively Manufactured Surface Model: Modal Decomposition Perspective

: As a consequence of rapid development of additive manufacturing (3D printing) methods, the academic/industrial demand has been continuously increasing. One ﬁeld of application is the manufacturing of heat exchanging devices using this promising method. In this regard, understanding the underlying mechanisms from a thermo-hydraulic viewpoint becomes important. Therefore, in this study, scale-resolving large eddy simulation (LES) is applied to reveal the ﬂow details in combination with a model of roughness topology occurring in additive manufacturing. To process the transient LES results, proper orthogonal decomposition (POD) is used to extract the coherent ﬂow structures, and the extended POD is used to rank the ﬂow modes based on thermal importance. The main aim of the present work is to go beyond the conventionally applied methodologies used for the evaluation of surface roughness, i.e., averaged numerical study or experimental overall performance evaluation of the ﬂow/thermal response of additively manufactured surfaces in heat exchangers. This is necessary to reveal the underlying ﬂow mechanisms hidden in the conventional studies. In this study, the behavior of the ﬂow over the micro-scale surface roughness model and its effects on heat transfer are studied by assuming cone-shaped roughness elements with regular placement as the dominant surface roughness structures. The major discussions reveal the footprint of ﬂow mechanisms on the heat transfer coefﬁcient spatial modes on the rough surface. Moreover, comparative study on the ﬂow/thermal behavior at different levels of roughness heights shows the key role of the height-to-base-diameter ratio of the roughness elements in thermal performance.


Introduction
The trending additive manufacturing 3D printing method has been recognized as a solution to manufacture light-weight/compact heat exchanging devices with minimized material waste and facilitation of freedom in designing complex elements. Generally, the technology applies the additive shaping principle. As a result, the three-dimensional body is made by successive addition of material [1]. The term 'additive' opposes the conventional 'subtractive' manufacturing methods, where selective removal of material is used to build the specific shape of bodies. Examples are milling, turning, and drilling [1]. In fact, in additive manufacturing, a part is built by layer-by-layer material addition, minimizing the overall waste of raw material [2]. Selective laser sintering, selective laser melting (very similar to selective laser sintering and differing primarily in the materials they use), electron beam melting, and direct metal laser sintering are the additive manufacturing methods dealing with metals [3,4].
One feature of 3D-printed heat exchanging devices that is important to thermal performance is the minimization of contact thermal resistance caused by conventional attachment of elements on surfaces [4]. Aside from the plus points, the challenges of the additive manufacturing must be considered; overhang printing, removal of supporting structures, compatible raw materials, and minimization of the overall costs can be named as a few [2]. In addition to the mentioned positive and challenging features, the processinduced surface roughness of additive manufacturing-enabled materials has become a topic from the perspective of fluid mixing and heat transfer consequences [5].
The growing trend of 3D printing technology in manufacturing heat exchanging and thermal management devices on the one hand, and the naturally rough surface of the 3D-printed parts on the other hand make the study of the surface roughness effect on the thermo-hydraulic performance of additively manufacturing heat exchangers a critical topic. Stimpson et al. [5] evaluated the surface roughness produced by the direct metal laser sintering process and measured its influence on flow losses and convective heat transfer through rectangular minichannels. Their work contains experimental data to study the overall performance of the rough surface in terms of the pressure loss and heat transfer coefficient. However, the detailed study of flow/thermal behavior is missing. The numerical study of rough surface effects on thermo-hydraulic performance has been conducted through both averaged models and detailed flow study methods, such as direct numerical simulation (DNS) and large eddy simulation (LES). Isaev et al. [6] provided a numerical study based on Reynolds averaged Navier Stokes (RANS) turbulence modelling on vortex heat transfer enhancement in the presence of oval-trench surface dimples. The focus of this study was on the mechanism of secondary flow restructuring and heat transfer enhancement due to the formation of vortex structures in the dimples and spiral vortices in behind. Despite the undeniable strength of RANS/URANS in providing desirably fast study of engineering cases in the presence of a turbulent flow regime, there are research items reporting unrealistic results obtained by such models in comparison to the scaleresolving LES model and DNS. Turnow et al. [7] provided a study on the flow structures and heat transfer in a staggered array of dimples on a surface in the presence of turbulent flow. They reported the flow structures determined by LES inside the dimple as chaotic, containing a wide range of scales where coherent structures are hard to detect. The applied POD showed tornado-like flow structures, and the different geometrical configurations of the dimples on the surface were shown to have considerable capacity to control the dimple-induced flow mixing.
Among the methods of analyzing the numerical/experimental data of turbulent flow thermal fields, POD has attracted considerable attention. Lumley [8] proposed this method in the field of fluid dynamics. The primary idea was to represent the complex turbulent flow organization using a set of deterministic functions. In this case, each of the functions contains a share of total turbulence kinetic energy. Eßl et al. [9] employed POD analysis of the gas-jet wiping process in hot-dip galvanizing. Using the decomposition analysis, they extracted the energetic flow modes of the flow. The identified energetic POD modes showed flapping and axial fluctuations of the jet core length. Extended-POD (EPOD) developed by Borée [10] is the method to correlate two quantities using the same temporal basis. In this method, the modes of one generic quantity are found by solving the eigenvalue problem for the other quantity. Lohrasbi et al. [11] provided a study on the correlation of flow and temperature fields for the domain with vortex generators of different angles of attack responsible for causing different signals within the computational domain. In this work, a spatial consideration was expressed to be adopted as the prerequisite for domain selection to ensure uncontaminated correlation study by EPOD.
In the present paper, an analysis of the thermal/flow behavior of a natural surface roughness topology arising in additive manufacturing from the perspective of modal decomposition analysis of flow and thermal fields is targeted. To the authors' knowledge, the detailed study of flow mechanisms and thermal consequences over 3D-printed surfaces is missing in literature, and with the growing trend of 3D printing technology, understanding the flow/thermal aspects of such surfaces becomes crucial. Therefore, in the present study, the scale-resolving turbulence model in the framework of LES along with POD/EPOD methods are used to extract the coherent structures and study the thermal-flow correlation. The dominant energetic flow modes and the thermally important modes are extracted, and the behavior of the natural roughness model from the thermo-hydraulic perspective is discussed. The surface roughness in this work is assumed as cone-shaped peak/valley structures at three levels of heights with regular placement on the surface ignoring the randomness of height/placement of the roughness elements in the real rough surface. The paper is structured as follows. Firstly, the employed method is described and validated. Secondly, model description is provided. Thirdly, the results of the velocity POD/EPOD study are provided for the roughness model, and the modes are described. Finally, POD is applied on the Nusselt field in order to examine the footprint of the flow modes on the modes corresponding to the heat transfer coefficient on the rough surface.

Methodology
ANSYS Fluent 19.1 is used to study turbulent flow and heat transfer. For this, Large Eddy Simulation is adopted whose governing equations are found by filtered Navier-Stokes equations in either Fourier space or configuration space, so that eddies smaller than the filter width or grid spacing are filtered out. Consequently, the obtained governing equations represent the dynamics of large eddies. For this, the filtered variable is defined by u = D G(x, x )u(x )dx , where overbar denotes the filtered field, D represents the fluid domain, and G is the filter function determined based on the volume of the computational cell [12]: ∂ρ ∂t where h and T represent filtered enthalpy and temperature, respectively. The subgridscale turbulence models for modelling subgrid-scale stresses in ANSYS FLUENT use the Boussinesq hypothesis as in the RANS models. The subgrid-scale turbulent stresses are calculated by: where µ t is the subgrid-scale turbulent viscosity and S ij is the rate-of-strain tensor for the resolved scale: Eddy viscosity is modelled by the so-called wall-adapting local eddy-viscosity (WALE) model: where κ represents the von Karman constant (0.41), d w represents the distance to the closest wall, C w is the WALE constant (0.325), and V is the volume of the computational cell. The working fluid is water/glycol (50-50) with the thermo-physical properties provided in Table 1. Validation of the LES calculation is done as shown in Figure 1. For this, a turbulent flow over a cylinder at Re = 3900 is adopted to check the validity of the method in predicting the separation of the shear layers in the presence of an obstacle, compared with the experimental data provided by Parnaudeau et al. [13]. closest wall, is the WALE constant (0.325), and V is the volume of the computational cell.
The working fluid is water/glycol (50-50) with the thermo-physical properties provided in Table 1. Validation of the LES calculation is done as shown in Figure 1. For this, a turbulent flow over a cylinder at Re = 3900 is adopted to check the validity of the method in predicting the separation of the shear layers in the presence of an obstacle, compared with the experimental data provided by Parnaudeau et al. [13].  Figure 1 indicates acceptable accordance of the LES-based simulation results considering the span-wise periodicity with the experimental data [13].
Proper orthogonal decomposition-developed by Lumley [8] as discussed beforehas the capability of extracting the coherent flow structures and ranking the flow modes based on the energy content. In fact, each one of these POD modes contains a share of the total turbulent kinetic energy (TKE) of the flow field. Besides, the modes are ranked based on the energy content. The input of POD is the collected transient field data in the domain of interest by a high-enough number of snapshots. To apply POD on the fluctuating part of velocity, the following decomposition is the starting point of the calculations [14]: where ⃗ ( ⃗ ) represents the velocity field ( = 1 to ), is the number of snapshots, and ⃗ is the spatial coordinates ( = 1 to ), where is the number of spatial positions. The fluctuating part of the velocity field can be represented by a linear set of space-  Figure 1 indicates acceptable accordance of the LES-based simulation results considering the span-wise periodicity with the experimental data [13].
Proper orthogonal decomposition-developed by Lumley [8] as discussed beforehas the capability of extracting the coherent flow structures and ranking the flow modes based on the energy content. In fact, each one of these POD modes contains a share of the total turbulent kinetic energy (TKE) of the flow field. Besides, the modes are ranked based on the energy content. The input of POD is the collected transient field data in the domain of interest by a high-enough number of snapshots. To apply POD on the fluctuating part of velocity, the following decomposition is the starting point of the calculations [14]: where → u i → x k represents the velocity field (i = 1 to N), N is the number of snapshots, and → x k is the spatial coordinates (k = 1 to M), where M is the number of spatial positions. The fluctuating part of the velocity field can be represented by a linear set of space-dependent spatial modes → φ j and the corresponding purely time-dependent temporal coefficients a ij : According to the method of snapshots [15], spatial POD modes and temporal coefficients can be calculated based on the eigenvalues (λ j ) and eigenvectors (e j ) of the N × N correlation matrix ( C) [9]: The temporal coefficients corresponding to the POD modes are calculated as follows: The spatial modes satisfy the orthonormality condition as follows: The spatial modes are calculated as follows: Note that each j th mode → φ j contains a share (λ j ) of total turbulent kinetic energy. Borée [10] proposed extended proper orthogonal decomposition (EPOD) to conduct correlation analysis between two vector fields or scalar fields, e.g., to study correlation between a flow field and temperature field. To achieve this, the eigenvalue problem for one quantity is solved. Following this, the obtained set of eigenvectors is used to find the modes of the other quantity (in this case, correlation is assessed between the temperature and flow fields). The correlated EPOD modes are hence found by: where g i is the field whose EPOD mode is to be calculated. The correlated field can be found using the EPOD mode [14]: Additionally, the uncorrelated field can be found by subtraction from the original field:

Model Description
Based on the roughness topology provided in [5] corresponding to a part manufactured by direct metal laser sintering, the roughness model for the current POD/EPODbased study is considered as cone-shaped peaks/valleys at three levels of heights (at each level, the roughness elements are of the same height) with a fixed base diameter It must be noted that the randomness of geometrical parameters corresponding to the roughness elements and their placement, which is naturally arising in the real surface roughness topology, is not considered in the model. In fact, the model takes into account the dominant roughness elements as cone-shaped peaks/valleys to extract the dominant flow mechanisms governing the heat transfer from the rough surface. The investigated cases in this study include roughness heights of 20 µm, 40 µm, and 60 µm, corresponding to a height-to-base diameter ratio of H/D = 0.175, 0.350, and 0.525, respectively. The span-wise and stream-wise distance between the cone peaks in the cone-shaped roughness model are assumed as 127.8 µm and 134.0 µm, respectively. The same is used for the stream-wise and span-wise distance between valleys. The computational domain and the applied boundary conditions are provided in Figure 2, and the three levels of roughness are illustrated in Figure 4. study include roughness heights of 20 μm, 40 μm, and 60 μm, corresponding to a heightto-base diameter ratio of / = 0.175, 0.350, and 0.525, respectively. The span-wise and stream-wise distance between the cone peaks in the cone-shaped roughness model are assumed as 127.8 μm and 134.0 μm, respectively. The same is used for the stream-wise and span-wise distance between valleys. The computational domain and the applied boundary conditions are provided in Figure 2, and the three levels of roughness are illustrated in Figure 4. As shown in Figure 2, stream-wise and span-wise periodicity is considered for the computational domain, and the mass flow rate (m ) in the stream-wise direction is specified along Z, with a specific inflow bulk temperature (T = 363.15 K). The bottom wall is subject to constant heat flux (10 W m ⁄ ). The top wall is taken to be adiabatic with a noslip condition. The mass flow rate over the periodic domain corresponds to Re = 1.26 × 10 , where the Re is calculated based on the height of the periodic domain [7].
To give an overview of the generated mesh and corresponding resolution, the threedimensional mesh for the periodic domain studying the roughness effect for the mid-level of roughness is shown in Figure 3. Cut-Cell Meshing technology is used in ANSYS-Meshing in order to achieve the desired near-wall mesh refinement with close-to-one aspect ratios. For the illustrated mesh, the total number of cells is 2.3 × 10 and the near-wall mesh resolution is ∆ ∆ 14.2 and ∆ 0.20, where the provided values for ∆ , ∆ , and ∆ are the facet maximum values representing the stream-wise, wallnormal, and span-wise resolutions, respectively. In more detail, ∆ represents the nondimensional wall distance for a wall-bounded flow defined as ∆ ≡ 2 ⁄ , where is the fluid density, is the friction velocity at the nearest wall, y is the distance to the nearest wall, and is the local dynamic viscosity of the fluid. Moreover, ∆ and ∆ are calculated as the square root of the multiplication of and face area magnitude, i.e., the magnitude of the face area vector for the non-internal faces, divided by the cell wall distance, i.e., the distribution of the normal distance of each cell centroid from the wall boundaries [16]. It is noteworthy that the difference in the resolution report of the wallnormal and stream-wise directions is due to the very thin inflation layers.   Figure 2, stream-wise and span-wise periodicity is considered for the computational domain, and the mass flow rate ( . m) in the stream-wise direction is specified along Z, with a specific inflow bulk temperature (T b = 363.15 K). The bottom wall is subject to constant heat flux (10 5 W/m 2 ). The top wall is taken to be adiabatic with a no-slip condition. The mass flow rate over the periodic domain corresponds to Re = 1.26 ×10 4 , where the Re is calculated based on the height of the periodic domain [7].

As shown in
To give an overview of the generated mesh and corresponding resolution, the threedimensional mesh for the periodic domain studying the roughness effect for the mid-level of roughness is shown in Figure 3. Cut-Cell Meshing technology is used in ANSYS-Meshing in order to achieve the desired near-wall mesh refinement with close-to-one aspect ratios. For the illustrated mesh, the total number of cells is 2.3 × 10 +6 and the near-wall mesh resolution is ∆x + ≈ ∆z + ≈ 14.2 and ∆y + ≈ 0.20, where the provided values for ∆x + , ∆y + , and ∆z + are the facet maximum values representing the stream-wise, wall-normal, and span-wise resolutions, respectively. In more detail, ∆y + represents the non-dimensional wall distance for a wall-bounded flow defined as ∆y + ≡ 2ρu τ y/µ, where ρ is the fluid density, u τ is the friction velocity at the nearest wall, y is the distance to the nearest wall, and µ is the local dynamic viscosity of the fluid. Moreover, ∆x + and ∆y + are calculated as the square root of the multiplication of y + and face area magnitude, i.e., the magnitude of the face area vector for the non-internal faces, divided by the cell wall distance, i.e., the distribution of the normal distance of each cell centroid from the wall boundaries [16]. It is noteworthy that the difference in the resolution report of the wall-normal and stream-wise directions is due to the very thin inflation layers.  Menter's SST k-ω [17] is conducted as a precursor RANS simulation to identify the suitable mesh resolution, and once this is done, the time step size is estimated. The bounded second-order implicit method is used for the transient formulation, and the iterative time advancement method is applied in the entire calculation. The LES calculations are initialized by the converged RANS steady solution with a small time step size corresponding to Courant numbers as low as 0.1, and a high number of iterations per time step (as high as 100). Following this, the time step size is ramped up to CFL 10, and the number of iterations per time step is also reduced to 20. During this time, the domain is washed at least once. The required time duration is roughly determined by the channel length and flow characteristic velocity. Finally, the time step size is ramped down to CFL around one up to the statistically averageable state. Once this is reached, the snapshots containing instantaneous thermal/flow fields are exported. The final time step used in this procedure is 5 × 10 s, and the maximum No. of iterations per time step is 100 at the beginning and 20 when exporting snapshots. The total No. of time steps for exporting snapshots is 12,000 with 4 as the data exporting frequency. The sub-grid scale model is the WALE model with = 0.325. The SIMPLE scheme is used for pressure-velocity coupling. The spatial discretization of the gradient, pressure, momentum, and energy is done by least squares cellbased, second order, bounded central differencing, and second order upwind, respectively.  Menter's SST k-ω [17] is conducted as a precursor RANS simulation to identify the suitable mesh resolution, and once this is done, the time step size is estimated. The bounded second-order implicit method is used for the transient formulation, and the iterative time advancement method is applied in the entire calculation. The LES calculations are initialized by the converged RANS steady solution with a small time step size corresponding to Courant numbers as low as 0.1, and a high number of iterations per time step (as high as 100). Following this, the time step size is ramped up to CFL ≈ 10, and the number of iterations per time step is also reduced to 20. During this time, the domain is washed at least once. The required time duration is roughly determined by the channel length and flow characteristic velocity. Finally, the time step size is ramped down to CFL around one up to the statistically averageable state. Once this is reached, the snapshots containing instantaneous thermal/flow fields are exported. The final time step used in this procedure is 5 × 10 −8 s, and the maximum No. of iterations per time step is 100 at the beginning and 20 when exporting snapshots. The total No. of time steps for exporting snapshots is 12,000 with 4 as the data exporting frequency. The sub-grid scale model is the WALE model with C w = 0.325. The SIMPLE scheme is used for pressure-velocity coupling. The spatial discretization of the gradient, pressure, momentum, and energy is done by least squares cell-based, second order, bounded central differencing, and second order upwind, respectively.  Menter's SST k-ω [17] is conducted as a precursor RANS simulation to identify the suitable mesh resolution, and once this is done, the time step size is estimated. The bounded second-order implicit method is used for the transient formulation, and the iterative time advancement method is applied in the entire calculation. The LES calculations are initialized by the converged RANS steady solution with a small time step size corresponding to Courant numbers as low as 0.1, and a high number of iterations per time step (as high as 100). Following this, the time step size is ramped up to CFL 10, and the number of iterations per time step is also reduced to 20. During this time, the domain is washed at least once. The required time duration is roughly determined by the channel length and flow characteristic velocity. Finally, the time step size is ramped down to CFL around one up to the statistically averageable state. Once this is reached, the snapshots containing instantaneous thermal/flow fields are exported. The final time step used in this procedure is 5 × 10 s, and the maximum No. of iterations per time step is 100 at the beginning and 20 when exporting snapshots. The total No. of time steps for exporting snapshots is 12,000 with 4 as the data exporting frequency. The sub-grid scale model is the WALE model with = 0.325. The SIMPLE scheme is used for pressure-velocity coupling. The spatial discretization of the gradient, pressure, momentum, and energy is done by least squares cellbased, second order, bounded central differencing, and second order upwind, respectively.   For the modal decomposition analysis, the snapshots of the flow/thermal fields are reported on the horizontal and wall-normal stream-wise planes as illustrated in Figure 5. The horizontal plane in all cases is located at = 2/3 , where is the peak height, and the wall-normal plane is located at the middle of the periodic domain along the X-axis at X = 0 as shown in the figure. In the present study, the shown 2D planes are used to export the flow/thermal data to conduct the modal decomposition analysis at a reduced computational cost in terms of input files. The limits of this approach must be considered; specifically, the mixing between the flow mechanisms observed in these two planes and the absence of the out-of-plane velocity component affects the analyses. Modal decomposition study of 2D planes extracted from the fluid domain has been used in numerous studies (see [18][19][20]) According to Meyer et al. [20], POD study conducted on planes would agree well with the "true" three-dimensional analysis as long as the plane contains the important phenomena in the flow. Accordingly, flow/thermal snapshots in this paper are exported on a horizontally mounted plane in addition to the wall-normal stream-wise plane so as to consider the span-wise mixing effect in the horizontal plane that is not observed in the stream-wise wall-normal plane.

Modal Decomposition Analysis of the Velocity Field on the Rough Surface Model
The TKE-ranked POD modes of the stream-wise velocity in the wall-normal plane for all three levels of roughness heights are shown in Figure 6. The modes corresponding to the high level of roughness indicate that among the low-order energetic flow modes, the so called Kelvin-Helmholtz instability [21] is observed (see modes No. 2, 3, 4, and 5). For the modal decomposition analysis, the snapshots of the flow/thermal fields are reported on the horizontal and wall-normal stream-wise planes as illustrated in Figure 5. The horizontal plane in all cases is located at Y = 2/3H, where H is the peak height, and the wall-normal plane is located at the middle of the periodic domain along the X-axis at X = 0 as shown in the figure. In the present study, the shown 2D planes are used to export the flow/thermal data to conduct the modal decomposition analysis at a reduced computational cost in terms of input files. The limits of this approach must be considered; specifically, the mixing between the flow mechanisms observed in these two planes and the absence of the out-of-plane velocity component affects the analyses. Modal decomposition study of 2D planes extracted from the fluid domain has been used in numerous studies (see [18][19][20]) According to Meyer et al. [20], POD study conducted on planes would agree well with the "true" three-dimensional analysis as long as the plane contains the important phenomena in the flow. Accordingly, flow/thermal snapshots in this paper are exported on a horizontally mounted plane in addition to the wall-normal stream-wise plane so as to consider the span-wise mixing effect in the horizontal plane that is not observed in the stream-wise wall-normal plane.  For the modal decomposition analysis, the snapshots of the flow/thermal fields are reported on the horizontal and wall-normal stream-wise planes as illustrated in Figure 5. The horizontal plane in all cases is located at = 2/3 , where is the peak height, and the wall-normal plane is located at the middle of the periodic domain along the X-axis at X = 0 as shown in the figure. In the present study, the shown 2D planes are used to export the flow/thermal data to conduct the modal decomposition analysis at a reduced computational cost in terms of input files. The limits of this approach must be considered; specifically, the mixing between the flow mechanisms observed in these two planes and the absence of the out-of-plane velocity component affects the analyses. Modal decomposition study of 2D planes extracted from the fluid domain has been used in numerous studies (see [18][19][20]) According to Meyer et al. [20], POD study conducted on planes would agree well with the "true" three-dimensional analysis as long as the plane contains the important phenomena in the flow. Accordingly, flow/thermal snapshots in this paper are exported on a horizontally mounted plane in addition to the wall-normal stream-wise plane so as to consider the span-wise mixing effect in the horizontal plane that is not observed in the stream-wise wall-normal plane.

Modal Decomposition Analysis of the Velocity Field on the Rough Surface Model
The TKE-ranked POD modes of the stream-wise velocity in the wall-normal plane for all three levels of roughness heights are shown in Figure 6. The modes corresponding to the high level of roughness indicate that among the low-order energetic flow modes, the so called Kelvin-Helmholtz instability [21] is observed (see modes No. 2, 3, 4, and 5).

Modal Decomposition Analysis of the Velocity Field on the Rough Surface Model
The TKE-ranked POD modes of the stream-wise velocity in the wall-normal plane for all three levels of roughness heights are shown in Figure 6. The modes corresponding to the high level of roughness indicate that among the low-order energetic flow modes, the so called Kelvin-Helmholtz instability [21] is observed (see modes No. 2, 3, 4, and 5). Due to the appearance among the low-order/high-energy POD modes, this instability is recognized as the dominant flow mechanism in stream-wise wall-normal planes, which is promoted and intensified by the peaks of the roughness profile as the wall-attached obstacles. In addition, the other point is the existence of strong active flow fluctuations in the core flow as well as in the near-wall flow as a result of the strengthened fluctuations by the high-rise peaks in the case of a high level of roughness. In spite of the positive aspects, the negative point for thermal transport is the minimal engagement of valleys in terms of the absence of strong velocity fluctuating structures due to the flow trap effect in the case of high roughness. This becomes clear when comparing the mode shapes with the corresponding ones at the low level of roughness height. This effect comes from the dead fluctuation zone in the valleys. To verify if the valleys are involved in thermally important modes, the EPOD modes of this case are discussed in the next sections. The POD modes of velocity corresponding to the medium level of roughness provided in Figure 6 show the same flow mechanism as the dominant/energetic as for the high roughness level. However, even among the energetic modes, a higher engagement of valleys in terms of strong flow fluctuations is observed. Active fluctuations in this area are important to avoid local dead areas/hot spots. The other point is that in the medium roughness level, a balance is observed between the strong Kelvin-Helmholtz structures and the engagement of valleys. This is in contrast to the case with the low level of roughness as shown later. Finally, when comparing the energy distribution and mode shapes of velocity eigenmodes at the low level of roughness height with the high roughness case, the first inferred point is the similar energy distribution and flow mechanism in the stream-wise wall-normal plane. In fact, the same flow mechanism is introduced by the wall-attached roughness structures. However, the main difference is where the active strong fluctuations take place. In this regard, valleys are more engaged so that the dead fluctuation zone is avoided. The other point is the shift of the coherent flow modes of the low level of roughness to the near-wall area, while the corresponding flow modes at the high level of roughness show strong fluctuating structures in the core area as well as near the wall. This is important to achieve core flow-near wall flow mixing. Due to the appearance among the low-order/high-energy POD modes, this instability is recognized as the dominant flow mechanism in stream-wise wall-normal planes, which is promoted and intensified by the peaks of the roughness profile as the wall-attached obstacles. In addition, the other point is the existence of strong active flow fluctuations in the core flow as well as in the near-wall flow as a result of the strengthened fluctuations by the high-rise peaks in the case of a high level of roughness. In spite of the positive aspects, the negative point for thermal transport is the minimal engagement of valleys in terms of the absence of strong velocity fluctuating structures due to the flow trap effect in the case of high roughness. This becomes clear when comparing the mode shapes with the corresponding ones at the low level of roughness height. This effect comes from the dead fluctuation zone in the valleys. To verify if the valleys are involved in thermally important modes, the EPOD modes of this case are discussed in the next sections. The POD modes of velocity corresponding to the medium level of roughness provided in Figure 6 show the same flow mechanism as the dominant/energetic as for the high roughness level. However, even among the energetic modes, a higher engagement of valleys in terms of strong flow fluctuations is observed. Active fluctuations in this area are important to avoid local dead areas/hot spots. The other point is that in the medium roughness level, a balance is observed between the strong Kelvin-Helmholtz structures and the engagement of valleys. This is in contrast to the case with the low level of roughness as shown later. Finally, when comparing the energy distribution and mode shapes of velocity eigenmodes at the low level of roughness height with the high roughness case, the first inferred point is the similar energy distribution and flow mechanism in the stream-wise wall-normal plane. In fact, the same flow mechanism is introduced by the wall-attached roughness structures. However, the main difference is where the active strong fluctuations take place. In this regard, valleys are more engaged so that the dead fluctuation zone is avoided. The POD modes illustrated in Figure 7 are the TKE-ranked POD modes for the stream-wise velocity in the horizontal plane. As shown in this figure for the high level of roughness height, according to the first few energetic modes, von Kármán instability is observed among the dominant flow mechanisms in the horizontal plane induced by the peaks of surface roughness elements. The von Kármán instability is clearer for a high level of roughness in modes No. 1, 2, 5, 6, 7, and 8. This mechanism is responsible for the intensification of span-wise flow mixing. One point is that this flow mechanism is exclusively produced by the peaks in the roughness profile, and the valleys do not contribute to it. Moreover, due to the desired flow mixing by this mechanism, it is beneficial for thermal transport. This is the other reason behind the beneficial effects of peaks in heat transfer. Finally, it is worth mentioning that this flow mixing is active for a wider wall-normal distance in the case of high-rise peaks compared to the low level of roughness. The mode shapes in this plane at a medium and low level of roughness indicate that the von Kármán mechanism is still observed at a low level of surface element roughness; however, at the low level of roughness height, it is observed in more streak-like shapes in the 2D plane. The fact that the dominance of this flow mechanism is lost is also observed in terms of losing the coherence as shown in the energy distribution diagram, where the trend of the energy distribution shows a strict decreasing behavior. The other point is that aside from the weaker mixing, in this case, the lower height of the peaks responsible for this flow mechanism reduces the overall wall-normal distance where this mixing mechanism exists, which is acting against the roughness-induced heat transfer enhancement.
The reason behind having two different dominant flow mechanisms in the two studied planes, i.e., stream-wise wall normal and horizontal planes, lies in the fact that in the wall-normal planes, flow is taking place over 'wall-attached obstacles'. In accordance with literature [21], for the case of flow over an obstacle with no gap between the obstacle and wall, i.e., fully wall-attached obstacle, the dominant flow mechanism is Kelvin-Helmholtz promoted by the obstacle. This is responsible for small-scale flow structures in the nearwall region promoting heat transfer from the surface. On the other hand, the low-order energetic flow modes in the horizontal plane indicate a different flow mechanism, i.e., von Kármán instability promoted by the peaks of the roughness elements. This is very similar to the simplified case of having pin-shaped obstacles attached to the walls, where von Kármán instability is generated by the pin structures. The POD modes illustrated in Figure 7 are the TKE-ranked POD modes for the stream-wise velocity in the horizontal plane. As shown in this figure for the high level of roughness height, according to the first few energetic modes, von Kármán instability is observed among the dominant flow mechanisms in the horizontal plane induced by the peaks of surface roughness elements. The von Kármán instability is clearer for a high level of roughness in modes No. 1, 2, 5, 6, 7, and 8. This mechanism is responsible for the intensification of span-wise flow mixing. One point is that this flow mechanism is exclusively produced by the peaks in the roughness profile, and the valleys do not contribute to it. Moreover, due to the desired flow mixing by this mechanism, it is beneficial for thermal transport. This is the other reason behind the beneficial effects of peaks in heat transfer. Finally, it is worth mentioning that this flow mixing is active for a wider wall-normal distance in the case of high-rise peaks compared to the low level of roughness. The mode shapes in this plane at a medium and low level of roughness indicate that the von Kármán mechanism is still observed at a low level of surface element roughness; however, at the low level of roughness height, it is observed in more streak-like shapes in the 2D plane. The fact that the dominance of this flow mechanism is lost is also observed in terms of losing the coherence as shown in the energy distribution diagram, where the trend of the energy distribution shows a strict decreasing behavior. The other point is that aside from the weaker mixing, in this case, the lower height of the peaks responsible for this flow mechanism reduces the overall wall-normal distance where this mixing mechanism exists, which is acting against the roughness-induced heat transfer enhancement.
The reason behind having two different dominant flow mechanisms in the two studied planes, i.e., stream-wise wall normal and horizontal planes, lies in the fact that in the wallnormal planes, flow is taking place over 'wall-attached obstacles'. In accordance with literature [21], for the case of flow over an obstacle with no gap between the obstacle and wall, i.e., fully wall-attached obstacle, the dominant flow mechanism is Kelvin-Helmholtz promoted by the obstacle. This is responsible for small-scale flow structures in the near-wall region promoting heat transfer from the surface. On the other hand, the low-order energetic flow modes in the horizontal plane indicate a different flow mechanism, i.e., von Kármán instability promoted by the peaks of the roughness elements. This is very similar to the simplified case of having pin-shaped obstacles attached to the walls, where von Kármán instability is generated by the pin structures. The EPOD modes representing thermally importance-based ranked flow modes are provided in Figure 8. According to the mode shapes corresponding to the high level of roughness, the EPOD modes verify the fact that the active/strong fluctuating structures are mainly introduced by peaks, and the valleys in the case of a high roughness level play a minimal role in thermal transport. The other point inferred from the EPOD modes is the similarity of the mode topologies to the energetic POD modes. This shows the importance of the energetic flow structures in heat transfer in the analyzed 2D stream-wise wall normal planes. The EPOD modes in Figure 8 corresponding to the medium level of roughness show the thermal importance of the large-scale Kelvin-Helmholtz-induced structures followed by small-scale structures distributed in the entire domain including valleys. This verifies the thermal importance of valley-involving flow modes, which is acting against the flow trap effect in these areas. With the reduction of the large-scale flow structures in the case of the low level of roughness, there is a change in the rank of thermally important flow modes compared to the energy-based ranked modes. In fact, the large-scale Kelvin-Helmholtz-induced structures still appear among the first few modes, but above them, the small-scale structures active all over the domain gain thermal importance. In more detail, comparison between Figures 6 and 8 indicates that for the high level of roughness, the same mode topologies appear as the first few energetic POD modes and thermallyimportant EPOD modes. However, when the same comparison is made for the low level of roughness, it is observed that the first few large-scale POD modes that appear as the most energetic among the TKE-ranked modes come after some small-scale flow structures among the thermally important EPOD modes. This is attributed to the fact that the energetic POD modes corresponding to the Kelvin-Helmholtz mechanism at a low level of roughness lose their thermal importance compared to the high level of roughness. In addition to the mentioned points, the engagement of valleys in the thermally important modes can also be observed in this figure.  The EPOD modes representing thermally importance-based ranked flow modes are provided in Figure 8. According to the mode shapes corresponding to the high level of roughness, the EPOD modes verify the fact that the active/strong fluctuating structures are mainly introduced by peaks, and the valleys in the case of a high roughness level play a minimal role in thermal transport. The other point inferred from the EPOD modes is the similarity of the mode topologies to the energetic POD modes. This shows the importance of the energetic flow structures in heat transfer in the analyzed 2D stream-wise wall normal planes. The EPOD modes in Figure 8 corresponding to the medium level of roughness show the thermal importance of the large-scale Kelvin-Helmholtz-induced structures followed by small-scale structures distributed in the entire domain including valleys. This verifies the thermal importance of valley-involving flow modes, which is acting against the flow trap effect in these areas. With the reduction of the large-scale flow structures in the case of the low level of roughness, there is a change in the rank of thermally important flow modes compared to the energy-based ranked modes. In fact, the large-scale Kelvin-Helmholtzinduced structures still appear among the first few modes, but above them, the small-scale structures active all over the domain gain thermal importance. In more detail, comparison between Figures 6 and 8 indicates that for the high level of roughness, the same mode topologies appear as the first few energetic POD modes and thermally-important EPOD modes. However, when the same comparison is made for the low level of roughness, it is observed that the first few large-scale POD modes that appear as the most energetic among the TKE-ranked modes come after some small-scale flow structures among the thermally important EPOD modes. This is attributed to the fact that the energetic POD modes corresponding to the Kelvin-Helmholtz mechanism at a low level of roughness lose their thermal importance compared to the high level of roughness. In addition to the mentioned points, the engagement of valleys in the thermally important modes can also be observed in this figure. The EPOD modes representing thermally importance-based ranked flow modes are provided in Figure 8. According to the mode shapes corresponding to the high level of roughness, the EPOD modes verify the fact that the active/strong fluctuating structures are mainly introduced by peaks, and the valleys in the case of a high roughness level play a minimal role in thermal transport. The other point inferred from the EPOD modes is the similarity of the mode topologies to the energetic POD modes. This shows the importance of the energetic flow structures in heat transfer in the analyzed 2D stream-wise wall normal planes. The EPOD modes in Figure 8 corresponding to the medium level of roughness show the thermal importance of the large-scale Kelvin-Helmholtz-induced structures followed by small-scale structures distributed in the entire domain including valleys. This verifies the thermal importance of valley-involving flow modes, which is acting against the flow trap effect in these areas. With the reduction of the large-scale flow structures in the case of the low level of roughness, there is a change in the rank of thermally important flow modes compared to the energy-based ranked modes. In fact, the large-scale Kelvin-Helmholtz-induced structures still appear among the first few modes, but above them, the small-scale structures active all over the domain gain thermal importance. In more detail, comparison between Figures 6 and 8 indicates that for the high level of roughness, the same mode topologies appear as the first few energetic POD modes and thermallyimportant EPOD modes. However, when the same comparison is made for the low level of roughness, it is observed that the first few large-scale POD modes that appear as the most energetic among the TKE-ranked modes come after some small-scale flow structures among the thermally important EPOD modes. This is attributed to the fact that the energetic POD modes corresponding to the Kelvin-Helmholtz mechanism at a low level of roughness lose their thermal importance compared to the high level of roughness. In addition to the mentioned points, the engagement of valleys in the thermally important modes can also be observed in this figure.

Modal Decomposition of the Transient Nusselt Field on the Rough Surface
Prior to applying modal decomposition study on the Nu field, the temporal mean of the local Nusselt on the heated wall is shown in Figure 9. The local Nu is calculated by Nu = (hD ) k ⁄ , where h, D , and k represent the convective heat transfer coefficient, hy-

Modal Decomposition of the Transient Nusselt Field on the Rough Surface
Prior to applying modal decomposition study on the Nu field, the temporal mean of the local Nusselt on the heated wall is shown in Figure 9. The local Nu is calculated by Nu = (hD h )/k f , where h, D h , and k f represent the convective heat transfer coefficient, hydraulic diameter, and fluid thermal conductivity, respectively. Besides, h is calculated as h = q /(T w − T b ), where q , T w , and T b represent the heat flux, wall temperature, and fluid bulk temperature, respectively. As illustrated in this figure, one plus point of the high-rise roughness peaks is observed as the impingement of the coolant on it. This effect is stronger at a higher level of heights. However, the negative aspect is also observed on the rear side of the peaks in terms of considerably low-Nu areas.

Modal Decomposition of the Transient Nusselt Field on the Rough Surface
Prior to applying modal decomposition study on the Nu field, the temporal mean of the local Nusselt on the heated wall is shown in Figure 9. The local Nu is calculated by Nu = (hD ) k ⁄ , where h, D , and k represent the convective heat transfer coefficient, hydraulic diameter, and fluid thermal conductivity, respectively. Besides, h is calculated as h = q (T − T ) ⁄ , where q , T , and T represent the heat flux, wall temperature, and fluid bulk temperature, respectively. As illustrated in this figure, one plus point of the high-rise roughness peaks is observed as the impingement of the coolant on it. This effect is stronger at a higher level of heights. However, the negative aspect is also observed on the rear side of the peaks in terms of considerably low-Nu areas. POD is applied on the Nu field over the rough surface to evaluate the footprint of the flow and thermal fields on the heat transfer characteristics of the roughness model. Figure  10 shows the mode variance distributions and spatial POD modes of the Nu field on the surface at the three levels of roughness heights. POD is applied on the Nu field over the rough surface to evaluate the footprint of the flow and thermal fields on the heat transfer characteristics of the roughness model. Figure 10 shows the mode variance distributions and spatial POD modes of the Nu field on the surface at the three levels of roughness heights.
Comparison between the spatial mode shapes of Nu on the rough surface between the low and high levels of roughness shows a streakier pattern of Nu eigenmodes when tending from the high level of roughness toward the more smooth surface. In fact, the dominant-high-variance Nu eigenmodes are mainly formed by the streak structures of low/high Nu elongated in the stream-wise direction. Note that even at the roughness height of 20 µm, the span-wise von Kármán mixing still exists, and the structure is not fully streaky. In fact, in the case of a fully smooth surface, the role of the abovementioned small-scale flow structures in the near-wall region is absent. In such cases, the large-scale sweeps/ejections in the system would govern the thermal response of the surface [21], and the small-scale near-wall vortices are almost absent due to the absence of rough elements responsible for such flow structures. Therefore, in such cases, the Nu eigenmodes of streak shapes are expected. However, in all of our cases, the mentioned flow instabilities cause near-wall vortices, but as mentioned, with a tendency toward a smooth surface, the more streak-like spatial mode shapes of Nu are observed. From the comparison between Nu mode shapes of the high and low levels of roughness (here H = 60 µm and 20 µm, respectively), it is inferred that in contrast with the streaky shapes of the POD modes of Nu at the low level of roughness, the pattern is strongly distorted with the increase of the roughness height. This is attributed to the intensification of span-wise mixing as a result of intensifying von Kármán induced by the high-rise peaks on the surface. Comparison between the spatial mode shapes of Nu on the rough surface between the low and high levels of roughness shows a streakier pattern of Nu eigenmodes when tending from the high level of roughness toward the more smooth surface. In fact, the dominant-high-variance Nu eigenmodes are mainly formed by the streak structures of low/high Nu elongated in the stream-wise direction. Note that even at the roughness height of 20 μm, the span-wise von Kármán mixing still exists, and the structure is not fully streaky. In fact, in the case of a fully smooth surface, the role of the abovementioned small-scale flow structures in the near-wall region is absent. In such cases, the large-scale sweeps/ejections in the system would govern the thermal response of the surface [21], and the small-scale near-wall vortices are almost absent due to the absence of rough elements responsible for such flow structures. Therefore, in such cases, the Nu eigenmodes of streak shapes are expected. However, in all of our cases, the mentioned flow instabilities cause near-wall vortices, but as mentioned, with a tendency toward a smooth surface, the more streak-like spatial mode shapes of Nu are observed. From the comparison between Nu mode shapes of the high and low levels of roughness (here = 60 μm and 20 μm, respectively), it is inferred that in contrast with the streaky shapes of the POD modes of Nu at the low level of roughness, the pattern is strongly distorted with the increase of the roughness height. This is attributed to the intensification of span-wise mixing as a result of intensifying von Kármán induced by the high-rise peaks on the surface.
Aside from the mode shapes, the mode variance distributions for the Nu numbers verify the intensification of the thermal role of span-wise mixing with the increase of the roughness height. In fact, with the increase of the roughness height to = 60 μm, the pattern of the Nu mode variance distribution (see Figure 10) becomes very similar to the mode variance distribution pattern of the corresponding height for the flow modes in the horizontal plane, which is dominated by the von Kármán mechanism (see Figure 7). This is attributed to the more dominant role of span-wise mixing caused by von Kármán instability with the increase of height level of high-rise peaks.
From the modal analysis viewpoint, it was indicated how the peaks intensify both observed flow mechanisms, how the valleys do not participate in the von Kármán mechanism, and how valleys can cause a dead fluctuation zone due to the flow trap effect. Accordingly, it seems necessary to provide proposals on how to manipulate additively manufactured surface roughness topology for the sake of thermal transport considera- Aside from the mode shapes, the mode variance distributions for the Nu numbers verify the intensification of the thermal role of span-wise mixing with the increase of the roughness height. In fact, with the increase of the roughness height to H = 60 µm, the pattern of the Nu mode variance distribution (see Figure 10) becomes very similar to the mode variance distribution pattern of the corresponding height for the flow modes in the horizontal plane, which is dominated by the von Kármán mechanism (see Figure 7). This is attributed to the more dominant role of span-wise mixing caused by von Kármán instability with the increase of height level of high-rise peaks.
From the modal analysis viewpoint, it was indicated how the peaks intensify both observed flow mechanisms, how the valleys do not participate in the von Kármán mechanism, and how valleys can cause a dead fluctuation zone due to the flow trap effect. Accordingly, it seems necessary to provide proposals on how to manipulate additively manufactured surface roughness topology for the sake of thermal transport considerations. To be pragmatic, the recommendation for the surface roughness control either during printing or by post-printing surface treatment is to consider the average roughness heightto-base-diameter ratio for a fair sample area on the rough surface as the surface design parameter. This is, because the mid-level of this parameter can balance the intensified Kelvin-Helmholtz mechanism, the von Kármán-induced mixing by peaks, and the reduction of dead fluctuation zones in valleys. The actual thermal response in terms of the maximum surface temperature and surface Nusselt number are provided in Figure 11. It is observed in this figure that the high level of the height-to-base-diameter causes local hot spots, and the close-to-perfectly-smooth surface would lack the desired flow mixing caused by the natural roughness of the surface. As indicated in this figure, the increase of the roughness height leads to an increase of the averaged Nu on the rough surface. However, due to the illustrated dead fluctuating zones in the valleys in the case of the high level of roughness, this increase in Nu is not seen in terms of a maximum temperature reduction. In other words, the dead fluctuating zones cause thermal performance saturation. hot spots, and the close-to-perfectly-smooth surface would lack the desired flow mixing caused by the natural roughness of the surface. As indicated in this figure, the increase of the roughness height leads to an increase of the averaged Nu on the rough surface. However, due to the illustrated dead fluctuating zones in the valleys in the case of the high level of roughness, this increase in Nu is not seen in terms of a maximum temperature reduction. In other words, the dead fluctuating zones cause thermal performance saturation.
(a) (b) Figure 11. Actual thermal performance: (a) maximum facet temperature on the rough surface, and (b) averaged heat transfer coefficient on the rough surface.

Conclusions
The conducted study investigated thermo-hydraulic behavior of a three-dimensional model representing a natural roughness profile observed in additive manufacturing. The surface roughness model is assumed as cone-shaped elements at three levels of height and regular placement of roughness elements as the major roughness structures playing a role in the flow/thermal performance of the system. LES was used to provide detailed flow/thermal fields in the transient form. The transient data were imported into modal decomposition analyses, i.e., POD/EPOD, to extract the energetic and thermally important flow modes. The obtained results indicated that in general, the flow mechanisms responsible for thermal transport are von Kármán in horizontal planes and Kelvin-Helmholtz in stream-wise wall-normal planes. A comparison was made between the POD/EPOD results of the three levels of roughness. Minimization of the dead areas corresponding to hot spots in the micro-scale roughness profile is of crucial importance. This is because regardless of the macro-scale flow mixing promoters, the high number of local dead fluctuation zones containing trapped flow acting as an insulator due to the lower thermal conductivity of common coolants in liquid/air cooling-based devices would reduce the overall performance. This is where roughness can act against the heat transfer purpose. Based on this, the surface roughness manipulation must act in a way that this effect is

Conclusions
The conducted study investigated thermo-hydraulic behavior of a three-dimensional model representing a natural roughness profile observed in additive manufacturing. The surface roughness model is assumed as cone-shaped elements at three levels of height and regular placement of roughness elements as the major roughness structures playing a role in the flow/thermal performance of the system. LES was used to provide detailed flow/thermal fields in the transient form. The transient data were imported into modal decomposition analyses, i.e., POD/EPOD, to extract the energetic and thermally important flow modes. The obtained results indicated that in general, the flow mechanisms responsible for thermal transport are von Kármán in horizontal planes and Kelvin-Helmholtz in stream-wise wall-normal planes. A comparison was made between the POD/EPOD results of the three levels of roughness. Minimization of the dead areas corresponding to hot spots in the micro-scale roughness profile is of crucial importance. This is because regardless of the macro-scale flow mixing promoters, the high number of local dead fluctuation zones containing trapped flow acting as an insulator due to the lower thermal conductivity of common coolants in liquid/air cooling-based devices would reduce the overall performance. This is where roughness can act against the heat transfer purpose. Based on this, the surface roughness manipulation must act in a way that this effect is minimized. The recommendation is to consider the average height-to-base-diameter ratio as the thermal consideration along with the arithmetic average roughness for the design of heat exchangers to be realized by additive manufacturing. While the high ratio roughness benefits from the plus points induced by the high-rise peaks, the thermally inactive dead valleys have a negative effect. At the mid values of this ratio, the high number of local dead areas is reduced, and it still benefits from the shown mixing mechanisms by peaks. It is noteworthy that in the case of a fully smooth surface, the role of von Kármán-induced mixing in horizontal planes and Kelvin-Helmholtz mixing in stream-wise wall normal planes fade. In this case, the transport would be on the shoulder of sweep/ejection mechanisms taking place on smooth surfaces. Finally, it must be noted that the surface roughness model in this work is based on cone-shaped peak/valley structures, so that compared to the real rough surface, the random geometrical parameters of each roughness element and the randomness of placement of roughness elements are absent. Future work is to be conducted on a real 3D-printed surface roughness sample, so that the CFD domain is reconstructed using the point cloud extracted by computed tomography scanning of