Fiber orientation interpolation for the multiscale analysis of short fiber reinforced composite parts

For short fiber reinforced plastic parts the local fiber orientation has a strong influence on the mechanical properties. To enable multiscale computations using surrogate models we advocate a two-step identification strategy. Firstly, for a number of sample orientations an effective model is derived by numerical methods available in the literature. Secondly, to cover a general orientation state, these effective models are interpolated. In this article we develop a novel and effective strategy to carry out this interpolation. Firstly, taking into account symmetry arguments, we reduce the fiber orientation phase space to a triangle in R2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\mathbb {R}}^2$$\end{document} . For an associated triangulation of this triangle we furnish each node with an surrogate model. Then, we use linear interpolation on the fiber orientation triangle to equip each fiber orientation state with an effective stress. The proposed approach is quite general, and works for any physically nonlinear constitutive law on the micro-scale, as long as surrogate models for single fiber orientation states can be extracted. To demonstrate the capabilities of our scheme we study the viscoelastic creep behavior of short glass fiber reinforced PA66, and use Schapery’s collocation method together with FFT-based computational homogenization to derive single orientation state effective models. We discuss the efficient implementation of our method, and present results of a component scale computation on a benchmark component by using ABAQUS ®.


State of the art
Composite materials are frequently used in engineering applications. However, the difference in sizes between the reinforcements and the parts generally prohibits a finite element analysis of the component where the underlying mesh resolves the heterogeneities. It is common in engineering practice to overcome this problem by resorting to effective models which characterize the composite's behavior on the component scale.
To identify these models, experimental methods are supported by mean field and computational upscaling techniques (cf. Zaoui [61] and Matouš et al. [35] for respective surveys), mathematically formalized by the theory of homogenization [6].
For nonlinear material behavior or large differences in the elastic properties of the constituents, numerical methods resolving the microstructure on representative volume elements (RVEs) for the mechanical behavior under considerations [28] are often the only accurate option, and a variety of numerical approaches specialized to homogenization problems have been developed [2,5,40]. Solving two-scale problems with nested volume element problems, known as FE 2 analysis or heterogeneous multiscale method [17], is limited by their formidable sizedespite impressive computational results like in Mosby-Matouš [38]. The micro-mechanical problems under consideration, however, share the same geometry, but involve different macroscopic loading conditions. Thus, a huge amount of very similar problems is solved. Based on this observation a variety of acceleration techniques have been developed, among them reduced order models [9,55,59], the (non-uniform) transformation field analysis [16,20,36], response surface models [4,12,46,50,60] and machine learning approaches like neural networks [21,24,33,34,52] giving rise to accurate surrogate models with significantly lower computational demands.
Surrogate models can than be used in a database concept [49,51].

Problem setting
Short fiber reinforced composite parts constitute an example for multiscale problems where the parameters determining the microstructure may vary continuously on the macro scale.
For illustration, Fig. 1 shows a quick release buckle socket 1 for which an injection molding simulation was conducted (see Sect. 4 for the exact simulation parameters). The colors indicate the occurring fiber orientation (FO) see Fig. 2: cyan corresponds to isotropic, yellow to planar and magenta to a unidirectional, i.e. an aligned, fiber orientation state. We see that large parts of the socket exhibit a rather aligned orientation. Yet, there are distinguished parts, like the flat stern where the orientation is mostly planar. Furthermore, there are regions where the melt fronts have collided (welding lines). Capturing the mechanical behavior of these regions is particularly important, as they correspond to weak spots of the structure.
The finite element mesh, cf. Fig. 1a, consists of more than 300,000 elements. This large number of elements is necessary for resolving the component. Notice that the fiber orientation is varying continuously, and large parts of the component are characterized by very similar fiber orientation states. Thus, it appears natural to employ an interpolation ansatz to cover the different orientations. Unfortunately, using representative volume elements (RVEs) for computing the effective mechanical response is, at first sight, incompat-ible with an interpolation scheme: the volume element does not depend continuously on the fiber orientation! The situation is even more delicate: one and the same fiber orientation state might be represented by two completely different volume elements. Thus, the association {FO} → {RVE} is not even a well-defined function. Similar problems arise if surrogate models of different type, in particular with different internal variables, are used for different orientations. This situation is characteristic for models obtained by model order reduction, for instance.
In this work, we propose a work-around for these problems, based on the insight that the effective stress is a well-defined function of the fiber orientation (and the applied load history).

Our contribution and organization of this article
Within the framework of generalized standard materials [23,25] we propose a method to interpolate arbitrary effective models. Suppose we have a collection of microstructures, indexed by {1, . . . , M}. For every microstructure whose constituents are generalized standard materials, the effective mechanical behavior can be described by a generalized standard material, see Suquet [48]. Thus, we obtain a free energy W i and a dissipation potential Φ i for each index i. Notice, however, that Φ i depends also on an internal variable vector Q i , which lives in some space X i .
Let the state to be interpolated be described by weights s 1 , . . . , s M forming a convex combination, i.e. satisfying the interpolated free energy and dissipation potential read where E denotes the applied macroscopic strain and Q = (Q 1 , . . . , Q M ) collects all internal variables. This construction naturally preserves convexity (and continuity) properties of the potentials, and yields the desired evolution equations with effective properties continuously varying between the microstructures. For a typical simplicial triangulation in the d-dimensional space describing the microstructure, each microstructure to be interpolated is surrounded by d + 1 nodes. Consequently, only d + 1 internal variables need to be tracked instead of the full M. This property contrasts with other methods, for instance the method of pseudo-grains [13,15].
In the context of short fiber reinforced composites we apply the proposed method to the interpolation of fiber reinforced structures described by the fiber orientation tensor of second order, see Sect. 2, corresponding-if objectivity is accounted for-to a two-dimensional configuration space.
To study the creeping behavior of short fiber reinforced composites, we equip the matrix with a Burgers model for PA66 [57,58], use FFT-based computational homogenization [39,40] to carry out the microstructure computations and identify a surrogate model by Schapery's collocation method [43], see Sect. 3. Last but not least we check the accuracy of the fiber orientation interpolation and investigate the behavior of the quick release buckle socket of Fig. 1, see Sect. 4. Care has been taken to ensure applicability of the scheme in state of the art engineering computations. For the micromechanical simulation, we use the Fraunhofer ITWM software FeelMath [18], distributed as the ElastoDict module of GeoDict [22]. The database containing the different effective models is filled by a script written in Python [42], and can be accessed by a user-defined material function, which is compatible with the UMAT interface of ABAQUS [1].

The fiber orientation tensor
Short fiber reinforced composites with fibers of equal length are locally characterized by the direction of their fiber reinforcements. Suppose a given volume V comprises a number of fibers with directions p 1 , . . . , p N ∈ S 2 , where S 2 denotes the the unit sphere, N can be very large. Then, the Advani-Tucker (second order) fiber orientation tensor A [3] computes via The fiber orientation tensor A is a symmetric positive semidefinite 3 × 3 tensor with tr (A) = 1 and represents a compact measure of the fiber orientation distribution in the volume V . For instance, A is degenerate precisely if all fibers { p i } N i=1 are contained in a single great circle, i.e. if the fiber orientation is planar. Furthermore, A is of the form A = p⊗ p for a direction p if and only if all fibers p i point in direction p (or, equivalently, in direction − p).
A carries a limited amount of information, as it comprises only five independent degrees of freedom. Still, it is the principal quantity of interest for commercial injec-tion and compression molding simulations, as it may vary for every Gauss point of the FEM Mesh for processed part and plausible higher order momenta are recovered by closure approximations (see Montgomery-Smith et al. [37] for a recent overview). For extensions to higher moment tensors compare Jack-Smith [29]. A principal component analysis of A reveals further information on the fiber orientation. Spectrally decomposing A leads to the expression with an orthogonal matrix R and a diagonal matrix Λ with sorted eigenvalues λ 1 ≥ λ 2 ≥ λ 3 . The columns r 1 , r 2 and r 3 of R represent the principal fiber directions, and the eigenvalues λ i describe the probability of finding fibers in direction r i , as the λ i are non-negative and sum to one. In particular, up to an orthogonal transformation, the fiber orientation tensor A can be described by by two positive real numbers λ 1 and λ 2 satisfying the two inequalities which geometrically corresponds to a planar triangle, cf. Fig. 2.

The framework of generalized standard materials
For the constitutive modelling we rely upon the two-potential framework of the generalized standard materials [23,25]. The constitutive relationships are derived from two thermodynamic potentials, the Helmholtz free energy density w(ε, q) and the dissipation potential φ(q) which are convex functions of the state variables, the infinitesimal strain ε and other internal variables q, and their time-derivative. The free energy w is the energy available in the system to trigger its evolution, whereas φ describes the evolution of the irreversible mechanisms.
In a first step, the thermodynamic forces stored in the materials are derived via Then, complementary laws relate the rate of the state variables and the forces acting upon them, i.e.
Combining these formulae leads to the equations

Effective properties of short fiber reinforced plastics and the invariance principle
From the work of Suquet [47,48] it is known that the macroscopic response of a microstructure with generalized standard material constituents is described by a generalized standard material, possibly with an infinite number of internal variables. Thus, we assume, that for every fiber orientation tensor A we have an effective free energy W A , depending on the macroscopic strain E and a vector of internal variables Q A , which, for simplicity of exposition, we assume to be a collection of two-tensors, and an effective dissipation potential Φ A , depending on the rate of change of the internal variablesQ A . Thus, the macroscopic constitutive relations according to (4) read and where Σ is the macroscopic stress. However, the material laws for different fiber orientation states are not independent. Due to objectivity and material frame indifference, both the free energy and the dissipation potential are invariant w.r.t. orthogonal transformations. More precisely, for any fiber orientation A, strain E, internal variable Q A and orthogonal transformation R we have the equalities and where we used that the fiber orientation is static,Ȧ = 0. Differentiating yields the equations and which are equivalent to (5) and (6) . The latter observations are illustrated on a unidirectional structure in Fig. 4. To compute the effective stress of the unidirectional structure on the left, which can be directed arbitrarily, we can use a unidirectional structure with a fixed direction (on the right), provided we transform the applied strain as well as the internal variables (not shown) accordingly, and do not forget to transform the effective stress back to the reference orientation. Thus, to derive effective models for all fiber orientation states A it is sufficient to consider the fiber orientations Λ = diag(λ 1 , λ 2 , 1 − λ 1 − λ 2 ) with λ 1 , λ 2 lying in the triangle (3). Effective constitutive laws for general orientation tensors are obtained by subsequent orthogonal transformations.
More precisely, suppose the family {W Λ , Φ Λ } Λ is given, and let A be arbitrary. Spectrally decomposing A = R T ΛR leads, for macroscopic strain E, to the constitutive law The free energy, the dissipation potential and the internal variables live in the diagonalized configuration, whereas both stresses and strains share the frame of the fiber orientation tensor A.

Fiber orientation interpolation
In the preceding section we have reduced the problem of finding the effective response of a fiber reinforced composite from a five-dimensional problem, corresponding to the full fiber orientation tensor A, to a two-dimensional problem, taking into account the fiber orientation triangle. Still, this triangle consists of an infinite number of points. From physical intuition it is clear that for fixed applied strain the free energy in equilibrium is continuous as a function of the orientation. Thus, if we discretize the fiber orientation triangle of Fig. 2 by a sufficient fine mesh, and associate a constitutive law (W Λ , Φ Λ ) to each nodal point Λ of the mesh, a subsequent interpolation should capture all effective materials law with sufficient accuracy, provided the triangulation is fine enough, see Fig. 3.
However, care has to be taken as to which quantities can be interpolated. In general, in contrast to macroscopic strains and stresses, it does not make sense to interpolate internal variables, as they live on the current configuration (and fiber orientation state). Thus, it is not clear how to transfer them between different states. Even more striking is the fact that upon discretization the number of internal variables becomes finite, and might differ for distinct fiber orientations. A solution to this dilemma is the following. Suppose Λ is a general fiber orientation state in the fiber orientation triangle. For a given triangulation let Λ 1 , Λ 2 and Λ 3 be the nodes of the triangle containing Λ, with constitutive laws with a unique real triple (s 1 , s 2 , s 3 ) satisfying We associate to Λ the free energy and the corresponding dissipation potential These two potentials include three sets of internal variables, and are continuous in Λ through the (s 1 , s 2 , s 3 )coefficients (12). Our construction preserves the convexity of both the free energies and the dissipation potentials.
Differentiating (14) according to (4) yields the constitutive law Thus, the nodal stresses are interpolated, and the nodal internal variables evolve independently of each other. Of course, if s i = 0 for some i = 1, 2, 3, the evolution for Q Λ i becomes irrelevant for the evaluation of the stress. In particular, if Λ = Λ i for some i = 1, 2, 3, we consistently recover the constitutive equations for Λ i .
linear interpolation using fiber orientation database stress response In this formulation, three sets of internal variables are stored. Each node receives the same strain, but uses its own constitutive law to produce a stress, and the evolution of the internal variables depends on the node only. The three stresses are averaged according to the convex representation (12), which is meaningful, see Fig. 5.
We have presented the method for a single triangle of the triangulation only. However, the method can be interpreted for the whole triangulation in a similar fashion, where only the adjacent nodal constitutive laws are active. Notice at this point that for every orientation only three material laws are required, regardless of the number of nodes used for the triangulation of the fiber orientation triangle. This property contrasts starkly with the method of pseudo-grains [13,15], which requires the full number of grains to be evaluated for calculating the effective response of the composite.
The fiber orientation interpolation method can be summarized as follows.
1. Offline phase: triangulate the fiber orientation triangle, cf. Fig. 3, and compute an effective model for each node in the triangulation 2. Online phase: for each macroscopic Gauss point and given strain E -Spectrally decompose the local fiber orientation ten- the action of R on two-tensors R i jkl = R ik R jl , and In this work, the offline phase is implemented in Python [42]: for a given triangulation of the fiber orientation triangle, the micromechanical solver (see Sect. 3.3) is called for each node of the triangulation and its results are used to identify the parameters of a surrogate model. As output, for each element in the fiber orientation triangulation, an ordered set of nodes with corresponding ABAQUS [1] UMATs and identified parameters for the UMAT are stored in an XML file. Once the component to be simulated is chosen, for each element of the component's finite element mesh, the eigendecomposition of the fiber orientation is computed, the Euler angles of the rotation matrix, the weights associated to the interpolation (12) and the corresponding element of the fiber orientation triangle are determined and stored.
For the online ABAQUS computation a dummy UMAT is used, which takes Euler angles, weights and the fiber orientation element as input parameters, and calls the three UMATs associated to the adjacent nodes of the fiber orientation triangle.

Identifying a viscoelastic surrogate model using FFT-based computational homogenization and Schapery's collocation method
To test the orientation interpolation technique we consider a fiber reinforced polymer with linearly viscoelastic matrix and linear elastic reinforcements. Linear viscoelasticity has the advantage that on the one hand the derivation of effective models is comparatively well-understood, but on the other hand due to the dependence on the material history complex constitutive behavior can be observed.

Generating fiber-filled volume elements
To generate periodic volume elements with prescribed fiber orientation and volume fraction we rely upon the Sequential Addition and Migration (SAM) method [44]. As input parameters, the fiber volume fraction φ, the second order fiber orientation tensor A, the fiber length L, the fiber diameter D and the edge length L x (= L y = L z ) of the cubical volume element needs to be specified. In a preliminary step, the required number N of fibers is calculated, s.t. the fiber volume fraction φ is matched with highest precision. Furthermore, the fourth order fiber orientation tensor A is computed from the A tensor with the help of the exact closure approximation [37].
In the first step of SAM, N overlapping cylindrical fibers are randomly placed in the volume. Then, in the second step, the fibers are displaced and rotated, s.t. the overlap is removed The SAM algorithm is unable to generate purely planar fiber orientation states, i.e. (up to orthogonal transformation) states corresponding to a fiber orientation tensor of the form However, these states do not appear in practice for industrial fiber volume fractions. Indeed, suppose a collection of fibers with directions p 1 , . . . , p N is given. This system is purely planar in the x-y-plane precisely if p i · e z = 0 for all i = 1, . . . , N , i.e. all fibers have vanishing z-component. In particular, no scatter is allowed in the z-direction.
For the study at hand, we force λ 3 to be at least 0.01. These orientation tensors are sufficiently close to planar for practical studies, and can easily be generated by the SAM method.

The Burgers model for the viscoelastic behavior of the matrix
As matrix material we consider a commercial polyamide 6.6 (PA66, DuPont Zytel 101), for which experimental creep curves and fitted viscoelastic model parameters are available in the literature, see Yang et al. [57,58]. The Burgers model, a classic in the literature [53], is widely used for the linear viscoelastic behavior of ideal thermoplastics. The rheological diagram is shown in Fig. 6 with a Maxwell and a Kelvin unit connected in series. The Burgers model exhibits a stationary creep rate after relaxation of the Kelvin element, includes an irreversible part of the viscous strain after unloading and the stress in a relaxation test approaches 0 as t ∞. According to Burgers' model, for an applied tensile stress of σ 0 , the tensile strain computes as where E M and η M are the Young's modulus and viscosity of the Maxwell spring and dashpot, respectively, and E K and η K are the Young's modulus and viscosity of the Kelvin spring and dashpot, respectively, whereas τ = η K /E K represents the retardation time.
Reciprocally, for an applied tensile strain of ε 0 , the tensile stress attains the form [19] Using the latter representation it is not difficult to exhibit Burgers' model as a generalized standard material. However, for practical purposes the creep function representation (16) is advantageous. Following Lai-Bakker [32] and Woldekidan [56] we obtain a 3D model by assuming Poisson's ratio ν = 0.38 to be constant over the complete time history. This is a common assumption for glassy or semi-crystalline thermoplastic polymers as the investigated matrix material. We arrive at the hereditary integral formulation relating the stress rateσ to the strain ε using the viscoelastic intensification tensor for t ≥ 0 and zero otherwise. Here, F denotes the creep compliance and D 0 is the instantaneous compliance. If F = 0, then D 0 + D 1 corresponds to the relaxed (elastic) compliance. In formulae, the action of the three compliances on a test stress S read where Id is the 3 × 3 identity matrix.
The chosen material parameters are summarized in Table  1, corresponding to the σ 0 = 20 MPa and T = 23 • C fitted parameters of Yang et al. [58].

Computational homogenization of linear viscoelasticity
Following Hashin [26,27] suppose Y ⊆ R 3 is a cuboid (for periodic boundary conditions) and let a heterogeneous linear viscoelastic medium be given in terms of the local relaxation function -compatibility -linear viscoelastic material law in relaxation function form -fixed macroscopic stress Here, ε is the local strain field, we write σ (t) ≡ σ (·, t) for the sake of notational clarity, ∇ s denotes the symmetrized gradient, and · Y stands for the averaging operator over Y.
Without loss of generality, u(t) Y = 0 can be enforced for every t ≥ 0. The quantity of principal interest is the macroscopic strain field history E. Due to the linearity and the time-invariance of the mapping Σ → E, there is a homogenized relaxation function J hom : [0, ∞) → R 3 ⊗ R 3 ⊗ R 3 ⊗ R 3 with both minor and major symmetries, s.t. we can write Indeed, J hom arises by choosing for six linearly independent Σ 0 ∈ Sym(d), as the corresponding macroscopic strain fields satisfy Thus, the effective creeping behavior completely determines the full linear viscoelastic effective behavior.
We solve (18)-(21) in the standard way [30] by resorting to a time discretization 0 = t 0 < t 1 < t 2 < · · · , s.t., for every time step, (18)-(21) becomes a linear elastic problem with eigenstrain, i.e. we seek, suppressing time indices, a strain E ∈ Sym(d) and a periodic displacement fluctuation u : Y → R 3 with zero mean, s.t. the equations div σ = 0, are satisfied, where D tan is a tangential compliance tensor and α stands for an eigenstrain depending on quantities computed at the previous time step. (23) constitutes a standard linear elastic homogenization problem with eigenstress and stress "boundary conditions", for which a variety of solution techniques are available.
FFT-based computational homogenization [39,40] constitutes our method of choice. For this method, a discretization on a regular grid or mesh decomposing Y is chosen. Due to the regularity of the grid and the periodic boundary conditions, problems of the form (23) with homogeneous D tan can be solved directly in terms of the discrete Fourier transform. The Operator of this auxiliary homogeneous problem serves to precondition the linear system (23).
For the problem at hand, we choose a discretization on a staggered grid [45], the stress-based formulation of Bhattacharya-Suquet [11] and the conjugate gradient method [62] for the resolution of (23). All mentioned methods are integrated into the Fraunhofer ITWM C++ code FeelMath [18].

On the resolution and the RVE size necessary for the precomputations
In this section we investigate the necessary resolution and representative volume element (RVE) size for solving (18)- (21), or, equivalently, (23), to engineering accuracy. We consider short E-glass fiber reinforcements with elastic parameters E = 72 GPa and ν = 0.22, length of 200 µm, aspect ratio r a = 20 and fiber volume fraction φ ≈ 17%, which corresponds to a fiber mass fraction of 30% in the PA66 matrix with the mechanical properties described in Sect. 3.2.
Recall from (17) that for creeping with stress σ 0 , Burgers' model predicts a strain in the matrix, whereas the glass fibers behave linear elastic (in particular, time independent) To investigate the viscoelastic behavior, we consider two extreme cases:

Case 1
The instantaneous response, corresponding to t = 0, i.e. (24) becomes simply the linear elastic relationship ε = D 0 : σ 0 Case 2 The creep rate at infinity. Differentiating (24) in time yieldṡ i.e., lim t→∞ε (t) = F : σ 0 , whereasε = 0 in the fibers. To study the creep rate at infinity, we study the linear elastic problem with local compliance D(x) = F, x in the matrix, 0, in the fiber.
As the fibers are rigid, this problem is much more difficult to solve than the instantaneous elastic problem.
As a first test we study the resolution necessary to obtain accurate effective properties, both for the instantaneous elastic case and the creep rate at infinity. For this purpose, we consider cubical elements with an edge length of L x = L y = L z = 600 µm, corresponding to thrice the fiber length L = 200 µm. We choose the extreme orientations, i.e. unidirectional (ud), planar isotropic and (3D) isotropic.
We have chosen to check four different resolutions, measured in voxels per fiber diameter. Recall that the fibers have a diameter of 10 µm, so that 5 voxels per diameter correspond to a voxel edge length of 2 µm, i.e. the (600 µm) 3 volume is discretized by 300 3 elements. Accordingly, 10, 15 and 20 For each of these resolutions, the homogenized instantaneous elastic tensor C hom was computed (Case 1). The relative error in the Frobenius norm of the corresponding Voigt matrices, measured relative to the highest resolution 1200 3 , is shown in Fig. 7a. Even for the crudest resolution, the relative error for all three considered orientations was below 1%.
The same computation was carried out for the long-term, i.e. viscous, response (Case 2). Upon replacing the strain ratė ε by the strain ε in (26) the problem is formally equivalent to an elasticity problem with rigid fibers and elastic matrix. The computed homogenized elastic tensor for the equivalent elastic problem is in fact a homogenized viscous tensor V hom , and the relative errors of the results are plotted in Fig. 7b, again relative to the highest resolution.
For the crudest resolution, the relative errors exceed 10% for all three considered orientations. The highest error of about 19.1% occurs for the planar fiber orientation state, and the lowest error (11.7%) is reached for the isotropic orientation. Doubling the resolution decreases the error for all three orientations below 2.5%. For 900 3 , the relative error is even smaller than 1%.
We see that in contrast to the instantaneous elastic case computing the long term viscous response requires a comparatively high resolution of 10 voxels per fiber diameter, i.e. a resolution of 1 µm. We fix the latter resolution for the succeeding. Next we determine the size of a volume element to be representative (see Bella and Otto [7] for recent mathematical findings in the context of linear elasticity) for the homogenized viscoelastic response, see Sect. 3.3, measured in terms of the fiber length L = 200 µm. With the resolution of 1 µm, the investigated volume elements comprise 400 3 , 600 3 , 800 3 , 1000 3 and 1200 3 voxels. To get an impression about the relative size of the volumes, compare Fig. 8, where volume elements with isotropic fiber orientation are shown. The number of fibers corresponding to each volume element size is listed in Table 2.
To study the size of an RVE, we compute the effective elastic/viscous tensors for the different sample sizes and compare these to the result obtained for the largest sample L x = 6× L. It should be mentioned that the variance of different realizations for a fixed unit cell length is negligible, as both the fiber volume fraction and the fourth order fiber orientation tensor are fixed with high accuracy by the SAM method, cf. Schneider [44]. Hence, only the systematic error is measured. For the instantaneous elastic case, see Fig. 9a, the relative error for all resolutions and orientations considered is well below 0.5%, and can be considered negligible. In particular, for the elastic response, twice the fiber length is a sufficient edge length for a representative volume element.
In the viscous case, see Fig. 9b, the relative errors are about an order of magnitude larger than for the elastic case. The unidirectional fiber orientation leads to the largest errors, whereas the isotropic state induces the lowest errors, with the planar case in between these extremes. For the smallest edge length, the relative error is still comparatively large, with about 4% for the unidirectional fiber orientation. We consider an edge length of thrice the fiber length as sufficient, since both the planar isotropic and isotropic fiber orientation states lead to a relative error of about 1%, and the unidirectional case is captured with < 2.5% relative error. Table 3 lists the number of iterations necessary for resolving the unit (a)  cell problems for different volume element sizes, both for the elastic and the purely viscous response, measured as the average of the six load cases. The dual formulation of Bhattacharya-Suquet [10] was solved by the conjugate gradient method up to a relative residual of 10 −5 . The iteration count mainly depends on the material contrast, is almost independent of the volume element size, but depends slightly on the orientation. For the unidirectional state, a lower number of iterations suffices compared to the planar and 3D isotropic structures. The viscous computations increase the iteration count almost by a factor of 20 compared to the elastic computations.

Schapery's collocation method
To obtain a surrogate model for the homogenized linear viscoelastic behavior of a heterogeneous Burgers medium we rely upon Schapery's collocation method [43]. Recall that the heterogeneous linear viscoelastic material law reads, for every microscopic point x ∈ Y ,  where, for Burgers' model, the creep tensor has the form Consider, for some fixed b > 0 and a positive integer n, the ansatz where D hom , F hom and B hom k are possibly anisotropic fourth-tensors to be determined. (28) are homogenized, giving rise to homogenized tensors K hom l . Finally, the tensors B hom k are eliminated from the 2n + 1 tensor equations

On the number of collocation points for Schapery's method
In this section we study the applicability of Schapery's method in our context.
14 Viscous orthotropic Young's moduli for different fiber orientations. a Young's modulus surface for the discretized fiber orientation triangle. b Young's moduli plotted along the corners of the fiber orienta-tion triangle (piso, iso, ud corresponds to the planar isotropic, isotropic, unidirectional fiber orientation) With the previously established parameters, an RVE edge length of three times the fiber length and a resolution of 1 µm, we have computed solutions to the Schapery problems (30) for a different number of collocation points, ranging from 1 to 7, and the basis b = 10.
To assess the quality of Schapery's method for different collocation points, we have conducted a creep simulation with 1 MPa stress applied in x-direction, with logarithmically distributed time steps up to 10 4 h. In Fig. 10 the creep strains computed from Schapery's model (28) are compared to the full resolution computation for the three extreme orientations.
For the isotropic fiber orientation state, see Fig. 10a, a single collocation point starts to strongly deviate from the reference curve at about 40 h, whereas 3 collocation points start to deviate at about 120 h. 5 collocation points slightly overestimate the creep strain, whereas the stress-strain curve for 7 collocation points is hard to distinguish from the reference solution.
The behavior for the unidirectional and the planar isotropic geometries are similar to the isotropic case.
To gain further insight into the quality of the Schapery approximation we consider the maximal relative error in the strain for the different collocation methods, computed  16 Viscous orthotropic shear moduli for different fiber orientations. a Shear modulus surface for the discretized fiber orientation triangle. b Shear moduli plotted along the corners of the fiber orienta-tion triangle (piso, iso, ud corresponds to the planar isotropic, isotropic, unidirectional fiber orientation) relative to the full-field solution. Here, maximal means the maximum error up to 10 4 h of applied stress. Figure 11 shows that there are strong differences in the error, depending on the direction of applied stress and the microstructure. The relative error is smallest for the unidirectional structure and tension in transverse direction, remaining below 1% even for a single collocation point. This is not unexpected, as this loading scenarios is dominated by the matrix response, and the matrix is described by Burgers' model, which features only a single collocation point. The situation in transverse direction strongly contrasts to extension in fiber direction for the unidirectional structure. For up to 3 collocation points, the relative error is larger than 10%, but decreases quickly for increased number of collocation points, and remains below 1% for seven collocation points. The behavior of the planar isotropic structure, loaded transversely, is somewhat similar to the unidirectional structure loaded in fiber direction. Also, there are similarities between the isotropic fiber orientation and the planar isotropic microstructure loaded in fiber direction: for these scenarios, the errors are largest, and decrease rather slowly. A relative error below 1% is reached for at least 9 collocation points. All in all we see that a relative error below 1% is guaranteed for more than 9 collocation points.

Fig. 17
Relative interpolation error of the instantaneous elastic orthotropic constants.
As a side remark notice that the error in Fig. 11 increases for more than 13 collocation points. This is no mistake, as these frequencies have little impact for creep up to 10 4 h, but change the optimum in the minimization (29) which covers also much larger times t ∈ [0, ∞). To further increase the accuracy, the base b would have to be changed from its current value 10. However, for the problem at hand, we consider the accuracy to be sufficient.
Last but not least we take a look at the computational costs of the simulations, which were carried out on 64 MPI compute nodes, equipped with a dual Intel Xeon E5-2670 and 64 GB RAM, on our cluster. The runtimes are plotted in Fig. 12 as a function of the frequency τ b k with τ = 14.2 h and varying integers k, see Table 1.
For k = 0 a runtime of 3.08 min is required. Then, as k decreases, the time slightly increases up to the elastic case "at infinity" requiring 3.5 min. For increasing k, the runtimes also increase compared to k = 0, at about 4 min for k = 1 and about 6 min for k = 2. For k = 4, already about 43 min are necessary for resolving the unit cell problem. Thus, the viscous computations are the most time-demanding, confirming the findings of Sect. 3.4.

Discussion of the results and errors
The database was filled with the previously determined parameters, i.e. an RVE length of thrice the fiber length, a resolution of 1 µm and 9 collocation points for Schapery's  Fig. 18 Relative interpolation error of the viscous orthotropic constants. a E v method, and a regular triangulation of the fiber orientation triangle with 15 nodes, as shown in Fig. 3. The entire computation took about 15 h on 64 MPI compute nodes, equipped with a dual Intel Xeon E5-2670 and 64 GB RAM, on our cluster. As the prescribed fiber orientation tensors of fourth order are orthotropic, the effective tensors K l , and thus, the tensors D hom , B hom k and F hom from (28) are orthotropic as well, see Schneider [44]. Hence it is reasonable to investigate the corresponding orthotropic engineering constants.
First, in Fig. 13, we take a look at the orthotropic Young's moduli of the instantaneous elastic response, i.e. those which can be read from D hom . Figure 13a shows the three Young's moduli E 1 , E 2 and E 3 as graphs plotted over the fiber ori-entation triangle, whereas Fig. 13b restricts the view to the edges of the fiber orientation triangle. The distances between ud (unidirectional), isotropic (iso) and planar isotropic (piso) on the x-axis of Fig. 13 correspond to the lengths of the edges connecting the corners of the fiber orientation triangle. We see that the graphs are smooth, but non-linear, and two of the three Young's moduli coincide on the edges connecting udiso and iso-piso, as they correspond to transversely isotropic fiber orientations. The contrasts between the Young's moduli is largest at the ud state (by about a factor of 2), much smaller at the piso state (about 4 3 ), and of course, minimal for the isotropic state.
For comparison, the "viscous" Young's moduli, i.e. those which can be read from F hom , are plotted in Fig. 14. The qual-itative behavior of the graphs is similar to the instantaneous Young's moduli. However, the contrasts are much larger. Indeed, for the ud case, the Young's modulus in fiber direction E v 1 and the transverse Young's modulus E v 2 = E v 3 differ by a factor of roughly 18, and the contrast for the piso case computes as about 5. Thus, the creeping behavior of a composite exhibits a much higher anisotropy than the elastic response.
Next we investigate the instantaneous shear moduli G 12 , G 13 and G 23 in Fig. 15. In this case, the differences between the different shear moduli is largest for the piso structure, with a contrast of about 1.5. Overall, the graphs appear to be at least quadratic.
Last but not least we take a look at the "viscous" shear moduli, which, as for the Young's moduli, exhibit qualitative similarities to their instantaneous counterparts, but differ strongly quantitatively. Indeed, the contrast between the shear moduli G v 12 and G v 13 = G v 23 at the piso state is about 5 (Fig. 16).
We have shown piecewise linear approximations of the orthotropic constants, considered as functions on the fiber orientation triangle. It remains to assess the quality of the approximation. For this purpose we have repeated the computation of effective elastic and viscous tensors D hom and F hom for the fiber orientations corresponding to the centers of mass of the triangles of our chosen triangulation and compared those tensors to their counterparts arising by interpolation. We investigate the relative error of the orthotropic elastic constants of the instantaneous elastic response in Fig. 17. The overall level of the relative error is rather small, bounded above by 0.7%.
For the orthotropic constants of the viscous tensor F hom the relative errors arising from interpolation are depicted in Fig. 18. The largest error of about 6% occurs for Poisson's ratio ν 31 . For the Young's moduli E v 1 , E v 2 and E v 3 the relative error is smaller than 2% except for the tip of the fiber orientation triangle and E v 2 , corresponding to the vicinity of the ud orientation, where the relative error rises to about 4%. This results from the strong contrast between the in fiber Young's modulus and the transverse Young's moduli in the viscous case, see Fig. 15b. The relative errors for the shear moduli are bounded by 4%.
All in all, the interpolation scheme is very accurate, except for some extreme cases (the Poisson's ratio in the viscous case in the vicinity of the ud state). Still, we consider the maximum error of about 6% acceptable since in most other cases, the relative errors are well below 5% relative error and thus within engineering accuracy.

A computational example
To prove our fiber orientation interpolation concept we investigate the viscoelastic creep behavior of a quick release  buckle, whose CAD geometry is publicly available, 2 see Fig. 1a We conducted a mold filling simulation with injection molding software InjectionMoldingFoam [41], which is based upon OpenFOAM's incompressible two-phase flow solver [54]. Parameters for mold-filling release buckles are available in the literature, see Bhat et al. [8]. The parameters we used are summarized in Table 4 partially corresponding to the following variant of the Carreau-WLF equation [31] η(T,γ )   Fig. 1a. Notice that we do not include a back-coupling of the fiber orientation on the viscosity by setting the particle number N p to zero. However, we do include the temperature and shear rate dependence of the viscosity. For the simulation at hand, we use the ORW3 closure approximation [14]. The resulting fiber orientation is shown in Fig. 1 in a bottom and side view with the color scale of Fig. 2. An almost aligned state is dominant throughout the component, with more isotropic and planar states visible at the weld lines and the flat parts, respectively. A quantitative analysis of the fiber orientation, see Fig. 19, reveals that almost the full fiber orientation triangle is indeed covered by occurring fiber orientations. Only the vicinity of the isotropic state and a neighborhood of the line connecting the unidirectional and the planar isotropic state are left blank. However, most of the fiber orientations are concentrated around (λ 1 , λ 2 ) = (7.7, 0.15).
The calculated fiber orientation tensors serve as input for a component analysis using the commercial finite element software ABAQUS [1]. To assess the influence of the fiber orientation, for comparison, the computation was repeated for the same component, but where the fiber orientation was assumed isotropic throughout the whole geometry.
In Fig. 21 the strain in x-direction (cf. Fig. 20 for the axes) is shown, comparing the constant isotropic to the distributed fiber orientations. For the initial elastic response (t = 0), the strain distributions are even qualitatively similar. However, after 200 h strong differences in the strain fields can be observed. Indeed, for the isotropic fiber orientation spots with much higher strains are visible (in yellow and orange).
To further investigate this effect we show the corresponding von Mises equivalent stress in Fig. 22, and investigate the local relaxation behavior. For the constant isotropic component, the von Mises stress apparently changes only slightly during the loading history. For the distributed orientation, however, there is a load rebalancing comparing the stress fields at t = 0 h and t = 200 h.
To gain deeper insight, we investigated the von Mises stress history up to 200 h in the most heavily loaded element for the distributed orientation, cf. Fig. 23. At the beginning of the loading, the stress level lies slightly above 25.5 MPa. After 200 h, the stress level increased to about 28 MPa, i.e. by about 10%.
We see that the constant isotropic computing overestimates the occurring local strains and does not account for the strongly anisotropic relaxation behavior of the component under consideration. In particular, an anisotropic component design necessitates larger safety factors than the computations with anisotropic relaxation behavior.

Conclusion
This work concerned the multiscale simulation of short fiber reinforced thermoplastics. Even though the fibers are much smaller than the typical component scale, the mechanical properties on the component scale are not homogeneous, but depend on the local fiber orientation. Thus, working with a single surrogate model for a fixed fiber orientation is insufficient, and the need to interpolate surrogate models corresponding to different fiber orientations arises.
We have proposed a versatile and robust approach to interpolate surrogate models, rigorously derived from the twopotential-framework [23,25]. The advantage of the approach is the ability to interpolate models of completely different type, in particular in terms of internal variables.
To test the framework on a non-trivial example, we investigated the viscoelastic creep behavior of PA66, reinforced by short glass fibers, and used Moulinec-Suquet's FFT-based computational homogenization scheme and Schapery's collocation method to identify anisotropic viscoelastic surrogate models for different fiber orientation states. It is shown that even for a rather rough triangulation of the fiber orientation triangle the interpolation errors are small.
Finally, the full identified viscoelastic model was used as a material model for an ABAQUS simulation of a benchmark component, exhibiting an anisotropy-induced strain rearrangement on the component scale.
This work contributes to bridging the gap between the insights provided by multiscale modelling of fiber reinforced composites achieved in recent years and fast computational analyses carried out by the constructing engineer. Indeed, even though the filling of the database is computationally demanding, the resulting material card enables extremely fast computations and is universal in the sense that it can be (re-)used for different components built from the same material.
As already mentioned, the proposed fiber orientation interpolation scheme is not restricted to viscoelastic material behavior. The availability of surrogate models for different fiber orientations is the only prerequisite. In particular, viscoplastic or damage effects could be investigated. Furthermore, the interpolation scheme is not limited to problems at small strains, but extends effortlessly to the finite strain framework. As the interpolation is based on the free energies and the dissipation potentials, the question which strains and stresses should be interpolated does not arise but follows directly from the energetic formulation. Also, the scheme could be extended to incorporate local variations in fiber volume fraction.