Numerical and experimental analysis of a RVE condition for heterogeneous materials using BEM and thermal images

This work focuses on the numerical and experimental study of the effective thermal conductivity in materials with heterogeneous composition. The numerical analysis relies in the Boundary Element method (BEM) formulation assisted by the sub-regions technique to simulate a two-dimensional square domain, with randomly distributed material inclusions. The boundary conditions are set to achieve a unidirectional heat flux condition. The Average Field Theory (AFT) is applied to determinate the representative volume element (RVE) condition for the studied cases. Finally, the effective heat conduction coefficient is obtained from the defined RVE’s. To verify the obtained numerical results, it is proposed an experimental procedure based on thermal images. The experimental assembly replicates the unidirectional heat flux condition over a steel plate. The experimental effective thermal conductivity is obtained from the analysis of the resultant temperature field on the surface of the plate taken by an infrared camera. The comparison between the numerical and experimental results showed a 3% difference between both results, pointing out for the successfulness of the proposed methodology.


INTRODUCTION
In face of the current fast-paced technological advances, the study and development of new high performance materials has been a requested field of research lately.The development of new "taylor-made" materials optimized for specific applications requires the understanding of the material composition and the effects of its microstructure on its macroscopic properties.
At a microstructural level, engineering materials are typically heterogeneous with many possible constitutional elements, each with its own physical properties.A volume element, defined by the observational scale, may contain just a few microconstituents.In these cases the effective physical properties of this element may be largely influenced by the individual properties of the microconstituents contained within it.
An industrial application of the effective properties concept is presented in [1].This work proposes a study of the effective thermal conductivities of porous materials applied on the fabrication of Vacuum Insulation Panels (VIPs) using simplified cell models.Also considering other heat transfer mechanisms, such as gaseous and radiative conductivities, the elaborated mathematical models provided relatively accurate results compared to experimental data for every studied type of material, except for powders.This evinced the limitation of simplified models and the need for other methods for modeling materials with complex microstructures.
In face of this, [2] presented a study on the development of correlational relations for determining the effective thermal conductivity of two-phase materials.Based on experimental data, it succeeded in developing correlational relations for three specific two-phased problems.This work points out the need for more accurate and efficient methods for prediction of the thermal properties of heterogeneous composition materials.
In another work, [3] evinced that most of the available methods for the prediction of material properties are based on idealized periodic constructions.It proposed a computational modeling of open-cell foams with random cell generation and growth, making possible the analysis of more complex and realistic porous structures.
As the observational scale is increased, the observed volume element effective properties gradually approach the material macroscopic physical properties.At these observational scales the influence of individual microconstituents tend to be less impacting on the element overall effective properties, as more microconstituents compose it.This variant behavior persists until a certain observational scale is reached, and the effective properties of the volume element approximate the material's macroscopic properties.
In this way, [4] proposed a homogenization method of materials with heterogeneous composition, the representative volume element method.It consists in determining a characteristic length for the observed volume element at which the effective properties measured over the volume element boundaries approach the macroscopic properties of the material.This volume element with the characteristic length is named as the representative volume element (RVE).This natural compatibility was explored in [5] and [6] where both methods, the BEM and RVE, were brought together to elaborate methodologies for the prediction of effective properties of micro-porous materials.The first one, [5], focused in the elastic effective properties of isotropic and orthotropic porous materials determining the RVE conditions through statistical analysis.In the same way, [6] applied a similar methodology for the analysis of permanent potential problems and the determination of the effective thermal conductivity coefficient for micro-porous materials.
Furthermore, [7] presented a modeling of a bidimensional heat transfer problem of plate with several random generated holes to simulate a RVE condition for a micro porous material.The numerical analysis proceeded using a "Fast Multipole" BEM formulation for increased computational efficiency.This work also presented an experimental analysis based on thermo graphical imaging to validate the obtained numerical results.The confrontation between numerical and experimental results showed a difference below 5%, indicating a good performance of the presented predictive methodology.
Following that idea, [8] presented a RVE study of the effective thermal conductivity property of materials composed by two different microconstituents with different overall properties.In this case, a BEM formulation with sub-regions was applied to consider the interaction between two different microconstituents that compose the domain.In this work, the RVE condition was obtained for different fractions of area occupied by less conductive circular inclusions randomly distributed over a square domain.The behavior observed on the graphical results showed that the RVE condition could be achieved at about 37 inclusions present within the studied domain.
The present work follows as a natural continuation of [8], proposing an experimental methodology to evaluate the RVE condition and validate the results shown in the aforementioned work.Aiming that, an experimental assembly is proposed to replicate the numerical case of a RVE represented by a square steel plate with randomly distributed circular inclusions subjected to a unidirectional heat flux condition.Thermal images are used to obtain the temperature fields over the plate surface.Finally, the effective thermal conductivity of the RVE is calculated based on these temperature fields.

REPRESENTATIVE VOLUME ELEMENT
In a microscopic scale, several different microconstituents compose typical engineering materials such as steel and aluminum.Despite the isotropic behavior of these materials, when observing a microscopic defined region, the mechanical properties measured at the edges of it may present different values for each direction.This occurs because, at this observational scale, the domain is composed by just a few micro constituents and hence, their geometry, positioning and own mechanical properties largely influence its effective properties.Increasing the observational scale also increases the number of microconstituents located within the observed domain, as illustrates Figure 1.As more microconstituents become present in the larger domain the individual influence of a single microconstituent on the domain's effective properties decreases.
Every observational scale is defined by a characteristic length, which corresponds to the length of the edges of the square domain (L).In [5] and [6] for microporous domains, and in [8] for domains composed by two different microconstituents, it is noticed that after successive increases on the observational scale, the effective properties measured tend to converge towards the material's macroscopic overall properties.Based on this behavior, [4] proposed a homogenization method for materials with heterogeneous composition, the RVE technique.It consists of defining the smallest observational domain at which the effective properties taken at its edges show sufficient convergence to the material's macroscopic properties.
In this sense, [8] presents a RVE determination methodology for materials composed by two different microconstituents.As done previously in [5], [6] and [7] it supposes a square observational domain, but with circular material inclusions instead of holes.Based on equation (1) and statistical analyses, the RVE condition for the effective thermal conductivity was obtained for different percentages of area occupied by the inclusions (R), being R given by: R = 100 At the end, despite the applied value of R, Oberg et al.
(2015) concluded that, for the studied cases, the RVE condition may be assumed at around 37 inclusions.

BOUNDARY ELEMENTS METHOD
The BEM formulation for potential heat transfer problems is, currently, one of the most well diffused BEM formulations.It has been featured in several works since the method was firstly developed and its detailed elaboration can be found in many books such as [9], [10] and [11] among others.
To briefly summarize the BEM formulation for potential problems using constant elements, an initial domain, as depicted on Figure 2, is assumed.The method derivation starts at the Laplace governing equation for 2D potential problems given as: Figure 2. Domain Ω and its boundary Γ.
For potential problems, the boundary conditions applied to its boundary may be of three types: Dirichlet, Neumann and/or Robin.This work considers only the first and the second types, respectively, imposed as: where u is the potential field in domain Ω, Γ is the boundary of Ω, and n is the outward normal.Note that the barred quantities are the values imposed by the boundary conditions.Given the boundary conditions, the solution of equation ( 3) is presented as: where and are the Green's functions for 2D problems while r represents the distance between the collocation point x and the field point y, as depicted in Figure 2. When x is taken to the boundary, the classic boundary integral equation (BIE) formulation of BEM is obtained as: Considering a smooth boundary at the collocation point x, the coefficient c(x) is assumed as 1/2.For the computational implementation of the method, the next step consists in discretizing the boundary Γ using N constant elements.The following equation presents the discretized BIE.
In this, and (j = 1,2,…,N) are the nodal values of u and q at the element ΔΓ j , respectively.Applying the boundary conditions at each node and switching the columns for grouping the unknown variables, equation ( 6) can be rearranged as linear system of equations.
The matricial representation of this is given as: Where A is the coefficient matrix, l is the unknown vector and B is the known right-hand side vector.

BEM WITH SUB-REGIONS
In the studied case, the sub-regions are delimited by closed boundary regions placed inside a main domain, as shown in Figure 3.In this example, a circular sub-region (Ω') is placed inside a square domain (Ω).The boundary elements applied to discretize the sub-regions are common to bothΩ' and Ω.To better illustrate this, Figure 4 separates the domains depicted on Figure 3 to evince the shared elements.In this figure, elements i (i = 1,2,3,4) and its respective elements i' are coincident when assembled together.Each of these shared elements adds two new unknown variables to the linear system.To deal with these new variables and couple the different regions it is necessary to impose a pair of equations to enforce the continuity condition.Therefore, for each shared element: where the variables u and q are the variables of the shared element related to one region and u' and q' are the same variables related to the equivalent element in the other region.Adding equation ( 8) to equation ( 6) satisfies the solving criteria, balancing the number of unknown variables and available equations.

EFFECTIVE THERMAL CONDUCTIVITY
The effective properties of a RVE are obtained based on the information gatherewwd over its boundaries.
Reached the RVE condition, the effective properties tend to represent the macroscopic properties of the material as a whole.In sight of this, the study of these properties and methodologies assoc iated to their determination have been the focus of many recent works.These effective properties are often related to mean values of distributions over the RVE edges, as observed in [5], [6] and [8].In [6] and  In the end, after obtaining the local heat conductivity coefficient for each pair of points, the RVE effective heat conductivity is calculated as the mean value of the local conductivities.This entire process is described by equation ( 9), where is the RVE effective heat conductivity coefficient, is the imposed heat flux, is the number of pair of points and L is the RVE characteristic length.

NUMERICAL ANALYSIS
In accordance with [8], the numerical procedure consists in numerically simulating the unidirectional heat flux condition using the BEM formulation with sub-regions.The object of analysis is a square steel plate with randomly distributed less conductive round inclusions.The size and quantity of inclusions are controlled by equation (1).In each performed study, R was kept constant, while the diameter of the inclusions (d) was calculated as a function of number of present inclusions (n).Therefore the act of simply increasing the value of n produces an effect similar to the increase of the observational scale depicted in Figure 1.
Established the simulation parameters, the next step is determining the influence of the observational scale over the resultant effective thermal conductivity.This is done by statistically evaluating the dispersion of the resultant effective thermal conductivity for each value of n as a standard deviation.Is this sense, as presented in [8], for each value of n, 34 randomly generated domains are analyzed to guarantee the statistical stability of the results.Figure 7 graphically demonstrates the variation of the calculated mean effective thermal conductivity for R = 0,20 and different values of n.In this figure, the drawn line conductivity and the vertical bars represent, respectively, the resultant mean values of the effective thermal conductivity and its associated dispersion intervals.From the graphical analysis, it is observed that, initially, the intervals tend to become smaller as more inclusions are inserted.This behavior persists until around n = 37.At this point the interval size becomes stabilized characterizing the RVE condition.
Repeating this same analysis for different values of R, [8] concluded that the RVE condition may be assumed at around 37 circular inclusions, regardless the value of R.

EXPERIMENTAL ANALYSIS AND RESULTS
Based on [7] and on the numerical results presented in [8], an experimental assembly was prepared to validate the presented numerical methodology.
The main goal of the experiment was to recreate the unidirectional heat flux condition pictured on Figure 5.
In order to reproduce the theoretical studied domains, it was used square SAE 1020 steel plates (155 x 155 x 9,52 mm) with 37 randomly 12,7 mm diameter holes.These holes were filled with same diameter nylon 6.6 cylinders.This distribution corresponds to an R value of 0,195.Being significantly less conductive than the steel, the nylon inclusions are expected to disturb the imposed unidirectional heat flux condition and produce results similar to the numerically observed.To assure the heat conduction between the inclusions and the plate, a high performance thermal compound was applied to every contact region.The top of the plates were painted black to increase the uniformity of the thermographic readings being taken.
A high watt density cartridge heater (220V) was used to uniformly heat one of the sides of the plate.The heater was build inside a refractory brick to prevent it from damaging the other experimental components in case of overheat.The resistance was built inside a refractory brick, leaving only the surface in contact with the experimental plate not insulated, reducing the heat loss in every other direction.This also prevents the heater from damaging other experimental components in case of overheat.The opposed plate side was cooled by contact with an aluminum recipient filled with water and ice kept at melting temperature.Thick 100mm Styrofoam walls were used to insulate the two remaining sides of the plate and its back to minimize the heat loss in these directions.The proposed experimental assembly scheme is shown on Figure 8.To evaluate the heat loss through the refractory brick holding the cartridge resistance, a solid steel plate (155 x 155 x 9,34 mm) is mounted on the experiment and heated until a permanent heat conduction condition is achieved.The thermal images of the top surface are taken using a thermographic camera (FLUKE Ti 125), allowing for relatively precise readings of the temperatures on the cooled and heated edges of the plate.Figure 9 shows the thermal image taken of the solid steel plate along with the maximum, minimum and mean temperatures taken at 3 specific drawn lines: (a), (b) and (c).The temperature distributions over lines (a) and (b) represent the temperature of the cooled and heated edges respectively, while (c) was traced to evaluate the expected linearity of the temperature gradient over the plate.Figure 10 displays the graphical analysis of the temperature distribution over the traced lines.Minor disturbances shown on the graphic lines may be associated with the thermal camera precision and possible non-uniformity of the applied black paint coat.
Given the steel heat conductivity coefficient () of 51.9 [W/m.k] and the mean temperatures at the heated () and cooled () edges, the heat flux through the steel plate () are calculated from equation (10).
where A corresponds to the plate transversal area of the plate while L to its side length.
The comparison between the total expected heat flux provided by the resistance and the calculated heat flux () determined a 20% heat loss through the refractory brick for the tested case.
The same experimental procedure is now repeated with steel plate with the randomly distributed nylon inclusions.Figure 11 presents the thermal image of the temperature field over top surface of the plate.
As done for the solid steel plate, 3 lines were traced over the thermal image to acquire the temperature distributions over the cooled edge (a), the heated edge (b) and the influence of the inclusions on the overall temperature field.The temperature distributions for this case are graphically shown on Figure 12.
The proposed experimental methodology resulted in an effective thermal conductivity of 35,9 ± 4 W/m.K.
The previously proposed numerical methodology based on BEM with sub-regions is applied to replicate the experimental RVE condition.Following the numerical procedure, 34 cases with the same experimental 0,195 R value and 37 randomly generated inclusions are simulated to obtain the RVE effective thermal conductivity with its associated statistical dispersion.The mean values of the heated and cooled edges measured with the thermal images are applied as boundary conditions.This analysis resulted in an effective thermal conductivity of 34,8 ± 0,2 W/m.K.
Confronting the numerical and experimental results it is noticed a 3% difference between both mean results, pointing out for a good performance of the numerical methodology.

Figure 1 .
Figure 1.Process of increasing the observation scale.

Figure 3 .
Figure 3. Example of a sub-region problem.

Figure 4 .
Figure 4. Sub-regions problem as two separate regions.

[ 8 ]
, a similar technique was applied to obtain the material's effective thermal conductivity.It consists of imposing a unidirectional heat flux over a square domain, illustrated below, and reading the temperature distributions over the heated and cooled edges.The temperature distributions are composed by temperature readings from several points equally distributed over the edges of interest.By horizontally pairing punctual temperatures on the heated (T ci ) and cooled (T hi ) edges it is possible to calculate the local heat conductivity coefficient between the pair of points.The following Figure depicts the pairing process for a better understanding.

Figure 7 .
Figure 7. Mean effective thermal conductivity vs. number of inclusions inserted for R = 20%.

Figure 9 .
Figure 9. Temperature distribution over the solid steel plate.

Figure 10 .
Figure 10.Temperature distributions over lines (a), (b) and (c) for the solid plates.

Figure 11 .
Figure 11.Temperature distribution over the steel plate with 37 randomly distributed nylon inclusions.