Homogenisation of the Local Thermal Conductivity in Injection-Moulded Short Fibre Reinforced Composites

This paper deals with predicting the effective thermal conductivity (ETC) of injection-moulded short fibre reinforced polymers (SFRPs) using two different homogenisation schemes: a scheme based on the dielectric theory for pseudo-oriented inclusions and a two-step homogenisation model based on the mean-field homogenisation approach. In both cases, the fibre orientation tensor (FOT) obtained from Autodesk Moldflow® simulation was used. The Moldflow FOT predictions were validated via structure tensor analysis of micro-computed X-ray tomography (micro-CT) scans of the part. In the dielectric-wise approach, the orientation of fibres was originally defined by a scalar parameter, which is related to the diagonal components of the FOT. In the two-step homogenisation approach, an interpolative model based on the Mori–Tanaka theory is used in the first step for calculating the ETC for the ideal case of unidirectional fibre alignment, followed by a second step in which orientation averaging based on the FOT inside each element is applied. The ETC was calculated using both schemes for the specific case of uniform fibre orientation distribution and at three different locations with non-identical FOTs of an injection-moulded SFRP part. The results are compared with each other and evaluated against the direct numerical simulation for the uniform fibre orientation and experimental measurements for the injection-moulded SFRP. This shows that while the two-step homogenisation can predict the ETC in the full range of orientations between the perfectly aligned and uniformly distributed fibres, the dielectric-wise approach is only capable of modelling the ETC when distributions are close to the two extreme ends of the orientation spectrum.


Introduction
Thanks to the large-scale manufacturing capacity of injection moulding, there has been dramatic growth in the production and applications of short fibre reinforced polymers (SFRPs) [1]. Within an injection-moulded SFRP part, however, the orientation of fibres may change considerably because of the non-uniform shear stress profile developed during injection [2]. This may result in anisotropic and heterogeneous mechanical, electric and thermal properties at different locations in injection-moulded SFRPs, which emphasises the importance and difficulties of developing efficient modelling approaches to predicting composite behaviour [3]. Thermal conductivity, which is an important parameter in the injection moulding simulations of SFRPs, has been regarded as a scalar constant in the existing software packages. Nevertheless, not only may it vary at different locations within injection-moulded SFRPs but it is a direction-dependent variable [4]. Therefore, thermal conductivity can be regarded as a second-order tensor that is influenced by the shape and orientation of fibres [5]. Since this can affect local cooling rates, part shrinkage and warpage, and hence the mechanical performance of the part, this study aimed to predict the effective thermal conductivity (ETC) by considering the local fibre orientation. To obtain a homogenised property in heterogeneous media, one can choose either asymptotic or analytical-based homogenisation approaches [5]. If the complicated microstructure is provided in detail, asymptotic homogenisation using the finite element method (FEM) may predict the effective properties of composites rather precisely [6][7][8]. However, generating and meshing a representative volume element (RVE) that minutely describes the fibre's placement and orientation in the injection-moulded SFRP is unfeasible. The computational process might be inefficient, time-consuming and costly [9]. Alternatively, analytical-based models can be used. For modelling the homogenised thermal conductivity of heterogeneous materials, two categories can be defined: Models based on the Eshelby theory of inclusions and methods based on Laplace's heat transfer equation [10]. Regarding the latter, which is also called Maxwell's homogenisation scheme [11], early works considered a dilute dispersion of ellipsoidal inclusions embedded in an isotropic matrix [12]. When the volume fraction of inclusions is high, differential methods based on the work of Bruggeman can be introduced to the workflow [13]. Using the Maxwell-Bruggeman solutions, Giordano [14] has proposed explicit formulae to predict the effective permittivity based on the dielectric theory of inclusions for unidirectional and uniform dispersions of ellipsoidal inclusions inside an isotropic matrix. Then, this model has been expanded to include all orientational distributions between the uniform and unidirectional ones by defining a scalar parameter that is dependent on the angular distribution of inclusions [15]. In this paper, the possibility of using Giordano s model for predicting the ETC of injection-moulded SFRPs will be studied.
With regard to the models based on the Eshelby theory of inclusions, which refers to the mean-field homogenisation, a method that has become popular is two-step homogenisation [16]. The procedure for implementing this model starts with decomposing the RVE into sets of pseudo-grains in which the orientations of fibres are identical. Homogenisation takes place within each pseudo-grain in the first step using a mean-field model, followed by homogenisation amongst all pseudo-grains in the second step. The first step in the homogenisation is commonly carried out via Mori-Tanaka (MT) or Lielens double inclusion (DI) models, whereas Voigt or Reuss models are normally used in the second step [17]. For homogenising the thermal conductivity, it has been reported that the MT-Voigt combination will predict the ETC most accurately compared to other possibilities [18]. Although adopting this approach may result in accurate predictions of effective properties, there are two things to keep in mind. First, the RVE decomposition necessitates analysing the orientation distribution function (ODF). Extracting data from the ODF could be time-consuming and tedious. Alternatively, the fibre orientation tensor (FOT) is a more practical orientation descriptor, since it can be derived directly from orientation equations solved in simulation software packages [19]. Moreover, when an orientation tensor is calculated for an assembly with a certain ODF, and then the ODF is reconstructed using the orientation tensor, some information (higher harmonics) is lost and the ODF is changed [20]. The second consideration is the necessity of transforming the homogenised tensorial property from the local coordinate system defined in each pseudo-grain to the global coordinate of the RVE [21].
It is possible to develop a more straightforward two-step homogenisation model which could be equivalent to what has been explained above, yet with fewer complexities. In the first step, a unit cell with a simple arrangement of unidirectional fibres with the same fibre volume fraction as in the original RVE is assumed. The unit-cell can be homogenised via different analytical models or the FEM. In the second step, an orientation averaging rule is used to take the fibre misalignment into account [22]. To the best of our knowledge, the efficiency of this method has already been demonstrated in the homogenisation of elastic properties [23], but not for the ETC prediction. Since the Voigt assumption is preferred in the second step of the ETC homogenisation [18], the very early orientation averaging model proposed by Advani and Tucker [24], which provides orientation averaged properties using the FOT and associated unidirectional properties by assuming the Voigt interaction is used, in this study. The ETC of the composite can thus be predicted without having to worry about the ODF and tensor transformation analyses.
Although various models have been developed to predict the ETC of multi-oriented inclusion-reinforced composites [5], developing an efficient method for predicting the ETC of injection-moulded SFRPs has been less investigated. This is explicitly addressed in this paper by comparing the developed two-step homogenisation scheme to the reformulated Giordano s model [15], as illustrated in Figure 1. Autodesk Moldflow ® is used to find the FOT which is required in both approaches. It has been shown that injection moulding software packages such as Autodesk Moldflow ® and Moldex3D ® can accurately predict the actual fibre orientation distribution in injection-moulded parts [25,26]. However, in this study the FOT obtained from the Moldflow simulations is also evaluated against experimental measurements of the injection-moulded part, performed via micro-computed X-ray tomography (micro-CT) scans. As mentioned previously, fibre orientation in Giordano s model is described by a scalar parameter whose relationship with the components of the FOT has not been stated explicitly. This will be clarified in this research by three basic equations. Furthermore, a two-step homogenisation method using orientation averaging is developed to find the local ETC of injection-moulded SFRPs. In the first step, an interpolative model based on the MT model is used for calculating the ETC in the hypothetical case of perfectly aligned fibres. The results are compared to those obtained from Mori-Tanaka, simplified Giordano, and FEM analyses, the latter which is supposedly the exact solution.
In the second step, the actual misaligned fibre orientation inside each element (FOT on elements) is taken into account using an orientation averaging approach. For validation, the results are evaluated against the experimental measurements and FEM homogenisation for two different SFRP systems. tation averaging model proposed by Advani and Tucker [24], which provides orientation averaged properties using the FOT and associated unidirectional properties by assuming the Voigt interaction is used, in this study. The ETC of the composite can thus be predicted without having to worry about the ODF and tensor transformation analyses.
Although various models have been developed to predict the ETC of multi-oriented inclusion-reinforced composites [5], developing an efficient method for predicting the ETC of injection-moulded SFRPs has been less investigated. This is explicitly addressed in this paper by comparing the developed two-step homogenisation scheme to the reformulated Giordano´s model [15], as illustrated in Figure 1. Autodesk Moldflow ® is used to find the FOT which is required in both approaches. It has been shown that injection moulding software packages such as Autodesk Moldflow ® and Moldex3D ® can accurately predict the actual fibre orientation distribution in injection-moulded parts [25,26]. However, in this study the FOT obtained from the Moldflow simulations is also evaluated against experimental measurements of the injection-moulded part, performed via microcomputed X-ray tomography (micro-CT) scans. As mentioned previously, fibre orientation in Giordano´s model is described by a scalar parameter whose relationship with the components of the FOT has not been stated explicitly. This will be clarified in this research by three basic equations. Furthermore, a two-step homogenisation method using orientation averaging is developed to find the local ETC of injection-moulded SFRPs. In the first step, an interpolative model based on the MT model is used for calculating the ETC in the hypothetical case of perfectly aligned fibres. The results are compared to those obtained from Mori-Tanaka, simplified Giordano, and FEM analyses, the latter which is supposedly the exact solution. In the second step, the actual misaligned fibre orientation inside each element (FOT on elements) is taken into account using an orientation averaging approach. For validation, the results are evaluated against the experimental measurements and FEM homogenisation for two different SFRP systems.   Within the model proposed by Giordano [15], which was originally developed for predicting effective electric permittivity, Maxwell s solution is expanded based on Bruggeman's differential approach to consider the orientational order/disorder of ellipsoidal particles in a composite by introducing a scalar parameter (S). The main results which have been reformulated to calculate the longitudinal (k xx ) and transversal (k zz ) ETC of SFRPs are

Homogenisation Models
. ( In Equations (1) and (2), ν f , k m and k f are the volume fraction of fibres, the thermal conductivity of the matrix and thermal conductivity of fibres, respectively. L is called the depolarisation factor, which is dependent on the aspect ratio of fibres (fibre length divided by its diameter, e = l/d), which can be computed from The depolarisation factor is calculated for an ellipsoidal inclusion with a x > a y = a z , where a x , a y and a z are semi-axes of the ellipsoid. Another assumption in developing the S parameter in this model is that the thermal conductivity of the composite is assumed to be transversely isotropic and can be represented by a tensor k: Concerning the orientational factor S, it can be shown that it is related to the diagonal components of the FOT (a ij ), as in When S = 1, the state of complete order is considered; i.e., all fibres are assumed to be aligned with x-direction (a xx = 1). In this case, Equations (1) and (2) will be simplified to calculate the ETC in the case of complete order (unidirectional fibre alignment): The case S = 0 refers to a random orientation in which fibres have no preferential direction (a xx = a yy = a zz = 1/3). The deviation from the complete ordered condition increases as S takes values from 1 to 0.
It should be noted that assuming the transverse isotropic behaviour requires two conditions to be satisfied in the FOT; a yy = a zz and a ij = 0 when i = j; that is, the offdiagonal components of the orientation tensor are zero. This leads to the major limitation of Giordano s model: Some particular distributions of fibres were considered, which are not exhaustive of all possible existing statistical distributions. In other words, there can be some distributions of inclusions, or equivalently FOT on some elements that cannot be represented by the parameter S. In this paper, we will investigate whether this limitation can cause major restrictions in predicting the local ETC of injection-moulded SFRPs.

The Mori-Tanaka (MT) Model and Lielens Interpolation
A family of models for predicting the effective properties of composite materials with the no-dilute concentration of inclusions has been developed from a proposal originally made by Mori and Tanaka [27], which is here reformulated for predicting the ETC of SFRPs.
The ETC of the composite can be expressed as a second-order tensor, k, which can be calculated from where A is the intensity-concentration tensor (similar to the strain-concentration tensor in elasticity), which relates the average temperature gradient (or strain) of fibres to that of the composite. An alternative concentration tensor,Â, is defined so as to link the average temperature gradient (or strain) of fibres to that of the matrix. The relationship between A andÂ is where I is the second-order identity tensor. Different models in the mean-field homogenisation scheme have suggested various formulae for obtaining the concentration tensors. Based on the MT, which is equivalent to the lower bound solution, we have: S is called the interior-point Eshelby thermal conductance tensor, which is a diagonal tensor depending only on the shape of inclusions. However, for the case of an ovary or elongated ellipsoidal inclusion with semi-axes a x > a y = a z , the diagonal components of S can be computed from [28]: The inverse Mori-Tanaka (IMT) formulation, which is identical to the upper bound solution, can be derived by assuming that elliptical inclusions of the polymer matrix are embedded in a continuous medium of the fibre material, resulting in a concentration tensor: As the MT prediction of the effective properties (the lower bound) may deviate from true values at high fibre volume fractions, Lielens et al. [29] proposed a formula to interpolate between the upper and lower bounds: f is the interpolating factor which can be obtained from the fibre volume fraction, as shown in Equation (15). As a result, the prediction of the effective properties is improved at higher fibre contents.

FEM Analysis: Unidirectional Fibre Alignment
Experimental validation of the predicted ETC in the case of complete order is problematic because fabricating samples with perfectly aligned short fibres is unfeasible. Instead, the FEM can be used for calculating the supposedly exact solution for this case. Based on [27], staggered arrays with square packing were chosen as the repetitive volume element, which is illustrated in Figure 2. Dimensions of the unit cell (a and c) are calculated using the volume fraction of the inclusions (v f ) and the fibre s length (l) and diameter (d). It should be noted that as long as we are dealing with unidirectional short fibres, increasing the size of the RVE to contain more fibres is useless, as it has been shown that increasing the size of the RVE, in this case, has almost no effect on the convergence of the values associated with the studied homogenised property [23].
Polymers 2022, 14, x FOR PEER REVIEW 6 of 17 element, which is illustrated in Figure 2. Dimensions of the unit cell ( and ) are calculated using the volume fraction of the inclusions ( ) and the fibre´s length ( ) and diameter ( ). It should be noted that as long as we are dealing with unidirectional short fibres, increasing the size of the RVE to contain more fibres is useless, as it has been shown that increasing the size of the RVE, in this case, has almost no effect on the convergence of the values associated with the studied homogenised property [23]. For the ETC calculation, the steady-state heat transfer equation is considered. Boundary conditions are imposed on the RVE as follows: For obtaining the ETC along each direction (conduction axes), a fixed temperature constraint is imposed on the faces of the RVE, whose normal vector is parallel to that direction. The other four faces of the unit cell are insulated; i.e., zero heat flux condition is assumed on the faces orthogonal to the conduction axis. The FEM solution provides the average heat flux along the conduction axis ( ). Knowing the imposed temperature gradient (∇ ) and on the basis of Fourier´s constitutive equation, the ETC along each axis can be computed using = | ⁄ |.

Orientation Averaging
In order to incorporate the effect of fibre orientation on the ETC, consider ( ) as a tensorial property of the composite associated with a unidirectional microstructure aligned in the direction of , which is a unit vector along each fibre´s axis. The orientation average of is denoted by 〈 〉 and is defined by where ( ) is the ODF. must be a transversely-isotropic tensor with as its axis of symmetry. Assume ( ) to be a second-order tensor, so that it can be expressed as: On the other hand, this tensor must also have the form: where and are two scalar constants, and is the Kronecker delta. If we take the orientation average of the property, 〈 〉 , we get: For the ETC calculation, the steady-state heat transfer equation is considered. Boundary conditions are imposed on the RVE as follows: For obtaining the ETC along each direction (conduction axes), a fixed temperature constraint is imposed on the faces of the RVE, whose normal vector is parallel to that direction. The other four faces of the unit cell are insulated; i.e., zero heat flux condition is assumed on the faces orthogonal to the conduction axis. The FEM solution provides the average heat flux along the conduction axis (Q i ). Knowing the imposed temperature gradient (∇T i ) and on the basis of Fourier s constitutive equation, the ETC along each axis can be computed using k i = |Q i /∇T i |.

Orientation Averaging
In order to incorporate the effect of fibre orientation on the ETC, consider k(p) as a tensorial property of the composite associated with a unidirectional microstructure aligned in the direction of p, which is a unit vector along each fibre s axis. The orientation average of k is denoted by k and is defined by where ψ(p) is the ODF. k must be a transversely-isotropic tensor with p as its axis of symmetry. Assume k(p) to be a second-order tensor, so that it can be expressed as: On the other hand, this tensor must also have the form: where c and d are two scalar constants, and δ ij is the Kronecker delta. If we take the orientation average of the property, k ij , we get: In other words, the orientation average of a second-order tensor is completely determined by the second-order orientation tensor, and by the underlying unidirectional property tensor [24]. Now we aim to find the relationships between the constants c and d, and unidirectional properties (k x and k z ). For this purpose, the first and second tensor invariants (I 1 and I 2 ) of two expressions of k based on Equations (17) and (18) are used to form a system of equations: Solving the system of equations results in: Substituting Equation (22) in Equation (19) results in Equation (23), in which the orientation average of the second-order tensorial property of SFRPs (for example, the ETC) can be calculated using the FOT and properties associated with unidirectional fibre alignment:

Autodesk Moldflow Simulation
The injection moulding process was simulated via Autodesk Moldflow ® 2019 using 3D meshing. Fibre orientation tensors on elements were exported from the software to be used in homogenisation models. The FOT on elements was calculated using the Moldflow rotational diffusion model with auto-calculated fibre interaction coefficient and coefficients of asymmetry. It has been reported that the simulated FOT is dependent on the mesh size, and by mesh refinement, the simulated FOT may converge to the experimentally measured FOT [30]. In our simulations, the mesh size gradually decreased, until the mesh sensitivity became negligible.

Case study 1: Uniform Distribution of Fibres
The FEM-homogenised ETC of a composite reinforced with 10 volume percent randomly distributed cylindrical fibres with an aspect ratio of 10 (l = 50 µm and d = 5 µm) was taken from [18]. The matrix and fibres were assumed to be isotropic with thermal conductivities of 72.0 and 32.0 Wm −1 K −1 , respectively. In order to predict the ETC of the composite with a uniform distribution of fibres via FEM homogenisation, the authors of [18] generated an RVE with the size of 100 µm 3 using the random sequential absorption (RSA) algorithm. After imposing periodic boundary conditions and solving the Fourier s constitutive equation numerically, they reported the ETC of 67.4 Wm −1 K −1 when the perfect interface was assumed between fibres and the matrix. This case was modelled using our two-step homogenisation mode and Giordano's approach, and the results are compared to those of FEM homogenisation.

Case Study 2: Fibre Orientation of Injection-Moulded SFRPs
DOMAMID 6LVG50H2BK, a commercial grade polyamide 6 (PA6) reinforced with 50 weight-percent short glass fibre (which is equivalent to around 30.8 volume percent), was injection-moulded to produce dog-bone shaped samples originally intended for tensile testing based on the ISO 527 1A standard. The part geometry and the injection moulding process parameters are represented in the Supplementary Material ( Figure S1 and Table S1).
Thermal conductivity of unfilled PA6 was measured at 0.26 Wm −1 K −1 , and the conductivity of glass fibres was assumed to be 1.04 Wm −1 K −1 [31]. The average length and diameter of glass fibres after injection moulding were found to be 260.8 and 11.4 µm, respectively, via micro-CT scans; Therefore, an average aspect ratio of around 23 can be considered for the glass fibres [32]. The 4 mm-thick tensile sample was injection-moulded via two gates located at two ends of the dog-bone sample in order to produce a weld-line at the middle of the part. This was done to make three zones in the injection-moulded sample in which the simulated FOTs would differ noticeably from each other, as illustrated in Figure 3. These regions are listed as A, B and C in Figure 3 corresponding to the wider section, gauge length and weld-line regions of the specimen, respectively. The detailed FOT analysis at these positions is discussed in Section 3.4.1. The flow direction was along the x-axis, and the thickness of the part was aligned with the z-axis. The local ETC was calculated using the two homogenisation approaches at three regions of the injection-moulded part and was evaluated against experimental measurements for assessment of the models. diameter of glass fibres after injection moulding were found to be 260.8 and 11.4 μm, respectively, via micro-CT scans; Therefore, an average aspect ratio of around 23 can be considered for the glass fibres [32]. The 4 mm-thick tensile sample was injection-moulded via two gates located at two ends of the dog-bone sample in order to produce a weld-line at the middle of the part. This was done to make three zones in the injection-moulded sample in which the simulated FOTs would differ noticeably from each other, as illustrated in Figure 3. These regions are listed as A, B and C in Figure 3 corresponding to the wider section, gauge length and weld-line regions of the specimen, respectively. The detailed FOT analysis at these positions is discussed in Section 3.4.1. The flow direction was along the x-axis, and the thickness of the part was aligned with the z-axis. The local ETC was calculated using the two homogenisation approaches at three regions of the injectionmoulded part and was evaluated against experimental measurements for assessment of the models.

Micro-Computed X-Ray Tomography
Micro-CT scanning acquisition was performed using a GE Phoenix Nanotom system which was equipped with a 180 kV/15 W nanofocus X-ray tube and a 2300 × 2300 pixel detector on a 12-bit Hamamatsu flat panel. The X-ray tube's accelerating voltage and the beam current were 100 kV and 135 μA, respectively. The installed target material consisted of tungsten on CVD synthetic diamond. There were no additional filter materials applied. The detector integration time was set to 0.5 s, and the number of projections was 3500. The specimen for micro-CT imaging was cut from the gauge length of the injectionmoulded tensile bar, as illustrated in Figure 3, with dimensions of 4 × 4 × 4 mm 3 . The scans have a resolution of 2.5 μm per pixel, and the average glass fibre diameter was 11.4 μm. Reconstruction of scans was performed in Datos|x software, and micro-CT images were imported into the VoxTex to extract the fibre orientation tensor in the scanned region. The structural tensor analysis was then used to assign local fibre direction data to the components of the voxel model [33]. The process for creating voxel models from the X-ray computed tomography data was thoroughly explained in [34]. Further details for the fibre orientation analysis in VoxTex are presented in Section 3.4.1.

Experimental Measurement of the ETC
A contact-based transient method using the C-therm Trident instrument was selected to measure the thermal conductivity of injection-moulded samples along different directions at room temperature. The Modified Transient Plane Source (MTPS) sensor was used to carry out the measurements according to ASTM D7984. The single-sided MTPS sensor was equipped with a guard ring to facilitate one-dimensional heat transfer inside the sample, which is required for the thermal conductivity measurement of anisotropic materials. When the density ( ) and specific heat ( ) of the test sample are known, the thermal

Micro-Computed X-ray Tomography
Micro-CT scanning acquisition was performed using a GE Phoenix Nanotom system which was equipped with a 180 kV/15 W nanofocus X-ray tube and a 2300 × 2300 pixel detector on a 12-bit Hamamatsu flat panel. The X-ray tube's accelerating voltage and the beam current were 100 kV and 135 µA, respectively. The installed target material consisted of tungsten on CVD synthetic diamond. There were no additional filter materials applied. The detector integration time was set to 0.5 s, and the number of projections was 3500. The specimen for micro-CT imaging was cut from the gauge length of the injection-moulded tensile bar, as illustrated in Figure 3, with dimensions of 4 × 4 × 4 mm 3 . The scans have a resolution of 2.5 µm per pixel, and the average glass fibre diameter was 11.4 µm. Reconstruction of scans was performed in Datos|x software, and micro-CT images were imported into the VoxTex to extract the fibre orientation tensor in the scanned region. The structural tensor analysis was then used to assign local fibre direction data to the components of the voxel model [33]. The process for creating voxel models from the X-ray computed tomography data was thoroughly explained in [34]. Further details for the fibre orientation analysis in VoxTex are presented in Section 3.4.1.

Experimental Measurement of the ETC
A contact-based transient method using the C-therm Trident instrument was selected to measure the thermal conductivity of injection-moulded samples along different directions at room temperature. The Modified Transient Plane Source (MTPS) sensor was used to carry out the measurements according to ASTM D7984. The single-sided MTPS sensor was equipped with a guard ring to facilitate one-dimensional heat transfer inside the sample, which is required for the thermal conductivity measurement of anisotropic materials. When the density (ρ) and specific heat (C p ) of the test sample are known, the thermal conductivity can be calculated via Equation (24) using thermal effusivity (ε), which is directly measured by the MTPS sensor: k = ε 2 /(ρC p ). The specific heat of the material was determined at 1140 Jkg −1 K −1 at room temperature using a Q200 TA instrument differential scanning calorimeter (DSC), and the density was taken from the material supplier at 1560 kgm −3 .
As the MTPS sensor requires specimens with a minimum diameter of 18 mm to cover the sensor's surface, samples with dimensions of 18 × 10 × 4 mm 3 were milled at three different locations of the injection-moulded sample shown in Figure 3. For each region, 10 samples were polished and then glued to each other using a thin layer of epoxy paste in order to form a block of 18 × 20 × 20 mm 3 , as illustrated in Figure 4. Two blocks were made for each region. This allowed us to measure thermal conductivities along each axis by rotating this block on the sensor. For instance, if the x-y face of the block was in contact with the sensor, thermal conductivity along the z-axis (thickness) was measured. Three measurements were performed at room temperature for each face of the block, and the average value is reported as the ETC along each direction. Glycol was used as the contact agent between the sample and the sensor.
conductivity can be calculated via Equation (24) using thermal effusivity (ε), which is directly measured by the MTPS sensor: The specific heat of the material was determined at 1140 Jkg −1 K −1 at room temperature using a Q200 TA instrument differential scanning calorimeter (DSC), and the density was taken from the material supplier at 1560 kgm −3 .
As the MTPS sensor requires specimens with a minimum diameter of 18 mm to cover the sensor's surface, samples with dimensions of 18 × 10 × 4 mm 3 were milled at three different locations of the injection-moulded sample shown in Figure 3. For each region, 10 samples were polished and then glued to each other using a thin layer of epoxy paste in order to form a block of 18 × 20 × 20 mm 3 , as illustrated in Figure 4. Two blocks were made for each region. This allowed us to measure thermal conductivities along each axis by rotating this block on the sensor. For instance, if the x-y face of the block was in contact with the sensor, thermal conductivity along the z-axis (thickness) was measured. Three measurements were performed at room temperature for each face of the block, and the average value is reported as the ETC along each direction. Glycol was used as the contact agent between the sample and the sensor.

Calculation of the Depolarisation Factor and Eshelby Tensor for the Cases Studies
Looking at Equations (3) and (12) shows that the depolarisation factor is identical to the second (or third) diagonal component of the Eshelby tensor. This can be proven by assuming = ⁄ and some other mathematical simplifications. Therefore: Originally, the depolarisation factor used in Equations (1) and (2) is a simplified notation of three depolarisation factors, , and , defined along each semi-axis of the ellipsoidal inclusion. Derivation of these depolarisation factors is explicitly expressed in terms of elliptic integrals in the reference [14], where it has been shown that + + = 1. However, for an ellipsoid of rotation in which = , we have = = and = 1 − 2 . This is exactly analogous to the concept of the Eshelby tensor, in which for the same inclusion we have = and = 1 − 2 . This means that the depolarisation vector corresponds to the diagonal Eshelby tensor, both having two independent

Calculation of the Depolarisation Factor and Eshelby Tensor for the Cases Studies
Looking at Equations (3) and (12) shows that the depolarisation factor is identical to the second (or third) diagonal component of the Eshelby tensor. This can be proven by assuming e = a x /a z and some other mathematical simplifications. Therefore: Originally, the depolarisation factor L used in Equations (1) and (2) is a simplified notation of three depolarisation factors, L x , L y and L z , defined along each semi-axis of the ellipsoidal inclusion. Derivation of these depolarisation factors is explicitly expressed in terms of elliptic integrals in the reference [14], where it has been shown that L x + L y + L z = 1. However, for an ellipsoid of rotation in which a y = a z , we have L y = L z = L and L x = 1 − 2L. This is exactly analogous to the concept of the Eshelby tensor, in which for the same inclusion we have S yy = S zz and S xx = 1 − 2S yy . This means that the depolarisation vector corresponds to the diagonal Eshelby tensor, both having two independent variables. Therefore, the construction of the Eshelby tensor for the electric/thermal problem is completely equivalent to the formulation based on the depolarisation factors, which was developed quite earlier in the literature [35]. The formalism based on the Eshelby tensor is simply more modern and useful for some generalisations, for instance, when the matrix properties are anisotropic. It is also interesting to notice that in Equation (25) lim e→∞ L = lim e→∞ S 22 = 1/2; that is, for fibres with high aspect ratios, we have: To understand what a "high aspect ratio" implies, Figure 5 shows how L or S yy changes with different values of aspect ratio (e). It can be seen that the curve reaches its limit at 0.5 rather quickly. For example, for the composites of cases 1 and 2 with aspect ratios of 10 and 23, L = S yy = 0.490 and 0.497, respectively. As a result, in both cases, Equation (26) was a reasonable approximation for further calculations. variables. Therefore, the construction of the Eshelby tensor for the electric/thermal problem is completely equivalent to the formulation based on the depolarisation factors, which was developed quite earlier in the literature [35]. The formalism based on the Eshelby tensor is simply more modern and useful for some generalisations, for instance, when the matrix properties are anisotropic. It is also interesting to notice that in Equation (25) To understand what a "high aspect ratio" implies, Figure 5 shows how or changes with different values of aspect ratio ( ). It can be seen that the curve reaches its limit at 0.5 rather quickly. For example, for the composites of cases 1 and 2 with aspect ratios of 10 and 23, = = 0.490 and 0.497, respectively. As a result, in both cases, Equation (26) was a reasonable approximation for further calculations.

ETC Prediction for Unidirectional Fibre Alignment
In the first phase of the two-step homogenisation approach, the ETC of the aforementioned case studies for unidirectional fibre alignment was obtained using Giordano's model (Equations (6) and (7)), the Mori-Tanaka model, and the Lielens interpolation, as described in Section 2.1.2. Unidirectional fibre reinforced composites are transversely isotropic; so the ETC is defined by longitudinal and transversal components, denoted by and respectively. In Table 1, the results are compared with each other and evaluated against the FEM homogenisation. Table 1 indicates that there is no difference in predicted based on three models, which is also consistent with the FEM results. On the other hand, the ETCs perpendicular to the fibre´s direction are not identical.
In case study 2, where the composite contained a relatively high fraction of fibres, Giordano and Lielens predictions turned out to be the most compatible with the FEM homogenisation, whereas all models produced similar results in case study 1 with lower fibre content. and predicted by three models are plotted versus the fibre volume fraction for PA6 composites reinforced with unidirectional short glass fibres (aspect ratio of around 23) in Figure 6. As can be seen, while all models predict the linear fluctuations of in accordance with the Voigt upper bound approach, this is not the case for over the whole range of . For low concentrations of fibres, all models provide the same predictions; but by increasing , they start to deviate from each other. Nevertheless, good agreement is shown in transversal conductivities obtained from the Giordano and Lielens

ETC Prediction for Unidirectional Fibre Alignment
In the first phase of the two-step homogenisation approach, the ETC of the aforementioned case studies for unidirectional fibre alignment was obtained using Giordano's model (Equations (6) and (7)), the Mori-Tanaka model, and the Lielens interpolation, as described in Section 2.1.2. Unidirectional fibre reinforced composites are transversely isotropic; so the ETC is defined by longitudinal and transversal components, denoted by k x and k z respectively. In Table 1, the results are compared with each other and evaluated against the FEM homogenisation. Table 1 indicates that there is no difference in predicted k x based on three models, which is also consistent with the FEM results. On the other hand, the ETCs perpendicular to the fibre s direction are not identical. homogenisation, whereas all models produced similar results in case study 1 with lower fibre content. k x and k z predicted by three models are plotted versus the fibre volume fraction for PA6 composites reinforced with unidirectional short glass fibres (aspect ratio of around 23) in Figure 6. As can be seen, while all models predict the linear fluctuations of k x in accordance with the Voigt upper bound approach, this is not the case for k z over the whole range of v f . For low concentrations of fibres, all models provide the same predictions; but by increasing v f , they start to deviate from each other. Nevertheless, good agreement is shown in transversal conductivities obtained from the Giordano and Lielens models up to v f = 0.5. For consistency, the data from Lielens interpolative model were hence taken to be used in the second step of the two-step homogenisation method.
Polymers 2022, 14, x FOR PEER REVIEW 11 of 17 models up to = 0.5. For consistency, the data from Lielens interpolative model were hence taken to be used in the second step of the two-step homogenisation method.

ETC Prediction for the Uniform Distribution of Fibres: Case Study 1
In this section, the FEM homogenisation result of case study 1 taken from [18] is compared to that obtained from the two-step homogenisation and Giordano´s model.
Solving Equation (27) by substituting the composite properties of case study 1 and = 0.49, calculated in Section 3.1, results in = 67.09 Wm −1 K −1 . Meanwhile, following the two-step homogenisation approach for this case, based on Equation (23) and Table 1, leads to an identical outcome of = 67.10 Wm −1 K −1 . Additionally, there is good agreement between these results and those of FEM-homogenised ET reported in [18]: = 67.40 Wm −1 K −1 . Therefore, it can be concluded that the two-step homogenisation method and the reformulated Giordano model could accurately predict the ETCs of composites in cases of unidirectional and uniformly distributed fibres.

ETC Prediction for the Uniform Distribution of Fibres: Case Study 1
In this section, the FEM homogenisation result of case study 1 taken from [18] is compared to that obtained from the two-step homogenisation and Giordano s model.  (1) and (2) results in a simplified Giordano s formulation where k xx = k zz = k: Solving Equation (27) by substituting the composite properties of case study 1 and L = 0.49, calculated in Section 3.1, results in k = 67.09 Wm −1 K −1 . Meanwhile, following the two-step homogenisation approach for this case, based on Equation (23) and Table 1, leads to an identical outcome of k = 67.10 Wm −1 K −1 . Additionally, there is good agreement between these results and those of FEM-homogenised ET reported in [18]: k = 67.40 Wm −1 K −1 . Therefore, it can be concluded that the two-step homogenisation method and the reformulated Giordano model could accurately predict the ETCs of composites in cases of unidirectional and uniformly distributed fibres.

ETC Prediction for Misaligned Fibres in Injection-Moulded SFRPs: Case Study 2 3.4.1. Simulated FOT Predictions and Validation
In order to predict the ETCs of injection-moulded SFRPs, the misaligned fibre orientations were scrutinised at three zones of an injection-moulded tensile bar via Autodesk Moldflow ® simulations, as illustrated in Figure 3. First, the FOTs on elements associated with each region were extracted and analysed to explain the differences in fibre orientation at these positions. Region A, which was located after the gate in the wider section of the tensile bar, showed a fairly misaligned in-plane orientation, a typical fibre orientation distribution in injection-moulded parts. While the average a xx in this region was around 0.63, a significant proportion of fibres were still aligned with the y-axis, the average of a yy being 0.35. As expected, a very small number of fibres tended to be out of plane, along z-direction in this region.
In contrast, due to converging flows in the gauge length of the sample, strong fibre alignment along the flow direction was seen in region B (average a xx = 0.84). Moreover, a small number of fibres were aligned along other two axes. Unlike these two samples, a considerable number of fibres were oriented out of the flow plane in region C because of the presence of a weld-line: on average, a zz = 0.20. The variation in a xx throughout the thickness of the injection-moulded part is shown in Figure 7. Despite the fluctuations in a xx throughout the thickness in region A, which are consistent with the skin-shell-core structure induced by the fountain flow in injection moulding [2], a xx in the gauge length of the samples remained almost unchanged due to the pronounced extensional flow.

Simulated FOT Predictions and Validation
In order to predict the ETCs of injection-moulded SFRPs, the misaligned fibre orientations were scrutinised at three zones of an injection-moulded tensile bar via Autodesk Moldflow ® simulations, as illustrated in Figure 3. First, the FOTs on elements associated with each region were extracted and analysed to explain the differences in fibre orientation at these positions. Region A, which was located after the gate in the wider section of the tensile bar, showed a fairly misaligned in-plane orientation, a typical fibre orientation distribution in injection-moulded parts. While the average in this region was around 0.63, a significant proportion of fibres were still aligned with the y-axis, the average of being 0.35. As expected, a very small number of fibres tended to be out of plane, along zdirection in this region.
In contrast, due to converging flows in the gauge length of the sample, strong fibre alignment along the flow direction was seen in region B (average = 0.84). Moreover, a small number of fibres were aligned along other two axes. Unlike these two samples, a considerable number of fibres were oriented out of the flow plane in region C because of the presence of a weld-line: on average, = 0.20. The variation in throughout the thickness of the injection-moulded part is shown in Figure 7. Despite the fluctuations in throughout the thickness in region A, which are consistent with the skin-shell-core structure induced by the fountain flow in injection moulding [2], in the gauge length of the samples remained almost unchanged due to the pronounced extensional flow. For validating the simulated FOT results generated by Autodesk Moldflow ® , experimentally measured FOT was obtained in the gauge length of the sample (outside of the weld-line region) through processing the micro-CT scans in the VoxTex software [34]. The micro-CT scans were imported to the software in the form of an image stack to set up a 3D model, the so-called voxel model. This is a 3D mesh of volume elements with the orientation vector and additional variables attached to each element. The next step was preprocessing the scans, which was performed by expanding the dynamic grey-scale range of the images and increasing the contrast by selecting the main part of the histogram. Then, a region of interest (ROI) with dimensions of 3125 × 3125 × 3125 μm 3 was cut from the original micro-CT scans in the VoxTex.
For generating the voxel model, two important parameters should be defined. The first one is the mesh density, which determines the size of the voxel model. Our voxel model was produced by setting the mesh density to 20 pixels, resulting in 226,981 voxels While the skin-shell-core structure is seen outside of the gauge length (region A), the fibre orientation throughout the thickness inside the gauge length is almost constant (regions B and C).
For validating the simulated FOT results generated by Autodesk Moldflow ® , experimentally measured FOT was obtained in the gauge length of the sample (outside of the weld-line region) through processing the micro-CT scans in the VoxTex software [34]. The micro-CT scans were imported to the software in the form of an image stack to set up a 3D model, the so-called voxel model. This is a 3D mesh of volume elements with the orientation vector and additional variables attached to each element. The next step was pre-processing the scans, which was performed by expanding the dynamic grey-scale range of the images and increasing the contrast by selecting the main part of the histogram. Then, a region of interest (ROI) with dimensions of 3125 × 3125 × 3125 µm 3 was cut from the original micro-CT scans in the VoxTex.
For generating the voxel model, two important parameters should be defined. The first one is the mesh density, which determines the size of the voxel model. Our voxel model was produced by setting the mesh density to 20 pixels, resulting in 226,981 voxels in the ROI. The second parameter is the window radius, which determines the computation accuracy of the voxel model. The general rule is that the larger the diameter of the fibre in the image (in pixels), the larger the window size should be set. For our case, the window radius was set to 8.
Each voxel in the model has data about the average grey scale (denotes material density), degree of anisotropy and its principal direction. When performing fibre orientation analysis in VoxTex, voxels containing exclusively the matrix material must be excluded. For this purpose, an anisotropy histogram of the voxel model should be plotted. Normally, this histogram has two peaks. The peak with a small degree of anisotropy represents the matrix, and the one with higher values corresponds to the fibres [33]. However, in our case only one peak was seen, which is related to the fibres. This means that there was no full-matrix voxel in our model. This could be due to the high fibre content (50 wt%) in the composite. As a result, no anisotropy filter is needed to be applied in the fibre orientation calculation. The principles of structure tensor analysis, which was implemented in VoxTex, for calculating the orientation distribution function and fibre orientation tensor, are explained in detail in Section 2.3 of [33]. Figure 8 shows examples of micro-CT scans in different planes, which were taken of the centre of the ROI, and the reconstructed 3D image. It is evident that the majority of fibres were oriented along the flow direction (x-axis), as predicted by the Moldflow. For quantitative analysis, after generating the voxel model in VoxTex, the diagonal components of the FOT calculated by the VoxTex were compared to those obtained from the Moldflow simulation in the same ROI, and are represented in Table 2. Overall, the Moldflow FOT predictions were found to be highly consistent with the experimentally measured FOTs. These findings, which are in accordance with the previous studies [25,26], demonstrate that the Moldflow fibre orientation predictions can be trusted for finding the homogenised ETC.
fibre in the image (in pixels), the larger the window size should be set. For our case, the window radius was set to 8.
Each voxel in the model has data about the average grey scale (denotes material density), degree of anisotropy and its principal direction. When performing fibre orientation analysis in VoxTex, voxels containing exclusively the matrix material must be excluded. For this purpose, an anisotropy histogram of the voxel model should be plotted. Normally, this histogram has two peaks. The peak with a small degree of anisotropy represents the matrix, and the one with higher values corresponds to the fibres [33]. However, in our case only one peak was seen, which is related to the fibres. This means that there was no full-matrix voxel in our model. This could be due to the high fibre content (50 wt%) in the composite. As a result, no anisotropy filter is needed to be applied in the fibre orientation calculation. The principles of structure tensor analysis, which was implemented in VoxTex, for calculating the orientation distribution function and fibre orientation tensor, are explained in detail in Section 2.3 of [33]. Figure 8 shows examples of micro-CT scans in different planes, which were taken of the centre of the ROI, and the reconstructed 3D image. It is evident that the majority of fibres were oriented along the flow direction (x-axis), as predicted by the Moldflow. For quantitative analysis, after generating the voxel model in VoxTex, the diagonal components of the FOT calculated by the VoxTex were compared to those obtained from the Moldflow simulation in the same ROI, and are represented in Table 2. Overall, the Moldflow FOT predictions were found to be highly consistent with the experimentally measured FOTs. These findings, which are in accordance with the previous studies [25,26], demonstrate that the Moldflow fibre orientation predictions can be trusted for finding the homogenised ETC.

Homogenisation of Thermal Conductivity
Regarding the use of Giordano s model to predict the ETC, a MATLAB script was written to select elements that comply with the transverse isotropy conditions. It is seen that not all elements associated with each zone will satisfy these conditions. The percentage of a number of elements satisfying the transverse isotropy conditions in each region can be used as a criterion for the validity of using Giordano s model for predicting the ETC in that region of the part. This is indicated by %Ele. in Table 3. After calculating the parameter S for each valid element, the longitudinal and transversal ETC are computed for each element using Equations (1) and (2), and then can be volume averaged for the region of interest. In the two-step homogenisation approach, firstly k x = 0.498 and k z = 0.386 Wm −1 K −1 were used in Equation (23) based on the Lielens prediction for unidirectional fibre alignment, according to Table 1. Afterwards, knowing the FOT for each element from the Moldflow analysis and using Equation (23), the orientation averaged ETC was computed along each direction for every element and then volume-averaged for the region. These data are reported in Table 3 and compared to the measured ETC obtained from experiments. First of all, Giordano's predictions can be implicitly trusted when the majority of elements associated with a region of interest satisfy the transverse isotropy conditions. This can be seen only in region B, since about two-third of elements comply with those conditions. Giordano's predictions in this region agree well with that of the two-step homogenisation approach and the experimental measurements. On the contrary, in region A, no element satisfies the conditions. Hence, Giordano's model failed to predict the ETC in this zone. Taking the values of diagonal components of the FOT in those regions into account, it can be deduced that Giordano's model can be used when at least two diagonal components of the FOT are identical; that is, the fibre orientation distribution is symmetric around one axis, which is normally the flow direction in injection moulding. This could occur in situations where the orientation of fibres slightly differs from that of the perfectly aligned fibres, as this is the case for region B. As a result, it can be concluded that Giordano's approach can be used in the ETC prediction of injection-moulded SFRPs when we are dealing with slightly-perturbed aligned orientations or a weakly expressed orientational preference [12]. These terms are used to describe orientation distributions that are close to perfectly aligned or completely uniform fibre distributions, respectively.
Overall, the predicted and measured ETC along each direction is proportional to the average values of the component of the FOT related to that direction. For instance, the average a xx is the highest in region B, followed by A and C. The same trend was also seen for modelled and measured k xx . Moreover, although there is generally good agreement between the predicted values of ETC along each direction and that of experimental measurements, it seems that the following trends existed for three zones: k xx Experiment < k xx modelling k yy Experiment > k yy modelling k zz Experiment > k zz modelling These trends could be explained by the fact that the MTPS method used for experimental measurement is not a purely one-dimensional thermal conductivity characterisation technique, even though it is an appropriate transient method for thermal conductivity measurement of an anisotropic material. In other words, when the ETC is measured along a particular axis normal to the plane which is in contact with the MTPS sensor, other in-plane components may play a role as well. For example, when measuring along the x-axis (the flow direction), along which the ETC is maximum in all three regions, interfering with the other two lower components of the ETC may decrease the measured k xx compared to actual values, so k xx Experiment < k xx modelling . Conversely, measured k yy and k zz could be overestimated because of the influence of k xx in measurements along y-and z-axes.

Conclusions
Two different homogenisation approaches based on the dielectric theory for pseudooriented inclusions (Giordano's model) and the mean-field homogenisation scheme (the two-step homogenisation method) were used to predict the ETC of an injection-moulded SFRP. When unidirectional fibre alignment or uniformly distributed fibres are assumed, both models predict nearly identical ETCs which are compatible with FEM homogenisation results. For predicting the ETCs in the case of fibres oriented according to a specific FOT, three regions with different fibre orientations were considered in an injection-moulded tensile bar. The injection-packing process was firstly simulated via Autodesk Moldflow ® in order to obtain the FOT of elements associated with each region. By taking micro-CT images and processing the scans in the VoxTex software, it was discovered that the Moldflow FOT predictions and the experimentally obtained FOTs are very compatible.
A MATLAB script was written to find elements whose FOTs satisfy the transverse isotropy criteria, as required for Giordano's model. Regarding the two-step homogenisation, an interpolative model based on Mori-Tanaka was used in the first step to predicting the ETC of unidirectional fibre alignment, followed by the orientation averaging in the second step. It was found that Giordano's model can be used only in the gauge length of the sample. In this region, the fibre orientation distribution differs slightly from the unidirectional alignment. However, in other regions where the fibre orientation distribution noticeably deviates from either unidirectional or random orientations, the transverse isotropy conditions cannot be satisfied. Consequently, only the two-step homogenisation predictions are valid in these regions. Overall, experimental measurements were found to be in good agreement with the predicted values. However, while the k xx measurements were marginally underestimated, the measured k yy and k zz seemed to be slightly higher than those of predictions. This observation can be explained by the fact that thermal conductivity measurements via the MTPS sensor are not exclusively unidirectional.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/polym14163360/s1. Figure S1: The geometry of the injectionmoulded part with the sprue and runners. Table S1: Injection moulding processing parameters.