1 Introduction

Due to increasing costs for energy and growing environmental awareness, energy and material efficiency are of ever-increasing importance in product development. One of the main causes for energy and material loss is friction in tribological contacts. In recent years, polymer-based materials have been used for various tribological applications in increasing quantities [9]. These materials combine the potential for lightweight design, dry running capability and cost-effective processing. Typically, a polymer-based tribological composite, produced by compounding, consists of a high-performance polymer which is reinforced with fibres to increase mechanical strength and wear resistance and modified with an internal solid lubricant to reduce friction [22, 30]. Further improvements of both friction and wear behavior can be achieved by the addition of hard micro- and nano-sized particles which improve the stiffness of the matrix material and interact with the carbon fibres in several ways [5, 6, 17, 29, 30].

The development of these hybrid composites is expensive and time consuming, due to the large number of experiments required to determine the tribological behavior of numerous formulations under different load conditions. Also, the tribological performance of polymer compounds is strongly dependent on the geometry of the tribological system. Consequently, results of early stage standardized tests may differ considerably from the results obtained in more application-oriented test setups. A reduction of the extend of experimental work and an improved transferability of results may be achieved using numerical methods to simulate tribological systems. This requires suitable homogenized models for the material properties, in particular for mechanical and thermal properties. In the literature, many mathematical methods for constructing effective media descriptions are available, generally called homogenization methods. These methods fall into two broad classes: Analytical homogenization methods are based on empirical observations or take advantage of physical understanding and mathematical derivations, usually asymptotic expansions. Numerical homogenization methods apply simulation techniques for the corresponding boundary value problems (e.g. for thermal conductivity or for mechanic deformation) on micro-scale Representative Volume Elements (RVE). Especially for fiber-reinforced materials, a number of research articles deal with homogenization techniques, mainly in the field of homogenization of the mechanical aspects of composites[2,3,4]. For the numeric homogenization of the Effective Thermal Conductivity (ETC) recent works deal with image-based models [10], laminates [19], cell-based models [16] and the FEM [15]. A more comprehensive review of recent work in this field was provided by Zhai et al. [28]. In this work, we propose a homogenization methodology implemented in two different numerical setups.

In the case of fillers with an aspect ratio of around 5 or above, the filler orientation inside the composite also needs to be considered. In this work a numerical homogenization method for the thermal conductivity of composites filled with fibres and particles is demonstrated. The main focus is to identify which aspects of the matrix-filler system are most relevant for the thermal properties.

In Section 2 the model for the orientation of non-spherical fillers in a composite will be discussed, since it requires some significant assumptions. The orientation distribution for any non-spherical filler phase is dominated by the production process and depends on the composition mainly via the viscosity.

In Section 3 and Section 4 the numerical homogenization approach is described and validated. The results are important to qualify the reliability of the studies performed thereafter.

In Section 5 several parameter studies are performed, to gain insight, which geometric and material properties have a large influence on the result. The investigated parameters encompass fundamental aspects, such as the gross composition of the composite or the thermal conductivity of the involved materials, as well as less well-structured properties, e.g. the length and orientation of the fibres in the composite. The gross composition of extruded composites is usually known rather accurately, but the production process may not guarantee a totally homogeneous dispersion of the filler (fibres) and other particles (agglomeration) in a polymer matrix. The thermal conductivity of the matrix material and the composite can generally be measured to a similar precision, but for the filler phases this is more difficult. E.g. for carbon fibres values ranging from 4.9 W/mK [26] to 17 W/mK [8] are reported. The length and orientation of fibres in a composite is generally difficult to assess. Even if the original length is known, the compounding and extrusion process often involves fiber cracking, so that the final length distribution in the composite is not directly available. At the end of Section 5 some minor variations of detail behavior like the size and thermal conductivity of filler particles, and issues with fiber-matrix bonding are investigated.

Section 7 concludes the numerical findings of this work on a structure-property relationship of the homogenized heat conductivity of the composites under investigation basically as a function of key characteristic parameters of the composite geometry. This relationship has been successfully tested against “worse case“ geometries not included in the fitting procedure.

2 Geometrical composite modeling

The orientation of the fibres in the composite is the most complex parameter to describe, because it cannot easily be reduced to a single number. The most direct way to obtain a RVE for homogenization, is to use an actual 3D image of the material. This approach however can only be applied to materials that actually exist (so it is not applicable in design of new materials) and it typically involves high costs (for a CT) and further difficulties if the contrast between the matrix material and the fillers is low (like carbon fibres in a polymer). Alternatively one could use a simulation of the intended production process (CFD simulation of polymer melt[21]) to predict a geometry for the system. Such a simulation would respect the influence of the production process on the geometry, but it requires a lot of material properties, and Hopmann [13] showed, that default values provided by typical CFD software packages often yield unsuitable results. The most feasible approach is to choose a suitable (parametrized) distribution for the directions of the fibres and adjust the parameters to the composite geometry. Distributions with many parameters should be avoided, because they would require a lot of data to actually estimate these parameters for a real (or even just planned) material.

A widely used assumption for injection molded parts made from ber reinforced plastics is the existence of three zones with different ber orientations along the thickness direction of the part. In the outer zones, the bres are predominantly oriented in the direction of molding due to shear of the polymer melt, while in the core they are mainly oriented perpendicular to the direction of molding due to extensional ow [11, 18]. Here, micro geometries are generated in a representative volume element (RVE) for each zone separately. As a basis for the geometry model, the direction distribution proposed by Schladitz [25] for non-woven materials is used. This model depends only on a single parameter. It represents distributions with one preferred direction while assuming symmetry in the two orthogonal directions. To relax the symmetry assumption, a second model, known as Angular Central Gaussian (ACG) [27] was also investigated. The covariance matrix dening the Gaussian distributions is assumed to be diagonal, which ts into the homogenization procedure as described in Section 3 below. In this case, the model depends on two parameters, including the model by Schladitz [25] as a special case. We note that the ACG distribution is an exact solution of the Fokker Planck equation describing the orientation density function dynamics, and represents an exact second order moment closure of the Folgar-Tucker equation, see [1, 12] for a discussion. The placement of already oriented bres and particles in the volume uses a random sequential adsorption process, see [23].

The parameters for both models could be estimated again from CFD-type simulations, see [24], or directly from material samples e.g. from microsection images. However, because only a very small number of parameters is involved, it is much simpler to quantify the influence of these parameters. Therefore the transfer of orientation information obtained for one existing compound to other compounds of moderately different composition is simpler.

For our study six different geometric parameter settings were considered, the parameters are listed in Table 1. In the last three columns the diagonal elements of the second order orientation tensor (FOT) [20] are added to give an idea how the geometries behave.

Table 1 Fibre orientation model parameters and second order orientation tensor diagonal elements for all 6 geometries investigated in this study. The Schladitz model is a special case of the ACG model, so that effective parameters for the other model can be given for B1-B4, but A1 and A2 can not be assigned a β parameter
Full size table

Using diagonalization and an appropriate coordinate transformation to the eigenbasis of the ACG tensor, we may assume, without loss of generality, that the coordinate axes of all RVEs used in this work represent the eigenvectors.

3 Mathematical homogenization procedure for the heat conductivity

In this work, a homogenized heat conductivity \(\lambda _{i}^{hom}\) in direction i, in a cubical RVE is calculated in the following setup:

For the sake of notation the unit cube (0, l)3 is assumed for the RVE. Inside the cube some parts are filled with fibres, and λ(x) denotes a constant symmetric temperature conductivity tensor with jumps on the interface I between the matrix and the fillers. The symbol \(\vec {n}\) shall denote the unit normal vector on I, always pointing towards one material (i.e. the matrix material), and x+ and x− denote limiting points towards a point x on I approaching x in the direction of \(\vec {n}\) and against the direction of \(\vec {n}\) respectively.

Then, the heat conduction equation and the transition and boundary conditions for the micro scale simulations used in this work for heat transfer in x-direction read:

$$ \begin{array}{@{}rcl@{}} \nabla \cdotp {\left( \lambda(x)\nabla{T(x)}\right)}&=&0, \quad (x,y,z)\in (0,l)^{3} \end{array} $$
(1)
$$ \begin{array}{@{}rcl@{}} T(x^{+})&=&T(x^{-}) \quad on \ I \end{array} $$
(2)
$$ \begin{array}{@{}rcl@{}} (\lambda(x^{+}) \nabla{T(x^{+})}) \cdotp \vec{n}&=& (\lambda(x^{-}) \nabla{T(x^{-})}) \cdotp \vec{n} \quad on \ I \end{array} $$
(3)
$$ \begin{array}{@{}rcl@{}} T(x=0,y,z)&=&T(x=1,y,z)+\delta T \quad (y,z) \in (0,l)^{2} \end{array} $$
(4)
$$ \begin{array}{@{}rcl@{}} T(x,y=0,z)&=&T(x,y=1,z) \quad (x,z) \in (0,l)^{2} \end{array} $$
(5)
$$ \begin{array}{@{}rcl@{}} T(x,y,z=0)&=&T(x,y,z=1) \quad (x,y) \in (0,l)^{2} \end{array} $$
(6)
$$ \begin{array}{@{}rcl@{}} (\lambda \nabla T) \cdot e_{x}|_{x=0}&=& (\lambda \nabla T) \cdot e_{x}|_{x=l} \end{array} $$
(7)
$$ \begin{array}{@{}rcl@{}} (\lambda \nabla T) \cdot e_{y}|_{y=0}&=& (\lambda \nabla T) \cdot e_{y}|_{y=l} \end{array} $$
(8)
$$ \begin{array}{@{}rcl@{}} (\lambda \nabla T) \cdot e_{z}|_{z=0}&=& (\lambda \nabla T) \cdot e_{z}|_{z=l} \end{array} $$
(9)

This set of equations describe a stationary heat conduction with transition conditions between materials respecting the first law of thermodynamics, namely continuity of temperature, energy conservation and periodic boundaries. The equations for the remaining direction are easily obtained by interchanging x, y, z.

From the solution T(x), a homogenized thermal conductivity is calculated as:

$$ \lambda_{i}^{hom}=\frac{\left\langle {\sum}_{j} \lambda_{ij}(x) \frac{\partial T(x)}{\partial x^{j}} \right\rangle}{\delta T / l} $$
(10)

All of the equations in Eq. 9 are linear in λ. Therefore a multiplication of λ with a single constant (the same for all directions and both materials) will not change the solution T(x). In Eq. 10 this will also multiply the result with the same value. For these reasons an increase in matrix conductivity by a factor c yields c times more conductivity than dividing the fiber conductivity by c.

In the same manner the equations are invariant when the length scale is changed, but now the term l in Eq. 10 changes as well, so that the homogenized conductivity is unchanged. Therefore the homogenized conductivity does not really depend on the fiber length, but only the aspect ratio.

Both properties are only strictly valid without particles in the composite. However since the sensitivity to the particle size and conductivity are rather low, as will be demonstrated later in this work, the approximation is viable.

4 Homogenization model

The matrix filler system is motivated from a polyetheretherketone (PEEK) polymer filled with HT short carbon fibres (SCF) to about 7% and up to 10% graphite particles.

The corresponding values are given in Table 2 in Section 5.

Table 2 Parameters used for the sensitivity study
Full size table

In all studies in this section use 5 different random geometry realizations based on the orientation model explained in Section 2. For studies on geometric parameters such as RVE size or fiber length separate geometries had to be created, but for parameters that effect the material properties the same geometries were used. The reported results are averages over 5 realizations of the stochastic fibre distribution, in many cases the standard deviations are displayed by error bars as well (Fig. 1).

Fig. 1
figure 1

RVE of a typical fibre reinforced polymer used in this work at a resolution of 1 μm. (Image generated by GeoDict)

Full size image

4.1 Numerical methods

The major part of the numerical studies was performed with the commercial tool GeoDict (Math2Market GmbH, Kaiserslautern, Germany. www.geodict.de). This tool implements the solution of the thermal diffusion equation for box-shaped domains with very complex interior structure (such as hundreds of cylinders and thousands of spheres surrounded by a matrix material). The method behind the tool resolves the geometry on a regular grid, which is called voxelization in the field; this approach requires a decent resolution to capture the fillers. GeoDict is specifically developed for numerical homogenization based on RVE calculations, and directly calculates the effective thermal conductivity from a given material geometry (Fig. 2).

Fig. 2
figure 2

RVE of a typical fibre reinforced polymer, including hard particles, used in this work at a resolution of 1 μm. (Image generated by GeoDict)

Full size image

Another commonly used method for simulations of thermal conductivities is the finite element method (FEM). Some simulations were performed with the FEM (implemented in the open source tool calculix, http://www.calculix.de [7]) to test aspects that cannot be handled by GeoDict. The FEM has several disadvantages with respect to calculations of RVEs. At first the generation of the mesh itself is rather time-consuming, and also suffers from instabilities of the algorithms used: At the boundary of the RVE, there are many fibres and some particles that are cut by that boundary. At least for the open source mesh generator NetGen(https://sourceforge.net/projects/netgen-mesher) used in this work, these cut objects posed a major problem, so that the boundary of the RVE had to be modified slightly in order to obtain a suitable mesh for the FE solver.

4.2 Numeric validation

As a first step for numerical homogenization an investigation of the dependence on size and resolution for the RVE was performed. For these two studies the maximum particle fraction (10%) and fiber length (280μm) were used, to guarantee numerical convergence for all the intended parameter studies. The remaining parameters were kept at their default values (see Table 2 in Section 5). In the voxelised representations of the geometries used, the volume fraction of the fibres was respected to up ± 0.01% (except for the box size of 250μm, where it was 0.03%) by the voxelization. The particle volume fraction differed slightly more, up to 0.2%. The studies afterwards all have similar uncertainty bounds.

The RVE for homogenization should likely be larger than the fiber length, therefore geometry samples of 2503, 3503, 5003 and 7003 μm were created and used for homogenization at a constant resolution (voxel size= 1μm). The results are shown in Fig. 3.

Fig. 3
figure 3

Influence of RVE size on thermal conductivity for a voxel size of 1μm. For each size 5 independent random realisations of the fibres have been generated. The standard deviation over these realisations is shown as error bars

Full size image

The dependence on the RVE size is rather small, this is likely because both the geometry and the homogenization were performed under complete periodic boundary conditions. Based upon these results a RVE size of 3503μm was considered sufficient for further calculations involving fibres.

As a second step the influence of the resolution was tested, for the now fixed RVE size of 3503μm. For this study only the first 2 geometry samples for the box size of 350μm were used. Different voxel sizes can be applied to one and the same geometry, such that there is no random variation involved in this step. The results for voxel sizes of 2, 1, 0.7 and 0.5 μm are shown in Fig. 4. There is a clear trend visible in these results. From these results a voxels size of 0.5μm is not significantly better than 1μm but it increases the numeric effort by a factor of 8.

Fig. 4
figure 4

Influence of voxel size on thermal conductivity. The first two of the random realisations used in Fig. 3 for a box size of 350μm are used again. In the limit voxel size towards 0, the results are seen to converge linearly. Since the uncertainties in the input parameters (heat conductivities of the materials comprising the composite, resolution of 3D images) are expected to contribute in a similar range to the effective heat conductivity, we use 1 micron as voxel size in subsequent calculations

Full size image

As a result of these studies all further calculations were performed at a box size of 350μm and a voxel size of 1 μm. A smaller voxel size or an extrapolation based upon multiple voxel sizes could be employed to obtain more accurate results, but we recommend to confirm at first that the all the parameters of the geometry and materials are known at a sufficiently high precision so that higher accuracy in discretizations are worth to be used at all.

As a separate test numerical homogenization based upon the FE method was performed for three out of the five geometries generated for the standard case (i.e. without particles and at fiber length 210μm). These FEM calculations used about 9 million 4-node tetraeder elements (with linear shape functions) built from about 1.5 million nodes. For the remaining two geometries the mesh generation failed, despite the modifications to the boundary. The results differ from those obtained from GeoDict by at most 7%. Due to the aforementioned manipulation of the geometry at the boundary and the different numerical solver techniques employed, differences in this range are to be expected.

5 Parameter study

In this section the influence of different input parameters referring either to the assumed geometry or the assumed input material properties is studied.

5.1 Length of the fibres

At first an analysis with respect to the fiber length was performed for four different fiber lengths of 70, 140, 210, and 280 μm. In all cases the fibre diameter was 7μm so the cases correspond to aspect ratios of 10, 20, 30, and 40. The remaining parameters are as given in Table 2.

The results are transferable to fibres of different diameter with the same aspect ratio, by scaling relations (see Section 3).

The corresponding results are shown in Fig. 5. The sensitivity is much higher for short fibres than for long ones.

Fig. 5
figure 5

Composite conductivity for carbon fibres of different length

Full size image

5.2 Volume fraction

The most basic parameter describing a composite is the composition, i.e. the fraction of volume occupied by fibre or particles. The results for fiber volume fractions from 5% to 9% and graphite volume fractions from 0% to 10% are shown in Fig. 6. Again the remaining parameters are kept as in Table 2. The results can be described by a linear model

$$ \begin{array}{@{}rcl@{}} \lambda_{i}^{hom}(v_{f},v_{g})= a_{i}\cdot (1 - v_{f} -v_{g}) + b_{i} \cdot V_{f} + c_{i} \cdot V_{g} \end{array} $$
(11)

quite well as decicted in Fig. 6.

Fig. 6
figure 6

Influence of fibre and graphite volume fraction for several selected values. *** Right Panel: Sensitivity test on variations of λg, consistent with Eq. 11, without significant effects

Full size image

The R-squared (R2) measure for the fit is above 99% and the maximum deviation between the linear fit and the results is below 1.6%. Additionally the offset values (ai) are very close to the matrix value of 0.23 W/mK and the discrepancies are comparable to the discrepancy between 1- and 2-direction. For this reason the linear fit can be expected to be valid without restrictions for smaller volume fractions of fibres as well. This result could be advantageous for material design even though the fit coefficients have been determined by numeric simulations (Table 3).

Table 3 Coefficients for model in Eq. 11
Full size table

5.3 Conductivity of input materials

Another very important parameter is the conductivity of the input materials, especially of the fibre and the matrix material.

In Fig. 7 the conductivity of the composite is shown as a function of the conductivity of the matrix for the base geometry described in Table 2. A variation of the matrix conductivity can be due to different matrix materials, a temperature dependence in the matrix material properties (quite common for typical materials) or additional micro particles that have been homogenized into the matrix in a separate step. There is only a small non-linear effect in the range shown, but for smaller matrix conductivity this can be expected to increase. The coefficients for the linear fits in Fig. 7 are about

$$ \lambda^{hom}_{1,2}=0.12+0.15\lambda^{M} $$
(12)
$$ \lambda^{hom}_{3}=0.5+0.2\lambda^{M} $$
(13)
Fig. 7
figure 7

Thermal conductivity of the composite dependent on the thermal conductivity of the matrix

Full size image

From the same data it is also possible to calculate the effect of the fiber material conductivity by a scaling law (see Eq. 10). For short HT carbon fibers, literature reports conductivity values from 4.9 [26] to 17 W/mK [8], therefore 5 - 20 W/mK seems to be a reasonable range of investigation. The results appear in Fig. 8, again a significant effect is seen.

Fig. 8
figure 8

Influence of fiber conductivity evaluated by scaling law

Full size image

In contrast, the precise value of conductivity for spherical inclusions has only low influence as long as it is significantly larger than the matrix conductivity. This effect is already included in the right part of Fig. 6. A change of thermal conductivity of graphite from 142 W/mK to only 14.2 W/mK lowers the conductivity in 1 and 2 direction only by about 2.5%, whereas a change of volume from 10% to 5% lowers the conductivity in that directions by 12%. For the 3 direction the numbers are lower but of similar ratio.

5.4 Microstructure effects

The thermal conductivity of fibres with higher aspect ratios has a strong influence on the composite properties, but at the same time there exist different values in the literature for supposedly very similar fibres. Given this difficulties in adequately determining the isotropic fiber conductivity, it seems rather unlikely that anisotropic behaviour could be measured reliably. In our numeric studies a reduction of the fibre conductivity in the fibre orthogonal directions by a factor of 10 led to a reduction of composite conductivity by less than 4% (2% in the preferred direction). This calculations were performed by the FEM, because GeoDict cannot deal with anisotropic thermal conductivity. Given the overall accuracy obtained by the FEM for homogenization this effect can be considered insignificant.

Additionally SCF are often coated to improve their contact to the matrix material, but this coating will be missing at the ends if the fibres are cracked during processing. Due to this the fibres may not actually be in contact with the matrix at their ends. This behavior was model by adding strongly isolating semi-spheres (below 1% of matrix conductivity) to the fibre ends. Again the effect on the thermal conductivity of the composite was rather small, below 1.5% in this case.

A further study was performed with the diameter of the graphite particles changed from 9μ m to 18 μ m at constant volume fraction. These effects are not significant as well. The total effect was below 1% (half for preferred fiber direction) and also significantly below the obtained standard deviation (Table 4).

Table 4 Thermal conductivity in 3-direction for different methods (GeoDict, FEM) and different values for the orthogonal thermal conductivity of the fibres
Full size table

5.5 Orientation sensitivity

The thermal conductivity of composites with strong orientation of fibres displays a significant amount of anisotropy. However, over all cases shown in Fig. 9, the isotropic average of the thermal conductivity differs by at most 4%. This can be helpful in comparison to experimental data, because it allows to separate incorrect orientation from other errors. Due to its complex structure, orientation tends to be known at very moderate accuracy.

Fig. 9
figure 9

Influence of fiber orientation on thermal conductivity. The orientations are labelled as in Table 1

Full size image

6 Regression model

The dependence of the thermal conductivity λhom on (i) the volume fractions (ii) fiber length and (ii) fiber orientation can be described by a single relation:

$$ \lambda^{hom}(v_{f},v_{g},a_{ii},l)=a(1-v_{f}-v_{g})+b v_{g}+(c a_{ii} (1-e^{-d l})+e) v_{f} $$
(14)

Here aii is the diagonalized Tucker orientation tensor, derived from the ACG model parameters, in the direction in which the heat transfer is simulated (Fig. 9).

This model actually has even less parameter then Eq. 11, because it handles all directions at once. The fit of the model, see table 5, is demonstrated in Fig. 10. Each individual plot only contains a part of the points such that the variables not on the axis are constant. In total the plots contain each point at least once, and demonstrate that the agreement of the fit is quite good (maximum 4% deviation).

Fig. 10
figure 10

Comparison of fit for all data

Full size image
Table 5 Parameter values for regression model in Eq. 14
Full size table

As a validation of the model two further geometries were evaluated, one with a high volume fraction (9%) of fibers of short length (70μm) and the other with a low volume fraction (5%) of long fibers (280μm). Both of these were created at maximum fiber alignment (configuration B1 from Table 1, i.e. β = 0.1) and evaluated for three different volume fractions of graphite (as in Table 2). For these results the dependence on the graphite volume fraction vg is shown in Fig. 11.

Fig. 11
figure 11

Results of test geometries for their dependence on graphite volume fraction vg

Full size image

For the geometry with fewer but longer fibers a nonlinear dependence is observed, potentially indicating the onset of a percolation effect.

The maximum relative deviation of the regression model Eq. 14 is 14% for the case with many short fibers and 8% for the long fibers. This is significantly more than for the original data used for to build the model (where it was 4%), but since this original data did by design not contain any correlation effects between length and orientation of fibers and the the test cases did deviate from the data in both at once the level of agreement is still acceptable.

7 Conclusion

Complex composite materials and their structure-property-relationships are an important area of research. In this paper, based on established structure modelling techniques, a theoretical sensitivity study of the homogenized heat conduction for double reinforced polymeric composites had been deducted and condensed into such a relationship. It turns out that linear models fit well findings in the composites under investigation. This offers the opportunity to use such structure-property relations in higher scale model approaches, for instance in tribological modeling approaches, in order to describe heat balances.

In this work, we used a ACG type geometrical model of double reinforced composites. Extensions of our approach to more sophisticated distributions with more complicated closure relations in the Folgar-Tucker equations may take advantage of the methods presented. Basically the only ingredient to be modified is an accordingly adapted stochastic fiber geometry generation method reflecting empirical closure relations, see [14]. Thus, the heat conductivity of a structure simulated by CFD-type models, like the ones used in [24], are predictable to a high accuracy using the structure-property relationship developed in this work.