A strain rate dependent and anisotropic failure material model for short fiber reinforced plastics based on energy density accounting for local fiber orientation

Abstract

A failure model for SFRP for FEM simulations is developed to describe the strain rate dependency, the influence of the local fiber orientation and of the stress state on the failure behavior. The material is considered as a continuum while internally calculating the micro-mechanics analytically. The described micro-mechanics are based on experimental observations and on analyzation with numerical studies. In particular the strain rate dependent delamination of fibers and matrix is incorporated in the model. The distortion energy density is defined as the driving value for failure and estimated by the model. This is achieved with the analytic solution by Eshelby for the stress field in the matrix and by introducing an additional phase for the plasticly deformed volume. The validation on characterization specimens as well as component test demonstrates that the influence of strain rate, fiber orientation, and stress state on the failure behavior can be described with only one material parameter, the critical distortion energy density.

Keywords

Composites micro-mechanics failure modelling strain rate effects debonding

Introduction

The importance of short fiber reinforced plastics (SFRPs) and their modelling in regards of crash safety has undergone a great increase in the last years in the automotive and aviation sector on account of lightweight aspects and cost efficiency.^1,2 A major hurdle for an increased application even in safety relevant structures is that the prediction of its behavior is complex due to the anisotropy and its dependency on the local fiber orientation, the strain rate, and the stress state. In Figure 1 this is demonstrated by characterization experiments on PP-GF30 displaying the anisotropy, and the influence of the strain rate and the stress state. The experiments were performed by Sun et al.,³ Sommer et al.⁴ and Lienhard.⁵

Figure 1.

The experimentally observed deformation and failure behavior of glass fiber reinforced polypropylene at different strain rates for uniaxial tensile test in fiber direction and transverse to the fiber direction, shear tests and a punch test for biaxial tension.

The individual representation of the fibers and the matrix material is numerically expensive for finite element simulations of larger crash scenarios. Hence, SFRP is often described as a continuum while only using the local fiber orientation as internal parameters.

Recently, advances have been made in modelling the deformation behavior of SFRP. For example, the model presented by Schulenberg et al.⁶ is capable of describing the reversible and irreversible deformation with their dependency on strain rate and the stress state while also considering the effect of the fiber orientation. Wherein homogenization methods are used, and the fiber orientation serves as weights for establishing the anisotropy. The model by Reithofer et al.⁷ applies the mean field theory of Mori and Tanaka and an incremental flow rule for the matrix material. Both approaches are efficient and flexible for predicting the deformation of SFRP. But for failure a phenomenological approach is chosen defining a separate parameter set for each possible situation (strain rate, stress state, and fiber orientation). Therefore, many parameters are needed which in turn results in high experimental expenses.

On the other hand, much effort has been made to find accurate failure models that are based on the underlying physics. A common approach is to calculate analytically the damage of the matrix in each time increment. Thereby the fibers are often split in unidirectional groups.^8–10 As will be discussed below, the fiber-matrix delamination plays an important role in the failure mechanics. To consider this, models have been developed for small cracks on fiber ends¹¹ or by adding a third phase for the cracks.¹² Other approaches^13–15 based on the work of Qu¹⁶ with an imperfect attachment between fiber and matrix allow a displacement on the surface and use a subsequent homogenization.

Although the listed models can predict the failure of SFRP within specific boundaries, there is still no model capable of describing its dependency of the strain rate or additional influences as the stress state and the fiber orientation. Such a model is developed and validated in this paper.

Material and mechanics

The here presented failure model is developed on the example of a glass fiber reinforced polypropylene with 30% fiber volume fraction (PP-GF30), a fiber diameter of 19 μm and a nominal length of 2 mm. A transfer to others short fiber reinforced thermoplastics is in principle possible.

To have a physically accurate basis the significant micro-mechanics must be described in a failure model. An apparent mechanic is the strain rate dependency of each component of the composite. Glass fibers deform linear elastically with a neglectable strain rate effect on the stiffness.¹⁷ They fail rapidly and can be describe as brittle whereby the failure strain is higher for dynamic loads.¹⁸

Pure polypropylene however deforms nonlinear early in a loading and can accumulate large irreversible deformations and is therefore usually described with plastic flow theory. This plasticity shows a significant strain rate effect just as the failure strains.^19,20 It is yet unknown how exactly the presence of fibers change the behavior of polypropylene, in terms of change of crystallinity. Wang and Liu²¹ detected a different nucleation rate close to carbon fibers. Zhao et al.²² were able to prove the existence of a 200–300 nm thick layer close to the fiber with a lower hardness. Though with analytic homogenization methods and applying the E modulus of not reinforced polypropylene one can predict the composite stiffness correctly and therefore show that this layer is neglectable. The nonlinear behavior cannot be analyzed in this way and is experimentally very challenging to characterize due to the small distance between the fibers. Consequently, in the following the behavior of the matrix material is assumed to be the same as the not reinforced polypropylene.

Equally as important is the interaction of the components and specifically the delamination. 1992 Kander and Siegmann²³ have discovered, that the delamination of glass fibers and polypropylene is strain rate dependent. Even before that the quasi-static detaching was characterized,²⁴ and more recently some tests were performed at higher strain rates.^25,26 The experiments prove the influence of the strain rate on delamination with higher critical stresses at higher strain rates.

As one can observe in the macroscopic experiments on PP-GF30 (Figure 1) the failure stress and strain increases with increasing strain rates. This holds for both loading direction in respect to the fiber orientation (tension in 0° and 90°) as well as for all applied stress states (uniaxial tension, shear, and biaxial tension). In other words, the absorbed energy until failure, the distortion energy

E_{d i s} = \int F d L,

increases with higher strain rates, where $F$ is the applied force and $L$ the displacement. Additionally, Huberth et al.²⁷ were able to develop a technique to test SFRP on a microscopic scale (with a specimen length of 1.1 mm). This can help to uncover the present micro-mechanics and the origin of specific influences on the failure behavior. As shown in Figure 2, the deformed material volume during a tensile test is higher at higher strain rates, analogical to the distortion energy. Therefore, the hypothesis can be made that the ratio between the energy absorption and the deformed volume $V$ , the distortion energy density

W_{d i s} = \int \frac{F d L}{V}

at failure is constant for all strain rates. The material model developed is based on this hypothesis and the distortion energy density is used as failure indicator.

Figure 2.

In the SFRP micro specimens just before failure the deformed material appears brighter. Strain rate of 0.0007/s (above) and 250/s (below) were applied.²⁷

Furthermore, this hypothesis was recently validated by Schweiger et al.²⁸ with a finite element study on the strain rate dependent microscopic failure mechanism of SFRP, modelling the fibers and the matrix material separately. The simulations showed that modelling the strain rate dependent delamination is necessary to capture the effect of the strain rate dependent deformed volume. And moreover, with this method it was possible to reproduce the results of the microscopic experiments based on the real fiber orientation. Figure 3 shows that the deformed volume is higher at high strain rates (lower) than at lower strain rates (upper). The analysis of the simulated distortion energy density showed a similar value at failure for both strain rates, providing another argument for describing failure with this value.

Figure 3.

Plastic strain distribution for the SFRP micro tensile test resolving the fibers and the matrix material modeling the strain rate dependent delamination and producing a strain rate dependent deformed volume. The same strain rate were applied as experimentally (0.0007/s above and 250/s below).²⁸

Failure model development

A predictive failure model for SFRP must describe the significant micro-mechanisms as discussed above. To achieve this, its core capability must be the prediction of the stress field in the matrix material. Eshelby²⁹ has presented a well-known analytic solution for ellipsoidal inclusions of different stiffness exposed to any outer deformation. As the solution is restricted to linear elastic behavior an extension for nonlinear behavior will be established later. For the Eshelby solution eigenstrains $ε_{k l}^{e i g e n}$ can be found that generate an equivalent strain field as the inclusion. This strain $ε_{i j}^{o u t e r} (x)$ in the matrix at a point $x$ is linear in the eigenstrains with the fourth order tensor $D_{i j k l} (x)$ as

ε_{i j}^{o u t e r} (x) = D_{i j k l} (x) ε_{k l}^{e i g e n},

which is in turn dependent on the geometry of the fiber, in particular the semi-axes

a_{n}

of the ellipsoid. For the specific expression of

D_{i j k l} (x)

see the work of Mura³⁰ or Jiang and Bhandari.³¹ With the matrix stiffness

C_{g h i j}

the local stress field caused by the fiber can be expressed as

σ_{g h} (x) = C_{g h i j} D_{i j k l} (x) ε_{k l}^{e i g e n} .

For prolate fibers ( $a_{1} > a_{2} = a_{3}$ ) the solution is of a closed form and can be implemented directly. This was done for glass fiber reinforced polypropylene ( $E_{g l a s s} = 73 G P a, ν_{g l a s s} = 0.23, E_{P P} = 2.1 G P a, ν_{P P} = 0.35$ ) applying a tension in fiber direction and compared with a numerical finite element (FE) simulation using the commercial software LS DYNA (see Figure 4). For this fully integrated hexahedral volume elements were used with a size of about 2.5 μm. 1/8 of the fiber was simulated with prohibited normal movement in the symmetry planes. The matrix material was modeled with an extent of 3 times the fiber length and 10 times the fiber diameter to approximate the infinite matrix assumed by Eshelby. To generate the following graphics of the stress field based on the Eshelby solution this analytic solution and its extensions were realized in python. The same was implemented in LS DYNA as a user material model for the validation simulations. As the difference of the stress fields does not exceed 4% of the local stress it can be concluded that it is negligible, and its origin is the numerical approximations.

Figure 4.

A comparison of the stress field close to an ellipsoidal fiber (black) between the analytic Eshelby solution (upper) and a numeric FE simulation (middle). Displayed is the 11-component for a tension in fiber direction as well as their relative difference (lower). The upper left quadrant is shown in the symmetry plane at x₃=0.

Even though the formulars are numerically efficient, the solution only holds with major assumptions. Firstly, only one inclusion in an infinite matrix is described. Secondly, the material behavior of the inclusion and the matrix is linear elastic. And thirdly, the fiber-matrix-contact is assumed to be perfect. As in SFRP the fibers are usually densely spaced, and the matrix deforms nonlinearly these assumptions do not hold and extensions to the solution must be made.

The influence of adjacent fibers on the stress field in the matrix can be accounted for by additive superposition of the stress fields caused by each fiber individually viewing the fibers as disturbance in the matrix (see Figure 5). The half fiber distance $d$ is used to compute the local fiber volume fraction

f = \frac{\frac{1}{6} π a_{1} {a_{2}}^{2}}{(a_{1} + d) {(a_{2} + d)}^{2}} .

f

is assumed to be known, this can be solved to

d = 2 Q \cos (\frac{\cos^{- 1} (R Q^{- 3})}{3}) - \frac{a_{1} + 2 a_{2}}{3} w i t h

Q = \frac{| a_{2} - a_{1} |}{3} a n d R = \frac{π a_{1} {a_{2}}^{2}}{12 f} - \frac{{(a_{1} - a_{2})}^{3}}{27} .

Figure 5.

A additive superposition of the stress field of adjacent fibers each calculate according to the Eshelby solution. The half fiber distance d is determined by the local fiber volume fraction. Displayed is the 11-component of the stress field for a tension in fiber direction in the symmetry plane at x₃=0.

The dashed line in Figure 5 shows the boundaries of a volume element that is repeated in all directions. The same behavior in all of them can be assumed. Therefore, only one of these representative volume elements is considered in the user material model for each element or integration point in the finite element simulation.

The local fiber orientation can be accounted for by rotating the adjacent fibers by a certain degree, as described in the following (see Figure 6). Firstly, to reduce the number of parameters necessary to describe the fiber orientation, the deformation is transferred to the local coordinate system in which the second order fiber orientation tensor³² $A$ is diagonal as

A = (\begin{array}{c} A_{1} & 0 & 0 \\ 0 & A_{2} & 0 \\ 0 & 0 & A_{3} \end{array}) w i t h A_{1} \geq A_{2} \geq A_{3} a n d A_{1} + A_{2} + A_{3} = 1 .

As the fibers oriented in the thickness direction for SFRP is known to be neglectable in common structures^4,33 the third eigenvalue is here assumed to be zero and

A

is simplified to

A = (\begin{array}{c} A_{1} & 0 & 0 \\ 0 & 1 - A_{1} & 0 \\ 0 & 0 & 0 \end{array}) .

Figure 6.

The overlapping stress field in the symmetry plane caused by rotated fibers around an angle of φ=13° based on the local fiber orientation.

The fiber orientation tensor represented in the analytic model $A^{a n a l y t i c}$ can be calculated by a closed integration over the distribution function $ψ (p, α)$ with the spatial direction $p$

A_{i j}^{a n a l y t i c} = \oint p_{i} p_{j} ψ (p, α) d p .

Thereby the distribution function $ψ$ depends on the rotation angle $α$ for the adjacent fibers. The rotation angle $α$ can be calculated, so that $A^{a n a l y t i c}$ matches the local fiber orientation in the component. For these calculations the local fiber volume fraction and the local fiber orientation are needed. Such can be provided for example by molding simulations or CT scans.

For the stress field a sum over all adjacent fibers is added with

σ_{g h} (x) = C_{g h i j} D_{i j k l} (x) ε_{k l}^{e i g e n} + \sum_{n}^{a d j . f i b e r s} R (α_{n}) C_{g h i j} D_{i j k l} (x_{r e l, n}) ε_{k l, r o t, n}^{e i g e n} .

Thereby for rotated fibers the applied load must be rotated resulting in different eigenstrains $ε_{k l, r o t, n}^{e i g e n}$ and the stress field must be rotated back into the global coordinate system by the rotation matrix $R (α_{n})$ . The rotation matrix is depending on the rotation angle $α$ which in turn depends on the local fiber orientation tensor $A$ . Also, the stress field must be evaluated at relative space coordinates $x_{r e l, n}$ considering the position of the adjacent fiber.

After a method has been established to describe the elastic behavior with the Eshelby solution, in the following an extension is developed for the plastic deformation of the matrix material and its effect on the stress within the matrix. Calculating a complete stress field as in Figures 4–6 consists of 10,000–100,000 points where the stress function has to be evaluated. Considering that this would only represent one finite element at one time step in the FE simulation, it is apparent that a more efficient method is needed to approximate the area, where plastic flow occurs.

Therefore, the stress is calculated point wise on straight lines equidistant on the fiber surface minimizing the difference to the yield stress by applying an extended Golden Section search.^34,35 Hereby the shape of the area with stresses above the yield stress can be approached as shown in Figure 7 (upper) for a tension 90° to the fiber direction.

Figure 7.

The optimal positioning of an ellipsoid within the plasticly deformed area by point wise optimization (upper) and the von Mises stress field caused by the plastic deformation of this volume (lower).

The effect of the plastic deformation of the matrix material on the stress field can again be calculated by the Eshelby solution if the plasticly deformed region is regarded as an additional inclusion with an E modulus deviating from the matrix material. This is done analogical to calculating the stress field caused by the glass fiber. But here the inclusion has a lower stiffness than the matrix material. To apply the Eshelby solution the plasticly deformed area is approximated by an ellipsoid of maximal volume, while not violating the border of the flow area (see Figure 7, upper).

To approximate the effect of the plasticity on the stress field with the Eshelby solution the load on the inclusion must be analyzed. A mean stress $σ^{p l V}$ in the so defined plasticly deformed volume can be calculated point wise. This stress exceeds the yield surface (see Figure 8). Beyond the yield surface the material deviates from the linear elastic behavior. Thus, only the stress above the yield stress must be considered for calculating the load on the inclusion. Under the assumption of proportional loading the specific yield stress $σ^{y i e l d}$ can be calculated by a normal projection to the yield surface.

Figure 8.

The relation between the present stress in plasticly deformed volume, the normal projection to the yield surface and the corresponding stress on the yield surface in the principal stress space.

Following the stress above the yield stress, here called overstress

σ^{o v e r} = σ^{p l V} - σ^{y i e l d},

is the effective stress acting on the plasticly deformed volume. Applying this in the Eshelby solution as outer load, the stress field caused by the plastic deformation can be calculated (see Figure 7, lower). With the accumulated plastic strain and the hardening curve, a tangential E modulus can be calculated and assigned as

C_{g h i j}^{p l V}

for the stiffness of the plasticly deformed volume. In the following time step the material model takes this stress field additively into account as well as the stress field caused by the plasticly deformed volumes at adjacent fibers:

σ_{g h} (x) = C_{g h i j} D_{i j k l} (x) ε_{k l}^{e i g e n} + \sum_{n}^{a d j . f i b e r s} R (α_{n}) C_{g h i j} D_{i j k l} (x_{r e l, n}) ε_{k l, r o t, n}^{e i g e n} + \sum_{n}^{p l a s t i c V o l u m e s} R (α_{n}) C_{g h i j}^{p l V} D_{i j k l}^{p l V} (x_{r e l, n}) ε_{k l, r o t, n}^{e i g e n, p l V}

A FE simulation was performed for comparison to the analytic results. Thereby the same boundary conditions were applied as above (see Figure 4) with a stress state dependent plasticity and the material parameters from Table 1. As depicted in Figure 9 in the numeric simulation the plastic flow lowers the stress close to the fiber. The analytic model does qualitatively the same. Therefore, the restricted form of the plastic region is in some cases a major simplification. But this is still necessary for the application of the Eshelby solution and can be improved in future works by implementing additional ellipsoids as plastic regions.

Table 1.

The material parameter based on literature which were not adjusted during the following simulations with the flow curves as $f (ε_{p l}) = σ_{y i e l d} + Q_{1} * (1 - e^{- C_{1} ε_{p l}}) - Q_{2} * (1 + e^{C_{2} ε_{p l}})$ .

Elastic		Yield stress		Flow curve		Delamination
E − glass	73 GPa	$σ_{y i e l d, t e n s i o n}^{q u a s i - s t a t i c}$	18 MPa	$Q_{1}$	6.2 MPa	$σ_{d e l a m i n a t i o n}^{q u a s i - s t a t i c}$	38 MPa
E − PP	1.3 GPa	$σ_{y i e l d, c o m p r e s s i o n}^{q u a s i - s t a t i c}$	19 MPa	$C_{1}$	140	$σ_{d e l a m i n a t i o n}^{d y n a m i c}$	75 MPa
$ν$ – glass	0.23	$σ_{y i e l d, s h e a r}^{q u a s i - s t a t i c}$	13 MPa	$Q_{2}$	0.5 MPa
$ν$ – PP	0.35	$σ_{y i e l d, b i a x i - t e n s i o n}^{q u a s i - s t a t i c}$	15 MPa	$C_{2}$	4
		$σ_{y i e l d, t e n s i o n}^{200 / s}$	31 MPa

Figure 9.

A comparison of the von Mises stress field between the analytic approximation (left) and the numeric validation simulation (right). Thereby the linear elastic for each modelling is displayed at the top. In the modelling accounting for plasticity (below) the regions of plastic deformation are marked.

A plastic material model for the matrix material must be implemented to predict incrementally the tangent E modulus and the accumulation of the plastic strain. Here a description analogical to the successful material model in LS DYNA “semi-analytical model for polymers” developed by Haufe et al.³⁶ is chosen. Thereby, the yield surface is based on two Drucker-Prager cones³⁷ while a dependency on the strain rate and on the multiaxiality of the load can be specified. By an analytic quadratic formulation, the $C_{0}$ -differentiability is ensured. The plastic strain $ε_{p l}$ is incrementally accumulated using a common plasticity algorithm where the flow condition is solved by a Newtons method with an associative flow rule and a radial return. This formulation is interchangeable, and not crucial to the main idea of the here presented failure model.

With an approximation of the extent of the plasticly deformed volume the distortion energy density

ρ = \int \frac{σ^{p l V}}{V_{p l} / V_{R V E}} d ε_{p l}

can be calculated. As

σ \cdot ε

already has the unit of an energy density it is divided by the fraction the plasticly deformed volume

V_{p l}

has on the volume of the whole RVE

V_{R V E}

. To convert the matrix

σ^{p l V}

to the scalar

σ^{p l V}

the von Mises stress was applied as norm. In the simulation failure occurs when a user defined critical value

ρ^{k r i t}

is reached.

As discussed above the debonding of fiber and matrix is a significant microscopic mechanism and is vital for describing the strain rate dependent failure of SFRP. The unmodified Eshelby solution assumes a perfect contact raising the need for another extension. The stress field of a partially delaminated fiber is analyzed with a numeric FE simulation in Figure 10 (upper) subjected to tension in fiber direction and with material parameters and elements as in the simulation of Figure 4. The delamination was modeled with a so called tied-break-contact in LS DYNA, enabling the nodes of fiber and matrix to separate if a critical stress is reached and with parameters discussed below (see Table 1). If one compares this with the stress field of a fiber which is shorter by the same length in Figure 10 (lower) the similarity is noticeable. The local difference (Figure 11, upper) is not negligible because the stress maxima is slightly shifted. By aligning the stress maxima one can show, that their shape and absolute values match (Figure 11, lower). The correct shape is more impactful than the specific location because by this the plastic regions will also be slightly shifted but the same size.

Figure 10.

A comparison of the von Mises stress field caused by a partially delaminated fiber (upper) and by a by the same length shorten fiber (lower) with numerical simulations.

Figure 11.

The relative difference of the stress fields of Figure 10 and the relative difference whereby location of the stress maxima were matched.

Consequently, the shortening of the effective fiber length is a good approximation of its debonding with regards to its stress field. This is implemented in the material model by checking the stress on the fiber surface point wise. If the delamination strength is exceeded the effective length of the fiber is reduced by that length. The delamination strength can be defined depending on the strain rate.

As the glass fibers are of cylindrical form and the Eshelby solution approximates the inclusions with ellipsoids, the influence of the fiber shape on the delamination must be discussed. In order to quantify this FE simulation are compared for quasi-static tension in fiber direction. First with cylindrical fibers and secondly with ellipsoids. In Figure 12 the stress field in the matrix material and the delamination is shown for three timesteps. For this simulation the linear elastic and delamination parameters discussed below (in Table 1) were applied. The delamination was again enabled by a tied-break-contact. At $t = 5 s$ the cylindrical fibers have already detached at the fiber ends because of the stress spike at the edge while with ellipsoidal fibers the contact is still intact. At $t = 15 s$ the ellipsoids have started detaching as well but the delamination with cylinders hast still advanced further. At $t = 30 s$ the delamination of the ellipsoid has caught up and is equal to the cylindric fibers.

Figure 12.

Finite element simulations for delamination with cylindrical and ellipsoidal fibers at three timesteps.

This exemplary simulation shows that the fiber shape has an influence on the delamination. The edge of the cylindric fibers causes an earlier detachment. But it also shows ellipsoidal fibers show qualitatively the same behavior and the delamination length even catches up in time. It can therefore be said that using ellipsoidal shape for the fibers is a simplification, but a justifiable and necessary simplification to find an analytical solution.

The failure model as described above was implemented as a user routine in LS DYNA. For the elastic and plastic material behavior the model of Schulenberg et al.⁶ was used. Thereby a strain rate dependent elasticity can be described and the local fiber distribution is taken into account with the homogenization methods of Mori and Tanaka³⁸ as well as from Advani and Tucker.³⁹ The plasticity is strain rate dependent as well and the yield surface consists of three Hill surfaces⁴⁰ weighted by the eigenvalues of the local fiber orientation tensor with an exponential flow rule.

The numerical performance of the user material model is noticeable different compared with standard material models in LS DYNA with anisotropic behavior. The simulation speed is about 5–10 times slower from the onset of plastic deformation. But with a revision of the structural setup of the failure model, greater speed would be possible, as its initial implementation placed more emphasis on clarity than numerical efficiency.

Validation on characterization experiments

In the following the established failure model for SFRP is applied to characterization experiments. The critical distortion energy density $ρ^{k r i t}$ is optimized and the capability to model the failure for a wide range of specific loads is validated. First the experimental and numerical set up is explained, where after the global stress–strain behavior is compared, followed by an analysis of the local prediction of the failure model.

Tensile specimens in molding direction and 90° to the molding direction, tensile shear specimens and punch specimens for biaxial loading were manufactured from molded plates with 2.8 mm thickness of a PP-GF30. The specific geometries and measurements are shown in Figure 13. Each specimen is subjected to a quasi-static and a dynamic load resulting in eight configurations. A more detailed description can be found in Sommer et al.⁴

Figure 13.

Geometries and dimensions of tensile, shear and punch specimens. Colored grey are firmly clamped areas.

The specimens are represented by a finite element mesh of volume elements of hexahedral form with an edge length of 0.5 mm and simulated in LS DYNA using the proposed failure model. The same deformation velocities as in the experiments were applied. The force per area and elongation normalized on $L_{0}$ was analyzed as stress–strain curves.

The failure model requires some material parameters that can be found in the literature and which were not fitted during the simulations (see Table 1). Firstly, the linear elastic parameter of fiber and matrix, which were direct manufacturer information. In addition, the plastic behavior of the matrix must be described including the influence of strain rate and stress triaxiality. These flow curves were obtained by simulating the characterization experiments on pure polypropylene in Lienhard and Schulenberg⁴¹ and Sommer et al.⁴ and fitting the material model parameters (see Figure 14). Here it was only necessary to fit the yield stress, but in general the whole flow curve can be adjusted.

Figure 14.

The strain rate and stress state dependent flow curves optimized to fit the Tension, Compression, Shear and Plane Strain Tension experiment from Lienhard and Schulenberg⁴¹ and Sommer et al.⁴

Most importantly the failure model features a strain rate dependent delamination. Whereas a critical stress for quasi-static loadings is well established in literature,^24,25 for high loading velocities the direct experimental characterization is still a challenge. With the experimental work of Tanoglu et al.²⁵ on the delamination of glass fiber reinforced epoxy-amine and the molecular dynamic simulations of Zhang et al.²⁶ on glass fiber reinforced polypropylene at very high strain rates the critical stress can be approximated to 75 MPa at 250/s with an uncertainty of about 10 MPa. Thus, this should be verified experimentally in the future.

The only parameter which was optimized was the critical energy density $ρ_{c r i t}$ at which failure occurs and the elements are deleted. A value of 135 Nmm/mm³ produced the overall best results.

In Figures 15–17 the global stress–strain curves of the characterization experiment are compared with the curves of the simulation. As the deformation model has been discussed in detail before,^6,42 here the focus will be on the failure exclusively. At first glance it can be seen that the failure can be predicted well by the failure model for all deformation types and strain rates. The tensile test in 0° and 90° to the molding direction in Figure 15 show that the anisotropy is captured correctly, even in combination with the strain rate dependency.

Figure 15.

Comparison of the experimental (solid lines) and the simulated results (dashed lines) for the tension in 0° and 90° to the molding direction. A quasi-static load at 0.001/s and a dynamic load at 200/s is applied.

Figure 16.

Comparison of the experimental (solid lines) and the simulated results (dashed lines) for shear test. A quasi-static load at 0.001/s and a dynamic load at 200/s is applied.

Figure 17.

Comparison of the experimental (solid lines) and the simulated results (dashed lines) for the punch test. A quasi-static load with a punch velocity of 0.01°mm/s and a dynamic load with 1400°mm/s is applied.

In Figure 16 one can see that the high deformations by applying shear can be reached as well. In the quasi-static case, the experiment shows a slight softening before failure, which cannot be described by the simulation. This suggests that implementing a stress-damage coupling would solve this divergence.

For the punch tests in Figure 17, additionally to the force–displacement curves the optically visible initial failure on the surface of the specimen is marked by points for the experiment and the simulation respectively. In the dynamic loading for these points a deviation is found. This could be explained in several ways. The exact time of the initial failure is often hard to determine as small cracks might occur which are not observed optically. Therefore, the initial failure might be experimentally detected too late, and the real initial failure happened earlier at a punch displacement closer to 2 mm. Another influence to consider is, that the punch is the only test with deformation components out of plane where bending occurs. This can lead to different behavior. Or these results might hint on a dependency of the critical energy density on the stress multiaxiality. Dynamic test on a biaxial tensile machine would lead to answers.

For a better understanding of the modeled microscopic effects in the following an analyzation is done of the calculation made by the failure models. Thereby, elements are selected where the failure in the specimen initiates. In Figure 18 the prediction for these critical elements for four simulations are displayed, the tension in 0° and 90° for the quasi-static and the dynamic load. As for the stress fields discussed previously, the upper left of a fiber is shown in the symmetry plane for $x_{3} = 0$ . Thereby the initial fiber with its full length is colored in grey and with black lines the delamination and analogical the shortening of the fiber is shown. In the (white) matrix ellipsoids of different brightness depict the prediction for the expanding plasticly deformed volume over time.

Figure 18.

The local predictions of the failure model in elements of first failure. The quasi-static case (upper) is compared with dynamic loading (lower) for the tension in fiber direction (left) and transverse (right). The fibers and their delamination are depicted (black and grey) as well as the development of the plasticly deformed volumes over time.

The failure model only calculates stresses within the borders of the symmetry element marked by the dashed line. Nonetheless for clarity neighboring fibers and their associated plasticly deformed volumes are displayed as well.

For the quasi-static tension in 0° a significant delamination is calculated and as a result of that the plasticly deformed volume stays small and localized close to the fiber ends. Together with the plasticly deformed volume of a neighboring fiber and the delaminated fibers, they form a band (from top to bottom) just as observed in the experiments (see Figure 2). A similar behavior is predicted in the quasi-static tension in 90°. Here the delamination leads to a small plasticly deformed volume stretching along the fiber length and forming a band in this direction (from left to right).

Applying a dynamic load for tension neither in 0° nor in 90° a delamination occurs with the failure model. Subsequently the size of the plasticly deformed volumes is higher covering almost the whole matrix and the critical energy density is reached later.

Validation on component experiments

After applying the failure model to characterization experiments, where a specific deformation is provoked, its transferability to more real crash scenarios is tested with complex deformations. The form of the molded component is a U-profile of 75 mm width, 245 mm length, and 60 mm height and a wall thickness of 2.4 mm. A full description of the experiments can be found in Sommer et al.⁴

For the simulations the local fiber orientation was obtained from a molding simulation and mapped to the FE model. Two load cases are tested. The first is a 4-Point-Bending, where one side is restricted to ensure a reproducible behavior (see Figure 19). As before in the characterization experiments, two strain rates are tested and simulated with the calibrated user model. For both velocities the material model is accurate in the initial stiffness and the first failure marked by the first rapid decline in the punch force. Dynamically the force drops to zero at 10 mm whereas quasi-statically a plateau is observed till about 20 mm. This can be predicted by the model as well, though the level of the force plateau is calculated slightly too low.

Figure 19.

A component test with 4-point-bending, quasi-static and with a punch velocity of 1250°mm/s, comparing the experiment (solid lines) and the simulations (dashed lines).

In a second experiment half of the component was crushed in its length direction, applying 1 mm/s and 1350 mm/s. As shown in Figure 20 the material model can predict the behavior for this test as well, reaching a maximal force of about ^8kN quasi-statically and ^13kN dynamically before the crack spans through the whole component at 6 mm and a sliding off with a slow decent in the punch force.

Figure 20.

The quasi-static and dynamic crush test on half of the component, comparing the experiment (solid lines) and the simulations (dashed lines).

The failure model is not only capable of describing the influences on failure (anisotropy, strain rate, and stress multiaxiality) separately as it occurs in characterization experiments. But by these component tests it has been proven, that complex load cases can be predicted with mixed stress states and general fiber orientation distributions on realistic components.

Discussion and outlook

A failure model for SFRP was developed for simulating crash scenarios with FEM by describing the SFRP as continuum and taking the local fiber orientation into account as an inner variable. Based on the micro-mechanics and the delamination of the fibers in particular, the model is highly efficient with only one fit parameter. The material behavior can be predicted in a wide range of load cases.

Motivated by the analysis of macroscopic and microscopic experiments the distortion energy normalized on the deformed volume was incorporated as the driving value for failure. Based on the Eshelby solution the stress field in the matrix was calculated, and with modifications the local fiber volume fraction and local fiber orientation were considered.

As an approximation for the plasticly deformed volume an additional ellipsoidal phase within the matrix was introduced. By that, the Eshelby solution can be applied once again. This represents a major simplification and could be improved for example by implementing multiple volumes to better fit the real shape of the deformed volume. Lastly the strain rate dependent delamination was modeled by reducing the effective fiber length.

The predictive nature of the model was validated on multiple experiments. Thereby, all material parameters were taken from literature and only the critical distortion energy density was optimized. For all tests the strain rate dependency was predicted correctly. By means of the tension in and transvers to the fiber direction the anisotropy due to the fiber orientation was tested successfully. Different stress states can be described as was demonstrated on tension, shear, and biaxial tension. And by two component tests, it was proven, that the model can also predict complex loadings.

A so fare not tested influence on the failure is the degree of alignment of the fibers. This can range from uniformly distributed in all directions to unidirectional aligned. With the failure model in the current state such tests would be possible, but (micro-) experiments on material with specific fiber alignment would be needed for validation.

Furthermore, the universality of the failure model could be tested by transferring it to other SFRP, for example to other fiber lengths, fiber volume fractions or other fiber materials as carbon or even plant-based fibers and thereby possibly improving the modeling of a wide range of materials.

Footnotes

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work was supported by AiF program of the BMWi for collaborative industrial research, a research association for vehicle engineering (Forschungvereinigung Automobiltechnik e. V. (FAT), Behrenstrasse 35, 10117 Berlin) grant no (IGF-No. 20418 N).

ORCID iD

Timo Schweiger

References

Jayan

. Fibre reinforced composites: constituents, compatibility, perspectives and applications. 1st ed. Duxford: Woodhead Publishing, 2021, pp. 1–24.

Lester

Nutt

. Composite materials: advantages and cost factors. Los Angeles, CA, USA: Elevated Materials, 2018.

Sun

. Entwicklung einer methode zur crashsimulation von langfaserverstärkten thermoplast (LFT) bauteilen auf basis der faserorientierung aus der formfüllsimulation. . Technical report. Berlin: Forschungsvereinigung Automobiltechnik e.V. (FAT), 2016.

Sommer

. Effiziente charakterisierung und modellierung des anisotropen versagensverhaltens von LFT für crashsimulation. Technical report. Berlin: Forschungsvereinigung Automobiltechnik e.V. (FAT), 2020.

Lienhard

. Auswirkungen von dehnrate und lokaler temperaturentwicklung auf das crashverhalten langfaser-verstärkter thermoplaste. PhD Thesis. Karlsruhe: Karlsruher Institut für Technologie, 2018.

Schulenberg

Seelig

Andrieux

, et al. An anisotropic elasto-plastic material model for injection-molded long fiber-reinforced thermoplastics accounting for local fiber orientation distributions. J Compo Mat 2017; 51: 2061–2078.

Reithofer

Fertschej

Hirschmann

, et al. *MAT_4A_MICROMEC - micro mechanic based material model. In: 14th German LS-DYNA Conference, Bamberg, Germany, 2016.

Nguyen

Khaleel

. A mechanistic approach to damage in short-fiber composites based on micromechanical and continuum damage mechanics descriptions. Compo Sci Tech 2004; 64: 607–617.

Nghiep Nguyen

Bapanapalli

Kunc

, et al. Prediction of the elastic-plastic stress/strain response for injection-molded long-fiber thermoplastics. J Compo Mat 2009; 43: 217–246.

10.

Kammoun

Doghri

Adam

, et al. First pseudo-grain failure model for inelastic composites with misaligned short fibers. Compo: Appl Sci Manufac 2011; 42: 1892–1902.

11.

Taya

Mura

. On stiffness and strength of an aligned short-fiber reinforced composite containing fiber-end cracks under uniaxial applied stress. J Appli Mech 1981; 48: 361–367.

12.

Meraghni

Benzeggagh

. Micromechanical modelling of matrix degradation in randomly oriented discontinuous-fibre composites. Composites Science and Technology 1995; 55: 171–186.

13.

Nguyen

Tucker

Khaleel

. A mechanistic approach to matrix cracking coupled with fiber-matrix debonding in short-fiber composites. J Eng Mat Tech 2005; 127: 337–350.

14.

Nghiep Nguyen

Kunc

. An elastic-plastic damage model for long-fiber thermoplastics. Int J Damage Mech 2010; 19: 691–725.

15.

Othmani

Delannay

Doghri

. Equivalent inclusion solution adapted to particle debonding with a non-linear cohesive law. Int J Solids Struct 2011; 48: 3326–3335.

16.

. The effect of slightly weakened interfaces on the overall elastic properties of composite materials. Mech Mat 1993; 14: 269–281.

17.

Zhang

Zou

Hao

, et al. Laboratory test on dynamic material properties of annealed float glass. Int J Protec Struc 2012; 3: 407–430.

18.

Arao

Taniguchi

Nishiwaki

, et al. Strain-rate dependence of the tensile strength of glass fibers. J Mater Sci 2012; 47: 4895–4903.

19.

Ferry

. Viscoelastic properties of polymers. New York and Chichester and Brisbane and Toronto and Singapore: John Wiley and Sons, 1980.

20.

Al-Maliky

. Strain rate behaviour of thermoplastic polymers. PhD thesis. Loughborough: Loughborough University, 1997.

21.

Wang

Liu

. Transcrystallization of polypropylene on carbon fibres. Polymer 1997; 38: 4715–4718.

22.

Zhao

Singh

Korach

. Effects of environmental degradation on near-fiber nanomechanical properties of carbon fiber epoxy composites. Compo: Appli Sci Manufac 2009; 40: 675–678.

23.

Kander

Siegmann

. The effect of strain rate on damage mechanisms in a glass/polypropylene composite. J Compo Mat 1992; 26: 1455–1473.

24.

Gupta

Mittal

Sharma

, et al. Some studies on glass fiber-reinforced polypropylene. Part II: mechanical properties and their dependence on fiber length, interfacial adhesion, and fiber dispersion. Polym Compos 1989; 10: 16–27.

25.

Tanoglu

McKnight

Palmese

, et al. A new technique to characterize the fiber/matrix interphase properties under high strain rates. Compo: Appli Sci Manufac 2000; 31: 1127–1138.

26.

Zhang

Jiang

Chen

, et al. The effect of temperature and strain rate on the interfacial behavior of glass fiber reinforced polypropylene composites: a molecular dynamics study. Polymers (Basel) 2019; 11(11): 1766.

27.

Huberth

Lienhard

Ragupathi

, et al. Dynamic micro-testing over a large range of strain rates for homogenous and heterogenous material characterization. In: 9th international Conference on computational Methods and Experiments in materials characterisation, Lisbon, Portugal, 2019: WIT Press, pp. 143–150.

28.

Schweiger

. Investigation of strain rate dependent microscopic failure mechanisms in short fiber reinforced plastics using finite element simulations. Submitted to. J Compo Mat 2022; 21. JCM-22-1076.

29.

Eshelby

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

30.

Mura

. Micromechanics of defects in solids. Mechanics of elastic and inelastic solids. Dordrecht: Kluwer Academic Publishers, 1998.

31.

Jiang

Bhandari

. Pressure variations among rheologically heterogeneous elements in Earth’s lithosphere: a micromechanics investigation. Earth Planet Sci Lette 2018; 498: 397–407.

32.

McGee

McCullough

. Characterization of fiber orientation in short-fiber composites. J Appl Phys 1984; 55: 1394–1403.

33.

Pötsch

Michaeli

. Injection molding: an introduction. München and Cincinnati. Hanser, 2008.

34.

Wiegand

Petrinic

Elliott

. An algorithm for determination of the fracture angle for the three-dimensional Puck matrix failure criterion for UD composites. Compo Sci Tech 2008; 68: 2511–2517.

35.

Kiefer

. Sequential minimax search for a maximum. Proc Am Math Soc 1953; 4: 502–506.

36.

Haufe

. A semi-analytical model for polymers subjected to high strain rates. In: 5th European LS-DYNA users conference, Birmingham, UK, 2005, pp. 25–26.

37.

Junginger

. Charakterisierung und Modellierung unverstärkter thermoplastischer Kunststoffe zur numerischen Simulation von Crashvorgängen. PhD Thesis. München: Univ. der Bundeswehr, 2002.

38.

Mori

Tanaka

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

39.

Advani

Tucker

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

40.

Hill

. Constitutive modelling of orthotropic plasticity in sheet metals. J Mech Phys Soli 1990; 38: 405–417.

41.

Lienhard

Schulenberg

. Strain rate dependent multiaxial characterization of long fiber reinforced plastic. Compo Part B: Eng 2018; 141: 164–173.

42.

Schulenberg

. Mehrskalenmodellierung langfaserverstärkter Thermoplaste - experimentelle und numerische Untersuchungen zum Deformations- und Versagensverhalten. PhD Thesis. Karlsruhe: Karlsruher Institut für Technologie, 2018.