A computational modeling approach based on random fields for short fiber-reinforced composites with experimental verification by nanoindentation and tensile tests

In this study a modeling approach for short fiber-reinforced composites is presented which allows one to consider information from the microstructure of the compound while modeling on the component level. The proposed technique is based on the determination of correlation functions by the moving window method. Using these correlation functions random fields are generated by the Karhunen–Loève expansion. Linear elastic numerical simulations are conducted on the mesoscale and component level based on the probabilistic characteristics of the microstructure derived from a two-dimensional micrograph. The experimental validation by nanoindentation on the mesoscale shows good conformity with the numerical simulations. For the numerical modeling on the component level the comparison of experimentally obtained Young’s modulus by tensile tests with numerical simulations indicate that the presented approach requires three-dimensional information of the probabilistic characteristics of the microstructure. Using this information not only the overall material properties are approximated sufficiently, but also the local distribution of the material properties shows the same trend as the results of conducted tensile tests.


Introduction
The compatibility of thermoplastic material for automated serial production like mold injection allows high production rates with a reasonable price per piece. Adding short fibers or nanoparticles to the base material of pure plastics leads to a significant increase of the stiffness and strength of the material without losing the ability to process the material by automated serial production. Hence, this combination results into a high interest in short fiber-reinforced composite (SFRC) not only in the automotive industry.
The main disadvantage due to the reinforcing elements is the spatial fluctuation of the material properties. The representation of the non-equally distributed material properties is hence, challenging and connected with high computational costs due to the use of probabilistic methods. This consequently is worthwhile to represent the probabilistic B Natalie Rauter natalie.rauter@hsu-hh.de 1 Helmut-Schmidt-University/ University of the Federal Armed Forces, Holstenhofweg 85, 22043 Hamburg, Germany characteristics of the microstructure on the component level by appropriate stochastic methods without extensive computational costs.
The spatial distribution of the material properties is mainly assigned to the finite fiber length as well as a production process induced variation of the fiber orientation due to different melt flow velocities and shear forces over the cross-section [35]. Since the fiber length is usually much smaller in comparison to the component the microstructural information must be transferred from the micro-and mesoscale to the component level. Therefore, a multi-scale approach appears advantageous for a sufficient modeling procedure for SFRC.
In the literature many different approaches are presented for the modeling of the probabilistic characteristics of reinforced materials. One approach is based on the orientation tensor introduced in [1]. In [3] the orientation tensor is implemented directly in the material description. An extended approach is used in [10], where the orientation tensor is combined with a material model based on a coupled micromechanical and phenomenological approach. The suitability is also discussed in [50]. A second approach is based on the extended finite element method (XFEM), which was developed for the numerical simulation of crack growth and interfaces [49]. In [32] the method is applied to reinforced concrete. It is also used in [58], where it is combined with a cohesive zone model. All these approaches have in common the lack of a scale transition. A first attempt for a scale transition is presented in [61] where XFEM is combined with a Monte Carlo simulation to estimate the size of a representative volume element (RVE) for random composites.
One technique for the representation of spatial varying parameters with a scale transition is the use of random fields [71] and hence, allow the modeling of inhomogeneous material properties [52]. For the discretization of random fields for the use in a numerical simulation the Karhunen-Loève expansion is a suitable tool. This allows a separation of the deterministic numerical modeling and the probabilistic representation of the material properties. In context of structural mechanics this approach is also referred to as the stochastic finite element method (SFEM) [47,55]. In [39] SFEM is combined with cohesive zone models for the analysis of crack growth in fiber reinforced cementitious composites. The fluctuation of the volume fraction is represented by random fields in [24]. The influence of the fiber mass fraction, the window size, and the contrast between the phases on the correlation structure is discussed in [62]. Finally, a combination of the SFEM with the concept of RVE is presented in [23,65].
The main objective of this research is to establish a robust computational modeling approach for components made of SFRC using homogeneous second-order random fields discretized by the Karhunen-Loève expansion. In this study a tensile test specimen is used exemplary. First the microstructural properties of the material are extracted from a two-dimensional micrograph as shown in [63]. Using the probabilistic characteristics of the fiber length, fiber orientation, fiber diameter, and the information about the layered structure of the tensile test specimen, cross-sections with different micromechanical representations of the specimen are analyzed and compared with each other. This is done by linear elastic numerical simulations that are validated by experimental investigations. For the numerical simulation artificial microstructures are generated on the basis of the derived microstructural characteristics. Using these artificial microstructures, a numerical model is established and the Young's modulus is determined by simulating tensile tests. In addition, the Young's modulus is also experimentally obtained by nanoindentation. This allows one to reduce the three-dimensional effects on the results and hence, to determine the near surface material properties [20]. Based on the results of the numerical and experimental investigations on the mesoscale the best approach is derived by taking also into account the aspect of computational costs. The selected approach is then used to expand the procedure to a representation of the material on the component level by taking into account the spatial distribution of the material proper-ties induced by the microstructural characteristics. This is achieved by representing the components of the elasticity tensor by homogeneous second-order random fields. Therefore, the material properties are still modeled probabilistically and not deterministically. This also indicates that there is no scale transition to the macroscale. For the discretization of these fields the Karhunen-Loève expansion is used, which requires the determination of correlation functions first.
Finally, the influence of the probabilistic microstructure characteristics on material properties obtained on the component level is analyzed. This is done by simulating tensile tests on a numerical basis. Again the numerically obtained values for the Young's modulus are compared with experimental results. Consequently, this study gives a comprehensive analysis of a computational modeling approach for components made of SFRC by numerical simulation, that is validated by experimental investigations.
The structure of the presented work is as follows. Sect. 2 gives a brief overview of the most important theoretical background including homogeneous second-order random fields, the Karhunen-Loève expansion for the discretization of random fields as well as the main aspects of the multi-scale modeling with respect to Hill's condition. In Sect. 3 the different specimens used in this study are presented. This is followed by the extraction of the probabilistic characteristics from a micrograph in Sect. 4. The numerical simulations and experimental investigations on the mesoscale and the impact of the resulting microstructural characteristics on the component level are shown in Sects. 5 and 6, respectively. Finally, Sect. 7 gives a summary and a conclusion of the presented work.

Second-order random fields
Among others probabilistic quantities are described by random variables Z (ω). A realized value, e.g. of a material parameter, is denoted by z. In context of random fields Z (ω, x) the random variables are also assigned to spatial coordinates x. For the synthesis of random fields by using the Karhunen-Loève expansion, it is necessary that the variance of the random field as well as the random variable are finite. In this case the following definition for random fields holds and one speaks of second-order random fields [4]. Their main properties are briefly presented below. Random fields are characterized by moments of their probability distribution [66,71]. In general, the n-th moment of a single random variable Z is defined as where f Z (z) is the probability density function. Based on this definition the first moment of a random field, also called expected value, is given by The mean square of Z (x) is given by the second moment. Besides these two moments central moments with respect to the expected value are introduced. The deviation of a value to the expected value is given by the second central moment, also known as the variance. Reformulating Eq.
(2) provides the definition for the variance of a random field Z (x) In addition, the standard deviation is often used, which is derived from the variance by Hence, the expected value as well as the variance are functions of the spatial coordinates x. However, in case of a homogeneous random field both the expected value and the variance become constants [71]. The observation of a random field at different locations x i is described by corresponding random variables Z = Z i . Besides the relation between an observation and the expected value of the random field the relation between observations at different points of a random field is of great interest. This relation is expressed by the covariance and reads for two random variables Z 1 and Z 2 Usually this expression is reduced to a dimensionless correlation parameter ρ 12 by dividing Eq. (6) by σ 1 and σ 2 Equations (6) and (7) are referred to as the auto-covariance and auto-correlation, respectively, if Z 1 and Z 2 are part of the same random field Z (x). Otherwise, Eqs. (6) and (7) are the cross-covariance and cross-correlation, respectively. Since the probability density function is unknown for most times the random field is represented by a discrete number of realizations ω i [56]. In this case the mean of the discrete values is used as expected value of the random field. In addition, the variance is rewritten as Finally, the dimensionless correlation coefficient for two random variables Z 1 and Z 2 is given by

Auto-correlated random fields
The Karhunen-Loève expansion is a generalization of the Fourier transform to probabilistic processes. A commonly used formulation for the expansion of a stochastic process X (ω, x) is given by [8] X where Z n are a set of uncorrelated random variables and λ n and φ n (x) are an eigenpair of a Fredholm integral equation of the second kind. The derivation of the integral equation is briefly summarized below by using the analogy to the Fourier transform as done in [45]. Hence, stochastic processes X (x) are represented by a series expansion of orthonormal functions φ n (x) as With this series expansion the position-dependent part is separated from the probabilistic component of the stochastic process, where the functions φ n (x) are orthonormal. Furthermore, it is beneficial to split the random field additively into a deterministic and a stochastic field. This separation can be written as where the mean of the random field μ(x) represents the deterministic part and α(ω, x) gives a random field with an expected value of zero. Due to this separation the expected value of the series expansion E[X (x)] and the expected value of the random variables Z n equal zero. In addition to the expected value the variance of the random variable is given by In context of random variables orthogonality is equivalent to independent and uncorrelated random variables, respectively. In this case the covariance reads for values of the same homogeneous random field X (x) at the two locations x and x Introducing μ = 0 and the independence of the random variables leads to In addition, combining the series expansion in Eq. (12) with the orthogonality provides for the random variables, where φ n (x) are the eigenfunctions of the auto-correlation. The corresponding eigenvalues are given by λ n . Therefore, before generating random fields by using the Karhunen-Loève expansion the auto-correlation must be determined. Expressing X (x) and X (x ) by using Eq. (12) and taking into account the independence of the random variables the auto-correlation is written as which leads finally to the following eigenvalue problem For a known correlation function the eigenvalues λ n and eigenfunctions φ n (x) can be calculated by solving the Fredholm integral equation of the second kind in Eq. (19). A closed solution is only possible for a few correlation functions. One example here is the exponential correlation function, which is given for homogeneous random fields by [17,53,66] Cov with the correlation length b i . An important property is the multiplicative decomposition of the different directions, which is used for the generation of two-dimensional random fields, as shown in Sect. 2.3. Introducing the correlation function the eigenvalue problem in Eq. (19) is rewritten for a one-dimensional field [−l, l] as The solution of this eigenvalue problem is given by [66] with for odd i and with for even i. The values of ω i are obtained by solving and for odd and even values of i, respectively. The computational solution of Eqs. (28) and (29) is not trivial since the range boundaries are poles with a transition from a positive to a negative sign. Usually, algorithms for the determination of zeroing are based on a sign transition. However, precision limitations may lead to an inaccurate evaluation of the boundary values. Therefore, these algorithms require a manually defined offset to ensure a correct calculation of the zeroing, see Fig. 1. In [8] a different approximation for ω i is given with This approach does not require a case differentiation between odd and even values of i. However, the problem regarding the sign transition at poles is still not solved, resulting in a missing sign transition. Therefore, in this work a further approach for a robust determination of zeroing is achieved by multiplying Eqs. (28) and (29). In this case there is no sign transition at the poles but within the interval for the determination of ω. Hence, a manually defined offset is no longer necessary. Figure 1 shows a comparison of the different approaches, for a correlation length of b = 5 mm and a field length of l = 200 mm.

Two-dimensional random fields
Applying random fields to problems in the field of mechanical engineering mostly requires the representation of twodimensional or even three-dimensional components. The Karhunen-Loève expansion is adapted to two-dimensional components by separating the correlation structure of the different directions following the multiplicative decomposition of the correlation functions. Hence, the two-dimensional eigenfunctions and eigenvalues are expressed by the product of the eigenfunctions and eigenvalues of each direction and For a rectangular component of the dimension l 1 × l 2 and an exponential correlation function this leads to [18] with the individual solution of each direction and respectively. The eigenvalues and eigenfunctions of these two individual equations are identical to the those presented in Sect. 2.2.

Apparent material properties
One technique to extract the relation between information of different locations x i is the moving window method [2,19]. Within this method a window of predefined size is used to extract segments of a larger microstructure. For these extractions the material properties are determined by numerical simulation. With the material properties available at different locations the correlation functions are obtained by a curve fit. The numerical simulations for the determination of the material properties are conducted in accordance with Hill's condition [33] σ : = σ : .
This procedure ensures the correct calculation of the effective material properties. For a RVE [34] the homogeneous material properties are written as and = S eff : σ (38) respectively. Here, C eff and S eff are the effective stiffness and compliance tensor. The strains and stresses σ hold information about the microstructure. In contrast to this · gives the volume average, defined by For the determination of the effective material properties among others the following two kinds of boundary conditions can be derived from Hill's condition. Based on the average strain theorem with linear elastic material properties as well as the average stress theorem [76] these boundary conditions read and respectively. Here Eq. (40) holds a boundary condition of Dirichlet type, because only pure displacement boundary conditions with the constant macroscopic strain 0 are defined on the complete surface. In contrast to this Eq. (41) holds a boundary condition of Neumann type, since only pure traction boundary conditions, with a constant macroscopic stress t 0 are defined on the complete surface. In case of a RVE the material properties obtained for these two boundary conditions are equal and hence, the material is assumed to be homogeneous. If the used volume is smaller than the RVE one speaks of apparent overall properties [30,31,37].
To determine these properties, the same procedure is used as for the effective material properties. However, in this case the results depend on the size of the extracted volume as well as the boundary conditions. Such kind of extraction on the mesoscale is referred to as a statistically volume element (SVE) [54]. Figure 2 gives an overview of all specimens and their orientations used in this research. The raw material of the conducted analyses is a polybutylenterephthalat (PBT) filled with glass fibers. The fiber mass fraction is specified with 30% (PBT GF 30). The initial component is a plate with the dimensions of 300 mm×300 mm×3 mm fabricated by an injection mold process (1). With this plate a coordinate system is introduced that is used throughout the whole investigation. The main objective of this research is the implementation of a linear-elastic numerical modeling procedure for the numerical simulation of components made of SFRC that provides not only the global material behavior but also information about the spatial distribution of the material properties on the mesoscale. This is done exemplary for a tensile test specimen. Therefore, in a first step tensile test specimens are cut out of the plate, while the orientation of the tensile test specimen coincides with the main flow direction of the melt front (2). The tensile test specimens are used to obtain the global Young's modulus on the component level in Sect. 6.2 of the material for a validation of the presented procedure. As the numerical model should also conclude information about the material microstructure a cross-section is defined (3), for which the microstructural properties are determined based on a micrograph (4). The selected cross-section of 3.37 mm × 3 mm is parallel to the main flow direction and hence, coincide with the load direction of the tensile test in Sect. 6.2 (2). By doing so, the cross-section is used to analyze the Young's modulus of the specimen's microstruc-ture on a numerical basis in Sect. 5.1 (5). Afterwards, the results are validated experimentally by nanoindentation in Sect. 5.2. Here, the specimens are small cuboids, where one side of the cuboid is perpendicular to the micrograph crosssection (6). The surface of the plate is the second side. This allows one to conduct indentation tests in the main flow direction as well as in thickness direction of the plate and specimens, respectively. The cuboids have an edge length of 20 mm × 5 mm × 3 mm. The numerical simulation on the component level assuming linear elastic material behavior refers to the same crosssection as the microstructural analysis. Here, the dimensions coincide with the standard measurement length of a tensile test specimen defined in [38]. Therefore, the edge length is 50 mm × 3 mm.

Specimen specification
The material properties of the two components used for the numerical simulations are provided in Table 1.

Theoretical description of the fiber orientation
The orientation of a single fiber is expressed by the unity vector p, that points in the direction of the fiber axis as shown in Fig. 3. The components of this vector are expressed by the angles introduced in Fig. 3 as [43,63] To derive the orientation from the micrograph the components of the orientation vector are formulated in terms of the geometrical properties of the fiber cross-section. In this case the minor and major axis of the ellipse as well as its horizontal and vertical dimensions with respect to the coordinate system are used, which leads to the following expressions for the trigonometrical functions of ϕ and θ Here Δ h and Δ v are the projection of the elliptic cross-section on the x− and y−axis, respectively. The parameters l min and l ma are introduced in Fig. 3. For a group of fibers, the orientation is given by the second order orientation tensor A introduced in [1] as Here Ψ (p) is the probability density distribution function for the fiber orientation defined in [13] and S the surface of the unit square. In context of a component made by mold injection and with respect to the introduced coordinate system the elements on the diagonal represent a fiber orientation in the melt flow direction (A 11 ), perpendicular to the melt front (A 33 ), and in thickness direction (A 22 ). The sum of these Angle and cross-section definition for the description of the fiber orientation and fiber length, see also [26,43,63] three elements equals one [73]. With the fiber orientation vector defined in Eq. (42) the elements of the orientation tensor is written as In terms of the geometrical properties of the elliptic crosssection this leads to [63]

Methodology
For the numerical simulation of SFRC the probabilistic characteristics of the microstructure are required. In [63] a procedure is presented, that allows one to derive the secondorder orientation tensor of the fibers from a micrograph. This procedure is adapted to obtain the probabilistic characteristics and their spatial distribution of the fiber length, diameter, and orientation as well as the second-order orientation tensor from a micrograph in this study. All presented results are The panorama image is generated by stitching 323 individual images with a magnification of 1500. This leads to a resolution of 558 px per 100 µm. Hence, the procedure implies ergodicity for the extracted microstructural properties with respect to the manufactured plate with an edge length of 300 µm.
The main steps for the extraction of the probabilistic characteristics are depicted in Fig. 4. Based on the original micrograph first the fibers are detected before the image is transformed to a binary image. This binary image is then further processed with the image processing toolbox of Matlab ® . Using the function regionprops the orientation, centroid as well as major and minor axes of each ellipse are determined. With these values the elements of the orientation tensor A i j given in Eq. (47) are calculated.

General results
For the determination of the probabilistic characteristics only fibers are taken into account that are within the boundaries of the micrograph. Fibers reaching over the micrograph edges are neglected since a correct determination of the major axis is not possible. In total 9721 fibers are detected within the micrograph cross-section. Below the results of the main fiber characteristics like fiber length (major axis), fiber diameter (minor axis), fiber orientation, and fiber volume fraction as well as their spatial distribution are discussed. Figure 5a shows the spatial distribution of the obtained fiber length. As the probability of detection for a single fiber depends strongly on the orientation of the fiber, in [63] a procedure is proposed that weights the fiber length with respect to its probability of detection. It is assumed, that fibers parallel to the micrograph cross-section are less likely to be detected in a micrograph than fibers perpendicular to it. In context of the presented framework this component of the fiber orientation is given by the angle θ as defined in Fig. 3 and Eq. (44). Using this angle the weighting factor is given by [63] The results given in Fig. 5a already include the weighting. Even then a profuse amount of short fibers with a length of less than 20 µm is indicated. This is also supported by a mean fiber length of 46.6 µm. Besides this, the results lead to the conclusion of a layered structure due to the mold injection process as shown in [5,12,14,16,48,60]. This is discussed more in detail in Sect. 4.3.2.
Next is the fiber orientation ϕ, that describes the angle between the major axis of the ellipse and the horizontal axis as shown in Fig. 3. The results are depicted in Fig. 5b. Again the distribution is plotted against the specimen thickness coordinate. The distribution clearly shows a preferred orientation of the fibers along the specimen and therefore, with the melt flow direction. In contrast to the fiber length, the orientation is independent of the location of the fiber. The preferred orientation of the fibers and therefore, an orientation dependence of the material properties of the compound with respect to the overall specimen orientation is further confirmed by the correlation analysis of the fiber length and the fiber orientation. With an increasing fiber length, the fiber orientation shows a decreasing spreading around an orientation angle of 0 • , see Fig. 5c.
The overall results of the fiber diameter are given in Fig. 6. The probability density function of the fiber diameter meets the expectation of a normal distribution with a mean of 9.6 µm. Compared to the results of other studies this value appears a little bit low. Usually, the diameter of glass fibers is approximated with 10 µm [25]. The fiber diameter corresponds with the minor axis of the pattern recognition results. If the fibers are cut perpendicular to the fiber length these two magnitudes coincide very well. Only if the fiber is almost parallel to the micrograph cross-section there might be a significant deviation. However, since the probability of detection of these fibers is small and there is no weighting considered for the fiber diameter the influence is neglectable. Another influence is more significant. As shown in Fig. 4 (left) there is a dark shadow between the fiber and the matrix material that cannot be clearly assigned to one of the two components. The difference between the Finally, the overall fiber volume fraction is 16.7%, which equals a fiber mass fraction of 28% based on the material properties provided in Table 1. For the raw material a fiber mass fraction of 30% was specified. Therefore, the overall fiber mass fraction appears a little bit below the specification. The fiber mass and fiber volume fraction are connected by the density of the matrix and fiber material, respectively. The precise density of the two components is unknown and the data provided in Table 1 are taken from literature, which usually gives only approximated values since the specific values vary between different charges. Taking these variations into account a fiber volume fraction of 16.7% leads to a fiber mass fraction in a range of 27% up to 30%, which meets the specifications.

Analysis of the layered structure
During the injection mold process, the fiber orientation varies over the specimen thickness due to different flow velocities and shear forces [35]. With respect to the fiber orientation the cross-section of the component made by mold injection is divided into different layers. Figure 7 gives an overview of the two commonly used orientation schemes [5,12,14,16,48,60].
The results of the fiber length distribution over the thickness indicates that the used specimen shows this kind of layered structure as well. To analyze the layered structure more in detail the components A 11 and A 33 of the orientation tensor introduced in Eq. (47) are calculated and plotted over the cross-section thickness. As discussed in Sect. 4.1 the probability of detection of a fiber orientation parallel to the molt flow direction is given by Perpendicular to the melt flow direction the probability is expressed by A 33 (ϕ = θ = 0 • ). The results depicted in Fig. 8 confirm the assumption based on the fiber length distribution. Five individual layers are recognizable.
Furthermore, the values of A 11 and A 33 indicate also that the orientation component in thickness direction A 22 is small compared to the other two components. This supports the overall assumption of mainly horizontally aligned fibers, since high values of A 22 would indicate vertically orientated fibers.
In the following sections a numerical modeling procedure for SFRC is presented. Since the main objective is the representation of the microstructural effects on the component level the layered structure must also be taken into account. Therefore, the layer thicknesses are determined based on the obtained orientation tensor as well as the distribution of the fiber length. The fiber length is an indirect measurement of the fiber orientation as short fibers with an almost circular cross-section are orientated perpendicular to the micrograph cross-section (θ = 0 • → A 33 ), whereas long fibers are parallel to the micrograph cross-section (θ = 90 • → A 11 ). Based on the results shown in Fig. 8 for the orientation tensor elements A 11 and A 33 as well as the results for the fiber length distribution Fig. 5a the cross-section consists of five layers. The thickness of the skin layers at the top and bottom of the cross-section is approximated with 100 µm, the core layer in the middle of the cross-section is 300 µm thick. Therefore, the shell layers between the skin and core layers have a thickness of 1250 µm each. Core layer Fig. 9 Histogram of the fiber length for the shell and core layers To determine the probabilistic characteristics of the numerical simulation in Sect. 5 the layers are analyzed individually and the results are compared with each other. However, due to the limited number of fibers within the skin layers these are not taken into account. In Fig. 9 the histogram of the fiber length is shown for the two shell layers (left) and the core layers (right). The results of the two shell layers (colored in red and blue) coincide very well. Therefore, these two layers can be modeled using the same probabilistic characteristics for the fiber length.
The fiber diameter distribution is independent of the layer scheme as the fiber diameter is only minimally affected by the mold injection process. The main influence is the initial fiber production [25]. Therefore, no further analysis is conducted regarding the fiber diameter.
For the fiber orientation both aspects the information about the layer as well as the correlation with the fiber length as depicted in Fig. 5b, c are taken into account. As already derived from Fig. 5b the fiber orientation distribution described by ϕ does not vary significantly between the different layers. However, as indicated in Fig. 5c the corre-lation between ϕ and the fiber length must be considered for the modeling of artificial SFRC microstructures.
The presented results seem to contradict the assumption of ergodicity of the micrograph at first, because the probabilistic characteristics of the fiber length depends on the spatial coordinates due to the layered structure of the cross-section. This means that the microstructural properties are no longer homogeneous or stationary, which implies that they cannot be ergodic. But as shown in Fig. 9 and discussed before the fiber length characteristics for the two analyzed shell layers show almost identical characteristics. The same holds for the orientation of the fibers as depicted in Fig. 8. Based on this and the fact that the fiber diameter characteristics are independent of the layered structure, each layer is assumed to be homogeneous and hence, layer-wise ergodic estimators appear appropriate.

Modeling on the mesoscale
In this section a procedure is introduced for the numerical modeling of SFRC on the mesoscale. The first approach is based on the micrograph used in Sect. 4 to extract the probabilistic characteristics of the microstructure. This is followed by the use of artificial microstructures. Here, the numerical model still represents the microstructure and hence, shows the same probabilistic properties on the microscale as the micrograph.

Model description
The numerical model on the mesoscale represents the micrograph used in Sect. 4 to extract the probabilistic properties of the microstructure. It has a rectangular shape with an edge length of 3370 µm × 3000 µm. The material definition and mesh generation follow the results presented in [59]. Therefore, the model is discretized by a mapped mash consisting of second order Lagrange elements with an edge length of 10 µm × 10 µm and the material properties are passed to the integration points by arrays that represent the underlying microstructure. In contrast to the results presented in [59] a plain strain state is assumed as the micrograph is a crosssection of a larger tensile test specimen. An overview of the specimen's geometry as well as the material properties of the two components are provided in Sect. 3, Fig. 2(3) and Table 1.
For the numerical simulation of the tensile test, the horizontal displacement is fixed at the left edge and the vertical displacement is fixed at the lower edge of the numerical model. Two load cases are analyzed with this model. First a load of 100 MPa is applied at the right edge, second is a where the indices h and v represent the horizontal and vertical component and the indices r and u the right and upper edge, respectively. All numerical simulations are based on linear theory. Finite deformations are not taken into account.

Numerical simulation based on the micrograph
In a first step the material properties of the micrograph crosssection itself are analyzed. This is done by generating three different arrays, one for each material parameter. The arrays represent the micrograph cross-section with a resolution of 1 µm × 1 µm. This is in accordance with the findings in [59]. With these arrays the material properties are passed to the integration points. Table 2 contains the results of the mean strain as well as Young's modulus in horizontal and vertical direction.

Numerical simulation based on the artificial microstructure
Below the analyzed configurations, the generation process of the microstructure as well as the numerical results are presented Analyzed configurations The analysis of the microstructure emanates from a layered structure of the cross-section. As presented in Sect. 4.3.2 the cross-section can be divided into five layers. Based on this finding three different configurations are investigated here. With this approach the influence of the layered structure on the material properties is analyzed. This is important for the representation of the microstructure by random fields in Sect. 6. The following, three different configurations are under investigation. (I) For the first configuration the probabilistic characteristics of the microstructural properties of the whole cross-section are used. This leads to some kind of homogenized material properties. (II) Next is the genera- tion of an artificial microstructure based on the probabilistic characteristics of the microstructural properties of the shell layers only. Therefore, the influence of the skin layers as well as the core layer are neglected. (III) In a last step the core layer is taken into account as well. Therefore, in this case the numerical model consists of three different layers. As discussed in Sect. 4.1 the probabilistic properties of the skin layers are not sufficient to derive the probabilistic characteristics of these layers. Due to this an artificial five-layer microstructure is not investigated. Furthermore, the representation of all layers contradicts a compact and cost saving approach of modeling SFRC.

Generation of the microstructure
The microstructure of the SFRC is generated in accordance with the procedure presented in [59]. In addition, the correlation between the fiber length and fiber orientation is taken into account. Hence, the generation of the microstructure comprises the following steps. First the fiber centroid is generated by a Monte Carlo sampling. If necessary, the coordinates are then assigned to one of the layers and the fiber length is sampled based on the corresponding layer characteristics. The layer identification is only required for configuration III, though. The overall characteristics are used to set a fiber diameter due to the missing correlation as shown in Sect. 4.3. For the assignment of the fiber orientation however, the correlation with the fiber length is essential. Finally, the fiber is added to the preset area if there is no overlap with an already existing fiber. The procedure is repeated until a predefined fiber volume fraction is reached. Since the fiber volume fraction of the different layers are not equal, the threshold varies for each configuration and the different layers. Table 3 gives an overview of the fiber mass fraction for each configuration and layer. In Fig. 10 the shell layers of the microstructure and an artificial microstructure are mapped one above the other. There are no major differences between the images recognizable.
The generated microstructures are analyzed in the same way as the numerical model of the micrograph. Instead of the micrograph cross-section the arrays represent now the artificial microstructure. For statistical reasons the simulation is done for 500 different microstructures of each configuration.  Table 4 gives an overview of the numerical results for the artificial microstructures. In addition to the mean value and standard deviation of the Young's modulus for the horizontal and vertical direction the deviation to the results of the micrograph are calculated. The results of the artificial microstructures meet the results based on the micrograph within a maximal deviation of 5% for the horizontal and 8% for the vertical direction. In both directions the highest deviations are observed for configuration I, where the different layers and their probabilistic characteristics are not distinguished. The deviation to the micrograph decreases for configuration II and III. The impact of modeling three different layers instead of only one is marginal compared to the computational costs. Hence, it appears to be convenient to reduce the modeling of the microstructure to the shell layers for the determination of the correlation structure, that is required for the numerical simulation on the component level in Sect. 6. This assumption is verified experimentally first. Furthermore, the numerical results indicate only a minor anisotropy of the material. This can be assigned to the derived aspect ratio of 4.9 based on the mean values of the fiber length and fiber diameter.

Experimental validation
The results of the numerical simulation presented in Sect. 5.1 are verified by an experimental investigation of the Young's modulus based on the material, that was characterized in Sect. 4. As the presented results and the probabilistic characteristics of the material are based on a two-dimensional analysis of the cross-section the Young's modulus is experimentally determined by nanoindentation. Compared to standard tensile tests this procedure allows one to measure the near surface material properties [20] and hence, in combination with a projected indentation area of approximately 7000 µm 2 characterize the material on the mesoscale.

Framework nanoindentation
For the experimental determination of the material properties by nanoindentation the slope of the force-displacement curve of the unloading process is calculated, see Fig. 11 (left). The indentation or reduced modulus is obtained by [51,57] where S is the stiffness of the contact, A p the projected area of the indentation tip for a contact depth h c , and β an empirical correction factor of the indenter uniaxial symmetry when using pyramidal indenter [29,69]. As done in [41] β is set to one in this analysis. This is possible as the value of β usually varies between 1.02 and 1.08, which leads to a deviation of the indentation modulus of up to 3% [51,57]. The projected area of a perfectly sharp Berkovich tip, which is used in this study, is usually approximated by [11] A p = π h 2 c tan 2 ϕ = 24.5h 2 c , with an effective semi-angle of ϕ = 70.32 • [6]. Due to imperfections, a calibration of the projected surface is mostly used. For this calibration the projected surface is approximated with the function [51] A p (h c ) = 24.56h 2 In Fig. 11 (right) the indentation mark of a Berkovich tip is depicted. Finally, between the indentation modulus and the Young's modulus of the specimen material the following relation holds Here, E i and ν i are the Young's modulus and Poisson's ratio of the indentation material, respectively. In this study the indentation tip is made of diamond with a Young's modulus of 1143 GPa and a Poisson's ratio of 0.0691 [40].

Setup and results
For the experimental verification of the numerical simulation on the mesoscale the indentation modulus is measured at the surfaces of the cuboid. First is the measurement in xdirection. The measurement points are arranged in a grid of 37 × 9 points, that cover an area of 9000 µm × 1600 µm. Consequently, the distance between the measurement points is 200 µm in z-and 250 µm in y-direction, respectively. For this configuration the indentation modulus is measured with a maximum indentation force of 500 mN and 1500 mN. The results are given in Table 5. The indentation modulus for both load cases match within the standard deviation. Hence, the indentation modulus is independent of the maximum indentation force in a range of 500 mN to 1500 mN. As in the numerical simulations, the indentation modulus is also determined experimentally in a second direction. Therefore, the indentation test is repeated in y-direction. Here only one load case is conducted with a maximum load of 1500 mN. The result is also provided in Table 5.
In Table 5 not only the results of the indentation tests are summarized, but also the calculated values of the Young's modulus are given. This is done be evaluating Eq. (53). For this the Poisson's ratio of the analyzed cross-section is obtained by using the analytical material model introduced by Halpin and Tsai [27,28,70] and the probabilistic character-  Fig. 11 Scheme load cycle nanoindentation [44,67] (left), micrograph of the indentation location after unloading (right) However, the obtained Poisson's ratio for the correction in x-direction cannot be used for a correction in y-direction, because the assumption of transversally isotropic material behavior is not applicable for the analyzed cross-section. The main reason is the varying orientation of the fibers. This holds not only when considering the layered structure of the cross-section, but also when reducing the representation of the cross-section to the shell layers only. This conclusion is supported by the mismatch of the Young's modulus in x-and y-direction, as shown in Table 4. Table 6 gives an overview of the experimentally and numerically obtained values of the Young's modulus on the mesoscale. For both directions the results of the numerical simulations meet the results of the experimental investigation within the standard deviation. Consequently, the presented modeling approach on the mesoscale is verified by the conducted experiments. Furthermore, the reduction of the cross-section representation to the shell layers, as suggested in Sect. 5.1.3, appears to be valid.

Modeling on the component level
So far the modeling process is based on the microstructural properties of the analyzed cross-section (see Sect. 4) given by the mean and standard deviation of the fiber length, fiber diameter, and fiber orientation. The aim of this section is to provide a numerical modeling approach on the component level including the microstructural properties of short fiberreinforced material. This is done by representing the material properties by homogeneous random fields. Since random fields also involve information about local correlation between the realization at each point the correlation structure and correlation function are determined first. In this context the artificial microstructures, introduced in Sect. 5 are used to derive the correlation functions of the elasticity tensor coefficients.
The procedure for the generation of the artificial microstructures for SFRC introduced in Sect. 5.1 leads to almost identical mechanical properties as the microstructure of the micrograph of Sect. 4.1. This is not only found on a numerical basis, but also verified experimentally by nanoindentation in Sect. 5.2. Hence, in this section the procedure is used to expand the modeling approach to the component level. A crucial issue for the representation of material properties is their positive nature. Since Gaussian random fields assume normal distributed underlying random variables, negative realizations are possible [17,66]. Therefore, due to the elliptic property of the elasticity partial differentiation operator the use of Gaussian random fields for the representation of material properties is not suitable in context of multiscale modeling of heterogeneous material [64]. Hence, the use of non-Gaussian random fields is required to ensure a stochastic solution of second-order and the positive nature of the elasticity coefficients, which is also state of the art [9,21,22,42,46,68].
A non-Gaussian random field M(ω, x) is given by where G is a non-linear mapping operator and X (ω, x) a centered, homogeneous Gaussian random field as introduced in Sect. 2 [22]. In this study discretized random fields are used. Although the stochastic solution obtained here (using a Gaussian field) is not of second-order, computational (sample-wise) results may not be strongly impacted given the choice of model parameters and are thus presented below.

Correlation structure
As shown in Sect. 2.1 to Sect. 2.3 random fields describe the spatial distribution of random variables. The main characteristic is the correlation function which gives an information about the dependence of values at different locations. Hence, the correlation functions of the elasticity tensor elements must be determined first. The correlation analysis is done in accordance with the procedure presented in [59] and hence, is based on the assumption of a SVE. As homogeneous second-order random fields are used to represent the material properties of a tensile test specimen, the correlation analysis as well as all remaining numerical simulations are based on the plain strain assumption. Using pure displacement and pure traction boundary condition as introduced in [74] the auto-and crosscorrelation lengths of the elasticity tensor coefficients are calculated, respectively. Taking into account that the elasticity tensor coefficients for a plain strain state depend not only on the in-plane material properties but also on ν 23 it is not possible to derive the engineering constants from the elasticity tensor coefficients without any further assumption. Due to this the correlation structure is obtained for the elasticity tensor coefficients instead of the engineering constants. In a next step these parameters are used to describe the material properties of the composite on the component level.

Fig. 12 Auto-and cross-correlation of the elasticity tensor elements
Using the moving window method, the auto-and crosscorrelation is analyzed for a window size of 250 µm, 500 µm, and 750 µm. As there are closed solutions for the Fredholm integral equation given in Eq. (19) for an exponential, Gaussian and triangle correlation function the use of one of these kernels is preferable. Table 7 gives the results of the correlation lengths for these three correlation functions and the corresponding root mean square error (RMSE). The results indicate, that all correlation functions approximate the numerically obtained correlation behavior of the elasticity tensor coefficients within a RMSE of less than 1 µm. Due to the common use of the exponential correlation function in the literature [17,53,66,72], this correlation function is selected for the following generation of random fields using the Karhunen-Loève expansion.

Generation of random fields
Independent of the correlation function the Karhunen-Loève expansion always assume a linear correlation between the parameters. For the numerical model of the elasticity tensor coefficients of SFRC this assumption is satisfied not only for the auto-correlation but also the cross-correlation, as shown in Fig. 12.
One main parameter for the accuracy and convergence of the expansion of random fields is the amount of considered eigenvalues and hence, the number of terms used in the Karhunen-Loève expansion. Moreover, the required amount of terms is strongly affected by the ratio of the correlation length and the dimension of the random field [36]. Low ratios coincide with highly correlated random fields, which allows one to reduce the terms of the Karhunen-Loève expansion to represent the random field. However, for large domains or small correlation lengths and hence, an increasing ratio, more terms should be considered. Since the ratio of the correlation length and observed area is considerably large in this study 3000 terms are used during the series expansion. This meets the magnitude of terms used in [75] for the numerical simulation of the wave propagation in thin-walled structures, which is also characterized by a weak correlated random field and hence, a large ratio of domain expansion to correlation length. Figure 13 shows examples of random fields for the elasticity tensor elements, that are not equal to zero, based on a window size of 250 µm and exponential correlation functions, when assuming a two-dimensional model and fiber reinforced material. As shown in [59] the symmetry property still holds.
As mentioned before, for the representation of the elasticity coefficients Gaussian random fields are used. This may lead to negative realization. To show, that in this case the use of Gaussian random fields is the generated fields are analyzed in detail. For this Table 8 provides the mean values, the standard deviation, as well as the minimal and maximal values. For the generation of these exemplary random field the underlying distribution is not modified to avoid negative real- izations. Due to the moderate standard deviation only positive results are calculated, even if a Gaussian distribution is used for the elasticity coefficients. Furthermore, the results of this analysis may be checked for plausibility. The minimal and maximal values can be explained by areas where almost no fiber is found and by areas with a relatively high fiber volume fraction, respectively. Cases where C 11 is smaller than C 22 the fiber can be assigned to a variation of the fiber orientation.

Numerical model
The numerical model for the analysis of the material properties of SFRC is based on a tensile test specimen. Since the standard observation length defined in [38] is 50 mm the specimen is represented by a rectangle with an edge length of 50 mm × 3 mm. Vertical and horizontal displacements are fixed on the left and lower edge of the rectangle, respectively.
As done in Sect. 5.1 100 MPa are applied on the right edge. Again all simulations are carried out assuming a plain strain state as well as linear theory.
For the determination of an adequate element size a convergence study is performed. This is done by increasing the number of elements in the vertical and horizontal direction until the variation of the resulting Young's modulus is less than 0.5%. The aspect ratio of the element is not affected by the convergence analysis and therefore, is fixed with a value of 1. Table 9 gives the results for the convergence analysis based on the Young's modulus. The element size is varied between an edge length of 200 mm up to 1 mm. Therefore, the smallest element size and hence, the distance between the integration points is comparable to the smallest correlation length. Based on the obtained mean value the results indicate that an element size of 1 mm × 1 mm is sufficient for the numerical simulation on the component level.
Compared to the obtained correlation length in Table 7 the material information at the integration points are independent of each other. This means that a Karhunen-Loève expansion is not stringently necessary to derive the random fields representing the material properties. The information can also be generated by a Monte Carlo sampling at each integration point. This circumstance allows one to neglect the cross-correlation between the elasticity tensor elements, that is provided in [59]. However, the auto-correlation and hence, the random fields used in this section to represent the material properties of the SFRC are still derived by Karhunen-Loève expansion. The simulations are carried out for all window sizes and boundary condition types used in Sect. 6.1.1 to analyze the correlation structure.

Results
The results of the numerical analysis based on the microstructural characteristics derived from a two-dimensional micrograph are summarized in Table 10. As expected the values based on an analysis with pure displacement boundary conditions (KUBC) decrease with an increasing window size whereas values based on pure traction boundary conditions (SUBC) increase with an increasing window size. This meets the multi-scale modeling approach with the Voigt and Reuss bounds [34] for a SVE. The very small standard deviation can  be assigned to the element size, which has the same magnitude as the obtained correlation length.

Experimental setup
The numerical results are validated experimentally by a standard tensile test. The tensile test specimens have the dimensions of 140 mm × 20 mm × 3 mm as shown in Fig. 2(2). As depicted in Fig. 2 the specimens are taken from a plate made by mold injection. The main orientation of the specimens coincides with the melt flow direction. In total 13 specimens are analyzed. To compare the results of the numerical simulation with the experimentally obtained values, the Young's modulus is determined for the linear elastic range. Therefore, the Young's modulus is calculated for a strain range of 0.05% to 0.25%. Based on the specimen geometry this equals a load range of 150 N to 750 N and is characterized by purely linear-elastic material behavior. The measurement length for the overall Young's modulus of a specimen is 45 mm.

Results and comparison
The results for the experimentally obtained Young's modulus is given in Table 10. For a comparison with the numerical analyses the results of the modeling on the component level and their deviation to the experimental results are included as well. The results of the numerical simulation based on the probabilistic characteristics of the microstructure determined in Sect. 4.1 show a significant deviation from the experimental results of 34.3% up to 58.4%. However, the numerical data matches the results of the simulation on the mesoscale very good. The Young's modulus in fiber direction is approximated on the mesoscale with 6 GPa in Sect. 5. This meets the numerical simulation conducted here with a minimal deviation of 2.2% for the correlation structure derived from a window size of 500 µm and pure displacement boundary condition. Therefore, the main reason for the significant deviation to the numerical simulation on the component level is assumed in the use of probabilistic characteristics of the microstructure obtained by a two-dimensional micrograph. This data is sufficient for numerical simulation of the material response on the mesoscale. However, for the representation of the material properties on the component level, which is done by the generation of homogeneous second-order random fields with the Karhunen-Loève expansion, requires a different set of input variables.

Application of the procedure on values taken from literature
As discussed in Sect. 6.2.2 the results obtained by numerical simulation in Sect. 6.1.3 and the results of tensile tests presented in Sect. 6.2 do not match. It is assumed that the main reason is the lack of information about the three-dimensional micromechanical properties. To confirm this conclusion, the presented procedure is applied to threedimensional microstructural data provided in [25]. In contrast to the analyses presented so far, tensile test specimens made of PA 6 are used in [25]. However, in both cases the fiber mass fraction is specified with 30%. In a first step the correlation functions as well as the three-dimensional microstructure information of the elasticity tensor components are identified by taking into account the slightly different material parameters of PA 6 compared to PBT. The material properties provided in Table 11 are assumed dry as molded.
To generate the microstructure for the correlation function analysis the probabilistic characteristics of the fiber length, fiber diameter, and fiber orientation are taken from [25]. The characteristics of the fiber length and fiber orientation are given for tensile test specimens directly produced by mold injection. The probability density function of the fiber length meets a two-parameter Weibull distribution [7,15] given by The Weibull parameters a and b are determined as The corresponding mean of the fiber length isl = 260 µm. The fiber orientation is described by an elliptic probability density function [25]. Therefore, holds. Here, h 1 and h 2 are the semi-minor and semi-major axis, respectively. The semi axes ratio for this probability density function measured by μCT is 22.1. The fiber diameter d is analyzed for tensile test specimens cut out of a preproduced plate. Since the fiber diameter is not influenced significantly by the mold injection process these values can be combined. As expected the fiber diameter shows a normal distribution. The corresponding probability density function is written as with a mean value μ = 10.9 µm and a standard deviation σ = 0.9 µm. The derived information of the correlation structure is summarized in Table 12.
The results of the numerical simulation of a tensile test using the presented data for the representation of the microstructure are provided in Table 13. In addition, the numerically obtained Young's modulus of the material is characterized experimentally in [25] as well. Therefore, Table 12 Results of the auto-correlation lengths of the elasticity tensor coefficients using an exponential correlation functions and a window size of 250 µm based on the microstructure characteristics provided in [25] Table 13 indicate that the introduced procedure for the numerical simulation of SFRC on the component level including information of the threedimensional microstructure leads to plausible results. The numerical simulations based on random fields characterized by correlation functions obtained for a window size of 250 µm and the mean values for pure displacement boundary conditions show a minimal deviation of 1.54%. The deviation starts to increase with an increasing window size. This indicates that the use of larger window sizes leads to a loss of information about the microstructure, that is necessary to obtain the correct material response on the component level. Due to the larger aspect ratio of 23.9 the results also reveal a more distinct anisotropy in contrast to the simulation based on the microstructure.
A further validation step is the comparison of an experimentally obtained local Young's modulus distribution with the spatial fluctuation of the Young's modulus based on the modeling with random fields. Unfortunately, a local distribution of the Young's modulus is not provided in [25]. Therefore, in a last step the presented characteristics are combined with the material properties of PBT GF 30.

Transfer of the microstructural characteristics to PBT GF 30
First a correlation analysis based on the new combination of microstructural and material properties is conducted, Table 14. As the previous analysis indicates that the representation of the microstructure with random fields for a correlation structure based on a window size of 250 µm shows the best fit, the final analysis is limited to this configuration. The results of the numerical and experimental analysis, given in Table 15, confirm the presented procedure. The deviation between the numerically and experimentally obtained values is less than 1%. Furthermore, the obtained minimal and maximal values show he same trend for the numerical simulation as well as the experiments. These values indicate the Young's modulus range for the analyzed specimens.
However, a further important property is the local distribution of the material property within one specimen. Therefore, the Young's modulus distribution over the observation length is measured by dividing the original measurement length into sections with a length of 5 mm. For each section the Young's modulus is obtained individually on a numerical and experimental basis. The results are given in Table 16. The mean values of the chosen data sets fit very well. In contrast to this, the standard deviation of the experimentally obtained data set is higher than the one obtained for the numerical results. This is also indicated by the range of the values. The corresponding minimal and maximal values are marked in Table 16. The reason for this is determination of the correlation length on the mesoscale in combination with the used element size for the numerical simulation.

Summary and conclusion
The presented work consists of three steps that build together the proposed modeling approach for a tensile test specimen made of SFRC. These are the analysis of the microstructure, the determination of the correlation structure of the elasticity tensor elements and the generation of random fields for the representation of the material properties by the Karhunen-Loève expansion on the component level.
The analysis of the microstructure is based on a twodimensional micrograph which reveals a layered structure of the cross-section. Based on the data the cross-section can be divided into five layers, namely a core layer, that is surrounded by a shell and a skin layer at the top and bottom. This meets the well-known orientation scheme due to mold injection process. Numerical simulation based on the derived probabilistic characteristics of the fiber length, fiber diameter, fiber orientation as well as fiber volume fraction show that the information about the layered structure influences the Young's modulus of the analyzed cross-section. A validation with experimentally obtained Young's modulus by nanoindentation indicates, that a homogenization based on the overall probabilistic characteristics underestimates the material properties of the cross-section. In contrast to this the modeling based on a layered structure consisting of the core and shell layers shows only a marginal deviation to an Table 15 Results of the numerically obtained Young's modulus based on the microstructural characteristics of PA 6 given in [25]  approach limited to the shell layers. Therefore, the representation of the microstructure by the shell layers only appears sufficient and is used for the determination of the correlation structure for the numerical simulation on the component level.
For the numerical simulation on the component level first the probabilistic characteristics based on a two-dimensional micrograph of the shell layers is used to derive the correlation structure. However, a comparison with conducted tensile tests show large deviations between the numerically and experimentally obtained results for the Young's modulus. In fact, the numerical results meet the Young's modulus obtained by nanoindentation instead. This leads to the conclusion that the numerical modeling on the component level requires a full three-dimensional analysis of the microstructure as a tensile test is no longer reduced to near surface measurements. This conclusion is proven by applying the presented approach to probabilistic characteristics of a microstructure of SFRC obtained by μCT. The resulting minimal deviation of the numerical simulation and the conducted experiment is about 1.5%. Finally, the local distribution of the material properties is investigated both numerically and experimentally. The results show the same trend. However, the spatial fluctuation seems to be more prominent experimentally than for the numerical simulations. One reason might be the used correlation lengths, that are smaller than the element size. This issue could be solved by deriving the correlation length of each elasticity tensor component not by the moving window method in the mesoscale. Instead the correlation lengths give additional independent parameters like the engineering constants. The correlation lengths are then determined in the same way as the material properties by fitting the parameters to experimental data. This approach is already successfully used in [74,75], where the correlation length of the Young's modulus is determined in such a way that the phenomenon of the continuous wave mode conversion, detected by the experimental observation of guided ultrasonic waves in carbon fiber reinforced plastics, can be observed in numerical simulations as well. This approach will be combined with the use of non-Gaussian random fields.
It is concluded, that the presented procedure is suitable for the numerical simulation of components made of SFRC. For an adequate representation on the mesoscale the probabilistic characteristics of a two-dimensional cross-section is sufficient. The numerical simulation on the component level, however, require a knowledge of the three-dimensional probabilistic properties of the microstructure.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.