Al composites

Advances in Engineering Software 77 (2014) 1–12

Contents lists available at ScienceDirect

Advances in Engineering Software journal homepage: www.elsevier.com/locate/advengsoft

An efficient numerical model for investigating the effects of anisotropy on the effective thermal conductivity of alumina/Al composites W. Leclerc ⇑, N. Ferguen, C. Pélegris, E. Bellenger, M. Guessasma, H. Haddad Université de Picardie Jules Verne, Laboratoire des Technologies Innovantes, 48 rue d’Ostende, IUT de l’Aisne, 02100 Saint Quentin, France

a r t i c l e

i n f o

Article history: Received 16 April 2014 Received in revised form 19 June 2014 Accepted 5 July 2014

Keywords: Composite materials Thermal conductivity Granular systems Fast Fourier transform method Numerical homogenization Discrete Elements

a b s t r a c t The paper describes an efficient numerical model for better understanding the influence of the microstructure on the thermal conductivity of heterogeneous media. This is the extension of an approach recently proposed for simulating and evaluating effective thermal conductivities of alumina/Al composites. A C++ code called MultiCAMG, taking into account all steps of the proposed approach, has been implemented in order to satisfy requirements of efficiency, optimization and code unification. Thus, on the one hand, numerical tools such as the efficient Eyre–Milton scheme for computing the thermal response of composites have been implemented for reducing the calculation cost. On the other hand, statistical parameters such as the covariance and the distribution of contact angles between particles are now estimated for better analyzing the microstructure. In the present work we focus our investigations on the effects of anisotropy on the effective thermal conductivity of alumina/Al composites. First of all, an isotropic benchmark is set up for comparison purposes. Secondly, anisotropic configurations are studied in order to direct the heat flux. A transversally isotropic structure, taking benefit of wall effects, is finally proposed for controlling the orientation of contact angles. Its thermal capabilities are related to the current issue of heat dissipation in automotive engine blocks. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction Thermal management of automotive engine blocks is a key issue on which depends the durability of certain polymer parts, electronic devices and embedded software which are more exposed to high temperatures. Recent European Union legislations restrict pollutant emissions of light vehicles. As a result, a reshape of engine block has become necessary in order to meet the new European standards. Engine downsizing is the most natural solution for reducing CO2 emissions but this requires a boosting device (such as a turbocharger or a supercharger) and a direct injection technology for maintaining the capacity of the engine block. Resulting heat stresses are sizeable and new materials are required for dissipating the heat flux. Recently, a manufacturing process was developed for elaborating alumina/Al composites with enhanced thermal capabilities [1]. The process is based on the method introduced by Descamps et al. [2,3] for the synthesis of organic porogen materials. Alumina/Al composites are composed of a ceramic phase which ensures the insulating character while the thermal conduction is ensured by aluminum spherical particles. Thus, such a composite could well be a potential candidate for the thermal ⇑ Corresponding author. Tel.: +33 323 503 697; fax: +33 323 503 698. E-mail address: [email protected] (W. Leclerc). http://dx.doi.org/10.1016/j.advengsoft.2014.07.004 0965-9978/Ó 2014 Elsevier Ltd. All rights reserved.

management of engine blocks. However, the manufacturing process leads to near-isotropic materials which are not able to direct the heat flux and manage the temperature gradient. That is why some solutions have to be found for modifying the elaboration of Alumina/Al composites in such a way as to generate directional effects. Different strategies are available for this purpose including functionally graded materials [4–7], and the control of the orientation of contact angles between particles according to specific loadings or modes of deposition [8]. This paper targets the last option. In a previous work, a numerical model was set up for predicting the effective thermal properties of alumina/Al composites [1]. On the one hand this approach uses the discrete element method (DEM) for reproducing the sedimentation process leading to the generation of granular systems composed of spherical aluminum particles. On the other hand, effective conductivities are assessed using the Finite Element Method (FEM). The model was initially used for investigating the influence of morphological and phenomenological parameters such as interconnection sizes, interfacial debonding and granulometry. Results confirmed the critical impact of the two first but no granulometry effect was exhibited. A C++ code called MultiCAMG is currently being developed for the generation and the analysis of granular systems using the numerical approach of Ferguen et al. [1]. The code satisfies some needs of efficiency, optimization and code unification, which are not met by the use of

2

W. Leclerc et al. / Advances in Engineering Software 77 (2014) 1–12

a number of commercial software applications. It takes into account all steps of the numerical approach, namely the generation of granular systems, the discretization and the evaluation of the Effective Thermal Conductivity (ETC). It also permits the connection with other codes or softwares such as R one for statistical analysis purposes [9]. The lack of efficiency of the FEM leads us to consider the well-established Fast Fourier Transform (FFT) based method [10,11] as an alternative way for evaluating effective properties. Efficient voxelization methods [12] have been adapted and used for handling such an approach using the fast Eyre–Milton scheme [13]. In the present work, our objectives are twofold. First, we extend the numerical approach of Ferguen et al. [1] to FFT-based methods and statistical analysis using MultiCAMG. Two kinds of systems are considered, namely random packings composed of unpenetrable hard spheres, and granular systems composed of partially penetrable spheres for which overlappings are enabled and controlled. We investigate effects of the volume fraction of spherical particles, the orientation of contact angles and the distribution of particles on the ETC of alumina/Al composites. An isotropic configuration is set up for comparison purposes. We assume the material to be isotropic when the distribution of particles and the orientation of contact angles are statistically independent of the orientation of chosen axes [8]. Furthermore, in the case of random packings, the volume fraction of spherical particles is set to 0.64 which approximately corresponds to the random close packing (or maximally random jammed packing) for monodispersed spherical particles [14–18]. Several numerical and statistical tools such as the covariance [19] and the radial plot [8,20] are set up and used for better characterizing the granular system in such a context. ETC are estimated by FFT and reference graphs are plotted as well. Secondly, anisotropic configurations are studied with the same numerical and statistical tools. Some lattice structures are investigated and compared to the isotropic system. A composite composed of chains of aligned particles is also described and analyzed. Such a structure is transversally isotropic and thus has a higher thermal conductivity along the preferred axis. However this is particularly difficult to design. That is why, finally, a more realistic methodology is proposed for controlling the orientation of contact angles using inner walls in order to take benefit of wall effects [21]. Parametrization, ETC and capability to direct the thermal flux in this framework are discussed, and related to the issue of thermal management in engine blocks. The paper is outlined as follows. First, we list the various steps in the manufacturing process. A third section is then dedicated to the numerical approach and MultiCAMG code. The fourth section describes the setting up and the investigation of isotropic granular systems. The last part deals with anisotropic configurations and gives the results of the proposed approach for controlling the orientation of contact angles. 2. Manufacturing process The process used for manufacturing a composite with graded microstructure is based on the porogen organic agents method developed by Descamps et al. [2,3] for the synthesis of macroporous b-tricalcium phosphate. The main idea of the technique consists of creating an alumina foam with controlled macroporosity, and infiltrate it with a melt aluminum alloy. Pore and interconnection sizes are critical in determining the properties of the final composite material. The main steps of the process are described below. 2.1. Producing of the organic frame A polymeric scaffolding composed of polymethylmethacrylate particles (PMMA) is used as precursor foam. The granular system is composed of spherical particles with fixed granulometric distri-

bution. Particles are chemically stuck within a metallic mold using a solvent and by squeezing the particle skeleton. The applied pressure slightly deforms the particles so that necks appear between PMMA particles [2,3]. Thus, the resulting network of spherical particles seems to be interconnected at contact points but no real overlapping arises here. The shrinkage of the granular system is measured by Scanning Electron Microscope (Hitachi S-3500N) in order to correlate the height of the bed to the pore interconnectivity of the organic frame. 2.2. Synthesis of the ceramic foam with controlled macroporosity Before starting the impregnation, an alumina slurry is made from mixed powder with varied concentrations in the range of 65–75 wt.%. When the organic frame and the alumina slurry are ready the impregnation process is carried out. After drying, a debinding treatment is then used at low temperature (1 °C/min, 30 h at 200 °C). This allows PMMA particles to be removed and generate macroporous foam with open pores whose dimensions depend on interconnection sizes of PMMA particles. Subsequently, the ceramic foam is sintered (5 °C/min, 1 h at 1670 °C) until consolidation of the ceramic walls occurs [2,3]. 2.3. Infiltration of the ceramic foam The infiltration process of the melt aluminum inside the ceramic foam is achieved in vacuum condition with temperature held at 720 °C. The fluidity of the melt aluminum inside the foam strongly depends on the wetting properties of the ceramic surfaces. Magnesium is added to the aluminum alloy in order to enhance the fluidity of the melt alloy and improve the filling of pores. Steps of the manufacturing process discussed above are summarized in a flowchart as shown in Fig. 1. 3. Numerical modeling A C++ based code called MultiCAMG (Multifunctional Code for the Analysis and the Modeling of Granular Packings) has been implemented for simulation purposes. This enables the generation and the study of random systems composed of spherical particles using the Discrete Element Method (DEM) on the one hand. The program also estimates the ETC of composites using either FEM or FFT approach. 3.1. Generation of granular systems The architecture of the alumina/Al composite is simulated considering a granular system composed of spherical particles. The particle skeleton is generated using the DEM which is a wellestablished approach to take into account both velocity of particles and inter-particle collisions [22,23]. A spring and dash-pot contact model is used to describe particle–particle and particle–wall collision. The normal force between two particles i and j in contact is given by Eq. (1) which allows small overlappings of rigid particles to occur at contact points. In Eq. (1) the shear component is neglected and no tangential displacement is considered.

pffiffiffiffiffiffiffiffiffiffi F n ¼ K n U n þ 2 mK n V n

ð1Þ

where n is the unit normal vector to the contact plane. K n and U n are respectively the normal stiffness and the displacement in the norpffiffiffiffiffiffiffiffiffiffi mal direction. 2 mK n is the critical damping constant and V n is the relative velocity at the contact between particles i and j. In the case of a particle–particle contact, m is given by Eq. (2).

m¼

mi mj mi þ mj

ð2Þ

3

W. Leclerc et al. / Advances in Engineering Software 77 (2014) 1–12

Fig. 2. Visualization of a 3D granular system: (a) at the end of sedimentation and (b) after vertical compression (contact force chains in black).

where Dab is the interparticle distance, and Ra and Rb are the radii of particles. The contact area of interconnected particles is characterized by the radius of the disk resulting of the intersection of two spheres. From now on, this radius is designated as the contact radius Rc and is estimated using the interconnection angle ha , as illustrated in Fig. 3:

"

D2 þ R2a  R2b ha ¼ arccos ab 2Dab Ra

Fig. 1. Manufacturing process flowchart.

# ð4Þ

Rc ¼ Ra sinðha Þ where mi and mj are respectively the mass of particles i and j. However, in the present context, DEM is only used for generating granular systems and not for evaluating thermal or elastic properties. That is why, K n is manually handled in order to control the overlap between each pair of particles [1]. The generation of 3D granular systems is performed in two steps. First, a set of particles with fixed radius are deposited by sedimentation under the effect of gravity in a rectangular container until equilibrium is reached. This latter is obtained when the ratio between the total kinetic energy and the mean mass of particles, is less than 0.001 m2/s2. The resulting granular system has a low volume fraction due to voids located near the walls of the container. The interconnection areas at contact points are, moreover, very small. That is why, in a second step, a movable top wall applies an external force in the z direction in order to increase the contact area and fill the voids. Fig. 2 illustrates the process and the state of the granular system before and after the compaction. The final interconnection areas strongly depend on the stiffness K n . In a previous work [1], charts were drawn for connecting K n to the overlap. 3.2. Parametrization Large samples of granular systems are efficiently generated with the help of the previous methodology. For each one, spatial positions and radii of spherical particles are recorded in a text file for further investigations. Indeed, several parameters are of interest and require specific studies. In the sequel, more emphasis is devoted to the contact radius, the volume fraction / of particles and the isotropy of the granular system. The volume fraction is defined as the ratio between the volume of aluminum phase and the total volume of RVE. The overlapped distance U n between two spheres A and B is determined as follows:

U n ¼ Ra þ Rb  Dab

ð3Þ

ð5Þ Rc

However, in what follows, a dimensionless contact radius (Eq. (5)) will be considered since no real dimension is considered in the present work.

Rc ¼

2Ra sinðha Þ Ra þ Rb

ð6Þ

Numerical and statistical tools have been implemented in the C++ code for better investigating the scope of parameters. Thus, the spatial distribution of particles, and the distribution of contact angles within the granular system are described using the covariance and the radial plot, respectively. These will be discussed in Section 4. 3.3. Evaluation of the ETC 3.3.1. Statistical RVE and thermal conductivities The numerical model is composed of two parts. The first part represents the system of aluminum spheres and is generated using the DEM approach. The second part represents the porous matrix and is obtained as a negative image of the particle skeleton. Both assembled parts form a cubic cell called the Representative Volume Element (RVE) of the composite [24]. Its dimensions are of crucial interest since they strongly affect the accuracy of the numerical calculations. Typically, the RVE has to be large enough to take into account a sufficient amount of information on the microstructure. This assumption ensures the cubic cell to be representative of the material. However, though larger dimensions lead to more accurate results, they also affect the calculation time. A more efficient approach consists in generating a sample of small Volume Elements (VE) the size of which is not sufficient to ensure the representative feature of the RVE but leads to more efficient calculations [12]. In this context, effective properties are asymptotically obtained by taking the average of apparent properties assessed for each VE. Thus a good compromise is found between

4

W. Leclerc et al. / Advances in Engineering Software 77 (2014) 1–12

Contacts

Fig. 3. Definition of the interparticle contact conditions in the neighborhood of particles.

efficiency and accuracy. However one has to know that a sample of too small VE provides erroneous predictions. That is why a statistical RVE with large enough dimensions for ensuring the validity of the sampling process has to be determined. In a previous work [1] in which the present numerical approach was proposed and validated, the size of the statistical RVE was investigated and determined as being 12 times the particle radius for monodispersed particles. This ratio corresponds to a cubic cell composed of about 300 spherical particles. From now on, unless otherwise specified, all calculations are performed with statistical RVE composed of 300 or more particles. The properties of each phase, alumina and aluminum, are summarized in Table 1. Furthermore, the term ‘‘statistical RVE’’ will be deliberately replaced by ‘‘RVE’’.

3.3.2. FFT The FFT approach is an efficient numerical technique for estimating effective properties. This is based on the Lippmann–Schwinger equation which is related to a classic linear problem via a Green’s operator. It deals with working in the Fourier space in which the convolution operator, which arises in the Lippmann–Schwinger equation, becomes a simple classical product. Numerous iterative schemes have been proposed for solving this kind of problem among which the pioneered algorithm of Moulinec and Suquet [10] and its improved version [11]. However, according to recent studies [26], the Eyre–Milton scheme [13] exhibits the best efficiency. From a practical point of view, an FFT-based approach has the great advantage to not require any mesh but is restricted to periodic boundary conditions. In the present work, three subroutines were implemented in our C++ code for evaluating thermal properties using the Eyre–Milton scheme [13] and the modified Green’s operator provided in [27]. The first subroutine discretizes the RVE into a regular voxel grid on which the FFT is based. An example of such a discretization with a 64  64  64 resolution, is depicted in Fig. 4. The illustrated periodical RVE is composed of a packing of 300 unpenetrable spheres and the volume fraction of particles (/) is 0.6331. The second one ensures the periodicity of the cell which is required by the FFT approach. Two cases have to be distinguished. In the first one, for which RVE are non-periodical, the unit cell is symmetrically duplicated in x, y and z-directions so that two opposite faces of the final cell perfectly coincide. In the second case of a periodical RVE, a direct discretization leads to isolated unmatched voxels on opposite faces due to the approximation of the geometry. A treatment is carried out for ensuring the

Table 1 Properties of each phase.

Alumina Al

Specific heat (kJ kg1)

Thermal conductivity (W/m K)

Mass density (kg m3)

850 875.2

23.4 116

3960 2780

Fig. 4. Example of a periodic discretized granular packing with a 64  64  64 resolution ð/ ¼ 0:6331Þ.

coincidence by giving a priority to voxels labeled as aluminum. The third one consists of the Eyre–Milton numerical scheme [13]. A major issue is related to the geometric approximation which leads to both a discontinuous interface and an approximate volume fraction of heterogeneities when the resolution is insufficient. This problem drastically affects the validity of predictions and is only solved using a sizeable resolution. A rigorous study has to be set up for determining the suitable resolution leading to accurate results. First, the apparent volume fraction of particles /a which is measured by dividing the number of voxels labeled as aluminum by the total number of voxels, must coincide with the real volume fraction /. Preliminary tests exhibit that 16  16  16 and 32  32  32 resolutions lead to very exaggerated values which cannot be considered. A 64  64  64 resolution is more suitable but the apparent volume fraction is still slightly too high (/a ¼ 0:6389 against / ¼ 0:6331 in the example depicted on Fig. 4). Thus, only 128  128  128 and finer resolutions are suitable for calculations. Fig. 5 illustrates the influence of both resolution and number of particles on the ETC. The Lubachevsky– Stillinger algorithm [28] is here used for generating a large sample of packings. Only monosized particles are considered and the real volume fraction of spherical particles is set to 0.63350:0005 for comparison purposes with FE results obtained in a previous work [1]. The choice of Lubachevsky–Stillinger algorithm is motivated by the inadequacy of the DEM approach to efficiently generate a large sample of patterns with an identical volume fraction. However one has to keep in mind that granular systems will only be

W. Leclerc et al. / Advances in Engineering Software 77 (2014) 1–12

5

and the sedimentation process described in 3.1. In this context, these ones are not periodical and prone to wall effects which affect the isotropy of the medium. First, several morphological features as the volume fraction of particles, the distribution of contact angles and the spatial distribution of heterogeneities are investigated in order to exhibit these drawbacks. Some solutions are tested for ensuring the isotropy of the medium. Second some assessments of ETC obtained from the method presented in 3.3 are provided and the study is extended to granular systems composed of penetrable spheres. 4.1. Morphological features

Fig. 5. Influence of both resolution and number of particles on the ETC ð/ ¼ 0:6335  0:0005Þ.

generated using the DEM from the Section 4. Each data is obtained from a set of 30 calculations performed with the same RVE size and resolution. Properties of each phase are those given in Table 1. The main observations are as follows. The ETC depends on both resolution and number of spheres. Thus, the convergence is only reached for a 256  256  256 discretization for which the predicted ETC is 66.360.15 W/m K for 300 particles. This value is in agreement with the FE prediction given in [1] for a particle in the range of radius from 900 to 1000 lm, namely 66.72 W/m K. One can also notice that, from some point, the greater the number of particles is, the more inaccurate the ETC is. This is due to the loss of geometric accuracy when considering a greater number of inclusions for a given resolution. From this range of results, we conclude that a resolution of 256 (or a finer one) is necessary for ensuring the exactness of numerical predictions. For this purpose, a number of spheres of 50 or 100 is sufficient for evaluating the ETC. However, we will consider a number of 300 particles in order to respect the RVE size determined in a FE approach. 3.3.3. FE An FFT approach is particularly efficient for evaluating effective properties of isotropic media. However, to our knowledge, the FFT is not easily usable when considering an anisotropic heterogeneous material. That is why, in such a framework, we consider the FE tool which is less efficient but much more flexible. FE calculations are performed under a steady-state heat flux. A temperature difference DT (K) is imposed between two opposite faces of the RVE and the other ones are assumed to be adiabatic. The total heat flow rate Q_ (W) resulting from heat conduction is numerically estimated and related to the temperature difference. Thus, the ETC keff is given by the Fourier’s law as follows,

keff ¼

Q_ l A DT

4.1.1. Volume fraction of aluminum spheres In a random configuration, many experimental and numerical measurements point to a volume fraction of around 0.64 [14,16,29]. This value is called the random close packing volume fraction and has to be checked out in a numerical model for validation purposes. In the present work, we consider systems composed of 300 spheres according to the discussion done in 3.3 about the issue of the RVE. Besides, in this subsection, interconnection between two or more particles is not permitted. The mean volume fraction is estimated to 0.582. This value turns out to be much lower than the random close packing volume fraction. This is because wall effects drastically affect the volume fraction of the granular system [21]. Such a drawback can be circumvented by considering a greater number of particles. Some tests were carried out up to 1500 particles and exhibit some improvements but, in the best case, the volume fraction was only 0.614. An efficient solution consists in considering free periodic boundaries [30] using fast algorithms [28]. However, such approaches affect the sedimentation process in such a way that both contact angles and contact radii are not well-controlled. Consequently, in order to be more consistent with the experimental process, we prefer another approach which consists in extracting a cubic RVE from an initial granular system. The extracted cubic cell is centered in the middle of the system, and its faces are parallel to the ones of the container. Fig. 6 illustrates the method which is done using a parameter q as follows:

q¼

Sc Sbox

ð8Þ

where Sc and Sbox respectively refer to the size of the cutting and the size of the box. q describes the rate of material removal in each main direction. Thus, voids located near the walls are eliminated and the volume fraction is enhanced. Fig. 7 shows the volume fraction estimated for q = 0, q = 0.05 and an extraction of 300 spheres from a denser granular system. For information purposes, the last case is performed considering an adapted q parameter which

ð7Þ

where l (m) and A (m2) are respectively the thickness and the crosssection of the material. More precisions about the numerical calculations and the determination of the RVE are provided in [1]. 4. Isotropic granular systems This section is devoted to the setting up of an isotropic benchmark for subsequent comparisons. Granular packings composed of unpenetrable spheres are generated using the DEM approach

Fig. 6. Scheme of the extraction method.

6

W. Leclerc et al. / Advances in Engineering Software 77 (2014) 1–12

depends on the initial number of particles. A random loose packing is obtained when no extraction is realized. Conversely, for q = 0.05, the volume fraction does not depend on the number of spheres and is quite close to the random close packing volume fraction, namely 0.64. In the third case, the volume fraction values are comparable to the previous ones as soon as the number of particles is greater than 1200 entities. 4.1.2. Contact angles between particles From a geometrical point of view, isotropy is confirmed when no directional effect arises. In such a configuration, morphological parameters such as the spatial positions and the contact angles between each pair of particles have to be randomly distributed. In the present work, the sedimentation process affects both the volume fraction of particles and the orientation of contact angles between particles. The contact angle of a pair of particles is defined as the angle the axis connecting their centroids makes with a reference plane. At first glance, the process is not well-suited for generating isotropic granular systems since wall effects affect the orientation of contact angles so that the medium is anisotropic [31]. However an extraction of a cubic cell from the granular system can be performed as well. Fig. 8 illustrates an example of distribution of contact angles obtained from a system composed of 300 particles without neither extraction nor interconnection. Three reference planes are considered, namely the xy, xz and yz ones. The radial plot describes the probability of occurrence for the contact angles arising in the granular system [20]. 36 classes are considered from 0° through 350°, in increments of 10°. A 0-degree orientation corresponds to a horizontal sphere-sphere contact, and a 90-degree one corresponds to a vertical sphere-sphere contact. In the present instance, both 0 and 90-degrees directions are over-represented. Walls affect the orientation of contact angles so that the resulting system is strongly anisotropic. Some extractions are considered in order to eliminate the wall effects. A first one is carried out for q ¼ 0:05. The corresponding results are presented in Fig. 9. The orientation distribution is almost identical to the previous one. This result was expected since the quantity of extracted material is insufficient to impact the number of contacts. A solution would consist in using a higher value of q but the final low number of spheres would conflict with the definition of the RVE. That is why, we consider an extraction of 300 spheres from a denser packing. Preliminary tests led us to choose an initial packing composed of 1500 particles. Results are depicted in Fig. 10. The

Fig. 7. Influence of the number of particles on the volume fraction for three kinds of extraction.

Fig. 8. Orientation distribution in xy, xz and yz plane, q ¼ 0.

orientation of contact angles turns out to be much more random and less prone to wall effects. Indeed, as exhibited in previous works [31], a sizeable number of spheres is necessary in order to obtain a near-random distribution of contact angles. 4.1.3. Covariance The covariance or two-point probability function [19,32] is a mathematical tool which is widely used in the field of heterogeneous media. Indeed, this one provides interesting information on both the microstructure and the spatial distribution of heterogeneities. In the present context, the covariance is used for detecting periodic and structured configurations as well. From a mathematical point of view, the parameter describes the probability that two random points belong to the same phase. Thus, this one directly depends on both the distance and the orientation of the vector connecting two random points. Typically, calculations are limited to the

Fig. 9. Orientation distribution in xy, xz and yz plane, q ¼ 0:05.

W. Leclerc et al. / Advances in Engineering Software 77 (2014) 1–12

7

Fig. 11. Covariances in x, y and z directions, case of a granular packing of 300 unpenetrable spheres without extraction ð/ ¼ 0:582Þ.

Fig. 10. Orientation distribution in xy, xz and yz plane, extraction of 300 spheres.

3 directions of the Cartesian coordinate system (x, y and z), and the distance is limited to the RVE size. Three values are particular. First, when the two points coincide, the covariance is equal to the volume fraction of spherical particles. Second, the covariance has a local minimum at a distance equal to the characteristic distance of inclusions, namely the sphere radius in the present context. Finally, the covariance is equal to the square of the volume fraction when both the distance between the two points tends to infinity and the heterogeneities are well-dispersed. This last point provides a good criterion for investigating whether the spatial distribution of particles is homogeneous or heterogeneous. Furthermore, in the context of a homogeneous state of dispersion, the covariance has to lead to equivalent results for each investigated direction. Figs. 11 and 12 show the evolution of the covariance according to the distance between two random points. The distance is normalized by the RVE size, and the calculations are carried out in x, y and z directions. In the first configuration, a granular packing of 300 unpenetrable spheres is considered and no extraction is used. In the second configuration, a granular packing of 1500 unpenetrable spheres is generated and an extraction of a cubic cell composed of 300 spheres is carried out. / refers to the real volume fraction of spheres which is respectively 0.582 and 0.635 for the first and the second configuration. First, in both cases, the first peak is noticeable for a normalized distance equal to the sphere radius (r). In addition, as expected, a normalized distance equal to 0 leads to the volume fraction of spherical particles. Second, the curves do not converge to the square of the volume fraction when no extraction process is employed. In addition, important disturbances, typical of poorly dispersed particles, are noticeable. In the second case for which an extraction of 300 spheres is carried out, the three curves converge to the square of the volume fraction. Small disturbances are visible but no structural defect is detected. Thus, the particles are homogeneously distributed. These findings confirm the previous ones about the benefit of an extraction of 300 spheres from a granular system of 1500 particles. This sole configuration is selected for further investigations. 4.2. ETC A numerical investigation of the ETC is set up. Calculations are carried out using the FFT approach described in 3.3, and the properties provided in the same subsection. In a first step, a complete

Fig. 12. Covariances in x, y and z directions, case of an extraction of a cubic cell composed of 300 spheres from a granular packing of 1500 unpenetrable spheres ð/ ¼ 0:635Þ.

sample of granular packings composed of unpenetrable spheres, is generated with a number of spheres ranging from 300 to 1500 particles. For each granular system, three kinds of extraction are considered, namely q ¼ 0, q ¼ 0:05 and an extraction of 300 spheres. Fig. 13 illustrates the predictions. The results are respectively depicted with a continuous black line and black circles, a dashed red line and red triangles, and a dotted green line and green plus symbols. One can observe that the number of spheres only impacts on the ETC when no extraction is done. Besides, in this case, numerical estimates are lower than the ones obtained for an extraction of 5% or 300 spheres. This observation is in agreement with the findings of Fig. 7 about the volume fraction of the granular packing when wall effects are taken into account. Indeed, in the context of a linear heat transfer problem without debonding, the volume fraction of heterogeneities is a predominant parameter. Thus, the lower the volume fraction value is, the lower the ETC becomes. The case of an extraction of 300 spheres leads to similar results to the ones obtained when considering an extraction of 5%. The predictions are moreover independent of the number of spheres even though the homogeneity of the distribution of

8

W. Leclerc et al. / Advances in Engineering Software 77 (2014) 1–12

Fig. 13. Influence of the number of spheres on the ETC for several kinds of extraction (/ < 0:614; /  0:64 and /  0:64).

contact angles is greatly dependent on the volume fraction of particles. These observations lead to the conclusion that the defects highlighted in 4.1.2 have little effect on the ETC. In a second step, interconnection areas are taken into account with a dimensionless contact radius Rc up to 0.25. In such a context, the granular system is composed of partially penetrable spheres and the more the particles overlap, the higher Rc is and the higher the volume fraction of spherical particles is. Thus, a priori interconnection areas should lead to higher ETC. A numerical investigation was set up for checking this point. For this purpose, first a large sample of granular systems composed of 1500 particles was generated, then a cubic RVE composed of 300 spheres was extracted from each one. Thus, the volume fraction of granular systems is close to the random close packing when Rc ¼ 0. Finally, ETC are predicted using the FFT approach described in 3.3.2. Fig. 14 illustrates the influence of the volume fraction of aluminum spheres on the ETC for a dimensionless contact radius Rc ranging from 0 to 0.25. A comparison is done with a third-order analytical estimate [33–35]. The model corresponds to the case of fully penetrable spheres for which overlaps are allowed. We talk about a system composed of soft spheres. In addition, effective conductivities are given for two kinds of structured systems, namely a simple cubic (SC) lattice and a centered cubic (CC) lattice. As for the case of random systems, interconnection areas are taken into account. Thus, Rc is equal to 0, 0.1 and 0.2 for the CC structure so that the volume fraction of particles is respectively 0.68, 0.691 and 0.72, and 0.4 for the simple cubic (SC) structure so that the volume fraction of particles is 0.685. The ETC is depicted with a dashed blue line and blue plus symbols which individually represent one single calculation. Error bars are not provided since these ones are too small to be visible. The third-order estimate is represented with a dash-dotted green line. SC and CC structures are respectively illustrated with black cube and pink triangle symbols. First, one can notice that all numerical results are located below the curve associated with the analytical predictions. This finding was expected since the third-order estimate is an upper bound of the present model for which the spheres are partially overlapped. Second, the ETC linearly increases with increasing volume fraction of spheres, and consequently increasing Rc . Thus, interconnection areas positively impacts on the effective conductivity of the composite. Finally, the structured granular systems, which have anisotropic orientation distributions, lead to similar results to the ones obtained with random granular systems. This confirms the relative

Fig. 14. Influence of the volume fraction of spherical particles on the ETC and comparison with a third-order estimate.

influence of contact angles on the ETC. This point will be much more developed in the next section.

5. Anisotropic granular systems The section is devoted to the investigation of anisotropic granular systems. ETC of lattice and aligned structures are predicted and compared to the results obtained in an isotropic configuration. Finally, a methodology is proposed for directing the heat flux using controlled granular systems.

5.1. Structured granular systems SC and CC structures are periodic crystal systems which are typically represented with a cubic unit cell called a Bravais lattice. The SC system is composed of eight points located on the corners of the cubic pattern. The CC structure is composed of nine points, namely the eight corners of the cube and its center. In both cases, sphere– sphere contacts have specific orientations. These ones can be described by radial plots as seen in 4.1.2 for random granular systems. Figs. 15 and 16 exhibit the radial plots of SC and CC structures. One can notice that the SC system leads to a plus pattern, and the CC system leads to a cross pattern. Each one is invariant by rotation of 90° and each structure is both anisotropic and symmetrical. Furthermore, volume fraction values associated to each system are respectively 0.523 and 0.680. Thus, a CC system is much denser than a SC one. Some predictions of ETC are provided in 4.2. For information purposes, contact radii are adjusted in order to fit the range of investigated volume fraction values. Thus, in the case of a SC system, the effective conductivity is estimated to 71.07 W/m K for a dimensionless contact radius Rc set to 0.4. In the case of a CC system, ETC are respectively 70.74, 71.76 and 75.04 W/m K for Rc set to 0, 0.1 and 0.2, respectively. As exhibited on Fig. 14, the predictions are quite close to the ones obtained for a random granular system. Relative differences with respect to the value given by the isotropic model are less than 0.2%. Consequently, even though SC and CC systems are anisotropic structures, these lead to isotropic ETC and no directional effect arises. However it is important to precise that this conclusion is only true for diffusion problems described by tensors of rank 2. For example, elastic problems for which tensors are of rank 4 lead to other findings [36].

W. Leclerc et al. / Advances in Engineering Software 77 (2014) 1–12

9

Fig. 17. Visual representation of a network of aligned spheres with x-axis preferential direction.

Fig. 15. Orientation distribution in the xy plane, SC structure.

of the total volume. Fig. 18 shows the ETC predicted in x, y and z directions for three preferential directions, namely 0°, 45° and 90° to the x direction. In each case, the results exhibit a strong anisotropy. Thus, at a 0-degree orientation, the effective conductivities are 59.81 W/m K, 50.24 W/m K and 50.47 W/m K for x, y and z directions, respectively. The associated relative gain in x direction with respect to values predicted in the xy plane, is about 19%. At a 90-degree orientation, results are quite similar except that the y-axis is now the preferential direction. In fact, a quasi-symmetry is obtained which highlights the rotational invariance of the present problem. When considering a 45-degree orientation, the effective conductivities are 54.58 W/m K, 54.58 W/m K and 48.02 W/m K for x, y and z directions, respectively. As a result, the predictions exhibit two preferential directions, namely x and y-axes. This finding is in agreement with the symmetric role of x and y-axes in such a configuration.

5.3. Controlled anisotropy

Fig. 16. Orientation distribution in the xy plane, CC structure.

5.2. Alignments of particles Transversally isotropic structures with a single preferential direction are studied. Chains of aligned spheres with interconnection are generated using MultiCAMG. Fig. 17 illustrates an example of structure in the xy plane. Spheres are aligned in a given direction, the x one in the present example, and no contact is permitted between two chains in perpendicular directions. The distance between two chains (e) is a tenth of the diameter (D) of spheres in order to minimize void areas. In this study, several directions are considered in the xy plane but one has to know that the system remains periodic in the z direction. Moreover, the volume fraction of spheres is much lower than the one of a random granular system, namely about 40% without interconnection. That is why, dimensionless contact radii are set to a high value of 40% in order to increase the volume fraction of aluminum particles to about 50%

5.3.1. Setting up A network of aligned spheres leads to a transversally isotropic thermal behavior of the composite. However, such a configuration has two major flaws. First, such a structure is difficult to manufacture. Second, the volume fraction of the particle skeleton is low in comparison with that of a random isotropic granular system. Hence the ETC is below our expectations even for a high dimensionless contact radius. We propose a methodology for generating anisotropic random granular systems for which the two previous drawbacks are avoided. It deals with controlling the anisotropy of random granular systems generated by sedimentation with the help of inner walls. Fig. 19 illustrates the process. Walls are regularly distributed within the container, and directed in the direction of sedimentation, namely the z direction. In the present example, 4 inner walls are used and the box is regularly subdivided into 3-by-3 columns. Thus, contact angles are preferentially directed towards the z direction too, and the composite material is anisotropic. The size of particles is described by the ratio between the column length L and the sphere diameter D. t denotes the wall thickness. Such a kind of structure leads to higher volume fractions of spherical particles than those obtained in the case of chains of aligned spheres. From an experimental point of view, this approach is feasible. The sole difficulty is that the walls must be removed after generating and drying the particle skeleton composed of

10

W. Leclerc et al. / Advances in Engineering Software 77 (2014) 1–12 Table 2 Influence of the number of columns on the ETC in x, y and z directions.

⁄

Configurations

/

kx (W/m K)

ky (W/m K)

kz (W/m K)

1-by-1 2-by-2 3-by-3 3-by-3⁄

0.4618 0.4603 0.4587 0.4625

50.67 51.04 51.41 51.79

50.75 50.75 51.35 51.74

52.86 53.66 53.86 54.03

Periodic boundary conditions.

Periodicity verifies the assumption of infinite medium and consequently leads to more accurate results according to some authors [25]. In the present study, results are in good agreement between the two FE calculations. In view of this numerical investigation, from now on, all calculations will be carried out using the 3-by-3 configuration.

Fig. 18. ETC as function of the orientation of aligned spheres (/ ¼ 0:498).

Fig. 19. Visualization of a structure with controlled anisotropy (3-by-3 subdivision, / ¼ 0:460).

PMMA Particles. Thus, voids left by the inner walls are filled during the immersion in the alumina slurry, and no bridge arises between two columns. 5.3.2. Validation We work with an infinite medium for which the number of spheres is itself infinite. Thus, in such a context, ETC has to be independent of the number of columns. A numerical study is proposed for checking out this point. Calculations are carried out using the FEM presented in 3.3.3. In a first step, L to D and t to D ratios are set to 2.5 and 0.1 respectively, and we carry out calculations for 1-by-1, 2-by-2 and 3-by-3 columns. Results are provided in Table 2. / is the mean volume fraction of particle skeletons. kx ; ky and kz denote the ETC in x, y and z directions, respectively. The volume fraction slightly fluctuates between 0.44 and 0.48 in the set of generated granular systems. That is why, only granular systems with a volume fraction of 0.46 ± 0.004 were selected for achieving consistent comparisons. Results highlight the anisotropic behavior of composites since kz is greater than kx and ky . One can notice that the number of columns slightly impacts on predictions. Thus, kz comes from 52.86 W/m K for a 1-by-1 configuration to 53.86 W/m K for a 3-by-3 one. In the last case, a comparison is realized with finite element (FE) simulations with periodic boundary conditions. This configuration is designated with an asterisk upperscript.

5.3.3. Numerical investigations A numerical study is set up in order to evaluate the ETC according to a set of relevant parameters. First L to D and t to D ratios are investigated for better understanding their impact on both volume fraction and thermal response of anisotropic composites. t to D ratio is chosen so that the wall thickness is minimal and the volume fraction of the particle skeleton is maximized. However care must be taken not to consider a too small value, since contact bridges could arise in the numerical simulation. Our choice is to set t to D ratio to 0.1. A study is then proposed for assessing the impact of the sphere diameter. Several samples of anisotropic granular systems with L to D ratios ranging from 2.5 to 4 are generated considering the 3-by-3 configuration, and ETC of composites are estimated using the FE approach described in 3.3.3. Results are provided in Table 3. The mean volume fraction / turns out to be sensitive to L to D ratio. At first sight, one would expect that the greater the L to D ratio is, the higher is the volume fraction. However, an L to D ratio equal to 4 leads to a lower volume fraction than an L to D ratio equal to 3.5. This result seems to be related to the formation of lattice structures within each column. For a nonwhole value, CC systems appear as observed on the third vertical column of Fig. 19 for which L to D ratio is 2.5. For a whole value, SC lattices constitute the majority of arrangements. This affects the volume fraction of spherical particles since an SC system has a low volume fraction, namely 0.523. Thus, the highest ETC are obtained when considering an L to D ratio equal to 3.5. On the other hand, kx and ky values remain close whatever is the chosen configuration. As expected, kz is greater than kx and ky but the relative gain greatly depends on the L to D ratio. The greater the L to D ratio is, the higher the relative gain is. with a maximum relative gain of 5.48% with respect to the average of kx and ky . In a second step, anisotropic granular systems composed of interconnected particles are generated using the proposed approach. Thus, the dimensionless contact radius Rc is varied from 0 to 0.35 while the volume fraction of granular systems ranges from about 0.56 to about 0.66. L to D and t to D ratios are respectively set to 3.5 and 0.1 in accordance with the previous discussions. Fig. 20 illustrates the estimated ETC in the z direction and the xy plane. kz is depicted with a continuous black line and black

Table 3 Influence of the L to D ratio on the ETC in x, y and z directions. L to D ratio

/

kx (W/m K)

ky (W/m K)

kz (W/m K)

2.5 3 3.5 4

0.4658 0.4898 0.5610 0.5327

51.85 53.52 60.72 58.72

51.82 53.70 60.74 58.77

54.74 57.70 63.96 61.99

W. Leclerc et al. / Advances in Engineering Software 77 (2014) 1–12

11

7. Funding The authors would like to gratefully acknowledge the European Union for the financial support under the INTERREG IV FranceWallonie-Vlaanderen Program PRISTIMAT2 (No. FW 1.1.27), and the région de Picardie for its CASIMAT Program funding (No. Reg. 13001).

Acknowledgment The authors thank LMPCA (UVHC), SIRRIS, and INISMA-CRIBC (members of EMRA) for supplying the experimental data about the manufacturing process.

References

Fig. 20. Influence of the volume fraction on the ETC in the z direction and the xy plane.

circles. The ETC in the xy plane is depicted with a dashed red line and red triangles. For information purposes, this one is measured by taking the average of kx and ky . A comparison is then done with the predictions obtained with random packings using the FFT approach, and provided in 4.2. These results are depicted with a dash-dotted blue line. The main observations are as follows. kz is greater than the ETC in the xy plane which highlights the anisotropic behavior of the composite. Thus, the relative difference with respect to predictions obtained in the xy plane comes from 5.32% for a volume fraction of about 0.561 to 11.75% for a volume fraction of about 0.660. In addition, predictions for a random system are quite close to those in the xy plane with a relative gap less than 1.5%. For example, when the volume fraction is 0.64, the ETC are 67.13 and 74.88 W/m K in the xy plane and the z direction, respectively, while the one obtained for the random system is 66.36 W/m K. As a result, a controlled anisotropy leads to an important gain of ETC in the z direction in comparison with the random configuration, up to about 14% for a volume fraction of 0.66. This benefit is, moreover, obtained without loss of ETC in the xy plane. 6. Conclusions In the present work, an extension of the numerical approach of Ferguen et al. [1] was developed for the statistical analysis of granular systems, and computational purposes. Thus, the wellestablished FFT-based method was implemented using the Eyre– Milton scheme for the fast assessment of ETC of alumina/Al composite materials. Statistical tools as the covariance and the radial plots were introduced and used for better analyzing the connection between microstructure and ETC. In a first step, an isotropic benchmark was set up. This one respects several criteria related to the volume fraction of particles and the homogeneity of the medium. These standards, which are not easily verified due to wall effects, led us to extract a cubic cell from granular systems. In a second step, anisotropic configurations were studied using the developed tools. Whether lattice structures lead to an isotropic thermal response, transversally isotropic materials are of interest for directing the heat flux. A technologically feasible anisotropic structure was proposed for this purpose. The concept takes advantage of wall effects for inducing a preferential orientation of contact angles using inner walls. Results highlight a 14% increase in ETC in comparison with a random system for a volume fraction of 0.66.

[1] Ferguen N, Cogné C, Bellenger E, Guessasma M, Pélegris C. A numerical model for predicting effective thermal conductivities of alumina/Al composites. J Compos Mater 2013;47:3311–21. [2] Descamps M, Richart O, Hardouin P, Hornez JC, Leriche A. Synthesis of macroporous b-tricalcium phosphate with controlled porous architectural. Ceram Int 2008;34:1131–7. [3] Descamps M, Duhoob T, Monchaud F, Lu J, Hardouin P, Hornez JC, et al. Manufacture of macroporous b-tricalcium phosphate bioceramics. J Eur Ceram Soc 2008;28:149–57. [4] Anné G, Vanmeensel K, Vleugels J, et al. Electrophoretic deposition as a novel near net shaping technique for functionally graded biomaterials. Key Eng Mater 2006;314:213–8. [5] Gandra J, Miranda R, Cilaaß P, et al. Functionally graded materials produced by friction stir processing. J Mater Process Technol 2011;211:1659–68. [6] Shin KH, Natu H, Dutta D, et al. A method for the design and fabrication of heterogeneous objects. Mater Des 2003;24:339–53. [7] Mishnaevsky Jr L, Derrien K, Baptiste D. Effect of microstructure of particle reinforced composites on the damage evolution: probabilistic and numerical analysis. Compos Sci Technol 2004;64:1805–18. [8] Cambou B. Behaviour of granular materials. CISM courses and lectures, vol. 385. New-York: Springer; 1998. [9] http://cran.r-project.org/. [10] Moulinec H, Suquet P. A numerical method for computing the overall response of nonlinear composites with complex microstructure. Comput Methods Appl Mech Eng 1998;157:69–94. [11] Michel J, Moulinec H, Suquet P. Effective properties of composite materials with periodic microstructure: a computational approach. Comput Methods Appl Mech Eng 2001;172:109–43. [12] Leclerc W, Karamian P, Vivet A. An efficient stochastic and double-scale model to evaluate the effective elastic properties of 2D overlapping random fibre composites. Comput Mater Sci 2013;69:481–93. [13] Eyre D, Milton G. A fast numerical scheme for computing the response of composites using grid refinement. Eur Phys J Appl Phys 1999;6:41–7. [14] Berryman JG. Random close packing of hard spheres and disks. Phys Rev A 1983;27:1053–61. [15] Jaeger HM, Nagel SR. Physics of granular states. Sciences 1992;255:1523–31. [16] Torquato S, Truskett TM, Debenedetti PG. Is random close packing of spheres well defined? Phys Rev Lett 2000;84:2064–7. [17] Donev A, Cisse I, Sachs D, Variano EA, Stillinger FH, Connelly R, et al. Improving the density of jammed disordered packings using ellipsoids. Sciences 2004;303:990–3. [18] Chiu SN, Stoyan D, Kendall WS, Mecke J. Stochastic geometry and its applications. Chichester: Wiley; 2013. [19] Lochmann K, Oger L, Stoyan D. Statistical analysis of random sphere packings with variable radius distribution. Solid State Sci 2006;8:1397–413. [20] Zhou C, Ooi JY. Numerical investigation of progressive development of granular pile with spherical and non-spherical particles. Mech Mater 2009;41:707–14. [21] Scott GD. Packing of spheres: packing of equal spheres. Nature 1960;188:908–9. [22] Cundall PA, Strack ODL. Discrete numerical model for granular assemblies. Geotechnology 1979;29:47–65. [23] Haddad H, Guessasma M, Fortin J. Heat transfer by conduction using DEM– FEM coupling method. Comput Mater Sci 2014;81:339–47. [24] Drugan WJ, Willis JR. A micromechanics-based nonlocal constitutive equation and estimates of representative volume element size for elastic composites. J Mech Phys Solids 1996;44:497–524. [25] Kanit T, Forest S, Galliet I, Mounoury V, Jeulin D. Determination of the size of the representative volume element for random composites: statistical and numerical approach. Int J Solids Struct 2003;40:3647–79. [26] Moulinec H, Silva F. Comparison of three accelerated FFT-based schemes for computing the mechanical response of composite materials. Int J Numer Methods Eng 2014;97(13):960–85.

12

W. Leclerc et al. / Advances in Engineering Software 77 (2014) 1–12

[27] Willot F, Abdallah B, Pellegrini YP. Fourier-based schemes with modified Green operator for computing the electrical response of heterogeneous media with accurate local fields. Int J Numer Methods Eng 2014;98(7):518–33. [28] Lubachevsky BD, Stillinger FH. Geometric properties of random disk packings. J Stat Phys 1990;60:561–83. [29] Scott GD, Kilgour DM. The density of random close packing of spheres. J Phys D 1969;2:863–6. [30] Finney JL. Fine structure in randomly packed dense clusters of hard spheres. Mater Sci Eng 1976;23:199–205. [31] André D, Iordanoff I, Charles JC, Neauport J. Discrete element method to simulate continuous material by using the cohesive beam model. Comput Methods Appl Mech Eng 2012;213:113–28.

[32] Torquato S. Random heterogeneous materials. New-York: Springer; 2002. [33] Beasley JD, Torquato S. Bounds on the conductivity of a suspension of random impenetrable spheres. J Appl Phys 1986;60:3576–81. [34] Torquato S. Effective stiffness tensor of composite media I. Exact series expansions. J Mech Phys Solids 1997;45(97):1421–48. [35] Torquato S. Effective stiffness tensor of composite media II. Application to isotropic dispersions. J Mech Phys Solids 1997;46(98):1411–40. [36] Nye JL. Physical properties of crystals their representation by tensors and matrices. Oxford: Clarendon; 1957.