Numerical Simulation of the Elastic–Ideal Plastic Material Behavior of Short Fiber-Reinforced Composites Including Its Spatial Distribution with an Experimental Validation

: For the numerical simulation of components made of short ﬁber-reinforced composites, the correct prediction of the deformation including the elastic and plastic behavior and its spatial distribution is essential. When using purely deterministic modeling approaches, the information of the probabilistic microstructure is not included in the simulation process. One possible approach for the integration of stochastic information is the use of random ﬁelds. In this study, numerical simulations of tensile test specimens were conducted utilizing a ﬁnite deformation elastic–ideal plastic material model. A selection of the material parameters covering the elastic and plastic domain are represented by cross-correlated second-order Gaussian random ﬁelds to incorporate the probabilistic nature of the material parameters. To validate the modeling approach, tensile tests until failure were carried out experimentally, which conﬁrmed the assumption of the spatially distributed material behavior in both the elastic and plastic domain. Since the correlation lengths of the random ﬁelds cannot be determined by pure analytic treatments, additionally numerical simulations were performed for different values of the correlation length. The numerical simulations endorsed the inﬂuence of the correlation length on the overall behavior. For a correlation length of 5mm, a good conformity with the experimental results was obtained. Therefore, it was concluded that the presented modeling approach was suitable to predict the elastic and plastic deformation of a set of tensile test specimens made of short ﬁber-reinforced composite sufﬁciently.


Introduction
The use of short fiber-reinforced composites (SFRCs) has increased significantly over the years. Due to the specific strength and stiffness in comparison to conventional polymers and the applicability for mold injection production processes, they are of special interest, e.g., in the automotive industry [1]. However, the mold injection process, as well as the finite length of the reinforcing elements result in a microstructure having a probabilistic nature, which significantly affects the mechanical properties [2] and, thus, the structural response. Following this, a probabilistic modeling approach appears to be promising to exploit the lightweight potential in the best possible way.
To reduce the computational costs, the main challenge here is the incorporation of microstructural information without explicit modeling. One stochastic technique to model spatial data is random fields [3][4][5], which are also used in the context of material modeling. Soize and Guilleminot developed a comprehensive framework starting with non-Gaussian positive-definite matrix-valued random fields [6] and tensor-valued random fields for a meso-scale stochastic model of anisotropic elastic microstructures [7]. Based on this paper by Soize, the approach was extended by Guilleminot et al., covering the stochastic fluctuations in fiber-reinforced composites on the mesoscale [8,9], random interphases from atomistic simulations of polymer nanocomposites [10], and a multiscale approach for heterogeneous materials with non-Gaussian random fields [11]. Beside this, random fields are widely used, e.g., for geosystems [12], thin-walled composite cylinders [13], three-dimensional concrete microstructures [14], and the representation of the continuous mode conversion observed by the propagation of guided ultrasonic waves in thin-walled structures made of fiber-reinforced composite [15]. A first application of random fields in the context of nonlinear behavior for isotropic material was provided by Zheng et al. [16]. The application to SFRCs has been limited to the linear elastic domain so far [2]. However, SFRCs show prominent nonlinear behavior, because the matrix material enters the plastic domain even at operation loads [17]. Hence, to predict the structural response correctly even at low stress levels, information about the plastic deformation must be considered. Since the plastic deformation is localized initially, the required nonlinear modeling approach needs to include information about the probabilistic nature of the microstructure causing spatially distributed material properties [18].
Beside the research on random fields for probabilistic material modeling, the modeling of SFRCs at the different scales mostly focuses on homogenization approaches and does not consider the spatial distribution of the material properties on the component level. The determination of a representative volume element (RVE) for random composites was discussed by Savvas et al. [19] by combining the extended finite element method with a Monte Carlo sampling. Breuer et al. [20,21] applied the RVE to SFRCs including artificial neural networks. Zhang et al. investigated the strain rate dependence of SFRCs based on RVEs [22], and Jia et al. applied the RVE concept to cyclic mechanical and thermal loading [23].
Hence, the main objective of this research was the incorporation of the probabilistic material characteristics of SFRCs into a modeling approach covering the elastic and plastic domain on the component level. The presented paper is linked to a correlation structure study of the elastic-ideal plastic material behavior [24], where the moving window method [19,25,26] was utilized for the elastic domain and a homogenization method was used to derive the local apparent plastic material parameters [27,28], respectively. Here, the obtained results are now employed for the numerical simulation of SFRCs on the component level, including an experimental validation. Therefore, first, cross-correlated random fields are generated to describe the spatial distribution and the probabilistic nature of the material properties introduced by the stochastic attributes of the microstructure. In a second step, a transversely isotropic elastic-ideal plastic material model including finite deformation is established and implemented in COMSOL Multiphysics ® to simulate tensile tests for specimens made of SFRCs. For validation of the presented modeling approach incorporating the spatial distribution and probabilistic characteristics of elastic and plastic material properties, uniaxial tensile tests are conducted experimentally.
Subsequently, the structure of the present paper is as follows. In Section 2, the experimental investigation is given. This includes a detailed description of the specimens and the experimental setup and procedure, respectively. This is followed by the generation of cross-correlated random fields in Section 3. After the implementation of an algorithm in Multiphysics ® covering elasto-plastic material behavior in Section 4, the random fields are used in Section 5 to incorporate spatially distributed material properties into the numerical modeling procedure. Finally, Section 6 gives a summary and conclusion.

Experiments
This section gives detailed information about the experimental investigation of the elastic and plastic properties of specimens made of SFRCs by tensile tests. First, the specimen specifications are given before presenting the experimental setup and procedure. The section is concluded with an overview of the experimental results covering the elastic and plastic domain.

Specimens
In this paper, the elastic and plastic deformation for specimens made of Ultradur B 4300 G6 [29], a polybutylene terephthalate (PBT) matrix material filled with glass fibers, were investigated. The fiber mass fraction was set to 30%, which is equal to a fiber volume fraction of 18.2%. The specimens were cut out of a larger plate that was manufactured by mold injection. The size of this plate was 300 mm × 300 mm and had a thickness of 3 mm. The geometry of the tensile test specimens was defined in accordance with DIN-ISO-527-1 [30]. However, due to the plate dimensions, the size of specimen type 1B was slightly adapted. Figure 1 gives an overview of the initial plate, the position of each specimen, and its exact dimensions.

Experimental Setup and Procedure
For the tensile tests covering the elastic and plastic deformation until failure, a Zwick-Roell Z050 tensile testing machine was used. With this testing machine, a load of up to 50 kN can be applied. As suggested by DIN-ISO-527-1 for plastics, the tensile test was performed with displacement driven at a speed of 1 mm min −1 . The strain was captured by a contactless laser extensometer with an initial measurement length of 50 mm. Furthermore, to analyze the spatial distribution of the material properties, additional measurement points within the original measurement length were added. They had a distance of 5 mm; see Figure 1.
For the elastic characterization of the material, Young's modulus was derived from the measured stress-strain data. It was determined by computing the slope of the curve between a strain level of 0.1% and 0.2%. For the plastic deformation, the strain at failure and the maximum stress level were obtained from the measurement data. This was performed based on the data for a measurement length of 50 mm for all specimens. In addition, to analyze the spatial distribution qualitatively, the data were obtained for sections of 15 mm for one specimen. Hence, three non-overlapping different sections of the initial measurement length were examined individually. In all cases, the engineering strain was measured. Figure 2 summarizes the results of the conducted tensile tests in the elastic and plastic domain. On the left-hand side, the stress-strain curve for each specimen based on the measurement length 50 mm is plotted. Furthermore, the diagram holds data provided by the manufacturer [29]. The stress-strain curves clearly revealed the prominent plastic deformation of the SRFC even at low stress levels. However, in comparison with the data sheet, the experimentally obtained stress and strain level at failure were significantly lower. This also held for the slope of the curve in the elastic domain. For a better overview, the corresponding values of Young's modulus, as well as the stress and strain level at failure are collected in Table 1. For each parameter, the mean value and the standard deviation are provided.

Results
The differences between the experimental data and the data provided by the manufacturer were most likely related to the manufacturing process of the specimen. It significantly influenced the microstructural characteristics such as the fiber length and fiber orientation of the specimen and the quality of the fiber-matrix bonding.
In the second diagram, the focus is set on the spatial distribution of the material properties. For one specimen, the stress-strain curve was not only determined for a measurement length of 50 mm, but also for three sections of 15 mm in length each. The sections did not overlap and were within the original measurement length. The different curves revealed that the plastic deformation was distributed over the specimen, confirming the assumption of spatially distributed material properties induced by the finite length of the reinforcing elements and the probabilistic characteristics of the microstructure.

Generation of Cross-Correlated Random Fields
In this section, first, random fields and their generation by numerical methods are discussed. Afterwards, the calculation of the cross-correlated random fields and their application to the problem at hand are addressed.

Random Fields
Random fields Z(ω, x) represent spatially distributed random variables Z(ω). To describe a random field completely, the function of the mean value µ(x), variance σ 2 (x), and correlation coefficient ρ(x, x ) are required.
Utilizing random fields in the context of material modeling, it is important to note that the material parameters are positive in nature. This does not necessarily apply to the coefficients of strain-energy density functions. However, it is indispensable for the yield strength, which contradicts the found normal distribution of the yield strength [24]. Therefore, when using Gaussian random fields and, hence, assuming normally distributed underlying random variables, negative realizations are possible [31,32]. In conclusion, the use of Gaussian random fields for the representation of the material properties in the context of multi-scale modeling of heterogeneous material is controversial [6]. With non-Gaussian random fields, negative realizations can be avoided, and hence, a stochastic solution of second-order and the positive nature of the elasticity coefficient are guaranteed [10,11,[33][34][35][36].
Non-Gaussian random fields M(ω, x) are defined by where G is a nonlinear mapping operator and Z(ω, x) a centered, homogeneous Gaussian random field [11]. With respect to the following investigations, where the concept of random fields was used within the finite element method framework, only discretized random fields were used, requiring a finite sampling. Consequently, as long as negative values are excluded, the difference between the results based on non-Gaussian fields and Gaussian fields may not differ significantly. Therefore, the approach presented here was based on homogeneous Gaussian random fields, which also allowed using the well-known techniques presented below.

Generation of Random Fields by Numerical Methods
Initially, random fields are continuous functions with respect to the spatial coordinates. However, to be able to utilize this concept within the framework of material modeling and numerical simulations by the finite element method (FEM), a discrete representation is necessary. Over the years, many different methods have been developed for this purpose. Examples are the midpoint method [37,38], spatial averaging [39], and the shape function method [40][41][42]. Commonly used approaches are based on the series expansion technique [43]. Spanos et al. [44] extended the approach by expressing random fields by a spatial decomposition of the correlation functions derived by the well-known Karhunen-Lòeve expansion (KLE) [45]. The spatial decomposition obtained by the KLE for a mean free random field with a standard deviation equal to one was given by Cho et al. [46]: Here, Z n (ω) are uncorrelated, standard normally distributed random variables and λ n and φ n (x) are the eigenvalues and eigenfunctions of the correlation function kernel, which can be obtained by solving a Fredholm integral equation of the second kind. Since the integral equation: has closed solutions only for a few types of correlation functions defined on a rectangular domain [47], e.g., the exponential correlation function [31,48], numerical methods are often required. This is especially the case for multidimensional fields [49]. By deploying numerical integration methods, the eigenfunctions φ n (x) as a solution of the Fredholm integral are approximated by a set of functions h i : where the parameters d n i are unknown and need to be determined. One technique here is the expansion optimal linear estimate (EOLE) developed by Li et al. [49]. It is based on the linear estimation theory and belongs to the group of series expansion methods. Furthermore, it can be shown that the EOLE is equivalent to the Nyström method [50] with uniformly distributed integration points [47]. Within the Nyström method, the integral eigenvalue problem is written as which can be reorganized in matrix notation as The N × N matrix C is symmetric positive semi-definite and holds the elements c n,i = Cov[x n , x i ]; the integration weights w i are stored in the diagonal matrix W of the size N × N, and y j is a column matrix with the entries y j,n =φ j (x n ). Assuming a uniform distribution of the points x i over the domain Ω and a equispaced structured grid, respectively, all integration weights w i are the same. In this case, the matrix W reads with I being the identity matrix and w = |Ω|/N. This leads to giving the equivalent of Equation (6) for the EOLE, with the eigenvaluesλ * j and eigenfunctions y j of the covariance matrix C. The eigenvalues of the EOLE are related to the eigenvalues of the Nyström method byλ * Finally, for normalized eigenvectors, the truncated KLE approximated by the EOLE reads

Cross-Correlated Random Fields
In the context of nonlinear material modeling, the correlation complex in nature, because the individual parameters are approximated best by different correlation functions and the parameters are strongly cross-correlated [24,51]. Therefore, to achieve a proper representation of the material properties, the generated random fields must include the information about the cross-correlation. Cho et al. introduced two different algorithms for the expansion of multi-correlated processes, namely the multiple uncorrelated Karhunen-Loève expansion (muKL) and the multiple correlated Karhunen-Loève expansion (mcKL) [46]. The main difference between these two procedures lies in the generation of the used random variables. In the context of the muKL, all correlated processes are described by a single set of uncorrelated random variables [46]. In contrast to this, the mcKL provides multiple sets of mutually correlated random variables for the discretization of correlated processes. Since the modeling of nonlinear material behavior requires a set of random variables for each parameter, the mcKL is applicable in this context.
Within this algorithm, the same equations as for uncorrelated random fields are used and solved. However, in an additional step, the uncorrelated random variables and eigenfunctions are transformed to obtain correlated random variables and eigenfunctions, respectively, which meet the given correlation structure. Here, first, the matrix K is introduced. The matrix is defined as where Z i k and Z j m are a set of zero-mean uncorrelated random variables for processes i and j, respectively. Using Equation (2), the cross-covariance between the processes i and j is written as [52] The result is an integral equation similar to the one for the auto-correlation shown in Section 3.1.2. Following Equation (3), the elements of K are obtained by solving In a final step, the matrix K is used to determine sets of correlated random variables for each process. It is assumed that K is positive definite, which allows a Cholesky decomposition given by [52] K = RR T .
Applying R on a set of random variables leads tõ where This means thatZ contains only uncorrelated random variables. In reverse, applying R on a set of uncorrelated random variablesZ gives a set of correlated random variables that represent cross-correlated random fields. The same procedure is also applicable to obtain cross-correlated eigenfunctions.

Application to the Elastic-Ideal Plastic Material Behavior of SFRCs
In this section, the presented framework for the generation of cross-correlated random fields is used to represent the material properties of the elastic-ideal plastic material behavior of SFRCs. The information of the correlation structure was derived and discussed in detail by Rauter et al. [24]. It is shown that the correlation structure of the material parameters is approximated best by a combination of the triangle correlation function, defined as and an exponential correlation function [32,48,53]: The parameter a indicates the strength of the correlation between two parameters; ξ 1 and ξ 2 hold the spatial coordinates, and b 1 and b 2 are the corresponding correlation lengths, respectively. In the case of auto-correlation, a is equal to one; in the case of cross-correlation, a can take values between −1 and 1.
The values of a were set in accordance with the results obtained by Rauter et al. [24]. In this initial study, the numerical simulation of the elastic-ideal plastic material behavior of SFRCs, only two parameters were described by spatially distributed material properties to reduce the complexity of the nonlinear numerical simulations. To ensure that both the elastic and plastic domain are captured by the stochastic multiscale approach, one parameter of the strain-energy density function, which is related to the elastic domain, and the yield strength are represented by random fields. Based on this, the routine can be reduced to two parameters capturing the cross-correlation of Λ and σ y . The matrix holding the factor a for each material parameter combination reads a = 1.00 0.53 0.53 1.00 . (19) In addition to the parameter a, the required mean values, standard deviations, and correlation lengths are provided in Table 2. Since the correlation length cannot be determined solely on numerical data, the values were normalized to the correlation length of the parameter Λ. For more details about the material model and the provided data, the reader is kindly referred to Section 4.1 and the correlation analysis presented by the authors [24], respectively. With this information at hand, realizations of random fields for these two material parameters were generated using the numerical methods for two-dimensional cross-correlated homogeneous second-order random fields presented in Section 3.1. In accordance with the numerical simulation of the tensile test in Section 4 and the experimental procedure in Section 2, the domain has a size of 50 mm × 3 mm. One realization of the parameter set is presented in Figure 3. The correlation length of Λ was set to l x = 20 mm and l y = 3 mm in the horizontal and vertical direction, respectively.
In addition, the quality of the procedure was analyzed. Therefore, in total, 1000 realizations for each parameter were generated. Based on the obtained data, the correlation structure was derived and compared to the results presented in the former study by the authors [24]. Figure 4 holds the results for Λ and σ y . The first two diagrams give the analytic auto-correlation functions of Λ and σ y , based on the results obtained in a previous correlation analysis [24] and the correlation function based on the 1000 realizations, respectively. In the last diagram, the cross-correlation function is plotted. For a better overview, the correlation function plots are limited to the horizontal direction. In addition to this, Figure 5 gives histograms of these two parameters based on all 1000 realizations. As discussed before, in this investigation, the material parameters are represented by Gaussian random fields. Since Gaussian random fields allow negative values in general, the distribution of the sampled yield strength was analyzed more in detail. Based on the 1000 realizations and a discretization truncated after 2000 terms, in total, 2 × 10 6 values were sampled for the yield strength. The corresponding mean value and standard deviation were 125.6 MPa and 14.7 MPa, which meet the values provided in Table 2. Furthermore, the minimal and maximal sampled values were 194.7 MPa and 64.0 MPa, respectively. Since the probability density function of the yield strength was not modified to avoid negative values, the moderate standard deviation of the Gaussian distribution resulted only in positive realizations of the random variable. Therefore, the use of second-order Gaussian random fields appears to be appropriate to incorporate spatially distributed material properties into the numerical simulation procedure.
Subsequently, the results for the correlation functions, as well as the material parameter distributions showed a very good agreement with the analytic values. Hence, the selected procedure is suitable to represent the locally varying material properties of the elastic-ideal plastic material behavior of SFRCs by random fields.

Numerical Simulation
In this section, the spatially distributed material properties described by cross-correlated second-order Gaussian random fields are combined with an algorithm for the elastic-ideal plastic material behavior. The main objective was the integration of the locally varying material properties, which were observed during the experimental investigations in Section 2, into a numerical modeling procedure. Therefore, first, the framework and constitutive equations for the transversely isotropic elastic-ideal plastic material behavior and its implementation in COMSOL Multiphysics ® are presented. The developed algorithm was then validated in detail. which need to satisfy the Kuhn-Tucker complementary conditions and the consistency condition [54]. The plastic constitutive material model introduced here was based on the research provided by Hashiguchi et al. [55] and Eidel et al. [56]. It was adapted to the transversely isotropic material by using the strain-energy density function by Bonet et al. [57]. Furthermore, hardening was not considered. Hence, the material was assumed to show the ideal plastic behavior [24].
Since the modeling approach covers finite deformation, the constitutive model needs to be formulated with respect to the nonlinear continuum mechanical framework. Therefore, the deformation is described by the multiplicative decomposition of the deformation gradient, also known as Kröner-Lee decomposition. The decomposition is given by [58][59][60][61] Here, F el holds the elastic and F pl the plastic deformation, respectively. Since the deformation gradient was used to map continuum mechanical quantities from the reference to the current configuration and vice versa, this decomposition revealed the definition of an additional configuration, the so-called intermediate configuration [54,58,60]. As shown in Figure 6, with F pl , the quantities were mapped from the reference to the intermediate configuration and with F el from the intermediate to the current configuration. Accordingly, holds, with C = F T · F being the right Cauchy-Green tensor, defined in the reference configuration. It also splits into an elastic and a plastic part written by the decomposition [62]: where Assuming J 2 -plasticity and, hence, plastic incompressibility leads to and respectively. With respect to the introduced configurations,C el , which is located in the intermediate configuration, contains the elastic deformation, whereas C pl was assigned to the reference configuration and includes the plastic deformation.
holds, with being the right Cauchy-Green tensor, defined in the reference configuration, which also splits into an elastic and a plastic part by Assuming 2 -plasticity and hence, plastic incompressibility, leads to and respectively. With respect to the introduced configurations pl contains the plastic deformation in the reference configuration and̄ el gives the elastic deformation in the intermediate configuration, which is indicated by the bar.
A crucial part of the plasticity framework is the plastic flow rule providing information about the irreversible deformation of a body. The fundamental equation is the second law of thermodynamic provided by the Clausian-Duhem-inequality for isothermal processes, which is written as 28 Figure 6. Configuration definition within the nonlinear continuum mechanics framework.
A crucial part of the plasticity framework is the plastic flow rule providing information about the irreversible deformation of a body. The fundamental equation is the second law of thermodynamics provided by the Clausian-Duhem inequality for isothermal processes, which is written as [63] −Ψ + S : where the dot indicates the derivative with respect to time, S is the second Piola-Kirchhoff stress tensor, and Ψ the Helmholtz free energy. Introducing the Helmholtz free energy as a function of C, the associated rate can be expressed bẏ Since, in this case, both kinematic and isotropic hardening were omitted, only the part Ψ el of the Helmholtz free energy Ψ, which was assigned to the elastic deformation, contributes to Equation (26). Following the multiplicative decomposition of F, the derivative of the right Cauchy-Green tensor with respect to time is obtained by [56] with the plastic velocity gradient in the intermediate configuration Combining the inequality in Equation (26) and the time derivative of C, which must hold for arbitrary thermodynamic processes, one can derive, first, the constitutive equation for elastic deformation, which defines the second Piola-Kirchhoff stress tensorS in the intermediate configuration: and, second, the reduced form of the Clausian-Duhem inequality [56,63]: again omitting isotropic and kinematic hardening. In this formulation, the Mandel stress tensorM is introduced. It is defined asM =C el ·S (32) and stands in the intermediate configuration.
The inequality in Equation (26) is satisfied by the associated flow rule [56,63]: assuming an isotropic yield function φ [56]. In this research, the framework of J 2 -plasticity was used. Therefore, the yield function is introduced as a von Mises criterion by [55] where σ yield is the reference yield stress for uniaxial loading and, hence, an independent material parameter and ||M dev || denotes the norm of the deviatoric part of the Mandel stress. Substituting the von Mises yield function in Equation (34) into Equation (33) leads to [55,63] L pl =λM dev ||M dev || .
With this framework at hand, the flow rule is written as [56] To solve this equation, an implicit backward Euler-type tensor exponential-based time integration scheme was used. Following this procedure, the updated plastic part of the deformation gradient is obtained by [56,64] F pl,n+1 = Q pl,n+1 F pl,n with the exponential function: The subscripts n and n + 1 denote the corresponding steps, where n indicates the initial and n + 1 the updated values.
The presented framework of elastic-ideal plastic deformation was applied to the transversely isotropic material behavior characterizing SFRCs in the next step. Therefore, the constitutive equations were used within the context of the finite element method. The resulting routine can be divided into four steps, which are presented in detail below.

Solution Procedure Elastic Predictor
In the initial step, the deformation of the body was assumed to be purely elastic, which is described in the intermediate configuration, as shown before. Following the research of Simo et al., who derived the framework of the finite strain elasto-plasticity by combining the multiplicative decomposition of the deformation gradient with hyperelastic strain-energy density functions [65][66][67], the elastic deformation is obtained by evaluating the hyperelastic strain-energy density function Ψ el with respect to the elastic part of the deformation.
In this study, the elastic material behavior of the SFRC is described by the potential [57] Ψ el (C el ) = Here, Λ, µ, α, β, and γ are the independent material parameters characterizing the transversely isotropic material behavior. In addition to the well-known invariants I 1 , I 2 , and I 3 for the description of the isotropic material behavior [68,69], the quantities I 4 and I 5 give the pseudo-invariants, which are defined for a symmetric second-order tensor B by [70,71] I 4 = a · B · a and The vector a gives the fiber orientation and, therefore, holds the transversely isotropic characteristics of the material. In this case, it was assumed that the fibers are mainly orientated along the tensile test specimen. Hence, a 1 = 1 and a 2 = a 3 = 0 hold.
With respect to the used transversely isotropic potential and the tensor derivatives, given by the second Piola-Kirchhoff stress tensor in the intermediate configuration reads

Checking Yield Condition
To be able to evaluate the yield condition, provided in Equation (34), in the intermediate configuration, the deviatoric part of the Mandel stress must be computed to derive the corresponding von Mises stress.
If the resulting von Mises stress is lower than the yield strength of the material, the body shows only elastic deformation. Therefore, no further action is required. The stress measure can be updated with the trial value obtained by the elastic predictor step. If, however, the von Mises stress is larger than the yield strength, the body shows not only elastic, but also plastic deformation. In this case, the determination of the resulting plastic deformation is covered by the plastic corrector step.

Plastic Corrector Step (Return Mapping)
Since the yield condition is given in the intermediate configuration, this also holds for the computation of the updated plastic part of the deformation gradient.
As shown before, an associative plastic flow rule introduced by Eidel et al. [56] was used to describe the evaluation of the plastic deformation. With this framework at hand, the return mapping procedure in the case of plastic deformation leads to the following system of equations (cf. Equation (37)): which needs to be solved to obtain the updated quantities at the step n + 1 for the plastic deformation. Since the plastic part of the deformation gradient depends solely on the plastic multiplier ∆λ (see Equations (37) and (38)), a Newton-scheme algorithm was used to determine the updated components of the plastic deformation gradient. Therefore, the yield criterion was iteratively solved with respect to ∆λ, and the procedure was stopped when a maximum number of iterations or a converged solution was reached by satisfying the condition: Furthermore, the required derivative was calculated numerically by deploying a one-sided difference scheme:

Updating
With the obtained value of ∆λ satisfying the yield condition, the updated plastic deformation gradient F pl,n+1 can be calculated by evaluating Equations (37) and (38). Next, the resulting second Piola-Kirchhoff stress tensor in the intermediate configuration was derived. Therefore, the updated elastic part of the deformation gradient was computed and introduced in Equation (42). COMSOL Multiphysics ® requires the computation of the second Piola-Kirchhoff stress tensor and the consistent tangent modulus tensor with respect to the reference configuration. Hence, a pull-back operation: was applied to obtain the second Piola-Kirchhoff stress tensor in the reference configuration. The last step was the computation of the consistent tangent modulus tensor. The corresponding procedure is provided in detail in the next section.

Variational Formulation and Consistent Tangent Modulus Tensor
For a body B, the local form of the balance of momentum with respect to the reference configuration reads Div P + ρ 0 b = ρ 0ü .
Here, P is the first Piola-Kirchhoff stress tensor, ρ 0 is the density, b are the volume forces, and u is the displacement field. To solve this differential equation, Neumann and Dirichlet boundary conditions were used. They are given by where ∂B is the boundary surface of B and N gives the normal to the boundary surface ∂B σ . To ensure that the boundary conditions cover the complete surface of the body, ∂B = ∂B ∪ ∂B σ and ∂B ∩ ∂B σ = ∅ hold. Applying standard variational calculus on Equation (47) and P = F · S leads to for the equilibrium state of the body B. Based on the definition: the virtual right Cauchy-Green deformation tensor, expressed with respect to the virtual deformation gradient δF = Grad δu, is given by To solve this nonlinear equation within the FEM framework, a standard Newton-Raphson iteration scheme was used. To ensure a quadratic convergence rate, a consistent linearization of the principal of virtual displacement in Equation (50) is crucial. Therefore, in this paper, the perturbation technique developed by Miehe [72] was used. This approach is widely used and has been extensively investigated with respect to convergence properties and the application within the finite strain framework [56,73,74]. Due to the symmetry of S, the linear increment of F can be formulated as where ∆δC = (∆F T · δF + δF T · ∆F) is the linearized virtual right Cauchy-Green deformation tensor given in terms of the incremental deformation gradient ∆F = Grad ∆u. To obtain the incremental second Piola-Kirchhoff stress tensor ∆S, the consistent tangent modulus tensor is required. In the reference configuration, it is expressed by which can be written as The incremental change of the right Cauchy-Green tensor was obtained by applying a perturbation technique to the deformation gradient. Following this, the incremental change of the right Cauchy-Green tensor is given by with Subsequently, the perturbed deformation gradient reads With this, the corresponding perturbed incremental right Cauchy-Green tensor can be calculated, which enables one to compute the stresses by applying the presented return mapping algorithm. This leads to the stress increment: and, hence, the consistent tangent modulus tensor by evaluating Equation (55) 4.2. Implementation in COMSOL Multiphysics ® 4.2.1. Algorithm Using an external material option in COMSOL Multiphysics ® requires the generation of a C-based dll, which can be loaded via the COMSOL Multiphysics ® interface. With respect to the presented elastic-ideal plastic material model introduced in Section 4.1, the general stress deformation routine was used. The input and output parameter for a general stress deformation subroutine for the implementation in COMSOL Multiphysics ® are provided in Table 3. In addition to these variables, optional ones can be stored utilizing state variables. Table 4 gives an overview of the state variables used in the presented algorithm for the elastic-ideal plastic material behavior of SFRCs. Output: Jacobian of the stress with respect to the deformation gradient as a 6-by-9 matrix of partial derivatives of the components of sPK with respect to the components of F With these definitions and the framework provided in Section 4.1, an algorithm was established to compute the second Piola-Kirchhoff stress tensor and the Jacobian matrix of the second Piola-Kirchhoff stress tensor with respect to deformation gradient. The structure of the developed routine is depicted in Figure 7. The main routine EVAL reads the input variables, calls the subroutine RETMAP, which is responsible for the return mapping procedure, and updates all output and state variables. A detailed overview of the routine is provided by Algorithm 1. In addition, Algorithm 2 and Table 5 hold all the information about the subroutine RETMAP. The details of all remaining supporting routines indicated in Figure 7 are given in Appendix A.
Algorithm 1 Main routine implemented in COMSOL Multiphysics ® for the calculation of S and J. Input: see Tables 3 and 4 Output: see Tables 3 and 4 1: procedure EVAL 2: read material parameters 3: read state variables 4: initialize Jacobian and identity matrix 5 compute Lagrangian consistent tangent modulus tensor Equation (55) 13: reorganize C in Voigt notation 14: compute Jacobian by calling COMSOL Multiphysics ® utilize function csext_jac_con 15: return 16: end procedure Algorithm 2 Return mapping algorithm implemented in COMSOL Multiphysics ® for the calculation ofS and F pl .

Validation of the Algorithm
For the validation process, a two-dimensional numerical model under the plane strain assumption representing the cross-section of a tensile test specimen was used. The size of the numerical model was 50 mm × 3 mm. In Figure 8, the numerical model and the corresponding boundary conditions are depicted. To reduce convergence issues during the plastic deformation, the simulation was displacement driven. Besides that, the numerical model was discretized by four-node elements with a size of 0.25 mm × 0.25 mm. Hence, over the thickness, the cross-section was discretized with 12 and along the cross-section with 200 elements. 16 RAUTER AND REINA 50 mm 3 mm  Tables 2 and 3) 2: read material parameters 3: read state variables 4: initialize Jacobian and identity matrix 5: computē +1 in intermediate configuration and pl, +1 calling Algorithm 3 ⊳ Eq. (40) 6: compute +1 by pull back of̄ +1 7: save +1 in Voigt notation 8: update state variable , +1 9: compute pertubed deformation gradient 10: computē +1 in intermediate configuration 11: compute , +1 by pull back of̄ +1

12:
compute Lagrangian algorithmic tangent moduli 13: transform ℂ to in Voigt notation 14: compute Jacobian by calling COMSOL Multiphysics utilize function csext_jac_con   Due to the fact that COMSOL Multiphysics ® experiences some issues when combining a user-defined hyperelastic strain-energy density function with plasticity, the provided algorithm was validated in two steps. In the first step, the external material algorithm was validated based on the well-known hyperelastic neo-Hookean material model [75]. Therefore, the strain-energy density function within the presented algorithm was adapted to the one implemented in COMSOL Multiphysics ® for the neo-Hookean material model given by Subsequently, two numerical simulations were performed, one with the implemented neo-Hookean hyperelastic material model and the second one with the provided external material algorithm. In both cases, the ideal plasticity and a plain-strain state were assumed. The corresponding material properties are provided in Table 6. Furthermore, a total displacement of u x = 1 mm corresponding to a strain of 2% was applied to ensure plastic deformation. In a second step, the strain-energy density function provided by Bonet et al. [57] was passed to COMSOL Multiphysics ® as a user-defined potential. Since this cannot be combined with the ideal plastic behavior due to convergence issues, this was only simulated assuming elastic deformation. The results were compared again to the results provided by Algorithm 1. The corresponding material parameters are also provided in Table 6. The results are depicted in Figure 9. The first row gives the result for the isotropic neo-Hookean simulation, and the second row holds the results for the transversely isotropic simulations. The diagrams on the left-hand side contain the stress-strain curves. The von Mises stress is plotted in black, and the horizontal component of the second Piola-Kirchhoff stress is plotted in gray, respectively. Following the theoretical framework provided in Section 4.1, it is important to note that the von Mises stress is given in the intermediate configuration, whereas the second Piola-Kirchhoff stress tensor stands in the reference configuration. The solid lines give the results of the implemented material model in COMSOL Multiphysics ® , and the small circles belong to the results obtained by the external material routine. As indicated by the diagrams on the right-hand side, which provide the deviation of the external material routine to the implemented material model in percent, both validation steps showed good agreement. Both comparisons experienced an increasing deviation at large displacements, which was still very small.

Application to Tensile Test Specimen
In this last step, the presented modeling approach was used to incorporate the spatially distributed elastic-ideal plastic material behavior of the SFRC into a numerical simulation of tensile tests. As mentioned before, the probabilistic characteristics of the material parameters were analyzed in a former study of the authors [24]. The main results are summarized in Table 2 for a window size of 750 µm, which showed the best approximation for the simulation on the component level [76]. In this section, the numerical model is briefly introduced, the numerical results are shown and the corresponding results discussed in detail, including a comparison with the experimental data presented in Section 2. This includes the analysis of the influence of the correlation length on the standard deviation of the obtained results and the selection of an optimal correlation length for which the numerical results and the experimental measurements coincided the best.

Numerical Model
The numerical model used to validate the novel modeling approach on the component level based on tensile test specimens was identical to the one used in Section 4.2.2. As already discussed in Section 3.2, only one parameter of the strain-energy density function and the yield strength were represented by the cross-correlated second-order Gaussian random fields. This reduced the complexity at this state of the investigation. Since, in this case, the standard deviation of the material parameter Λ did not cover the total standard deviation of elasticity tensor element C 11 , which is significantly related to the structural behavior under uniaxial tension, the corresponding value was set to 1.74 GPa. Furthermore, the experimental data clearly revealed that local plastic deformation occurred even at low stress levels. This was related to the yield strength of the matrix material and to the plastic deformation induced by the stress peaks at the edges of the reinforcing elements. In addition, the maximal stress level of the data sheet and, hence, the results of the homogenization were not reached. Therefore, the mean value and the standard deviation were adjusted and set to 100 MPa and 31 MPa, respectively. The remaining parameters were assumed to be homogeneous and were set to the mean values obtained by a previous correlation analysis [24]. Furthermore, a plain strain state was employed. This was eligible because the microstructure of the SFRC consisted of layers with different main fiber orientations [77,78]. The core layer, which was located in the center of the cross-section, was characterized by a fiber orientation perpendicular to the main flow direction of the injection mold process. With respect to the specimens used in this study, the fibers within the small core layer were perpendicular to the load direction and, hence, inhibited the out-of-plane deformation.
Further crucial parameters for the numerical simulation were the correlation lengths of the random fields. Since the correlation length significantly influences the local distribution of the material properties, it was strongly linked to the standard deviation of the numerical results. Therefore, in this study, different values were used to analyze the influence more in detail. In Figure 10, realizations for the strain-energy density function parameter Λ are depicted for a correlation length in the horizontal direction of l x = 5 mm, l x = 10 mm, and l x = 20 mm. In the vertical direction, the correlation length was set to l y = 3 mm. The effect of the correlation length on the random field is clearly observable. In analogy to the description of the oscillations, the correlation length can be interpreted as the wavelength of the dominant eigenfunctions. In this paper, the numerical simulations were performed for horizontal correlation lengths of l x = 5 mm, l x = 10 mm, and l x = 15 mm. The correlation length in the vertical direction was set to 0.3 mm for all simulations. This value was chosen to allow plastic deformation at different locations, as was observed during the experimental investigations. In combination with the ideal plastic modeling approach, the resulting stress-strain relation approximated the experimental data for the reinforced material. Were the correlation length be identical to the specimen's thickness of 3 mm, the plastic deformation was limited to the region where the yield strength was reached first, because the surrounding material cannot take higher stress levels. The obtained results would then be characterized purely by the ideal plastic behavior of the matrix material.
The numerical simulations were carried out for eight different realizations per correlation length. Hence, in total, 24 simulations were conducted. The results are presented below.

Results
The results of the numerical simulations are summarized in Figure 11. Furthermore, a comparison of the numerical data with the experimental results is included. The diagrams on the left-hand side contain the stress-strain curve (black), taken from the data sheet [29], the experimental data (gray) (see Section 2), and the numerically obtained results (colored) for correlation lengths of l x = 5 mm, l x = 10 mm, and l x = 15 mm from top to bottom. The diagrams on the right hold again the stress-strain curve (black), taken from the data sheet [29] and the numerical data (gray). One dataset is highlighted. The corresponding stress-strain curve was based on a measurement length of 50 mm. For these specimens, the spatial distribution of the material properties was qualitatively investigated by deriving the results for three sections of 15 mm each.  Figure 11. Results of the numerical simulation in comparison with the experimentally obtained results and a stress-strain curve provided by the manufacturer [29].
For a better comparability, Table 7 holds some additional information. For strain levels of 0.2 %, 0.5 %, 1.0 %, 1.5 %, and 2.0 %, the mean stress, as well as the standard deviations were computed for the numerical and experimental data. For all numerical results, the stress refers to the horizontal component of the second Piola-Kirchhoff stress tensor.

Discussion
The diagrams provided in Figure 11 reveal that the numerical obtained stress-strain curves clearly showed the expected elasto-plastic behavior. Despite the fact that the material behavior was modeled based on the assumption of an elastic-ideal plastic material model, the resulting stress-strain curves showed an overall good agreement with the experimental data. This was related to the localized plastic deformation and the rearrangement of the applied loading. In addition, the diagrams disclose, which the numerical results showed, a better agreement with the experimental data for a decreasing correlation length. This is also indicated by the values provided in Table 7. For a decreasing correlation length, the mean stress state converged to the experimentally obtained values; for a correlation length of 5 mm, the results showed an overall good agreement. The deviation between the mean stress states at the different strain levels for the experimental and numerical results was covered by the standard deviation. However, despite the good overall conformity, the results indicated two shortcomings. First, the experimentally obtained stress range at strain levels up to 1 % was not covered by the numerical simulation. This might be due to the reduction of the parameters, which were represented by random fields. Second, even though the composite material did not show explicitly ideal plasticity, the plastic deformation within the numerical simulations started at higher strain levels, as indicated by the experimental results. This might be resolved by reducing the yield strength locally or taking into account hardening effects.
The modeling approach was further endorsed by the results on the right-hand side by two important aspects. First, the modeling approach led to a plastic deformation covering the whole specimen and, hence, did not show localized plastic deformation. Second, when dividing the measurement length into small sections, each section showed a slightly different behavior. Both phenomena were confirmed by the experimental investigation.
All in all, it can be stated that the presented modeling approach is a very appropriate method to incorporate the spatial distribution of the elastic and plastic material properties of SFRCs.

Summary and Conclusions
In this manuscript, a probabilistic modeling approach for the elasto-plastic material behavior at finite strain based on second-order Gaussian random fields was introduced. The description of the material behavior was based on the Kröner-Lee decomposition. Following this, the strain-energy density function introduced by Bonet et al. [57] was combined with a classical J 2 -plasticity formulation assuming ideal plastic behavior. Since this material model was not available in commercially available finite element codes, first, an algorithm covering the transversely isotropic elastic-ideal plastic material behavior was developed and implemented in COMSOL Multiphysics ® . The validation based on the well-known neo-Hookean material model in the elastic and plastic domain, as well as utilizing a user-defined strain-energy density function in the elastic domain showed an excellent conformity.
Therefore, in the second step, based on the correlation analysis performed in a previous study [24], cross-correlated second-order Gaussian random fields were generated to represent the spatial distribution, the strain-energy density function parameter Λ, and the yield strength induced by the probabilistic characteristics of the microstructure. The remaining parameters were assumed to be homogeneous at this state of the study. With the spatially distributed material properties, tensile tests for specimens made of PBT-GF-30 were simulated. Even though the initial modeling approach was based on the elastic-ideal plastic material behavior, the spatial distribution of the two material parameters led to localized plastic deformation covering the whole measurement length. This caused a strain-dependent reduction of the material stiffness.
To validate the numerical results and the presented modeling approach, uniaxial tensile tests were conducted to characterize the elastic and plastic deformation of PBT-GF-30. The obtained stress-strain curves were compared to the numerical results. It was shown that, with an decreasing correlation length, the numerical data matched the experimental results quite well. The deviation at large strain can be assigned to the assumption of ideal plasticity. Furthermore, this initial investigation did not show a significant deviation between the results of different realizations at a low strain level. The reasons were the reduction of the complexity by representing only one parameter of the strain-energy density function by a random field and the limitation of Gaussian random fields. To ensure only positive values for the yield strength, low values were underestimated.
Summarizing, it can be concluded that the presented initial probabilistic approach is suitable for the numerical simulation of the elasto-plastic material behavior of SFRCs at finite strain including the spatial distribution on the component level without the need for an explicit microstructural modeling. However, the approach is not limited to SFRCs. It can be used to represent any material system with a probabilistic microstructure, which is the key advantage of the presented approach.
Future research will focus on the implementation of specimen unloading, hardening, non-Gaussian random fields, and the representation of additional parameters by stochastic methods.

Funding:
The financial support of the German Academic Exchange Service (Grant No. 91789376) is gratefully acknowledged.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available upon request from the corresponding author.
Acknowledgments: I want to acknowledge and extend my sincere thanks to Celia Reina for the valuable ideas and tips during our discussions, which significantly improved and shaped the research.

Conflicts of Interest:
The author declares no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of the data; in the writing of the manuscript; nor in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript:

Appendix A
For the implementation of the elastic-ideal plastic routine in COMSOL Multiphysics ® , several subroutines were used to calculate the second Piola-Kirchhoff stress tensor in the intermediate and reference configuration, respectively. Furthermore, the subroutine to calculate the exponential functions matexp provided by Hashiguchi et al. [55] written in Fortran was transferred into C. Besides these subroutines, COMSOL Multiphysics ® provides a library of utility functions for various tensor operations, computing principle and equivalent stresses. These routines are not explained here; a detailed overview is given in the documentation [79].