AN ITERTIVE DESIGN METHOD TO REDUCE THE OVERALL THERMAL RESISTANCE IN A CONJUGATE CONDUCTION-FREE CONVECTION CONFIGURATION

A design approach is proposed and demonstrated to identify desirable two-dimensional solid geometries, cooled by natural convection, that offer superior thermal performance in terms of reduced overall (conduction-convection) thermal resistance. The approach utilizes (i) heat transfer modeling in conjunction with (ii) various novel shape evolution rules. Predictions demonstrate the evolution of the solid shape and associated reduction of the overall thermal resistance. Parametric simulations reveal the dependence of the predicted solid shape on the evolution rule employed, the thermal conductivity of the solid material, and the strength of advection within the fluid.


INTRODUCTION
The design of various thermal systems has long been a point of research with a common goal being to reduce the overall conduction-convection thermal resistance while minimizing the weight or volume of the solid. Relative to heat sink design, early on this was achieved using size and geometrical optimization of heat sinks that were comprised of geometrically well-defined sub-components (e.g. straight, circular, or pin fins of various cross sectional shapes). For example, Bar-Cohen (1979) developed expressions for optimum fin thicknesses, spacing between fins, and fin widths, that maximize the heat dissipation per unit fin cross sectional area. The correlations were developed based on the assumption of a uniform heat transfer coefficient and fin efficiency. In a later study by Bar-Cohen et al. (2003), the work was extended to incorporate experimentally validated correlations for the heat transfer coefficient (Nusselt number), including for non-isothermal plates. In a similar study, Kim (2012) optimized the size and shape of plate fin heat sinks that had a trapezoidal cross-section and found better thermal performance when compared to a similar optimized plate fin heat sink of uniform thickness.
While size and shape optimization provided good designs of traditional heat sink layouts comprised of, for example, plate fin or pin fin arrays, a truly optimal heat sink shape that achieves the best thermal performance is often complex and might consist of unanticipated geometrical features. To this end, topology optimization (TO) is an evolving design methodology which provides the best solid shapes that will achieve a specified goal while satisfying imposed constraints (Bendsøe and Sigmund, 2003).
The general TO process involves: (i) discretization of a computational domain consisting of a solid and perhaps an adjoining fluid, (ii) specification of an initial spatial distribution of solid material, (iii) solution of the appropriate equations that describe the physical process or processes of interest, and (iv) optimization of the solid material shape and configuration using an appropriate algorithm. Steps * Corresponding author. Email: cdsevart@ku.edu (iii) and (iv) are repeated until the final solid shape is achieved. A comprehensive review of the different variations of TO and their characteristics can be found in Sigmund and Maute (2013).
TO was first introduced by Bendsøe and Kichuchi (1988) and used as a method to most effectively distribute stresses in a structural member of unknown shape but has since been adapted to the fields of heat transfer and fluid flow. To this end, Li et al. (1999) applied TO to heat conduction problems; they used the principles of evolutionary structural optimization (ESO) to design solids that resulted in the most uniform distribution of a surface heat flux. Other works that apply TO to heat conduction problems include Gersborg-Hansen et al. (2006), Gao et al. (2008), Marck et al. (2012), Dirker and Meyer (2013) and Xia et al. (2018). Subsequently, several studies included convection heat transfer occurring at the boundary of solid structures into TO problems by use of a constant convection coefficient and Newton's law of cooling (Yin and Ananthasuresh, 2002;Bruns, 2007;Ahn and Cho, 2010).
To the authors' knowledge, the first investigators to apply TO to fluid flow configurations were Borrvall and Petersson (2003) who sought to minimize the pressure drop experienced by the fluid in low Reynolds number flows. A study by Yoon (2010) involved the application of TO to forced convection problems accounting for nonuniform convective conditions that were determined by solution of an advection-diffusion model. Similar studies which incorporate the effects of fluid flow by solving the advection-diffusion equations include Koga et al. (2013), Haertel and Nellis (2017), Qian and Dede (2016), and Subramaniam et al. (2019). Topology optimization for natural convection problems has been recently studied by Alexandersen et al. (2014) where the authors used density-based TO to design complex 2D heat sinks as well as micropumps driven by natural convection. In a later work, the same method was used to design a 3D heat sink for light emitting diodes (Alexandersen et al.,2018). In a similar study by Coffin and Maute (2016), a level-set TO method was introduced to determine optimal solid shapes experiencing steady-state and transient natural convection. Joo et al. (2017) used a simplified surrogate model for natural convection which incorporated a shape-dependent convection coefficient in order to design a 3D natural convection heat sink using density based TO.
Although TO is a powerful design methodology, it can be computationally expensive, especially when considering a complex physical problem such as conjugate heat transfer with conduction in the solid region and natural convection in the fluid. The objective of this work is to propose and demonstrate an overall evolutionary design method (EDM) for 2D solids cooled by natural convection. The EDM consists of a heat transfer (HT) sub-model and one of four proposed shape evolution method (SEM) sub-models. The proposed EDM is comparable to bidirectional evolutionary structural optimization (BESO), in which a solid structure evolves by the simultaneous removal of unnecessary material that is less stressed and addition of material to regions that are more stressed. Through this process the structure gradually approaches the maximum overall stiffness per unit material volume. The difference between BESO and the design method of this work is that the solid is reconfigured based on local temperatures or heat fluxes. Similar evolutionary design methods have been applied to heat conduction problems, such as (Li et al. 1999 and2004). The proposed design method is relatively simple compared to many TO methods because (i) the shape evolution does not require a formal optimization analysis and (ii) the governing equations are solved using a finite volume method, making the solid redistribution straightforward. The predictions of the proposed EDM will be compared to an optimized heat sink design (Coffin and Maute, 2016) operating under identical conditions. A parametric study will also be conducted to observe the influence of solid thermal conductivity and the domain size (Rayleigh number) on the final solid shape.

PHYSICAL AND NUMERICAL MODELS
The situation of interest is depicted in Fig. 1. A two-dimensional, square domain is composed of a fixed amount of solid and an adjacent fluid. As shown, the solid is arbitrarily specified to have an initial, semicircular cross-section of radius R and cross-sectional area At = R 2 /2. A heat rate per unit length, q', is applied at the bottom center of the domain, while the top boundary is isothermal at a reference temperature To. All remaining boundaries are adiabatic. The temperature of the solid at the location where q' is applied is to be minimized by allowing the solid shape to evolve under the constraint that the total amount of solid remains constant, affecting both conduction in the solid and free convection in the fluid. The conduction and convection processes are coupled at the solid-fluid interface.
The overall EDM simulation is initiated with the domain shown in Fig. 1. After the steady-state temperature distributions (including the maximum solid temperature at x = H/2, y = 0) are calculated with the HT sub-model (Section 2.1) the solid cross-sectional shape is modified according to a SEM (Section 2.2), the conjugate conduction-convection problem is re-solved, and the maximum solid temperature is recalculated. The heat transfer prediction -shape evolution -heat transfer prediction process is continually repeated with the goal to reduce the solid temperature at the location where q' is applied.

Heat Transfer Sub-Model
Heat transfer within the fluid is described by (i) the conservation of mass equation, (ii) the Navier-Stokes equations and (iii) the conservation of energy equation. A Boussinesq, Newtonian fluid is considered, and viscous dissipation is neglected. The fluid flow and heat transfer processes are steady-state, and it is assumed that all thermophysical properties are constant. Conditions will be specified so that the fluid flow is laminar. Therefore, the governing equations for the fluid are where u and v are the x and y components of the velocity respectively, p is the pressure, µf is the viscosity, f is the density, f is the coefficient of thermal expansion, cp,f is the specific heat, and kf is the thermal conductivity of the fluid. Heat transfer in the solid is governed by an energy equation similar to Equation (4), but with u = v = 0 and ks specified instead of kf. Radiation heat transfer is neglected, so at the solid-fluid interface the heat flux in the solid normal to the interface is equal to the heat flux in the fluid normal to the interface, and the temperature of the solid is equal to the temperature of the fluid. No-slip conditions are applied at the boundaries of the domain. Local temperatures and heat fluxes are obtained using the finite volume method (Patankar 1980), with each control volume being entirely solid or entirely fluid. A staggered grid was used to solve the discretized forms of the Navier-Stokes equation, the power-law differencing scheme was employed, and harmonic mean thermal conductivities were calculated to properly determine thermal conditions at the control volume surfaces that separate the solid and fluid phases. The equations were solved using the SIMPLE algorithm. Because it is not known a priori how the solid and fluid sub-domains will evolve, the entire computational domain was populated with control volumes of uniform size. The computational model was validated as discussed in the Appendix.

Shape Evolution Methods
A flow chart of the overall EDM is shown in Fig. 2. The simulation begins with specification of an initial solid geometrical shape (Fig. 1).
The governing heat transfer equations are then solved using the HT submodel, and all local heat fluxes and temperatures along the solid-fluid interface are calculated. Within each iteration of the EDM, one control volume is switched from solid to fluid, while a second control volume is concurrently switched from fluid to solid in adherence to a specified SEM.
As discussed in Section 1, the solid reallocation associated with the SEM is similar to that of structural evolutionary methods such as BESO, where solid is removed from locations where it is less beneficial and added to locations where it is more useful. Unlike structural problems, however, evolutionary solid reallocation schemes for conjugate conduction-convection heat transfer problems are less obvious because (i) addition or subtraction of solid along the interface will affect the fluid flow in a nonlinear manner and (ii) the solid is of relatively high thermal conductivity and therefore poses a small thermal resistance relative to that posed by the fluid. Four SEMs are considered here.
In SEM I, all solid (fluid) control volumes along the solid-fluid interface are checked, and the solid (fluid) is converted to fluid (solid) in the control volume that experiences the largest (smallest) heat flux to the fluid (from the solid). In SEM II, all solid (fluid) control volumes along the solid-fluid interface are checked, and solid (fluid) is converted to fluid (solid) in the control volume with the lowest (highest) surface temperature. SEM III is the opposite of SEM II, that is, all solid (fluid) control volumes along the solid-fluid interface are checked, and solid (fluid) is converted to fluid (solid) in the control volume with the highest (lowest) surface temperature. SEM IV is the opposite of SEM I, that is, all solid (fluid) control volumes along the solid-fluid interface are checked, and the solid (fluid) is converted to fluid (solid) in the control volume that experiences the smallest (largest) heat flux to the fluid (from the solid).
The overall EDM simulation is curtailed when either of two criteria is achieved. The first curtailment criterion is reached when part of the solid becomes disconnected from the rest of the solid and a boundary of the domain, resulting in some of the solid "floating" unrealistically within the surrounding fluid. The second curtailment criterion is met when the predicted solid shape oscillates from iteration-to-iteration, marking the end of the solid shape evolution.

EDM PREDICTIONS
Results are presented for the following cases. Prior to employing the full EDM, the HT sub-model (without utilizing a SEM) is used to both (i) replicate the heat transfer and fluid flow predictions and (ii) quantify the overall thermal resistance of the conjugate conduction/convection system, ′ ≡ [ ( 2 ⁄ , 0) − ]/ ′, associated with the benchmark solid shape predicted by Coffin and Maute (2016). Second, using the same conditions specified by Coffin and Maute (2016), final solid shapes and the associated overall thermal resistances are reported for predictions generated by the overall EDM model using each of the four SEM sub-models described in Section 2.2. After these predictions are reported, an additional constraint is imposed in the EDM model of Fig.  2, and results are discussed. Finally, the influence of (i) the solid phase thermal conductivity and (ii) the domain size (Rayleigh number) on the final solid shape and overall thermal resistance is examined.
The optimal shape of the benchmark study is shown in Fig. 3a, along with the predicted benchmark temperature and streamline distributions. As evident, thermal and velocity conditions are not symmetric about x = H/2 despite the symmetric conditions of the problem and consist of a large clockwise circulation of the fluid (air) that sculpts an irregular, nearly-isothermal sloped solid that is thicker on the LHS of the computational domain. Note that no solid is in contact with the left adiabatic wall of the domain; this is attributed to the complexity of the conjugate conduction and convection heat transfer processes as well as inclusion of an additional constraint in the benchmark model. Specifically, the total perimeter of the solid was constrained by Coffin and Maute (2016); a similar constraint was not incorporated in this study due to the use of an orthogonal computational mesh. The consequence of this additional constraint will be discussed further in Section 3.3.
To demonstrate the veracity of the HT sub-model, the solid shape of Fig. 3a was replicated and is shown in Fig. 3b. The HT sub-model was then applied to the combined solid-fluid domain resulting in the predicted streamline (Fig. 3c) and temperature (Fig. 3d) distributions shown, which are in excellent qualitative agreement with those of the benchmark (Fig. 3a). The maximum temperature reported by Coffin and Maute (2016) is T(H/2,0) = 2.16 K while the maximum temperature predicted by the HT sub-model is T(H/2,0) = 2.19 K. The predicted thermal resistance associated with Fig. 3d is ′ = 23.888 K ⋅ m/W.

Shape Evolution using the Four SEM Sub-models
With the accuracy of the HT sub-model demonstrated, the overall EDM model, composed of the HT and SEM sub-models, was used to predict the topology of the solid of total area At = 39.5 mm 2 . Solid shapes at various stages of iteration of the overall EDM model, corresponding to a cumulative amount of displaced solid AD, and number of design iterations, iter, are shown in Fig. 4 for each of the four SEM sub-models. For SEM Method IV, the final shape was reached after only AD = 1 mm 2 of solid was displaced, which corresponds to 4 design iterations, and is shown in the upper right corner of the figure. With use of SEM IV, curtailment of the EDM was triggered by the oscillation between two solid shapes on subsequent iterations resulting in ′ = 28.008 K ⋅ m/W. Shapes corresponding to AD = 1 mm 2 are also shown for the other three SEMs in the top row of the figure.
When AD = 7 mm 2 , the solid shape predicted using SEM III becomes separated from itself, resulting in part of the solid being artificially suspended in the fluid, and the simulation was curtailed. Just prior to solid separation, ′ = 27.375 K ⋅ m/W.
Although the solid shapes predicted with SEM I and SEM II are similar when AD = 7 mm 2 , the shape predicted using SEM II is noticeably more symmetric. Shortly thereafter, at AD = 8 mm 2 , the simulation associated with SEM II was curtailed due to oscillation between solid shapes in subsequent iterations with ′ = 27.041 K ⋅ m/W. In contrast, the solid shape predicted with use of SEM I continues to evolve, and eventually reaches its final configuration only after a significant amount of solid is displaced. (Interestingly, the amount of displaced solid exceeds the amount of solid present, AD = 43 mm 2 > At). The final shape predicted using SEM I bears similarity to that of the benchmark (Fig. 3a), with the sloped solid being thicker on the LHS of the domain. Unlike the benchmark, however, the solid makes contact with the vertical LHS insulated surface of the enclosure. The thermal resistance associated with the final shape is ′ = 26.126 K ⋅ m/W Curtailment of the simulation was triggered by the oscillation between solid shapes on subsequent iterations of the EDM. Hence, the solid that is shown is the final shape generated by the overall EDM model. A magnified view of the final (AD = 43 mm 2 ) predicted solid shape using SEM I, along with the corresponding streamline and isotherm distributions are presented in Fig. 5. Similar to the benchmark results of Fig. 3, the system displays highly asymmetrical behavior with the solid displaced in the direction of the adjacent fluid velocities, resulting in sculpting effects similar to those noted for the benchmark result shown in Fig. 3a. The streamline and temperature distributions of Figs. 3 and 5 are also qualitatively similar. Unlike the solid of the benchmark study (Fig. 3a), however, no constraints were employed to generate the results of Fig. 5 with respect to the solid's peripheral length, so the disparity between the details of final solid shape predicted here and that of the benchmark is not surprising.
The evolution (open symbols) of the predicted thermal resistance associated with each SEM, as well as that of the benchmark geometry (dashed line), is reported in Fig. 6. The filled symbols correspond to the four final geometries of Fig. 4. As is evident for all of the SEM submodels employed, the thermal resistance generally decreases as more solid is displaced. This trend is attributed to the fact that the total area of the solid-fluid interface generally increases as adjustments are made to the solid shape. Compared to SEM II through SEM IV however, SEM I displaced over four times as much solid and has the smallest  final thermal resistance. Although the solid geometry predicted by SEM I ( Fig. 5a) is similar in shape to that of the benchmark (Fig. 3a) its thermal resistance ( ′ = 26.126 K ⋅ m/W) is closer in value to that of the initial geometry (Fig. 1, ′ = 28.062 K ⋅ m/W) than to that of the benchmark ( ′ = 23.888 K ⋅ m/W). This result is unexpected, and is attributed to the buildup of solid adjacent to the insulated left face of the enclosure that is promoted by SEM I. The solid adjacent to the insulated vertical wall is relatively inactive thermally, and might be better utilized if it had been moved from the adiabatic boundary to the upper surface of the solid.

IMPLEMENTATION OF A CUTOFF CONSTRAINT IN SEM I
SEM I operates by moving solid to the solid-fluid interface location that experiences the smallest heat flux. Hence, the solid will accelerate toward an adiabatic boundary, eventually becoming thermally inactive as it makes contact with the boundary. In order to prevent the solid from reaching the insulated LHS of the domain, SEM I was modified by introducing an additional criterion that, as will become evident, prohibits the solid from making contact with the insulated vertical surface. This modification is based on a cutoff heat flux, ′′ , that is quantified in terms of a cutoff ratio, , and the maximum local heat flux along the solid-fluid interface, ′′ . If the local heat flux along the solid-fluid interface is below the cutoff heat flux, solid will not be added to that location. Therefore, the modified version of SEM I is as follows.
(1) Remove solid from the location of ′′ and (2) Add the same amount of solid to the location of the smallest q ′′ that is greater than ′′ where To demonstrate the modified SEM I, a value of RC = 0.05 was specified and the EDM predictions were compared to those using the un-modified version of SEM I (RC = 0). The final solid shape, streamline distribution and temperature distribution are presented in Fig. 7. Clearly, implementation of the cutoff criterion prevents the solid from contacting the adiabatic wall, with the final solid shape (with RC = 0.05) bearing more similarity to the benchmark shape than the shape predicted with the un-modified version of SEM I (RC = 0). Despite being in better qualitative agreement with the benchmark shape, the thermal resistance associated with the final shape and RC = 0.05 is ′ = 26.119 K ⋅ m/W which is approximately the same as for the solid of Fig. 5a. The evolution history of the thermal resistance generated by use of SEM I with and without the cutoff ratio criterion is shown in Fig. 8. The cutoff criterion using RC = 0.05 comes into play only when AD = 24 mm 2 . For AD > 24 mm 2 , the thermal resistance is smaller for RC = 0.05 than for RC = 0, although the differences in the values of ′ predicted with and without the cutoff ratio imposed are small.

Parametric Simulations
As demonstrated in Section 3.2 and 3.3, the EDM model is able to predict solid shapes that are in qualitative agreement with those of the benchmark, resulting in reduced thermal resistances across the computational domain. Because heat transfer across the domain is due to conduction in the solid phase and convection in the fluid, parametric simulations were performed in order to assess the sensitivity of the predicted solid shape (and the corresponding thermal resistance) to changes in the conduction process (i.e. variations in the thermal conductivity of the solid), and changes in the convection process (i.e. the size of the computational domain, or Rayleigh number). conductivity of the solid used in the benchmark study and in the preceding simulations is approximately 5 orders of magnitude greater than that of the fluid. Therefore, the solid exhibits nearly isothermal conditions throughout. Figure 9 shows the final predicted solid shape, thermal resistance and amount of displaced solid using SEM I with RC = 0 for ks = 237, 23.7, and 2.37 W/m•K. Slight differences in the solid shapes corresponding to ks = 237 and 23.7 W/m•K can be noted upon close inspection. Decreasing the solid's thermal conductivity would generally lead to a higher overall thermal resistance, but as evident, the thermal resistance is slightly lower for the ks = 23.7 W/m•K case than for the ks = 237 W/m•K benchmark case. This unexpected result is attributed to offsetting effects associated with the modest differences in the solid (and fluid) domain shapes of Fig. 9a and Fig. 9b. As the value of the solid phase thermal conductivity is further reduced to ks = 2.37 W/m•K (Fig. 9c), the solid shape exhibits a notably different shape relative to the two higher thermal conductivity simulations, and the thermal resistance begins to increase. Additional reductions in ks (not shown) result in even greater values of the overall thermal resistance. Predictions using SEM I and RC = 0.05 are shown in Fig. 10 and exhibit the same trends as noted in Fig. 9. Again, the ks = 23.7 W/m•K case (Fig. 10b) yields the lowest overall thermal resistance.

Influence of the strength of convection (Rayleigh number)
Uniformly increasing the size of the overall domain, while holding all other conditions constant, will strengthen convection heat transfer and is expected to reduce the overall thermal resistance. Figure 11 shows the final solid shape, thermal resistance, and amount of solid displaced using SEM I with RC = 0 for H = 30, 45, and 60 mm. Note that the Rayleigh number associated with each domain size is calculated using the maximum predicted temperature difference, = ∆ 3 ⁄ , where = , ⁄ . As evident in the figure, (i) the overall resistance is decreased by approximately 40 percent as a result of doubling H and (ii) the solid shapes associated with larger Ra exhibit better qualitative resemblance to the optimized shape of Coffin and Maute (2016). The effect of the increased strength of convection on the final solid shape is as expected. Figure 12 illustrates the sensitivity of the final solid shape, thermal resistance, and amount of solid displaced using SEM I and = 0.05 on the domain size. As for the RC = 0 case, the solid shapes associated with higher Rayleigh numbers have lower thermal resistances and are in better qualitative agreement with those of the benchmark prediction.
The evolution history of the overall thermal resistance using SEM I and H = 60 mm is shown in Fig. 13. Unlike the trends evident in Fig.  8 (H = 30 mm), the predictions exhibit a more substantial increase in thermal resistance (for 10 ≲ ≲ 25) followed by a decrease in the thermal resistance to its final value. For 10 ≲ ≲ 25, the solid is being moved toward the LHS vertical adiabatic boundary. As the boundary is approached, the solid becomes less thermally active, increasing the overall thermal resistance. The increase in the thermal resistance is not as pronounced for the lower Rayleigh number cases (Fig. 8) because the solids for those cases retain somewhat semi-circular shapes as they approach the LHS boundary, whereas the solids associated with Fig. 13 transition from semi-circular shapes to flatter shapes in the range 10 ≲ ≲ 25. The flattening solids become further removed from the cold top boundary, increasing the overall thermal resistances to a more significant degree than for the lower Ra cases. For either case, for ≳ 30 the solids become less flat and evolve upward as shown in Figs. 11a and 11c, increasing the exposed solid surface area and reducing the thermal resistance.
The percentage reductions in thermal resistance for the three different domain sizes, using SEM I with and without implementation of the cutoff ratio, are presented in Table 1. The largest reduction in    thermal resistance is associated with the lowest Rayleigh number. Implementation of the cutoff ratio has a greater effect for the larger Rayleigh number cases

SUMMARY, CONCLUSIONS, AND RECOMMENDATIONS
A novel EDM model, consisting of heat transfer (HT) and shape evolution method (SEM) sub-models was developed and its use was demonstrated. Four SEMs were considered, and predictions were compared to a benchmark optimized design. The best performance was achieved using SEM I in which high thermal conductivity material is removed from the location of the largest interfacial heat flux and added to the location of the smallest flux. An additional constraint was proposed and examined, with the objective of prohibiting the solid from contacting the adiabatic side boundary of the domain. This constraint lead to a thermal resistance similar in value to that predicted without the constraint. Parametric simulations were conducted to assess the influence of conduction (solid thermal conductivity) and convection (Rayleigh number). Perhaps unexpectedly, a solid thermal conductivity value (23.7 W/m•K) lower than that used in the benchmark led to a final solid shape with a lower overall thermal resistance. Increasing the strength of convection decreased the thermal resistance and resulted in final solid shapes that bear a closer resemblance to the benchmark study. The decrease in ′ and qualitative similarity of the solid shape generated by SEM I relative to a benchmark solid shape suggest there is merit to using SEM I. However, the predicted thermal resistance is greater than that of the optimized design. Hence, although the EDM proposed here is straightforward and easily implemented, it lacks the accuracy of more complex optimization routines.
General recommendations for future research include: (i) improving the relatively straightforward evolutionary methods proposed here to achieve the high accuracy of TO methods or (ii) developing new computational techniques to reduce the expense of high accuracy TO methods. More specifically, improving the evolutionary methods proposed here might entail: (i) specification of different material redistribution rules, (ii) incorporation of additional or new constraints to guide the evolution of the solid shape, and (iii) developing new or additional curtailment criteria to allow the solid to more fully evolve.

APPENDIX
The HT and SEM sub-models were validated to the extent possible as follows. First, predictions of the convection heat transfer aspects of the HT model were compared to the classic benchmark solution provided by de Vahl Davis (1983) involving a square cavity containing air with adiabatic top and bottom surfaces and isothermal side walls. Values of the average Nusselt number were predicted for Rayleigh numbers ranging from 10 3 to 10 6 . Using a 60 × 60 uniform mesh, the largest difference between predicted average Nusselt numbers was 3.1 percent at Ra = 10 6 . The difference between the maximum local Nusselt numbers was 7.1 percent at Ra = 10 6 . The HT sub-model was then used to predict a benchmark solution (Costa, 2012) for natural convection of air in a partitioned cavity involving conjugate conduction-free convection heat transfer, as shown in Fig. 14. Two cases were considered: (i) hotter wall on the right and (ii) hotter wall on the left. Both a 50 × 50 mesh and a 100 × 100 mesh were employed in the comparison exercise. The average Nusselt numbers for each case are shown in Table 2. The predictions generated by the HT sub-model are in good agreement with the benchmark solutions, and there is only a slight improvement in agreement with the benchmark as the mesh is refined from 50 × 50 to 100 × 100. Based on the preceding discussion, a 60 × 60 mesh was Fig. 14 Domain for the conjugate conduction-free convection benchmark solutions (Costa, 2012). A third comparison was made to predict the evolution of a geometrical shape in a conduction scenario similar to that considered by Li et al. (1999). In this comparison, an initially square domain is (i) completely filled with a conducting material and a second isothermal material at its square center and is (ii) exposed to isothermal exterior boundaries. The solid shape (defined by the interface between white and black areas in Fig. 15) evolves by continually removing control volumes adjacent to the boundaries of the solid that experience the smallest average heat flux. The boundary of the solid shape is maintained at the same temperature as the original boundary of the square domain. As evident in Fig. 15, the qualitative agreement between the two predicted solid shapes is excellent. The iteration history, showing evolution of the minimum and maximum local heat fluxes along the solid boundary is shown for both studies in Fig. 16. As in Fig. 15, the predictions and benchmark results are in good qualitative agreement, especially at later stages of the iteration when the solid shapes approach their final configurations. Differences between the current and benchmark results are attributed to the different numerical techniques used (finite volume versus finite element of Li et al., 1999).  A final verification exercise was conducted to determine the sensitivity of the final predicted solid shape to the initial solid shape specified in the overall EDM model of Fig. 2. Specifically, predictions based on the semicircular initial solid considered in this study were compared to those generated by specifying an initially nearly rectangular solid containing the same amount of material as the semicircle. Figure 17 includes the initial shape (top row) and final shape without (RC = 0, middle row) and with (RC = 0.05, bottom row) the cutoff criterion applied. As evident, the final shapes associated with RC = 0 are nearly independent of the initial specified shape. For RC = 0.05, the final shapes are similar but exhibit more noticeable yet minor variations; the initially rectangular shape yielded a slightly lower (2.5 percent difference) thermal resistance than the initially semicircle shape simulation as reported in Table 3.

Fig. 17
Influence of the initial solid shape on the final solid shapes without (RC = 0) and with (RC = 0.05) the cutoff criterion.