A tessellated continuum approach to thermal analysis: discontinuity networks

Tessellated continuum mechanics is an approach for the representation of thermo-mechanical behaviour of porous media on tessellated continua. It involves the application of iteration function schemes using affine contraction and expansion maps, respectively, for the creation of porous fractal materials and associated tessellated continua. Highly complex geometries can be produced using a modest number of contraction mappings. The associated tessellations form the mesh in a numerical procedure. This paper tests the hypothesis that thermal analysis of porous structures can be achieved using a discontinuous Galerkin finite element method on a tessellation. Discontinuous behaviour is identified at a discontinuity network in a tessellation; its use is shown to provide a good representation of the physics relating to cellular heat exchanger designs. Results for different cellular designs (with corresponding tessellations) are contrasted against those obtained from direct analysis and very high accuracy is observed.


Introduction
In the recently developed method of tessellated continuum mechanics, physics defined on fractal porous media is mapped to a corresponding continuum [1,2]. The classical Galerkin finite element method (GFEM) can then be combined with tessellated continuum mechanics to describe thermal physics in cellular structures. The drawback with the Galerkin approach, however, is that discontinuous behaviour arising from the closing of the pores during the mapping (which is pervasive on tessellations) cannot readily be accommodated. This paper extends the work in references [1,2] by incorporating a network of lines/surfaces on a tessellation at which discontinuous behaviour occurs. The discontinuities can be viewed as enforcing internal boundary conditions where (say) heat fluxes on different parts of a pore's surface are mapped to a single line. The new feature is called a discontinuity network and is designed to capture more completely the physics on cellular materials.
Many natural materials possess hierarchical microstructures extending over length scales that cover many orders of magnitude [3]. Those microstructures observed with geometric similarity over a range of scales are denoted as self-similar [4]. Self-similar objects possessing interesting mathematical properties have received the attention of researchers in many fields including mathematics, biology, chemistry and engineering. Analysis of these types of structures falls under the heading of multi-scale modelling [5], where physical self-similarity can in principle be used to describe scaling relationships connecting behaviours across a range of length scales.
Of particular interest in this paper is the behaviour of self-similar porous structures subject to fluid-induced thermal loading.
The performance of heat exchangers involving thermo-fluids and cellular materials is an area of practical interest, designs of which can be found in everyday use for the cooling of printed circuit boards [6]. Self-similar porous materials that maximise fluid-solid contact areas offer potentially very high rates of heat transfer. This, when combined with high thermal conductivities and good enhancement of fluid mixing, may be exploited for heat exchanger designs [7]. Examples of commercial heat exchangers which incorporate a complex geometrical structure include miniature shell and tube heat designs [8]. Alternative configurations are honeycomb metallic structures [9], metallic foams [10] and plate-fin designs [8]. These designs play an important part in applications requiring high-density energy transfers, and growth in their use and application is expected [8]. Experimental evidence showcasing the benefits of cellular designs can be found in reference [11].
The complex geometries involved in porous/cellular heat exchangers raise particular challenges for analysis. These challenges include: the accurate representation of the geometry; the application of boundary conditions; meshing problems that arise with large and thin sections; and any computational constraints imposed by the large number of elements involved. The standard approach to the analysis of flow through porous heat exchangers would typically not include a direct description of the geometry of the porous media (see  Fig. 1 Pre-fractal and corresponding tessellated domains containing discontinuity networks for the Sierpinski gasket references [12][13][14]). Classical approaches founded on continuum equations, and involving coarse-grained geometrical parameters such as porosity and permeability [15] provide limited accuracy. Recent alternatives to these classical treatments include discrete techniques such as the lattice Boltzmann method (LBM), which is particularly attractive for low Reynolds number flows involving complex boundary conditions [16]. An example of pore-scale turbulent flow modelled using LBM is provided in [11], where four types of porous structures are considered. Unfortunately, despite the LBM offering some advantages, it still suffers in a similar manner to continuum approaches as a consequence of the complex geometries involved [17].
In order to represent better the types of complex geometries found in cellular material, the authors have investigated the use of fractals [2,18], which arise in many other areas of science [19] and play a part in describing a variety of behaviours in nonlinear systems [20]. Fractals can be constructed in a number of ways, but a particularly attractive procedure is by application of a small number of affine contraction maps [2,20,21]. Although complex geometries can readily be represented using fractals, their subsequent use in practical analysis is beset with difficulties. The principal concern is the lack of meaningful measures on the fractal, where traditional notions of length, area and volume are undefined as a consequence of non-integer Hausdorff fractal dimensions [22]. Consequently, traditional physical quantities such as fluxes (rates per unit area) and densities (per unit volume) are likewise ill-defined. Attempts to overcome these difficulties can be found in the literature [23,24]. One possible method to avoid the problem altogether is to use pre-fractals, which are Pre-fractal and corresponding tessellated domains containing discontinuity networks for the Sierpinski carpet created with fractal construction procedures, but do not suffer from non-integer Hausdorff fractal dimensions. A popular indirect approach is through the use of fractional derivatives [25,26] in the governing equations describing the physics. An example is the fractional model considered in reference [27] for describing the flow of two incompatible liquids through homogenous porous media. The drawback of approaches of this type is that the governing equations are not directly derivable from the underlying physics. An alternative approach is to consider transport equations in their integral (Euler) form. A variant of this approach is considered by Tarasov [28,29] and underpins the fractional derivative approach outlined by Ostoja-Starzewski [30,31]. In this paper, the physics on the cellular material-whose geometry is represented by a pre-fractal-is linked to the physics on a tessellation by means of a hole-fill map. The map is identified as hole-filling because when applied to a pre-fractal it has the effect of closing the holes within the porous pre-fractal structure. Holefill maps can be constructed by means of function composition or by identifying corresponding elements in the pre-fractal and tessellation. Function composition is more suitable for analysis, so focus is directed towards the latter, which is preferred for the type of numerical work considered here.
A novel aspect of the work is the creation of a discontinuity network, which is a network of lines/surfaces where discontinuous behaviour is possible. The discontinuity network is simultaneously created during the formation of the fractal and associated tessellation, resulting in an elegant procedure controlled through a Pre-fractal and corresponding tessellated domains containing discontinuity networks for a finger-like porous fractal relatively small number of maps. The physics on the tessellations/pre-fractals is represented using weighted transport equations, which can be formulated directly without recourse to partial differential equations (PDEs). This aspect is important particularly if fractals (rather than pre-fractals) are involved, as the governing PDE in this latter case may not be obvious; indeed, it may not even be definable using traditional derivatives. The paper as a whole investigates how the finite element method can be applied to a tessellated mesh (where tiles double up as elements) incorporating a discontinuity network. A critical aspect of this work is that the physics on the porous structure are represented exactly in the continuum. This facilitates very precise analysis and is particularly important when high accuracy is required. The procedures for construction of pre-fractals, tessellations and discontinuity networks for three classical fractals are presented in Sect. 2. The transport theory linking the tessellated space with the physical space is described in Sect. 3, with particular focus on heat transfer. This is followed by weighted transport equations in Sect. 4, which can be formulated directly and lead immediately to the finite element method. The role a discontinuity network plays in the finite element method is discussed. Transient and steady-state analytical solutions are given in Sect. 5 on a simple dust fractal and contrasted with numerical predictions. The role and importance of discontinuity networks are investigated in Sect. 6. In the absence of analytical results, numerical predictions-both with and without a discontinuity network-are compared with results obtained for a pre-fractal via the commercial package ABAQUS. The influence of different tessellation refinements is investigated in Sect. 7, where an efficient procedure for tessellation refinement is examined. Section 8 examines the influence of alternative expansion maps and confirms that they do not influence the accuracy of the numerical predictions.
tessellationT k is exactly the same as the number of pre-fractal elements inÊ k ; there is thus a one-to-one correspondence between tiles and pre-fractal elements, achieved by using an identical number of expansion and contraction maps. The one-to-one correspondence means that analysis results obtained on a tessellation T k can be immediately lifted and returned to the corresponding pre-fractalÊ k .
Typical tessellations for some classical fractals are depicted in Figs. 1, 2 and 3 with affine contraction and tessellation maps provided in Appendices A-C. An aspect deduced on examination of Figs. 1, 2 and 3 is that E 0 (not shown in figures) is supplied with an initial tessellation. The initial domainÊ 0 has six tiles for the pre-fractals in Fig. 1 and eight for those in Figs. 2 and 3, and each tile is purposely triangular. The initial tiling has two distinct functions: first it serves to form a piecewise linear hole-fill map and second it provides a means to modify the number of tiles inT k for the purpose of numerical analysis.
The creation of pre-fractals and tessellations is achieved independently using bespoke iteration function schemes involving the (non-unique) maps given in "Appendices A-C". A hole-fill map-on application to a pre-fractal-has the effect of closing holes and producing the corresponding tessellation. Although it is possible to form a hole-fill map explicitly through composition of contraction and expansion maps (and their inverses), there is little merit in so doing. The independent creation of pre-fractals and tiles means that triangular elements in each domain are related. This relationship is shown explicitly in Figs. 1, 2 and 3 by means of arrows connecting a selection of tiles in each domain. A hole-fill map for a tile is simply a linear relationship which maps points from a pre-fractal element to a tile; the map is linear as a consequence of the triangular domains involved.
Figures 1, 2 and 3 also highlight a discontinuity networkD k in each tessellationT k . It is important to appreciate that likeT k andÊ k a discontinuity networkD k is created recursively and satisfies the relationship Thus a network is formed by forming a union withD 1 and images of the previous network under the maps P i . Examination of Figs. 1, 2 and 3 reveals that networks are associated with hole edges in pre-fractals. The closing of holes has the effect of bringing together the edges which appear on the discontinuity network in a tessellation. This makes clear why discontinuous behaviour is ubiquitous on any tessellation-there is little expectation that temperature and heat transfer rates are identical on opposite edges of a hole. The networkD 1 inT 1 , in Fig. 1 for example, is obtained from the boundary ofT 0 =Ê 0 , and the networkD 2 inT 2 includesD 1 and the three expansion maps {P 1 , P 2 , P 3 } (see "Appendix A") applied toD 1 Similarly, the networkD 3 inT 3 again includesD 1 and the three expansion maps {P 1 , P 2 , P 3 } applied to the networkD 2 inT 2 (i.e.D 3 =D 1 ∪ P 1 D 2 ∪ P 2 D 2 ∪ P 3 D 2 ). It is thus apparent that creation of discontinuity networks parallels the recursive procedure for creation of a tessellation.

Transport theory
The ability to match physics on pre-fractals and tessellation was established in reference [2]. The basic idea is to describe the physics of interest in each domain on control volumes, which identify continuous spaces through which matter can flow. Control volumes are employed because they are concerned with the description of physics on a continuum (a continuous space) rather than on a discontinuous set. Thus, the difficulty posed by fractal geometry can in principle be overcome with physics described using traditional measures. Underpinning the approach adopted here is the continuum hypothesis (advanced by Cantor in 1878 [32]), which asserts that any infinite set (in 1-D) forms a bijection with either the set of natural numbers or the set of real numbers. The fact that the Cantor dust (for example) forms a bijection with the real line was proved by Cantor in 1883 [33]. The concept underpins the strategy adopted in this paper where pre-fractals (and in the limit, fractals) are mapped to continuous tessellations. Consider two control volumes s and r , with s located in the physical domain in which the cellular material resides (represented byÊ k ); r is located in the tessellated space. For clarity, attention is restricted to heat transfer problems, although more complex physics is possible via the inclusion of additional transport equations. The transport equations for transient heat transfer through stationary control volumes (the latter again chosen for simplicity) are given by d dt where h is the specific enthalpy,Q is a heat source and the heat flux is defined asq · n =ĥ (T − T ∞ ) withĥ a heat transfer coefficient and T ∞ is the bulk temperature of an external cooling/heating medium. The differentials dV s and dV r are measures of volume on the two control volumes and are thus well defined. The tessellation is continuous and dV r = dV tess , where dV tess is a volume measure on the tessellated set. Defining dV frac as a measure of volume on a pre-fractal, dV frac = μdV s , where μ is a support function taking the value of zero or one. It follows from the discussion in Sect. 2, however, that hole-fill maps (being a material map) relate dV frac to dV tess with a relationship of the form dV tess = |F| dV frac . This expression happens to be valid for precisely the situation where dV r = dV tess and dV s = dV frac , and it follows that dV r = |F| dV s . The interpretation placed on Eq. (1) therefore is that the integrands involved may not be defined for points off the pre-fractal, which is in fact an interpretation typical to classic control volume theory. In view of the continuous tessellation and the validity of the relationship dV r = |F| dV s almost everywhere on the tessellation, any support function is removable, as it takes the value of one. Substitution of Nanson's identities [34] which can be compared with Eq. (1).   Recognising that a mass-conserving map satisfies ρ r dV r = ρ s dV s , it follows that matching physics between the tessellated continuum and physical spaces is achieved with ρ r = |F| −1 ρ s , h r = h s ,Q r =Q s anḋ q r = |F| −1 F ·q s . The convective heat transfer boundary condition can be provided in the form oḟ which is satisfied for T r = T s andĥ r d r =ĥ s d s giving ∇ s T s = F T · ∇ r T r , and (in view of the identitẏ q r = |F| −1 F ·q s ) leads toq r = −K r · ∇ r T r . The relationship of thermal conductivity between tessellation and pre-fractal is where K r is the associated orthotropic conductivity tensor on the tessellated domain and K s is the isotropic conductivity on the pre-fractal.

The finite element with discontinuity networks
In the absence of a governing partial differential equation, the finite element method can be applied directly to a transport equation via the inclusion of a weighting function [35]; such an approach applied to Eq. (1) gives d dt where W s is a weighting function which is assumed continuous at least up to the first derivative.
Substitution into Eq. (7) of the relationships established in Sect.
which matches Eq. (6) on setting W r = W s . The finite element method describing the physics on a tessellationT k (and henceÊ k ) is obtained by setting where N i is a shape function and e r is a tile inT k ; the tiles by choice double up as elements in the proposed numerical procedure.
Three-noded triangular elements are used in the analysis of a tessellation, which reflects the restriction imposed by the retention of linear hole-fill maps. Those elements that share no edges with a discontinuity network are continuous in both W r and T r and are consistent with a standard Galerkin finite element method. For elements connected by an edge on a discontinuity network, continuity cannot be assumed. Each element in this case is subject to a mapped convective boundary condition from the associated hole in the corresponding pre-fractal in accordance with Eq. (4).   10 Temperature plots along x = y for pre-fractalsÊ k using UoMFEC and ABAQUS-direct Table 6 Maximum temperature difference atD k for the Sierpinski carpet The commercial finite element package ABAQUS (version 6.13) has also been used for the purpose of confirming the validity of the theory. The methods proposed in the preceding sections do not directly analyse the pre-fractals but rather their tessellations. The ABAQUS code is able to analyse directly the pre-fractal and thus provide benchmark solutions. The elements selected for the analysis with ABAQUS are also linear three-noded elements.
It should be appreciated that the proposed work is aimed at very precise analysis on a continuum and takes advantage of the manner in which tessellations and pre-fractals are generated. The process is efficient in this regard involving recursive processes; tessellations, pre-fractals and discontinuity networks are all produced recursively. It is recognised, however, that the use of a tessellation, which is recursively produced, as a mesh, does place a constraint on the possible distribution of elements forming the mesh. This can possibly result in too few elements being located in regions of greatest need or equally too many in regions where there is no need.
Prior to testing out the full numerical theory on realistic fractals, however (as depicted in Figs. 1, 2, 3), it is of interest to test out the ideas initially on a relatively simple fractal on which analytical solutions can be developed.

Exact solutions on Cantor dusts
This section is concerned with the development of the analytical solution on the Cantor dust in one dimension, with internal heat source and convection boundary conditions. The exact solution is compared with the finite element method developed in the earlier sections. The governing partial differential equation arising from the heat transport Eq. (2) is of the form whereQ 0 is a heat source loading density on tiles in any tessellation.

Pre-fractal and tessellation construction
The pre-fractalsÊ k for the Cantor dust and corresponding tessellationsT k with highlighted discontinuity networks (appearing as discontinuous points) are shown in Fig. 4. The contraction maps that produceÊ k in Fig. 4 are along with the expansion maps forT k are with the original setÊ 0 = [0, 0 ] used for both maps. The expansion maps are not unique although those shown in Eq. (12) are a natural choice, producing tiles of equal size. When viewed as a mesh it is evident the maps in Eq. (12) provide the coarsest mesh possible originating from only one element onÊ 0 . Mesh refinement is possible and one approach is to supply an initial mesh toÊ 0 , where, for example, settingÊ 0 = [0, 0 /2] [ 0 /2, 0 ] provides the two-element refinement in Fig. 5a, b. Equivalently, composition of the expansion maps can be used where, for example, which, when appliedÊ 0 produces four elements onT 1 in Fig. 5c, d).
Examples of refined tessellations and corresponding pre-fractals are shown in Fig. 5. Observe also the discontinuity networkD k highlighted at particular points in the tessellations, which correspond to pre-fractal holes. The correspondence of elements means it is relatively straightforward to stipulate an hole-fill map betweenÊ k andT k ; where and it is assumed here that nodes are enumerated from left to right, so that with k being the length of an element inÊ k [1]. For this relatively simple case, the hole-fill map can be written in terms of a support function μ k as [1] (i) (ii) (iii) Fig. 13 Contour temperature plots forÊ 3 using ABAQUS-direct. i Mesh onÊ 3 with eight elements per pre-fractal element. ii Temperature contours onÊ 3 with eight elements per pre-fractal element. iii Temperature contours onÊ 3 with 152 elements per pre-fractal element where D 1 is the Hausdorff fractal dimension (see Ref. [36] for derivation).
The relationship between material on pre-fractals and tessellations is readily obtained with knowledge of the hole-fill map (Eq. (15)). Conservation of mass onÊ k andT k gives ρ r = |F| −1 ρ s = k Note that in view of the fact that k is constant for each element of the kth pre-fractal, the density ρ r and thermal conductivity K r are homogenous for the corresponding tessellation.

Analytical solution for a 1-D Cantor dust model
In view of the ultimate objective to describe the behaviour of cellular heat exchangers, it is fitting to assume that coolant flows in a normal direction through the voids in the pre-fractals for the Cantor dust depicted in Table 10 Maximum temperature difference atD k for the finger fractal Figs. 4 and 5. A uniform heat source is applied to each pre-fractal element, defined to ensure that the total rate of energy supplied is constant onÊ k (and consequently onT k ). The loading appears as a constant value ofQ 0 in Eq. (10) as a consequence of the total volume remaining constant on a tessellation. At nodes in the discontinuity network, cooling is set as a boundary condition for each adjoining element and is of the forṁ q r =ĥ x (T r − T wat ), where T wat is the bulk coolant temperature andĥ x is a heat transfer coefficient associated with a pre-fractal hole. In this case, Eq. (4) givesq r =q s , i.e. the cooling experienced at a pre-fractal hole is also experienced at the discontinuity network. External surfaces are exposed to a heat flux of the forṁ q 0 =ĥ 0 (T r − T ∞ ), where T ∞ is the ambient temperature of the surrounding medium andĥ 0 is the associated heat transfer coefficient. To arrive at an analytical solution on the tessellated structure requires the union of solutions obtained on individual tiles in the tessellation.

Steady-state solution for a 1-D Cantor dust
According to Eq. (10), the governing equation onT k is where w is the width of the bar and is part of an additional term to capture convective heat transfer from the top and bottom faces of the bar.
The form the solution takes on a tile can be assumed to be T r = T cf + T source , where the complementary function is of the form T cf (x) = a cosh(αx) + b sinh(αx) and satisfies K r T cf − 2ĥ 0 w −1 T cf = 0 with α = 2ĥ 0 /wK r 1 2 , and where a and b are integration constants. The temperature T source takes the form T source = T ∞ + 2 −1ĥ−1 0 wQ 0 (on assumptionĥ 0 is nonzero) and satisfies K r T source +Q 0 − 2ĥ 0 w −1 (T source − T ∞ ) = 0, which is readily verified on substitution. The solution is where a i and b i are determined by two boundary conditionsĥ i T x=x i i+1 , which applies at the edges of a tile.

Steady-state validation tests
The analytical solution is compared with results from ABAQUS and UoMFEC. The material properties for the Cantor dust material (copper) are: thermal conductivity K s = 400 W/mK, density ρ = 8930 kg/m 3 and specific heat capacity c p = 385 J/kgK. The dimensions for the bar are taken to be: edge length l 0 = 1 m and width w = 1 m. The ambient temperature of the surrounding medium is taken to be T ∞ = 323 K, with the heat transfer coefficientĥ 0 = 200 W/m 2 K. A uniform heat-loading rateQ 0 = 400 W/m 3 is applied to the bar as an internal heat source to account for thermal loading. The required heat transfer coefficientsĥ wat s for the pre-fractal holes are determined using an empirical relationship obtained on combination of the Dittus-Boelter and the Darcy-Weisbach equations [37], which provideŝ   where d hole s = 4Â/P is the hydraulic diameter,Â is cross-sectional area,P is the wetted perimeter (in this case the perimeter of the cooling channel), k is the thermal conductivity of the bulk fluid, μ is the viscosity, ρ is the density and c p is the specific heat capacity at constant pressure.
For water flowing through the hole at temperature T wat = 293 K, with pipe length L = 1 m and pressure drop p = 50 kPa, the heat transfer coefficientĥ wat s with different sizes of d hole s on pre-fractals for the dust is calculated using Eq. (18) and given in Table 1.
Temperature plots on pre-fractalsÊ k for the analytical solution, for solutions from UoMFEC with different initial meshes onÊ 0 and solutions from the commercial package ABAQUS are depicted in Fig. 6, with different values of k. The pre-fractal temperatures from the analytical solution and from UoMFEC are mapped from temperatures obtained on the corresponding tessellated structureT k . To obtain similar results from ABAQUS, two options are available for the Cantor dust: (i) indirect temperature determination via a tessellation with results lifted to the corresponding pre-fractal and (i.e. a procedure similar to that used by UoMFEC and termed ABAQUS-lifted) and (ii) temperatures obtained directly on the pre-fractal (termed ABAQUS-direct). Although not an issue for the Cantor dust (but certainly so for more complex fractals), the use of ABAQUS for the two approaches has particular drawbacks: with ABAQUS-lifted, complex orthotropic material properties for the tessellation have to be determined and input as a separate exercise; with ABAQUS-direct, complex geometries have to be loaded into the software.
The various approaches are tested on the pre-fractals and tessellation depicted in Fig. 5, where temperatures lifted from tessellations are of principal interest. The results are provided in Fig. 6 and Table 2 which, in general, establish the credibility of the tessellation strategy via the high accuracy returned. The averagesD andD % given in Table 2 are obtained with the relationships where superscripts (a) and (b) distinguish the methods (analytical solution, UoMFEC, ABAQUS-lifted and ABAQUS-direct) and n is the number of the data points. Examination of Table 2 reveals relatively small errors, with the largest of these for ABAQUS-direct (i.e. for ABAQUS applied directly to the pre-fractals); this is a consequence of a modelling error arising from the approximate enforcement of continuity on a disconnected set.

Transient solution for a 1-D Cantor dust
The transient equation on the tessellated constructionT k with a discontinuity network for each tile is given by

coefficients. This assumption is an approximation but allows for an analytical solution of the form T r (x, t) = -T (x, t)+-S (x), where -S (x) is the steady-state solution of Eq. (17) and -T (x, t) satisfies
which, on assumption that -T (x, t) = X (x) τ (t) presents as This equation provides solutions with the transient part where β = K r /ρ r c r and α = 2ĥ 0 /wK r 1 2 .
The tile considered here is x i , x i+1 , so it is convenient to set the spatial part of the solution to To satisfy the specific temperatures as end conditions for each element, it is necessary for A i = 0 and ω n = 2 k nπ/ 0 , with B i unspecified, where use is made of the tile length l k = 0 /2 k . It follows that solution on the interval where in the limit t → ∞ the temperature T i r (x, t) → -S(x), i.e. steady-state conditions are obtained after a sufficient period of time.
For convenience, the initial condition along the bar is set to zero, i.e. T (x, 0) = 0, which means on which is a Fourier sine series representation of −-S(x) on the interval x i , x i+1 , where B i n are obtained from the Euler formulae The complete solution is obtained on joining solutions on tiles

Transient validation tests
Transient results from the analytical solution are compared with results from UoMFEC and ABAQUS. The assumption adopted in Sect. 5 Fig. 7, with average differences provided in Table 3; averages are again obtained using Eqs. (18) and (19). Errors in results from UoMFEC and ABAQUS on the tessellationsT k and pre-fractalsÊ k are larger with lower values of k but reduce significantly as k increases. This trend is apparent on examination of Fig. 7 and Table 3. Increasing the value of k has the effect of increasing the number of elements involved, which appears to override deleterious effects coming from an increase in geometric complexity.

The importance of discontinuity networks
The results provided in reference [2] were performed without the use of a discontinuity network, which meant that temperature and heat transfer discontinuities could not be accommodated. It is of interest therefore to investigate the influence of discontinuity networks on the accuracy of predictions on realistic pre-fractals. In the absence of a network, cooling can be achieved along the edge of tiles by application of an appropriate heat sink as described in reference [2].
For the results presented in this section, the material selected for the heat exchangers is copper as used in Sect. 4 and the types of exchanger designs considered are limited to prismatic designs with cross sections represented by pre-fractals. The heat transfer coefficientĥ air associated with flowing air with temperature T air = 323 K around the pre-fractals is taken to be 100 W/m 2 K. Each pre-fractal has water flowing through the voids in order to match typical working conditions pertaining to heat exchangers and is subjected to a uniform heat loading with internal heat sourceQ 0 = 500 kW/m 3 .

Thermal analysis of a Sierpinski carpet heat exchanger
The Sierpinski carpet is selected to represent a closed-pore cellular material for the heat exchanger. The contraction and expansion maps for the pre-fractals and corresponding tessellations are provided in "Appendix B". It is assumed that the water coolant flows through the voids in the pre-fractal depicted in Fig. 2 with temperature T wat = 293 K. The heat transfer coefficientsĥ wat s for each hole in the pre-fractals (with different hole sizes and hydraulic diameters d hole s ) are given in Table 4. Note that heat transfer coefficients in Table 4 are transferred to the corresponding edges in a tessellation but first must be scaled to satisfy Eq. (6).
Two different meshes are defined in ABAQUS as follows: (i) one is a coarse mesh with exactly the same number of tiles as used in UoMFEC and (ii) the other is a sufficiently fine mesh to provide a converged result. Temperature distributions obtained with UoMFEC on tessellationT 3 with and without discontinuity networks, and with results subsequently mapped to pre-fractalÊ 3 are depicted in Fig. 8. Similarly, results obtained directly onÊ 3 are obtained using ABAQUS (on two meshes) which are presented in Fig. 9. To provide a quantitate basis of comparison, the average errors (using the converged ABAQUS solution as a benchmark) are calculated using Eqs. (18) and (19) along the diagonal line x = y onÊ k . These errors are given in Table 5   and corresponding temperature distributions are presented in Fig. 10. It can be seen that errors increase with an increase in k for the tessellated approach without discontinuity networks as a consequence of the unrealistic enforcement of continuity at cooling channel holes. Although in this case reasonable accuracy is obtained in the absence of a discontinuity network, higher accuracy is achieved with its inclusion. Continuity assumptions are enforced on tessellations without discontinuity networks, which is equivalent to temperatures and heat transfer rates on opposite sides of pre-fractal holes having identical values. This lack of realism gives rise to problem-dependent modelling errors. It is of interest, however, to quantify the maximum temperature difference T diff at a discontinuity networkD k to provide a direct indication of the significance of this aspect; T diff is formally defined to be where T + − T − is the jump in temperature across the discontinuity networkD k . In the case of the Sierpinski carpet heat exchanger, the values for T diff are presented in Table 6. The magnitude is significant onD 3 andD 2 but zero onD 1 as a consequence of symmetric boundary conditions and the one symmetrically located hole.
Transient temperature plots mapped onÊ k for the different points shown in Fig. 2 are provided in Fig. 11 for UoMFEC on two tessellations, which are compared with convergent ABAQUS-direct results (ABAQUS applied directly to the pre-fractals with fine mesh). In addition, average differences calculated by Eqs. (19) and (20) are given in Table 7 at two points inÊ k . Similar to the steady-state solution, the discontinuity network has little impact on transient results onÊ 1 . However, the influence of the network becomes obvious with increase in k, where its inclusion provides greater accuracy. Errors tend to be greatest at or close to boundaries of holes, and an increase in k provides a greater number of holes.

Finger-like fractal
It is of interest to investigate the open-pore structure depicted in Fig. 3 to represent a cellular heat exchanger. The open-pore structure is assumed to be placed in a close-fitting container to enforce the flow of coolant through the structure. The contraction and expansion maps for the pre-fractals and corresponding tessellations shown in Fig. 3 are provided in "Appendix C". The heat transfer coefficientsĥ wat s associated with each hydraulic diameter d hole s in the pre-fractals are given in Table 8. Examination of Fig. 3 and Table 8 highlights a peculiar feature with higher heat transfer coefficients appearing at smaller holes on different tessellations. This is an unfortunate consequence of the rather irregular wetted perimeter used in the calculation for the hydraulic diameter.
Temperature distributions for the finger-like porous fractal are provided Figs. 12 and 13. Equations (19) and (20) are again used to determine the average errors (contrasted against the benchmark converged ABAQUSdirect results) along the diagonal line x = y onÊ k . The temperature differences are presented in Table 9 with 24 Principal thermal conductivity ratio (K r /K s ) distribution for Vicsek fractal with hole-fill map formed with expansion map (a). i Contours onT 1 . ii Contours onT 2 . iii Contours onT 3 . iv Graph onT 1 . v Graph onT 2 . vi Graph onT 3 corresponding temperature plots in Fig. 14. It is evident on examination in Table 9 and Fig. 14 that greater accuracy is achieved with the inclusion of a discontinuity network. Errors are observed to be particularly large for the continuous Galerkin approach at extreme corners and remain high despite higher values of k. Large error at the corners is an indication that insufficient elements are being placed there, which is a constraint imposed by tiles doubling up as elements with tessellations being created recursively. Equation (29) is again employed to quantify the extent of the temperature differences onD k for the finger fractal with the results given in Table 10. The results are similar to the results obtained for the Sierpinski gasket (see Table 6) with non-trivial differences onD 3 andD 2 but zero onD 1 because of symmetric boundary conditions and the location of the (in this case) four holes.
Using the points illustrated in Fig. 3 onÊ k , transient temperatures are depicted in Fig. 15, obtained both fromT k (UoMFEC) and directly obtained with ABAQUS-direct with fine mesh. The average differences calculated by Eqs. (19) and (20) are provided in Table 11. The importance of including a discontinuity network is clear on examination of the results, with maximum errors appearing at corners. The temperatures produced by the continuous Galerkin method close to the corners are higher than those with a discontinuity network.
Comparing errors produced on the closed-pore structure (Sierpinski carpet) against the open-pore structure (finger-like fractal) reveals the latter structure's relative sensitivity to the inclusion of a discontinuity network.  This is not unexpected, since enforcing continuity between points separated across holes is difficult to justify and is unrepresentative for open-pore constructions of the type depicted in Fig. 3.

The influence of mesh
It is apparent from the results in previous sections that UoMFEC provides improved accuracy in its mapped results when compared with results obtained from ABAQUS applied directly to pre-fractals. A feature of the tessellated approach, however, is that the tessellated mesh used in any analysis is related to the initial tessellation onÊ 0 ; it is thus of interest to investigate the influence the initial tessellation has on the accuracy of predicted results. Consider again the Sierpinski gasket with an initial tessellation of 6 elements onT 0 =Ê 0 . Re-examination of Fig. 1 reveals tessellations onT 1 ,T 2 andT 3 consisting of 18 = 3 × 6, 54 = 3 2 × 6 and 162 = 3 3 × 6 tiles. The tile numbers arise from the number of maps involved (i.e. 3) and the initial number of tiles (i.e. 6). Each of these tessellations is produced by the expansion maps provided in "Appendix A". These tessellations could all in principle be used as alternative initial tessellations onÊ 0 . The expansion maps now associated with these initial tessellations onT 0 =Ê 0 are obtained by function compositions of the maps in "Appendix A", but there is no requirement to formulate these explicitly because all tessellations can be produced from the original set of maps. This aspect is demonstrated in Fig. 16, which shows tessellations produced solely from the maps in "Appendix A", without recourse to composite functions. The different initial tessellations onT 1 andÊ 1 in 16 are obtained from different initial tessellations onÊ 0 , themselves formed from a tessellation consisting of six tiles. Thermal analysis performed on the Sierpinski carpet in Sect. 5 is repeated but on the gasket with an internal heat loadingQ tot = 300kW/m 3 and cooling-channel heat transfer coefficients given in Table 12, for the various hydraulic diameters involved. Temperature distributions onT 3 mapped toÊ 3 as determined by UoMFEC are shown in Fig. 17. These are compared with results from ABAQUS-direct applied directly toÊ 3 for a similar numbers of elements as presented in Fig. 18. All the temperatures from UoMFEC (all UoMFEC results are lifted from the tessellation) provide very close agreement with ABAQUS-direct results, providing further confidence in the tessellation method and discontinuity networks. The average errors determined using Eqs. (19) and (20) on the edge y = 0 are given in Table 13. The results correspond to temperatures presented in Fig. 19.
The average error (contrasting UoMFEC against a convergent ABAQUS-direct solution) obtained using Eq. (19) is recorded to give an indication on how error reduces with mesh refinement. Mesh refinement is achieved through the increasing the number of tiles on the original setÊ 0 . The recorded errors onÊ k are depicted in Fig. 20, where convergent behaviour is observed.
Transient temperature profiles are shown in Fig. 21 along with average differences given in Table 14. Both steady-state and transient results confirm that accuracy improves with tessellated mesh refinement.

The influence of different hole-fill maps
One of the features of the tessellated approach is the non-uniqueness of the sub-expansion maps P i j in the formation of on an expansion map P i . This variability is reflected in changes in the tessellation onT 1 and changes in P i change the hole-fill map x : s → x and influence the precise manner in which holes are closed. Recall that hole-fill maps relate points s on a pre-fractal to points x in a tessellation. Any temperature T (x, t) calculated on a tessellation is related to a temperature on a pre-fractal through the hole-fill map x : s → x, i.e. T (s, t) = T (x (s) , t). Predictions on a pre-fractal are required to be invariant with respect to the non-unique choice of hole-fill map x : s → x. The explicit form the hole-fill map takes for a fixed set of contraction maps S i depends on the choice of expansion maps P i . It is of interest therefore to demonstrate explicitly the invariance of predictions for two significantly different hole-fill maps. It is also of interest to examine the effect distorted tiles have on the accuracy of predictions. These aspects are examined on pre-fractals for the classic Vicsek fractal, on tessellations formed with two distinct hole-fill maps.

Pre-fractal and tessellation construction
The Vicsek fractal is a non-product fractal, constructed by the recursive application of the affine contraction maps shown in Table 23 Thermal conductivity contours using Eq. (5) for point-wise maximum thermal conductivities are shown in Figs. 24 and 25. Examination of the two figures reveals significant differences in the distributions; it should be appreciated that differences in the tessellated space can be expected but of principal concern is the retention of isotropy in the physical space.

Thermal analysis of a Vicsek fractal heat exchanger
The open-pore Vicsek structure is assumed to be placed in a close-fitting container to enforce the flow of coolant through the structure. Consider then water coolant flowing through the voids shown in Figs heat exchanger is manufactured from the same copper material used in Sect. 5 with the exact same surrounding conditions. A uniform heat loading ofQ tot = 500 kW/m 3 is applied on each pre-fractal as an internal heat source. The cooling-channel heat transfer coefficientsĥ wat s are provided in Table 15 for the hydraulic diameters met on the Vicsek fractal.
Temperature distributions on pre-fractalÊ 3 and corresponding tessellationT 3 determined by UoMFEC from the two different tessellation maps are shown in Fig. 26. These can compared against results obtained with ABAQUS-direct for two meshes as depicted in Fig. 27. Examination of Fig. 26 reveals stark differences in temperatures on the tessellation, yet almost identical results in the physical space. This result provides good evidence that the manner in which tessellations are created does not affect the final outcome as seen on prefractals in the physical space. To quantify the results more precisely, temperature distributions along y = 0.5 are shown in Fig. 28 for both ABAQUS-direct and UoMFEC. It is clear that the two sets of tessellations return high accuracy. The errors calculated via Eqs. (19) and (20) are also given in Table 16 and also confirm high accuracy.
Transient plots for Point 1 in Fig. 22 with coordinates (0.5, 0.5) onÊ k are depicted in Fig. 29 along with the average differences calculated by Eq. (19) and (20) in Table 17. The results provided further evidence that different hole-fill maps have little impact on results viewed in the physical space.   It is worth noting that highly skewed tiles also have little impact on the accuracy of the results with the tessellated approach. This is expected because changes in thermal conductivities compensate to provide an analysis which is akin to doing calculations with the regular meshes appearing on the corresponding pre-fractals.

Conclusions
This paper tests the hypothesis that heat transfer analysis on porous media for cellular heat exchangers can be performed on a tessellation with new form of tessellated continuum mechanics using the finite element method. Thermal analysis of a relatively complex structures replicated by pre-fractals is demonstrated on corresponding tessellations incorporating discontinuity networks. In establishing the tessellated finite element approach for heat transfer on tessellations with discontinuity networks, the following conclusions can be drawn: • Temperatures on tessellated structures have been shown to return the temperatures on pre-fractals in support of the contention that heat transfer analysis on pre-fractals can be achieved through analysis on tessellations with high accuracy. • The accuracy of the tessellated approach with the finite element method can be significantly improved by incorporating discontinuity networks in tessellations. • The use of highly skewed tiles and alternative expansion maps has little impact on the accuracy of the procedure. • Increased accuracy can be obtained by refining the initial tessellation on the initial domainT 0 =Ê 0 , which can be achieved very efficiently through multiple application of affine expansion maps.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Appendix A: Maps for the Sierpinski gasket
The Sierpinski gasket is an example of a non-product fractal set and is constructed by the recursive application of the following three affine contraction maps: where each S i maps a triangular domain into a smaller triangular as deduced on inspection of Fig. 1.
The initial six tiles on the triangular domainÊ 0 provide a means to define the expansion maps for the construction ofT k . Each expansion map consists of a number of sub-maps with a number equal to the number of initial tiles onÊ 0 (6 in this case). Explicitly, the expansion maps for the Sierpinski are illustrated in Table 18.

Appendix C: Maps for the finger-like fractal
The finger-like fractal is constructed by the five affine contraction maps in Table 21; the associated expansion maps are defined in Table 22.    If y ≥ x and y ≤ −x + 1