Homogenization of elastic properties of short-fiber reinforced composites based on measured microstructure data

Abstract

Mechanical properties of short-fiber reinforced composites are crucially influenced by their microstructure. The microstructure itself is mainly governed by the manufacturing process like injection or compression molding. The main contribution of this paper lies in the homogenization of linear elastic properties using experimental microstructural information. For this purpose, the microstructure of injection-molded specimens made of polypropylene reinforced with 30wt.% of short glass fibers are analyzed through micro-computer tomography (μCT) measurements. Applying a recently developed segmentation algorithm, the spatial position, the orientation distribution and the length of the fibers are determined. This data is evaluated in terms of orientation tensors and length distribution, and is used within three mean field approaches: a self-consistent homogenization method, the interaction direct derivative estimate, which is based on the three-phase model, and a two-step bounding method. All methods account for the orientation, the length and the diameter distribution. The numerical results are compared to experimental tensile tests.

Keywords

Micro-mechanics self-consistent method interaction direct derivative fiber-reinforced composite micro-computer tomography

Introduction

Due to the growing effort in building lightweight parts, materials such as short-fiber reinforced composites (SFRCs) have taken over a substantial market share. This class of material is relatively cheap, comparatively easy to process, and provides advantageous specific mechanical properties. However, to design parts made of SFRCs is still a challenging task, since their mechanical behavior strongly depends on the microstructure. The microstructure of SFRCs is characterized by the spatial distribution, the orientation distribution, the distribution of the aspect-ratio, and the volume fraction of the fibers.¹ Particularly in shell-like parts, the fibers are oriented in layers.^2,3 In the boundary layer, the fibers are mostly oriented in the filling direction, and in the core layer in the transversal direction, which is perpendicular to the filling direction, see Figure 1. Thus, a reliable prediction of the mechanical properties depends on the particularities of the microstructure. In order to obtain this microstructural information, X-ray observations by means of computed tomography are, meanwhile, a common approach. Micro-computer tomography (μCT) in laboratories or synchrotron tomography have been used not only to observe and evaluate microstructures but also to characterize damage through ex situ and in situ observations.^4,5 A detailed review on X-ray tomography can be found in Maire and Withers.⁶ This methodology produces a three-dimensional voxel-based picture of the specimen in grayscale, and the question arises of how this information could be used to predict the mechanical properties. As a consequence of the complex microstructure and material behavior, phenomenological approaches are not quite suitable for this task. Alternatively, full field approaches based on solutions of finite element and finite cell models or fast Fourier transformation methods may be applied partially relying on mean microstructure data like orientation tensors, which have been previously extracted from CT data.^7–9 Compared to the details of the microstructure, efforts to calculate full field solutions directly based on CT scans are still handicapped by rough meshes since the computational time is extremely high for these methods.^10,11 Due to their efficient and straight forward algorithms, mean field approaches can be applied to calculate the effective material behavior regarding experimentally determined microstructure data. Especially in the case of SFRCs, the microstructure information, such as the fiber orientation distribution, can be extracted from the CT images in two different ways: In the first approach, the mean fiber orientation distribution in terms of fabric or orientation tensors is calculated directly using voxel data.² The second approach consists of a segmentation of the fibers, and a calculation of the mean orientation distribution based on these segmented data.

Figure 1.

Three-layer model of injection-molded shell-like parts.

Generally, the set of mean field methods can be divided into bounding and estimating methods. The former specify an admissible range of possible effective properties for given microstructural information. Since the first order or simple bounds, known as the Voigt and Reuss bounds,^12,13 only take into account the volume fractions, they enclose a wide range of admissible effective properties. Provided by Hashin and Shtrikman,^14–16 the second-order bounds are based on a variation principle. This approach was extended and applied to anisotropic materials¹⁷ and also nonlinear material behavior¹⁸ as well. Contrary to the bounding methods, estimating approaches give approximations of the effective properties. The Mori-Tanaka,¹⁹ the self-consistent²⁰ (SC), the generalized self-consistent²¹ and the interaction direct derivative (IDD)²² estimates are prominent representatives of such mean field schemes, which take the interaction of the inclusions into account. Based on the Eshelby solution,²³ the general idea of SC is to embed each inclusion in an infinite matrix with the properties of the effective material. Due to this assumption, the SC approach owns an inherently implicit character. The SC method was applied to granular and also to particulate materials with multiple phases and nonlinear properties.^24–26 The major shortcoming of SC is that the interaction between the inclusions and the surrounding matrix is not considered directly. This drawback is the main motivation to apply methods like IDD, which are based on the three-phase model. Here, one inclusion interacts with the matrix directly and with the other inclusions through the effective medium. Additionally, two-step (TS) methods may be applied to homogenize particulate materials with arbitrarily oriented inclusions.²⁷ Within these methods, the microstructure is decomposed into as many domains as there exist different inclusions. Each domain is homogenized individually, and then all domains are homogenized to determine the effective properties. This method can be applied to a much larger class of microstructures compared to second-order bounding schemes²⁸ and also in conjunction with strength prediction.²⁹

In the present work, a composite material consisting of polypropylene reinforced with 30 wt.% of short glass fibers (PPGF30) is examined. The microstructure of this composite material is analyzed and segmented using μCT and a new segmentation algorithm. The segmented microstructure data consist of non-aligned fibers with varying aspect-ratios. It is used within three homogenization approaches directly, namely the SC method, the IDD and a bounding TS approach. In all applied methods, each fiber is considered in the homogenization procedure without a need of orientation averaging of transversal isotropic stiffnesses or compliances for unidirectional structures like it is often done in conjunction with the Mori-Tanaka approach.³⁰

Notation

A direct tensor notation is preferred throughout the text. If tensor components are used, then Latin indices are used and Einstein’s summation convention is applied. Vectors and second-order tensors are denoted by lowercase and uppercase bold letters, e.g., $a$ and $A$ , respectively. The composition of two second-order or two fourth-order tensors is formulated by $AB$ and $AB$ . A linear mapping of second-order tensors by a fourth-order tensor is written as $A = ℂ [B]$ . The scalar product is denoted by $A \cdot B$ . We define the composition operator □ via $(A$ □ $B) [C] = ACB$ , the dyadic product operator ⊗ as $(A \otimes B) [C] = (B \cdot C) A$ and the contraction operator $[[\cdot]]$ with $(a \otimes b) \cdot (ℂ [[a \otimes b]]) = (a \otimes a) \cdot (ℂ [b \otimes b])$ . Arbitrary vectors $a$ and $b$ , second-order tensors $A, B$ and $C$ , and the fourth-order tensor $ℂ$ are used in the foregoing definitions. The identity on symmetric second-order tensors is denoted by $I B s$ . The brackets $〈 \cdot 〉$ indicate ensemble averaging which for ergodic media can be identified with volume averages in the infinite volume limit. A superimposed bar, e.g., $\bar{ɛ}$ , indicates that the quantity refers to the macroscale. The quantity 30 wt.% denotes a mass fraction of 30%.

Experimental methods

Fabrication of specimens

The materials for the test specimens were compounded at the Fraunhofer-LBF to ensure control of all constituents. The polypropylene (PP) Moplen HP500N (Lyondell-Basell) was used as matrix material. In the reinforced configuration (PPGF30), a fraction of 30 wt.% of short glass fibers was added to the polymeric bulk material. To assure stability, durability, and an adequate performance of the reinforced system, three additives were appended to the compound. A fraction of 0.5 wt.% of Irganox B225 (BASF) served as long-term thermal stabilizer. Moreover, 0.07 wt.% of the calcium stearate acid scavenger Ceasit AV (Baerlocher) and 3.0 wt.% of the coupling agent Exxelor PO 1015 (Exxon Mobil) were added. To ensure the determination of the effect of the reinforcing glass fibers on the mechanical properties, these additives were added not only to the reinforced material but also to the pure PP material.

Plates of the dimensions $80 mm \times 80 mm \times 2.5 mm$ were produced by means of injection molding from each of the compounded materials. In order to achieve a homogeneous filling of the plate cavity with a parallel melt front, a triangular gating system was used as infeed.^31,32 The specimens required for mechanical testing were machined from these plates. The specimen geometry has been chosen in accordance with Becker.³¹ For the determination of anisotropic material behavior, the mechanical tests were performed at different orientations. Considering the SFRC, the specimens were extracted at an angle of $0 \circ$ in reference to the filling direction during injection molding, and angles of $45 \circ$ and $90 \circ$ , respectively. For the examination of the non-reinforced material, the specimens were extracted at $0 \circ$ and $90 \circ$ . Details concerning the specimen geometry of the tension rod as well as the different extraction angles are shown in Figure 2(a).

Figure 2.

(a) Locations of specimen extraction from plate (circles distinguish position of samples for microstructural analysis) and (b) example of 2D surface strain field calculated by the correlation software, and averaging area on a tension specimen.

Quasi-static testing methods

The quasi-static tensile tests were conducted on a Z020 testing machine by Zwick-Roell with a velocity of 1 mm/min. Temperature and humidity during the testing procedure were kept at constant values of 23℃ and 50%, respectively. To avoid influences of long-term post-crystallization processes taking place in the PP matrix phase, all samples were tested within a two week time frame.

The deformation of the specimen was observed using 2D digital image correlation.^33,34 Therefore, all specimens were spray-painted with a stochastic grayscale pattern. During the testing procedure, different stages of the sample deformation were captured with a monochromatic CCD-camera equipped with a resolution of $1280 px \times 1024 px$ , a color depth of $8 bit$ and a pixel-size of 8 µm/px. Subsequently, local strain measurement was performed using the software package Vic2D (Limes GmbH). This software calculates the strain field of the grayscale pattern for an area of interest, by analyzing the deformation of predefined grayscale facets between consecutive deformation steps and a reference image. As the correlation results refer to pixel differences, they have to be calibrated to the actual geometry of the specimen. In order to get stress–strain curves, the local strain values in a rectangular area of $12 mm \times 2 mm$ over the region of maximum deformation were averaged and combined with the data from the machine integrated force sensor. The surface strain field and the averaging area are shown in Figure 2(b).

For each test configuration, the final result was generated from average values of five single tests. In Figure 3(a), these results are shown for the pure PP matrix material for specimen oriented in flow ( $0 \circ$ ) and transversal ( $90 \circ$ ) direction. Since the difference between the two curves is in the range of scatter of the single experiments, the matrix material was considered to be isotropic and the results of both orientations were averaged.

Figure 3.

(a) Stress–strain curve of PP matrix material showing details on derivation of Young’s modulus and (b) stress–strain curve of SFRC with 30 wt.% of glass fibers.

The stress–strain response of the SFRC shows a significant deviation depending on the orientation of sample extraction from the sample plate, see Figure 3(b). This mechanical anisotropy is a result of the material orientation induced by the fiber distribution inherent in the injection-molded plates.

The linear elastic properties for the subsequent homogenization procedure have been calculated according to DIN EN ISO 527 1-4.³⁵ Since no distinct definition of a linear elastic region is possible, the data were evaluated in the small strain intervals 0.05–0.25%, 0.05–0.5% and 0.05–1.0%, as shown in Table 1. Any influence of viscous material behavior was neglected. Based on Hooke’s law, Young’s modulus was calculated as the secant gradient in the examined strain range as shown in Figure 3(a). The values of Young’s modulus and Poisson’s ratio for both materials are given in Table 1. The Poisson’s ratios were calculated as the secant gradient in the interval of interest of the lateral strain versus longitudinal strain curve. The elastic properties of the glass fibers were taken from literature and are also shown in Table 1.

Table 1.

Experimental values of Young’s modulus E and Poisson’s ratio ν.

	Strain interval
	0.0005 − 0.0025		0.0005 − 0.005		0.0005 − 0.01
Material	$E [MPa]$	$ν [-]$	$E [MPa]$	$ν [-]$	$E [MPa]$	$ν [-]$
PP 0°	1665.5 ± 6.6%	0.364 ± 6.4%	1608.1 ± 3.3%	0.360 ± 4.4%	1450.4 ± 2.1%	0.371 ± 2.0%
PP $90 \deg$	$1748.4 \pm 15.3 %$	$0.346 \pm 3.0 %$	$1663.7 \pm 6.0 %$	$0.357 \pm 3.4 %$	$1501.3 \pm 3.1 %$	$0.368 \pm 2.2 %$
mean PP	1706	0.355	1636	0.359	1476	0.37
PPGF30 $0 \deg$	$4482.4 \pm 2.9 %$	$0.271 \pm 11.2 %$	$4137.4 \pm 4.8 %$	$0.273 \pm 10.1 %$	$3671.4 \pm 2.9 %$	$0.294 \pm 6.3 %$
PPGF30 $45 \deg$	$3540.9 \pm 5.4 %$	$0.304 \pm 4.7 %$	$3250.6 \pm 5.1 %$	$0.300 \pm 7.4 %$	$2947.7 \pm 4.3 %$	$0.310 \pm 8.0 %$
PPGF30 $90 \deg$	$3452.4 \pm 2.2 %$	$0.217 \pm 2.1 %$	$3227.7 \pm 1.8 %$	$0.222 \pm 3.5 %$	$2938.8 \pm 1.2 %$	$0.238 \pm 5.0 %$
Glass fibers	73000	0.22	–	–	–	–

Analysis of microstructure

Micro-computer tomography measurement

The cylindrical specimens for micro-computer tomography (μCT) measurements were milled from the aforementioned plate at the positions defined in Figure 2(a). To achieve the maximum resolution needed to accurately resolve the reinforcement fibers, the diameter of the specimen should have a maximum diameter of $4 mm$ . An exemplary μCT specimen is shown in Figure 4(a) on the planar mount.

Figure 4.

(a) Specimen for analysis of microstructure on μCT mount and (b) μCT measurement setup.

The samples have been analyzed with a SkyScan 1072-100 CT apparatus, using its maximum resolution of 1.8 µm. As shown in Figure 4(b), the specimen is exposed to a low intensity X-ray while being rotated around its vertical axis. The absorption of the X-ray depends on the density of the constituents in the sample and, therefore, the amount of the X-ray that reaches the detector varies with the sample orientation. For each full rotation, the sample is analyzed in small steps on the vertical axis, resulting in layer-wise voxel information of the density distribution. The detected images are transformed to slice images as shown in Figure 5 and to a three-dimensional geometrical representation of the sample.

Figure 5.

μCT slice image of PPGF30.

Evaluation of microstructure

With the new evaluation method, it is possible to determine the fiber length, the diameter, the orientation and the local position distributions by means of an iterative analysis of μCT images.³⁶ This analysis method is based on model assumptions considering geometrical and symmetrical properties of the fiber objects, such as a constant diameter and a small curvature without turning points. From these assumptions, heuristic approaches for the identification attempts of single fibers are derived. Furthermore, the model assumptions allow for assigning voxels of the μCT images to individual fibers in spite of the inherent image noise. As shown in Figure 6, the new approach consists of four functional units: in the image initialization phase, the μCT images are prepared for the evaluation process by means of filtering and binarization (see Figure 5). To achieve a reasonable binarization, the proper fiber volume fraction has to be defined. For the present work, the fiber volume fraction, known from the plate fabrication process, was confirmed by a thermogravimetric analysis. In the next step, the iterative Monte-Carlo pattern recognition process is started. Each of the parallel executed recognition processes evaluates integrals over a spherical region at random positions to identify centers of separate fiber sections.³⁶ After the identification of these fiber centers, first approximations of the fiber orientation $n$ are determined by using the main eigenvalues of second-order moments. Exact orientation and length of the fiber are then obtained by the calculation of integrals over a cylindrical region for different directions $n = \pm n \pm δ n$ . Parallelization allows multiple scans of fibers at variable starting points, whereby fiber recognition can be improved. Each recognition process is followed by the control process that merges information from parallel calculations, finds duplicate fibers, deletes voxels of reasonably well detected fibers from the input image and adjusts the Monte-Carlo pattern parameters before starting the next iteration of the fiber recognition process. In case of reaching a predefined detection accuracy, the control process terminates the recognition procedure and initiates the geometric analysis of the detected reinforcement fibers. The resulting data consisting of the position, orientation, length and diameter of each recognized fiber is called segmented μCT data in the following and can be used to reconstruct the microstructure, see Figure 7.

The algorithm was evaluated in three ways: first, artificially generated fiber data with three different lengths were mapped on a regular mesh, which correspond to a 3D voxel picture. The fiber positions were chosen randomly, however, the length and orientation were correlated. An analysis of this volume-image reveals an overestimation of the fiber length of up to 8% for the short fibers. Secondly, the polymer material of the composite was removed via combustion and the residual fibers were measured via microscope. Qualitatively, similar length distributions were found. Thirdly, a visual comparison of the volume-image and the detected fibers was undertaken. This confirmed visually the mentioned small overestimation of the fiber length. Additionally, the reproducibility of the algorithm was checked by repeated measurement of three different specimens. A standard deviation of the length of approximately 5% was found. Further details on the image analysis method and its properties can be found in Glöckner et al.³⁶

Figure 6.

Fiber recognition process.

Figure 7.

Segmented μCT data for specimen at the positions (a) $(- 10)$ , (b) (00), and (c) (10).

Properties of microstructure

In the reconstructions of the microstructure of the evaluated specimen shown in Figure 7, the characteristically graded mesostructure of injection-molded plates made of SFRCs can be recognized. This phenomenon is well known in the literature.³² Since the segmented μCT data allow a detailed analysis of the properties of the microstructure, the μCT datasets were partitioned into 20 equal-sized layers in through-thickness direction of the plate. The mean aspect-ratio of the fibers $\bar{a} γ$ and the mean orientation distribution in terms of the second-order moment tensor $N γ$ of the fiber orientations have been determined for each layer. The second-order moment tensor, also called fabric tensor³⁷ or orientation tensor,³⁸ is generally defined by

N = \int S f (n) n \otimes n d S .

(1)

Herein,

f (n)

is the fiber orientation distribution function,

n

the fiber axis orientation and

d S

a surface element of the unit sphere

S : = {n \in ℝ 3 : | | n | | = 1}

. Utilizing Dirac’s delta distribution

δ (n - n β)

on unit vectors, the empirical orientation distribution function³⁷ for N orientations is given by

f (n) = \frac{1}{N} \sum_{β = 1}^{N} δ (n - n β) .

(2)

Inserting the last equation in the definition of the second-order tensor given in equation (1) results in a discrete form of the second-order fabric tensor

N = \frac{1}{N} \sum_{β = 1}^{N} n β \otimes n β .

(3)

Since the available microstructure data does not contain only the orientations of the fiber axes, but also their length and diameter, it is reasonable to consider these quantities during the calculation of the fabric tensor. Thus, in this work, fabric tensors are considered in a weighted manner. For each of the 20 layers, the weighted fabric tensor

N γ

is calculated as

N γ = \sum_{β = 1}^{N γ} ω β n β \otimes n β .

(4)

In the last equation, $N γ$ is the number of fibers in layer γ, and $ω β$ is the volume fraction of the fiber with the axis direction $n β$ . In Figure 8, the mean aspect ratio of the fibers $\bar{a} γ$ , the diagonal components (N₁₁, N₂₂, and N₃₃), and the eigenvalues (λ₁, λ₂, and λ₃) of the second-order fabric tensor $N γ$ are shown for all three segmented μCT datasets. The indices $(\cdot) 11, (\cdot) 22$ , and $(\cdot) 33$ correspond to the $e 1, e 2$ , and $e 3$ -directions, respectively, introduced in Figure 2(a), whereby, the $e 2$ -direction is the filling direction. From Figure 8(a), (d), and (g), it is obvious that the mean aspect ratio of the fibers in each layer is not constant throughout the thickness of the plate. Especially in the case of the μCT-specimens (−10) and (00), see Figure 2(a), the mean aspect ratio of the fibers in the layers near the bounding sections is slightly larger than those in the core section. The overall mean aspect ratio and the overall mean length for each specimen are approximately 26.3 and 330 μm, respectively. The overall mean aspect ratio is marked in Figure 8(a), (d), and (g).

Figure 8.

(a), (d), and (g) Mean aspect ratio of the fibers in each layer $\bar{a} γ$ and the whole dataset $\bar{a}$ ; (b), (e), and (h) Diagonal components of $N γ$ for each layer; (c), (f), and (i) eigenvalues of $N γ$ for each layer.

The fiber orientation distribution is also not constant throughout the thickness. In Figure 8(b), (e), and (h), the diagonal components N₁₁, N₂₂, and N₃₃, and in Figure 8(c), (f), and (i), the eigenvalues of the fabric tensor $N γ$ , are shown. In particular, the former three diagrams highlight the sectional structure of injection-molded plates made of SFRCs: in the bounding sections, layers 1–2 and 19–20 in Figure 8, respectively, the majority of the fibers are oriented in the filling direction. In the layers 11 and 12 at the core section of the specimen, the fibers are oriented in the transversal direction.

In the literature,^39,40 the fiber orientation distribution is usually characterized only through the diagonal components N₁₁, N₂₂, and N₃₃ of the second-order fabric tensor $N γ$ . The comparison of the diagonal components with the eigenvalues of $N γ$ of the 8th and 14th layer for the (00) data in Figure 8(e) and (f) highlights that the diagonal components do not specify completely the orientation distribution. According to the diagonal components of $N γ$ in the aforementioned layers, the orientation distribution seems to be planar isotropic. The eigenvalues for the same cases indicate, however, a planar isotropic distribution for layer 14, not for layer 8. These results are confirmed through the stereographic projections of all fiber axes orientations in each of the layers 8, 11, and 14 in Figure 9. Here, the out-of-plane direction coincides with the filling direction ( $e 2$ -direction) and the horizontal direction with the transversal direction ( $e 1$ -direction). Thus, the vertical direction in these figures corresponds to the through thickness direction, which is perpendicular to the injection-molded plate. In Figure 9(a), a preferred orientation of the fibers between the out-of-plane direction and the transversal direction in layer 8 can be seen. In the pole figure for layer 11, Figure 9(b), the preferred orientation of the fibers corresponds to the transverse direction. In the case of layer 14 in Figure 9(c), an approximately planar isotropic fiber orientation distribution can be noticed. The through-thickness orientation distribution in the other segmented μCT datasets is qualitatively and quantitatively similar to the section discussed before.

Figure 9.

Pole figures for the layers 8, 11, and 14 of μCT data at position (00).

The stereographic projections of the entire segmented μCT data are shown in Figure 10. Since each μCT dataset contains more than 6000 fibers, continuous representations of the stereographic projections are shown. In all examples, the fiber orientation distribution is only approximately planar isotropic. There exists rather a preferred orientation different from the filling direction. This can be affiliated with the manufacturing process and the flow conditions during manufacturing. The mean length $\bar{l}$ , the diagonal components, and the eigenvalues of $N$ for the entire datasets are summarized in Table 2.

Figure 10.

Stereographic projection of fiber axes for the μCT data $(- 10)$ , (00), and (10).

Table 2.

Number of fibers, mean length with standard deviation, diagonal components and eigenvalues of $N$ for three μCT datasets.

Data	N	$\bar{l}$ [μm]	$\bar{N} 11; \bar{N} 22; \bar{N} 33$	$eig (\bar{N})$
$(- 10)$	6329	330.28 ± 165.9	$(0.381; 0.593; 0.026)$	$(0.594; 0.380; 0.026)$
(00)	6355	330.18 ± 167.4	$(0.392; 0.584; 0.024)$	$(0.635; 0.341; 0.024)$
(10)	6200	332.66 ± 164.3	$(0.404; 0.571; 0.025)$	$(0.613; 0.361; 0.025)$

Homogenization of linear elastic properties

Modeling preliminaries

In the present work, the polymeric matrix and the arbitrarily oriented glass fibers are assumed to be linear elastic and isotropic with piecewise constant properties. The microstructure of this composite consists of a matrix reinforced with N fibers. The matrix phase is specified by the stiffness tensor $ℂ B M$ and the volume fraction $c B M$ . Each fiber is characterized with the stiffness $ℂ α$ , the volume fraction $c α$ , the length $l α$ , the diameter $d α$ , and the orientation of its axis $n α$ . The axis is parameterized in spherical coordinates with the polar angle θ and azimuthal angle ϕ, see Figure 11. The fibers are approximated by spheroids of equal length and volume. Thus, the half-axes of the spheroid are correlated with the length and diameter of the corresponding fiber according to the following prescription: $2 a 1 α = l α$ and $2 a 2 α = \sqrt{3 / 2} d α$ , whereby $a 1 α$ and $a 2 α$ are the first and the second half-axis of the spheroid related to fiber α. The total fiber volume fraction is $c B F = 1 - c B M = \sum_{α = 1}^{N} c α$ . The material behavior of the fibers is assumed to be uniform, i.e., $ℂ α = ℂ B F \forall α$ .

Figure 11.

Parametrization of the fiber axis $n$ in spherical coordinates.

Self-consistent homogenization

The effective stiffness tensor $\bar{ℂ}$ of a $N + 1$ -phase particulate composite is generally given by

\bar{ℂ} = ℂ B M + \sum_{α = 1}^{N} c α (ℂ α - ℂ B M) A α,

(5)

where

c α

denotes the volume fraction of fiber α and

A α

the fourth-order strain localization tensor.⁴¹ N corresponds to the number of fibers inside the considered volume element, see Figure 7. The index

(\cdot) N + 1

is replaced with the index

(\cdot) B M

, which denotes the matrix phase. The special case of an spheroidal approximation of the fibers and the assumption that each fiber is embedded in an infinite homogeneous matrix with the properties equal to the effective material

\bar{ℂ} = ℂ B SC

, lead to the following expression for the strain localization tensor

A α

A α = A (ℂ B SC, ℂ α, Z α) = (I B s + ℙ_{α}^{B SC} (ℂ α - ℂ B SC)) - 1 .

(6)

In this representation, the strain localization tensor $A α$ depends on the effective stiffness $ℂ B SC$ , the fiber stiffness $ℂ α$ and the spheroidal approximation of the fiber $Z α$ , where the orientation of the fiber axis is included

| | x | | Z 2 = x \cdot (Z_{α}^{B T} Z α x) \leq 1 .

(7)

In the last equation, $x$ specifies a position vector in three-dimensional space. Since, an spheroidal geometry of the fibers has been assumed, the polarization tensor $ℙ_{α}^{B SC}$ can be determined explicitly⁴¹

ℙ_{α}^{B SC} (ℂ B SC, Z α) = \frac{1}{4 π det (Z α)} \int S ℍ (n) (n \cdot (Z_{α}^{- B T} Z_{α}^{- 1} n)) - 3 / 2 d S

(8)

with

ℍ (n) =

I B s (K - 1

□

(n

⊗

n)) I B s

, and

K = ℂ B SC [[n \otimes n]] . ℙ_{α}^{B SC}

depends on the stiffness

ℂ B SC

and the spheroidal geometry of the fibers

Z α

. In the last equation,

d S

is a surface element of the unit sphere

S : = {n \in ℝ 3 : | | n | | = 1}

and

det (Z α)

represents the determinant of

Z α

. The integral in equation (8) has been calculated numerically using Gauss quadrature after a decomposition of the surface of the sphere according to Weber.⁴²

Combining equations (5) and (6), and using the self-consistent approximation $\bar{ℂ} = ℂ B SC$ gives an implicit equation for the stiffness $ℂ B SC$

ℂ B SC = ℂ B M + \sum_{α = 1}^{N} c α (ℂ α - ℂ B M) (I B s + ℙ_{α}^{B SC} (ℂ α - ℂ B SC)) - 1 .

(9)

This implicit equation for the unknown tensorial quantity $ℂ B SC$ was solved numerically using a Newton-Raphson algorithm combined with a line search procedure. Despite the isotropic modeling of the matrix and fiber material, the SC homogenization delivers a non-isotropic effective stiffness tensor, which reflects the anisotropic fiber orientation distribution.

Interaction direct derivative estimate

The interaction direct derivative (IDD) estimate, proposed by Zheng and Du²² is based on the generalized self-consistent scheme²¹ (GSCS), which for its part is based on the three-phase model. In the three-phase model, the inclusions are embedded in a finite matrix material region, directly. This inclusion-matrix cell itself is embedded in the unbounded initially unknown effective medium. The difference between the three-phase model and IDD can be identified in the estimation of the stresses in the inclusions: contrary to the three-phase model, Zheng and Du²² assume an unbounded medium with the properties of the matrix for the embedding of the inclusion in the matrix material. Du and Zheng⁴³ have proven that the stresses in the inclusions are well approximated by this assumption with an error of second order of the inclusion volume fraction. The advantage of the IDD is its explicit structure, which is valid for multi-phase composites with different material symmetries and distributions.

The following prescription gives the IDD estimation of the properties of the effective medium

ℂ B IDD = ℂ B M + (I B s - \sum_{β = 1}^{N} c β (ℂ β - ℂ B M) ℕ β ℙ_{β}^{B D}) - 1 \times \sum_{α = 1}^{N} c α (ℂ α - ℂ B M) ℕ α,

(10)

with

ℕ α = (I B s + ℙ α (ℂ α - ℂ B M)) - 1

. Here,

ℙ α = ℙ (ℂ B M, Z α)

is Hill’s polarization tensor as defined in equation (8) but

ℂ B SC

replaced with

ℂ B M

. If the matrix-inclusion cell takes on an ellipsoidal shape, then,

ℙ_{α}^{B D} = ℙ (ℂ B M, Z_{α}^{B D})

is the polarization tensor corresponding to an ellipsoidal inclusion with geometry of the matrix-inclusion cell

Z_{α}^{B D}

, which is embedded in an infinite matrix with the stiffness

ℂ B M

. The shape of the matrix-inclusion cell defines the inclusion distribution in the composite.²² In the present work, this shape is assumed to be equal to the shape of the spheroidal approximation

Z α

of the corresponding inclusion itself.

Particularly, if the distribution is constant for all inclusions ( $Z_{α}^{B D} = Z B D$ ), the IDD estimate is equivalent to the Hashin–Shtrikman estimate of Ponte Castañeda and Willis,⁴⁴ which is based on the Hashin–Shtrikman variational structure in the form developed by Willis.^45,46 A detailed discussion of the relation of IDD to the estimate of Ponte Castañeda and Willis⁴⁴ and, e.g., Mori-Tanaka can be found in Zheng and Du²² and Du and Zheng.⁴³

A two-step bounding method

The SC and IDD approaches deliver estimates of the effective elastic properties. The two-step bounding method (TS) is a simple method providing an admissible range of possible effective properties. The first step consists of the decomposition of the microstructure into as many domains as there are different fibers. Fibers of equal shape and equal orientation are combined into unidirectional domains. In each domain, the fiber volume fraction corresponds to the total fiber volume fraction $c B F$ . For each combination of fiber and matrix, the effective elastic properties are calculated using the unidirectional (UD) special case of the second-order Hashin–Shtrikman (HS) bounds⁴¹

ℂ_{α}^{B UD} = c B F ℂ α A α + (1 - c B F) ℂ B M A B M .

(11)

Here, the strain localization tensors

A α

and

A B M

do not depend on the effective material, but on a reference material

ℂ 0

A β = M β á M ñ - 1, M β = (I B s + ℙ_{0}^{B UD} (ℂ β - ℂ 0)) - 1

(12)

with

á M ñ = c B F M α + (1 - c B F) M B M

, where β denotes the fibers α or the matrix

B M

, respectively. In the case of unidirectional aligned ellipsoidal inclusions, the polarization tensor

ℙ_{0}^{B UD}

is known explicitly.⁴¹ Specifying the softer matrix material or the stiffer fiber material as the reference material

ℂ 0

gives a lower

ℂ_{α}^{B UD -}

and an upper bound

ℂ_{α}^{B UD +}

, respectively, for the quasi-coated fiber domains. The resulting bounds exhibit a transversely isotropic behavior.

Within the second step, again HS bounds are calculated assuming an isotropic two-point correlation function for the domains.⁴⁷ Herein, corresponding bounds are combined: the lower (upper) HS bound of the domains is homogenized with the lower (upper) HS bound for the granular system. Thus, this procedure results in the stiffnesses tensors $ℂ B TS -$ and $ℂ B TS +$ , respectively

ℂ B TS \pm = \sum_{α = 1}^{N} \frac{c α}{c B F} ℂ_{α}^{B UD \pm} A_{α}^{\pm} = \sum_{α = 1}^{N} \frac{c α}{c B F} ℂ_{α}^{B UD \pm} M_{α}^{\pm} á M \pm ñ - 1

(13)

with

M_{α}^{\pm} = (I B s + ℙ 0 (ℂ_{α}^{B UD \pm} - ℂ_{0}^{\pm})) - 1, á M \pm ñ = \sum_{β = 1}^{N} \frac{c β}{c B F} (I B s + ℙ 0 (ℂ_{β}^{B UD \pm} - ℂ_{0}^{\pm})) - 1 .

(14)

In the last equation, the definition of

ℙ 0

depends on modeling assumptions concerning the shape and distribution of the domains. Assuming an isotropic two-point correlation of the domains,

ℙ 0

corresponds to the spherical polarization tensor. Another possible assumption would be an ellipsoidal two-point correlation. Due to simplicity, the isotropic two-point correlation has been assumed. In the case of the upper (lower) HS bound, for

ℂ_{0}^{\pm}

the maximum (minimum) isotropic part of all stiffness tensors of the domains

ℂ_{α}^{B UD +}

(

ℂ_{α}^{B UD -}

) is taken.

Results and discussion

The mean values of Young’s moduli and Poisson’s ratios of the isotropic polypropylene material in the strain region 0.0005–0.0025, as shown in Table 1, were used as input parameters for the homogenization of the composite. The fiber mass fraction is 30 wt.%, which corresponds to a fiber volume fraction of $c B F = 0.13$ .

Combining these data with the segmented μCT data, the effective elastic properties were calculated according to the methods described in the foregoing section. To compare the experimentally measured Young’s moduli to the computational results, for each homogenized stiffness, the direction-dependent Young’s modulus⁴⁸ $E (d)$ , given in equation (15), was determined and evaluated in the appropriate direction

\frac{1}{E (d)} = d \otimes d \cdot S [d \otimes d] .

(15)

In combination with the direction-dependent bulk modulus $K (d)$

\frac{1}{3 K (d)} = I \cdot S [d \otimes d]

(16)

the linear elastic behavior is completely defined. In equation (15), the direction

d

is parametrized with spherical coordinates, and

S = ℂ - 1

is the compliance tensor.

In the following first, the numerical results based on the three segmented μCT datasets are compared to each other, and second, the numerical estimations are compared to the experimentally measured Young’s moduli. In Figure 12, the contours of the direction-dependent Young’s moduli in the x–y plane are shown in dependence of the azimuthal angle as depicted in Figure 11. It can be seen that the SC and IDD estimates are located between the bounding TS results for all three datasets. SC predicts a stiffer behavior than IDD. The direction-dependence of Young’s modulus is more pronounced by SC and IDD than by the TS methods. The TS bounding methods exhibit a rather large range of admissible Young’s moduli. Furthermore, in Figure 12, it can be observed that the direction-dependence of Young’s moduli is very similar for the (00) and the (10) data, see Figure 12(b) and Figure 12(c). The data at the position $(- 10)$ yield a different direction-dependence of the Young’s modulus, see Figure 12(a). This can be attributed to the flow conditions during the manufacturing process. The same phenomenon is visible in the stereographic projections of all fibers in Figure 10. The maximum Young’s modulus is located apart from the filling direction ( $e 2$ -direction), which is also an effect of the inhomogeneous flow process.

Figure 12.

Comparison of numerical Young’s moduli for the three μCT datasets.

Since the volume element analyzed with the μCT measurements is small compared to the size of the tensile specimen, the three μCT datasets are combined to an overall dataset for the comparison of the numerical and experimental results. In Figure 13, the shapes of the direction-dependent Young’s moduli are shown on the x₁–x₃, x₂–x₃ and x₁–x₂ plane for the overall data. As the specimens for the tensile test were prepared in the x₁–x₂ plane, the experimental results are added to Figure 13(c). It can be observed that the direction-dependent Young’s modulus of the overall data exhibits a defined symmetry in the x₁–x₃ and the x₂–x₃ plane. In the x₁–x₂ plane, no characteristic symmetry is obvious. In Table 3, the experimental and the numerical results for all homogenization methods for the overall data are given.

Figure 13.

Orientation dependence of Young’s modulus in the (a) x₁–x₃ plane, (b) x₂–x₃ plane and (c) x₁–x₂ plane for the overall data.

Table 3.

Comparison of experimental and numerical Young’s modulus $E [MPa]$ for the overall dataset.

	$E B \exp$	$E B TS -$	$E B IDD$	$E B SC$	$E B TS +$
$0 \circ$	4.482	3.798	4.809	5.463	7.754
$90 \circ$	3.452	3.035	3.741	4.167	6.842
$45 \circ$	3.541	2.986	3.661	4.087	6.912

In addition to the qualitative evaluation of the numerical results compared to experimental measurements in Figure 13 and Table 3, it is interesting to notice that the predicted anisotropy of the direction-dependent Young’s modulus essentially depends on the homogenization procedure. Especially, the upper bound of the two-step approach does not represent the anisotropy accurately. It predicts a rather isotropic material behavior.

To observe this more quantitatively, the ratios $E 90 / E 0, E 45 / E 0$ and $E 90 / E 45$ between the Young’s modulus values in different directions were calculated for experimental and numerical results separately. In Figure 14, the comparison of these ratios is given for the overall μCT data.

Figure 14.

Anisotropy ratios of Young’s modulus for the overall μCT data.

The TS methods predict a smaller anisotropy for $E 90 / E 0$ and $E 90 / E 45$ and $B TS -$ also for $E 45 / E 0$ . SC predicts a greater anisotropy for $E 90 / E 0$ and $E 45 / E 0$ and IDD only for $E 45 / E 0$ . The lower bound of TS reflects the experimental anisotropy ratios better than the upper bound. Generally, all ratios of the upper bound are on a similar elevated level. This conveys a nearly isotropic material behavior.

Summary and conclusions

The present paper discusses the examination of the thermoplastic composite PPGF30 through tensile tests in three different directions, and its microstructure through μCT measurements in combination with a stochastic fiber analysis approach. The microstructure was analyzed by means of the resulting segmented μCT data. Three different mean field homogenization methods were applied to estimate the effective elastic properties directly based on segmented microstructural data.

Relying on the experimental and numerical results, the following conclusions can be made:

Anisotropy: The evaluated injection-molded SFRC specimens show an anisotropic material behavior. The experiments in different directions indicate the importance of attaining detailed knowledge of the microstructure in order to dimension parts made of such a material.

Inhomogeneity: Important properties of the microstructure, like the fiber orientation distribution and fiber length distribution, can be analyzed by using the segmented μCT data. Thus, a layer-wise analysis reveals the dependence of these properties on the position in the specimen. Thereby, a more accurate observation of the microstructure of SFRCs is possible.

Microstructural information: All presented mean field homogenization approaches consider the segmented μCT data and take advantage of the information about the distribution of the fiber axes, fiber lengths and fiber radii.

Two-step homogenization: In order to get an admissible range for the effective properties, two-step bounding methods have been applied: In the first step, the Hashin–Shtrikman bounds for unidirectional aligned fiber domains and, in the second step, the Hashin–Shtrikman bounds for granular systems have been calculated.

SC and IDD estimates: Both methods, SC and IDD, allow for a direct consideration of the segmented μCT data. IDD uses the matrix material as reference material. SC is an implicit approach using the effective properties as reference material. In addition to the shape, the orientation, and the material properties, IDD accounts for a statistical description of the distribution of the inclusions. Here, a spheroidal distribution equal to the shape of the spheroidal approximation of the fibers has been used.

Experimental vs numerical results: The self-consistent homogenization method delivers a stiffer material behavior compared to experimental measurements. Due to the inherent characteristics of the SC approach,^49,50 this method tends to overestimate the elastic properties for the case of hard inclusions embedded in a soft matrix similar to the problem at hand. The Young’s modulus values predicted by IDD are located between the lower TS bounds and SC. Like SC, this approach delivers also a stiffer material behavior compared to experiments. The IDD predictions of the Young’s modulus values show the closest agreement with experimental measurements. In terms of the anisotropy ratios, SC and IDD show similar trends. Both methods predict a close agreement with experiments for $E 90 / E 0$ , a larger anisotropy for $E 45 / E 0$ and a smaller for $E 90 / E 45$ . The TS methods, except $B TS -$ for $E 45 / E 0$ , show larger anisotropy ratios compared to experiments. Two reasons are responsible for that: first, in the second step, the reference material has to be isotropic. Therefore, the minimal and maximal isotropic part of the domain stiffnesses, respectively, has been taken as reference material. Second, an isotropic two-point correlation function has been assumed by applying the spherical polarization tensor for $B TS -$ and $B TS +$ .

The present investigation shows the modeling of the elastic behavior of SFRCs, using homogenization techniques in combination with detailed microstructural data from a micro-computer tomography image analysis method. The homogenization results demonstrate the need of realistic microstructural data for the simulation of parts made of SFRCs. In engineering practice, the fiber configuration of complex SFRC components is derived from injection molding simulations. Further investigations should, thus, be dedicated to the validation of these simulations, taking advantage of the new possibility of extracting detailed μCT fiber information from real parts. Additionally, further experimental measurements like bending tests should reveal properties of the composite material, which are due to the inhomogeneous fiber orientation distribution in through-thickness direction.

Footnotes

Funding and Acknowledgment

The research documented in this manuscript was partially initiated by the German-Canadian research group “Integrated engineering of continuous-discontinuous long fiber reinforced polymer structures”. The support of the Institute of Engineering Mechanics at KIT by the German Research Foundation (DFG) is gratefully acknowledged.

Conflict of interest

None declared.

References

Thomason

. Micromechanical parameters from macromechanical measurements on glass reinforced polypropylene. Compos Sci Technol 2002; 62: 1455–1468.

Bernasconi

Cosmi

Dreossi

. Local anisotropy analysis of injection moulded fibre reinforced polymer composites. Composites Science and Technology 2008; 68: 2574–2581.

Bernasconi

Cosmi

Hine

. Analysis of fibre orientation distribution in short fibre reinforced polymers: A comparison between optical and tomographic methods. Compos Sci Technol 2012; 72: 2002–2008.

Bull

Helfen

Sinclair

. A comparison of multi-scale 3D X-ray tomographic inspection techniques for assessing carbon fibre composite impact damage. Compos Sci Technol 2013; 75: 55–61.

Bull

Spearing

Sinclair

. Observations of damage development from compression-after-impact experiments using ex situ micro-focus computed tomography. Compos Sci Technol 2014; 97: 106–114.

Maire

Withers

. Quantitative X-ray tomography. Int Mater Rev 2013; 59: 1–43.

Demirci

Acar

Pourdeyhimi

. Finite element modelling of thermally bonded bicomponent fibre nonwovens: tensile behaviour. Computat Mater Sci 2011; 50: 1286–1291.

Duschlbauer

Böhm

Pettermann

. Computational simulation of composites reinforced by planar random fibers: homogenization and localization by unit cell and mean field approaches. J Compos Mater 2006; 40: 2217–2234.

Müller V, Böhlke T, Kabel M, et al. Homogenization of linear elastic properties of short-fiber reinforced composites – A comparison of mean field and voxel-based methods (submitted. 2014).

10.

Düster

Sehlhorst

Rank

. Numerical homogenization of heterogeneous and cellular materials utilizing the finite cell method. Computat Mech 2012; 50: 413–431.

11.

Węglewski

Bochenek

Basista

. Comparative assessment of Young’s modulus measurements of metal–ceramic composites using mechanical and non-destructive tests and micro-CT based computational modeling. Computat Mater Sci 2013; 77: 19–30.

12.

Voigt

. Ueber die Beziehung zwischen den beiden Elasticitätsconstanten isotroper Körper. Annalen der Physik 1889; 274: 573–587.

13.

Reuss

. Berechnung der Fließgrenze von Mischkristallen auf Grund der Plastizitätsbedingung für Einkristalle. ZAMM J Appl Math Mech/Zeitschrift für Angewandte Mathematik und Mechanik 1929; 9: 49–58.

14.

Hashin

Shtrikman

. On some variational principles in anisotropic and nonhomogeneous elasticity. J Mech Phys Solid 1962; 10: 335–342.

15.

Hashin

Shtrikman

. A variational approach to the theory of the elastic behaviour of polycrystals. J Mech Phys Solid 1962; 10: 343–352.

16.

Hashin

Shtrikman

. A variational approach to the theory of the elastic behaviour of multiphase materials. J Mech Phys Solid 1963; 11: 127–140.

17.

Böhlke

Lobos

. Representation of Hashin–Shtrikman bounds of cubic crystal aggregates in terms of texture coefficients with application in materials design. Acta Materialia 2014; 67: 324–334.

18.

Talbot

DRS

Willis

. Some simple explicit bounds for the overall behaviour of nonlinear composites. Int J Solid Struct 1992; 29: 1981–1987.

19.

Mori

Tanaka

. Average stress in matrix and average elastic energy of materials with misfitting inclusions. Acta Metallurg 1973; 21: 571–574.

20.

Kröner

. Bounds for effective elastic moduli of disorderd materials. J Mech Phys Solid 1977; 25: 137–155.

21.

Christensen

. Solutions for effective shear properties in three phase sphere and cylinder models. J Mech Phys Solid 1979; 27: 315–330.

22.

Zheng

. An explicit and universally applicable estimate for the effective properties of multiphase composites which accounts for inclusion distribution. J Mech Phys Solid 2001; 49: 2765–2788.

23.

Eshelby JD. The determination of the elastic field of an ellipsoidal inclusion, and related problems. In: Proc Royal Soc London Ser A, Math Phys Sci 1957; 241: 376–396.

24.

Budiansky

. Thermal and thermoelastic properties of isotropic composites. J Compos Mater 1970; 4: 286–295.

25.

Gillman

Matouš

Atkinson

. Microstructure-statistics-property relations of anisotropic polydisperse particulate composites using tomography. Phys Rev E 2013; 87: 022208–022208.

26.

Miled

Sab

Le Roy

. Effective elastic properties of porous materials: Homogenization schemes vs experimental data. Mech Res Commun 2011; 38: 131–135.

27.

Pierard

Friebel

Doghri

. Mean-field homogenization of multi-phase thermo-elastic composites: a general framework and its validation. Compos Sci Technol 2004; 64: 1587–1603.

28.

Doghri

Friebel

. Effective elasto-plastic properties of inclusion-reinforced composites. Study of shape, orientation and cyclic response. Mech Mater 2005; 37: 45–68.

29.

Kaiser

Stommel

. Micromechanical modeling and strength prediction of short fiber reinforced polymers. J Polym Eng 2012; 32: 43–52.

30.

Ghossein

Lévesque

. A comprehensive validation of analytical homogenization models: the case of ellipsoidal particles reinforced composites. Mech Mater 2014; 75: 135–150.

31.

Becker F. Entwicklung einer Beschreibungsmethodik für das mechanische Verhalten unverstärkter Thermoplaste bei hohen Deformationsgeschwindigkeiten. PhD Thesis, Martin-Luther-Universität Halle-Wittenberg, Zentrum für Ingenieurwissenschaften, 2009.

32.

Mlekusch

Lehner

Geymayer

. Fibre orientation in short-fibre-reinforced thermoplastics I. Contrast enhancement for image analysis. Compos Sci Technol 1999; 59: 543–545.

33.

Tarigopula

Hopperstad

Langseth

. A study of large plastic deformations in dual phase steel using digital image correlation and FE analysis. Exp Mech 2008; 48: 181–196.

34.

Palanivelu

De Pauw

Van Paepegem

. Validation of digital image correlation technique for impact loading applications 2009; vol. 1, Brussels, Belgium: EDP Sciences, pp. 373–379.

35.

DIN EN ISO 527 1-4. Plastics – Determination of tensile properties, DIN German Institute for Standardization, 2012.

36.

Glöckner R, Heiliger C and Kolling S. A Monte-Carlo algorithm for 3D fibre detection from microcomputer-tomography. (submitted).

37.

Kanatani

. Distribution of directional data and fabric tensors. Int J Eng Sci 1984; 22: 149–164.

38.

Advani

Tucker

III CL

. The use of tensors to describe and predict fiber orientation in short fiber composites. J Rheol 1987; 31: 751–784.

39.

Laspalas

Crespo

Jimenez

. Application of micromechanical models for elasticity and failure to short fibre reinforced composites. Numerical implementation and experimental validation. Comput Struct 2008; 86: 977–987.

40.

Vincent

Giroud

Clarke

. Description and modeling of fiber orientation in injection molding of fiber reinforced thermoplastics. Polymer 2005; 46: 6719–6725.

41.

Willis JR. Variational estimates for the overall response of an inhomogeneous nonlinear dielectric. In: Ericksen JL, Kinderlehrer D, Kohn R, Lions JL (eds) Homogenization and effective moduli of materials and media. No. 1 in The IMA Volumes in Mathematics and its Applications. New York: Springer, 1986, pp.247–263.

42.

Weber

Kenmeugne

Clement

. Improvements of multiaxial fatigue criteria computation for a strong reduction of calculation duration. Computat Mater Sci 1999; 15: 381–399.

43.

Zheng

. A further exploration of the interaction direct derivative (IDD) estimate for the effective properties of multiphase composites taking into account inclusion distribution. Acta Mech 2002; 157: 61–80.

44.

Ponte Castañeda

Willis

. The effect of spatial distribution on the effective behavior of composite materials and cracked media. J Mech Phys Solid 1995; 43: 1919–1951.

45.

Willis

. Bounds and self-consistent estimates for the overall properties of anisotropic composites. J Mech Phys Solid 1977; 25: 185–202.

46.

Willis

. Variational principles and bounds for the overall properties of composites. Continuum Model Discrete Syst 1978; 2: 185–215.

47.

Ponte Castañeda

Suquet

. Nonlinear composites. Adv Appl Mech 1997; 34: 171–302.

48.

Böhlke

Brüggemann

. Graphical representation of the generalized Hooke’s law. Technische Mechanik 2001; 21: 145–158.

49.

Hill

. A self-consistent mechanics of composite materials. J Mech Phys Solid 1965; 13: 213–222.

50.

Torquato

. Random heterogeneous materials: microstructure and macroscopic properties, New York: Springer Science & Business Media, 2002.