Two‐scale conjugate heat transfer solution for micro‐structured surface

The primary challenges for simulating a turbulent flow over a micro‐structured surface arise from the two hugely disparate spatial length scales. For fluid–solid coupled conjugate heat transfer (CHT), there is also a time‐scale disparity. The present work addresses the scale disparities based on a two‐scale framework. For the spatial scale disparity, a dual meshing is employed to couple a global coarse‐mesh domain with local fine‐mesh blocks around micro‐structures through source terms generated from the local fine‐mesh and propagated to the global coarse‐mesh domain. The convergence and robustness of the source terms driven coarse‐mesh solution is enhanced by a balanced eddy‐viscosity damping. In this work, the two‐scale method previously developed only for a fluid‐domain is extended to a solid domain so that thermal conduction around micro‐elements can now be resolved accurately and efficiently. The fluid–solid timescale disparity is dealt with by a frequency domain approach. The time‐averaged (zeroth harmonic) is effectively obtained in the same way as steady CHT. And remarkably, wall temperature unsteadiness can be simply obtained from the fluid temperature harmonics through a wall fluid–solid temperature harmonic transfer‐function at minimal computational cost. The developed CHT capability is validated for an experimental internal cooling channel with multiple surface rib‐elements. For a test configuration with 100 micro‐structures, the fluid domain‐only, the solid domain‐only and the fluid–solid coupled CHT solutions are analyzed respectively to examine and demonstrate the validity of the present framework and implementation methods. Some of the results also serve to illustrate the primary underlying working of the methodology.


INTRODUCTION
Fluid flow and convective heat transfer for micro-structured surfaces are important to many engineering applications. Impact of general stochastic roughness on friction drag and heat transfer performance is a subject of long-standing interest. There have been many efforts in correlating increased friction and heat transfer to the equivalent roughness parameters. 1 The current consensus is that flow velocities in outer-flow region can be reasonably well correlated to a general similarity profile for smooth and roughened surfaces, as first postulated by Townsend. 2 The near-wall flow as well as the main scaling parameters are however dependent on more detailed roughness geometry. 3 Comprehensive understanding of how roughness affects the boundary layer aerothermal performance is still lacking, particularly for relatively large roughness scales. Direct numerical simulations 4 help the understanding of roughness characteristics, and indicate that resolving detailed roughness micro-structure geometries may also be needed in many cases for better characterization of roughness behavior and performance impact than the conventional largely empirical treatment. 5 Enhanced prediction of surface micro-structures is also motivated by the potential performance gains with manufacturable "regular roughness." It has been indicated that detailed shape and pattern of surface micro-structures can considerably affect the aerothermal performance. [6][7][8] The major challenges in the method development for resolving, rather than modeling, surface micro-structures arise from huge scale disparities, in space given typical sizes of micro-structures, as well as in time when fluid-solid thermal coupling is included.
In terms of the spatial disparity, a starting point for micro-structured surfaces is to recognize a scale-dependent solvability. On the one hand, high temporal and spatial gradients within each micro-structure require high resolution locally. On the other hand, the geometrical similarity among the micro-elements leads to a similarity in steady (time-averaged) flow patterns among the micro-elements. The consideration leads to the intent to develop a two-scale approach which couples a local fine-mesh solution for a small set of micro-elements with a global coarse-mesh solution, through source terms generated by upscaling the space-time averaged local fine-mesh solution. 9 Verifications and validations of the two-scale method have been carried out for effusion cooling flow 9 and for a dimple surface. 7 The validity and effectiveness of the two-scale methodology for scale-resolving turbulence solutions has been underpinned for canonical channel flow and tripped turbulence boundary layers respectively compared to well-established full LES/DNS at different Reynolds numbers. 10,11 In the context of convective heat transfer, the applicability of the conventional approach with wall temperatureinvariant heat transfer coefficient (HTC) may be questioned in general. 12 Considerable influence of wall temperature on near-wall flow is indicated, 13,14 as well as on the transition to turbulent flow. 15,16 These observations support the fluid-solid coupled conjugate heat transfer (CHT) method in general. For steady CHT problems, the temporal scale disparity should not have a significant bearing as the fluid and solid domains can adopt very different time steps. For unsteady CHT problems, the physical time-consistency requirement for both fluid and solid parts would lead to a direct conflict. A very small fluid time-step required for sufficient temporal resolution will have to accommodate the long time-scale required for disturbances to propagate through a solid domain. If a direct time-domain CHT solution is performed, the time scale ratio between the two domains can be up to 10 4 for URANS based CHT. 17 Effectively, the number of the time steps for fluid-only solution may have to be multiplied by 10 4 if a direct time domain CHT is performed simply to satisfy a time consistency as required by the solid part. For a scale-resolving turbulent solution based CHT, the challenge will be further exacerbated due to a much smaller time step for resolving turbulence, albeit even just for resolving large-scale eddies.
There are a few attempts in addressing the challenges of unsteady CHT. A widely adopted approach for circumventing the time scale disparity is simply to neglect all unsteadiness in the CHT as commonly adopted in the so-called "loosely coupled" CHT methods. [18][19][20] It would be desirable to have some unsteady CHT capabilities for multiscale unsteady disturbances in a relatively low frequency range (e.g. for low frequency thermo-acoustic dynamics in a gas turbine combustor). Although in many practical situations, wall temperature fluctuations tend to be small, it has also been indicated that magnitudes of wall temperature fluctuations can be considerably affected by solid material properties, particularly thermal conductivity. 21,22 Therefore, there is a need for a consistent CHT framework methodology to deal with unsteady wall temperature fluctuations, which, to the best knowledge of the present author, has been lacking for unsteady CHT in general.
In the present work, a full unsteady variable is decomposed into the time-mean and unsteady fluctuation parts, and the two parts are treated in the two domains differently according to their distinctive characteristics and scales respectively. The time-mean part can be easily treated in the same way as for a steady CHT. The unsteady part in a harmonic form is also effectively solved as a steady-like problem, He and Oldfield. 17 A moving-average based LES-CHT has been developed and applied to a smooth surface of a turbine blade. 22 Given the background and related recent progress as briefly reviewed above, the present work is aimed to develop an efficient and accurate LES based conjugate heat transfer method for micro-structured surfaces, which appears non-existent, to the present author's knowledge. In the following sections, the main framework methodology and implementation methods will be first introduced for both fluid and solid domains, and for the CHT interface. The case studies will then be described for validation and demonstration. Finally, some concluding remarks will be made.

Overall approach and rationale
A unified (monolithic) field solver is preferred for both fluid and solid domains for simplicity, even though the two domains face different challenges arising from the two distinctive scales in space and in time. We thus need to consider, firstly how to deal with the interior parts of the two domains in a unified manner, and secondly how to treat the fluid-solid interface to effectively capture the interactions between the two domains as well as between the different scales. For a fluid domain, the two spatial scales of relevance are the small ("micro") scales with a typical characteristic length of a micro-structure element and the large ("macro") scales of main flow path. A computational domain containing large number micro-structures can be discretised in a coarse base mesh, and it is then further divided into a large number of local fine-mesh blocks around the micro-structures. A common basic conflicting scenario should now be recognized. A local fine-mesh solution should provide a high resolution, but it however can only be useful if the local solution is properly conditioned by its surrounding. On the other hand, the global coarse-mesh domain may be well conditioned, but its solution can only be useful if it is adequately resolved. Here we resort to a two-scale method to address the conflict between the low resolution of the global coarse-mesh domain and the poor conditioning of the local fine-mesh domain. The idea is to solve only a small number of local fine-mesh domains which can be used to set a target solution for the coarse mesh. Source terms can thus be generated and propagated to the global coarse-mesh domain to drive the coarse-mesh solution. In the process for solving the coupled system, the global coarse-mesh and the local fine-mesh are made to correct each other, so that the converged solution for the coupled system will lead to a better conditioned local fine-mesh solution and an effectively better resolved global coarse-mesh solution.
In terms of time-scale disparity, the fluid-solid interface is where a close attention should be paid. The interface is of prime relevance to unsteady CHT in twofold. Firstly, within whole solid domain, the unsteady temperature fluctuations tend to be the highest at a solid wall boundary surface interfacing with the fluid domain. Note that the principal (if not the only) pathway for wall temperature fluctuations to influence the time-averaged heat transfer is to alter the time-averaged near-wall flow and thermal field through local nonlinear turbulent interactions. Secondly the interface is also where the challenge of time scale disparity manifests itself most acutely.
It should be recognized that treating unsteady solid temperature in frequency domain would effectively realign the mismatched time scales, as instigated for periodic unsteadiness by He and Oldfield. 17 For a scale-resolving turbulent flow solver (LES) based CHT, a general question of interest is, can we use a seemingly deterministic set of harmonics at their corresponding frequencies in a Fourier spectrum to represent seemingly "random" turbulence unsteadiness? To respond to the question, it should be first observed that turbulence energy is most commonly measured in the form of a Fourier spectrum (in computations and in experiments). If the spectrum is time-invariant, it should mean that each and every harmonic of spectrum should be time invariant. A harmonic in the context of a frequency spectrum of seemingly "random" turbulence may simply be interpreted as a statistically mean unsteady component at the frequency captured by Fourier transform over a long enough time period. In the same vein, it should not be surprising that a discrete Fourier spectrum is used for generating synthetic inflow conditions for LES, for example, Reference [23]. There have also been increasing applications of reduced order modeling (ROM) representations of turbulence, for example, POD 24 and DMD. 25 These ROMs are shown to be able to describe the majority of kinematics and dynamics of seemingly stochastic turbulence by using a relatively small number of deterministic spatiotemporal modes. The successful applications of the ROMs for turbulence flow analyses in themselves provide extra support for using a discrete Fourier spectrum in the interface treatment for LES based CHT. Once a continuous updating procedure of the Fourier spectrum is established, it can be applied locally at those fluid mesh cells adjacent to an interface, leading to a harmonic-balance interface treatment completely for the time-averaged part (zeroth harmonic), as well as for each and every unsteady harmonic retained in the Fourier spectrum, as shown by He. 22

2.2
Spatial two-scale framework for fluid and solid domains 2.2.1 Dual meshing by local embedding Figure 1 shows a surface with a large number of micro-structures. We start with a global coarse mesh which will be significantly under-resolved in the near-wall flow region around the micro-structures. If we locally refine mesh sufficiently for each and every one of all micro-structures, the adequate resolution will be achieved both locally and globally. This reference direct solution will be expensive, prohibitively so if we have a large number of micro-structures. The question of interest then is, are we able to get a global solution with an equivalent accuracy to the direct solution by solving only a small number of local fine-mesh blocks?
Here we adopt a dual mesh strategy ( Figure 2) where a small number of fine-mesh blocks are generated locally. The embedded fine-mesh blocks are simply generated by subdividing the corresponding coarse mesh cells. The intent is that the solutions of these few locally embedded fine-mesh blocks will provide a stepping-stone approaching to the target solution where all micro-structures are subject to their own local fine-mesh blocks respectively.

2.2.2
Upscaling by space-time averaging In the dual mesh system, the target solution for a coarse mesh cell is the locally spatial-averaged solutions of those fine-mesh cells embedded in this coarse-mesh cell. Furthermore, when scale-resolving turbulent solution is sought for the fluid domain, a time-averaging is introduced. As illustrated in Figure 3, the instantaneous flow field shows hardly any similarity among micro-structures ( Figure 3A). However, for the time-averaged flow ( Figure 3B), firstly note a spanwise similarity among different rows of the microstructures (in the vertical direction). Secondly, the local flow field around each micro-structure also becomes similar streamwise (in the horizontal direction) after the two upstream micro-structures. Thus, we take a working path through a time-averaged flow where the smooth variation of flow parameters makes it easier to propagate the fine-mesh solution effect from one micro-structure to its neighbors. Then the tasks become (a) to set a target for a coarse-mesh, and (b) to find a means of propagating the corresponding fine-mesh solution to drive the coarse-mesh solution toward the target. We consider a general governing equation for a fluid or solid domain.
R D denotes the discretized flux residuals which can result from different spatial discretization schemes of different orders of accuracy. Thus Equation (1) is in a semi-discrete form, to facilitate the introduction to the space-time averaging also in a discrete form. U is the vector of five conservative flow variables for a fluid domain. For a solid domain, the heat conduction equation is in the same form as the fluid energy equation with zero velocity and given solid material properties, U is then solid temperature.
Considering the space-time average for a coarse mesh cell M in a discrete form, we first take a local volume-weighted average. The instantaneous spatial averageŨ f for coarse-mesh cell M of total volume V M at time step n becomes, N m is the number of fine-mesh cells with volume of v m , (m = 1, … , N m ) contained in coarse-mesh cell M. The subscript "f " denotes that the variable is originated from a fine-mesh solution. The corresponding space-time averageŪ f for coarse-mesh cell M is thenŪ N t is the total number of time steps with a constant physical time step size being assumed. Equation (3) effectively defines the time-averaged target for a coarse-mesh when it is subject to a locally embedded fine-mesh. Then we need to find a way to drive the coarse-mesh solution toward the target when the coarse mesh has no local fine-mesh embedded.
In the context of multi-scale methods, the term "upscaling" is referred to as a means to formulate an augmented equation for a coarse-scale of the same form as that for a fine-scale. 26 In the present work, we upscale the governing equation To that end, the equation for coarse mesh variable U C can only be balanced with extra source terms: Equation (4) is effectively the upscaled semi-discrete coarse-mesh equation (fluid or solid). The source term S st (x) is time-invariant, consisting of two parts: The first part is generated in the fine-mesh domain by applying the space-time averaging to the fine-mesh solution: Equation (6) means that for a given fine-mesh solution, we can obtain this source term simply from the discrete flux residual directly computed using the space-time averaged variables on the corresponding coarse-mesh. The second part in Equation (5) is due to the nonlinear time-averaging effect of the unsteady solution on the coarse mesh which needs to be accounted for in order to balance the upscaled equation for the targeted space-time averaged solution: Given the source terms as defined above (Equations 5-7), we can now get the final time-averaged solution on the coarse mesh by time-averaging the upscaled unsteady governing equation (Equation 4). The resultant equation in a discrete form will be, Thus, the time averaged solution U C on the coarse-mesh should approach the target solution originated from the corresponding fine-meshŨ f as intended.
It may help the interpretation of the space-time averaging by noting two relevant scenarios. In the first one, when there is no spatial refinement, the space-time averaging is reduced to a time-averaging only. Equation (8) then becomes, For a fluid domain with a scale-resolving turbulent flow solution, this is equivalent to the Reynolds-averaging and the source term should reduce to the lumped Reynolds stress terms for the three momentum equations, and the extra heat and shear-stress work fluxes arising from nonlinearity for the energy equation. For an incompressible flow model, the source terms simply become a divergence of the Reynolds stresses for the three momentum equations respectively. Now the flow equations are effectively driven by the given Reynolds stresses toward a target time-independent solution of the time-averaged flow field, for example, by Ning and He in reference [27].
In the second scenario, when the problem is steady, the space-time averaging reduces to a spatial-averaging only for the coarse-mesh with an embedded fine mesh. Equation (8) then becomes, Consider more specifically steady conduction in solid. For given material properties, the conduction equation is purely linear. Still for the discretised fluxes, we will haveR D ( The non-zero source term is generated for this linear equation similarly to the time-averaging of the nonlinear flow equations when the Reynolds stresses are generated by nonlinearity. In the steady conduction problem however, the source is generated due to discretization truncation errors, so that the spatial-averaged solution from the fine mesh will not satisfy the steady conduction equation discretised on the coarse mesh. The source term (RHS of Equation 10) will thus be needed to balance out the truncation errors and drive the coarse-mesh solution U C toward the targetŨ f .
Having set up a dual-mesh system and identified a target solution for the coarse mesh, we can now think of a global target solution as a result of coupling the two sides of the two-scale system. The coarse-mesh part, which on its own would be poorly resolved, can be improved by the source terms originated from the fine-mesh blocks. On the other hand, the fine-mesh blocks, which on their own would be poorly conditioned, can become better conditioned with the improved global coarse-mesh solution. A remaining task is to propagate the source terms from those local fine-mesh blocks to the corresponding global coarse-mesh domain. The task is carried out by the block spectral mapping, to be described in Section 3.3.

2.3
Temporal two-scale treatment for fluid-solid interface between the two can be rearranged as 28 : For an order of magnitude estimate, take a typical stainless steel as the solid, air of Prandtl number Pr f at the ambient condition as the fluid, . We then have: , respectively. That means the total number of time steps for a fluid-only or a solid-only solution are, respectively: Given the time steps ratio based in the numerical stability and resolution is 17 we can then make following observations.
(i) For a steady CHT, the time scale disparity does not really matter. As there is no time consistency constraint, each domain can be time-marched completely in its own time and by using its own time step. Consequently a steady CHT will only need a number of time-steps comparable to its counterpart for a separate decoupled fluid-only solution or solid-only solution respectively: (ii) For an unsteady CHT, there will be a profound dual requirement 28 : (a) Time accuracy: small enough step size for fluid disturbances, Δt CHT = Δt f (b) Time consistency: large enough time scale to cover for solid domain, τ CHT = τ s . Without the long enough time period covered, initial transients may not be cleared, and convergence of the solid conduction part (thus the overall CHT solution) may not be warranted.
Therefore, the total number of time steps needed for the unsteady CHT becomes N unsteady CHT can thus be drastically increased (by four orders of magnitude) compared to that for a fluid-only or solid-only solution. Therefore, a direct time-consistent time-accurate unsteady CHT will be prohibitively costly.

Fourier-spectral CHT method
Given the modeling considerations discussed above, our priority is to avoid a direct time-consistent and time-accurate unsteady CHT solution. The overall approach taken here is to treat the unsteadiness at a fluid-solid interface ( Figure 4) in frequency domain. This effectively converts an unsteady CHT to two equivalently steady problems tasked with: (i) ensuring the continuity of time-averaged (zeroth harmonic) temperature and heat flux across the interface in the same way as for steady CHT; (ii) ensuring the harmonic continuity of unsteady temperature and heat flux across the interface by balancing each and every harmonic retained in the Fourier spectrum.
For variable ϕ (temperature T or heat flux q), a discrete temporal Fourier spectrum with N F harmonics is expressed as:

Coupling time averaged fluid with steady solid (zeroth Fourier mode)
For solid conduction, the temperature fluctuations in response to fluid unsteadiness tend to be small, so unsteady solid temperatures can be treated as linear perturbations about the steady temperature field. As such a much larger solid time step can be taken. In the cases of the LES subject to steady boundary conditions, the solid domain can simply be solved with the maximum local solid time step dictated by numerical stability. We first time-average the related variables on the fluid side for the mesh cells adjacent to the solid wall ( Figure 4) and then couple the time-averaged fluid parameters with those of a steady state on the solid side. The time-averaged steady flux continuity, when discretized with one-sided finite difference, should lead to an explicit expression of the wall temperature in terms of the time averaged parameters on the fluid side and the steady ones on the solid side, The time-averaging of the variables at the fluid mesh cell adjacent to the interface can be simply and effectively carried out in a moving-average, updated from the one at the previous time step and the current solution. 22

Temperature transfer function and harmonic balance for unsteadiness
For unsteadiness cast in the frequency domain, the flux continuity (energy conservation) can be achieved by balancing all N F harmonic components respectively.
For a local 1D unsteady conduction behavior for a semi-infinite solid domain, similar to that in well-established data processing methods in transient heat transfer experiments, 30 there is an analytical relation between the wall temperature and the heat flux: The complex number  17) then leads to a semi-analytical expression of the wall temperature harmonic response: The wall temperature harmonics (Equation 19), and the time-mean (zeroth harmonic) wall temperatures (Equation 16) are then used to construct the instantaneous unsteady wall temperature in time to proceed with the time-domain fluid LES solution. As such, the Fourier spectral interface coupling effectively provides a Dirichlet condition at every time step of the nonlinear time-domain scale-resolving flow solution.
Remarkably, thanks to the semi-analytical harmonic transfer function, we can now get unsteady wall temperatures but in so doing without even solving the unsteady conduction equation in the solid domain at all. This is simply because the T-q transfer function (Equation 18) already analytically includes the unsteady interaction between the wall flux and temperature in the solid domain. Consequently, the solid domain can be solved completely as a steady conduction problem with a local time step size several orders of magnitude larger than that of the time-domain.
More details for the frequency domain unsteady interface treatment, particularly those detailed formulations of the T-q harmonic transfer function and the moving-average discrete Fourier transfer of near-wall fluid temperatures can be found in Reference [22].

METHOD IMPLEMENTATION
The flow solver is based on an in-house multi-block code which the author has developed and updated over many years. 9,29,31 The equations are spatially discretized using a cell center based finite volume scheme. The spatial discretization for the advective fluxes is based on the pressure-convective flux split upwinding discretization, AUSM, Liou. 32 The solution method uses a 4-stage Runge-Kutta scheme to explicitly integrate the spatially discretised equations in a pseudo time , accelerated by the local time stepping and the multi-grid techniques. The physical temporal terms for unsteady flows are discretized by a 2nd order backward difference and solved in a dual time stepping.
For the solid domain, the temperature field is governed by the steady heat conduction which effectively is in the same form as the fluid energy equation, with zero velocities and corresponding material properties. Again, the equation is in a similar simple form to the equations for the flow domain and solved by the same time marching method in the pseudo time. Further details of the baseline methods can be found in Reference [22]. A few distinctive aspects specific to the present work are described in the following.

Scale-dependent coarse-fine mesh interface treatment
For a locally embedded fine-mesh block as shown in Figure 2. Given that the fine-mesh is embedded in a coarse one, a coarse-fine mesh interface will be subject to hanging-mesh notes. A key requirement is the total flux conservation. This is met by overwriting the total fluxes across a coarse-mesh face of a boundary mesh cell by the sum of total fluxes of the corresponding fine-mesh cell faces. The interface treatment needs to capture the physical interactions between the two mesh domains, as enabled by the respective mesh resolutions respectively. More specifically, we recognize that the global coarse-mesh domain provides the flow conditioning at its low mesh resolution level. The fine mesh block needs to be directly interfacing with and influenced by the flow variables of the coarse mesh side of the interface. But at this same time, the fine-mesh side should avoid being completely fixed by the coarse-mesh side with a much lower resolution. A scale-dependent interface treatment 9 is thus adopted.
Express a variable at a fine-coarse mesh interface ϕ as a base coarse-mesh variable ϕ C and a local fine-mesh disturbance Δϕ f on top of the base value.
Consider a pair of interface mesh points, for example, "1" and "2", between the two spanwise side boundaries of the fine-mesh block ( Figure 2). First, the coarse mesh variables ϕ C (x, t) will be directly interfacing locally between the two sides of the local interface. Thus values of ( ϕ C ) 1 and ( ϕ C ) 2 are determined by the direct interfacing at the local interfaces respectively. For the fine-mesh disturbance part however, a direct periodic condition is applied: Thus the full variables at the two paring points are updated in Equations (22) and (23), respectively.
The scale-dependent interface treatment can also be applied to the streamwise pairing points (1 ′ and 2 ′ ) if the local flow is fully developed, as illustrated in a locally embedded two-scale LES solution for canonical turbulent channel flows. 10

Balanced eddy viscosity (BEV) damping for source-term driven flow equations
It is recognized that the solutions to flow equations with added source terms may experience convergence issues.
Recently, Wu et al. 33 have carried out a stability analysis of the steady RANS equations of an incompressible flow, and identified that the equations with explicitly added Reynolds stresses can be ill-conditioned. Much of the discussion in Reference [33] is on contrasting an implicit formulation of the Reynolds stresses with an explicit one for the improved conditioning. From the present author's viewpoint, the conditioning improvement might perhaps have more to do with the turbulent eddy viscosity being introduced in the reformulation than the Reynolds stresses being added implicitly or explicitly. Closely relevantly, the role of the eddy viscosity in stabilizing Reynolds stress driven solution is also reported earlier by He 9 for a reformulated mixing-length eddy-viscosity model, and Zhang and He 34 for a reformulated one-equation S-A model. The common link is believed to arise from the basic property that the eddy-viscosity tends to dampen out disturbances (physically or numerically originated) which may otherwise grow into instabilities.
For the present source term driven solutions, the balanced eddy viscosity method (BEV), originally introduced and tested for unsteady Reynolds stresses source term driven compressible RANS solutions 9 is followed. The preference is made not just because of the author's own previous implementation and operational experience of the method. Compared to the one-equation model of Zhang and He, 34 the explicit mixing-length model is simpler and more efficient computationally. The implicit eddy viscosity formulation of Wu et al. 33 needs an inversion of the mean flow strain for a simple adoption, which may cause numerical difficulties for complex flow with nearly zero local strains. The difficulty can be avoided by fitting an eddy viscosity field through the minimization process as adopted in Reference [33]. A field minimization nevertheless is not a trivial operation. For a pure data-driven problem with a given Reynolds stresses field, the minimization for fitting the eddy viscosity for the whole flow field can be pre-processed prior to the RANS solution. In the present work however, we have a moving target as the source terms needs to be updated during the solution to capture the interaction/interdependence between the two scales/domains. Thus, a locally updated eddy viscosity formulation which can be simply updated with a minimal overhead during the solution is preferred. Moreover, the BEV should have no effect on the final time-averaged (steady or unsteady) solution when approaching the space-time averaged target originated from the fine-mesh solution.
Consider a simple eddy-viscosity model based on the target time-averaged flow propagated from the fine mesh.
In the previous work, 9 the maximum mixing-length l was prescribed as a fraction of the characteristic length of flow domain, and the insensitivity to the empiric input was demonstrated for steady compressible RANS solutions driven by unsteady stress source terms. In the present work, a length limit is taken relative to the local mesh-size similarly to a sub-grid scale model, thus: where k is the von Karman constant, 0.41, y is the distance to the wall. l SGS is the local length scale for a sub-grid scale model. For the Smagorinsky model, l SGS = C SGS Δ mesh (a nominal value, C SGS = 0.16). Δ mesh is the local mesh size, C ev is the coefficient to control the eddy-viscosity damping. In this way, the user can conveniently select the mixing-length limit relative to the local mesh size, commonly used as the length scale in a sub-grid scale model for LES. We now add the eddy viscosity to the upscaled coarse mesh flow equation (Equation 4).
where LD denotes the discrete lumped linear dissipative terms. For the three momentum equations the frictional stresses terms included in LD all have the eddy-viscosity as a time-invariant coefficient. For the energy equation, the linear dissipative terms only include those for heat conduction with a given constant Prandtl number and specific heat, while the shear work-done terms by frictional stresses, being nonlinear, should not be included in LD.
Time-averaging Equation (26) and noting LD being time-linear, thus LD (U C ) = LD , we will then have a resultant discrete equation: Thus, the discrete equation (Equation 27) can only be satisfied (numerically balanced) when the converged time-averaged discrete solution U C on the global coarse-mesh approaches the discrete target solutionŨ f , as intended. As such, the eddy-viscosity acts to stabilize the coarse-mesh field during the solution process through LD (U C ). The net-effect of the added eddy-viscosity on the final time-averaged solution will nevertheless be balanced out. It is thus called "balanced eddy viscosity" (BEV). This important property will also be further tested and demonstrated in the result Section 4.2.1.

Block-spectral mapping
A key enabler in facilitating the mutual interactions between the coarse-mesh and fine-mesh scales is to propagate the equivalent high-resolution target from the small number of locally embedded fine-mesh blocks to the global coarse-mesh domain. This involves a construction of a spectral variation and then map the target to the global domain accordingly. For a surface with large number micro-structures (Figure 1), the similarity in the geometry and the pattern of time-averaged (instead of instantaneous) flow among microstructures ( Figure 3) means a smooth variation can be obtained if a spectrum is established at the same position relative to each micro-structure. Thus, when fine-mesh blocks of the same topology and connectivity are used for those micro-structures, a mesh-pointwise variation among blocks can be accurately represented by a block-to-block spectrum for the same mesh point in all blocks, which is thus called "Block Spectrum." Furthermore, for a micro-structured surface, a block spectral propagation will only need to be carried out in the two wall parallel directions.
The detailed formulations for a block-spectral mapping based on a double Fourier series can be found in Reference [9], which are not to be repeated here. A key point to note is that once the coarse-mesh target solution is identified from the locally embedded fine-mesh blocks, the block spectral method can used to map any variables from the local fine-mesh to the global coarse-mesh. The first option is to propagate the space-time averaging source terms (S st ) f (Equation 6) generated from the fine-mesh solution.
Alternatively, we can first directly map (Ũ f ), the target space-time averaged flow variables from the local fine-mesh blocks, to the global coarse-mesh, and then work out the corresponding source terms locally on the coarse-mesh with the mapped (Ũ f ). Although for the simple mesh tested, both mapping options have similar results, the second option can potentially make it easier to extend the block-spectral mapping to non-uniform (or even unstructured) meshes.
In addition, it should be noted that the block spectral mapping is applicable to the surface fluxes. The coarse-mesh boundary fluxes overwritten by those of the corresponding embedded fine-mesh can also be mapped to other unembedded coarse-mesh boundary surfaces. The flow chart of the time-marching solution process in one time-step for the coupled fluid-solid two-scale system is given in Figure 5.

RESULTS
Computational case studies are carried out principally for two main tasks: 1. To validate the conjugate heat transfer capability against experimental data for relevant configurations with extended surface structures; 2. To validate the two-scale method with fine-mesh blocks only for a small number of micro-structures against the direct full solution with fine-mesh blocks for all micro-structures.
In addition, the results of some cases can also further help elaborate several aspects relevant to the present framework and method implementation: (i) The interaction between the local fine-mesh and global coarse-mesh, how can we tell? (ii) The wall temperature unsteadiness in CHT, how does it manifest? (iii) The BEV, does it affect the time-averaged coarse-mesh solution result?

4.1
Validation of conjugate heat transfer for internal ribbed channel

Case setup
The VKI internal flow experiment extensively studied by Casarsa and Arts 35 and Cukurel et al. 36 is judged to be a suitable case to validate the CHT capability of the present development. It is a simplified gas turbine blade internal cooling configuration of a rib-roughened cooling flow channel of a square cross-section. The configuration consisting of six ribs is taken ( Figure 6A). In the experiment, the ribs are heated from the bottom surface with a constant flux heater, so a constant heat flux condition is applied. For the top and side walls, the adiabatic condition is applied. The key non-dimensional parameter to be matched is Reynolds number, Re Dh = 40,000 based on the bulk velocity and the cross-section hydraulic diameter D h , The rib height H is 0.3 D h . and the rib-rib (pitch) distance is 3 D h . The computations are carried out at Mach number around 0.2, to avoid the convergence difficulty for a density based compressible flow solver like the present one, similarly to that of Scholl et al. 37 At the inlet to the domain, stagnation pressure profile is specified to produce a turbulent velocity profile (1/7th law based on the local distance from wall). No turbulence fluctuations are introduced at the inlet, as a fully turbulent flow pattern is seen to be quickly reached at the second rib, insensitive to the detailed inlet velocity profile as discussed in Scholl et al. 37 The low sensitivity of local aerothermal behavior to the detailed flow conditioning of ribs further upstream should also support locally refining only the mesh for the 3rd, 4th, and 5th ribs where the aerothermal parameters for the 4th rib will be closely examined with its immediate neighboring 3rd and 5th ribs providing required up−/down-stream conditioning. The base coarse-meshes for these three ribs are shown for the fluid and solid domains in Figure 6B,C, respectively.
All the solutions are run several throughflows when well established unsteady flow patterns emerge. It is observed that time-averaged parameters of interest tend to become sufficiently constant after being time-averaged over 5-6 throughflow periods. For all heat transfer results, the local Nusselt number (Nu) is the local heat flux nondimensionalised by local fluid thermal conductivity, hydraulic diameter D h and the temperature difference between the wall and the bulk flow linearly interpolated between the inlet and exit. Furthermore, the effectiveness of internal ribbed channels is measured by the effective enhancement (EF), Nu∕Nu 0 , (Nu 0 = 0.023Re Mesh sensitivity and validation for flow velocity and heat transfer Given the importance of mesh resolution as part of the primary intent of the present two-scale method development, we shall first get some indication of the mesh sensitivity. With the base coarse mesh ( Figure 6B,C), we now look at the different levels of refinement by embedding different numbers of fine-mesh cells in a coarse mesh cell. Table 1 lists the mesh cell counts for the base mesh only case and those with different embedded fine-meshes. Here the cases denoted as "Embed8," "Embed12," "Embed18" are for those with 8, 12 and 18 fine-mesh cells embedded in one base coarse-mesh cell respectively. Figure 7 shows the Nu distributions along the wall surface coordinate (s) around the 4th rib. Note that s = 0 is taken at the middle point of the top surface of the 4th rib. Therefore, −1.5 < s/H < −0.5 is the left sidewall of the rib and 0.5 < s/H < 1.5 is the right sidewall, as viewed in Figure 6A.
It is clear that the base mesh is very poorly resolved. Embed8 shows a reasonable mesh convergence for large part of the whole surface around the rib, but it is still clearly under resolved around the rib top surface (−0.5 < s/H < 0.5). The results for the finer mesh embedding, Embed12 and Embed18 indicate the range of the refinement required for an adequately mesh independent solution.
The mesh sensitivity for the corresponding flow field is illustrated in Figure 8A for the mean velocity and Figure 9A for the unsteady fluctuations in RMS, where velocity profiles in the vertical direction (y) on the mid-span plane are taken at different streamwise positions (s). Overall, we can ascertain that the Embed12 and Embebd18 provide far less mesh dependency than Embed8 for both the time mean and fluctuation parts than Embed8, as expected. More remarkable however, is that the least mesh sensitivity for velocity is observed in the region right above the top surface of the rib ( Figures 8A and 9A), which happens to have the highest mesh sensitivity for the heat transfer ( Figure 7). The observed difference in mesh sensitivity between the heat transfer and velocity field may seem unusual from the viewpoint that fluid flow should dictate convective heat transfer. It can be argued that a mesh resolution for getting a good velocity field (easier to compute and to measure) may not be good enough for getting a right gradient/differential on the wall for the local heat flux. The corresponding comparisons with the experimental data for the time-mean and fluctuation velocities are given in Figures 8B and 9B, which can be regarded as in good agreement.
The present heat transfer result in comparison with the experimental data may be regarded as reasonably satisfactory ( Figure 10). For this case, there are a few published computational heat transfer (Nu) results with a quite considerably TA B L E 1 Mesh count for VKI internal ribbed channel (embedded fine mesh for ribs 3,4,5).

Mesh-cell counts
Base mesh 1.5 million varying degree of differences from the experimental results. A very comprehensive case study with detailed comparisons for this case is made by Scholl et al. 37 Their CHT heat transfer result is also included in Figure 10. It is fair to say the present CHT result follows closely that of Scholl et al. 37

Wall temperature unsteadiness sensitivity to interface treatment
We now take a look at the impact of the fluid-solid interface treatment on wall temperature unsteadiness. The time scale disparity has been seen to cause an extra solid-domain mesh dependence of unsteady wall temperature fluctuation, which can be over-predicted by an order of magnitude. 21 For the present experimental case, the solid conducting part is made of AISI304 steel. We compare the unsteady wall temperatures from two solutions: the first is the present developed Fourier interface treatment where 100 harmonics (from 100 Hz to 10,000 Hz) are retained in case of any significant unsteadiness may emerge in this frequency range. The second is the time-domain direct coupling as commonly available in commercial CFD codes and often used for CHT simulations. It must be noted that the time step used in the solid domain is much larger, so the time-domain direct coupling is not time-accurate. The much larger solid time step is necessarily needed for the solution convergence in the solid domain, as discussed in Section 2.3.1. Should a time-accurate and time-consistent CHT solution be adopted, the solid domain would have had to be time marched in the same fluid time step size, requiring a factor of 10 3 ∼ 10 4 more time-steps over that for a fluid-domain only solution (see Equation 14), practically out of reach. If on the other hand, the CHT system is time-marched in the small fluid time step through only a short time period like that typical for a fluid-only solution, the solid domain would still be subject to initial transients, a long way away from a credibly converged solution. Figure 11 shows the time traces of instantaneous wall temperature relative the time mean value at the middle point on the rib top surface (s/H = 0, Figure 6A). The unsteady wall temperatures are normalized by the mean temperature difference between the inflow and wall, which is roughly the mean driving temperature difference. We can see that the commonly available time-domain interface treatment can produces very large wall temperature fluctuations, up to 15%-20% of the fluid driving temperature difference. On the other hand, the present Fourier interface treatment picks up very small unsteady wall temperature fluctuations. The time-domain direct coupling interface overpredicts the unrealistically high wall temperature unsteadiness by ∼50 times over that by the present Fourier-spectral method. Table 2 compares the time mean and unsteady RMS heat fluxes corresponding to 10,000 time-steps. Here, we can see that the time-domain interface method underpredicts the heat flux RMS, in line with that expected from the wall temperature fluctuations, as shown in Figure 11. However, there is little difference in the time-averaged values. Thus, for the present case, if we TA B L E 2 Time-mean and unsteady (rms) heat flux predicted by two interface methods. are only after the time-averaged heat transfer, the time-average (zeroth harmonic) based on the simple moving-average should suffice. In many cases of practical interest, the time-averaged part should be of primary importance, and thus should serve as the baseline CHT solution. The Fourier spectral method provides a simple and consistent way to add on the unsteadiness with very small extra computational cost. For this case, the Fourier-spectral unsteady CHT with 100 harmonics adds only a small overhead (∼20%) over the baseline steady (time-averaged) CHT, while providing the capability to be able to capture the wall temperature unsteadiness, which is generally completely missed in other CHT methods. For practical CHT analyses, one may retain only a small number of harmonics initially to just probe the magnitude of the unsteady wall temperature with minimal extra computational effort. The number of harmonics can be expanded during the solution, should it be needed. It is worth noting that although in many cases, the wall temperature unsteadiness tends to be small, it may vary considerably with solid material. 21,22

Two-scale solutions for micro-structured surface
The test configuration considered is a steel slab with 100 micro-structures on its top surface subject to an inflow in the x direction ( Figure 12). We will first look at the computational tests for the fluid-domain only and the solid-domain only, respectively, to assess the validity of the two-scale block-spectral methods with locally embedded fine-mesh blocks for a subset of micro-structures. The two-scale solutions will be compared with the direct full solution where the fine meshes are applied to all micro-structures. In addition, the coupling between the local fine-mesh blocks and the global coarse-mesh domain will be further illustrated by comparing the block-spectral two-scale solutions with those solutions with locally embedded fine-mesh blocks but without the source term coupling (denoted as "One-scale" solutions).

Turbulent flow solution over micro-structures (fluid-domain only)
For the flow domain, most boundary conditions are the same for this case and the fluid-solid coupled conjugate heat transfer case to be presented later. The key difference is that for this case, the adiabatic wall condition is applied to the bottom surface of the domain subject to micro-structures. The inflow is in the x direction with specified stagnation temperature of 300 K, a bulk flow stagnation pressure of 100,000 pa. The exit static pressure of 98,000 pa is specified to keep the flow in an incompressible flow regime. All 100 micro-structures are of the same simple shape with a square base and a height/width of 0.5. The upper halves of the two side surfaces are tapered to the top surface by 10% of the width, as shown in the closeup of the micro-structure element ( Figure 12). The Reynolds number based on free stream velocity and micro-element height is 2800. At the inlet, the stagnation pressure is reduced toward the bottom wall to mimic an endwall boundary layer profile measured in a wind tunnel 38 where the test section is after a 7.5:1 contraction, and thus is deemed to be of low free-stream turbulence. The boundary layer thickness is taken to be about seven times of the height of the micro-element. In addition, a spanwise nonuniformity is introduced in stagnation pressure distortion in the form of, where Yp is the total spanwise length of the domain. The amplitude A p01 corresponds to a spanwise stagnation pressure distortion of 15% of the inlet dynamic head in the bulk flow region outside the inlet boundary layer, linearly reduced to zero at the bottom wall.
Given the efficient mesh refinement as the central theme of the work, we first need to get some sense of the sensitivity for direct solutions before assessing the two-scale method. The base mesh is very coarse with only three mesh cells to cover the micro-element width and two cells for the height, so it clearly will be poorly resolved. It is of interest to see if we can achieve a two-orders of magnitude refinement around each micro-structure element. Two fine-mesh embedding in coarse mesh are thus examined: one with 5 × 5 × 5 fine mesh cells in each coarse cell (denoted as "Embed125") and the another finer one in the x and y directions with 6 × 7 × 5 (denoted as "Embed210"). The time averaged velocity traverses are taken at the location downstream of the final row of micro-elements (90% domain length) in the streamwise x direction and at the element height in the vertical z direction. The results for Embed125 and Embed210 are in good agreement, and both are in a clear contrast to the poorly resolved base mesh result as shown in Figure 13. Consequently, the embedded refinement with the mesh density ratio of 1: 125 is deemed sufficient for the present case and will be used in the following tests.
Now we examine the two-scale solution capability. Figure 14 indicates the scale resolving capability of the all fine-mesh direct solution compared to the two-scale solutions. In the direct solution ( Figure 14A), all 20 spanwise rows of micro-elements are refined by the fine mesh blocks (Embed125). In the two-scale solution however, the fine-mesh blocks are only embedded in a region covering three spanwise rows for the 3rd, 4th, 5th row, for the 10th, 11th, 12th row, and for the 17th, 18th, 19th row, respectively ( Figure 14B). Here the instantaneous flow patterns are given in iso-Q plot colored by the x velocity contours. We can see that the flow passing around the first streamwise row of the micro-element become unstable with shed vortices. These vortices develop into quite coherent hairpin vortical structures around the second streamwise row and in turn start to interact more strongly spanwise. Further downstream, the interactive flow structures become less coherent toward the end of the domain. Note that the locally embedded three-row fine-mesh solution appears to capture well the whole streamwise development of the unstable vortical flow structures. But more relevantly, we need to know if these partial fine-mesh solutions can serve to fill the coarse-mesh gap regions in between properly. The three-row embedded fine mesh block is taken to reduce the unideal interface influence on the source terms, as observed previously in analyzing multiple jets. 9 At each time step, the source terms are first generated from row 4, 11, and 18. The mesh pointwise spectrum is constructed and then propagated to other unembedded coarse-mesh regions. The impact of the source term coupling is shown to be very effective (Figures 15 and 16).
A final test case for this section is on the sensitivity of the time-averaged flow solution to the BEV to enhance convergence of the source term driven coarse-mesh solution, as introduced in Section 3.2. In the present BEV model, this mixing-length limit is measured in terms of the local mesh size similarly to the standard sub-grid scale model. The effect of the eddy viscosity is controlled by a multiplier coefficient Cev (Equation 25).
For most cases tested, the coarse-mesh solutions driven by the source terms seem to converge well to the target solution. There are few cases where certain local flow regions exhibit poor convergence behavior, drifting away from the space-time-averaged target flow, although the corresponding coarse-mesh solutions do not show any full-blown instabilities. For these cases, introducing the BEV is shown to be very effective in improving the convergence. Figure 17 shows an example of the influence of the BEV damping. In the case, the base solution without the BEV (Cev = 0) indicates a poor convergence where the flow entropy variable is drifting persistently about 30% different from the space-time averaged target solution. The solution is then restarted with the BEV introduced with Cev = 1, such that the mixing-length is limited to the local mesh size as in the Smagorinsky sub-grid model (Equation 25). Figure 17 shows how the entropy variable E changes with time steps when the BEV damping is introduced. The variable is expressed in |E − E av |/E 0 − E av |, where E 0 is the value at the restart (N = 0), E av is taken from the final converged time-averaged solution. We can see that the BEV improves the convergence markedly in just a few hundred time-steps.
Next, we assess the sensitivity of the final solutions to the BEV. Two cases with two damping values are compared. The first one is for Cev = 1. In the second case (Cev = 4), the maximum mixing-length is capped at four times of that of the SGS model. It means that the eddy viscosity for the second case is 16 times higher than that of the first (see Equations 24 and 25). The time-averaged velocity traverses in spanwise (y) direction at the 90% domain length and the micro-element height are shown in Figure 18. The excellent agreement between the two solutions confirms that the final time-averaged

Conduction solution for micro-structures (solid-domain only)
An examination on the solid domain solution is conducted before moving on to the CHT solution when both fluid and solid domains are coupled. This is needed as a verification given the new extension of the two-scale blockspectral mapping to the solid conduction solution. Consider the solid material of stainless steel with density 8000 kg/m 3 , specific heat capacity 500 J/kg/K, and thermal conductivity 15.7 W/m/K. The solid domain configuration and boundary condition setup are shown in Figure 19. The 100 solid micro-structure elements are of the same shape as the in previous case (Section 4.2.1) over the top boundary of the solid domain. The outer surfaces of the 100 micro-elements are subject to a specified wall temperature with a spanwise variation in a form of, where T 01 is the mean temperature over the micro-element surfaces, also taken as the reference temperature (T 01 = 280 K). A Tw is taken to correspond to 10 K amplitude in the spanwise temperature variation. The bottom boundary of the solid domain is subject to a specified low temperature T Wbottom = 0.6 T 01 . Effectively the micro-elements act as 100 discrete heaters. The adiabatic condition is applied to all side surfaces of the domain and to the unheated flat parts of the top surface.
The temperature contours on the top surface plane cutting through the bottom faces of the micro-elements ("A-A" plane, the zoom-in view in Figure 19) are shown in Figure 20. For the direct solution ( Figure 20A), the fine-mesh blocks are embedded for all 100 micro elements. For the two-scale solution ( Figure 20B), only are nine micro-elements are locally embedded with the fine-mesh blocks of a fine to coarse mesh density ratio of 1:125 (Embed125). The source terms are generated from these nine fine-mesh blocks and mapped to other micro-structures regions with only the base coarse mesh. Also included is the "one-scale" solution ( Figure 20C), which is directly solved with the same mesh as the two-scale one but without the source-term coupling.
We can see that the coupled two-scale solution ( Figure 20B) is in good agreement with the direct full fine-mesh solution ( Figure 20A). The large discrepancies between the two-scale solution and the one-scale solution ( Figure 20C) underscore the inadequate resolution of the base mesh and the marked improvement that the two-scale coupling can make. Particular attention is drawn to those embedded fine-mesh regions of the nine micro-elements (indicated by the solid lines in Figure 20B,C). Here both the two-scale and the one-scale solutions have the same mesh. But the differences in the local temperature distribution between the two solutions are still quite marked. The contrasting results underscore the principal working of the present two-scale methodology in harnessing the important interaction between the local fine-mesh blocks and the global coarse-mesh domain. Without the source terms coupling, the fine and coarse domains can only communicate at the fine-coarse mesh interfaces. Then the coarse-mesh solution, though well-conditioned, will be poor due to its under-resolution. Similarly but also in contrast, the fine-mesh solution, though well-resolved, would also be poor due to the inadequate conditioning from its surrounding ( Figure 20C). When the two domains are made to interact through the source terms coupling, the two domains correct each other which in turn help correct themselves, respectively ( Figure 20B), as intended.

4.2.3
Conjugate heat transfer solution (fluid-solid coupled) Having examined the fluid-only and solid-only cases respectively, we now put the two domains together for conjugate heat transfer of a micro-structured surface. The fluid domain retains the same boundary conditions as the flow-only case except for the bottom boundary surface, which is now subject to a fluid-solid conjugate heat transfer interface. Similarly, the solid domain retains the same boundary conditions except for its top boundary face now subject to a solid-fluid CHT interface.
For the two-scale solution, the fluid domain has the fine-mesh "Embed125" embedded in three regions each covering three spanwise rows as in the flow-only case. The source terms are generated from the middle row: the 4th, the 11th and the 18th row for the three fine-mesh regions respectively. Correspondingly, the fine-mesh blocks are embedded in the 4th, the 11th and the 18th row of the solid domain respectively. The source terms from the three rows in both domains are used to construct the mesh-pointwise block-to-block (row-to-row) spectra which propagate the source terms to other un-embedded coarse-mesh regions. In addition, the time-averaged fluxes on the wall boundary faces of the embedded coarse mesh cells are also mapped to other corresponding coarse-mesh boundary faces to drive the local time-averaged fluxes to the targets.
For the fluid thermal field, we compare the stagnation temperature contours on the cut plane at the middle height of the micro-elements, as shown in Figure 21. The two-scale solution ( Figure 21B) is again compared to the direct solution with all rows subjected to the embedded fine mesh ( Figure 21A). Also plotted is the "one-scale" solution without source term coupling ( Figure 21C). The corresponding stagnation temperature traverses at 90% domain length are shown in  Figure 22. The source term coupling is demonstrated to be very effective in driving the coarse-mesh solution to the target, without which the coarse mesh regions are clearly very poorly resolved ( Figures 21C and 22).
The time-averaged temperature contours on the solid surface cutting through the bottom planes of the micro-structures ("A-A" cut plane, Figure 19) are compared in Figure 23. The results consistently demonstrate the validity of the present two-scale block-spectral CHT method ( Figure 23B) in driving the coarse-mesh aerothermal field to the target solution ( Figure 23A). Without the source terms coupling, not only is the coarse-mesh region poorly resolved, even the fine-mesh regions (the 4th, 11th and 18th row) also have large errors due to the inadequate conditioning  Figure 23C). The result comparison between the two-scale and the uncoupled one-scale solution again highlights the effectiveness of the underlying working of the present method in addressing the basic conflict between a poorly resolved global domain and a poorly conditioned local one.
Finally, a general issue of interest is the interaction between heat transfer and flow. This may manifest in terms of the local Reynolds number, based on local flow parameters and the micro-element height h (Re local = ρuh/μ). The traverses of Re local at 90% domain length and mid micro-element height are shown in Figure 24. The conjugate heat transfer and the adiabatic solutions are contrasted. The local Reynolds number of the CHT solution is shown to be on average higher than that of the adiabatic case by ∼30%. The increase in the local Reynolds number due to a wall cooling is also observed to correspond to a rise in heat transfer coefficient 13 and an earlier laminar-turbulent transition. 15,16 Interactions between flow and convective heat-transfer in general deserve further investigations.

CONCLUDING REMARKS
The main objective is to develop an efficient and accurate two-scale conjugate heat transfer method for micro-structured wall surfaces. The key challenges faced are the scale disparities in space as well as in time.
To address the spatial scale disparity issue, a dual meshing is employed to facilitate the coupling between a global coarse-mesh and a small number of local fine-mesh blocks. The pathway taken is first to set up a target solution for the coarse-mesh by space-time averaging the counterpart fine-mesh solution. Then, propagate the target from the small subset of local blocks to the global coarse-mesh domain. This is achieved efficiently and accurately by a mesh pointwise block-spectral mapping of time-invariant source terms harnessing the similarity of time-averaged local flow patterns around micro-structures. The convergence of the source-term driven coarse-mesh solution is further enhanced by a BEV damping. This source-term based approach provides an effective two-way coupling. The global coarse-mesh solution corrects the local fine-mesh solution by improving its surrounding environment/conditioning, while the local fine-mesh solution corrects the low resolution of the global coarse-mesh by providing the source terms. In so doing, both domains effectively correct themselves. As a key enabler for the present work, the two-scale method is extended to the solid domain, so that the detailed local high gradient temperature distributions around large number micro-structure elements can now be resolved efficiently and accurately.
For the time scale disparity in CHT, a present work is based on a temporal Fourier spectral framework to avoid time marching the fluid and solid domains with hugely mismatched time scales. The time-averaged (the zeroth harmonic) is effectively obtained in the same way as a steady CHT aided with a moving-average. Wall temperature unsteadiness is simply and explicitly obtained from the fluid temperature harmonics through a semi-analytical wall temperature transfer function. As a result, the solid domain can be solved completely as a steady conduction problem in its own time step which is by several orders of magnitude larger than the physical time step applicable to the fluid domain.
The developed CHT capability is validated for an internal cooling channel with a ribbed surface. For a test configuration with 100 micro-structures, a fluid domain-only, a solid domain-only and a coupled CHT solution are examined respectively. These test cases with two orders of magnitude local embedded mesh refinement consistently demonstrate the validity and potential effectiveness of the present framework and implementation methods. Some of the results also clearly illustrate the principal interactive working of the methodology.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request. ORCID L. He https://orcid.org/0000-0002-6791-809X