© ESO, 2016

1. Introduction

Meteorites preserve in their structure and composition rather detailed information on the processes that were acting during planetary formation and, in particular, on processes active in their respective parent bodies. It is possible to recover this information by modelling the structure, composition, and thermal history of meteorite parent bodies and by comparing this modelling with results of laboratory investigations of the composition and structure of meteorites, and of their thermal evolution as inferred from closure temperatures and ages of radioactive decay systems from cooling rates. The basis of such investigations are models for the thermal evolution of their parent bodies.

The parent bodies of the ordinary chondrites are suspected to be undifferentiated planetesimals with diameters ranging between at least 100 km and at most a few hundred kilometers. There also exists the possibility that part of the parent bodies of chondrites is differentiated in their interior region and only bear an outer crust of chondritic material; see Elkins-Tanton et al. (2011). Irrespective of such possible structural differences the thermal evolution of these bodies is characterised by a brief initial heating period due to decay of short-lived radioactive nuclei (mainly ²⁶Al), lasting for a period of a few million years. This heating period is followed by a long-lasting cooling period of at least 100 Ma duration, where the heat liberated within the body is transported to its surface where it is radiated away. A realistic modelling of this rise and fall of temperature requires the use of a value of the heat conduction coefficient, K, that is as accurate as possible to solve the heat conduction equation.

The present paper is the fourth in a series in which we attempt to improve the model construction for the thermal evolution of planetesimals (Henke et al. 2012a, 2013; Gail et al. 2015). Here we aim to improve the modelling of the heat conductivity, K, of the material. A determination of K is complicated because of the rather complex structure of the chondrite material. First, it is a mixture of many different materials with very different individual properties and, second, it is of granular structure that consists of a two-component porous mixture of chondrules of 0.1 to 1 mm size and a very fine-grained matrix material with particle sizes of the order of micrometers. The proportions of chondrules and matrix strongly vary between different meteorite classes. This granular material is found in meteorites of the petrologic classes 1 to 3 and is assumed to represent the initial structure of the material from which the parent bodies formed.

During the course of the thermal evolution of the parent body the initially porous granular material is compacted once the temperature is raised to sufficiently high values such that creep processes in the solids are activated. This results in a gradual deformation of the granules and finally in a closure of all pores. This process is witnessed by the meteorites of petrologic classes 4 to 6 (or 7), which represent different preserved stages of the compaction process at increasingly higher temperature.

In a completely compacted state, the material still is a heterogeneous mixture of several mineral compounds and iron and/or iron sulphide. The next thermal evolutionary step is the onset of melting. Meteorites where the material has once been molten, at least partially, are achondrites or primitive achondrites. For the ordinary chondrite classes, no achondrites seem to exist that could clearly be related through compositional similarities to one of the parent bodies of the ordinary chondrite classes. Either the parent bodies of the ordinary chondrites never were heated to such an extent in their interior that their central region was molten during part of their thermal evolution, or they never suffered very strong collisions with other bodies by which material from their interior region is excavated. Generally it is assumed that the first alternative holds, but this is not really proven.

To study the thermal evolution of parent bodies of ordinary chondrites, one has to determine the heat conductivity of the granular material from its likely initial stage as a loosely packed powder or sand-like material through its various compaction stages to the final completely compacted stage of a strongly heterogeneous mixture of several components. If melting also occurs, one also needs the heat conductivity of the melt and of the material after its solidification, but this is a much simpler problem than that before melting; we do not consider bodies for which this is relevant.

One way to determine the heat conductivity is to measure the heat conductivity coefficient for a number of chondrites and temperatures, as carried out by Yomogida & Matsui (1983), Opeil et al. (2010, 2012). Then one can use such results in model calculations. This is how K is determined in most published model calculations of the thermal evolution of meteoritic parent bodies, including our own. The shortcoming of this approach is that K depends on a number of parameters that change during the evolution of the bodies and that not all stages of evolution of all the bodies are represented in meteorite collections. Furthermore, a lot of cracks formed in mineral components via shock waves, for example from the impact that excavated the meteorite. Such cracks reduce the heat conductivity of the material in an unpredictable way, such that measured heat conductivities of meteoritic material are not necessarily representative for the same material during the early evolution of the body. This makes it desirable to develop an approach that allows us to calculate K from the properties of the individual components from which the material is composed.

Another approach is to calculate the heat conductivity from properties of the mixture components via approximate mixing rules, which have been developed for technical applications (see, e.g., the review Berryman 1995). This approach has many advantages because this often combines a simple calculation method with acceptable accuracy. We compare such methods with our more detailed calculations of the heat conductivity.

In this paper we aim to calculate K for composite materials encountered in the parent bodies of ordinary chondrites by numerically constructing samples of such composite materials and solving the equation of heat transfer for the boundary value problem of a sample held at different fixed temperatures on opposite sides of a sample. From the calculated stationary heat flux through the sample one readily determines an effective heat conductivity of the composed material. This allows us to determine a value of K for given composition and structure of the composite material. Although this procedure is numerically demanding, for applications it is sufficient to generate a table for the dependence of K on temperature and porosity for the mixture of interest; from this table K can be determined by interpolation during the course of model calculations.

The heat conductivity is determined by three processes: the transport via phonons in electrically insulating minerals, by conduction electrons in the metallic iron component of the mixture, and by radiative transfer through the transparent minerals and voids of the mixture. These have different dependencies on temperature and their contribution to the total heat conduction coefficient is of a different order of magnitude. At low temperatures, conduction by the solid material dominates and the contribution by radiative transfer is negligible. At high temperatures (T> 1000 K), the radiative contribution steeply increases and finally becomes the dominating mode of heat transfer. Both modes of heat transfer, by phonons and conduction electrons, on the one hand, and by radiation, on the other hand, are independent of each other and can simply be superposed.

The problem of heat conduction through minerals by radiative transfer was already studied for modelling the structure of planets. Here one can appeal to studies in planetary physics (e.g. Hofmeister 1999, 2014) and we do not elaborate on this part of the problem. We attempt to calculate the heat conduction via solid-state heat conduction for a material that is an intimate mixture of different mineral and iron grains and, before compaction, of voids.

The voids in the granular material are filled with gas from the accretion disk and possibly with outgassed volatiles. This also contributes to the heat conductivity, but in earlier test calculations for the models in Henke et al. (2012a) we found this to be negligible compared to the other processes. However, this may be different for carbonaceous chondrites (not considered here) during the phase of ice melting and evaporation.

The main dependencies of interest for our problem are the temperature dependence of K, its dependence on the composition and granular structure of the material, and how this is determined from the properties of the different components. A possible pressure dependence is of minor importance because central pressures in the bodies that we have in mind do not much exceed 1 kbar or stay well below that value. Effects resulting from the compressibility of the material remain negligible in that case. We do not consider these effects.

The plan of this paper is as follows: in Sect. 2 we discuss the composition of chondrites, in Sect. 3 we describe the calculation of heat conductivity of composite materials, and in Sect. 4 we describe the calculation of heat conductivity for granular materials. In Sect. 5 we apply our results to the parent body of H chondrites. The paper closes with some final remarks in Sect. 6. The appendices contain descriptions of our calculational methods for solving the heat conduction equation, generating random multicomponent mixtures, and generating the packing of spheres.

2. Composition of chondritic material

2.1. Structure of chondrites

The refractory solids of the precursor material from which asteroids form is composed of two morphologically different components: matrix and chondrules. The matrix is a very fine grained material consisting of particles with diameters ≲1 μm, while chondrules are glassy rounded beads (sometimes broken) with diameters usually in the range between 0.1 mm and 1 mm. In meteorites, these two components form a granular material, and a significant fraction of the total space filled by this binary mixture are empty pores between the particles. In the precursor material of carbonaceous chondrites part or all of this pore space was once filled with ice.

Table 1 shows some characteristic data for the different known chondrite groups. The relative abundances of matrix material and chondrules is very different between different groups, ranging for carbonaceous chondrites from almost pure matrix in CI chondrites to a slightly higher volume fraction of chondrules than matrix material in CO chondrites, while for ordinary chondrites the volume fraction of matrix is small to very small, except perhaps for the rare K grouplet.

The properties of such binary matrix-chondrule mixtures were discussed in Gail et al. (2015). They can be characterised by the ratio, f_ma, of matrix volume fraction, V_ma, to total volume fraction of matrix and chondrules (without the pore space) ${\mathit{f}}_{\mathrm{ma}}\mathrm{=}\frac{{\mathit{V}}_{\mathrm{ma}}}{{\mathit{V}}_{\mathrm{ma}}\mathrm{+}{\mathit{V}}_{\mathrm{ch}}}\mathrm{·}$(1)On average this quantity (see Table 1) ranges for ordinary chondrites between f_ma ~ 0.1 and f_ma ~ 0.15 such that the matrix-chondrule mixture in ordinary chondrites of type H, L, LL, EH, and EL is chondrule dominated as defined in Gail et al. (2015), where matrix fills part of the space not occupied by chondrules. In carbonaceous chondrites this value varies on average between f_ma ~ 0.42 and f_ma = 1, such that the matrix-chondrule mixture in carbonaceous chondrites of type CI, CM, and CK is matrix-dominated, where chondrules are interspersed within a matrix ground material. The CR, CO, and CV chondrites are intermediate cases.

The chondritic material is porous (see Consolmagno et al. 2008,for a detailed review). The volume fraction of voids, i.e. the porosity, is given by $\mathit{\phi }\mathrm{=}\frac{{\mathit{V}}_{\mathrm{vo}}}{{\mathit{V}}_{\mathrm{ma}}\mathrm{+}{\mathit{V}}_{\mathrm{ch}}\mathrm{+}{\mathit{V}}_{\mathrm{vo}}}\mathit{,}$(2)where V_vo is the volume of void space, which varies in meteorites between a few percent and up to about 30 percent (Consolmagno et al. 2008). Typical values are shown in Table 1; for individual meteorites of a group, there is significant scatter of φ. These numbers have to be considered with some caution for two reasons. First, shock waves experienced in the course of asteroid collisions, including those ejecting a meteorite from its parent body, produce cracks in the material (cf. Consolmagno et al. 2008). The porosity measured for meteorite specimens is in part impact generated. In particular, low values of φ can be suspected to be of impact origin. Second, an impact into porous material causes a partial compaction of the target material (Beitz et al. 2014). The excavated material may show a higher degree of compaction than before the impact. For these reasons, the measured porosity of meteoritic material cannot simply be taken as an indigenous property of the parent asteroid material.

The existence of pores in a material obviously has considerable influence on the heat conductivity. For modelling the thermal evolution of a planetesimal, one has to reconstruct the initial degree of porosity of the material from which the parent body was formed because one cannot simply take the porosity measured for meteorites of low petrologic type as the initial state. If one disregards the possible existence of ice in the initial state of the parent bodies of carbonaceous chondrites, the initial material is composed of a mixture of chondrules and dust. The diameters of the chondrules vary in principle over a wide range, but most of them are from a rather narrow range of diameters with size ratios typically ≲2 (e.g. Rubin & Keil 1984). The same most likely holds for the dust particles (Wozniakiewicz et al. 2012, 2013). The random densest packing of granular material of two components with size ratios of less than 2 deviates only very little from that of a mono-disperse granular medium (Fiske et al. 1994). We ignore such minor effects and assume for simplicity that the granular units of the chondrule component all have the same diameter d_ch, and also that the granular units of the dust component all have the same diameter d_du.

This two-component mixture is discussed in the appendix of Gail et al. (2015). It is shown that for such a two-component granular medium the porosity in the chondrule dominated mixture is $\mathit{\phi }\mathrm{=}\mathrm{1}\mathrm{-}\frac{\mathrm{1}\mathrm{-}{\mathit{\phi }}_{\mathrm{ch}}}{\mathrm{1}\mathrm{-}{\mathit{f}}_{\mathrm{ma}}}\mathrm{·}$(3)Here φ_ch is the porosity of the chondrule component if the matrix were absent. In the matrix-dominated case it is $\mathit{\phi }\mathrm{=}\mathrm{1}\mathrm{-}\frac{\mathrm{1}\mathrm{-}{\mathit{\phi }}_{\mathrm{ma}}}{\mathrm{1}\mathrm{-}{\mathit{\phi }}_{\mathrm{ma}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{f}}_{\mathrm{ma}}\mathrm{)}}$(4)where φ_ma is the porosity of the matrix component. For a granular medium of equal-sized grains there are two critical packings with porosities φ ≈ 0.56 and φ ≈ 0.64. The lower value corresponds to the random loose packing, which is just stable under weak external pressure; the higher value corresponds to the random densest packing of equal-sized spheres (cf. Jaeger & Nagel 1992). Since during the growth process of planetesimals the material is joggled again and again, it is likely that the porosities φ_ch in the chondrule-dominated case and φ_ma in the matrix-dominated case equal the random closest packing that is usually achieved in a granular medium after vigorous tapping and joggling. Table 1 shows the porosities φ₀ of the granular mixture predicted by such considerations. They are estimates of the likely initial porosity of the material in planetesimals before the onset of sintering.

For carbonaceous chondrites, except for CK, the highest observed porosities are close to the theoretical values. For ordinary chondrites the predicted value is much higher. This may have a number of reasons, for example sintering owing to higher peak temperatures or compaction by shocks. The decrease of highest measured porosity of ordinary chondrites with increasing shock stage (Fig. 3 of Consolmagno et al. 2008) may be a hint that compaction by impact shock wave played a significant role for observed porosities. For the lowest shock stage 1, the highest observed porosities agree with our model. We take the two-component model as initial structure of chondritic material (corresponding to petrologic type 3).

For the parent bodies of ordinary chondrites, the initially porous chondritic material is compacted by sintering if the body is sufficiently heated by radioactive decays. The pore space vanishes by sintering and the final structure of the material after compaction corresponds to what is found in meteorites of petrologic type 6. The structure of the compacted material is still complicated because the initial chondrule and matrix material is already a heterogeneous mixture of different minerals and iron. In order to determine the heat conductivity of the material we, thus, have to study how, for the initially porous and heterogeneous material, the heat conductivity varies from the initial state with decreasing porosity to the final fully compacted but still heterogeneously composed material. This is accomplished by constructing a numerical model for a granular material with varying compositions of the granular units and by solving numerically the heat conduction equation for this model.

For the carbonaceous chondrites, one has the additional complication that the parent body initially contains ice which, once molten, reacts with the silicates to form phyllosilicates. The corresponding compositional changes (ice bearing → water bearing → phyllosilicates) result in changes in the heat conductivity, which have to be included in model computations. Here we concentrate on ordinary chondrites.

2.2. Mineral components

The composition of the mineral mixture and the mass fractions, X_min, or volume fractions, f_i, of the components varies between the parent bodies of chondrites. The mass-fractions X_min are essential parameters that determine the properties of the material. Because of the marked hierarchy of abundances in the cosmic element mixture, only a few elements can form abundant components in the mixture and a rather limited suite of minerals needs to be considered to determine the gross properties of the material. We only consider the more important components, which contain almost the whole inventory of the abundant rock-forming refractory elements (Si, Mg, Fe, Al, Ca, Ni, Na, K). The pure solid substances and their solid solutions that typically are found in a mineral mixture at temperatures below the melting point of the Fe-FeS eutectic (T = 1 261 K, where differentiation possibly starts; see Fei et al. 1997) and in the absence of water are given in Table 2. The mixture of minerals in ordinary chondrites is dominated by these components.

In the applications we consider in particular the mixture of minerals and of the iron-alloy found in H and L chondrites, which are taken to be equal to the modal mineralogy of H and L chondrites as given in Yomogida & Matsui (1983). This is listed in Table 3 together with X_min, the mass-fractions of the different components in the mixture, which are also taken from Yomogida & Matsui (1983).

The density of the solid solutions given in the table is simply calculated as average of the densities of the pure components, weighted with their mole fractions in the solid solution. This assumes that non-ideal mixing effects are small and do not need to be considered for our calculations.The average density ϱ_bulk of the mixture is calculated from ${\mathit{\varrho }}_{\mathrm{bulk}}\mathrm{=}{\left(\sum _{\mathit{j}}\frac{{\mathit{X}}_{\min \mathit{,j}}}{{\mathit{\varrho }}_{\mathit{j}}}\right)}^{-1}\mathrm{·}$(5)This is the density of the consolidated meteoritic material, i.e. without pores. The volume fractions f_min,i of the components are ${\mathit{f}}_{\min \mathit{,i}}\mathrm{=}\frac{{\mathit{X}}_{\min \mathit{,i}}}{{\mathit{\varrho }}_{\mathit{i}}}\hspace{0.17em}{\mathit{\varrho }}_{\mathrm{bulk}}\mathit{.}$(6)

2.3. Thermal conductivity

This work aims to calculate the heat conductivity K of the mixture of several minerals, metal-alloy, and pores in chondrites from the heat conductivities of the constituents of the mixture. Most of the constituents are not pure substances but are solid solutions with varying concentrations of the solution components in different meteorite classes and with varying heat conductivities depending on the concentrations. Table 2 lists the heat conductivities of the pure components at room temperature (300 K) for the pure substances forming the mineral mixture of ordinary chondrites. We refrain for the moment from the temperature dependence of K and consider only the concentration dependencies. The composition of the solid solutions relevant for ordinary chondrites are listed in Table 3.

Figure 1 shows the heat conductivity of the solid solutions of interest. A detailed explanation of each solid solution as follows:

Nickel-iron alloy:

the data are from Ho et al. (1978) where the heat conductivity of the nickel-iron alloy is tabulated for a large number of concentrations (mole fractions) between x = 0.0048 and x = 0.9947 of the Ni component and for a wide range of temperatures between 4 K and 1100 K. The figure shows the concentration variation of K at room temperature (300 K). The experimental values can be fitted in the range x< 0.34 by the following expression: $\begin{array}{ccc}& & {\mathit{k}}_{\mathrm{1}}\mathrm{=}\mathrm{67.863}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\mathrm{17.232}\hspace{0.17em}\mathit{x}}\\ & & {\mathit{k}}_{\mathrm{2}}\mathrm{=}\mathrm{33.542}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\mathrm{2.0312}\hspace{0.17em}\mathit{x}}\\ & & {\mathit{k}}_{\mathrm{3}}\mathrm{=}\mathrm{361.44}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\mathrm{10.612}\hspace{0.17em}\mathit{x}}\\ & & \mathit{K}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}{\left\{{\mathit{k}}_{\mathrm{1}}^{\mathrm{4}}\mathrm{+}{\mathrm{(}{\mathit{k}}_{\mathrm{2}}^{-4}\mathrm{+}{{\mathit{k}}_{\mathrm{3}}^{-4}}^{\mathrm{)}}}^{-1}\right\}}^{\mathrm{1}\mathit{/}\mathrm{4}}\mathit{,}\end{array}$(7)where x is the mole fraction of nickel in the solid solution. Table 3 shows the interpolated value of K for the alloy in H and L chondrites.

The concentration range of interest for nickel in the nickel-iron alloy in ordinary chondrites is around x = 0.1, see Table 3. The accuracy of the data for K in this range is ± 15% according to Ho et al. (1978). Additionally the dependence of K on x is rather strong such that variations of the Ni fraction in individual iron grains in chondrites may result in noticeable variations of K. Since we do not consider such variations in our calculations of the effective heat conductivity of the mixed material, the value of the heat conductivity of the Ni, Fe alloy in chondrites has some uncertainty.

Olivine:

this is the solid solution of forsterite and fayalite. The data for heat conductivity of the solid solution at a number of fayalite concentrations, shown in Fig. 1, are from Horai & Baldridge (1972). In that paper, values are given for a number of composition intervals between the pure end members. The dots in Fig. 1 correspond to the midpoints of the intervals. The data are for room temperature. The experimental values can be fitted by the following expression: $\mathit{K}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\mathrm{5.3548}\mathrm{-}\mathrm{5.7704}\hspace{0.17em}\mathit{x}\mathrm{+}\mathrm{3.6216}\hspace{0.17em}{\mathit{x}}^{\mathrm{2}}\mathit{,}$(8)where x is the mole fraction of fayalite in the solid solution. Table 3 shows the interpolated value of K for the olivine component in H and L chondrites.

Orthopyroxene:

this is the solid solution of enstatite and ferrosilite. Data for the concentration dependence of the heat conductivity shown in Fig. 1 are again from Horai & Baldridge (1972), where values for a number of composition intervals between the pure end members are given. Only the values for low concentrations of ferrosilite (x ≤ 0.3) are experimental data; they are estimates for higher concentrations, but fortunately only data for low Fs concentrations are needed for chondritic material. The data from Horai & Baldridge (1972) can be fitted by the following expression: $\mathit{K}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\mathrm{4.9094}\mathrm{-}\mathrm{5.0649}\hspace{0.17em}\mathit{x}\mathrm{+}\mathrm{3.5078}\hspace{0.17em}{\mathit{x}}^{\mathrm{2}}\mathit{,}$(9)where x is the mole fraction of ferrosilite in the solid solution. Table 3 shows the interpolated value of K for the orthopyroxene component in H and L chondrites.

Plagioclase:

this is the solid solution of albite, anorthite, and a small amount of orthoclase. Data for heat conductivity of the solid solution of anorthite and albite for different concentration values are given by Horai & Baldridge (1972). No data for the ternary mixture could be found. Since the heat conductivity of pure orthoclase is very similar to that of albite, and the orthoclase concentration is small, we treat the orthoclase component as albite and calculate the heat conductivity for the modified binary mixture. The data from Horai & Baldridge (1972) for the albite-anorthite solid solution can be fitted by the following expression: $\mathit{K}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\mathrm{2.2673}\mathrm{-}\mathrm{3.1000}\hspace{0.17em}\mathit{x}\mathrm{+}\mathrm{2.7768}\hspace{0.17em}{\mathit{x}}^{\mathrm{2}}\mathit{,}$(10)where x is the mole fraction of anorthite in the solid solution. Table 3 shows the interpolated value of K for the plagioclase component in H and L chondrites.

Clinopyroxene:

the composition of the clinopyroxene component in chondrites corresponds to diopside. Table 3 shows the value of K for diopside from Clauser & Huenges (1995).

Troilite:

the value given in Table 3 is the value for pyrrhotite (Fe_1−xS) from Clauser & Huenges (1995). A somewhat lower value of 3.53 ± 0.05 is given in Diment & Pratt (1988). Because troilite is a rare component in the mixture, the modification of the effective heat conductivity using the lower value of K for troilite is negligible.

3. Heat conductivity of composite materials

3.1. Basic model

In principle, it is possible to solve for a composed material, such as found in meteorites, the heat conduction equation, $\mathit{\varrho }{\mathit{c}}_{\mathit{v}}\frac{\mathit{\partial }\hspace{0.17em}\mathit{T}}{\mathit{\partial }\hspace{0.17em}\mathit{t}}\mathrm{=}\mathrm{-}\mathrm{\nabla }\left(\mathit{K}\mathrm{\nabla }\mathit{T}\right)\mathit{,}$(11)to determine the temperature field inside a specimen of that material. Here, c_v is the specific heat capacity of the material, ϱ its mass density, and K the heat conductivity. The spatial variation of the functions ϱ, c_v, and K is very complicated, corresponding to the distribution of different mineralic and metallic particles and voids. To solve the heat conduction equation, in this case one has to decompose the volume for which the temperature field is to be determined into subvolumes corresponding to the individual particles of the material, and to solve the heat conduction equation inside each particle with the appropriate boundary conditions at the contact areas or at the border to adjacent voids. The corresponding solution shows significant local variations in temperature and heat fluxes, in particular, if heat conductivities of the mixture components are significantly different. Such local fluctuations on length scales of the individual particle sizes are without interest and it is appropriate to apply a homogenisation procedure where one only considers averages performed over regions such that, for all relevant mixture components, a big number of the individual particles are contained in that region.

For this purpose we consider for the following an infinitely extended sheet of material of thickness L filled with a multi-component medium and a cube of side-length L within the sheet in which two of its surfaces coincide with the surfaces of the sheet (cf. Fig. 2). The two surfaces of the sheet are held at fixed but slightly different temperatures, T₁ and T₂. If two opposite sides of such a body are held at different temperatures, one has a net heat flux from the hot to cold side. If we solve the heat conduction equation for this cube, the heat fluxes through the individual particles of the mixture and their temperature are found to be slightly different because of different heat conductivities of different components of the material. It can be expected that if the local heat fluxes are averaged over areas with extensions much bigger than the sizes of the biggest particles, there results a net heat flux vector that has constant direction and modulus, independent of the volume over which the averaging was performed. Also, if one averages the temperature over areas that are much more extended than the size of the largest particles and are perpendicular to the average heat flux vector, one finds a constant average temperature over length scales that are much bigger than the particles and independent of the particular size of the area.

Therefore instead of calculating the heat flux through a composed medium by considering each single particle and the corresponding locally fluctuating heat flux and temperature, we only consider the (locally) averaged temperature distribution and (locally) averaged heat flux, and an effective heat conductivity of the mixture. This effective heat conductivity is defined as follows: We solve the heat conduction equation for our cube for fixed temperatures at the two front faces and with appropriate boundary conditions at the four cube faces inside the sheet, which provides us with the total heat flux through the front faces in the stationary state. Then, following the standard method to define a heat conductivity, the average heat conductivity K_eff is defined by ${\mathit{K}}_{\mathrm{eff}}\mathrm{=}\frac{\mathit{L}\hspace{0.17em}\mathit{F}}{{\mathit{T}}_{\mathrm{1}}\mathrm{-}{\mathit{T}}_{\mathrm{2}}}\mathit{,}$(12)where F is the heat flux per unit area through the front faces and L the length of the cube edge.

The method for how the heat conduction equation is solved numerically is described in Appendix A.

3.2. Two-component mixture of forsterite and nickel-iron

As one important application the heat conductivity of a non-porous mineral mixture is considered. This kind of structure corresponds to the structure of the compacted material found inside an asteroid for most of its volume and for most of the time during its thermal evolution.

The first task for determining an effective heat conductivity of a composite material is to generate a distribution of the different components of the mixture within the cube for which we then solve the heat conduction equation in a next step. The method used to generate a random distribution of components for given volume fractions is described in Appendix B.

As a case study, we consider a mixture of olivine and nickel-iron. We chose these two components because they are two major components of chondritic material and because their thermal conductivity is considerably different from each other such that the effective conductivity of the mixture is significantly different from that of the components.

A two-component mixture, with volume fractions of the components varying between 0 and 1 in steps of 0.05, is generated by the procedure described in Appendix B. For the heat conductivity K of the components its value at 300 K is used. The temperature at one of the front surfaces of the cube was set to T₂ = 310 K and the temperature at the opposite front surface was set to T₁ = 300 K. A possible variation of K with temperature was neglected in this model calculation; the temperature dependence is considered elsewhere. The cube was decomposed into cells of different sizes by choosing n = 100,200,300. The size of the balls for generating the random mixture is chosen as 1/20 of the cube length L. With this ball size, the prescribed volume fraction of the mixture components is met by our algorithm with an accuracy better then 6 × 10^-4 in all cases. The heat conduction equation was solved for all cases until a stationary temperature distribution has evolved. For this solution, the heat flux through the cube is calculated and, via Eq. (12), the effective heat conductivity K_eff is determined.

The results are shown in Fig. 3 and in Table 4 for resolutions of n = 100 and n = 200. As is evident from the figure and tabular values, the differences between the effective heat conductivities, K_eff, determined for the solutions with n = 100 and n = 200 is less than one percent and the difference between n = 300 and n = 200 of the order of 2%0. This difference is less than the accuracy with which the heat conductivities K of the mixture components are known. This suggests that it suffices to use a resolution of n = 100 for calculations of the effective heat conductivity for chondritic material, as it is needed for model calculations.

Some scatter of the results may be introduced with different sequences of random numbers for constructing the spatial distribution of the mixture components. Table 5 shows as an example for one case the numerical values of K_eff obtained for ten different realisations of the sequences of random numbers using different seed numbers for the random number generator. The average value of the heat conductivity is K_eff = 13.67 ± 0.02 W/mK and K_eff = 13.73 ± 0.019 W/mK for a resolution of n = 100 and n = 200, respectively. The typical relative deviation of the calculated values of a single calculation for the heat conductivity to the average value of many different calculations, using different seed numbers for the random number generator, is ~0.1% for both resolutions.

The differences resulting for different realisations of the distribution of the components with prescribed volume fractions are small such that there is no need to determine averages for a large collection of realisations. Just a few suffice.

3.3. Analytical approximations for heterogeneous media

For technical applications a number of semi-empirical mixing rules have been developed to calculate the heat conductivity of a composite material from the properties of its components. An overview of the problem can be found, for example in Berryman (1995). In the following, we compare some of these approximations with the results of our numerical calculations to check their applicability to thermal models of asteroids.

One simple approximation is the geometric mean (see Clauser & Huenges 1995) $\begin{array}{ccc}{\mathit{K}}_{\mathrm{eff}\mathit{,}\mathrm{g}}\mathrm{=}\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathit{N}}{\mathrm{􏽙}}}\hspace{0.17em}({\mathit{K}}_{\mathit{i}}{\mathrm{)}}^{{\mathit{f}}_{\mathit{i}}}\mathit{,}& & \end{array}$(13)where f_i are the volume fractions of the components and K_i the individual heat conductivities of the components. The geometric mean is frequently applied in numerical model calculations to determine the heat conductivity of a mixture. Figure 3 shows the effective heat conductivity calculated from this approximation for the olivine-iron mixture as a short-dashed line. The deviation between this and the numerically calculated effective heat conductivity of the mixture is significant. For the mixture considered here the geometric mean is only an approximation of low accuracy that should be avoided.

We define, following Berryman (1995), the function $\mathrm{\Sigma }\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}{\left(\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}\frac{{\mathit{f}}_{\mathit{i}}}{\mathrm{2}\mathit{z}\mathrm{+}{\mathit{K}}_{\mathit{i}}}\right)}^{-1}\mathrm{-}\mathrm{2}\mathit{z}\mathit{.}$(14)An upper and lower bound for the heat conductivity of the mixture are given by (cf. Berryman 1995) $\begin{array}{ccc}{\mathit{K}}_{\mathrm{eff}\mathit{,}\max }& \mathrm{=}& \\ {\mathit{K}}_{\mathrm{eff}\mathit{,}\min }& \mathrm{=}& \end{array}$where K_min and K_max are the lowest and highest values of the heat conductivities of the mixture components, respectively. Figure 3 shows both limits to K_eff as long-dashed lines. The computed values of K_eff are between these limits.

Another approach are effective medium theories (cf. Berryman 1995) resulting in certain mixing rules. A frequently used rule is the Bruggeman mixing rule (Bruggeman 1935), which follows from Eq. (14) by letting Σ(z) = K_eff,B be the effective heat conductivity. This defines the effective heat conductivity K_eff,B as the solution of the equation ${\mathit{K}}_{\mathrm{eff}\mathit{,}\mathrm{B}}\mathrm{=}{\left(\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{N}}\frac{\mathrm{3}{\mathit{f}}_{\mathit{i}}}{\mathrm{2}{\mathit{K}}_{\mathrm{eff}\mathit{,}\mathrm{B}}\mathrm{+}{\mathit{K}}_{\mathit{i}}}\right)}^{-1}\mathrm{·}$(17)In Fig. 3 the effective heat conductivity calculated from this equation is shown as a solid line. For our problem, the effective heat conductivity predicted from the Bruggeman mixing rule is rather close to what is found from numerically solving the heat conduction equation. An inspection of Table 4 confirms this. This suggests that for compacted chondritic material the effective heat conductivity can be determined with reasonable accuracy from Eq. (17).

3.4. Chondritic mixtures

Next the effective heat conductivity is calculated for the six-component mixtures of ordinary chondrites as defined in Table 3. The numerical calculation is carried out for resolutions of n = 100 and n = 200 and for ten different seed numbers for the random number generator to obtain different realisations of the mixture with prescribed average composition. The results are shown in Table 6 (and also in Table 3).

The results for vanishing porosity φ = 0 correspond to a completely compacted material, for example as found in chondrites of petrologic type 6. This cannot be directly compared to laboratory determined values for the heat conductivity of chondrites since the material of chondrites is always somewhat porous because these objects were subject to impact-generated shocks, which generate fractures in the material even if it was initially completely compacted. Therefore we have to extrapolate the measured data to zero porosity. Despite strong scatter, the data of Yomogida & Matsui (1983) show a clear systematic trend in the heat conductivity with porosity (cf. Fig. 4), which was used in Henke et al. (2012a) to fit the data for H and L chondrites with an analytic fit formula. Unfortunately, the number of data points is too small and their scatter too high as to obtain separate fits for H and L chondrites. The average heat conductivity was fitted by $\mathit{K}\mathrm{=}\mathrm{4.3}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\mathit{\phi }\mathit{/}\mathrm{0.08}}\hspace{0.17em}\mathrm{W}\hspace{0.17em}{\mathrm{m}}^{-1}{\mathrm{K}}^{-1}\mathit{.}$(18)The average heat conductivity extrapolated to zero porosity is K₀ = 4.3 W m^-1 K^-1. The uncertainty of the least-squares fit for the value of K₀ is ± 0.67 W m^-1 K^-1. In view of the paucity and strong scatter of available data, the agreement between the extrapolated value of measured heat conductivities to zero porosity and the results of the theoretical calculation of K = 4.89 W m^-1 K^-1 for H chondrites and K = 4.19 W m^-1 K^-1 for L chondrites can be considered as satisfactory. Therefore it is possible to reliably predict the average heat conductivity of the mineral mixture of chondrites by our numerical model.

The numerical results are also compared in Table 6 with results of the effective heat conductivities calculated by the Bruggeman mixing rule, Eq. (17). The result of the Bruggeman mixing rule is close to the numerically calculated value. Since the accuracy of the input data for the heat conductivity of the pure components is low (cf. the tables in Clauser & Huenges 1995), the accuracy of the mixing rule probably suffices in practical applications to calculate the average heat conductivity of a multi-component medium.

4. Heat conductivity of porous material

In the following, we study the dependence of the heat conductivity of a granular material on the porosity.

4.1. Porous material with randomly distributed pores

A porosity of the chondritic material is simulated by introducing a new component in the mixture with a heat conductivity coefficient of K = 0.01. This value is sufficiently small such that the volume elements filled with this fictitious component can be considered as effectively insulating in comparison with the more than a factor of 100 times higher conductivity of the other components. On the other hand, the jump in the value of K between adjacent cells of different conductivity is sufficiently small that the finite difference solution method still works and no special treatment at internal interfaces to vacuum is required.

The first problem we consider is a porous material generated by the method described above. This results, for low porosities, in voids interspersed in an otherwise compact material. This can be assumed to describe the situation encountered in a sintered material shortly before and after closure of the pore space.

Table 6 shows the results of our calculations for the average heat conductivity at different porosities. For comparison, the results obtained by the Bruggeman mixing rule are also shown. The Bruggeman mixing rule yields results that agree reasonable well with the numerical results and may also be used for porous material.

We compare in Fig. 4 the results of our numerical model with measurements of the heat conductivity of sandstone (data from Desai et al. 1974). Sandstone is chosen because this material is formed by compaction of a granular material. In this respect, one could expect a close similarity between our model and sintered chondritic material with sandstone, except that the granules in the latter case are mainly quartz particles. The data points from Desai et al. (1974) for φ > 0 correspond in most cases to a material with a volume fraction of >95% quartz and some admixture of other minerals (e.g. feldspar). The data point with φ = 0 corresponds to pure quartz (from Clauser & Huenges 1995). This is not the same mineral composition as in chondrites, but the heat conductivity of the compact material has a value of 6.15 W K^-1 m^-1 very similar to chondritic material such that the variation with porosity can be directly compared.

An inspection of Fig. 4 shows that our numerical model predicts a variation of K with φ, which is very similar with that observed for sandstone. Only the decrease of heat conductivity with increasing porosity is slightly steeper for sandstone as in our model. Hence, the numerical model with a random distribution of pores is able to reproduce the heat conduction properties of a compacted material. There is only a slight deviation at increasing porosity, which may be attributed to different microscopic structures between our random distribution of the components and the decisively granular structure of sandstone. This is considered in the next section.

4.2. Granular material

We construct a numerical model for a granular medium formed by a packing of balls to study a material with a structure close to chondritic meteorites or sandstone. For this purpose, we simulate the following process: In a box filled with a viscous medium, a large number of balls with prescribed distribution of radii is inserted at the top at random positions (generated with a random number generator). The balls fall under the action of a downwards directed force. The balls interact with the wall with a repulsive force, and mutually with each other, if the distance of the centre of a ball to the wall becomes smaller than the radius of the ball or if the distance between the centres of two balls becomes smaller than the sum of their radii, respectively. The bottom of the box and one of the sidewalls (this suffices) can be in vibrational motion for some period of time to simulate joggling of the box. The equations of motions of the balls are solved by a simple numerical procedure; the details of the model are described in Appendix C.

After switching off the vibrations, the balls settle to a random packing of spheres because of the frictional force. As an example, the random packing of 2800 equal-sized spheres in a box obtained by this procedure is shown in Fig. 5. It is also possible to use some distribution of radii, but for the moment we consider this somewhat simplified model for the chondrules existing in the parent body of chondritic meteorites before the onset of sintering. We cut out the central region from our packing by omitting about two layers adjacent to the wall to get rid of wall effects where the packing fraction is lower than in an unlimited medium. The porosity of the arrangement of balls obtained in this way is found to be 0.35, i.e. the packing corresponds to the random close packing of balls (cf. Gail et al. 2015). This packing of balls can be considered a simple model for the packing of chondrules in chondritic material before sintering.

One can homologously shrink this package by reducing the coordinates of the centres of the balls with respect to the centre of gravity by a fixed factor for all three directions of space. This generates a certain degree of overlap between the balls and correspondingly reduces the total void volume. The model of Arzt et al. (1983; see also the discussion in Henke et al. 2012a) also assumes that sintering a package of balls can be modelled in this way, except that by us the ball radii are held fixed and ball distances are scaled while in the sintering model ball radii and distances both are scaled. In this way, one can simulate the sintering of the ensemble of balls and generate numerical models for different stages of the sintering process.

We decomposed the ball packing into a grid of cubes by counting a cube as belonging to a ball if for at least one ball the midpoint of the cube is located inside a ball. The remaining cubes are counted as void space. We calculated the effective heat conductivity of such a model by assuming the heat conductivity of H chondrite material from Table 6 for the ball material and assigning a small conductivity artificially to the void space, as in the preceding model. Figure 4 and model 1 from Table 7 show the resulting variation of K with φ for this model. As a second case, a packing of 3800 balls with two different radii with ratio 2:1 is also shown in Fig. 4 and by model 2 in Table 7. The bigger number of balls is required in this case to fill about the same volume. We obtain now an improved representation of the porosity dependence of the average heat conductivity of sandstone material with voids filled by air. Thus, the numerical method used here is suited to handle such type of porous materials.

We consider also an empirically determined size distribution of chondrules from H chondrites as given in Friedrich et al. (2015, which give also distributions for other meteorite classes). The data are taken from their Fig. 1 for Hammond Downs (H4). The distribution of chondrule diameters is peaked around a diameter of 0.46 mm with half-width of ~± 0.15 mm. The results for the heat conductivity calculated for such a mixture of spheres is shown in Table 7 as model 3. The results are closer to that of the mono-dispersed case (model 1) than to the binary mixture with a 2:1 size ratio (model 2). In view of the relatively sharply peaked size distribution of the chondrules in H chondrites this was to be expected.

Finally, we obtained from F. Möhlmann, E. Beitz, and J. Blum (working group planetary formation at the Institute for Geophysics and Extraterrestrial Physics at the University of Braunschweig, Germany) data for a sintered cube of originally equal-sized glass beads for which the spatial structure was determined by tomography, resulting in a data file for the decomposition of the body into a grid of cubes. For this, packings of different porosities were generated by shrinking or extending the original spheres. The results of the calculation of the heat conductivity of such granular material are shown in Table 7, as model 4, and in Fig. 4 with a green line with circles. The results are comparable to the other cases considered and the variation of K with φ rather closely follows that of natural sandstone.

The theoretical modelling of heat conductivity of sintered granular material with spheres as primitive granular units and the natural specimens of sandstone show a high conformity in the conductivity variation with porosity. This obviously results from the similarity of the geometric structure of the void space in such materials. Since the material of ordinary chondrites is dominated by nearly spherical chondrules, it is obvious to assume that the variation of heat conductivity of chondritic material with porosity during sintering in an up-heating parent body is the same as found in this calculation.

4.2.1. Approximation for heat conductivity

For the purpose of model calculations it is advantageous to use an analytic expression for the porosity dependence of the heat conductivity. The variation of our results for the mean heat conductivity of partially sintered granular material with varying porosity can be approximated by $\mathit{K}\mathrm{\left(}\mathit{\phi }\mathrm{\right)}\mathrm{=}{\mathit{K}}_{\mathrm{b}}\left(\mathrm{1}\mathrm{-}\mathrm{2.216}\hspace{0.17em}\mathit{\phi }\right)\mathit{,}$(19)where K_b is the bulk heat conductivity of the compact material. This expression goes to zero for φ = 0.451 such that we better clip negative values by letting ${\mathit{K}}_{\mathrm{1}}\mathrm{\left(}\mathit{\phi }\mathrm{\right)}\mathrm{=}{\mathit{K}}_{\mathrm{b}}\begin{array}{c}\max \end{array}\left(\left(\mathrm{1}\mathrm{-}\mathrm{2.216}\hspace{0.17em}\mathit{\phi }\right)\mathit{,}\mathrm{0}\right)\mathit{.}$(20)A zero heat conductivity is unrealistic, however. For high porosity the heat conductivity should tend to the heat conductivity of highly porous material as studied, for example in Krause et al. (2011). We have approximated their results for very porous materials by ${\mathit{K}}_{\mathrm{2}}\mathrm{\left(}\mathit{\phi }\mathrm{\right)}\mathrm{=}{\mathit{K}}_{\mathrm{b}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{1.2}\mathrm{-}\mathit{\phi }\mathit{/}\mathrm{0.167}}\mathit{.}$(21)As described in that paper, we join both approximations into a single prescription for calculating K(φ) in the form $\mathit{K}\mathrm{\left(}\mathit{\phi }\mathrm{\right)}\mathrm{=}{\mathrm{\left(}\mathrm{\right(}{\mathit{K}}_{\mathrm{1}}^{\mathrm{4}}\mathrm{\left(}\mathit{\phi }\mathrm{\right)}\mathrm{+}{\mathit{K}}_{\mathrm{2}}^{\mathrm{4}}\mathrm{(}\mathit{\phi }{\mathrm{\left)}}^{\mathrm{\right)}}}^{\frac{\mathrm{1}}{\mathrm{4}}}\mathit{.}$(22)This interpolates smoothly between the limit cases of low and high porosity. The bulk heat conductivity K_b is calculated for composed material, as described in this work.

4.3. Comparison with meteoritic data

The results of the model calculations for the heat conductivity of chondritic material are compared in Fig. 4 with meteoritic data of some H and K chondrites (Yomogida & Matsui 1983). The dependence of the heat conductivity on porosity found in our numerical models agrees rather well with the experimental results for the case of sandstone with voids filled by air. Our results strongly deviate, however, from the experimental results for meteorites. Obviously there is a strong discrepancy between the numerical results for our model of a porous material and the empirical data. This hints to a fundamental difference between the structure of void space in compacted spherical granules and the void space that is encountered in chondritic material.

In the material that we considered, the voids are interstitials between some kind of granular units; however the structure of void space in the meteoritic material is different. For meteorites, even if the shock stage is low, the material contains numerous microcracks, which are mainly responsible for the observed porosity of ordinary chondrites (Consolmagno et al. 2008). They result from impact induced shock waves that traversed the meteoritic material and fractured the material in a certain volume around the impact site. These cracks are essentially planar structures that act as efficient barriers for the heat flux, even if the corresponding void volume associated with them is small. The most plausible explanation for the different dependencies on φ in materials with a void-space structure such as that in sandstone, on the one hand, and that for meteoritic material, on the other hand, seems to be the strongly different geometry of the voids in both cases. This also means that the pore space measured for meteorites is in large part not pristine but a secondary product of later bombardment of the parent body with boulders since the formation of the body.

4.3.1. Crack-like voids

In order to check the hypothesis that the experimental values obtained for meteorites reflect the large number of cracks, we model such cracks by insulating inclusions with heat conductivity K_i and with a shape corresponding to strongly flattened rotational ellipsoids, which are embedded in an otherwise homogeneous material with heat conductivity K_m. For such a model, the effective heat conductivity K_eff can be calculated from the Bruggeman mixing rule for non-spherical inclusions (cf. Berryman 1995), i.e. $\sum _{\mathit{n}\mathrm{=}\mathrm{1}}^{\mathit{N}}{\mathit{f}}_{\mathit{n}}\mathrm{(}{\mathit{K}}_{\mathit{n}}\mathrm{-}{\mathit{K}}_{\mathrm{eff}}\mathrm{)}\frac{\mathrm{1}}{\mathrm{9}}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathrm{3}}\frac{\mathrm{1}}{{\mathit{L}}_{\mathit{j}}{\mathit{K}}_{\mathit{n}}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{L}}_{\mathit{j}}\mathrm{)}{\mathit{K}}_{\mathrm{eff}}}\mathrm{=}\mathrm{0}\mathit{,}$(23)where the sum n runs over all components in the mixture (inclusions and matrix in our case), f_n is their volume fraction, and K_n their heat conductivity; the sum j runs over the principal axis of the ellipsoid. This is a non-linear equation for K_eff that has to be solved numerically. The L_j are the depolarisation coefficients. These are given for an oblate rotational ellipsoid with principal axis A = B>C by (Bohren & Huffman 1983) ${\mathit{L}}_{\mathrm{1}}\mathrm{=}\frac{\mathit{g}}{\mathrm{2}{\mathit{e}}^{\mathrm{2}}}\left(\frac{\mathit{\pi }}{\mathrm{2}}\mathrm{-}\arctan \mathit{g}\right)\mathrm{-}\frac{{\mathit{g}}^{\mathrm{2}}}{\mathrm{2}}\mathit{,}$(24)with L₂ = L₁ and L₁ + L₂ + L₃ = 1, where ${\mathit{e}}^{\mathrm{2}}\mathrm{=}\mathrm{1}\mathrm{-}{\mathit{C}}^{\mathrm{2}}\mathit{/}{\mathit{A}}^{\mathrm{2}}\mathit{,} \mathit{g}\mathrm{=}\frac{\sqrt{\mathrm{1}\mathrm{-}{\mathit{e}}^{\mathrm{2}}}}{\mathit{e}}\mathrm{·}$(25)We assume for the matrix K_m = 4.3 W m^-1 K^-1 as derived by extrapolation to zero porosity for the chondrites considered in Yomogida & Matsui (1983) and K_i = 10^-4 W m^-1 K^-1 for the insulating inclusions, and we vary the axis ratio from A:C = 1:1 to A:C = 1000:1. The resulting run of K_eff with porosity is shown in Fig. 4 as the lilac lines. Most experimental values are within a factor of about two near the line with A:C = 100:1. This suggests that, indeed, the porosity variation of the heat conductivity observed for chondritic material is determined by crack-like voids in the meteorites.

5. Application to H chondrite parent body

5.1. Porosity of surface layers

According to our present knowledge, planetesimals form by some process that assembles free floating dust and chondrules from the accretion disk into planetesimals. The details of this process are not important for our problem, only the fact that initially the material should be a loosely packed bimodal granular medium consisting of about mm-sized granules and less than μm-sized dust.

The growth of bodies to the size of the parent bodies of meteorites involves, according to present models, low relative velocities of ≲1 km s^-1 (see Johansen et al. 2014, 2015, for a review on asteroid formation). These velocities are much too low to shatter the dust particles and chondrules. It is only if Jupiter has formed that relative velocities in the region of the asteroid belt are pumped up to the order of 5 km s^-1 where collisions become destructive. At that time obviously bodies with sizes up to several hundred km had already formed from which the parent bodies of the chondrites are the survivers. The initial structure of the material from which the bodies are built, thus, is generally assumed to be a granular material, which is to some extent compressed by the action of self-gravitation of the body and by collisions during the growth process.

Heating by decay of radioactive nuclei raises the temperature within less than 1 Ma after formation above the critical temperature of ~900 K, where sintering removes the pores from the material (Gail et al. 2015) resulting in compact rocky material. Under certain circumstances the central temperature may become even high enough for melt migration of iron such that they become partially differentiated (Elkins-Tanton et al. 2011). However, there remains an outer layer where temperatures remain insufficient for sintering. This layer preserves the initial granular structure. Such incompletely compacted but thermally metamorphosed material has been found to form the material of a number of ordinary chondrites and has been studied with micro-tomography in some detail by Friedrich et al. (2008), Sasso et al. (2009), Friedrich & Rivers (2013), and Friedrich et al. (2014). Such chondrites have large porosities and most of their pore space is intergranular; abundant microcracks are not found. A detailed discussion of such chondrites is given in Friedrich et al. (2014) who argued that they are “rare materials that survived the earliest stages of impact processing on ordinary chondrite parent bodies”. The existence of such material shows that the initial texture of the material of the parent bodies resembles the structure of granular material as studied in the preceding section.

After formation of Jupiter (at some still unknown instant, but probably later than 2 Ma after CAIs) the relative velocities of bodies in the asteroid belt are stirred up and hypervelocity impacts begin, which result in cratering of the surface of the bodies and formation of an outer regolith layer. The formation of the regolith layer is obviously accompanied by extensive formation of cracks in the material. If one models the thermal evolution of the parent bodies of meteorites, one has to discriminate between the structure of the void space of the initial sand-like mixture of chondrules and matrix from which the body formed and the void space structure induced in material compacted by impacts. The heat conductivity and how it depends on void space is very different in both cases.

It is not possible to use the values measured for meteorites for model calculations of the early thermal evolution of meteoritic parent bodies, as has been done up to now in our previous papers and by others (Yomogida & Matsui 1984; Bennett & McSween 1996; Akridge et al. 1998; Harrison & Grimm 2010; Henke et al. 2012a,b; Henke et al. 2013). For the earliest phase of evolution one has to use instead the results for compacted porous material for which sandstone seems to be an analog. Later the properties of the surface material are modified by impacts and abundant cracks are formed. It is then necessary to consider how the properties of the outer few km of the body develop under the action of a bombardment from outside. In particular, how the time to form an impact caused regolith layer is related to the characteristic heating and cooling time of the outer layers is important.

The formation of a regolith layer on the surface of asteroidal bodies has been studied by Housen & Wilkening (1979), Ward (2002), and Warren (2011). Explicit results for the time required to build up such a layer are only given in the paper by Housen & Wilkening (1979). These authors considered impacts on basaltic and sand-like material, but they only give results for impacts into basaltic material for bodies of 100 km and 300 km in diameter. From their figures one infers that the formation of a full-fledged regolith layer lasts 1 ... 2 Ga. This is significantly longer than the ~100 Ma period for the essential part of the thermal evolution of the parent bodies. During this initial period, only a thin regolith cover of the surface forms, which means within the frame of their model that only a small fraction of the surface suffered cratering by impacts. Since their model deals with average values, the real covering of the surface by craters and their ejecta during the initial growth phase should be rather patchy and large parts of the surface are likely not yet affected (cf. also Ward 2002,Fig. 11). By a simple test calculation, we show later (see Fig. 6) that a thin early regolith layer does not strongly affect the initial thermal evolution of a parent body and can be neglected in zero order approximation, but it is clear that the porosity evolution by impacts should be included in future model calculations.

5.2. Evolution model

5.2.1. Model calculation

Our results for the heat conductivity of the chondritic material are used to calculate thermal evolution models of planetesimals. The method for the model construction is described in Henke et al. (2012a). The numerical method for solving the heat conduction equation is now changed from a method based on finite differences to a method based on finite elements. The improved treatment of sintering described in Gail et al. (2015) is included. The instantaneous formation approximation is used since it was found that this fits best to the observed cooling properties of H chondrites (Henke et al. 2013). Additionally, melting of the Fe,Ni-FeS complex and of silicates is implemented (the method for modelling is described elsewhere), but no melt migration is considered. Melting is included in the modelling because the parent body of ordinary chondrites could have been partially molten in their innermost core regions, although presently no chondrites are known that show indications of partial melting and, at the same time, seem to originate from one of the parent bodies of the ordinary chondrites.

The basic parameter values used for the following model calculations are listed in Table 8. We consider parent bodies of ordinary chondrites where chondrules strongly dominate over matrix material. The basic properties of this mixture are listed in Table 1. It is assumed that the properties of the initial chondrule matrix mixture do not vary within the body. The initial porosity is calculated as described in Sect. 2.1 (see Table 1). The sinter properties of the material of ordinary chondrites are dominated by the chondrules. The method for calculating the sintering of such material and the set of parameters describing the properties of the material are described in Gail et al. (2015).

For the model calculations, the dependence of the average heat conductivity on porosity is calculated according to Eq. (22). For the bulk heat conductivity K_b, we use the value calculated for material with the average composition of H chondrites given in Table 6. For the temperature dependence of K_b, we assume a provisional dependency ∝T^−1/2, which is a reasonable approximation for rock material (e.g. Xu et al. 2004). A more detailed investigation of the temperature dependence will be published elsewhere.

The heat production in the parent bodies is determined by the decay of short- and long-lived radioactive nuclei. The heat production by ²⁶Al and ⁶⁰Fe governs the initial period of thermal evolution. Their initial abundances at time of CAI formation used in our calculations are given in Table 8. The reasons for this choice are discussed in Henke et al. (2013). The calculation also includes heating by long-lived radioactive nuclei as described in Henke et al. (2012a). They are important for the late thermal evolution.

5.2.2. Sample calculation

We compare models for the thermal evolution of a planetesimal of radius 150 km and formation time 2 Ma after CAI calculated with two different assumptions with respect to the porosity dependence of the heat conduction. The surface temperature is set in all cases to 200 K. We denote the models using the porosity dependence of heat conductivity as for granular material as model G, and with porosity dependence as derived from the meteoritic data of Yomogida & Matsui (1983) and Opeil et al. (2012) as model Y.

Figure 6 shows the temperature evolution at depths of 3, 10, and 30 km beneath the surface and the temperature evolution at the centre of the body. The temperature evolution in both cases is significantly different as a result of the strongly different blocking of heat flux by pores in a granular material and in cracked material. Although the heat conductivity is identical after all voids are closed by sintering, the structure of the outer layer where temperatures never became sufficient for sinter compaction of the material is different for both cases. The surviving porous layer extends much deeper in the case where we apply the porosity dependence of K for granular material than for the case of cracked meteoritic material. Correspondingly, the thermal shielding effect of the outer porous layer is less marked for the granular layer because of its higher heat conductivity. This affects the temperature structure in the model in the sense that, except for the shallowest layers, the descending part of the temperature curve at a certain depth is shifted in the model to earlier times for the granular case by about 30 Ma compared to the case of impact fractured material.

Figure 7 shows the resulting maximum temperature achieved at some radius during the thermal evolution of the body. This determines the compaction and thermal metamorphism of the material as it is reflected by the different petrologic classes of meteorites. One recognises the significant differences in the depth of the surviving porous outer region between the two kinds of models. The larger extent in model G of the porous outer layer from which meteorites of petrologic type 3 originate appears to be more realistic than the very thin layer in model Y because the outermost 1 km surface layer of a small body is probably significantly eroded by impacts during the 4.5 Ga evolution of the asteroid belt, as may be inferred from the models in Housen & Wilkening (1979).

One conclusion we can draw from Fig. 6 is that the temperature evolution inside the body, as far as it is of interest for the thermochronometers, is completed within about 100 Ma. If during this period a disruptive collision occurs, this completely changes the temperature history of the remaining fragments, even if they re-assemble to a rubble pile. If no such catastrophic event occurs within the first 100 Ma, strong impacts may result in strong local disturbances of the thermal evolution of the matter, but they do not result in a global change of the thermal history of the material. If disruptive impacts occur after the critical initial evolutionary phase of about 100 Ma, the bigger fragments should preserve the information on the initial temperature evolution in their interior. Then at least some of the meteorites derived from the bigger fragments preserve the record on the thermal history within the crucial initial ~100 Ma period.

For a crude estimation of the effect that a growing impact-modified outer layer has for the thermal evolution, we consider a very simple approximation. We follow the results of Housen & Wilkening (1979) for a body with 300 km diameter and strength scaling and assume that a regolith layer grows linear in time to 2 km thickness within 200 Ma. An initial porosity in this layer is assumed to equal 0.25 as for the granular material. The porosity dependence of the heat conductivity in this impact-generated covering layer is assumed as for meteoritic matter (model Y). Figure 6 shows the temperature evolution for the centre and at three selected depths in such a model. The temperatures in this model are close to those in model G for the initial ~100 Ma of evolution and only start to deviate somewhat for later times. In all cases, the deviations remain ≲30 K for temperatures >400 K, which is less than the accuracy with which the closure temperatures of thermochronometers are known (see Fig. 8 or Table 5 of Gail et al. 2015, for closure temperatures of the thermochronometers used in this work). The closure times derived from thermochronometers, therefore, do not seem to be significantly influenced by a gradually growing impact-generated regolith layer. For the moment we feel this finding justifies ignoring the whole subject, but the impact evolution of the surface layers clearly needs a detailed study, which is out of the scope of this paper.

5.3. Fit of thermochronological data of H chondrites

We apply our calculation of the heat conductivity of chondritic material to the problem of reconstructing the parent body of the H chondrites from thermochronological data of a number of selected chondrites of class H. We search for a model that reproduces the thermal history of a number of H chondrites for which their individual cooling history is sufficiently constrained by empirical data. We have discussed this problem in our earlier papers (Henke et al. 2012a,b, 2013; Gail et al. 2015) and will not repeat the details. We base our calculation on the same set of nine meteorites and 37 data points for closure ages, as described in Gail et al. (2015). There are still only nine meteorites available that satisfy the requirement that their individual cooling history is determined by closure ages for at least three different thermochronometers.

The basic physical parameters (G, φ₀, etc.) assumed for the models are given in Table 8. The basic parameters that determine the properties of an evolution model of a planetesimal (radius R_p, formation time t_form, surface temperature T_srf) are determined by comparing cooling curves with empirical cooling histories of meteorites. For a given model, the unknown burial depths within the parent body from which the meteorites originate are determined for each meteorite via an optimisation method such that the temperature evolution at the burial depth fits as close as possible to the thermochronological data points of that meteorite. The parameters determining the evolution model of the parent body are varied by applying an evolution algorithm such that an overall best fit for the set of all meteorites is obtained in such a way that the calculated deviation χ² between the model and data is lowest.

The result of the optimisation process is shown in Table 9 and in Fig. 8 as model G. For comparison we ran the same optimisation with the porosity dependence of the heat conductivity as determined for meteorites (model Y). The first group of data in the table show the optimised values of the model parameters (R_p, t_form, T_srf) and the resulting maximum central temperature and thickness of the residual porous layer. We also show the best-fit model H from Gail et al. (2015) for comparison; this model differs from model Y only in so far as some estimated value is used rather than the calculated bulk heat conductivity K_b of the material of H chondrites (Table 6).

The radius and formation time of the models are very similar for all the models. They are essentially constrained by the cooling history at large depth, i.e. by the data for the meteorites of high petrologic type. The temperature evolution is almost the same for the deep layers of the body, as is readily seen in Fig. 8. The differences are more pronounced at shallow outer layers where the structure of the model is different for model G and Y. As already shown by the sample calculations, the residual porous layer in model G is significantly thicker than in model Y as result of the higher value of the heat conduction coefficient for sintered granular material than for a material where heat flux is efficiently blocked by cracks. The temperature lapse rate close to the surface therefore is shallower in model G than in model Y and the temperature for efficient sintering of the chondrules at ≳900 K is only encountered at a higher depth in model G than in model Y.

The strong difference in heat conductivity in the porous outer layer is also the reason for the higher surface temperature required in model G than in model Y because otherwise cooling would be too rapid in the layers where the H4 chondrites come from. The cooling history of the low petrologic classes mainly constrains the surface temperature.

On the whole, the model using the new model for heat conductivity of pristine chondritic material results in a somewhat better fit quality, expressed by a lower value of χ²/ (N−p) (N = number of data, p = number of parameters). To some extent, this supports our hypothesis that the heat conductivities measured for meteoritic material cannot be applied to the early evolutionary phase of the parent bodies because the microstructure of the material is strongly impact modified over the long period between completion of the initial heating and cooling phase and the investigation of the meteorite in the laboratory.

6. Concluding remarks

We have studied the heat conductivity of chondritic material. It is shown that it is possible to determine, for ordinary chondrites, the effective heat conductivity coefficient from the individual properties of the components of the complex mixture of minerals and metal in chondrites by solving the heat conduction equation for a test body. We found that the effective heat conductivity of the compact material can be approximated with sufficient accuracy by the Bruggeman mixing rule. This reproduces the result of the numerical calculation with an accuracy of better than 1 percent. Hence the value of the heat conductivity coefficient of compact material can be calculated equally well from the Bruggeman mixing rule.

We have also studied the effective heat conductivity of porous material by numerically solving the heat conduction equation for a test body. We considered a model with a stochastic distribution of voids, and a model in which a granular material, formed from spheres with some different mixtures of sphere radii, is numerically generated and artificially sintered. It is found that the variation of effective heat conductivity strongly deviates from what is found in laboratory measurements of heat conductivity of meteorites, while the calculations can reproduce measured values for natural sandstone. It is argued that the discrepancy can be explained by the fact that the pore space in meteorites is dominated by impact-generated cracks, which have a sheet-like structure that differs completely from the geometry of voids in compacted granular material formed from spheres. We demonstrate, by applying the Bruggeman mixing theory for disk-like voids in a porous material, that the difference between calculated and measured heat conductivities can be readily explained by microcracks in meteoritic material.

We argue then that meteoritic material does not represent the pristine material present in the earliest evolutionary phase of the parent bodies of asteroids because, according to theory, the build up of a heavily impact modified surface layer, where numerous cracks are generated in the material, requires timescales longer than the critical period for the thermal evolution of the parent bodies of chondrites. Laboratory data for heat conductivity of chondritic material are only representative for present day asteroids. The pristine material present during the initial thermal evolution seems to exist no more, and one has to determine its properties from calculations, as we did, or one could produce a crack-free material in the laboratory and measure its properties.

Our calculations also show that the porosity dependence of the granular material may be calculated with no more than about 10 percent deviation from the numerical results from the Bruggeman mixing rule. On the other hand, the behaviour of the meteoritic material may be modelled by the Bruggeman mixing rule for disk-like inclusions. We also give an analytic approximation to our numerical results, which may be applied in model calculations.

We repeated our earlier fits of thermal evolution models for the H chondrite parent body with thermochronological data for the set of nine H chondrites with rather well-defined cooling histories and determined the putative radius and formation time of the body. The results are close to our earlier results based on an effective heat conductivity as measured for meteorites. With the new data on heat conductivity, the quality of the fit is, however, somewhat improved.

The discrepancy between the results of laboratory measurements and theoretical predictions for the effective heat conductivity strongly suggests that attention has to be paid to the question of what are the properties of the pristine material in planetesimals. It is obviously necessary to study in detail the impact evolution of the surface layers before one can come to a final conclusion on the thermal history of meteorites and their parent bodies.

The model for porosity of meteoritic material used in this paper cannot be applied to the uppermost few meters below the surface where space weathering is active. There one probably encounters material properties for which our model does not make allowances (e.g. a very high porosity; Emery et al. 2006; Vernazza et al. 2012) and that requires a different treatment.

Acknowledgments

We greatly acknowledge that E. Beitz, F. Möhlmann, and J. Blum shared their data on the structure of the sintered body of glass beads with us. This work was performed as part of a project of “Schwerpunktprogramm 1385”, supported by the “Deutsche Forschungsgemeinschaft (DFG)”. M.T. acknowledges support by the Klaus Tschira Stiftung gGmbH. This research has made use of NASA’s Astrophysics Data System.

References

Appendix A: The numerical model

The heat conduction equation is numerically solved for a cube of edge-length L. We chose for this purpose an implicit finite differences method, which is about an order of thousand times faster on a single cpu than an explicit method since there is no restriction on the time step length in implicit methods to guarantee stability. The cube is divided into a grid of n³ equal-sized subcubes of edge length h = L/n, where n is the number of grid cells in each dimension.

Writing Eq. (11) in implicit finite differences yields ${\mathit{T}}_{\mathit{ijk}}^{\mathit{n}\mathrm{+}\mathrm{1}}\mathrm{-}{\left(\frac{\mathrm{1}}{\mathit{\rho }{\mathit{c}}_{\mathit{p}}}\right)}_{\mathit{ijk}}\mathrm{·}\left(\mathrm{\nabla }\mathit{K}\mathrm{·}\mathrm{\nabla }\mathit{T}\right)\begin{array}{c}\mathit{n}\mathrm{+}\mathrm{1}\\ \mathit{ijk}\end{array}\mathrm{·}\mathrm{\Delta }\mathit{t}\mathrm{=}{\mathit{T}}_{\mathit{ijk}}^{\mathit{n}}\mathit{,}$(A.1)where Δt is the time step size and $\begin{array}{ccc}{\left(\mathrm{\nabla }\mathit{k}\mathrm{·}\mathrm{\nabla }\mathit{T}\right)}_{\mathit{ijk}}& \mathrm{=}& \\ & & {\left(\frac{{\mathit{K}}_{\mathit{i}\mathrm{+}\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{T}}_{\mathit{i}\mathrm{+}\mathrm{1}}\mathrm{+}\left({\mathit{K}}_{\mathit{i}\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{-}{\mathit{K}}_{\mathit{i}\mathrm{+}\mathrm{1}\mathit{/}\mathrm{2}}\right){\mathit{T}}_{\mathit{i}}\mathrm{+}{\mathit{K}}_{\mathit{i}\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{T}}_{\mathit{i}\mathrm{-}\mathrm{1}}}{{\mathit{h}}^{\mathrm{2}}}\right)}_{\mathit{jk}}\\ & & \mathrm{+}{\left(\frac{{\mathit{K}}_{\mathit{j}\mathrm{+}\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{T}}_{\mathit{j}\mathrm{+}\mathrm{1}}\mathrm{+}\mathrm{(}{\mathit{K}}_{\mathit{j}\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{-}{\mathit{K}}_{\mathit{j}\mathrm{+}\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{)}{\mathit{T}}_{\mathit{j}}\mathrm{+}{\mathit{K}}_{\mathit{j}\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{T}}_{\mathit{j}\mathrm{-}\mathrm{1}}}{{\mathit{h}}^{\mathrm{2}}}\right)}_{\mathit{ik}}\\ & & \end{array}$Here the index n denotes the nth time step of length Δt, and i,j,k the indices of the spatial coordinates in each spatial direction, and h = L/n denotes the spatial step size (grid constant). Equation (A.2) results in a large and sparse coupled system of equations with n³ unknowns that is solved by a conjugate gradient method.

The quantities K_{i ± 1/2} denote the heat conductivity at the midpoint of the interval between index i and that of index i ± 1. One possibility is to define this as the arithmetic mean ${\mathit{K}}_{\mathit{i}\mathrm{\pm }\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{=}\mathrm{(}{\mathit{K}}_{\mathit{i}\mathrm{\pm }\mathrm{1}}\mathrm{+}{\mathit{K}}_{\mathit{i}}\mathrm{)}\mathit{/}\mathrm{2}\mathit{,}$(A.6)and another possibility is the harmonic mean ${\mathit{K}}_{\mathit{i}\mathrm{\pm }\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{=}\mathrm{2}\mathit{/}\mathrm{(}{\mathit{K}}_{\mathit{i}\mathrm{\pm }\mathrm{1}}^{-1}\mathrm{+}{\mathit{K}}_{\mathit{i}}^{-1}\mathrm{)}\mathit{.}$(A.7)The harmonic mean is a more appropriate treatment if the heat conductivities of the different components used vary significantly but cannot be used if one component has a heat conductivity of zero. In this paper, the harmonic mean was used in all cases where the heat conductivity of chondritic material is calculated. For the calculations of the heat conductivity of systems of spheres in a vacuum, the arithmetic mean is used.

Figure A.1 shows a comparison of the resulting effective heat conductivites resulting from the two means for a two-component material with different component mixing ratios. Plotted are results of both means at grid resolutions of n = 100, 200, and 300. It is easily seen that the results for both means approach each other for increasing resolution, and this approaching is mainly seen for the solution based on the arithmetic mean, whereas the K_eff for the harmonic one stays basically the same. From this behaviour we conclude that using the harmonic mean in our models is the better choice if we do not consider areas with zero heat conductivity.

As boundary conditions for solving the equation, we have constant temperatures T₁ and T₂ with T₁ ≠ T₂ on two opposite sides of the cube, i.e. a Dirichlet boundary condition. On the other four sides, we describe Neumann boundary conditions. This boundary condition assumes that there is no effective heat flux through the side walls of the cube. This is certainly not correct for an individual grain located at the boundary because of certain local fluctuations of the heat flux. Strictly speaking this boundary condition can only apply on the average over larger pieces of the boundary area, i.e. after we have already performed the averaging procedure. In practice, however, we wish to consider a situation where averaged quantities vary over length scales, which are much larger (kilometers) than typical sizes of individual particles (μm to mm). In that case local deviations of heat flux from the mean must be very small such that a Neumann boundary condition appears to be a reasonable approximation.

As an initial condition, any temperature distribution throughout the cube can be prescribed that fulfills the boundary conditions. Here we prescribe the mean between T₁ and T₂ for each grid cell not belonging to the front surfaces with constant boundary temperature.

To determine the effective heat conductivity, this model is then run until a stationary state establishes for the temperature distribution. We assume a stationary state to be achieved, if the maximum temperature change in any grid cell for a time step is less then 10^-7 K/s. In that stationary state the heat flow does not change any more and can be used to determine the effective heat conductivity k_eff Eq. (12). The total heat flow F through the front sides can be determined by summing over the heat flows of each grid cell that belongs to a certain boundary surface of constant temperature.

The model can be run for any resolution, whereas a resolution of 200 × 200 × 200 is a good compromise between a high enough resolution and a reasonable computation time. For a cube with L = 1 cm and an effective heat conductivity of the order of 10 W/mK, the physical time until a stationary state establishes is of the order of 1 min.

Appendix B: Random distribution of materials

A random distribution of components with given volume fractions is generated via the following procedure:

• 1.

  The material is supposed to consist of N components. These components are to be distributed randomly within the cube with a prescribed volume fraction for each component.

• 2.

  In the first step, the cube is decomposed into subcubes. The number of grid cells, n, for each direction is chosen and the set of volume elements of the decomposition is generated. All grid-cells are supposed to be filled with one of the components, now called component 1.

• 3.

  In the second step, centres of equal-sized balls with a prescribed radius R are generated at random using a random number generator. For each grid cell for which its centre lies within such a ball, its component is changed to component 2. The balls are generated one by one until the number of cells filled with component 2 attains the prescribed volume fraction of this component. This results in a cube with a matrix component (component 1) containing some balls consisting of the second component. In this procedure overlapping of balls is not excluded, resulting in arbitrary shapes of component 2 aggregates within the cube if the number of balls and their radius is such that there is significant overlap.

• 4.

  If there are three or more components, the preceding step is repeated for each of the further components. For each grid cell where its centre is located in one of these new balls and its component is still component 1, its component is changed to the new component.

This procedure allows us to generate a random spatial distribution of a number of different materials inside the cube with a prescribed volume fraction of each component. By means of construction, the prescribed volume fraction can only be fitted approximately because of the final increments of component volume during the construction process. Different constructions starting with different seeds of the random number generator give slightly different volume fractions around its nominal value. This resembles what one observes in the compacted material of chondrites of petrologic types 5 and 6, and we take such artificially generated distributions as a model for compacted chondritic material.

The distribution of components generated in this way does not depend on the grid resolution of the finite difference method but only on the assumed radius R of the balls and the volume fractions of the components. Thus, starting with the same sets of balls for the different components, the same spatial distribution of components with different decompositions of the cube into n³ grid cells can be generated for calculating the effective heat conductivity by varying n. In addition, the radius R of those balls can be varied to define component structures of varying sizes.

Appendix C: Construction of a granular medium

A box filled with balls either of a fixed radius R or with radii from some assumed radius distribution is considered.

A rectangular Cartesian coordinate x, y, z is used. The coordinates of the walls perpendicular to the x, y, and z axis are X₀, X₁ and Y₀, Y₁ and Z₀, Z₁, respectively. The width D of the box

in the x and y direction have to be chosen as a suited multiple of the (average) ball radius. In case of fixed ball radii, D must not be a multiple of R because the system tends to form an unwarranted cubic sphere packing in this case. The height of the box in z direction must be at least twice the width D.

It is assumed that there is some downwards acceleration, ${\mathit{A}}_{\mathit{z}}^{\mathrm{\left(}\mathrm{g}\mathrm{\right)}}\mathrm{=}\mathrm{-}\mathit{g,}$(C.1)acting on the balls and that there is some frictional acceleration, ${\mathit{A}}_{\mathit{x}}^{\mathrm{\left(}\mathrm{fr}\mathrm{\right)}}\mathrm{=}\mathrm{-}{\mathit{u}}_{\mathit{x}}\mathit{/}\mathit{\tau ,}$(C.2)where u_x is the x component of the velocity and τ a given stopping time, and analogously for the other two directions. The balls experience a repulsive acceleration on contact given by ${\mathit{A}}_{\mathit{x,ij}}^{\mathrm{\left(}\mathrm{r}\mathrm{\right)}}\mathrm{=}\{\begin{array}{c}\\ & \mathrm{if} {\mathit{\rho }}_{\mathit{ij}}\mathit{<}{\mathit{R}}_{\mathit{i}}\mathrm{+}{\mathit{R}}_{\mathit{j}}\\ & \mathrm{else}\end{array}$(C.3)and analogously for the other two components, where ${\mathit{\rho }}_{\mathit{ij}}\mathrm{=}{\mathrm{\left(}\mathrm{\right(}{\mathit{x}}_{\mathit{i}}\mathrm{-}{\mathit{x}}_{\mathit{j}}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{(}{\mathit{y}}_{\mathit{i}}\mathrm{-}{\mathit{y}}_{\mathit{j}}{\mathrm{)}}^{\mathrm{2}}\mathrm{+}\mathrm{(}{\mathit{z}}_{\mathit{i}}\mathrm{-}{\mathit{z}}_{\mathit{j}}{\mathrm{)}}^{\mathrm{2}}\mathrm{)}}^{\mathrm{1}\mathit{/}\mathrm{2}}$(C.4)and x_i, y_i, ... are the coordinates of the centres of the balls, R_i and R_j the radii of the balls, and α_b is a suited constant. If the distance of a ball centre to one of the walls is less than the radius, the same kind of repulsive acceleration ${\mathit{A}}_{\mathit{x,i}}^{\mathrm{\left(}\mathrm{w}\mathrm{\right)}}\mathrm{=}{\mathit{\alpha }}_{\mathrm{w}}\frac{{\mathit{x}}_{\mathit{i}}\mathrm{-}{\mathit{x}}_{\mathit{j}}}{{\mathit{R}}_{\mathit{i}}}$(C.5)is applied.

The following equations are solved for the balls $\begin{array}{c}\mathrm{¨}\\ \mathit{x}\end{array}\mathrm{=}{\mathit{A}}_{\mathit{z}}^{\mathrm{\left(}\mathrm{g}\mathrm{\right)}}\mathrm{+}{\mathit{A}}_{\mathit{x}}^{\mathrm{\left(}\mathrm{fr}\mathrm{\right)}}\mathrm{+}{\mathit{A}}_{\mathit{x,ij}}^{\mathrm{\left(}\mathrm{r}\mathrm{\right)}}\mathrm{+}{\mathit{A}}_{\mathit{x,i}}^{\mathrm{\left(}\mathrm{w}\mathrm{\right)}}$(C.6)with a simple leap-frog method. The balls are introduced one after another at the ceiling of the box with random x, y coordinates at short intervals.

One of the sidewalls and the bottom are vibrating, e.g. the X₀ and Z₀ coordinates are periodically moved with slight amplitude. This simulates a shaking of the box, which in turn induces a chaotic motion of the balls. The shaking is gradually switched off some time after the last ball is introduced. At the same time the position of the ceiling is gradually lowered until it fits to the upper layer of balls.

The system evolves by friction to a final state where all balls are at rest. Rather big time steps are taken for integrating the equations of motion because it is not necessary to achieve an accurate orbit of the balls; in any case, the assumed forces are artificial. Only the final stage is of interest after all balls have come to a rest. This procedure results in a random packing of balls. The parameter values typically used in a run are given in Table C.1.

All Tables

All Figures