Abstract
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.
Similar content being viewed by others
1 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 [54, 55], see also the work of Germain [30], Toupin [76, 77] as well as Eringen [29]. A general non-linear framework on higher-order theories has been proposed by Javili et al. [40]. Based on these most general formulations, specific mechanical problems like elastic nets and woven fabrics have been addressed in Steigmann [70, 71], Steigmann and dell’Isola [72] and Steigmann and Pipkin [73], see also dell’Isola et al. [21] for the application on panthographic structures. Schulte et al. [62] and Duong et al. [27] , a model of a woven fabric as presented in Steigmann [71] 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 [69] and Soldatos [68]. Note that we will show here, that the formulation as proposed by dell’Isola et al. [21] 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 [32, 33]. 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 [78], Tvergaard and Needleman [79], Needleman and Tvergaard [57], Leblond et al. [44], Nahshon and Hutchinson [56], Xue et al. [81], Li et al. [49] and Huespe et al. [38] 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. Hütter et al. [39] and Reusch et al. [60, 61], 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. [66] and Cayzac et al. [16].
A model for initiation and propagation of ductile fracture using the damage plasticity theory has been proposed in Bai and Wierzbicki [9]. Moreover, a large number of purely phenomenological approaches exist in literature describing ductile fracture in the context of continuum damage mechanics, see Lemaitre [46, 47], Lemaitre and Chaboche [48], Steinmann et al. [74], de Borst et al. [19], Besson [10], Enakoutsa et al. [28], Larsson et al. [43], Seabra et al. [65] and Brünig and Brenner [14] and [13]. To calculate fracture with complex crack topologies within an efficient computational environment, phase-field methodologies have been developed and applied to multiphysical environments, see, among many other, Miehe et al. [53], Hesch et al. [35, 36], Borden et al. [11], Kuhn et al. [42], Verhoosel and De Borst [80], Paggi and Reinoso [58], Teichtmeister et al. [75], Zhang et al. [83], Dittmann et al. [25, 26], Heider and Markert [34], Bryant and Sun [15] and Aldakheel et al. [2] for brittle fracture.
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 [1, 3, 22, 26, 41, 50] based on the variational principle. 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 [33], 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 [31]. 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. [26] and Miehe et al. [52]. An extension towards finite strain thermo-porous-plasticity based on the phase-field approach has been recently developed in Dittmann et al. [23].
The paper is structured as follows: the governing equations for the coupled problem are outlined in Sect. 2, whereas algorithmic issues are addressed in Sects. 3 and 4. In Sect. 5, a variety of representative numerical examples is presented including a verification of the implementation. Finally, conclusions are drawn in Sect. 6.
2 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 \(\nabla (\bullet )\) and \(\nabla _x(\bullet )\) of a vector field \(\varvec{a}\) as
and of a second-order tensor field \(\varvec{A}\) as
The divergence operator with respect to the reference configuration \(\nabla \cdot (\bullet )\) of a second-order tensor field \(\varvec{A}\) and third-order tensor field \(\varvec{\mathfrak {A}}\) is defined as
respectively. Moreover, the double contractions of a second-order and third-order tensor, i.e. \(\varvec{a}=\varvec{A}:\varvec{\mathfrak {A}}\) and \(\varvec{b}=\varvec{\mathfrak {A}}:\varvec{A}\), are defined as \([\varvec{a}]_k=[\varvec{A}]_{ij}:[\varvec{\mathfrak {A}}]_{ijk}\) and \([\varvec{b}]_i=[\varvec{\mathfrak {A}}]_{ijk}:[\varvec{A}]_{jk}\), respectively.
2.1 Primary fields and state variables
We consider a fiber reinforced composite as a three-dimensional continuum body which occupies the domain \(\mathcal {B}_0\subset \mathbb {R}^3\) referred to as reference configuration. Assuming that the deformation of the fibers coincides with deformation of the matrix material, we introduce
as a common field mapping at time \(t\in \mathcal {T}=[0,T]\) points \(\varvec{X}\in \mathcal {B}_0\) onto points \(\varvec{x}\in \mathcal {B}\) of the current configuration. The material deformation gradient is defined by \(\varvec{F}=\nabla \varvec{\varphi }(\varvec{X},t)\) with its determinant \(J=\det (\varvec{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 \({\mathfrak {s}}=0\) refers to the undamaged and \({\mathfrak {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 \(\varvec{L}\)-direction and \(\varvec{M}\)-direction, respectively, where \(\varvec{L}\) and \(\varvec{M}\) are constant, orthogonal unit vector fields within the body 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 \(\varvec{\varphi }\), the absolute temperature field \(\theta \), the equivalent plastic strain field \(\alpha \), the dissipative plastic resistance force \(r^{\mathrm {p}}\), the crack phase-field \({\mathfrak {s}}\) of the matrix material and the dual crack phase-field \([{\mathfrak {s}}_{\mathrm {L}},{\mathfrak {s}}_{\mathrm {M}}]\) of the fiber material. Moreover, the Lagrangian plastic deformation map \(\varvec{F}^{\mathrm {p}}\) and the void volume fraction f will be condensed within the balance equations.
2.2 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
which can also be defined in terms of the elastic parts of the principal stretches \(\lambda _a^{\mathrm {e}}\) with \(a=\{1,2,3\}\) and the principal directions of the left and right stretch tensors \(\varvec{n}_a\) and \(\varvec{N}_a\) as
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 and Weinberg [37], fracture insensitive parts of the elastic principal stretches are given as
where \(g=a_{\mathrm {g}}((1-{\mathfrak {s}})^3-(1-{\mathfrak {s}})^2)-2(1-{\mathfrak {s}})^3+3(1-{\mathfrak {s}})^2\) is an adjustable degradation function via the modeling parameter \(a_{\mathrm {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, 24], 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 \(\varvec{l}=\lambda _{\mathrm {L}}\tilde{\varvec{l}}\) and \(\varvec{m}=\lambda _{\mathrm {M}}\tilde{\varvec{m}}\) 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. \(\nabla \varvec{l}=\nabla \varvec{F}\varvec{L}\) and \(\nabla \varvec{m}=\nabla \varvec{F}\varvec{M}\), have to be taken into account. In particular, we consider
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 \(\varvec{L}\)-direction as
and for the fiber initially aligned in \(\varvec{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
where \(g_{\mathrm {L}}=a_{\mathrm {gL}}((1-{\mathfrak {s}}_{\mathrm {L}})^3 -(1-{\mathfrak {s}}_{\mathrm {L}})^2)-2(1-{\mathfrak {s}}_{\mathrm {L}})^3+3(1-{\mathfrak {s}}_{\mathrm {L}})^2\) and \(g_{\mathrm {M}}=a_{\mathrm {gM}}((1-{\mathfrak {s}}_{\mathrm {M}})^3-(1-{\mathfrak {s}}_{\mathrm {M}})^2) -2(1-{\mathfrak {s}}_{\mathrm {M}})^3+3(1-{\mathfrak {s}}_{\mathrm {M}})^2\) are degradation functions with modeling parameters \(a_{\mathrm {gL}}\) and \(a_{\mathrm {gM}}\), 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
2.3 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.
2.3.1 Energetic response
The stored thermoelastic energy density of the composite material is defined by the functional
where \(\zeta \in [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
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 \(\mu _b\) and \(\alpha _b\) with \(b=\{1,\hdots ,N\}\) are related to the shear modulus and the parameters \(\kappa \) and \(\beta \) are related to the bulk modulus. Moreover, \(\theta _0\) is a reference temperature and the parameters \(\epsilon \) and \(\gamma \) 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 \(\upsilon \) 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_{\perp }\), which can be interpreted as the in-plane and out-of-plane bending stiffness, see Section 5.1.1, [7, 62], and [20] 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_{\mathrm {mat}}\) and \(c_{\mathrm {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 \(\varvec{\tau }=\varvec{\tau }_{\mathrm {mat}}+\varvec{\tau }_{\mathrm {fib}}\) as
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 \(\varvec{d}^{\mathrm {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 \(\nu _{\mathrm {p_{mat}}}=\nu _{\mathrm {f}_{\mathrm {mat}}}=\nu _{\mathrm {f}_{\mathrm {fib}}}=1\). Note, however, that the plastic dissipation factor \(\nu _{\mathrm {p}_{\mathrm {mat}}}\) is typically chosen in the range of 85–95% in the context of thermoplasticity, see e.g. [45, 67, 82]. In addition, based on experimental observations it is reasonable to set fracture dissipation factors to \(\nu _{\mathrm {f}_{\mathrm {mat}}}<1\) and \(\nu _{\mathrm {f}_{\mathrm {fib}}}<1\), see the discussion related to an energy transfer into the thermal field in [26, 63] and the references therein.
To model the plastic and fracture mechanical response, we introduce an auxiliary functional as
The plastic contribution \(\widehat{\Psi }_{\mathrm {mat}}^{\mathrm {p}}\) describes the response of isotropic strain-gradient hardening related to the matrix material. To be specific, we focus on the equivalent plastic strain \(\alpha \) and its gradient \(\nabla \alpha \) with the particular form
Here, \(l_{\mathrm {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(\alpha ,\theta )\) is an isotropic local hardening function taken from Cayzac and Saï [16] and Selles [66] and adapted to thermoplasticity following Simo and Miehe [67], Reis et al. [59] and Da Costa Mattos et al. [18]. In particular, we use the saturation-type function
with the three temperature-dependent material parameters \(y_0>0\), \(y_1\ge 0\) and \(y_{2}\ge 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, \(y_{2}(\theta )(1-\mathrm {exp}[-\omega _{\mathrm {p2}}\alpha ])\) describes an initial hardening stage and \(y_{1}(\theta )\mathrm {exp}[\omega _{\mathrm {p1}}\alpha ]\) allows for the simulation of large stretches of fibrils which leads to an abrupt increase of stress. This phenomenon is often called rheo-hardening. Moreover, \(\omega _{\mathrm {p1}}\) and \(\omega _{\mathrm {p2}}\) are saturation parameters and \(\omega _{\mathrm {t}0}\), \(\omega _{\mathrm {t}1}\) and \(\omega _{\mathrm {t}2}\) 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, 66]. Next, we formulate fracture contributions for the matrix as well as the fiber material. Therefore, we approximate a sharp crack surface \(\Gamma _{\bullet }\) by a regularized functionalFootnote 1
based on a specific crack regularization profile \(\widehat{\gamma }_{\bullet }\) defined per unit volume of the reference configuration and the fracture length scale \(l_{\mathrm {f}_{\bullet }}\) which controls the regularization. Concerning ductile fracture of the matrix material, we require that \(l_{\mathrm {p}}\ge l_{\mathrm {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
Here, \(g_{\mathrm {c}_{\bullet }}\) 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 \(\omega _{\mathrm {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 \(\widehat{\Psi }\) with respect to \(\alpha \) and \({\mathfrak {s}}_{\bullet }\) as
and
2.3.2 Dissipative response
Regarding the porous elastoplastic material behavior, we consider a GTN type function [33, 57, 78] which implicitly defines the effective scalar stress \(\bar{\sigma }:=\bar{\sigma }(\varvec{\sigma }_{\mathrm {mat}},f)\) in terms of the Cauchy stress tensor \(\varvec{\sigma }_{\mathrm {mat}}=\varvec{\tau }_{\mathrm {mat}}/J\) and the void volume fraction f
Here, \(\sigma _{\mathrm {eq}}=\sqrt{3/2}\, \Vert \varvec{\tau }^{\mathrm {dev}}_{\mathrm {mat}} / J \Vert \) denotes the von Mises equivalent stress, \(p=\frac{1}{3}\mathrm {tr}[\varvec{\tau }_{\mathrm {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. \( \bar{\sigma }=\sigma _{\mathrm {eq}}\). With the effective stress \(\bar{\sigma }\) and the dissipative resistance force \(r^{\mathrm {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 \(\dot{f}=(1-f)\mathrm {tr}[\varvec{d}^{\mathrm {p}}]\). Following Miehe et al. [51], the current void volume fraction is given in terms of the plastic deformation as
A plastic Lagrange multiplier \(\lambda ^{\mathrm {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 \(\mathcal {H}_{\bullet }\) are bounded by crack resistance forces \(r^{\mathrm {f}}_{\bullet }\) dual to the crack phase-field variables \({\mathfrak {s}}_{\bullet }\). Similar to plasticity, we introduce fracture Lagrange multipliers \(\lambda ^{\mathrm {f}}_{\bullet }\) to enforce the Karush–Kuhn–Tucker conditions of the respective crack phase-field
Based on the concept of maximum dissipation and the set \(\mathfrak {C}=[\varvec{\sigma }_{\mathrm {mat}},r^{\mathrm {p}},\mathcal {H}-r^{\mathrm {f}},\mathcal {H}_{\mathrm {L}}-r^{\mathrm {f}}_{\mathrm {L}},\mathcal {H}_{\mathrm {M}}-r^{\mathrm {f}}_{\mathrm {M}},\lambda ^{\mathrm {p}},\lambda ^{\mathrm {f}},\lambda ^{\mathrm {f}}_{\mathrm {L}},\lambda ^{\mathrm {f}}_{\mathrm {M}}]\), we define an extended dissipation potential and obtain a constrained optimization problem as
where the Lagrange multipliers \(\lambda ^{\mathrm {p}}\) and \(\lambda ^{\mathrm {f}}_{\bullet }\) 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 as
A penalty regularization of the Lagrange multipliers can be utilized as followsFootnote 2
where \(\eta _{\mathrm {p}}, \mathrm {n}_\mathrm {p} \) and \(\eta _{\mathrm {f}_{\bullet }}\) 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 \(\zeta \) and \((1-\zeta )/2\), respectively.
2.3.3 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 \(\varvec{K}\) in terms of the phase-field parameter \({\mathfrak {s}}\). Moreover, \(K=\zeta K_{\mathrm {mat}}+(1-\zeta )K_{\mathrm {fib}}\) is a average conductivity parameter related to the composite material and \(K^{\mathrm {conv}}\) is a convection parameter.
2.3.4 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 yieldsFootnote 3
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 \(\partial ^2\mathcal {B}_0\) and \(\partial ^3\mathcal {B}_0\) which represent curves and points, see e.g. Schulte et al. [62] and Javili et al. [40]. Assuming that the principle of virtual work \(\delta W^{\mathrm {e,int}}-\delta W^{\mathrm {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 \(\varvec{B}\) is a given body force per unit volume of the reference configuration.
Eventually, we obtain the local form of the mechanical problem as
supplemented by boundary conditions
with prescribed fields \(\bar{\varvec{\varphi }}\) and \(\nabla \bar{\varvec{\varphi }}\varvec{N}\) at the mechanical Dirichlet boundaries \(\Gamma _0^{\varphi }\) and \(\Gamma _0^{\nabla \varphi }\). As usual for fourth-order boundary value problems, the entire boundary is decomposed twice, i.e. \(\Gamma _0=\Gamma _0^{\varphi }\cup \Gamma _0^{T}\) with \(\Gamma _0^{\varphi }\cap \Gamma _0^{T}=\emptyset \) and \(\Gamma _0=\Gamma _0^{\nabla \varphi }\cup \Gamma _0^{M}\) with \(\Gamma _0^{\nabla \varphi }\cap \Gamma _0^{M}=\emptyset \). For details related to the enforcement of the gradient condition given in (72)\(_2\) see e.g. Schuß et al. [64].
2.3.5 Coupled problem
Based on the derivations concerning the mechanical field within the previous section, the set of admissible test functions related to \(\mathfrak {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 as
included within the Sobolev functional space of square integrable functions and derivatives \(H^k\) with \(k\ge 0\). Then, the weak form of the coupled multifield problem reads
where \(\mathcal {R}\) is a given heat supply per unit volume of the reference configuration and \(\bar{Q}\) is a heat supply across the thermal Neumann boundary \(\Gamma _0^{Q}\). For each other field, homogeneous Neumann boundary conditions are applied and appropriate thermal Dirichlet boundary conditions are formulated in terms of prescribed temperature \(\bar{\theta }\), 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 \(\chi _{\mathrm {f}_{\bullet }}\) in (75), we demand \(\dot{{\mathfrak {s}}}_{\bullet }\ge 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 \(\chi _{\mathrm {f}_{\bullet }}\equiv 1\) and restrict only the fully broken state, i.e. we allow for healing until the respective crack phase-field reaches the value one.
3 Isogeometric discretization
Concerning the spatial discretization, the domain \(\mathcal {B}_0\) is subdivided into a finite set of non-overlapping elements \(e\in \mathbb {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 \(\varvec{\varphi }^{\mathrm {h}},\delta \varvec{\varphi }^{\mathrm {h}}\in \mathcal {H}^2(\mathcal {B}_0^{\mathrm {h}})\). To meet this continuity requirement an isogeometric analysis approach which employs Non-Uniform Rational B-splines (NURBS) of order \(p_{\alpha }\ge 2\) can be applied. Accordingly, rational approximations of the deformed geometry \(\varvec{\varphi }\) and its variation \(\delta \varvec{\varphi }\) are defined as
respectively, where \(\varvec{q}_{A}\in \mathbb {R}^3\) and \(\delta \varvec{q}_{A}\in \mathbb {R}^3\). Moreover, the approximations of the crack phase-fields \({\mathfrak {s}}_{\bullet }\) and the temperature field \(\theta \) as well as their variations \(\delta {\mathfrak {s}}_{\bullet }\) and \(\delta \theta \) read
and
where \({\mathfrak {s}}_{\bullet ,A},\delta {\mathfrak {s}}_{\bullet ,A},\theta _{A},\delta \theta _{A}\in \mathbb {R}\). Introducing global shape functions \(R^A:\mathcal {B}_0^{\mathrm {h}}\rightarrow \mathbb {R}\) associated with control points \(A\in \mathcal {I}=\{1,\,\hdots ,\,\mathfrak {N}\}\), NURBS based shape functions read
where \(B^{i_{\alpha }}\) are univariate non-rational B-splines defined on a parametric domain which is subdivided by the knot vector \([\xi _{1}^{\alpha },\xi _{2}^{\alpha },\hdots ,\xi _{\mathfrak {N}_{\alpha }+p_{\alpha }+1}^{\alpha }]\), \(\mathfrak {N}=\mathfrak {N}_1\mathfrak {N}_2\mathfrak {N}_3\). The recursive definition of a single univariate B-spline is given as follows
beginning with
Moreover, \(w_{\varvec{j}}\) are corresponding NURBS weights. For further details on the construction of NURBS based shape functions as well as the construction of local refinements related to the IGA concept, see e.g. Cottrell et al. [17], Bornemann and Cirak [12], Hesch et al. [36] and Dittmann [22].
Next, the hardening variable \(\alpha \) and its variation \(\delta \alpha \) are approximated as
and the dual driving force \(r^{\mathrm {p}}\) to the hardening variable and its variation \(\delta r^{\mathrm {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\in \mathcal {J}=\{1,\,\hdots ,\,\mathfrak {n}\}\) and the corresponding number of nodes \(\mathfrak {n}\).
Remark: The more natural choice using the same NURBS shape functions for the approximation of the hardening variable \(\alpha \) and the dual driving force \(r^{\mathrm {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. [24].
Inserting (93)–(95) along with (99) and (100) into (75) yields the semi-discrete set of coupled equations
Therein, \(\varvec{\tau }^{\mathrm {h}}\), \(\mathfrak {P}^{\mathrm {h}}\), \(\mathcal {\eta }^{\mathrm {h}}\) and \(\mathcal {H}_{\bullet }^{\mathrm {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). \(\mathcal {D}_{\mathrm {int}}^{\mathrm {h}}\) and \(\varvec{Q}^{\mathrm {h}}\) are semi-discrete definitions of the dissipation density and heat flux, cf. (40), (62) and (63). Moreover, the semi-discrete external contributions in (101)\(_1\) and (101)\(_2\) are formulated as
and
The coefficients of the matrices in (101)\(_3\) and (101)\(_4\) take the form
whereas the matrices in (101)\(_5\) are given by
Eventually, the semi-discrete functions \(\widehat{y}^{\mathrm {h}}\), \(\Phi ^{\mathrm {p,h}}\) and \(g_{\mathrm {c}}^{\mathrm {h}}\) denote the local hardening, the plastic yield and the critical fracture energy density, cf. (43), (53) and (47).
4 Temporal discretization
In a final step, the semi-discrete coupled problem (101) 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 \(\mathcal {T}\) into a sequence of times \(t_{0},\hdots ,t_{n},t_{n+1},\hdots ,T\), where \((\bullet )_{n}\) and \((\bullet )_{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 \(\{\varvec{q}_{A,n},\theta _{A,n},\alpha _{i,n}, r^{\mathrm {p}}_{i,n},{\mathfrak {s}}_{A,n},{\mathfrak {s}}_{\mathrm {L},A,n},{\mathfrak {s}}_{\mathrm {M},A,n}\}\) and the local plastic deformation variable \(\varvec{F}_{n}^{\mathrm {p,h}}\) at time \(t_{n }\) are known and the time step size \(\Delta 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 (101) defined as
Therein, a full-discrete definition of the internal dissipation is given by
Using small values for the plastic viscosity parameter \(\eta _{\mathrm {p}}\), we obtain
such that the internal dissipation can be recast as
for practical reasons.
In-plane bending test. Problem setting. The lines illustrate the fiber structure
In-plane bending test. Strain energy distribution of the Kirchhoff–Love shell formulation (left) and higher-gradient continuum formulation (right)
In-plane bending test. Problem setting. The lines illustrate the fiber structure
To solve the above multifield problem, we apply a staggered solution scheme, i.e. the displacement field along with the plastic and hardening fields \(\{\varvec{q}_{A,n+1}, \alpha _{i,n+1}, r^{\mathrm {p}}_{i,n+1}, \varvec{F}_{n+1}^{\mathrm {p,h}}\}\), the crack phase-fields \({\mathfrak {s}}_{\bullet ,A,n+1}\) and the temperature field \(\theta _{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 asFootnote 4
assuming that no further plastic deformation occurs within the time step. Based on this trial state, we evaluate the yield criteria (53). If \(\Phi ^{\mathrm {p}}_{\mathrm {tr}}\le 0\), then the process is purely elastic and the elastic trial state is the solution. If on the other hand \(\Phi ^{\mathrm {p}}_{\mathrm {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
where \(\sigma _{\mathrm {mat},a,n+1}=(\tau _{\mathrm {mat},a,n+1}^{\mathrm {dev}}+\tau _{\mathrm {mat},a,n+1}^{\mathrm {vol}})/J_{n+1}\). Note that in contrast to standard von Mises plasticity \(\varvec{n}_{n+1} \ne \varvec{n}_{\mathrm {tr}}\) and \(\Vert \varvec{n}_{n+1}\Vert \ne 1\), i.e. the plastic correction has to be performed by the Lagrange multiplier \(\lambda ^{\mathrm {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. [23].
5 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 \(\alpha _1=2\) in (28).
5.1 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. [62], 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.
Four point bending test. Boundary conditions of the four point bending test
Four point bending test. Force-displacement curves for bending tests
Tensile test (unidirectional). Problem setting. The lines illustrate the fiber structure
Tensile test (unidirectional). Load deflection results for unidirectional fiber reinforcements with different orientations
Tensile test (unidirectional). Results of the fiber crack phase-field (first row), 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 Fig. 7
Tensile test (bidirectional). Problem setting. The lines illustrate the fiber structure
Tensile test (bidirectional). Load deflection results for bidirectional, orthotropic fiber reinforcements with different orientations
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 Fig. 10
Thermal investigations. Load deflection results for isothermal simulations at different temperatures of \(\theta = [253,\, 273,\, 293]\,\mathrm {K}\) and a fiber orientation of \(\vartheta =30^{\circ }\)
Thermal investigations. Results of isothermal simulations at different temperatures of \(\theta = [253,\, 273,\, 293]\,\mathrm {K}\) (each from left to right) and a fiber orientation of \(\vartheta =30^{\circ }\). Results are shown for the fiber crack phase-field (first block), the plastic strain field (second block) and crack phase-field of the matrix material (third block) at the last deformation state marked in Fig. 12
5.1.1 In-plane bending test
We consider a benchmark example from dell’Isola et al. [62], where the left edge of a Kirchhoff-Love shell is clamped while the right edge is subjected to an external in-plane torque, chosen to match a reference analytical solution. To verify the proposed formulation in terms of in-plane bending stiffness parameterization, we take the shell deformation result obtained in [62], extrude it to the corresponding 3D geometry and calculate the energy. The plate is of size \(L \times W \times H =10\,\mathrm {mm} \times 1\,\mathrm {mm} \times 0.5\,\mathrm {mm}\) and is discretized by \(8 \times 2\times 1\) quadratic B-spline based elements, see Figure 1. Furthermore, we assume that the plate consists of a single fiber bundle with a cross section of \(A=HW=0.5\,\mathrm {mm}^{2}\) and a tensile stiffness of \(E_{\mathrm {fib}}=79,000\,\mathrm {N}/\mathrm {mm}^2\). The area moments of inertia of the fiber bundle with respect to the \(\varvec{e}_3\)-axis and the \(\varvec{e}_2\)-axis are given by \(I_{\varvec{e}_3}=HW^{3}/12=0.0417\,\mathrm {mm}^{4}\) and \(I_{\varvec{e}_2}=WH^{3}/12=0.0104 \,\mathrm {mm}^{4}\), respectively. Using these quantities we calibrate the bending stiffness parameters as \(c_{\#}=E_{\mathrm {fib}}I_{\varvec{e}_3}/A=6583.3333\,\mathrm {N}\) and \(c_{\perp }=E_{\mathrm {fib}}I_{\varvec{e}_2}/A=2212\,\mathrm {N}\).
In Fig. 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 in-plane bending stiffness parameter \(c_{\#}\). Note that the parameter \(c_{\perp }\) does not contribute to the simulation result, but will be investigated within the next example.
5.1.2 Four point bending test
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 \times W \times H =125\,\mathrm {mm} \times 25\,\mathrm {mm} \times 0.5\,\mathrm {mm}\) discretized by \(50 \times 10\times 2\) quadratic B-spline based elements. The bidirectional composite material has a matrix volume ratio of \(\zeta =0.53\) and the fibers are aligned in the \(\vartheta =0^{\circ }\) configuration, see Fig. 3 for the details on the fiber orientation. The four point bending test as shown in Fig. 4 leads to a pure out-of-plane 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 \(\mu =1630.4\,\mathrm {N/mm}^{2} \) and \(\alpha _1=2\) for the deviatoric part and \(\kappa = 6250\,\mathrm {N/mm}^{2}\) and \(\beta =-2\) for the volumetric part, which corresponds to a Young’s modulus of \(E_{{\text {mat}}}=4500 \, \mathrm {N/mm}^{2}\) and a Poisson’s ratio of \(\nu = 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_{{\text {fib}}}=79000\,\mathrm {N/mm}^{2}\) and the bending stiffness to \(c_{\perp }=0\,\mathrm {N}\). Secondly, we set the tensile stiffness of the fibers to zero and adjust the bending stiffness as \(c_{\perp }=E_{{\text {fib}}}H^2/12=1645.83\,\mathrm {N}\).
The applied bending stiffness of the continuum fiber model correlates to the out-of-plane 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 must be removed.
5.2 Tension Test
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\times W\times H =125\,\mathrm {mm} \times 25\,\mathrm {mm} \times 2\,\mathrm {mm}\). Figures 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\,\mathrm {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\,\mathrm {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_{\mathrm {fib}}=0.0025\,\mathrm {mm}^{2}\) and obtain a bending stiffness of \(c_{\perp }= c_{\#} = E_{{\text {fib}}} A_{\mathrm {fib}}/12 = 16.46\,\mathrm {N}\).
5.2.1 Unidirectional fiber reinforcement
We first analyze the behavior of a unidirectional reinforced composite material, with fiber orientations of \(\vartheta =[0^{\circ },10^{\circ },20^{\circ },30^{\circ },40^{\circ },65^{\circ },90^{\circ }]\), see Fig. 6. The load deflection results for isothermal simulations at \(\theta =293\,\mathrm {K}\) are shown in Fig. 7. Therein, crack initialization and final rupture of the fiber material are indicated by \(\square \) and \(\circ \), respectively. In addition, crack initialization and final rupture of the matrix material are indicated by \(\diamond \) and \(\varvec{\times }\), respectively. In Fig. 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 \(\vartheta =0^{\circ }\), 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^{\circ }\) 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 \(\vartheta =20^{\circ }\), 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 \(\vartheta =30^{\circ }\), 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 \(\vartheta =[40^{\circ },65^{\circ },90^{\circ }]\), 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 \(\vartheta =90^{\circ }\) 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\,\mathrm {mm}\).
5.2.2 Bidirectional fiber reinforcement
Next, we investigate the same tension test using a bidirectional reinforced material with fiber orientations of \(\vartheta =[0^{\circ },10^{\circ },20^{\circ },30^{\circ },45^{\circ }]\) as shown in Fig. 9. The load deflection results for isothermal simulations at \(\theta =293\,\mathrm {K}\) are depicted in Fig. 10. Again, crack initialization and final rupture of the fiber material are indicated by \(\square \) and \(\circ \), respectively. The final rupture of the matrix material is indicated by \(\varvec{\times }\). Crack phase-field results of the fiber and matrix material as well as results of the plastic strain are depicted in Fig. 11. Note that only phase-field results of the fiber aligned in \(\vartheta \)-direction is plotted.
For fiber orientations of \(\vartheta =[0^{\circ },10^{\circ },20^{\circ }]\), 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 stiffer 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 \(\vartheta =[30^{\circ },45^{\circ }]\), this stiffening effect becomes more pronounced as can be observed in Fig. 10. As already discussed, the additional, orthogonal oriented 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 ➈–➌.
5.2.3 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 Fig. 6 and apply a fiber orientation of \(\vartheta =30^{\circ }\).
Figure 12 shows the load deflection result for isothermal simulations using temperatures of \(\theta = [253,\, 273,\, 293]\,\mathrm {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 Fig. 13. As already observed previously, for \(\theta =293\,\mathrm {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 \(\theta =[ 253,\,273]\,\mathrm {K}\) the matrix material fails before any fiber cracks occur.
6 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 second-gradient theory, a novel hybrid phase-field model and a temperature dependent GTN-type 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.
Notes
The \(\bullet \) indicates the matrix material or the respective fiber direction.
The Macaulay brackets are defined by \(\langle x\rangle =(x+|x|)/2\).
In index notation the internal virtual work yields
$$\begin{aligned} \delta W^{\mathrm {e,int}}= & {} \int \limits _{\mathcal {B}_0} -(\tau _{ij}F^{\mathrm {-T}}_{jK})_{,K}\delta \varphi _{i} + (\tau _{ij}F^{\mathrm {-T}}_{jK}\delta \varphi _{i})_{,K} \nonumber \\&-\,\mathfrak {P}_{iJK,K}\delta \varphi _{i,J} + (\mathfrak {P}_{iJK}\delta \varphi _{i,J})_{,K} \,\mathrm {d}V \end{aligned}$$(67)For the sake of readability, we neglect the labeling of the spatial approximation in the following.
References
Aldakheel F (2016) Mechanics of nonlocal dissipative solids: gradient plasticity and phase field modeling of ductile fracture. PhD thesis, University of Stuttgart
Aldakheel F, Hudobivnik B, Hussein A, Wriggers P (2018) Phase-field modeling of brittle fracture using an efficient virtual element scheme. Comput Methods Appl Mech Eng 341:443–466
Aldakheel F, Mauthe S, Miehe C (2014) Towards phase field modeling of ductile fracture in gradient-extended elastic-plastic solids. Proc Appl Math Mech 14:411–412
Aldakheel F, Wriggers P, Miehe C (2018) A modified Gurson-type plasticity model at finite strains: formulation, numerical analysis and phase-field coupling. Comput Mech. 62:815–833
Alessi R, Ambati M, Gerasimov T, Vidoli S, De Lorenzis L (2018) Comparison of phase-field models of fracture coupled with plasticity. Advances in computational plasticity, vol 46. Springer, pp 1–21
Amor H, Marigo JJ, Maurini C (2009) Regularized formulation of the variational brittle fracture with unilateral contact: numerical experiments. J Mech Phys Solids 57:1209–1229
Asmanoglo T, Menzel A (2017) A finite deformation continuum modelling framework for curvature effects in fibre-reinforced nanocomposites. J Mech Phys Solids 107:411–432
Asmanoglo T, Menzel A (2017) A multi-field finite element approach for the modelling of fibre-reinforced composites with fibre-bending stiffness. Comput Methods Appl Mech Eng 317:1037–1067
Bai Y, Wierzbicki T (2008) A new model of metal plasticity and fracture with pressure and Lode dependence. Int J Plast 24:1071–1096
Besson J (2010) Continuum models of ductile fracture: a review. Int J Damage Mech 19:3–52
Borden MJ, Hughes TJR, Landis CM, Verhoosel CV (2014) A higher-order phase-field model for brittle fracture: Formulation and analysis within the isogeometric analysis framework. Comput Methods Appl Mech Eng 273:100–118
Bornemann PB, Cirak F (2013) A subdivision-based implementation of the hierarchical b-spline finite element method. Comput Methods Appl Mech Eng 253:584–598
Brünig M, Brenner D, Gerke S (2015) Stress state dependence of ductile damage and fracture behavior: experiments and numerical simulations. Eng Fract Mech 141:152–169
Brünig M, Gerke S (2011) Simulation of damage evolution in ductile metals undergoing dynamic loading conditions. Int J Plast 27:1598–1617
Bryant EC, Sun W (2018) A mixed-mode phase field fracture model in anisotropic rocks with consistent kinematics. Comput Methods Appl Mech Eng 342:561–584
Cayzac H-A, Saï K, Laiarinandrasana L (2013) Damage based constitutive relationships in semi-crystalline polymer by using multi-mechanisms model. Int J Plast 51:47–64
Cottrell JA, Hughes TJR, Bazilevs Y (2009) Isogeometric analysis: toward integration of CAD and FEA. Wiley
Da Costa Mattos HS, Brandao JF, Amorim F, De Souza Araujo PV, Reis J (2019) A unified expression to estimate the stress-strain curve of polyamides at different temperatures. Mater Res Express 6:015304
de Borst R, Pamin J, Geers MGD (1999) On coupled gradient-depended plasticity and damage theories with a view to localization analysis. Eur J Mech A/Solids 18:939–962
dell’Isola F, Giorgio I, Pawlikowski M, Rizzi NL (2016) Large deformations of planar extensible beams and pantographic lattices: heuristic homogenization, experimental and numerical examples of equilibrium. Proc R Soc Lond A: Math, Phys Eng Sci 472(2185):20150790
dell’Isola F, Seppecher P, Spagnuolo M, Barchiesi E, Hild F, Lekszycki T, Giorgio I, Placidi L, Andreaus U, Cuomo M, Eugster SR, Pfaff A, Hoschke K, Langkemper R, Turco E, Sarikaya R, Misra A, De Angelo M, D’Annibale F, Bouterf A, Pinelli X, Misra A, Desmorat B, Pawlikowski M, Dupuy C, Scerrato D, Peyre P, Laudato M, Manzari L, Göransson P, Hesch C, Hesch S, Franciosi P, Dirrenberger J, Maurin F, Vangelatos Z, Grigoropoulos C, Melissinaki V, Farsari M, Muller W, Abali BE, Liebold C, Ganzosch G, Harrison P, Drobnicki R, Igumnov L, Alzahrani F, Hayat T (2019) Advances in pantographic structures: design, manufacturing, models, experiments and image analyses. Cont Mech Thermodyn 31:1231–1282
Dittmann M (2017) Isogeometric analysis and hierarchical refinement for multi-field contact problems. PhD thesis, University of Siegen
Dittmann M, Aldakheel F, Schulte J, Schmidt F, Krüger M, Wriggers P, Hesch C (2020) Phase-field modeling of porous-ductile fracture in non-linear thermo-elasto-plastic solids. Comput Methods Appl Mech Eng 361:112730
Dittmann M, Aldakheel F, Schulte J, Wriggers P, Hesch C (2018) Variational phase-field formulation of non-linear ductile fracture. Comput Methods Appl Mech Eng 342:71–94
Dittmann M, Hesch C, Schulte J, Aldakheel F, Franke M (2017) Multi-field modelling and simulation of large deformation ductile fracture. In: Proceedings of the XIV international conference on computational plasticity. Fundamentals and applications, pp 556–567
Dittmann M, Krüger M, Schmidt F, Schuß S, Hesch C (2019) Variational modeling of thermomechanical fracture and anisotropic frictional mortar contact problems with adhesion. Comput Mech 63(3):571–591
Duong TX, Khiem VN, Itskov M, Sauer RA (2020) A general theory for anisotropic Kirchhoff-Love shells with embedded fibers and in-plane bending. Mater Sci (submitted). arXiv:2101.03122
Enakoutsa K, Leblond JB, Perrin G (2006) Numerical implementation and assessment of a phenomenological nonlocal model of ductile rupture. Comput Methods Appl Mech Eng 196:1946–1957
Eringen AC (1999) Microcontinuum field theories I: foundations and solids. Springer
Germain P (1973) The method of virtual power in continuum mechanics. Part 2: microstructure. SIAM J Appl Math 25:556–575
Gurland J, Plateau J (1963) The mechanism of ductile rupture of metals containing inclusions. Trans ASM 56:443–454
Gurson AL (1975) Plastic flow and fracture behavior of ductile materials incorporating void nucleation, growth and coalescence. PhD thesis, Division of Engineering, Brown University
Gurson AL (1977) Continuum theory of ductile rupture by void nucleation and growth, part I—yield criteria and flow rules for porous ductile media. J Eng Mater Technol 99:2–15
Heider Y, Markert B (2017) A phase-field modeling approach of hydraulic fracture in saturated porous media. Mech Res Commun 80:38–46
Hesch C, Gil AJ, Arranz Carreno A, Bonet J, Betsch P (2014) A Mortar approach for Fluid-Structure Interaction problems: immersed strategies for deformable and rigid bodies. Comput Methods Appl Mech Eng 278:853–882
Hesch C, Schuß S, Dittmann M, Franke M, Weinberg K (2016) Isogeometric analysis and hierarchical refinement for higher-order phase-field models. Comput Methods Appl Mech Eng 303:185–207
Hesch C, Weinberg K (2014) Thermodynamically consistent algorithms for a finite-deformation phase-field approach to fracture. Int J Numer Methods Eng 99(12):906–924
Huespe A, Needleman A, Oliver J, Sánchez, PJ (2012) A finite strain, finite band method for modeling ductile fracture. Int J Plast 28:53–69
Huetter G, Linse T, Mühlich U, Kuna M (2013) Simulation of ductile crack initiation and propagation by means of a non-local Gurson-model. Int J Solids Struct 50:662–671
Javili A, dell’Isola F, Steinmann P (2013) Geometrically nonlinear higher-gradient elasticity with energetic boundaries. J Mech Phys Solids 61(12):2381–2401
Krüger M, Dittmann M, Aldakheel F, Härtel A, Wriggers P, Hesch C (2020) Porous-ductile fracture in thermo-elasto-plastic solids with contact applications. Comput Mech 65:941–966
Kuhn C, Schlüter A, Müller R (2015) On degradation functions in phase field fracture models. Comput Mater Sci 108:374–384
Larsson R, Razanica S, Josefson BL (2015) Mesh objective continuum damage models for ductile fracture. Int J Numer Methods Eng 106(10):840–860
Leblond JB, Perrin G, Devaus J (1995) An improved Gurson-type model for hardenable ductile metals. Eur J Mech—A/Solids 14:499–527
Lehmann T, Blix U (1985) On the coupled thermo-mechanical process in the necking problem. Int J Plast 1:175–188
Lemaitre J (1985) A continuous damage mechanics model for ductile fracture. J Eng Mater Technol 107:83–89
Lemaitre J (1992) A course on damage mechanics. Springer
Lemaitre J, Chaboche JL (1990) Mechanics of solid materials. Cambridge University Press
Li H, Fu MW, Lu J, Yang H (2011) Ductile fracture: experiments and computations. Int J Plast 27:147–180
Miehe C, Aldakheel F, Teichtmeister S (2017) Phase-field modeling of ductile fracture at finite strains: a robust variational-based numerical implementation of a gradient-extended theory by micromorphic regularization. Int J Numer Methods Eng 111(9):816–863
Miehe C, Kienle D, Aldakheel F, Teichtmeister S (2016) Phase field modeling of fracture in porous plasticity: a variational gradient-extended Eulerian framework for the macroscopic analysis of ductile failure. Comput Methods Appl Mech Eng 312:3–50
Miehe C, Schänzel L, Ulmer H (2015) Phase field modeling of fracture in multi-physics problems. Part I. Balance of crack surface and failure criteria for brittle crack propagation in thermo-elastic solids. Comput Methods Appl Mech Eng 294:449–485
Miehe C, Welschinger F, Hofacker M (2010) Thermodynamically consistent phase-field models of fracture: variational principles and multi-field FE implementations. Int J Numer Methods Eng 83(10):1273–1311
Mindlin RA (1964) Micro-structure in linear elasticity. Arch Ration Mech Anal 16:51–78
Mindlin RA (1965) On the equations of elastic materials with microstructure. Int J Solids Struct 1:73–78
Nahshon K, Hutchinson JW (2008) Modification of the Gurson model for shear failure. Eur J Mech A/Solids 27:1–17
Needleman A, Tvergaard V (1984) An analysis of ductile rupture in notched bars. J Mech Phys Solids 32:461–490
Paggi M, Reinoso J (2017) Revisiting the problem of a crack impinging on an interface: a modeling framework for the interaction between the phase field approach for brittle fracture and the interface cohesive zone model. Comput Methods Appl Mech Eng 321:145–172
Reis JML, de Santanna T, Coutinho PP, Monteiro AH, Teixeira SCS, Chaves EG, da Costa Mattos HS (2016) Coupled effect of ageing and temperature in the mechanical behaviour of a polyamide. Polym Test 53:264–275
Reusch F, Svendsen B, Klingbeil D (2003) A non-local extension of Gurson-based ductile damage modeling. Comput Mater Sci 26:219–229
Reusch F, Svendsen B, Klingbeil D (2003) Local and non-local Gurson-based ductile damage and failure modelling at large deformation. Eur J Mech A/Solids 22:779–792
Schulte J, Dittmann M, Eugster SR, Hesch S, dell’Isola F, Hesch C (2020) Isogeometric shell elements for fiber reinforced composites. Comput Methods Appl Mech Eng 362:112845
Schulte J, Krüger M, Dittmann M, Hesch C (2018) Multi-field modeling of thermomechanical coupled fracture problems. Proc Appl Math Mech 18(1):e201800348
Schuß S, Dittmann M, Klinkel S, Wohlmuth B, Hesch C (2019) Multi-patch isogeometric analysis for Kirchhoff-Love shell elements. Comput Methods Appl Mech Eng 349:91–116
Seabra MRR, Andrade JMA, Cesar de Sa FXC, Pires FMA (2011) Continuous-discontinuous formulation for ductile fracture. Int J Mater Form 4:271–281
Selles N, King A, Proudhon H, Saintier N, Laiarinandrasana L (2018) Time dependent voiding mechanisms in polyamide 6 submitted to high stress triaxiality: experimental characterisation and finite element modelling. Mech Time-Depend Mater 22:351–371
Simo JC, Miehe C (1992) Associative coupled thermoplasticity at finite strains: formulation, numerical analysis and implementation. Comput Methods Appl Mech Eng 98(1):41–104
Soldatos KP (2010) Second-gradient plane deformations of ideal fibre-reinforced materials: implications of hyper-elasticity theory. J Eng Math 68:99–127
Spencer AJM, Soldatos KP (2007) Finite deformations of fibre-reinforced elastic solids with fibre bending stiffness. Int J Non-linear Mech 42:355–368
Steigmann DJ (2012) Theory of elastic solids reinforced with fibers resistant to extension, flexure and twist. Int J Non-linear Mech 47:734–742
Steigmann DJ (2018) Equilibrium of elastic lattice shells. J Eng Math 109:47–61
Steigmann DJ, dell’Isola F (2015) Mechanical response of fabric sheets to three-dimensional bending, twisting, and stretching. Acta Mech Sin 31(3):373–382
Steigmann DJ, Pipkin AC (1991) Equilibrium of elastic nets. Philos Trans R Soc Lond 335:419–454
Steinmann P, Miehe C, Stein E (1994) Comparison of different finite deformation inelastic damage models within multiplicative elastoplasticity for ductile materials. Comput Mech 13:458–474
Teichtmeister S, Kienle D, Aldakheel F, Keip M-A (2017) Phase field modeling of fracture in anisotropic brittle solids. Int J Non-linear Mech 97:1–21
Toupin RA (1964) Theories of elasticity with couple stress. Arch Ration Mech Anal 17:85–112
Toupin RA (1962) Elastic materials with couple-stresses. Arch Ration Mech Anal 11:385–414
Tvergaard V (1982) On localization in ductile materials containing spherical voids. Int J Fract 18:237–252
Tvergaard V, Needleman A (1984) Analysis of the cup-cone fracture in a round tensile bar. Acta Metall Mater 32:157–169
Verhoosel CV, de Borst R (2013) A phase-field model for cohesive fracture. Int J Numer Methods Eng 96:43–62
Xue Z, Pontin MG, Zok FW, Hutchinson JW (2010) Calibration procedures for a computational model of ductile fracture. Eng Fract Mech 77:492–509
Zdebel U, Lehmann Th (1987) Some theoretical considerations and experimental investigations on a constitutive law in thermoplasticity. Int J Plast 3:369–389
Zhang X, Vignes C, Sloan SW, Sheng D (2017) Numerical evaluation of the phase-field model for brittle fracture with emphasis on the length scale. Comput Mech 59(5):737–752
Acknowledgements
Support for the present research was provided by the Deutsche Forschungsgemeinschaft (DFG) under grant HE5943/8-1 and DI2306/1-1. The authors C. Hesch and M. Dittmann gratefully acknowledge this support.
Funding
Open Access funding enabled and organized by Projekt DEAL.
Author information
Authors and Affiliations
Corresponding authors
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Dittmann, M., Schulte, J., Schmidt, F. et al. A strain-gradient formulation for fiber reinforced polymers: hybrid phase-field model for porous-ductile fracture. Comput Mech 67, 1747–1768 (2021). https://doi.org/10.1007/s00466-021-02018-0
Received:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s00466-021-02018-0