A strain-gradient formulation for fiber reinforced polymers: Hybrid phase-field model for porous-ductile fracture

A novel numerical approach to analyze the mechanical behavior within composite materials including the inelastic regime up to final failure is presented. Therefore, a second-gradient theory is combined with phase-field methods to fracture. In particular, we assume that the polymeric matrix material undergoes ductile fracture, whereas continuously embedded fibers undergo brittle fracture as it is typical e.g. for roving glass reinforced thermoplastics. A hybrid phase-field approach is developed and applied along with a modified Gurson-Tvergaard-Needelman GTN-type plasticity model accounting for a temperature-dependent growth of voids on microscale. The mechanical response of the arising microstructure of the woven fabric gives rise to additional higher-order terms, representing homogenized bending contributions of the fibers. Eventually, a series of tests is conducted for this physically comprehensive multifield formulation to investigate different kinds and sequences of failure within long fiber reinforced polymers.


Introduction
In the past decade, lightweight designs of composite materials have gained increasing attention in research including the field of computational engineering. This is primarily due to the wide range of industrial applications of composites. In the last years, novel machines for additive manufacturing have been designed and introduced to the market, which are able to construct tailor made composites with controllable mechanical properties by adding fibers during the manufacturing process into the produced parts. New challenges arise for the implementation of this class of materials in a suitable simulation environment to obtain realistic predictions of the material behavior. In particular, numerical investigations of complex structures made of composite materials, which undergo large deformations with regard to thermomechanical effects are demanding and of high interest at the same time. In the present work, we focus on long fiber reinforced thermoplastics which are complex materials with pronounced deformations and temperature dependent material properties. In addition, various damage and fracture mechanisms have to be investigated to formulate realistic simulation models.
Fiber reinforced materials can be understood as materials with a detailed microstructure since the length scales of the embedded fibers are small compared to the surrounding continuum. A most general framework to incorporate the effects of microstructures within a continuum formulation has been presented in the pioneering and most fundamental work on higher-gradient theories of Mindlin [52,53], see also the work of Germain [28], Toupin [75,74] as well as Eringen [27]. A general non-linear framework on higher-order theories has been proposed by Javili et al. [38]. Based on these most general formulations, specific mechanical problems like elastic nets and woven fabrics have been addressed in Steigmann et al. [71,68,70,69], see also dell'Isola et al. [20] for the application on panthographic structures. In Schulte et al. [60], a model of a woven fabric as presented in Steigmann [69] has been embedded into the Kirchhoff-Love shell theory. Therein, the in-plane flexural resistance of the fibers is taken into account in addition to first gradient anisotropic effects. Within this previous work, we were able to demonstrate on the basis of experimental results that a classical Cauchy continuum theory without higher-gradient contributions can only be adapted to a specific fiber orientation and load case, but never independently on the orientation. In contrast, the proposed formulation as generalized continuum with higher-gradient contributions allows for an independent modeling without recalibration of the material for the specific fiber direction. Further discussions on higher-order contributions for the constitutive modeling of composites have been presented in Asmanoglo and Menzel [7,8], using the framework as provided in the preliminary work of Spencer and Soldatos [67] and Soldatos [66]. Note that we will show here, that the formulation as proposed by dell'Isola et al. [20] can be recast into the formulation of Asmanoglo and Menzel [7,8].
Fiber reinforced polymers are exposed to various damage mechanisms. A large number of phenomenological and micromechanical approaches exists in literature for the modeling of damage within polymers used here as matrix material. To describe such phenomena, the material behavior at the microscale must be incorporated within the continuum formulation. In particular, the growth of microvoids prior to final rupture at the macroscale needs to be considered as rooted in the pioneering works of Gurson [30,31]. Therein, a macroscopic yield surface has been developed by homogenization of a porous representative volume element with assumed rigid plastic flow, which degrades with increasing void fraction. Later, this model was modified by Tvergaard [76], Tvergaard and Needleman [77], Needleman and Tvergaard [55], Leblond et al. [42], Nahshon and Hutchinson [54], Xue et al. [79], Li et al. [47] and Huespe et al. [36] to account for damage growth, where the yield criterion function has been extended by introducing new material parameters to account for nucleation and coalescence effects. In Hütter et al. [37] and Reusch et al. [58,59], non-local Gurson-models to overcome the nonphysical mesh sensitivity in the softening materials have been presented. Actual work on the application of this model towards polyamide, as often used for fiber reinforced thermoplastics, can be found e.g. in Selles et al. [64] and Cayzac et al. [16].
Phase-field methodologies have been successfully extended towards ductile fracture by coupling the gradient damage mechanism with models of elastoplasticity. The extension to finite strains is considered in [3,1,48,21,25,39] based on a rigorous variational principle. In Alessi et al. [5] a comparative study between different phase-field models of fracture coupled with plasticity is outlined. As already stated in the context of damage for polymers, see Gurson [31], pressure effects should be included in the modeling of failure in ductile materials to account for complex phenomena at the microscale, such as nucleation, growth and coalescence of microvoids. This has been observed experimentally in Gurland and Plateau [29]. To this end, Aldakheel et al. [4] extend the phase-field modeling of fracture towards porous finite plasticity to account for this complex phenomena at the microscale as well as for the final rupture at the macroscale. Thermomechanical extensions of the Gurson type model with brittle crack propagation in thermoelastic solids have been shown in Dittmann et al. [25] and Miehe et al. [50]. An extension towards finite strain thermo-porous-plasticity based on the phase-field approach has been recently developed in Dittmann et al. [22].
The paper is structured as follows: The governing equations for the coupled problem are outlined in Section 2, whereas algorithmic issues are addressed in Section 3 and 4. In Section 5, a variety of representative numerical examples is presented including a verification of the implementation. Finally, conclusions are drawn in Section 6.

Governing equations
In this section we give a brief summary of the fundamental equations for the modeling of thermomechanical damage in fiber reinforced composites. In order to provide a clear representation of mathematical operations used therein we define the gradient with respect to the reference and current configuration ∇(•) and ∇ x (•) of a vector field a as and of a second-order tensor field A as The divergence operator with respect to the reference configuration ∇ · (•) of a secondorder tensor field A and third-order tensor field A is defined as

Primary fields and state variables
We consider a fiber reinforced composite as a three-dimensional continuum body which occupies the domain B 0 ⊂ R 3 referred to as reference configuration. Assuming that the deformation of the fibers coincides with deformation of the matrix material, we introduce ϕ(X, t) : as a common field mapping at time t ∈ T = [0, T ] points X ∈ B 0 onto points x ∈ B of the current configuration. The material deformation gradient is defined by F = ∇ϕ(X, t) with its determinant J = det(F ) > 0. Moreover, the absolute temperature is introduced as a further common field representing the thermal state of the matrix as well as the fiber material.
Regarding the different damage behavior, the above common variables are supplemented by variables describing porous plasticity and ductile fracture of the matrix material and brittle fracture of the fiber material. For the matrix material we introduce the equivalent plastic strain and its dual, the dissipative resistance force along with the plastic deformation map and the void volume fraction where f 0 denotes the initial porosity. In addition, the crack phase-field is described by an order parameter where the value s = 0 refers to the undamaged and s = 1 to the fully ruptured state of the matrix material.
Assuming that the fiber reinforcement exhibits a woven structure, we introduce a dual crack phase-field for the fiber material. To be specific, the order parameters and describe the crack phase-field of the fiber aligned in L-direction and M -direction, respectively, where L and M are constant, orthogonal unit vector fields in the reference configuration.
The above introduced variables characterize a multifield setting for the formulation of temperature-dependent micro-and macromechanical damage in fiber reinforced composites based on seven independent fields the finite deformation map ϕ, the absolute temperature field θ, the equivalent plastic strain field α, the dissipative plastic resistance force r p , the crack phase-field s of the matrix material and the dual crack phase-field [s L , s M ] of the fiber material. Moreover, the Lagrangian plastic deformation map F p and the void volume fraction f will be condensed within the balance equations.

Kinematics and deformation measures
In a first step we derive the required deformation measures related to the matrix and fiber material. To this end, we apply a multiplicative split of the deformation gradient and its determinant as usual in non-linear elastoplasticity and obtain the elastic parts as F e = F (F p ) −1 and J e = J(J p ) −1 (13) which can also be defined in terms of the elastic parts of the principal stretches λ e a with a = {1, 2, 3} and the principal directions of the left and right stretch tensors n a and N a as F e = a λ e a n a ⊗ N a and J e = a λ e a .
Since the elastoplastic response of the matrix material relies on different mechanisms for the deviatoric and volumetric contributions, it is convenient to introduce the isochoric elastic parts of the principal stretches Following the ansatz proposed in Hesch & Weinberg [35], fracture insensitive parts of the elastic principal stretches are given as λ e a = (λ e a ) g(s) andλ e a = (λ e a ) g(s) , is an adjustable degradation function via the modeling parameter a g . Assuming that fracture requires a local state of tensile/shear deformation as considered, e.g. in Amor et al. [6] and Dittmann et al. [23,22], we define the elastic fracture insensitive part of the isochoric deformation gradient and the Jacobian determinant as Concerning the fiber material, we introduce as stretch of the respective fiber and as change of the angle between both fibers. Here l = λ Ll and m = λ Mm are deformed fiber configurations decomposed into fiber stretches and normalized fiber directions. To describe fiber bending, the gradients of the deformed fiber vectors, i.e. ∇l = ∇F L and ∇m = ∇F M , have to be taken into account. In particular, we consider ∇lL = λ L ∇lL + (∇λ L · L)l and ∇mM = λ M ∇mM + (∇λ M · M )m (20) which are projections of the fiber configuration gradients onto the initial fiber direction, cf. Asmanoglo and Menzel [7]. These expressions include terms related to stretch gradients of the fibers as well as fiber curvatures. We introduce the curvature measure for the fiber initially aligned in L-direction as and for the fiber initially aligned in M -direction as Assuming that the fiber material is brittle compared to the matrix material and that fiber rupture requires a local tensile state, we formulate fracture insensitive parts of the stretches as are degradation functions with modeling parameters a g L and a g M , cf. (16). Next, we formulate a corresponding measure related to the shear deformation of the fiber material as Note that this deformation measure is completely degraded in case of single fiber rupture even if the remaining fiber is undamaged. Eventually, fracture insensitive measures of the fiber curvatures read

Variational formulation
Next, we propose the constitutive framework for thermomechanical damage in fiber reinforced composites. To be specific, we introduce constitutive energetic and dissipative response functions based on the above definitions and derive the required relations and evolution laws to formulate the multifield variational problem.

Energetic response
The stored thermoelastic energy density of the composite material is defined by the functional where ζ ∈ [0, 1] is the volume fraction of the matrix material. The elastic contribution to the stored energy function of the matrix material is decomposed into volumetric and deviatoric parts Ψ e mat = Ψ e,iso mat F e (λ e 1 , λ e 2 , λ e 3 , s), θ + Ψ e,vol mat J e (λ e 1 , λ e 2 , λ e 3 , s), θ .
As a representative non-linear constitutive law, a modified Ogden material model with the associated strain energy density function and is used for the numerical examples. The parameters µ b and α b with b = {1, . . . , N } are related to the shear modulus and the parameters κ and β are related to the bulk modulus. Moreover, θ 0 is a reference temperature and the parameters and γ are related to the thermal expansion coefficient. Assuming that the fiber portion in both directions is identical, the corresponding elastic contribution of the fiber material is defined by where a and b are stiffness parameters related to stretch and shear of the fiber material and υ denotes the thermal expansion coefficient. Moreover, the stiffness tensor related to fiber curvature is given as taking into account a geometric dependency via the stiffness parameters c # and c ⊥ , see Section 5.1.1 for details. Next, the purely thermal contributions to the stored energy of the matrix and the fiber material are defined by and respectively. Therein, c mat and c fib are constant parameters representing specific heat capacities of the respective material.
The evolution of the stored thermoelastic energy is given by Regarding the partial derivatives therein, we introduce first relations related to the Kirchhoff stress τ = τ mat + τ fib as ∂λ e a n a ⊗ n a (35) and the higher-order stress of the fiber material as the driving force of the respective crack phase-field as and the specific entropy as Moreover, we introduce a dissipation function to account for a transfer of dissipated energy due to plastification and fracture into the thermal field, where d p denotes the Eulerian plastic rate of deformation tensor. The above relations are derived in a thermodynamically consistent manner by assuming that the dissipated energy is completely transfered into the thermal field, i.e. by setting ν pmat = ν fmat = ν f fib = 1. Note, however, that the plastic dissipation factor ν pmat is typically chosen in the range of 85% to 95% in the context of thermoplasticity, see e.g. [65,80,43]. In addition, based on experimental observations it is reasonable to set fracture dissipation factors to ν fmat < 1 and ν f fib < 1, see the discussion related to an energy transfer into the thermal field in [25,61] and the references therein.
To model the plastic and fracture mechanical response, we introduce an auxiliary functional as The plastic contribution Ψ p mat describes the response of isotropic strain-gradient hardening related to the matrix material. To be specific, we focus on the equivalent plastic strain α and its gradient ∇α with the particular form Here, l p is a plastic length scale related to a strain-gradient hardening effect and accounts for size effects to overcome the nonphysical mesh sensitivity of the localized plastic deformation in softening materials, as outlined in [1]. Moreover, y(α, θ) is an isotropic local hardening function taken from [16,64] and adapted to thermoplasticity following [65,57,18]. In particular, we use the saturation-type function with the three temperature-dependent material parameters y 0 > 0, y 1 ≥ 0 and y 2 ≥ 0 defined as Note that this formulation is typically applied for polyamide which is often used as matrix material of composite structures. Therein, the initial yield stress y 0 + y 1 determines the threshold of the effective elastic response, describes an initial hardening stage and y 1 (θ)exp[ω p1 α] allows for the simulation of large stretches of fibrils which leads to an abrupt increase of stress. This phenomenon is often called rheohardening. Moreover, ω p1 and ω p2 are saturation parameters and ω t0 , ω t1 and ω t2 are thermal hardening/softening parameters. Note that since we are only interested in the mean mechanical effects of the semi-crystalline matrix material, we consider a unified elastoplastic model with averaged material parameters taken from the multimechanism model in [16,64]. Next, we formulate fracture contributions for the matrix as well as the fiber material. Therefore, we approximate a sharp crack surface Γ • by a regularized functional † based on a specific crack regularization profile γ • defined per unit volume of the reference configuration and the fracture length scale l f• which controls the regularization. Concerning ductile fracture of the matrix material, we require that l p ≥ l f such that the regularized crack zone lies inside of the plastic zone. Using the regularization given in (45), the approximated fracture energy of the composite material reads (46) † The • indicates the matrix material or the respective fiber direction.
Here, g c• denotes the Griffith-type critical energy density required to create fracture within the respective material. For the matrix material, the critical energy density is decomposed additively into elastic and plastic contributions as where ω f is a modeling parameter. Summarized, the phase-field fracture contributions are given in terms of crack density functions as and Eventually, we obtain dissipative resistance forces of the plastic field and the respective crack phase-field via the variational derivatives of Ψ with respect to α and s • as and

Dissipative response
Regarding the porous elastoplastic material behavior, we consider a GTN type function [31,76,55] which implicitly defines the effective scalar stressσ :=σ(σ mat , f ) in terms of the Cauchy stress tensor σ mat = τ mat /J and the void volume fraction f Here, σ eq = 3/2 τ dev mat /J denotes the von Mises equivalent stress, p = 1 3 tr[τ mat /J] the pressure and q 1/2 are fitting parameters. Note that for q 1 = 0 the influence of the pressure and the void volume fraction vanishes, i.e.σ = σ eq . With the effective stressσ and the dissipative resistance force r p we define the plastic yield function as Focusing on void growth and thereby neglecting other influences such as void nucleation or void softening due to shear, the evolution form of the void growth readsḟ = (1 − f )tr[d p ]. Following [49], the current void volume fraction is given in terms of the plastic deformation as A plastic Lagrange multiplier λ p is introduced to enforce the Karush-Kuhn-Tucker conditions For the incorporation of the fracture mechanical behavior, we define crack threshold functions as where the energetic driving forces H • are bounded by crack resistance forces r f • dual to the crack phase-field variables s • . Similar to plasticity, we introduce fracture Lagrange multipliers λ f • to enforce the Karush-Kuhn-Tucker conditions of the respective crack phase-field Based on the concept of maximum dissipation and the set C = [σ mat , r p , , we define an extended dissipation potential and obtain a constrained optimization problem as where the Lagrange multipliers λ p and λ f • control the non-smooth evolution of plasticity and fracture, respectively. Then, the associated plastic evolution equations follows as and the evolution equation of the respective crack phase-field aṡ A penalty regularization of the Lagrange multipliers can be utilized as follows ‡ where η p , n p and η f• are material parameters which characterize the viscosity of plastification and crack propagation. Note that in the sense of continuum setting as defined in (26) and (41), the rates obtained in (59) and (60) are weighted by the respective volume fraction ζ and (1 − ζ)/2, respectively. ‡ The Macaulay brackets are defined by x = (x + |x|)/2.

Heat conduction
Regarding the heat transfer within the composite material, we introduce a relation for the Piola-Kirchhoff heat flux vector as which is known as Duhamel's law of heat conduction. The thermal conductivity tensor is defined as In case of fracture, the conduction degenerates locally such that we achieve a pure convection problem and the heat transfer depends mainly on the crack opening width of the matrix material. Here, we formulate the conductivity tensor K in terms of the phase-field parameter s. Moreover, K = ζK mat + (1 − ζ)K fib is a average conductivity parameter related to the composite material and K conv is a convection parameter.

Localization
To identify and collect the internal contributions to the boundaries of the mechanical field, i.e. the resulting bending moments and normal stress contributions, we derive the internal virtual work as The usage of the transformation along with integration by parts yields (66) and a second integration by parts related to the third term yields In a last step, we apply the divergence theorem for the second and third term such that we obtain Note that also contributions in tangential direction at the boundaries can be considered such that further integrations by parts incorporates the boundaries ∂ 2 B 0 and ∂ 3 B 0 which represent curves and points, see e.g. Schulte et al. [60] and Javili et al. [38]. Assuming that the principle of virtual work δW e,int − δW e,ext = 0 is valid with respect to the corresponding functional spaces of admissible solution and test functions, the external contribution can be formulated as where B is a given body force per unit volume of the reference configuration. Moreover, T andM are a surface traction and a surface torque acting at the mechanical Neumann boundaries Γ T 0 and Γ M 0 . Eventually, we obtain the local form of the mechanical problem as supplemented by boundary conditions with prescribed fieldsφ and ∇φN at the mechanical Dirichlet boundaries Γ ϕ 0 and Γ ∇ϕ 0 . As usual for fourth-order boundary value problems, the entire boundary is decomposed twice, i.e.

Coupled problem
Based on the derivations concerning the mechanical field within the previous section, the set of admissible test functions related to U is given as i.e. variations of the deformation, the absolute temperature, the equivalent plastic strain, the dual plastic resistance force, the crack phase-field of the matrix material and the variables of the dual crack phase-field of the fiber material, where their spaces are defined included within the Sobolev functional space of square integrable functions and derivatives H k with k ≥ 0. Then, the weak form of the coupled multifield problem reads where R is a given heat supply per unit volume of the reference configuration andQ is a heat supply across the thermal Neumann boundary Γ Q 0 . For each other field, homogeneous Neumann boundary conditions are applied and appropriate thermal Dirichlet boundary conditions are formulated in terms of prescribed temperatureθ, see Table  1. Note that we neglect inertia terms within the mechanical balance equation, i.e. we consider only quasi static problems. Additionally, internal conditions for the crack phase-field equations are given by ensuring that a fully broken state remains broken. The Karush-Kuhn Tucker conditions in (55) and (57), are evaluated by inserting local variables defined as Using local variables χ f• in (74), we demandṡ • ≥ 0 for thermodynamical consistency, avoiding a transfer of dissipated energy back into the mechanical field. This prevents healing effects, which may be taken into account as well. We can also set χ f• ≡ 1 and restrict only the fully broken state, i.e. we allow for healing until the respective crack phase-field reaches the value one.

Isogeometric discretization
Concerning the spatial discretization, the domain B 0 is subdivided into a finite set of non-overlapping elements e ∈ N such that Due to the incorporation of curvature contributions into the fiber material, the mechanical part of the variational problem requires approximation functions which are globally at least C 1 -continuous to satisfy ϕ h , δϕ h ∈ H 2 (B h 0 ). To meet this continuity requirement an isogeometric analysis approach which employs Non-Uniform Rational B-splines (NURBS) of order p α ≥ 2 can be applied. Accordingly, rational approximations of the deformed geometry ϕ and its variation δϕ are defined as respectively, where q A ∈ R 3 and δq A ∈ R 3 . Moreover, the approximations of the crack phase-fields s • and the temperature field θ as well as their variations δs • and δθ read and where s •,A , δs •,A , θ A , δθ A ∈ R. Introducing global shape functions R A : B h 0 → R associated with control points A ∈ I = {1, . . . , N}, NURBS based shape functions read
Next, the hardening variable α and its variation δα are approximated as and the dual driving force r p to the hardening variable and its variation δr p are approximated as where we make use of linear shape functions N i defined on the physical mesh representation of the NURBS geometry with nodes i ∈ J = {1, . . . , n} and the corresponding number of nodes n.
Remark: The more natural choice using the same NURBS shape functions for the approximation of the hardening variable α and the dual driving force r p leads to oscillations within both fields, indicating stability issues. The above described scheme using quadratic shape function R A and linear shape functions N i has shown to be stable and numerically robust within our numerical examples, cf. Dittmann et al. [23].
Inserting (92)-(94) along with (98) and (99) into (74) yields the semi-discrete set of coupled equations Therein, τ h , P h , η h and H h • are semi-discrete versions of the Kirchhoff stress tensor, the higher-order stress tensor, the local entropy and the phase-field driving forces obtained via the partial derivatives of the semi-discrete stored energy density cf. (26)- (39). D h int and Q h are semi-discrete definitions of the dissipation density and heat flux, cf. (40), (62) and (63). Moreover, the semi-discrete external contributions in (100) 1 and (100) 2 are formulated as and The coefficients of the matrices in (100) 3 and (100) 4 take the form whereas the matrices in (100) 5 are given by Eventually, the semi-discrete functions y h , Φ p,h and g h c denote the local hardening, the plastic yield and the critical fracture energy density, cf. (43), (53) and (47).

Temporal discretization
In a final step, the semi-discrete coupled problem (100) has to be discretized in time to obtain a set of non-linear algebraic equations to be solved via a Newton-Raphson method. Therefore, we subdivide the considered time interval T into a sequence of times t 0 , . . . , t n , t n+1 , . . . , T , where (•) n and (•) n+1 denote the value of a given physical quantity at time t n and t n+1 , respectively. Assume that the discrete set of state variables at t n given by {q A,n , θ A,n , α i,n , r p i,n , s A,n , s L,A,n , s M,A,n } and the local plastic deformation variable F p,h n at time t n are known and the time step size ∆t = t n+1 − t n is given. Then, the goal is to determine the corresponding fields at time t n+1 via the algorithmic approximation to the weak formulation (100) defined as (106) Therein, a full-discrete definition of the internal dissipation is given by Using small values for the plastic viscosity parameter η p , we obtain such that the internal dissipation can be recast as To solve the above multifield problem, we apply a staggered solution scheme, i.e. the displacement field along with the plastic and hardening fields {q A,n+1 , α i,n+1 , r p i,n+1 , F p,h n+1 }, the crack phase-fields s •,A,n+1 and the temperature field θ A,n+1 are solved successively. For the time integration of the plastic evolution equations, the construction of a return mapping algorithm is most crucial. Therefore, we define a trial state as § F e tr = F n+1 (F p n ) −1 (110) § For the sake of readability, we neglect the labeling of the spatial approximation in the following.
assuming that no further plastic deformation occurs within the time step. Based on this trial state, we evaluate the yield criteria (53). If Φ p tr ≤ 0, then the process is purely elastic and the elastic trial state is the solution. If on the other hand Φ p tr > 0, then the trial state is not admissible and a plastic correction is required. Therefore, we apply an exponential integration scheme regarding (59) which leads to λ e a,n+1 = λ e a,tr exp [−∆tλ p n+1 n a,n+1 ] with n a,n+1 = where σ mat,a,n+1 = (τ dev mat,a,n+1 + τ vol mat,a,n+1 )/J n+1 . Note that in contrast to standard von Mises plasticity n n+1 = n tr and n n+1 = 1, i.e. the plastic correction has to be performed by the Lagrange multiplier λ p n+1 as well as the components n a,n+1 which can be obtained by solving the non-linear relations via an internal Newton-Raphson iteration. In addition, the void volume fraction f n+1 is locally calculated by For further details on the return map algorithm see Dittmann et al. [22].

Numerical examples
In this section we investigate the accuracy and performance of the proposed formulation for endless fiber reinforced polymers. We start with a verification of the higher-order contributions of the fiber material by the means of two simple bending tests. Subsequently, a series of tensile tests demonstrates the capability of the proposed hybrid phase-field model to investigate different failure mechanisms for a prototypical fiber reinforced composite depending on the fiber configuration. This study is completed by thermal investigations on the damage behavior of the model and its impact on final failure. Without loss of generality we apply a Neo-Hookian model for the matrix material within all examples, i.e. we set b = 1 and α 1 = 2 in (28).

Bending Test
This first examples is dedicated to the verification of the higher-order, bending contributions of the fiber material. In particular, we investigate the in-plane bending behavior using a benchmark from Schulte et al. [60], originally used for the verification of gradient shell formulations, as well as the out-of-plane bending behavior using a four point bending test. Therefore, we consider a purely elastic behavior of the material, i.e. we neglect thermoplastic effects and fracture. We consider a rectangular geometry of size L × W × H = 10 mm × 1 mm × 0.5 mm discretized by 8 × 2 × 1 quadratic B-spline based elements. The left edge of the plate is clamped, whereas the right edge is subject to an external in-plane torque µ, see Figure 1.

In-plane bending test
A similar problem setting is applied to a Kirchhoff-Love shell formulation in [60], where the external torque is chosen to match an analytical reference solution. Here, we adopt the deformation result obtained in [60] to verify the parameter setting related to the bending stiffness. In particular, we assume that the plate consists of a single fiber bundle with a cross section of A = HW = 0.5 mm 2 and a tensile stiffness of E fib = 79000 N/mm 2 . The area moments of inertia of the fiber bundle with respect to the e 3 -axis and the e 2axis are given by I e 3 = HW 3 /12 = 0.0417 mm 4 and I e 2 = W H 3 /12 = 0.0104 mm 4 , respectively. Using these quantities we calibrate the bending stiffness parameters as c # = E fib I e 3 /A = 6583.3333 N and c ⊥ = E fib I e 1 /A = 1645.8333 N.
In Figure 2, the strain energy density is depicted for both the Kirchhoff-Love shell formulation as well as the proposed higher-order continuum formulation. Therein, we can observe the same homogeneous distributions which verifies the calibration of the bending stiffness parameter c # . Note that the parameter c ⊥ does not contribute to the simulation result, but will be investigated within the next example.  Next, the out-of-plane bending behavior of the fiber material is investigated using a four point bending test. Therefore, we consider again a rectangular geometry of size L × W × H = 125 mm × 25 mm × 0.5 mm discretized by 50 × 10 × 2 quadratic B-spline based elements. The composite material has a matrix volume ratio of ζ = 0.53 and the fibers are aligned in the ϑ = 0 • configuration, see Figure 3 for the details on the fiber orientation. The four point bending test as shown in Figure 4 leads to a pure out-ofplane bending deformation of the structure. In particular, we prevent the displacement in upward direction for the outer support points and prescribe a displacement in downward direction for the inner contact points. Additionally, the left support point is horizontally fixed, whereas we allow sliding for the other contact points. The material setting of the matrix material reads µ = 1630.4 N/mm 2 and α 1 = 2 for the deviatoric part and κ = 6250 N/mm 2 and β = −2 for the volumetric part, which corresponds to a Young's modulus of E mat = 4500 N/mm 2 and a Poisson's ratio of ν = 0.38. Two different settings of the fiber material properties are applied assuming a single layer of fibers over thickness direction. Firstly, we set the tensile stiffness of the fibers to a = E fib = 79000 N/mm 2 and the bending stiffness to c ⊥ = 0 N. Secondly, we set the tensile stiffness of the fibers to zero and adjust the bending stiffness as c ⊥ = E fib H 2 /12 = 1645.83 N.

Four point bending test
The applied bending stiffness of the continuum fiber model correlates to the out-ofplane bending stiffness for a shell model with the same high. As shown in the previous example, the proposed strain-gradient continuum formulation match the contributions of a gradient shell formulation, provided that the stiffness is chosen properly. Thus, if we resolve the thickness of sufficiently flat geometry with in the continuum model to obtain the same deformation as expected for the shell theory, a coincident bending behavior of the structure should result. Figure 5 shows the load deflection result for the investigated material settings. As expected, both results match in a good agreement, i.e. the tension/compression behavior of the continuum fiber model in this bending example can be described by the bending terms themselves. This is an important and well known result, as strain-gradient contributions emanate from a length-scale dependent microstructure and if this microstructure is already resolved by the first order continuum framework, the second-order contributions have to removed.  In this next example, we conduct a serious of tension tests to investigate the crack behavior of a prototypical roving glass composite material with different fiber configurations. Therefore, we consider a flat specimen of size L × W × H = 125 mm × 25 mm × 2 mm. Figure 6 and 9 show the geometry in the reference configuration along with the applied boundary conditions and the fiber configurations. The outer areas of length 20 mm are subject to Dirichlet boundary conditions. To be specific, one flap is fixed and the other flap is moved by a displacement rate of 0.5 mm/s within a quasi-static simulation setting neglecting inertia effects. The computational mesh consists of 2432 quadratic NURBS elements. The material setting of the composite is summarized in Table 2. We assume a quadratic cross section of the fibers with A fib = 0.0025 mm 2 and obtain a bending stiffness of c ⊥ = c # = E fib A fib /12 = 16.46 N.

Unidirectional fiber reinforcement
We first analyze the behavior of a unidirectional reinforced composite material, with Figure 6. The load deflection results for isothermal simulations at θ = 293 K are shown in Figure 7. Therein, crack initialization and final rupture of the fiber material are indicated by and •, respectively. In addition, crack initialization and final rupture of the matrix material are indicated by and ×, respectively. In Figure 8, the crack phase field results of the fiber material and the matrix material are depicted along with results of the plastic strain field for the marked points. Note that black marker indicate states without fiber fracture, as the  specimen is already fully broken.
For a fiber orientation of ϑ = 0 • , the fibers account for most of the load transfer due to the different Young's modulus and fracture abruptly in the center of the specimen . Subsequently, the matrix material undergoes plastification and ductile fracture due to an abrupt load rearrangement . Note that the resulting high strain rates lead to a pronounced viscoplastic behavior within the matrix material which is controlled by the viscous regularization parameter.
Concerning the unidirectional 10 • fiber configuration, the fibers start to crack near the clamping zones which is additionally driven by the bending contribution to the crack diving force. This process is slowed down due to the hardening behavior of the matrix material . At the state , the fibers are fully ruptured and the matrix material starts to fracture in the same region.
For a fiber orientation of ϑ = 20 • , brittle fracture of the fibers starts again near the clamping zones . However, a more pronounced plastification and thus hardening of the matrix material occur such that the fiber and matrix material undergo final rupture nearly at the same deformation state .
Applying a fiber orientation of ϑ = 30 • , a direct load transfer between both boundaries by the fibers is not possible since fibers which are clamped at the lower end do not reach the upper clamping zone. Hence, the load has to be transferred towards the matrix material leading to higher plastification and ductile fracture at the center of the specimen . Note that the fibers begin to fracture only in small areas near the clamping zones .
For fiber orientations of ϑ = [40 • , 65 • , 90 • ], the fiber material does not fracture due to small loads acting in fiber direction -. Instead, the matrix material undergoes suitable plastification and subsequently ductile fracture leading to failure orthogonal to the fiber orientations. Concerning the ϑ = 90 • fiber configuration, the fibers controls the necking which can be observed by comparing the deformation with the results obtained for pure matrix material . This can also be observed by a slightly higher stiffness within the load deflection results before crack initiation, whereas the results are nearly identical up to a displacement of u = 40 mm.

Bidirectional fiber reinforcement
Next, we investigate the same tension test using a bidirectional reinforced material with fiber orientations of ϑ = [0 • , 10 • , 20 • , 30 • , 45 • ] as shown in Figure 9. The load deflection results for isothermal simulations at θ = 293 K are depicted in Figure 10. Again, crack initialization and final rupture of the fiber material are indicated by and •, respectively. The final rupture of the matrix material is indicated by ×. Crack phase-field results of the fiber and matrix material as well as results of the plastic strain are depicted in Figure  11. Note that only phase-field results of the fiber aligned in ϑ-direction is plotted.
For fiber orientations of ϑ = [0 • , 10 • , 20 • ], the bidirectional reinforced material shows a similar behavior compared to the corresponding unidirectional reinforced counterparts. The additional orthogonal fiber merely accounts for a higher necking resistance -. Moreover, the orthogonal fiber configuration restrict a relative movement between the respective fibers leading to a slightly more stiff material behavior. This has to be investigated in terms of experimental measurements which is out of the scope of present work.
Applying fiber orientations of ϑ = [30 • , 45 • ], this effect becomes more pronounced as can be observed in Figure 10. As already discussed, the additional, orthogonal oriented Figure 11: Tensile Test (bidirectional). Results of the fiber crack phase-field (first row) as well as the plastic strain field (second row) and crack phase-field of the matrix material (third row). The results are shown for the different deformation states and fiber configurations marked in Figure 10. fiber counteract the necking behavior due to the Poisson effect of the matrix material such that fractures within the matrix material emerges near the clamping zones and not in the center of the specimen -.

Thermal investigation
Eventually, we investigate the temperature dependency of the proposed model. Therefore, we reuse the tension test with a unidirectional fiber reinforcement as shown in Figure 6 and apply a fiber orientation of ϑ = 30 • . Figure 12 shows the load deflection result for isothermal simulations using temperatures of θ = [253, 273, 293] K. The corresponding crack phase-field results of the fiber and matrix material as well as results of the plastic strain are depicted for the last deformation step in Figure 13. As already observed previously, for θ = 293 K the matrix material undergoes large plastic deformations followed by fiber fracture in small areas near the clamping zones and finally the matrix material undergoes ductile fracture at the center of the specimen. Lower temperatures increase the yield stress of the matrix material leading to a higher elastic energy and thus an earlier, less ductile fracture behavior of the matrix material. Note that for isothermal simulations with θ = [253, 273] K the matrix material fails before any fiber cracks occur.

Conclusions
The non-linear framework presented in this work allows for a comprehensive investigation of damage and fracture in fiber reinforced polymers. The combination of a secondgradient theory, a novel hybrid phase-field model and a temperature dependent GTNtype plasticity model provides a numerical framework which is able to describe different failure mechanisms in detail. This approach allows for improvements in the design of such composite materials since we are able to predict fiber and matrix failure and their sequence dependent on the fiber orientation. Moreover, due to the fully-coupled, thermomechanical approach we can optimize the fiber orientation for specific loads and thermal states. Several numerical tests conducted within this work have demonstrate the capability of the proposed framework to investigate such a complex behavior including the growth of microvoids, plasification and necking, crack initiation and propagation within the composite material and its components, respectively.