A thermodynamically consistent elastoviscoplastic phase-field framework for structural damage in PTFE

Abstract

Deformation in polymers is highly dependent on molecular structures and motion and relaxation mechanisms, which are highly influenced by temperature and mechanical load history. These features imply that some models can fit for specific classes of polymers and not for others; moreover, these models also include several non-linearities, which turns out to be challenging for computational simulation. This work develops and simulates a thermal-structural phase-field model for the polytetrafluorethylene (PTFE) polymer. The constitutive multimechanism model used considers a non-isothermal non-linear elastoviscoplastic model, able to represent elastic molecular interactions, and viscoplastic flow from polymer segments. Material parameters for complex rheological models are addressed, through a genetic algorithm, to adjust curves from simulated models to stress-strain experiments available in literature. Results of stress-strain curves, followed by rupture, for a temperature ranging from -50° C to 150° C are plotted in comparison with experimental results, presenting a reasonable fit.

Keywords
Phase-field; Finite Element Method; Polymer Mechanics

Graphical Abstract

1 INTRODUCTION

Nowadays, polymeric materials have been used largely in engineering applications. Due to advances in reinforcement fibres and combined molecular structures, polymeric materials can offer many advantages when compared to others materials even in most precise or under extremely severe usage conditions.

Considering the significant increasing of polymeric products and components in many applications taking into account different safety and quality requirements, the appropriate description of the mechanical behaviour of polymers and their structures is necessary for design purposes. Usually, polymers have quite different structural and thermal responses when compared to metals, which can impact their manufacturing process and also usage conditions. For this reason, many material constitutive models were developed in last decades and still represent a growing field.

Deformation in polymers are highly dependent on their molecular structures. Motion and relaxation mechanisms are highly influenced by temperature and mechanical load history. These features imply that some models can fit well for specific classes of polymers and not for others. There are still several non-linearities, which turn out to be challenging for computational simulation.

Most common failure modes in mechanical components include fracture and, consequently, are very relevant in terms of simulation objectives, in order to avoid damage to structures and hazard to people. Due to this fact and new material models, failure simulation is applied in very early stages of product development to evaluate risks and, consequently, ensure safety to products. In this context, phase-field models have emerged as an efficient and promising technique for evaluating fracture evolution, due to continuous treatment of discontinuities and sharp surfaces that challenge simulations.

This paper aims to develop and simulate a thermal-structural phase-field model for the polytetrafluorethylene (PTFE) semi-crystalline polymer taking into account relevant constitutive models in the literature. The response of this model will be compared to experiments in order to evaluate its behaviour, explore possibilities for internal variables and, finally, contribute to the current state of the art for this material.

2 LITERATURE REVIEW

For describing properly the behavior of polymers, a large number of works, started decades ago, aimed to develop constitutive models in order to reproduce different aspects of mechanical response. When properly adjusted, these models can, with reasonable accuracy, reproduce the macroscopic response of polymers under multi-axial loading conditions over a wide range of strain rates and temperatures.

A starting point for many others authors describing amorphous polymers response, Haward and Thackray (1968) presented an one-dimensional rheological model, composed by two mechanisms (Eyring dashpot and Langevin spring). Stress response is built from different contributions, both from intermolecular interactions and also from preferred orientation of the entangled molecular network. The model reproduced features of both glass-like and non-linear rubber-like behaviors, for materials like cellulose acetate, cellulose nitrate and PVC.

Subsequently, and counting on development of rubber elasticity by Treloar (1975), many others authors extended constitutive models for including more dimensions and complex behaviors. Parks, Argon and Bagepalli (1984) implemented a three-dimensional constitutive model for PMMA, taking into account in the rheological model the macromolecular structure and micro-mechanisms of plastic flow as well. Boyce, Parks and Argon (1988) included the effect of deformation rate, pressure, true strain softening and temperature on the plastic resistance. The result of this work is the well known BPA model. Wu and Giessen (1993), based on constitutive description of rubber elasticity, described an isothermal three-dimensional model for elastomers and also for polycarbonate, taking into account large strain of molecules through a non-Gaussian statistical mechanics. Tervoort et al. (1997) introduced a compressible-Leonov model, in which the elastic volume response is separated from the elasto-viscoplastic isochoric deformation and making possible to extend the model in a straightforward way to include a spectrum of relaxation times. In Govaert, Timmermans and Brekelmans (2000), the compressible-Leonov model is extended to describe the phenomenological aspects of the large strain mechanical behavior of glassy polymers and applied to polycarbonate response simulation. Anand and Gurtin (2003), in the sequence of including complex mechanisms of polymers, modelled a three-dimensional continuum isothermal theory for elasto-viscoplastic deformation of amorphous solids such as polymeric and metallic glasses. This model was able to capture highly non-linear stress-strain behaviour that precedes the yield-peak and gives rise to post-yield strain softening. Buckley et al. (2004) showed constitutive responses of three glassy thermoset polymers at high and low strain rates.

All these models mentioned above were used to describe, primarily, isothermal deformation of polymers below their glass transition temperature. To include temperature evolution, Anand et al. (2009) presented a detailed continuum mechanic development of a thermo-mechanically coupled elasto-viscoplastic theory to model the strain rate and temperature dependent large deformation response of amorphous polymeric materials. This work was applied to large-strain compression experiments on three amorphous polymeric materials: polymethylmethacrylate (PMMA), polycarbonate (PC), and a cyclo-olefin polymer (Zeonex-690R). Temperature range spanning room temperature to slightly below the glass transition temperature of each material, and strain rates ranging from 1×10⁻⁴ s⁻¹ to 1×10⁻¹ s⁻¹ were considered in the work of Ames et al. (2009), and also for temperature range spanning room temperature to approximately 50°C over the glass transition temperature are showed by Srivastava et al. (2010).

In addition to amorphous polymers, multi-mechanism models were also used for simulating semi-crystallyne polymers and, in special, fluoropolymers. In this area, several models are proposed, like the Dual Network Fluoropolymer Model (DNF) and the Three Network Model (TNF), as described by Bergstrom (2015), in an evolution of models described earlier by Bergström and Boyce (1998), Arruda, Boyce and Jayachandran (1995) and Kletschkowski, Schomburg and Bertram (2002). In many cases, due to complexity of constitutive models, composed by several branches and internal variables that cannot be obtained experimentally, additional methods are used in order to characterize the material. One of the most used method is based on genetic algorithms, as presented by Fernández et al. (2018) for fitting parameters in elastomeric material model, and by Colak and Cakir (2019) in modelling strain and temperature dependent behavior of epoxy resin.

Many methods for evaluating damage mechanisms are based on Griffith (1921) theory, which introduced the concept of energy release rate as the responsible for the onset of crack.

This approach, with limitations in predicting crack paths and new crack nucleations, was later improved by others researchers in order to understand damage evolution, fatigue and fracture in materials, such as Lemaitre and Desmorat (2006), Kanninen, McEvily and Popelar (1986) and Anderson (2017).

An important variant of Griffith method was proposed by Francfort and Marigo (1998), in which the crack initiation, propagation and ramification process is driven by minimization of an energy functional, in terms of displacement field and crack surface. Considering that crack surface is unknown in the beginning of the process, it is replaced by a continuous scalar field, named phase-field, which interpolates smoothly between undamaged and broken material.

Due to these new techniques, more attention has been given to the phase-field methodology as a promising way to deal with modelling of material damage. Originally applied to problems related to separation of fluids by Cahn and Hilliard (1958), these models have been also used in a great variety of multiphase problems such as in the study of tumor growth by Silva (2009), vesicle-fluid iteration by Du, Li and Liu (2007), solidification by Allen and Cahn (1972) and fluid-structure interation problems by Sun, Xu and Zhang (2014). In addition to these categories, many fracture and crack initiation problems are solved using the Ginsburg-Landau free-energy functional in order to predict the phase-field evolution in continuum mechanics and solid materials, as described by Fabrizio, Giorgi and Morro (2006) and Amendola, Fabrizio and Golden (2016). In the sequence of these advances in phase-field for fracture and fatigue, Boldrini et al. (2016) presented a general thermodynamically consistent, non-isothermal, nonlocal framework under the hypothesis of small deformation.

Concerning papers considering aspects more closely related to the ones of the present work, we mention the following. SCHÄNZEL (2015) and Miehe, Schaenzel and Ulmer (2015) published a phase-field model for evaluating crazing-induced brittle fracture of glassy polymers. The use of driving force for damage evolution, however, is not thermodynamically consistent. In order to improve this aspect, Narayan and Anand (2021) developed a new gradient-damage theory, able to describe deformation and failure of amorphous polymers in thermodynamic consistent way. However, Narayan and Anand (2021) did not consider the important infuence of the temperature on the polymer behavior. Srivastava et al. (2010) considered this infuence of temperature, but did include the possibility of material damage.

To advance the knowledge about the infuence of temperature and damage on polymer behavior, in this work we present a new phase-field model by extending the work of Srivastava et al. (2010) and including in a thermodynamically consistent way the possibility of material damage. This done by using a framework similar to the one of Boldrini et al. (2016). Moreover, the results of numerical simulations with our model were compared to experiments with the polytetrafluorethylene (PTFE) polymer in a rather large range of temperatures with good results. We think that this shows the relevance of the presented model.

3 PTFE STRUCTURE

Polymers containing fluorine present advantages related to improved mechanical, electrical, chemical and friction wear resistance. According to Brown et al. (2004), the most used fluorocarbon polymer in applications is the Polytetrafluorethylene (PTFE) commonly known under the trade name Teflon, which has a complex microstructure with different crystalline phases in terms of temperature at ambient pressure. These phases influence the mechanical behavior, including the resistance to fracture.

PTFE is a semi-crystalline polymer with amorphous and crystalline phases whose proportions depend on the parameters of the production process, such as cooling rate of the sintering temperature.

Amorphous phase, according to Calleja et al. (2013), is composed by two distinct regions. The first one is called mobile amorphous fraction (MAF), whose relaxation process occurs at low temperatures (T = -103°C). The other one is at the interface between crystalline and MAF domains and called rigid amorphous fraction (RAF). Its relaxation process occurs at highers temperatures (T = 116°C) and strongly depends on the degree of crystallinity.

As discussed by Brown and Dattelbaum (2005), the crystalline phase of PTFE depends on temperature and pressure. Experimentally, it is observed four different phases. Three of them (I, II and IV) exist at ambient pressure with transition temperatures at 19°C and 30°C, which makes predictions on the PTFE behaviour in this range more important as they are closer to the ambient temperature. In terms of microstructure, transition at 19°C (from phase II to phase IV) represents an unraveling of helical conformation, given by a well ordered triclinic structure (13 atoms/180 turn) to a partial ordered hexagonal structure (15 atoms/turn). After 30°C, more disorder is included in crystalline structures and phase I has a pseudohexagonal arrange.

4 DERIVATION OF THE GENERAL BALANCE EQUATIONS

The main purpose of this section is to derive a one-dimensional non-local damage and fracture phase-field model for polymers, adapted from Boldrini et al. (2016), based on continuum mechanics and thermodynamics balance laws. In this context, a review of basic concepts of kinematics and stress-strain theories under large displacement is first considered. The multiplicative decomposition of the deformation gradient $F$ is used, as in Anand and Ames (2006).

4.1 Kinematics

Consider a bar of length $L_0$ in a given one-dimensional reference configuration at time $t=0$, whose arbitrary material point is denoted by $X$ and boundaries, denoted by $\Gamma$, given by the two points at the bar ends, with normals with values -1 and 1. The motion of the bar, denoted by $\chi \left(X,t\right)$, maps the undeformed to the current configuration with spatial points indicated by $x$ and length $L$, as illustrated in Figure 1.

Figure 1
Bar in reference configuration at t = 0 and arbitrary material point X mapped to the current (deformed) configuration at point x.

The motion of each particle is described by a smooth one-to-one mapping $x=\mathbf{\chi }\left(X,t\right)$, in which $\mathbf{\chi }\left(X,t\right)=X+u\left(X,t\right)$, where $u$ is the displacement field. The deformation gradient, which maps the material infinitesimal fibre $dX$ at point $X$ to the spatial infinitesimal fibre $dx$ at point $x$, velocity and velocity gradient in spatial description are given, respectively, by

$\begin{array}{l}F\left(X,t\right)=\frac{\partial \chi }{\partial X},\\ v\left(x,t\right)=\dot{\mathbf{\chi }}\left({\chi }^{-1}\left(x,t\right),t\right),\\ L\left(x,t\right)=\frac{\partial v}{\partial x}=\dot{F}{F}^{-1}.\end{array}$ (1)

For an elasto-viscoplastic finite strain model, used in this work, the Kroner-Lee multiplicative decomposition (Lee (1969)) is considered in terms of the elastic $F_e$ and plastic $F_p$ parts of $F$ as

$F=F_e{F}_p.$ (2)

This hypothesis considers the complete separation between the elastic deformation due to stretch and the plastic deformation due to flow of defects, both in the microscopic structure.

Due to the microscopical complexity, a general and accurate polymer response including thermo-structural behaviour, requires multiple structural mechanisms. In this sense, the classical Kroner-Lee decomposition is replaced by the generalized Zener multi-mechanism model, with the following decomposition of the deformation gradient:

$F=F_e^{\left(\alpha \right)}{F}_p^{\left(\alpha \right)},$ (3)

where $\alpha$ indicates a mechanism.

Figure 2 illustrates the generalized Zener model with $M$ mechanisms. A more detailed material description including, for example, stress-relaxation and creep, requires larger number of mechanisms and, consequently, more experimental parameters must be obtained. For this reason, the lowest number of mechanisms $\alpha$ should be considered, since finding the correct material parameters become an important aspect of these phenomenological continuum-level models.

Figure 2
Generalized Zener model. Adapted from Anand and Ames (2006).

Each branch of the rheological model comprises an elastic part described by a spring in series with a set composed of a dissipative dashpot term (driven by plastic deformation $F_p^{\left(\alpha \right)}$ and internal variables ${\xi }^{\left(\alpha \right)}$) in parallel with an energy storage mechanism associated to plastic deformation (driven by $A^{\left(\alpha \right)}$).

Physically, $F_e^{\left(\alpha \right)}$ represents an elastic deformation, such as stretching and rotation of intermolecular bonds and polymer chains. $F_p^{\left(\alpha \right)}$ is a local inelastic deformation at $X$, due to mechanisms such as slippage of polymer molecules. This last transformation takes material coordinates to a coherent structure in the relaxed space $M_X$, being $M_X$ the range of $F_p^{\left(\alpha \right)}$, as shown in Figure 3.

Figure 3
Reference spaces including structural space given by the Kroner-Lee decomposition.

This relaxed space, also called sometimes intermediate or structural space, is defined as the output of the linear transformation $F_p^{\left(\alpha \right)}$ and the input of the linear transformation $F_e^{\left(\alpha \right)}$. Thus, $M_X$ is not related to a real physical space, but to a mathematical construction based on the Kroner-Lee transformation. Quantities in this space will be represented with a bar over the variable symbol.

Applying the Kroner-Lee decomposition for the velocity gradient, Eq. (1)(iii) is then rewritten as

$\begin{array}{c}L=L_e^{\left(\alpha \right)}+F_e^{\left(\alpha \right)}\overline{L}{\mathrm{}}_p^{\left(\alpha \right)}{\left(F_e^{\left(\alpha \right)}\right)}^{-1},\end{array}$ (4)

with

$\begin{array}{l}{L}_e^{\left(\alpha \right)}={\dot{F}}_e^{\left(\alpha \right)}{\left(F_e^{\left(\alpha \right)}\right)}^{-1},\\ {\overline{L}}_p^{\left(\alpha \right)}={\dot{F}}_p^{\left(\alpha \right)}{\left(F_p^{\left(\alpha \right)}\right)}^{-1},\end{array}$ (5)

where ${\overline{L}}_p^{\left(\alpha \right)}$ is the inelastic velocity gradient in the intermediate configuration and $F_e^{\left(\alpha \right)}{\overline{L}}_p^{\left(\alpha \right)}{\left(F_e^{\left(\alpha \right)}\right)}^{-1}$ its counterpart in the current configuration.

In an one-dimensional configuration, there is no rotation. Consequently, the decomposition of velocity gradient in a symmetric (stretching) and skew-symmetric (spin) parts is given by

$\begin{array}{l}{L}_e^{\left(\alpha \right)}=D_e^{\left(\alpha \right)},\\ {\overline{L}}_p^{\left(\alpha \right)}={\overline{D}}_p^{\left(\alpha \right)}.\end{array}$ (6)

Analogously, polar decompositions of $F_e^{\left(\alpha \right)}$ and $F_p^{\left(\alpha \right)}$ provide

$\begin{array}{c}{F}_e^{\left(\alpha \right)}=U_e^{\left(\alpha \right)}=V_e^{\left(\alpha \right)},\\ F_p^{\left(\alpha \right)}=U_p^{\left(\alpha \right)}=V_p^{\left(\alpha \right)},\end{array}$ (7)

with $U_e^{\left(\alpha \right)}$, $U_p^{\left(\alpha \right)}$, $V_e^{\left(\alpha \right)}$ and $V_p^{\left(\alpha \right)}$ being variables describing deformations.

The right and left Cauchy-Green deformations, for both elastic and plastic cases, are given, respectively, by

${\overline{C}}_e^{\left(\alpha \right)}={\left(F_e^{\left(\alpha \right)}\right)}^T{F}_e^{\left(\alpha \right)},\mathrm{}{B}_e^{\left(\alpha \right)}=F_e^{\left(\alpha \right)}{\left(F_e^{\left(\alpha \right)}\right)}^T,$ (8)

$C_p^{\left(\alpha \right)}={\left(F_p^{\left(\alpha \right)}\right)}^T{F}_p^{\left(\alpha \right)},\mathrm{}{\overline{B}}_p^{\left(\alpha \right)}=F_p^{\left(\alpha \right)}{\left(F_p^{\left(\alpha \right)}\right)}^T.$ (9)

We define also

$J\left(X\right)=\text{det}\left(F\left(X\right)\right),$ (10)

which for 1D problems is equal to $F\left(X\right)$ in the case of compressibility and equal to one in case of incompressibility.

4.2 Governing equations

Consider, initially, the spatial reference frame and the respective variables for obtaining the governing equations through the principle of virtual power (PVP).

Following the arguments given in Frémond and Nedjar (1996), there is in addition to the standard velocity $v$, a microscopic velocity ($c$) referring to the time rate of change of phase-field $\phi$ used here to describe the failure of the material. In this context, the phase-field is considered a dynamical variable with microscopic bulk forces and stresses. In this work, analogously to Boldrini (2018), it is assumed that microscopic acceleration forces are zero, which implies in a purely dissipative evolution for the structures described by $\phi$.

Applying the PVP, consider that the external power comes from the macroscopic stresses $t$ on the bar ends and the axial distributed loads $b\left(x\right)$ per unit of mass along the bar length. Moreover, it is considered a density of energy, $a$, supplied by exterior to evolving structures (such as phase-field describing damage in material) and a surface density of energy given by flux $t_h$. Virtual macroscopic and microscopic velocities are denoted by $\widehat{v}$ and $\widehat{c}$, respectively. Based on the previous definitions, the expression of the external virtual power is given by

${\mathcal{P}}_e\left(\widehat{v},\widehat{c}\right)={\int }_0^L\rho b\widehat{v}dx+t\left(0\right)\widehat{v}\left(0\right)+t\left(L\right)\widehat{v}\left(L\right)+{\int }_0^L\rho a\widehat{c}dx+t_h\left(0\right)\widehat{c}\left(L\right)+t_h\left(L\right)\widehat{c}\left(L\right).$ (11)

Considering the internal power, due to the existence of elastic and plastic stresses ($T_e^{\left(\alpha \right)}$ and $T_p^{\left(\alpha \right)}$) for each branch, the expression must be split in terms of virtual rates ${\widehat{L}}_e$ and ${\overline{\widehat{L}}}_p$. Analogously to the external power, microscopic terms $b_{\phi }$ and $h$ are considered, representing, respectively, the volumetric density of energy exchanged by variation of a unit of $\phi$ and the flux of energy associated to the spatial variation of a unit of $\phi$ in a unit of time. Based on previous definitions, the expression of the internal virtual power for a material with $M$ mechanisms is given by

$\begin{array}{c}{\mathcal{P}}_i\left({\overline{\widehat{L}}}_p^{\left(\alpha \right)},{\widehat{L}}^{e\left(\alpha \right)},\widehat{c}\right)={\int }_0^L[\sum _{\alpha =1}^M\left(-T_p^{\left(\alpha \right)}\left(F_e^{\left(\alpha \right)}{\overline{\widehat{L}}}_p^{\left(\alpha \right)}{\left(F_e^{\left(\alpha \right)}\right)}^{-1}\right)-T_e^{\left(\alpha \right)}{\widehat{L}}_e^{\left(\alpha \right)}\right)\\ -\left(b_{\phi }\widehat{c}+h\frac{d\widehat{c}}{d{x}_t}\right)]dx.\end{array}$ (12)

The virtual power of the inertia (acceleration) loads is expressed as follows

${\mathcal{P}}_a\left(\widehat{v},\widehat{c}\right)={\int }_0^L\rho \dot{v}\widehat{v}dx.$ (13)

Given the previous expressions, the PVP is written as follows:

${\mathcal{P}}_a\left(\widehat{v},\widehat{c}\right)={\mathcal{P}}_i\left({\overline{\widehat{L}}}_p^{\left(\alpha \right)},{\widehat{L}}^{e\left(\alpha \right)},\widehat{c}\right)+{\mathcal{P}}_e\left(\widehat{v},\widehat{c}\right).$ (14)

Substituting Eqs.(11), (12) and (13) in Eq.(14), we have

$\begin{array}{c}{\int }_0^L\rho \dot{v}\widehat{v}dx={\int }_0^L\rho b\widehat{v}dx+t\left(0\right)\widehat{v}\left(0\right)+t\left(L\right)\widehat{v}\left(L\right)+{\int }_0^L\rho a\widehat{c}dx+t_h\left(0\right)\widehat{c}\left(L\right)+t_h\left(L\right)\widehat{c}\left(L\right)+\\ +{\int }_0^L[\sum _{\alpha =1}^M\left(-T_p^{\left(\alpha \right)}\left(F_e^{\left(\alpha \right)}{\overline{\widehat{L}}}_p^{\left(\alpha \right)}{\left(F_e^{\left(\alpha \right)}\right)}^{-1}\right)-T_e^{\left(\alpha \right)}{\widehat{L}}_e^{\left(\alpha \right)}\right)-\left(b_{\phi }\widehat{c}+h\frac{d\widehat{c}}{d{x}_t}\right)]dx.\end{array}$ (15)

To obtain the dynamic equations, consider that $\widehat{V}=\left(\widehat{v},{\widehat{L}}^e,{\overline{\widehat{L}}}_p^{\left(\alpha \right)},\widehat{c}\right)$ is arbitrary and set ${\overline{\widehat{L}}}_p^{\left(\alpha \right)}$ and $\widehat{c}$ to zero. Therefore,

$\begin{array}{c}{\int }_0^L\rho \dot{v}\widehat{v}dx=-{\int }_0^L\sum _{\alpha =1}^M\left(T_e^{\left(\alpha \right)}\frac{d\widehat{v}}{dx}\right)dx+t\left(0\right)\widehat{v}\left(0\right)+t\left(L\right)\widehat{v}\left(L\right)+{\int }_0^L\rho b\widehat{v}dx.\end{array}$ (16)

Applying integration by parts in the first term on the right hand side and rearranging terms, the following equations are obtained:

$\begin{array}{l}\rho \dot{v}=\sum _{\alpha =1}^M\frac{d{T}_e^{\left(\alpha \right)}}{dx}+\rho b,\\ t\left(0\right)=-\sum _{\alpha =1}^M{T}_e^{\left(\alpha \right)}\left(0\right),\\ t\left(L\right)=\sum _{\alpha =1}^M{T}_e^{\left(\alpha \right)}\left(L\right).\end{array}$ (17)

It can be noticed that $T_e=\sum _{\alpha =1}^M(T_e^{\left(\alpha \right)})$ works as Cauchy stress for the problem.

Next, to find an expression for $T_p=\sum _{\alpha =1}^M{T}_p^{\left(\alpha \right)}$, set $\widehat{v}$ and $\widehat{c}$ to zero, and from Eq. (4)

${\widehat{L}}_e^{\left(\alpha \right)}=-F_e^{\left(\alpha \right)}{\overline{\widehat{L}}}_p^{\left(\alpha \right)}{\left(F_e^{\left(\alpha \right)}\right)}^{-1}=-{\overline{\widehat{L}}}_p^{\left(\alpha \right)}.$ (18)

Equation (15) reduces to

$\begin{array}{c}{\int }_0^L\sum _{\alpha =1}^M\left(-T_p^{\left(\alpha \right)}{\overline{\widehat{L}}}_p^{\left(\alpha \right)}-T_e^{\left(\alpha \right)}{\widehat{L}}_e^{\left(\alpha \right)}\right)dx=0,\end{array}$ (19)

and therefore,

$\begin{array}{c}{T}_e^{\left(\alpha \right)}=T_p^{\left(\alpha \right)}.\end{array}$ (20)

Similarly to the cases above, let $\widehat{v}$, ${\overline{\widehat{L}}}_p^{\left(\alpha \right)}$ and ${\widehat{L}}_e^{\left(\alpha \right)}$ be zero. In this case, applying integration by parts and rearranging the terms, we obtain

$\begin{array}{l}\frac{dh}{dx}-b_{\phi }+\rho a=0,\\ t_h\left(0\right)=-h\left(0\right),\\ t_h\left(L\right)=h\left(L\right).\end{array}$ (21)

Next, the principle of energy conservation is considered to obtain the balance of energy in the bar. This principle describes the evolution of internal energy, given that a thermomechanical system can store kinetic ($\mathcal{K}$) and internal ($e$) energies.

Balance of energy states that the total energy rate is equal to the sum of the external power (${\mathcal{P}}_e$) and the heat flux. Therefore,

$\frac{d}{dt}{\int }_0^L\rho edx+\frac{d}{dt}\mathcal{K}\left(v\right)={\mathcal{P}}_e\left(v,c\right)+{\int }_0^L\rho gdx-\left[q\left(L\right)-q\left(0\right)\right],$ (22)

where the kinetic energy $\mathcal{K}$ does not consider microscopical velocities and is given by $\mathcal{K}\left(v\right)=\frac{1}{2}{\int }_0^L\rho vvdx.$

Moreover, $g$ is the specific heat source density; $e$ is the specific internal energy density, and $q$ is the heat flux (all in Eulerian coordinates).

Previous expression of the balance of energy, combined with the balance of mechanical work, which is obtained from Eq.(15) by taking $\widehat{v}=v$ and $\widehat{c}=\dot{\phi }$, gives the reduced form of the balance of energy in the integral form, given by

$\frac{d}{dt}{\int }_0^L\rho edx=-{\mathcal{P}}_i\left({\overline{L}}_p^{\left(\alpha \right)},L_e^{\left(\alpha \right)},\dot{\phi }\right)+{\int }_0^L\rho gdx-\left[q\left(L\right)-q\left(0\right)\right].$ (23)

Due to the conservation of mass, we have $\frac{d}{dt}{\int }_0^L\rho edx={\int }_0^L\rho \dot{e}dx$, and, using also the divergence theorem, we obtain

$\begin{array}{l}\rho \dot{e}=-\frac{dq}{dx}+\rho g+\sum _{\alpha =1}^M\left(T_p^{\left(\alpha \right)}{\overline{L}}_p^{\left(\alpha \right)}+T_e^{\left(\alpha \right)}{L}_e^{\left(\alpha \right)}\right)+b_{\phi }\dot{\phi }+h\frac{d\dot{\phi }}{dx}.\end{array}$ (24)

Next, the second law of thermodynamics is given in local form of the Clausius-Duhem inequality by

$\rho \dot{\eta }\ge -\frac{d{q}_{\eta }}{dx}+\rho g_{\eta },$ (25)

with $\eta$ the possible specific entropy density, $q_{\eta }$ the entropy flux and $g_{\eta }$ the specific density of sources and sinks of entropy, which is also called entropy production term. Entropy flux and production terms are given, respectively, by

$\begin{array}{l}{q}_{\eta }=\frac{q}{\theta },\mathrm{}{g}_{\eta }=\frac{g}{\theta },\end{array}$ (26)

where $\theta$ is the absolute temperature.

Inequality (25) used with the Coleman-Noll procedure is fundamental for obtaining constitutive relations for stresses in terms of the Helmholtz free-energy, $\Psi$, which is given in Eulerian coordinates by

$\Psi =e-\theta \eta .$ (27)

Using the balance of energy and the Helmholtz free energy definition, dissipation inequality can be rewritten as

$-\rho \dot{\psi }-\rho \eta \dot{\theta }+\sum _{\alpha =1}^M\left(T_p^{\left(\alpha \right)}{\overline{L}}_p^{\left(\alpha \right)}+T_e^{\left(\alpha \right)}{L}_e^{\left(\alpha \right)}\right)+b_{\phi }\dot{\phi }+h\frac{d\dot{\phi }}{dx}-\frac{d\theta }{dx}{q}_{\eta }\ge 0.$ (28)

Collecting the previous results, the basic governing equations in the spatial reference frame are

$\left\{\begin{array}{l}\dot{u}=v,\\ \rho \dot{v}=\sum _{\alpha =1}^M\frac{d{T}_e^{\left(\alpha \right)}}{dx}+\rho b,\\ \frac{dh}{dx}-b_{\phi }+\rho a=0,\\ \rho \dot{e}=-\frac{dq}{dx}+\rho g+\sum _{\alpha =1}^M\left(T_p^{\left(\alpha \right)}{\overline{L}}_p^{\left(\alpha \right)}+T_e^{\left(\alpha \right)}{L}_e^{\left(\alpha \right)}\right)+b_{\phi }\dot{\phi }+h\frac{d\dot{\phi }}{dx},\\ -\rho \dot{\psi }-\rho \eta \dot{\theta }+\sum _{\alpha =1}^M\left(T_p^{\left(\alpha \right)}{\overline{L}}_p^{\left(\alpha \right)}+T_e^{\left(\alpha \right)}{L}_e^{\left(\alpha \right)}\right)+b_{\phi }\dot{\phi }+h\frac{d\dot{\phi }}{dx}-\frac{d\theta }{dx}{q}_{\eta }\ge 0,\\ F^{\left(\alpha \right)}=F_e^{\left(\alpha \right)}{F}_p^{\left(\alpha \right)},\\ T_e^{\left(\alpha \right)}=T_p^{\left(\alpha \right)}.\end{array}\right.$ (29)

The dynamic equations obtained are now written in terms of material coordinates. The terms in this configuration are denoted with the subscript $(.{)}_0$. However, their physical meaning remain unchanged.

The balance of linear momentum in the material configuration is expressed by

$\begin{array}{l}{\rho }_0{\ddot{u}}_0=\sum _{\alpha =1}^M\frac{d{P}^{\left(\alpha \right)}}{dX}+{\rho }_0{b}_0,\\ T_0^{\left(\alpha \right)}\left(0\right)=-P^{\left(\alpha \right)}\left(0\right),\\ T_0^{\left(\alpha \right)}\left(L_0\right)=-P^{\left(\alpha \right)}\left(L_0\right),\end{array}$ (30)

with the Cauchy stress $T_e^{\left(\alpha \right)}$ replaced by the first Piola-Kirchhoff stress $P^{\left(\alpha \right)}$. They are related as

$\begin{array}{l}{P}^{\left(\alpha \right)}=J{T}_e^{\left(\alpha \right)}{F}^{-T},\\ {\rho }_0=\rho J.\end{array}$ (31)

The microscopical balance of Eq. (21) is given now by

$\begin{array}{l}\frac{d{h}_0}{dX}-b_{\phi 0}+{\rho }_0{a}_0=0,\\ t_{h0}\left(0\right)=-h_0\left(0\right),\\ t_{h0}\left(L\right)=h_0\left(L\right),\end{array}$ (32)

with

$h_0=Jh{F}^{-T}.$ (33)

Rewriting the balance of energy of Eq. (24) in material reference leads to following equation:

$\begin{array}{c}{\rho }_0{\dot{e}}_0=-\frac{d{q}_0}{dX}+{\rho }_0{g}_0+\sum _{\alpha =1}^M(J{T}_p^{\left(\alpha \right)}{\overline{L}}_p^{\left(\alpha \right)}+J{T}_e^{\left(\alpha \right)}{L}_e^{\left(\alpha \right)})+b_{\phi 0}\dot{\phi }+h_0\frac{d\dot{\phi }}{dX},\end{array}$ (34)

With

$\begin{array}{c}{q}_0=Jq{F}^{-T}.\end{array}$ (35)

For describing the dissipation inequality it is convenient to use the Clausius-Duhem and the Helmholtz free-energy definitions both in material reference. Therefore,

$\begin{array}{l}{\rho }_0{\dot{\eta }}_0\ge -\frac{d{q}_{\eta 0}}{dX}+{\rho }_0{g}_{\eta 0},\\ {\Psi }_0=e_0-\theta {\eta }_0,\end{array}$ (36)

and again the classical choices for the entropy flux and production terms are considered

$\begin{array}{c}{q}_{\eta 0}=\frac{q_0}{\theta },g_{\eta 0}=\frac{g_0}{\theta },\end{array}$ (37)

resulting in

$-{\rho }_0{\dot{\psi }}_0-{\rho }_0{\eta }_0\dot{\theta }+\sum _{\alpha =1}^M(J{T}_p^{\left(\alpha \right)}{\overline{L}}_p^{\left(\alpha \right)}+J{T}_e^{\left(\alpha \right)}{L}_e^{\left(\alpha \right)})+b_{\phi 0}\dot{\phi }+h_0\frac{d\dot{\phi }}{dX}-\frac{d\theta }{dX}\frac{q_0}{\theta }\ge 0.$ (38)

It is useful to rewrite the term $J{T}_p^{\left(\alpha \right)}{\overline{L}}_p^{\left(\alpha \right)}+J{T}_e^{\left(\alpha \right)}{L}_e^{\left(\alpha \right)}$ of Eq.(38) without variables in the spatial reference. For this purpose, partial right Cauchy-Green ${\overline{C}}_e^{\left(\alpha \right)}$ and $C_p^{\left(\alpha \right)}$ and partial Green-Lagrange ${\overline{E}}_e^{\left(\alpha \right)}$ and $E_p^{\left(\alpha \right)}$ strains are used and defined as

$\begin{array}{c}{\overline{E}}_e^{\left(\alpha \right)}=\frac{1}{2}\left({\overline{C}}_e^{\left(\alpha \right)}-1\right),\mathrm{}{E}_p^{\left(\alpha \right)}=\frac{1}{2}\left(C_p^{\left(\alpha \right)}-1\right).\end{array}$ (39)

The derivatives of the previous expressions are given by

$\begin{array}{c}{\overline{\dot{C}}}_e^{\left(\alpha \right)}=2{\overline{\dot{E}}}_e^{\left(\alpha \right)}=2{\left(F_e^{\left(\alpha \right)}\right)}^T{D}_e^{\left(\alpha \right)}{F}_e^{\left(\alpha \right)},\\ {\dot{C}}_p^{\left(\alpha \right)}=2{\dot{E}}_p^{\left(\alpha \right)}=2{\left(F_p^{\left(\alpha \right)}\right)}^T{\overline{D}}_p^{\left(\alpha \right)}{F}_p^{\left(\alpha \right)}.\end{array}$ (40)

Using these definitions and remembering that $D_e^{\left(\alpha \right)}$ and ${\overline{D}}_p^{\left(\alpha \right)}$ are the symmetric parts of $L_e^{\left(\alpha \right)}$ and ${\overline{L}}_p^{\left(\alpha \right)}$, respectively, and using also the fact that $T_e^{\left(\alpha \right)}$ is symmetric, results in the following identity:

$T_e^{\left(\alpha \right)}{L}_e^{\left(\alpha \right)}=T_e^{\left(\alpha \right)}{D}_e^{\left(\alpha \right)}=T_e^{\left(\alpha \right)}\left({\left(F_e^{\left(\alpha \right)}\right)}^{-T}{\dot{E}}_e^{\left(\alpha \right)}{\left(F_e^{\left(\alpha \right)}\right)}^{-1}\right)={\tilde{T}}_e^{\left(\alpha \right)}{\dot{E}}_e^{\left(\alpha \right)}=\frac{1}{2}{\tilde{T}}_e^{\left(\alpha \right)}{\dot{C}}_e^{\left(\alpha \right)},$ (41)

with

${\tilde{T}}_e^{\left(\alpha \right)}={\left(F_e^{\left(\alpha \right)}\right)}^{-1}{T}_e^{\left(\alpha \right)}{\left(F_e^{e\left(\alpha \right)}\right)}^{-T}=J^{-1}{S}_e^{\left(\alpha \right)},$ (42)

where $S_e^{\left(\alpha \right)}$ is the elastic second Piola-Kirchhoff stress.

Considering the Kroner-Lee decomposition of Eq.(2), the first Piola-Kirchhoff stress is expressed by

$P^{\left(\alpha \right)}=J{F}^{\left(\alpha \right)}{\left(F_p^{\left(\alpha \right)}\right)}^{-1}{\tilde{T}}_e^{\left(\alpha \right)}{\left(F_p^{\left(\alpha \right)}\right)}^{-T},$ (43)

or, equivalently,

$P^{\left(\alpha \right)}=J{F}_e^{\left(\alpha \right)}{\tilde{T}}_e^{\left(\alpha \right)}{\left(F_e^{\left(\alpha \right)}\right)}^T{\left(F^{\left(\alpha \right)}\right)}^{-T}.$ (44)

Thus, the motion equation in Lagrangian coordinates is given by

${\rho }_0{\ddot{u}}_0=\sum _{\alpha =1}^M\frac{d\left(J{F}_e^{\left(\alpha \right)}{\tilde{T}}_e^{\left(\alpha \right)}{\left(F_e^{\left(\alpha \right)}\right)}^T{\left(F^{\left(\alpha \right)}\right)}^{-T}\right)}{dX}+{\rho }_0{b}_0.$ (45)

After simplifications, we obtain for the plastic term that

$T_p^{\left(\alpha \right)}{F}_e^{\left(\alpha \right)}{L}_p^{\left(\alpha \right)}{\left(F_e^{\left(\alpha \right)}\right)}^{-1}=T_p^{\left(\alpha \right)}{L}_p^{\left(\alpha \right)}.$ (46)

Given that the skew-symmetric part of the gradient of velocity $W^p$ is zero because there are no rotations, we have that

$T_p^{\left(\alpha \right)}{D}_p^{\left(\alpha \right)}=T_p^{\left(\alpha \right)}\left({\left(F_p^{\left(\alpha \right)}\right)}^{-T}{\dot{E}}_p^{\left(\alpha \right)}{\left(F_p^{\left(\alpha \right)}\right)}^{-1}\right)={\tilde{T}}_p^{\left(\alpha \right)}{\dot{E}}_p^{\left(\alpha \right)},$ (47)

with

${\tilde{T}}_p^{\left(\alpha \right)}={\left(F_p^{\left(\alpha \right)}\right)}^{-1}{T}_p^{\left(\alpha \right)}{\left(F_p^{\left(\alpha \right)}\right)}^{-T}.$ (48)

It is possible, then, to rewrite the dissipation inequality as

$-{\rho }_0{\dot{\psi }}_0-{\rho }_0{\eta }_0\dot{\theta }+\sum _{\alpha =1}^M(J{\tilde{T}}_p^{\left(\alpha \right)}{\dot{E}}_p^{\left(\alpha \right)}+\frac{1}{2}J{\tilde{T}}_e^{\left(\alpha \right)}{\dot{C}}_e^{\left(\alpha \right)})+b_{\phi 0}\dot{\phi }+h_0\frac{d\dot{\phi }}{dX}-\frac{d\theta }{dX}\frac{q_0}{\theta }\ge 0.$ (49)

Finally, from Eqs.(42) and (48), the relation $T_e^{\left(\alpha \right)}=T_p^{\left(\alpha \right)}$ obtained in the deformed configuration is then redefined in the Lagrangian frame of reference as follows

$F_p^{\left(\alpha \right)}{\tilde{T}}_p^{\left(\alpha \right)}{\left(F_p^{\left(\alpha \right)}\right)}^T=C_e^{\left(\alpha \right)}{\tilde{T}}_e^{\left(\alpha \right)}.$ (50)

Collecting the previous results, the basic governing equations in the material reference frame are

$\left\{\begin{array}{l}{\dot{u}}_0=v,\\ {\rho }_0{\ddot{u}}_0=\sum _{\alpha =1}^M\frac{d\left(J{F}_e^{\left(\alpha \right)}{\tilde{T}}_e^{\left(\alpha \right)}{\left(F_e^{\left(\alpha \right)}\right)}^T{\left(F^{\left(\alpha \right)}\right)}^{-T}\right)}{dX}+{\rho }_0{b}_0,\\ \frac{d{h}_0}{dX}-b_{\phi 0}+{\rho }_0{a}_0=0,\\ {\rho }_0{\dot{e}}_0=-\frac{d{q}_0}{dX}+{\rho }_0{g}_0+\sum _{\alpha =1}^M(J{\tilde{T}}_p^{\left(\alpha \right)}{\dot{E}}_p^{\left(\alpha \right)}+\frac{1}{2}J{\tilde{T}}_e^{\left(\alpha \right)}{\dot{C}}_e^{\left(\alpha \right)})+b_{\phi 0}\dot{\phi }+h_0\frac{d\dot{\phi }}{dX},\\ -{\rho }_0{\dot{\psi }}_0-{\rho }_0{\eta }_0\dot{\theta }+\sum _{\alpha =1}^M(J{\tilde{T}}_p^{\left(\alpha \right)}{\dot{E}}_p^{\left(\alpha \right)}+\frac{1}{2}J{\tilde{T}}_e^{\left(\alpha \right)}{\dot{C}}_e^{\left(\alpha \right)})+b_{\phi 0}\dot{\phi }+h_0\frac{d\dot{\phi }}{dX}-\frac{d\theta }{dX}\frac{q_0}{\theta }\ge 0,\\ F^{\left(\alpha \right)}=F_e^{\left(\alpha \right)}{F}_p^{\left(\alpha \right)},\\ F_p^{\left(\alpha \right)}{\tilde{T}}_p^{\left(\alpha \right)}{\left(F_p^{\left(\alpha \right)}\right)}^T=C_e^{\left(\alpha \right)}{\tilde{T}}_e^{\left(\alpha \right)}.\end{array}\right.$ (51)

4.3 Constitutive model

The balance laws obtained in the previous sections are not enough to fully characterize the physical response of a given body. The material properties must be considered from the constitutive equation. As mentioned earlier, due to complex elasto-viscoplastic response of polymers, a set of internal variables are used and different free energies are applied to each branch of the rheological model.

In order to obtain thermodynamically consistent expressions for constitutive relations, we follow in this section arguments similar to the ones introduced by Truesdell and Noll (1965).

From the dissipation inequality and free-energy density, a set of constitutive equation is defined as

$\left\{\begin{array}{l}\begin{array}{l}\psi =\sum _{\alpha =1}^M{\psi }^{\left(\alpha \right)}\left(\Gamma \right),\mathrm{}\text{with}\mathrm{}\Gamma =\left(\theta ,\phi ,\nabla \theta ,\nabla \phi ,C_e^{\left(\alpha \right)},A^{\left(\alpha \right)}\right),\\ {\tilde{T}}_e^{\left(\alpha \right)}={\tilde{T}}_e^{\left(\alpha \right)}\left(C_e^{\left(\alpha \right)},\theta ,\phi \right),\\ T_p^{\left(\alpha \right)}=T_p^{\left(\alpha \right)}\left(D_p^{\left(\alpha \right)},{\Lambda }^{\left(\alpha \right)}\right),\mathrm{}\text{with}\mathrm{}{\Lambda }^{\left(\alpha \right)}={\Lambda }^{\left(\alpha \right)}\left(C_e^{\left(\alpha \right)},B_p^{\left(\alpha \right)},A^{\left(\alpha \right)},{\mathbf{\xi }}^{\left(\alpha \right)},\theta \right),\\ {\dot{\xi }}_i^{\left(\alpha \right)}=h_i^{\left(\alpha \right)}\left(D_p^{\left(\alpha \right)},{\Lambda }^{\left(\alpha \right)}\right)-R_i^{\left(\alpha \right)}\left({\Lambda }^{\left(\alpha \right)}\right),\\ \dot{\overline{A^{\left(\alpha \right)}}}=D_p^{\left(\alpha \right)}{A}^{\left(\alpha \right)}+A^{\left(\alpha \right)}{D}_p^{\left(\alpha \right)}-G^{\left(\alpha \right)}\left({\Lambda }^{\left(\alpha \right)}\right)d_p^{\left(\alpha \right)}-G_{static}^{\left(\alpha \right)}\left({\Lambda }^{\left(\alpha \right)}\right),\mathrm{}\text{with}\mathrm{}{d}_p^{\left(\alpha \right)}:=\left|D_p^{\left(\alpha \right)}\right|,\end{array}\end{array}\right.$ (52)

where $A^{\left(\alpha \right)}$ is an internal variable and $\xi =\left({\xi }_1,{\xi }_2,\dots ,{\xi }_m\right)$ represents the internal and hidden state variables, whose evolution depends on the dynamic term $h_i^{\left(\alpha \right)}$ and the static recovery $R_i$. These variables are related to the microstructural resistance to plastic flow.

For metals, the Bauschinger-effect is related to decrease in the yield strength when a body is submitted to a certain amount of axial extension into the plastic range and, subsequently, compressed. Due to different microstructural mechanisms in polymers, the same effect can be noticed (even before chain locking in large strains, as described by Hasan and Boyce (1995)).

For taking into account energy storage mechanisms due to plastic deformation, an internal back-stress describes Bauschinger-like phenomena on unloading and reverse loading. For this reason, Anand et al. (2009) introduced variable $A^{\left(\alpha \right)}$ and its evolution equation, which is the same as Eq.(52.v). This variable represents a dimensionless squared stretch-like quantity, whose evolution considers static ($G_{static}^{\left(\alpha \right)}$) and dynamic ($G^{\left(\alpha \right)}$) recovery terms.

Similarly to Frémond and Shitikova (2002), ${\tilde{T}}_e$, $T_p$, $b_{\phi 0}$, $h_0$ and $q_0$ are split in their reversible (non-dissipative) and irreversible (dissipative) parts, which are indicated, respectively, by the superscripts $(.{)}^{\left(r\right)}$ and $(.{)}^{\left(ir\right)}$ as below

$\begin{array}{c}{\tilde{T}}_e={\tilde{T}}_e^{\left(r\right)}+{\tilde{T}}_e^{\left(ir\right)},\\ T_p=T_p^{\left(r\right)}+T_p^{\left(ir\right)},\\ b_{\phi 0}=b_{\phi 0}^{\left(r\right)}+b_{\phi 0}^{\left(ir\right)},\\ h_0=h_0^{\left(r\right)}+h_0^{\left(ir\right)},\\ q_0=q_0^{\left(r\right)}+q_0^{\left(ir\right)}.\end{array}$ (53)

Non-dissipative (reversible) parts may, in general, depend on $\Gamma$, while the dissipative (irreversible) parts may depend both on $\Gamma$ and some of their derivatives in time or space. The non-dissipative terms are expected to give no contribution to entropy increasing, while dissipative terms necessarily contribute to it.

For simplicity, as Frémond and Shitikova (2002), we assume that dissipation (irreversibility) appears only due to $\dot{\phi }$ and $\nabla \theta$ in Eq.(49), and that the heat flux is purely dissipative (irreversible). Thus, we consider

$h_0^{\left(ir\right)}\equiv 0,$ (54)

$q_0^{\left(r\right)}\equiv 0.$ (55)

The expressions for terms in Eqs.(53) must be found such that the entropy condition is satisfied for any admissible process.

Considering the chain rule for ${\psi }^{\left(\alpha \right)}$ and using the evolution law of $A^{\left(\alpha \right)}$ given in Eq.(52.v), we obtain

$\begin{array}{rr}{\dot{\psi }}^{\left(\alpha \right)}\left(\Gamma \right)& =\frac{\partial {\psi }^{\left(\alpha \right)}}{\partial C_e^{\left(\alpha \right)}}{\dot{C}}_e^{\left(\alpha \right)}+\frac{\partial {\psi }^{\left(\alpha \right)}}{\partial A^{\left(\alpha \right)}}{\dot{A}}^{\left(\alpha \right)}+\frac{\partial {\psi }^{\left(\alpha \right)}}{\partial {\theta }^{\left(\alpha \right)}}{\dot{\theta }}^{\left(\alpha \right)}\\ & +\frac{\partial {\psi }^{\left(\alpha \right)}}{\partial \nabla {\theta }^{\left(\alpha \right)}}\nabla {\dot{\theta }}^{\left(\alpha \right)}+\frac{\partial {\psi }^{\left(\alpha \right)}}{\partial {\phi }^{\left(\alpha \right)}}{\dot{\phi }}^{\left(\alpha \right)}+\frac{\partial {\psi }^{\left(\alpha \right)}}{\partial \nabla {\phi }^{\left(\alpha \right)}}\nabla {\dot{\phi }}^{\left(\alpha \right)},\end{array}$ (56)

with

${\frac{\partial {\mathrm{\psi }}^{\left(\mathrm{\alpha }\right)}}{\partial A^{\left(\alpha \right)}}\dot{A}}^{\left(\alpha \right)}=\mathrm{}2\frac{\partial {\mathrm{\psi }}^{\left(\mathrm{\alpha }\right)}}{A^{\left(\alpha \right)}}{A}^{\left(\alpha \right)}{D}_{\mathrm{p}}^{\left(\alpha \right)}-\mathrm{}\frac{\partial {\mathrm{\psi }}^{\left(\alpha \right)}}{\partial A^{\left(\alpha \right)}}\mathrm{}{G}^{\left(\alpha \right)}{\mathrm{d}}_{\mathrm{p}}^{\left(\alpha \right)}-\frac{\partial {\mathrm{\psi }}^{\left(\alpha \right)}}{\partial A^{\left(\alpha \right)}}{G}_{\left\{\mathrm{s}\mathrm{t}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{c}\right\}}^{\left(\alpha \right)}$ (57)

The entropy inequality (written in terms of the free-energy), after some manipulation and collecting similar terms is rewritten as

$\begin{array}{l}-{\rho }_0\left({\eta }_0+\frac{\partial {\psi }^{\left(\alpha \right)}}{\partial {\theta }^{\left(\alpha \right)}}\right)\dot{\theta }+\left(-{\rho }_0\frac{\partial {\psi }^{\left(\alpha \right)}}{\partial {\phi }^{\left(\alpha \right)}}+b_{\phi 0}^{\left(r\right)}+b_{\phi 0}^{\left(ir\right)}\right)\dot{\phi }-\left({\rho }_0\frac{\partial {\psi }^{\left(\alpha \right)}}{\partial \nabla {\phi }^{\left(\alpha \right)}}-h_0^{\left(r\right)}\right)\nabla \dot{\phi }\\ +\left(\frac{1}{2}J{\tilde{T}}_e^{\left(r\right)}+\frac{1}{2}J{\tilde{T}}_e^{\left(ir\right)}-{\rho }_0\frac{\partial {\psi }^{\left(\alpha \right)}}{C_e^{\left(\alpha \right)}}\right){\dot{C}}_e^{\left(\alpha \right)}+\left(J{T}_p^{\left(r\right)}-2{\rho }_0\frac{\partial {\psi }^{\left(\alpha \right)}}{\partial A^{\left(\alpha \right)}}{A}^{\left(\alpha \right)}\right)D_p^{\left(\alpha \right)}\\ -{\rho }_0\frac{\partial {\psi }^{\left(\alpha \right)}}{\partial \nabla {\theta }^{\left(\alpha \right)}}\nabla \dot{\theta }-{\rho }_0\left(\frac{\partial {\psi }^{\left(\alpha \right)}}{\partial A^{\left(\alpha \right)}}G\right)d_p^{\left(\alpha \right)}+J{T}_p^{\left(ir\right)}{D}_p^{\left(\alpha \right)}-\frac{1}{\theta }{q}^{\left(ir\right)}\frac{d\theta }{dX}\ge 0.\end{array}$ (58)

Next, we choose the reversible terms of the last inequality such that they do not contribute to the increase of entropy for any admissible process, that is,

$\begin{array}{c}-{\rho }_0\left({\eta }_0+\frac{\partial {\psi }^{\left(\alpha \right)}}{\partial {\theta }^{\left(\alpha \right)}}\right)\dot{\theta }+\left(-{\rho }_0\frac{\partial {\psi }^{\left(\alpha \right)}}{\partial {\phi }^{\left(\alpha \right)}}+b_{\phi 0}^{\left(r\right)}\right)\dot{\phi }-\left({\rho }_0\frac{\partial {\psi }^{\left(\alpha \right)}}{\partial \nabla {\phi }^{\left(\alpha \right)}}-h_0^{\left(r\right)}\right)\nabla \dot{\phi }\\ +\left(\frac{1}{2}J{\tilde{T}}_e^{\left(r\right)}-{\rho }_0\frac{\partial {\psi }^{\left(\alpha \right)}}{\partial C_e^{\left(\alpha \right)}}\right){\dot{C}}_e^{\left(\alpha \right)}+\left(J{T}_p^{\left(r\right)}-2{\rho }_0\frac{\partial {\psi }^{\left(\alpha \right)}}{\partial A^{\left(\alpha \right)}}{A}^{\left(\alpha \right)}\right)D_p^{\left(\alpha \right)}-{\rho }_0\frac{\partial {\psi }^{\left(\alpha \right)}}{\partial \nabla {\theta }^{\left(\alpha \right)}}\nabla \dot{\theta }=0.\end{array}$ (59)

Since in Eq.(59) we can give any values to $\dot{\theta }$, $\nabla \dot{\theta }$, $\dot{\phi }$, $\nabla \dot{\phi }$, ${\dot{C}}^e$ and ${\dot{D}}^p$ their coefficients can be set to zero. Thus,

${\eta }_0=-\frac{\partial {\psi }^{\left(\alpha \right)}}{\partial {\theta }^{\left(\alpha \right)}},$ (60)

$\frac{\partial {\psi }^{\left(\alpha \right)}}{\partial \nabla {\theta }^{\left(\alpha \right)}}=0,$ (61)

$b_{\phi 0}^{\left(r\right)}={\rho }_0\frac{\partial {\psi }^{\left(\alpha \right)}}{\partial {\phi }^{\left(\alpha \right)}},$ (62)

${\tilde{T}}_e^{\left(r\right)}=2J^{-1}{\rho }_0\frac{\partial {\psi }^{\left(\alpha \right)}}{\partial C_e^{\left(\alpha \right)}},$ (63)

$T_p^{\left(r\right)}=2J^{-1}{\rho }_0\frac{\partial {\psi }^{\left(\alpha \right)}}{\partial A^{\left(\alpha \right)}}{A}^{\left(\alpha \right)},$ (64)

$h_0=h_0^{\left(r\right)}={\rho }_0\frac{\partial {\psi }^{\left(\alpha \right)}}{\partial \nabla {\phi }^{\left(\alpha \right)}}.$ (65)

By using the previous results, Eq.(58) is reduced to

$\begin{array}{c}{b}_{\phi 0}^{\left(ir\right)}\dot{\phi }+\frac{1}{2}J{\tilde{T}}_e^{\left(ir\right)}{\dot{C}}_e^{\left(\alpha \right)}+{\rho }_0\left(\frac{\partial {\psi }^{\left(\alpha \right)}}{\partial A^{\left(\alpha \right)}}G\right)d^p+J{T}_p^{\left(ir\right)}{D}_p^{\left(\alpha \right)}-\frac{1}{\theta }{q}^{\left(ir\right)}\frac{1}{\theta }\frac{d\theta }{dX}\ge 0.\end{array}$ (66)

The dissipation inequality above is, then, split into four inequalities, being the first related to the mechanical dissipation

$\begin{array}{c}{Y}_p^{\left(\alpha \right)}{D}_p^{\left(\alpha \right)}+{\rho }_0\left(\frac{\partial {\psi }^{\left(\alpha \right)}}{\partial A^{\left(\alpha \right)}}G\right)d_p^{\left(\alpha \right)}\ge 0,\end{array}$ (67)

with

$\begin{array}{c}{Y}_p^{\left(\alpha \right)}=T_p^{\left(\alpha \right)}-2J^{-1}{\rho }_0\frac{\partial {\psi }^{\left(\alpha \right)}}{\partial A^{\left(\alpha \right)}}{A}^{\left(\alpha \right)},\end{array}$ (68)

representing the dissipative part of $T_p^{\left(\alpha \right)}$. Thus, it can be noticed that the plastic stress $T_p^{\left(\alpha \right)}$, conjugated to $D_p^{\left(\alpha \right)}$, is composed by a dissipative and a reversible (energetic) terms as

$T_p^{\left(\alpha \right)}\left(D_p^{\left(\alpha \right)},{\Lambda }^{\left(\alpha \right)}\right)=Y_p^{\left(\alpha \right)}\left(D_p^{\left(\alpha \right)},{\Lambda }^{\left(\alpha \right)}\right)+\left(2J^{-1}{\rho }_0\frac{\partial {\psi }^{\left(\alpha \right)}}{\partial A^{\left(\alpha \right)}}{A}^{\left(\alpha \right)}\right).$ (69)

The energetic part is also called back-stress, and given by

$M_{back}^{\left(\alpha \right)}=2J^{-1}{\rho }_0\left(\frac{\partial {\psi }^{\left(\alpha \right)}\left(C_e^{\left(\alpha \right)},A^{\left(\alpha \right)},\theta \right)}{\partial A^{\left(\alpha \right)}}{A}^{\left(\alpha \right)}\right).$ (70)

Proceeding with the analysis of the dissipation inequality, the last term in Eq.(66) is the heat conduction inequality, whose positiveness is imposed by

$q^{\left(ir\right)}=-\tilde{c}\frac{1}{\theta }\frac{d\theta }{dX},$ (71)

with $\tilde{c}$ a non-negative coefficient.

Analogously, the first term in Eq.(66) corresponds to the microscopic energy inequality, whose positiveness is imposed by

$b_{\phi 0}^{\left(ir\right)}=\tilde{\gamma }\dot{\phi },$ (72)

with $\tilde{\gamma }$ a non-negative coefficient.

For this work, the irreversible part of ${\tilde{T}}_e$ to ensure positiveness is

${\tilde{T}}_e^{\left(ir\right)}=2{\tilde{b}}_d{\dot{C}}_e^{\left(\alpha \right)}.$ (73)

An important hypothesis, to be used for the specific constitutive models in next sections, is that the free energy can be separated in three independent terms for each branch $\alpha$, being the first one an elastic energy (isotropic with relation to $C_e^{\left(\alpha \right)}$), the second a defect energy (isotropic with relation to $A^{\left(\alpha \right)}$), associated to plastic flow, and the third related to damage. Therefore,

$\psi ={\psi }_e^{\left(\alpha \right)}\left(C_e^{\left(\alpha \right)},\theta ,\phi \right)+{\psi }_p^{\left(\alpha \right)}\left(A^{\left(\alpha \right)},\theta \right)+{\psi }_d^{\left(\alpha \right)}\left(\phi ,\nabla \phi \right).$ (74)

For evaluation of stresses in plastic flow, considering that $M_e^{\left(\alpha \right)}$ is the Mandel stress and defined as $M_e^{\left(\alpha \right)}=J{C}_e^{\left(\alpha \right)}{\tilde{T}}_e^{\left(r\right)\left(\alpha \right)}$, the effective Mandel stress $M_{\left(eff\right)}^{\left(\alpha \right)}$ is given by

$M_{\left(eff\right)}^{\left(\alpha \right)}=M_e^{\left(\alpha \right)}-M_{back}^{\left(\alpha \right)}.$ (75)

From the constitutive relation given in Eq.(69) follows the definition:

$M_{\left(eff\right)}^{e\left(\alpha \right)}=Y^{p\left(\alpha \right)}\left(D_p^{\left(\alpha \right)},{\Lambda }^{\left(\alpha \right)}\right),$ (76)

with $Y$ being a scalar flow strength of the material for each $\alpha$.

Given the relation for $D_p^{\left(\alpha \right)}$, using $N_p^{\left(\alpha \right)}$ as the direction for plastic flow for each $\alpha$ (remembering that in 1D domain it is the sign of driving stress), the flow rule is given by

$\begin{array}{l}{\dot{F}}_p^{\left(\alpha \right)}=D_p^{\left(\alpha \right)}{F}_p^{\left(\alpha \right)},\\ \\ D_p^{\left(\alpha \right)}=d_p^{\left(\alpha \right)}{N}_p^{\left(\alpha \right)},\\ \\ N_p^{\left(\alpha \right)}=\frac{\left(M_{eff}^{\left(\alpha \right)}\right)}{\left(\left|M_{eff}^{\left(\alpha \right)}\right|\right)}=\text{sign}\left(M_{eff}^{e\left(\alpha \right)}\right).\end{array}$ (77)

5 CONSTITUTIVE MODEL FOR POLYMERS

A large number of constitutive models for polymers, specially for amorphous, can be found in the literature. Most of these models, for simplicity reasons, make assumptions in order to reduce material parameters, many of them being difficult to be predicted or collected through experimental tests.

Considering the PTFE features, with its large dependence on temperature due to internal micromechanisms, it is crucial to evaluate constitutive models in order to be able to predict mechanical behaviour under a broader range of temperature, despite the elevated number of constants to be found.

Earlier studies considering polymer hardening in finite deformation were performed by Boyce, Parks and Argon (1988) and Boyce, Parks and Argon (1989) , which considered the same behaviour originally modelled for rubbers. These models use an analogy of entropic spring as polymer branches.

Models for amorphous polymers take into consideration important changes in mechanical properties next to $T_g$, considering that from this temperature on, constrained polymer chains acquire rotational and translational movement. Being PTFE (semi-crystalline material) our focus, however, adapting and understanding the validity of these models turns out to be essential.

This work uses a constitutive model able to reproduce polymer behaviour in a broad range of temperatures, proposed by Srivastava et al. (2010). It will be described briefly here in the uniaxial context, trying to expose only fundamental topics for their understanding. For a deeper explanation on the derivation of all equations, it is suggested to read the original paper. This model uses the multimechanism concept, where each one has a specific deformation nature, aiming to represent different sources of resistance in polymers (intramolecular and intermolecular), as suggested earlier by Haward and Thackray (1968). The dependence of the free energy on the phase-field used to describe damage is not included yet; it will be added to the model in next sections.

Srivastava et al. (2010) proposed a model composed by three branches, as shown in Fig. 4. In the first, there is a non-linear spring representing elastic resistance to intermolecular energetic bond-stretching and a dashpot representing thermally-activated plastic flow due to inelastic mechanisms, such as chain-segment rotation and relative slippage of the polymer chains between neighboring mechanical cross-linkage points. There is also another non-linear spring in parallel with the dashpot, representing energy storage mechanism caused by the viscoplastic flow mechanisms.

Figure 4
Esquematic model of the amorphous constitutive model by Srivastava et al. (2010).

Branches 2 and 3 represent responses under and above $T_g$, respectively, given the abrupt changes in the behaviour of the material. Both branches are composed by non-linear springs in series with dashpots. Non-linear springs represent resistances due to changes in the free energy upon stretching of the molecular chains between the mechanical cross-links and dashpots represent thermally-activated plastic flow due to slippage of the mechanical cross-links.

The free energies of the three branches, and their specificities, for this model are described next.

5.1 Free energy and stresses for branch 1

Consider a bar subjected to an incompressible deformation. The focus of one-dimensional analysis is on effects in longitudinal direction of the domain and, thus, only $F_{11}=\lambda$ is considered in the model. Moreover, given the Kroner-Lee decomposition (Eq.(3)) with simplification given in Eq.(7), we have

$\begin{array}{rr}{F}^{\left(\alpha \right)}=& {\lambda }_e^{\left(\alpha \right)}{\lambda }_p^{\left(\alpha \right)}.\end{array}$ (78)

In the first branch, the elastic free energy is given by

$\begin{array}{l}{\overline{\psi }}^{e\left(1\right)}=\frac{1}{2}E(E_e^{\left(1\right)}{)}^2-E{\alpha }^{th}(\theta -{\theta }_0\left)\right(E_e^{\left(1\right)}),\end{array}$ (79)

with $E_e^{\left(1\right)}$ being the elastic logarithmic strain measure, $E$ is the Young modulus, ${\alpha }^{th}$ is the coefficient of thermal expansion and ${\theta }_0$ is the reference temperature.

The dependence on temperature of $E$ (following Dupaix and Boyce (2007)) is given by $E\left(\theta \right)=\frac{1}{2}\left(E_{gl}+E_r\right)-\frac{1}{2}\left(E_{gl}-E_r\right)\text{tanh}\left(\frac{1}{\Delta }\left(\theta -{\theta }_g\right)\right)-M\left(\theta -{\theta }_g\right),$with $E_{gl}$ and $E_r$ being the Young modulus in glassy and rubbery regions, respectively. Parameter $\Delta$ is related to size of transition close to ${\theta }_g$, and $M$ represents slope of the shear modulus in terms of temperature, being divided in two regions (pre and post ${\theta }_g$).

The defect free energy in the first branch is

$\begin{array}{l}{\overline{\psi }}^{p\left(1\right)}=\frac{1}{4}B\left(\theta \right)\left[(\text{ln}{a}_1{)}^2+(\text{ln}{a}_2{)}^2+(\text{ln}{a}_3{)}^2\right],\\ \\ B\left(\theta \right)=X\left(\theta -{\theta }_0\right)\text{for}\theta <{\theta }_0\end{array}$ (80)

with $a_i\left(i=\mathrm{1,2},3\right)$ given under uniaxial conditions by $a_1=A^{\left(1\right)}$, and $a_2=a_3={\left(A^{\left(1\right)}\right)}^{-1/2}$, in accordance with Anand et al. (2009). Parameter $B\left(\theta \right)$ is the backstress modulus, assumed to be a linearly decreasing function of temperature with slope $X$ (to be defined as material coefficient), when $\theta <{\theta }_g$. This backstress modulus is assumed zero when $\theta \ge {\theta }_g$.

The Mandel stress is defined as

$M_e^{\left(1\right)}=E\left(E_0^{\left(1\right)}\right)-E{\alpha }^{th}\left(\theta -{\theta }_0\right).$ (81)

The constitutive relation of Eq.(70) is used to define the backstress and the defect free-energy expression of Eq.(80.(i)). After some simplification and rearranging terms, we obtain

$M_{back}=\frac{3}{2}B\text{ln}\left(A^{\left(1\right)}\right).$ (82)

Finally, the driving stress for plastic flow is given by the difference between Mandel stress and backstress as $M_{eff}=M_e^{\left(1\right)}-M_{back}.$

5.2 Free energies and stresses for branches 2 and 3

The free energies for branches 2 and 3 are similar and based on the free energy. Therefore,

$\begin{array}{l}{\overline{\psi }}^{e\left(2\right)}=-\frac{1}{2}{\mu }^{\left(2\right)}{I}_m^{\left(2\right)}\text{ln}\left(1-\frac{I_1^{\left(2\right)}-3}{I_m^{\left(2\right)}}\right),\end{array}$ (83)

$\begin{array}{l}{\overline{\psi }}^{e\left(3\right)}=-\frac{1}{2}{\mu }^{\left(3\right)}{I}_m^{\left(3\right)}\text{ln}\left(1-\frac{I_1^{\left(3\right)}-3}{I_m^{\left(3\right)}}\right),\end{array}$ (84)

with ${\mu }^{\left(\alpha \right)}$ representing the ground state rubbery shear modulus of the material and $I_m^{\left(\alpha \right)}$ is the maximum value of $\left(I_1^{\left(\alpha \right)}-3\right)$, associated with the limited extensibility of the polymer chains. The term $I_1^{\left(\alpha \right)}$ is given by

$I_1^{\left(\alpha \right)}=({\lambda }_e^{\left(\alpha \right)}{)}^2+\frac{2}{{\lambda }_e^{\left(\alpha \right)}},$ (85)

where ${\lambda }_e^{\left(\alpha \right)}$ is given by Eq.(78).

The stresses for branches 2 and 3 can be simplified to

$T^{\left(\alpha \right)}={\mu }^{\left(\alpha \right)}{\left(1-\frac{I_1^{\left(\alpha \right)}-3}{I_m^{\left(\alpha \right)}}\right)}^{-1}\left(({\lambda }_e^{\left(\alpha \right)}{)}^2-({\lambda }_e^{\left(\alpha \right)}{)}^{-1}\right),$ (86)

with $\alpha$=2 or 3.

5.3 Flow rules

The evolution equation in the first branch is given by

${\dot{F}}_p^{\left(1\right)}=D_p^{\left(1\right)}{F}_p^{\left(1\right)},$ (87)

with $D_p^{\left(1\right)}$ given by

$D_p^{\left(1\right)}={\nu }^{p\left(1\right)}=\left\{\begin{array}{ll}0& \text{if}{\overline{\tau }}^{e\left(1\right)}\le 0,\\ {\nu }_0^{\left(1\right)}\text{exp}\left(-\frac{Q}{k_b\theta }\right){\left[\text{sinh}\left(\frac{{\tau }_e^{\left(1\right)}V}{2k_b\theta }\right)\right]}^{1/m_1}& \text{if}{\overline{\tau }}^{e\left(1\right)}>0,\end{array}\right.$ (88)

with ${\nu }^{p\left(1\right)}$ the equivalent plastic shear strain rate; ${\overline{\tau }}^{e\left(1\right)}$ is the equivalent shear stress, given by

$\begin{array}{l}{\overline{\tau }}^{e\left(1\right)}=M_{eff}-\left(S_b+{\alpha }_p\overline{p}\right),\\ \overline{p}=-\frac{1}{3}\text{tr}{M}_e^{\left(1\right)}.\end{array}$ (89)

In Eq.(88), several material parameters must be calibrated, since their physical meaning are not easily extracted from experiments. ${\nu }_0^{\left(1\right)}$ is a pre-exponential factor, $Q$ is an activation energy, $k_b$ is the Boltzmann constant, $V$ is activation volume, $m_1$ is a strain-rate sensivity parameter and ${\alpha }_p$ is the pressure sensivity of plastic flow.

The term $S_b$ in Eq.(89) is the dissipative resistance to plastic flow, used to model hardening at large strains, considering the friction among polymeric chains. This term avoids an excessive elastic recovery, after the unloading in large strains. Its evolution is given by the differential equation

$\begin{array}{l}{\dot{S}}_b=h_b\left(\overline{\lambda }-1\right)\left(S_b^{\mathrm{*}}-S_b\right){\nu }^{p\left(1\right)},\\ \\ S_b^{\mathrm{*}}=\frac{1}{2}\left(S_{gl}+S_r\right)-\frac{1}{2}\left(S_{gl}-S_r\right)\text{tanh}\left(\frac{1}{\Delta }\left(\theta -{\theta }_g\right)\right),\end{array}$ (90)

which has others material parameters ($h_b$, $S_{gl}$, $S_r$) to be defined and $\overline{\lambda }$ given by

$\overline{\lambda }=\sqrt{\frac{1}{3}\left({\lambda }_e^2+{\frac{2}{\lambda }}_e\right)}.$ (91)

Still in branch 1, the evolution of $A$ is given by

${\dot{A}}^{\left(1\right)}=D_p^{\left(1\right)}{A}^{\left(1\right)}+A^{\left(1\right)}{D}_p^{\left(1\right)}-\gamma A^{\left(1\right)}\left(\text{ln}{A}^{\left(1\right)}\right){\nu }_p^{\left(1\right)},$ (92)

with $\gamma$ being a parameter driving the dynamic recovery of $A$.

The flow rules for branches 2 and 3 are given by the followin simple power laws, respectively:

${\nu }^{p\left(2\right)}=\frac{1}{2}{\nu }_0^{\left(2\right)}{\left(\frac{{\overline{\tau }}^{\left(2\right)}}{S^{\left(2\right)}}\right)}^{1/m2}\left[1+\text{tanh}\left(\frac{1}{{\Delta }_x}\left(\theta -{\theta }_g\right)\right)\right],$ (93)

${\nu }^{p\left(3\right)}=\frac{1}{2}{\nu }_0^{\left(3\right)}{\left(\frac{{\overline{\tau }}^{\left(3\right)}}{S^{\left(3\right)}}\right)}^{1/m3}\left[1+\text{tanh}\left(\frac{1}{{\Delta }_x}\left(\theta -{\theta }_g\right)\right)\right],$ (94)

with ${\nu }_0^{\left(\alpha \right)}$ being the reference plastic shear strain rate and $m\left(\alpha \right)$ the positive-valued strain-rate sensitivity parameter. $S^{\left(2\right)}$ is a positive-valued shear resistance. In branch 3, which is linked to behavior in temperatures above $T_g$, the internal variable $S^{\left(3\right)}$ is related to a dissipative resistance during molecules movement, known as reptation movement. Its evolution is given by the following rule:

${\dot{S}}^{\left(3\right)}=h_3\left(\overline{\lambda }-1\right){\nu }^{p\left(3\right)},$ (95)

with $\overline{\lambda }$ given in Eq.(91) and the initial conditions for $S^{\left(3\right)}$ and parameter $h_3$ given by

$\begin{array}{l}{S}^{\left(3\right)}=S_g^{\left(3\right)}\text{exp}\left(-Y\left(\theta -{\theta }_g\right)\right),\\ h_3=h_{3g}\text{exp}\left(-Z\left(\theta -{\theta }_g\right)\right),\end{array}$ (96)

where $Y$ and $Z$ are parameters to be found in the calibration process.

6 FREE ENERGY DENSITIES INCLUDING DAMAGE

The free energy of the material including the phase-field damage variable is

$\psi \left(\phi ,\nabla \phi ,C_e,A,C\right)={\psi }^{e\left(1\right)}\left(\phi ,C_e^{\left(1\right)}\right)+{\psi }_p^{\left(1\right)}\left(A\right)+{\psi }_e^{\left(\alpha \right)}\left(C\right)+\mathcal{J}\left(\phi ,{\nabla }_p\phi \right),$ (97)

with $\alpha$ being 2 or 3, depending on the constitutive model considered.

The free energy related to elasticity corresponds to the first term on the right hand side of Eq. (97), expressed by a degradation term ($D\left(\varphi \right)$) and an energetic term (${\psi }_0^{e\left(1\right)}$) as

${\psi }_e^{\left(1\right)}\left(\phi ,C_e\right)=D\left(\phi \right){\psi }_0^{e\left(1\right)}\left(\theta ,C_e^{\left(1\right)}\right).$ (98)

The energetic terms were already defined for the two constitutive models described in previous sections. The degradation function, for reproducing the experimental behaviour, must induce the failure in an abrupt manner. Its graph is shown in Fig. 5, and it is defined as

$D\left(\phi \right)=({\phi }^2-1{)}^2-a\phi (\phi -1{)}^2,$ (99)

where $a>0$ is related to the damage rate for the virgin material that is, when $\phi =0$. In Fig. 5, we show the graph of this function for several values of $a$; in this work, we take $a=0.001.$

Figure 5
Degradation function for different values of parameter a.

Further explanation about why we choose (99) as the degradation function can be found in Section 9.

The free energy related to damage is

$\mathcal{J}\left(\phi ,\nabla \phi \right)=g_c\left(\frac{\gamma }{2}{\left|\frac{d\phi }{dX}\right|}^2+\frac{c}{\gamma }\mathcal{H}\left(\phi \right)\right),$ (100)

where $g_c$ is the Griffith’s critical energy release rate, $\gamma$ a positive constant related to fracture layer width and $c$ a (small) positive constant equal to 0.001 in this work. The expression for the potential $\mathcal{H}$ is given by

$\mathcal{H}\left(\phi \right)=\left\{\begin{array}{ll}\frac{1}{2}{\phi }^2& \text{for}0\le \phi \le 1,\\ 1& \text{for}\phi >1,\\ 0& \text{for}\phi <0.\end{array}\right.$ (101)

Further explanation about the choice of expression (101) for the damage energy can be found in Section 10.

The considered degradation function and the damage free energy are able to correctly describe an abrupt rupture of the material under tension, as noticed, for instance, in experimental results at -10°C by Gamboni et al. (2016).

7 FINITE ELEMENT APPROXIMATION

7.1 Approximation of damage equation

Using the constitutive relations for variables $h_0$, $b_{\phi 0}$ and $a_0=0$ from Eq.(51.(iii)) and applying the free energy defined in Eq. (97), the damage equation is rewritten as

$\begin{array}{rr}\dot{\phi }& =\frac{\gamma g_c}{\tilde{\lambda }}\Delta \phi -\frac{1}{\tilde{\lambda }}\left(\phi -1\right)\left(a+4\phi -3a\phi +4{\phi }^2\right){\psi }_0^{e\left(1\right)}-\frac{g_{c}c}{\tilde{\lambda }\gamma }{\mathcal{H}}^{\mathrm{\text{'}}}\left(\phi \right).\end{array}$ (102)

The previous expression is multiplied by a test function $w$ at both sides and integrated in the bar length $L_0$. After integrating by parts, the first term on the right-hand side, we have

$\begin{array}{l}{\int }_{L_0}\dot{\phi }wd{L}_0=-{\int }_{L_0}\frac{\gamma g_c}{\tilde{\lambda }}\frac{d\phi }{dX}\frac{dw}{dX}d{L}_0-{\int }_{L_0}\frac{1}{\tilde{\lambda }}\left(\phi -1\right)\left(a+4\phi -3a\phi +4{\phi }^2\right){\psi }_0^{e\left(1\right)}wd{L}_o\\ -{\int }_{L_0}\frac{g_{c}c}{\tilde{\lambda }\gamma }{\mathcal{H}}^{\mathrm{\text{'}}}\left(\phi \right)wd{L}_0+\frac{\gamma g_c}{\tilde{\lambda }}\left(\frac{d\phi }{dX}\left(L_0\right)w\left(L_0\right)-\frac{d\phi }{dX}\left(0\right)w\left(0\right)\right).\end{array}$ (103)

The functions with non-linear terms in $\phi$, indicated generically as $f\left(\phi \right)$, and the test functions $w$ are approximated in a generic element $e$ of the mesh, respectively, by

$\begin{array}{l}{f}^{\left(e\right)}\left(\phi \right)=\sum _{i=1}^{P+1}{N}_i^{\left(e\right)}\left(\xi \right)f_i^{\left(e\right)}\left(\phi \right)=\left[{\mathbf{N}}^{\left(e\right)}\right]{\{\mathbf{\phi }}^{\left(e\right)}\},\\ w^{\left(e\right)}=\sum _{i=1}^{P+1}{N}_i^{\left(e\right)}\left(\xi \right)w_i^{\left(e\right)}=\left[{\mathbf{N}}^{\left(e\right)}\right]{\{\mathbf{w}}^{\left(e\right)}\},\end{array}$ (104)

where $N_i^{\left(e\right)}$ are the element shape functions, $P$ is the polynomial order and $P+1$ is the number of nodes of the element. The approximation of the derivatives of $\phi$ and $w$ are given, respectively, by

$\begin{array}{l}\nabla {\phi }^{\left(e\right)}=\sum _{i=1}^{P+1}{N}_{i,\xi }^{\left(e\right)}\left(\xi \right){\phi }_i^{\left(e\right)}=\sum _{i=1}^{P+1}{B}_i^{\left(e\right)}\left(\xi \right){\phi }_i^{\left(e\right)}=\left[{\mathbf{B}}^{\left(e\right)}\right]{\{\mathbf{\phi }}^{\left(e\right)}\},\\ \nabla w^{\left(e\right)}=\sum _{i=1}^{P+1}{N}_{i,\xi }^{\left(e\right)}\left(\xi \right)w_i^{\left(e\right)}=\sum _{i=1}^{P+1}{B}_i^{\left(e\right)}\left(\xi \right)w_i^{\left(e\right)}=\left[{\mathbf{B}}^{\left(e\right)}\right]{\{\mathbf{w}}^{\left(e\right)}\},\end{array}$ (105)

where $B_i^{\left(e\right)}\left(\xi \right)$ are the global derivatives of the element shape functions.

Using the Galerking method, we have

$\begin{array}{l}{\int }_{L{e}_0}[{\mathbf{N}}^{\left(e\right)}{]}^T\left[{\mathbf{N}}^{\left(e\right)}\right]\left\{{\dot{\mathbf{\phi }}}^{\left(e\right)}\right\}d{L}_0=-{\int }_{L{e}_0}\frac{\gamma g_c}{\tilde{\lambda }}\left[{\mathbf{B}}^{\left(e\right)}{]}^T\right[{\mathbf{B}}^{\left(e\right)}\left]\right\{{\mathbf{\phi }}^{\left(e\right)}\}d{L}_0\\ -{\int }_{L{e}_0}\frac{1}{\tilde{\lambda }}(-a[{\mathbf{N}}^{\left(e\right)}{]}^T+(4a-4)\left[{\mathbf{N}}^{\left(e\right)}{]}^T\right[{\mathbf{N}}^{\left(e\right)}\left]\right\{{\mathbf{\phi }}^{\left(e\right)}\}\\ -3a\left[{\mathbf{N}}^{\left(e\right)}{]}^T\right[{\mathbf{N}}^{\left(e\right)}\left]\right\{{\mathbf{\phi }}^{\left(e\right)}{\}}^2+4\left[{\mathbf{N}}^{\left(e\right)}{]}^T\right[{\mathbf{N}}^{\left(e\right)}\left]\right\{{\mathbf{\phi }}^{\left(e\right)}{\}}^3){\psi }_0^{e\left(1\right)}d{L}_o-{\int }_{L{e}_0}\frac{g_{c}c}{\tilde{\lambda }\gamma }\left(\left[{\mathbf{N}}^{\left(e\right)}{]}^T\right[{\mathbf{N}}^{\left(e\right)}\left]\right\{{\mathbf{\phi }}^{\left(e\right)}\}\right)d{L}_0.\end{array}$ (106)

Grouping terms and rearranging, the previous equation is rewritten as

$\begin{array}{rr}\left[{\mathbf{M}}^{\left(e\right)}\right]{\dot{\mathbf{\phi }}}^{\left(e\right)}=& \left[{\mathbf{K}}_1^{\left(e\right)}\right]\left\{{\mathbf{\phi }}^{\left(e\right)}\right\}+\left[{\mathbf{K}}_2^{\left(e\right)}\right]\{{\mathbf{\phi }}^{\left(e\right)}{\}}^2+[{\mathbf{K}}_3^{\left(e\right)}\left]\right\{{\mathbf{\phi }}^{\left(e\right)}{\}}^3+\left[{\mathbf{K}}_4^{\left(e\right)}\right],\end{array}$ (107)

with

$\left[{\mathbf{M}}^{\left(e\right)}\right]={\int }_{L{e}_0}[{\mathbf{N}}^{\left(e\right)}{]}^T\left[{\mathbf{N}}^{\left(e\right)}\right],$ (108)

$\left[{\mathbf{K}}_1^{\left(e\right)}\right]=-{\int }_{L{e}_0}\frac{\gamma g_c}{\tilde{\lambda }}\left[{\mathbf{B}}^{\left(e\right)}{]}^T\right[{\mathbf{B}}^{\left(e\right)}]+\frac{\left(4a-4\right)}{\tilde{\lambda }}[{\mathbf{N}}^{\left(e\right)}{]}^T\left[{\mathbf{N}}^{\left(e\right)}\right]{\psi }_0^{e\left(1\right)}+\frac{g_{c}c}{\tilde{\lambda }\gamma }\left[{\mathbf{N}}^{\left(e\right)}{]}^T\right[{\mathbf{N}}^{\left(e\right)}]d{L}_0,$ (109)

$\left[{\mathbf{K}}_2^{\left(e\right)}\right]={\int }_{L{e}_0}\frac{3a}{\tilde{\lambda }}\left[{\mathbf{N}}^{\left(e\right)}{]}^T\right[{\mathbf{N}}^{\left(e\right)}]{\psi }_0^{e\left(1\right)}d{L}_0,$ (110)

$\left[{\mathbf{K}}_3^{\left(e\right)}\right]=-{\int }_{L{e}_0}\frac{4}{\tilde{\lambda }}\left[{\mathbf{N}}^{\left(e\right)}{]}^T\right[{\mathbf{N}}^{\left(e\right)}]{\psi }_0^{e\left(1\right)}d{L}_0,$ (111)

$\left[{\mathbf{K}}_4^{\left(e\right)}\right]={\int }_{L{e}_0}\frac{a}{\tilde{\lambda }}[{\mathbf{N}}^{\left(e\right)}{]}^{T}d{L}_0.$ (112)

The global equations are obtained using the standard assembly procedure of the element operators given in Eq.(107).

The time integration uses the implicit Euler algorithm. Thus, the Eq.(107) is rewritten as follows:

$g\left(\phi \right)=-{\mathbf{\phi }}_{k+1}^{\left(e\right)}+{\mathbf{\phi }}_k^{\left(e\right)}+\Delta t{\left[{\mathbf{M}}^{\left(e\right)}\right]}^{-1}\left(\left[{\mathbf{K}}_1^{\left(e\right)}\right]\left\{{\mathbf{\phi }}_{k+1}^{\left(e\right)}\right\}+\left[{\mathbf{K}}_2^{\left(e\right)}\right]\{{\mathbf{\phi }}_{k+1}^{\left(e\right)}{\}}^2+[{\mathbf{K}}_3^{\left(e\right)}\left]\right\{{\mathbf{\phi }}_{k+1}^{\left(e\right)}{\}}^3+\left[{\mathbf{K}}_4^{\left(e\right)}\right]\right),$ (113)

with $\Delta t$ the time step and $k$ the step number.

The starting value for ${\phi }_{k+1}^{\left(0\right)}$ is given by converged values in last step ${\phi }_k$. The updated values in step $k+1$ are given iteratively by

${\phi }_{k+1}^{\left(n+1\right)}={\phi }_{k+1}^{\left(n\right)}-{\left[\frac{dg\left({\phi }_{k+1}^{\left(n\right)}\right)}{d{\phi }_{k+1}^{\left(n\right)}}\right]}^{-1}.g\left({\phi }_{k+1}^{\left(n\right)}\right).$ (114)

7.2 Approximation of motion equation

Based on the motion equation of Eq.(17.i) and given that the second Piola-Kirchhoff stress $S$ is work conjugate to the virtual Green-Lagrange strain $\delta E$, the following expression for the virtual internal work is valid

$\delta W_i={\int }_0^{L_0}{S}^{\left(\alpha \right)}\delta E{A}_{0}dX.$ (115)

The expression for the virtual external work, for any virtual material axial displacement $\delta u$, is

$\delta W_e=P^{\left(\alpha \right)}\left(0\right)\delta u\left(0\right)+P^{\left(\alpha \right)}\left(L_0\right)\delta u\left(L\right)+{\int }_0^{L_p}{\rho }_0{b}_0\delta udX.$ (116)

The virtual work $\delta W_a$ for the inertia load is given by

$\delta W_a=-{\int }_0^{L_0}{\rho }_0\ddot{u}\delta u{A}_{0}dX.$ (117)

The Principle of Virtual Work, given by $\delta W_e+\delta W_a=\delta W_i$ is, then written as

${\int }_0^{L_0}{S}^{\left(\alpha \right)}\delta E{A}_{0}dX=-{\int }_0^{L_0}{\rho }_0\ddot{u}\delta u{A}_{0}dX+{\int }_0^{L_0}{\rho }_0{b}_0\delta udX+P^{\left(\alpha \right)}\left(0\right)\delta u\left(0\right)+P^{\left(\alpha \right)}\left(L_0\right)\delta u\left(L_0\right)=0.$ (118)

The Newmark method employs the following time approximations for velocity and acceleration:

$\begin{array}{c}{v}_{n+1}={\chi }_4\left(u_{n+1}-u_n\right)+{\chi }_5{v}_n+{\chi }_6{\dot{v}}_n,\\ {\dot{v}}_{n+1}={\chi }_1\left(u_{n+1}-u_n\right)-{\chi }_2{v}_n-{\chi }_3{\dot{v}}_n,\end{array}$ (119)

with

$\begin{array}{l}{\chi }_1=\frac{1}{{\beta }_N\Delta t^2},{\chi }_2=\frac{1}{{\beta }_N\Delta t},{\chi }_3=\frac{1}{2{\beta }_N}-1,\\ {\chi }_4=\frac{{\gamma }_N}{{\beta }_N\Delta t},{\chi }_5=1-\frac{{\gamma }_N}{{\beta }_N},{\chi }_6=\left(1-\frac{{\gamma }_N}{2{\beta }_N}\right)\Delta t,\end{array}$ (120)

where ${\gamma }_N$ and ${\beta }_N$ are the Newmark constants.

The expression for the global residue of the motion equation, $r_m$, is obtained as

$\begin{array}{c}{r}_m={\int }_0^{L_0}{\rho }_0\left({\chi }_1\left(u_{n+1}-u_n\right)-{\chi }_2{v}_n-{\chi }_3{\dot{v}}_n\right)\delta u{A}_{0}dX+{\int }_0^{L_0}{S}_{n+1}\delta E_{n+1}{A}_{0}dX\\ -{\int }_0^{L_0}{\rho }_0{b_0}_{n+1}\delta udX-P^{\left(\alpha \right)}\left(0\right)\delta u\left(0\right)-P^{\left(\alpha \right)}\left(L_0\right)\delta u\left(L_0\right).\end{array}$ (121)

The solution of the equation of motion requires the linearization of the PVW for the application of the Newton-Raphson method. The directional derivative of residue $D_{\Delta u}\left[r_m\right]$ in the direction of displacement increment, $\Delta u$, considering that the forces are not dependent on the displacement field, is given by

${\mathrm{D}}_{\mathrm{\Delta }\mathrm{u}}\left[{\mathrm{r}}_{\mathrm{m}}\right]={\int }_0^{{\mathrm{L}}_0}{\mathrm{\rho }}_0{\mathrm{\chi }}_1{\mathrm{u}}_{\mathrm{n}+1}\mathrm{}{\mathrm{A}}_0\mathrm{}\mathrm{d}\mathrm{X}+{\int }_0^{{\mathrm{L}}_0}{\mathrm{D}}_{\mathrm{\Delta }\mathrm{u}}\left[{\mathrm{S}}_{\mathrm{n}+1}\right]\mathrm{\delta }{\mathrm{E}}_{\mathrm{n}+1}{\mathrm{A}}_0\mathrm{}\mathrm{}\mathrm{d}\mathrm{X}+{\int }_0^{{\mathrm{L}}_0}{\mathrm{S}}_{\mathrm{n}+1}{\mathrm{D}}_{\mathrm{\Delta }\mathrm{u}}\left[\mathrm{\delta }{\mathrm{E}}_{\mathrm{n}+1}\right]{\mathrm{A}}_0\mathrm{}\mathrm{}\mathrm{d}\mathrm{X}.$ (122)

We have the following identities for the Green-Lagrange strain:

$\delta E_{n+1}=F_{n+1}\frac{d\left(\delta u\right)}{dX},$ (123)

$D_{\Delta u}\left[\delta E_{n+1}\right]=\frac{d\left(\delta u\right)}{dX}\frac{d\left(\Delta u\right)}{dX}$,(124)

$D_{\Delta u}\left[S_{n+1}\right]={\mathcal{C}}_{n+1}{D}_{\Delta u}\left[E_{n+1}\right],$ (125)

$D_{\Delta u}\left[E_{n+1}\right]=F_{n+1}\frac{d\left(\Delta u\right)}{dX}.$ (126)

$\mathcal{C}$ is the constitutive function, also called tangent modulus, given by

${\mathcal{C}}_{n+1}=\frac{\partial S_{n+1}}{\partial E_{n+1}}.$ (127)

The residual expression, therefore, is rewritten as

$\begin{array}{l}{D}_{\Delta u}\left[r_m\right]={\int }_0^{L_0}{\rho }_0{\chi }_1{u}_{n+1}\delta u{A}_{0}dX+{\int }_0^{L_0}{F}_{n+1}\frac{d\left(\delta u\right)}{dX}{\mathcal{C}}_{n+1}{F}_{n+1}\frac{d\left(\Delta u\right)}{dX}{A}_{0}dX\mathrm{}+{\int }_0^{L_0}{S}_{n+1}\frac{d\left(\delta u\right)}{dX}\frac{d\left(\Delta u\right)}{dX}{A}_{0}dX.\\ \end{array}$ (128)

The final expression after the linearisation of the PVW is

$\begin{array}{l}{\int }_0^{L_0}{\rho }_0{\chi }_1{u}_{n+1}\delta u{A}_{0}dX+{\int }_0^{L_0}{F}_{n+1}\frac{d\left(\delta u\right)}{dX}{\mathcal{C}}_{n+1}{F}_{n+1}\frac{d\left(\Delta u\right)}{dX}{A}_{0}dX\\ +{\int }_0^{L_0}{S}_{n+1}\frac{d\left(\delta u\right)}{dX}\frac{d\left(\Delta u\right)}{dX}{A}_{0}dX={\int }_0^{L_0}{\rho }_0\left({\chi }_1\left(u_{n+1}-u_n\right)-{\chi }_2{v}_n-{\chi }_3{\dot{v}}_n\right)\delta u{A}_{0}dX\\ +{\int }_0^{L_0}{S}_{n+1}{F}_{n+1}\frac{d\left(\delta u\right)}{dX}{A}_{0}dX-{\int }_0^{L_0}{\rho }_0{b_0}_{n+1}\delta udX-P^{\left(\alpha \right)}\left(0\right)\delta u\left(0\right)-P^{\left(\alpha \right)}\left(L_0\right)\delta u\left(L_0\right).\end{array}$ (129)

Analogously to the damage equations, the approximation of the material coordinate, in an element $e$, using the local normalized coordinate system $\xi$ is given by

$X\left(\xi \right)=\sum _{i=1}^{P+1}{N}_i^e\left(\xi \right)X_i=\left[{\mathbf{N}}^{\left(e\right)}\left(\xi \right)\right]{\{X}^{\left(e\right)}\},$ (130)

where $X_i$ are the nodal coordinates of the initial mesh, $N_i^{\left(e\right)}\left(\xi \right)$ are the element shape functions, $P$ is the polynomial order and $P+1$ is the number of nodes of the element.

For an isoparametric element, the material axial displacement and strain are interpolated as

$u\left(\xi \right)=\sum _{i=1}^{P+1}{N}_i^{\left(e\right)}\left(\xi \right)u_i=\left[{\mathbf{N}}^{\left(e\right)}\left(\xi \right)\right]{\{u}^{\left(e\right)}\},$ (131)

$\frac{du\left(\xi \right)}{dX}=\sum _{i=1}^{P+1}{N}_{i,\xi }^{\left(e\right)}\left(\xi \right)u_i=\sum _{i=1}^{P+1}{B}_i^{\left(e\right)}\left(\xi \right)u_i=\left[{\mathbf{B}}^e\left(\xi \right)\right]\left\{u_i^{\left(e\right)}\right\},$ (132)

where $B_i^{\left(e\right)}\left(\xi \right)$ are the global derivatives of the element shape functions.

The virtual axial displacement and the infinitesimal normal strain are obtained deriving the previous expressions to each nodal coordinate $u_i^{\left(e\right)}$, as follows

$\begin{array}{l}\delta u\left(\xi \right)=\left[{\mathbf{N}}^{\left(e\right)}\right(\xi ){]}^T,\\ \frac{d\left(\delta u\left(\xi \right)\right)}{dX}=\left[{\mathbf{B}}^{\left(e\right)}\right(\xi ){]}^T,\end{array}$ (133)

Finally, the displacement increment $\Delta u$ and its derivative are interpolated as

$\begin{array}{l}\Delta u\left(\xi \right)=\left[{\mathbf{N}}^{\left(e\right)}\left(\xi \right)\right]\left\{\Delta u^{\left(e\right)}\right\},\\ \frac{d\left(\Delta u\left(\xi \right)\right)}{dX}=\left[{\mathbf{B}}^e\left(\xi \right)\right]\left\{\Delta u_i^{\left(e\right)}\right\},\end{array}$ (134)

where $\left\{\Delta u_i^{\left(e\right)}\right\}$ is the vector of nodal displacement increments in the material configuration.

Replacing the approximations in the Eq. (129), we have the residue and Jacobian for each element given, respectively, by

$\begin{array}{l}{r}_m={\int }_{L_{e0}}{\rho }_0\left[{\mathbf{N}}^{\left(e\right)}{]}^T\right[{\mathbf{N}}^{\left(e\right)}]A_{0}dX\left({\chi }_1\left(u_{n+1}^{\left(e\right)}-u_n^{\left(e\right)}\right)-{\chi }_2{v}_n^{\left(e\right)}-{\chi }_3{\dot{v}}_n^{\left(e\right)}\right)+\\ {\int }_{L_{e0}}[{\mathbf{B}}^{\left(e\right)}{]}^T{F}_{n+1}^{\left(e\right)}{S}_{n+1}^{\left(e\right)}{A}_0\mathrm{}dX-{\int }_{L_{e0}}{\rho }_0\left[{\mathbf{N}}^{\left(e\right)}{]}^T\right[{\mathbf{N}}^{\left(e\right)}]{b_0}_{n+1}^{\left(e\right)}dX-\{P_0^{\left(e\right)}{P}_L^{\left(e\right)}{\}}^T,\end{array}$ (135)

$\begin{array}{c}{J}_m^{\left(e\right)}={\chi }_1{\int }_{L_{e0}}{\rho }_0\left[{\mathbf{N}}^{\left(e\right)}{]}^T\right[{\mathbf{N}}^{\left(e\right)}]dX+{\int }_{L_{e0}}\left[\left[{\mathbf{B}}^{\left(e\right)}{]}^T{S}_{n+1}^{\left(e\right)}\right[{\mathbf{B}}^{\left(e\right)}]+[{\mathbf{B}}^{\left(e\right)}{]}^T{F}_{n+1}^{\left(e\right)}{\mathbb{C}}_{n+1}^{\left(e\right)}{F}_{n+1}^{\left(e\right)}\left[{\mathbf{B}}^{\left(e\right)}\right]\right]dX.\end{array}$ (136)

Again, the global equations are obtained from the assembling procedure of the element operators.

8 RESULTS

8.1 Parameter fitting

Defining material parameters for complex rheological models requires an extensive effort, since most of the terms are not obtained directly from experiments. For this reason, it is necessary to apply optimization techniques, in order to adjust the material model response to the results of stress-strain curves available in the literature.

The parameter fitting requires hundreds of runs of simulation. For this reason and considering that experimental data gives stress-strain values, it considers the constitutive model applied to a single integration point, in order to run the model faster. Moreover, since the damage occurs abruptly and just before the fracture, it is not included in this fitting procedure.

An genetic algorithm, based on the principle of natural selection, was used for parameter fitting. This technique, as shown in Fig. 6, allows the evolution of an initial population of by several individuals to a final population which maximizes (or minimizes) a given cost function after various generations.

Figure 6
Steps of the inverse problem solved by genetic algorithm.

For the current problem, each individual is the simulation of the constitutive model for a given initial state of parameters. After each generation, the cost function is evaluated comparing the result of the simulation to the experimental stress-strain curve of the PTFE. The cost function to be minimized is given by

For the current problem, each individual is the simulation of the constitutive model for a given initial state of parameters. After each generation, the cost function is evaluated comparing the result of the simulation to the experimental stress-strain curve of the PTFE. The cost function to be minimized is given by

$Cost=\sum _{i=0}^{Steps}\sqrt{({T_{model}}_{\left(i\right)}-{T_{experiment}}_{\left(i\right)}{)}^2},$ (137)

with $T$ being the stresses evaluated at each time step of the simulation and from experiment.

The starting parameters were taken from Marin et al. (2007) and Srivastava et al. (2010). The experimental results for the stress-strain were obtained from Rae and Brown (2005). In this reference, the values of the Young modulus for several temperatures, obtained directly from the experiments, were used in the model without the need to be fitted with the genetic algorithm.

The genetic algorithm used in this work is called eaSimple, from the Python library DEAP (Distributed evolutionary algorithms in Python). The population size was set to 200, evaluated in 50 generations. Others internal parameters required for this algorithm are the crossover, individual mutation and gene mutation, whose probability values were set to 90%, 20% and 5%, respectively.

The optimization process is used for each one of the 7 values of temperature available. Consequently, an optimum solution is found for each temperature and 7 completely different sets of parameters are determined. Due to this, a second and tedious part of analysis is needed, aiming to keep the largest number of parameters constant through the entire range of temperatures (and, if possible, eliminating unnecessary coefficients). After this manual stage, certain parameters still cannot be kept constant.

The parameters with constant values for all temperatures are: $E_{gl}$ = 650 MPa, $E_r$ = 550 MPa, $M_{gl}$ = 31 MPa.K${}^{-1}$, $M_r$ = 4.3 MPa.K${}^{-1}$, $\Delta$ = 11 K, $h_b$ = 10.9, ${\alpha }_p$ = 0.3, ${\nu }_0^{\left(1\right)}$ = $8.20\times {10}^{12}$ s${}^{-1}$, $d$ = $1.50\times {10}^{-2}$K${}^{-1}$, ${\zeta }_{gl}$ = 0.17, $V$ = $8.00\times {10}^{-28}$ m${}^3$, $m^{\left(1\right)}$ = 0.15, $\gamma$ = 100, $r$ = $4.50\times {10}^4$, $S_r$ = 0.1 MPa, $I_m^{\left(2\right)}$ = 24, ${\nu }^{\left(2\right)}$ = $3.70\times {10}^{-4}$ s${}^{-1}$, $S_{gl}^{\left(2\right)}$ = 25 MPa, $S_r^{\left(2\right)}$ = $1.10\times {10}^4$ Pa, $m^{\left(2\right)}$ = 0.3, $I_m^{\left(3\right)}$ = 24, ${\nu }_0^{\left(3\right)}$ = $9.00\times {10}^{-3}$ s${}^{-1}$, $m^{\left(3\right)}$ = 0.18, $S_g^{\left(3\right)}$ = 6 MPa, $h_{3g}$ = $9.50\times {10}^2$ MPa, $Y$ = 0 K${}^{-1}$ and $Z$ = 0 K${}^{-1}$.

The parameters listed in Table 1 are dependent on temperature.

Table 1
Final values of parameters of the non-isothermal material model.

Parameter	Temperature
	-50°C	-15°C	15°C	25°C	50°C	100°C	150°C
$Q_{gl}\left[J\right]$	$1.49\times {10}^{-19}$	$1.41\times {10}^{-19}$	$1.30\times {10}^{-19}$	$1.22\times {10}^{-19}$	$1.22\times {10}^{-19}$	$1.22\times {10}^{-19}$	$1.22\times {10}^{-19}$
$Q_r\left[J\right]$	$1.49\times {10}^{-19}$	$1.41\times {10}^{-19}$	$1.30\times {10}^{-19}$	$1.22\times {10}^{-19}$	$1.22\times {10}^{-19}$	$1.22\times {10}^{-19}$	$1.22\times {10}^{-19}$
$S_{gl}\left[Pa\right]$	$9.00\times {10}^7$	$8.20\times {10}^7$	$1.00\times {10}^7$	$1.00\times {10}^5$	$1.00\times {10}^5$	$1.00\times {10}^5$	$1.00\times {10}^5$

The parameters ${\mu }_g^{\left(2\right)}$ and ${\mu }^{\left(3\right)}$ not listed in Table 1 are also temperature dependent, since their values are given by

$\begin{array}{c}{\mu }_g^{\left(2\right)}={\mu }^{\left(3\right)}=W-Q\text{tanh}\left[\frac{1}{0.5(T_2-T_1)}\left(\theta -\frac{T_2+T_1}{2}\right)\right]\text{exp}\left[a\left|\left(\theta -\frac{T_2+T_1}{2}\right)\right|\right],\\ \end{array}$ (138)

where W = 3×10⁶, Q = 3.15×10⁵, a = 0.001, T₁ = 19^°C (transition from phase II to IV), and T₁ = 30^°C (transition from phase IV to I).

Eq.(138) plotted in Fig.(7) looks reproducing typical shape of storage modulus ($E\text{'}$) of PTFE, with different behaviour due to phase change, between 19°C and 30°C, which represent transitions in crystallyne phases.

Figure 7
Behaviour of ${\mathrm{\mu }}_{\mathrm{g}}^{\left(2\right)}$ and ${\mathrm{\mu }}^{\left(3\right)}$, inspired in the storage modulus from Tóth, Baets and Szebényi (2020).

8.2 Constitutive model without damage

The first simulation uses the previously adjusted parameters in a single integration point of the constitutive model of a non-damaged PTFE model, under temperatures from -50°C to 150°C, and the resulting stress-strain curves are shown in Fig. 8. Strain rates for model and experiment are both set to $5\times {10}^{-3}{s}^{-1}$.

Figure 8
Adjusted stress-strain curves for the 3-branches constitutive model for the PTFE at several temperatures.

Experimental curves were obtained from Rae and Brown (2005). There is a good fit for all temperatures, except for 25°C, where there is overestimation of results for mid-strains, representing a higher stiffness of material in this region.

It is important to highlight that the last strain value considered in the simulation for each temperature is set in accordance with the end values from experimental results.

Fig. 9 shows the influence of the branches according to the temperature. In negative temperatures, the main influence is given by branch one. After increasing from 25°C, there is higher influence of branch 3 as noticed from Fig. 9.

Figure 9
Influence of stresses from each branch for PTFE at several temperatures.

8.3 Relaxation test model

In order to verify the behaviour of PTFE under relaxation cycle, the 3-branch constitutive model is now simulated under the strain loading profiles shown in Fig. 10. Taking a proportional value to the maximum strain in each temperature under analysis as reference, a sequence of five steps is considered, composed by alternate cycles of strain loading followed by steady strain.

Figure 10
Alternate cycles of strain loading followed by steady strain for PTFE.

The results for these load profiles are presented in Fig. 11. As can be noticed, after a period of steady strain, the stresses present relaxation and, after a loading restart, the stresses resume their evolution from the point where the relaxation started. This behaviour is in accordance with experimental results from Khan and Zhang (2001).

Figure 11
Relaxation results in terms of strain for different temperatures.

8.4 Fracture damage evolution

After calibrating the constitutive model, uniaxial simulation of damage using the finite element model is performed.

According to Rae and Brown (2005), it is noticed an increase of the failure strain from -50°C to 25°C C. From 25°C on, a constant failure strain is reached, with a small decrease near 150°C. On the other hand, the failure stresses reach a limit value in 25°C, where crystalline structure of PTFE is in phase IV.

Considering this profile of fracture depending on temperature, the model should present an abrupt failure in accordance with the stress-strain curves.

To simulate the fracture of the bar under different temperature conditions, the constitutive model described in Section 5 is used, in order to capture all the differences in microstructural behavior. In addition, parameter $\tilde{\lambda }$ representing fracture onset is adjusted in order to represent different responses to failure in all temperature range.

Results of the stress-strain curves for temperatures from -50°C to 150°C are shown in Fig.12. All strain rates are set to 5 × 10⁻³ to be compliant with experiments described by Rae and Brown (2005).

Figure 12
Stress-strain curve for PTFE, including fracture at different temperatures.

Another influence on polymer response is the strain-rate. For PTFE, however, this variable has a lower influence when compared to temperature, with an almost insensibility to different strain rates. For the evaluation of the proposed adjusted model, a simulation was run for different strain rates at 25°C. The simulation results are shown in Fig.13. The results are quantitatively similar to experiments described in Nunes, Dias and Mattos (2011) and Rae and Brown (2005), mainly for lower strain rates.

Figure 13
Stress-strain curves from model response for of different strain rates, T = 25 C for all cases.

9 CONSIDERATIONS ABOUT THE DEGRADATION FUNCTION $D\left(\phi \right)$

The degradation function $D\left(\phi \right)$ describes the influence of the phase-field on the possible decreasing of the elastic response of the material and, therefore, its choice is very important.

Since degradation involves changes in the microstructure, $D\left(\phi \right)$ may depend on the material being considered because different materials may have several different degradation mechanisms; for instance, the appearance of voids and microcracks are important degradation mechanisms in metals and ceramics; degradation in polymers are more complex and may occur due to sliding, tangling or breaking of polymeric chains and also due to crazing.

Also, these several degradation mechanisms may be described in different details according to the material model being developed. This means that the effects in the degradation of the elastic response of the material, described by $D\left(\phi \right)$, may also depend on how the the model being developed describes the degradation mechanisms presented by the material being considered.

In most of the cases, the several different degradation mechanisms are not individually described. The phase-field considered in the model is viewed as the mean damaged state of the material, due to changes in the material microstructure resulting maybe from several degradation mechanism; the corresponding degradation function must then give the mean effect in the corresponding elastic properties. This is so because it is in general difficult to describe in detail and in a realist way each of the degradation mechanisms and their cross influence; moreover, when this is done, several material parameters, which are difficult to measure or estimated, are introduced in the material model making it more difficult to validate by comparison with experiments.

Models that view the phase-field as a mean damage state of the material usually obtain the corresponding degradation function in a phenomenological way; that is, $D\left(\phi \right)$ is chosen to fit to the observed material response. The importance of the choice of $D\left(\phi \right)$ to obtain correct results has been discussed along the years, as one can see, for instance, in Borden et al. (2016), Haveroth et al. (2020) and Costa-Haveroth et al. (2022) and the references cited there in. The model of the present work is of this kind, and so the expression given in (99) was chosen to fit the experimental PTFE material response.

However, the previous considerations may raise the question whether a simpler expression, like the more standard $D\left(\phi \right)={(1-\phi )}^2$ could be used also to fit the experimental results if the material parameters were suitably chosen. The answer is no, as we show in the following numerical simulations varying the material parameters and using $D\left(\phi \right)={(1-\phi )}^2$ instead of (99).

Figure 14 plots the stress-strain curves considering, respectively, the variation of one of the parameters $c,\gamma ,\frac{1}{\tilde{\lambda }},g_c$ of the damage equation. The values of these parameters when they not changed are $c={10}^{-2},\gamma ={10}^{-2},\frac{1}{\tilde{\lambda }}=\frac{1}{{10}^5},g_c=1000.$

Figure 14
Stress-strain curves of the 1D finite element model at -15 ^oC using D(ϕ) = (1− $\phi$)². For each plot above, only one of the parameters of the damage equation $(c,\gamma ,\frac{1}{\tilde{\lambda }},g_c)$ is changed, to verify its impact in the response.

As we can see from the previous numerical simulations, by using $D\left(\phi \right)={(1-\phi )}^2$ in our model, we did not obtain an abrupt failure. We think that is because $D\left(\phi \right)={(1-\phi )}^2$ degrades too much the elastic response at rather small values of $\phi$; so, the material softens and the failure becomes slower. Thus, to obtain the correct failure behavior, $D\left(\phi \right)$ should have a slow decay for rather small values of $\phi$ as in Fig. 5.

A final remark is that Narayan and Anand (2021) were able to obtain abrupt failure by using the degradation function $D\left(\phi \right)={(1-\phi )}^2$ because they include an explict modelling of the degradation mechanism of crazing which is very rele-vant to the material they are considering. In the present work, the degradation function (99) implicitly includes the effects of all degradation mechanisms.

10 CONSIDERATIONS ABOUT THE DAMAGE ENERGY $\mathcal{I}(\phi ,\nabla \phi )$

One may wonder why the expression (100) used in the present work for the damage energy $\mathcal{I}(\phi ,\nabla \phi )$ is quite different from the classical AT1 and AT2 phase-field formulations. Before we answer this, we must briefly explain the differences between the two phase-field methodologies.

Both AT1 and AT2 formulations are derived as follows: starting by assuming that the phase-field is a diffusive approximation (diffusive representation) of a discontinuous crack, a specific form for such approximation is chosen around a crack. Next, a suitable equation, having that approximation as solution is considered and it is assumed that the phase-field variable must satisfy a variant of that chosen equation.

For instance, in the AT1 one dimension formulation, the chosen approximation around a crack located at the origin is $d\left(x\right)={\left(\right|x|/(2l_c)-1)}^2$, for $\left|x\right|\le 2l_c$ and $d\left(x\right)=0$ for $\left|x\right|>2l_c$, with $l_c$ a given length scale parameter, while in the AT2 formulation the corresponding approximation is $d\left(x\right)=e^{-\left|x\right|/l_c}$. For a brief description of these formulations, see for instance, Molnár et al. (2022). In those cases, it is then observed that there is a “damage energy” expression associated to the proposed equation in the sense that it is the corresponding Euler-Lagrange equation associated to the minimum of that damage energy.

This approach to the derivation of the phase-field equation is known as the diffusive approximation methodology.

However, there are several difficulties with this approach: (1) the diffusive profile around the crack, which gives the damage distribution around it, is a priori chosen, which raises the question about why one approximation is preferable than others; (2) it is the “damage energy” expression that depends on the chosen form of the approximation, while the opposite is physically expected; that is, the basic physical laws and the expressions of the several energies, including a meaningful damage energy, should imply the damage profile around the crack (the diffusive approximation); (3) the view that the phase-field equation is the Euler-Lagrange equation associated to the minimum of the damage energy usually leads to a stationary (elliptic) equation, which forces the use of the notion of the pseudo-time to obtain the damage evolution; (4) it is doubtful whether models obtained using this approach are really thermodynamically consistent in the strict sense; in particular, this is certainly so in the cases where the temperature is involved.

To avoid such difficulties, the entropy methodology was developed; it reverses the arguments of the diffusive approximation methodology. That is, instead of trying to obtain an a priori equation for the phase-field variable, it is treated at the same level as the other physical state variables; that is, it is not just seen as an approximation, but a state variable that must be coherent with the basic physical laws and the thermodynamic restrictions.

Then, to obtain the model equations, one proceeds as follows: the principle of virtual powers is used to derive the basic dynamical equations for all the state variables; next, assuming that the expressions for the Helmholtz free energy and the pseudo-potential of dissipation of the system are given, the Coleman-Noll approach is used to obtain the necessary constitutive equations in terms of the free energy and the pseudo-potential of dissipation to guarantee that the entropy inequality (dissipation inequality) is actually satisfied. The derived basic dynamical equations, the balance of energy equation and the derived constitutive equations give then the general form of the system to be solved.

In any particular case, one then proposes reasonable specific expressions for all the involved energies, including the damage energy, and also for the pseudo-potentials of dissipation to specialize those equations and finally obtain the mathematical model associated to that case. This includes the equation for the phase-field, which is now a dynamical equation, depending thus on the true time, not on a pseudo-time; then, it can be used to obtain the time evolution of the damage. Moreover, the damage profile around cracks is now not a priori chosen, but something that is derived from the involved energies, as physically expected. Also, the coupling with temperature appearing in the model is thermodynamically sound.

This approach is much more general and physically sound than the diffusive approximation methodology and does not have the previously mentioned difficulties. The entropy approach is the one used to build the framework proposed in Boldrini et al. (2016) and in the present work.

We remark that, for sure, the quality of the final mathematical model derived using the entropy approach depends on the quality of the chosen free-energy and pseudo-potential of dissipation, but this is also true for the other approaches.

About the damage energy proposed in this work, we observe that, since our approach does not a priori choose the form of the diffusive approximation, which fixes the form of the associated damage energy, we now have more freedom to choose its expression. For this, we observe that the damage energy may accumulate in the bulk of the material or on the faces of cracks. In terms of the phase-field, the bulk damage energy depends on the phase-field itself, while the damage energy accumulated on the faces of cracks depend on the gradient of the phase-field; that is, on the transition layers of the phase-field.

This explains why the expression (100) for the total damage energy $\mathcal{I}(\phi ,\nabla \phi )$ has two parts: the bulk damage energy whose expression is $g_c\frac{c}{\gamma }\mathcal{H}\left(\phi \right)$ and the damage energy at the faces of the smeared cracks: $g_c\frac{\gamma }{2}{\left|\frac{d\phi }{dX}\right|}^2$. Here $c$ is a material parameter that weights how much of the damage energy is accumulated in the bulk instead of on the faces of the cracks. This means that, when it is small, most of the damage energy tends to accumulate on the faces of the cracks, which implies that $c$ has a rather large influence on the velocity of the crack propagation. The value of the material parameter $c$ used in this work was the most adequate to fit the experimental results.

We finally observe that, since the diffusive approximation methodology a priori chooses the form of the diffusive approximation, fixing thus the form of the associated damage energy, it also a priori fixes the relation between the bulk and crack faces damage energies.

11 CONCLUSIONS

In this work, it was presented an one-dimensional finite element model along with a phase-field approach for describing fracture behaviour of polymers, specifically for the polytetrafluoroethylene (PTFE). The model was able to accurately capture the material response over a wide temperature range (-50°C to 150°C), with good agreement between simulation and experimental data including the abrupt behaviour of fracture evolution.

One contribution of this work is the parameter fitting that was performed to model the PTFE through a genetic algorithm. This is particularly important given the amount of material parameters. In addition, the model was able to reproduce relaxation test and the influence of strain rate in stress-strain test, with good results for low and medium strains.

The implementation of this model provides a tool for understanding and predicting the fracture behaviour of PTFE and potentially other polymeric materials. The phase-field approach used in this work can be extended to more dimensions, enabling more accurate and efficient modelling of fracture in a variety of applications.

Overall, this work demonstrates the potential of the phase-field approach for modelling fracture in polymer materials and can help to inform the design of more resilient and fracture-resistant polymer materials, with potential applications in a range of industries.

References

• ALLEN, S. M.; CAHN, J. W. Ground state structures in ordered binary alloys with second neighbor interactions. Acta Metallurgica, Elsevier, v. 20, n. 3, p. 423-433, 1972.
• AMENDOLA, G.; FABRIZIO, M.; GOLDEN, J. Thermomechanics of damage and fatigue by a phase-field model. Journal of Thermal Stresses, Taylor & Francis, v. 39, n. 5, p. 487-499, 2016.
• AMES, N. M. et al. A thermo-mechanically coupled theory for large deformations of amorphous polymers. part ii: Applications. International Journal of Plasticity, Elsevier, v. 25, n. 8, p. 1495-1539, 2009.
• ANAND, L.; AMES, N. On modeling the micro-indentation response of an amorphous polymer. International Journal of Plasticity, Elsevier, v. 22, n. 6, p. 1123-1170, 2006.
• ANAND, L. et al. A thermo-mechanically coupled theory for large deformations of amorphous polymers. part i: Formulation. International Journal of Plasticity, Elsevier, v. 25, n. 8, p. 1474-1494, 2009.
• ANAND, L.; GURTIN, M. E. A theory of amorphous solids undergoing large deformations, with application to polymeric glasses. International Journal of Solids and Structures, Elsevier, v. 40, n. 6, p. 1465-1487, 2003.
• ANDERSON, T. L. Fracture mechanics: fundamentals and applications. [S.l.]: CRC press, 2017.
• ARRUDA, E. M.; BOYCE, M. C.; JAYACHANDRAN, R. Effects of strain rate, temperature and thermomechanical coupling on the finite strain deformation of glassy polymers. Mechanics of Materials, Elsevier, v. 19, n. 2-3, p. 193-212, 1995.
• BERGSTROM, J. S. Mechanics of solid polymers: theory and computational modeling. [S.l.]: William Andrew, 2015.
• BERGSTRÖM, J. S.; BOYCE, M. C. Constitutive modeling of the large strain time-dependent¨ behavior of elastomers. Journal of the Mechanics and Physics of Solids, Elsevier, v. 46, n. 5, p. 931-954, 1998.
• BOLDRINI, J. et al. A non-isothermal thermodynamically consistent phase-field framework for structural damage and fatigue. Computer Methods in Applied Mechanics and Engineering, Elsevier, v. 312, p. 395-427, 2016.
• BOLDRINI, J. L. Phase-field: A methodology to model complex material behavior. In: Advances in Mathematics and Applications [S.l.]: Springer, 2018. p. 67-103.
• BORDEN, M. J. et al. A phase-field formulation for fracture in ductile materials: Finite deformation balance law derivation, plastic degradation, and stress triaxiality effects. Computer Methods in Applied Mechanics and Engineering, Elsevier, v. 312, p. 130-166, 2016.
• BOYCE, M. C.; PARKS, D. M.; ARGON, A. S. Large inelastic deformation of glassy polymers. part i: rate dependent constitutive model. Mechanics of Materials, Elsevier, v. 7, n. 1, p. 15-33, 1988.
• BOYCE, M. C.; PARKS, D. M.; ARGON, A. S. Plastic flow in oriented glassy polymers. International Journal of Plasticity, Elsevier, v. 5, n. 6, p. 593-615, 1989.
• BROWN, E. N.; DATTELBAUM, D. M. The role of crystalline phase on fracture and microstructure evolution of polytetrafluoroethylene (ptfe). Polymer, Elsevier, v. 46, n. 9, p. 3056-3068, 2005.
• BROWN, E. N. et al. Phase dependent fracture and damage evolution of polytetrafluoroethylene (ptfe). 2004.
• BUCKLEY, C. et al. Deformation of thermosetting resins at impact rates of strain. part 2: constitutive model with rejuvenation. Journal of the Mechanics and Physics of Solids, Elsevier, v. 52, n. 10, p. 2355-2377, 2004.
• CAHN, J. W.; HILLIARD, J. E. Free energy of a nonuniform system. i. interfacial free energy. The Journal of chemical physics, American Institute of Physics, v. 28, n. 2, p. 258-267, 1958.
• CALLEJA, G. et al. Where is the glass transition temperature of poly (tetrafluoroethylene)? a new approach by dynamic rheometry and mechanical tests. European Polymer Journal, Elsevier, v. 49, n. 8, p. 2214-2222, 2013.
• COLAK, O.; CAKIR, Y. Material model parameter estimation with genetic algorithm optimization method and modeling of strain and temperature dependent behavior of epoxy resin with cooperative-vbo model. Mechanics of Materials, Elsevier, v. 135, p. 57-66, 2019.
• COSTA-HAVEROTH, T. C. da et al. A damage phase-field model for fractional viscoelastic materials in finite strain. Computational Mechanics, Springer, v. 69, n. 6, p. 1365-1393, 2022.
• DU, Q.; LI, M.; LIU, C. Analysis of a phase-field navier-stokes vesicle-fluid interaction model. Discrete & Continuous Dynamical Systems-B, American Institute of Mathematical Sciences, v. 8, n. 3, p. 539, 2007.
• DUPAIX, R. B.; BOYCE, M. C. Constitutive modeling of the finite strain behavior of amorphous polymers in and above the glass transition. Mechanics of Materials, Elsevier, v. 39, n. 1, p. 39-52, 2007.
• FABRIZIO, M.; GIORGI, C.; MORRO, A. A thermodynamic approach to non-isothermal phase-field evolution in continuum physics. Physica D: Nonlinear Phenomena, Elsevier, v. 214, n. 2, p. 144-156, 2006.
• FERNÁNDEZ, J. R. et al. A genetic algorithm for the characterization of hyperelastic´ materials. Applied Mathematics and Computation, Elsevier, v. 329, p. 239-250, 2018.
• FRANCFORT, G. A.; MARIGO, J.-J. Revisiting brittle fracture as an energy minimization problem. Journal of the Mechanics and Physics of Solids, Elsevier, v. 46, n. 8, p. 1319-1342, 1998.
• FRÉMOND, M.; NEDJAR, B. Damage, gradient of damage and principle of virtual power. International Journal of Solids and Structures, Elsevier, v. 33, n. 8, p. 1083-1103, 1996.
• FRÉMOND, M.; SHITIKOVA, M. Non-smooth thermomechanics. Applied Mechanics Review, v. 55, n. 5, p. B99-B100, 2002.
• GAMBONI, O. C. et al. On the formation of defects induced by air trapping during cold pressing of PTFE powder. Polymer, Elsevier, v. 82, p. 75-86, 2016.
• GOVAERT, L.; TIMMERMANS, P.; BREKELMANS, W. The influence of intrinsic strain softening on strain localization in polycarbonate: modeling and experimental validation. Journal of Engineering Materials and Technology, v. 122, n. 2, p. 177-185, 2000.
• GRIFFITH, A. A. The phenomena of rupture and flow in solids. Philosophical transactions of the Royal Society of London. Series A, containing papers of a mathematical or physical character, The Royal Society London, v. 221, n. 582-593, p. 163-198, 1921.
• HASAN, O.; BOYCE, M. A constitutive model for the nonlinear viscoelastic viscoplastic behavior of glassy polymers. Polymer Engineering & Science, Wiley Online Library, v. 35, n. 4, p. 331-344, 1995.
• HAVEROTH, G. et al. A non-isothermal thermodynamically consistent phase field model for damage, fracture and fatigue evolutions in elasto-plastic materials. Computer Methods in Applied Mechanics and Engineering, Elsevier, v. 364, p. 112962, 2020.
• HAWARD, R.; THACKRAY, G. . The use of a mathematical model to describe isothermal stress-strain curves in glassy thermoplastics. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences, The Royal Society London, v. 302, n. 1471, p. 453-472, 1968.
• KANNINEN, M. F.; MCEVILY, A.; POPELAR, C. H. Advanced fracture mechanics. The American Society of Mechanical Engineers (ASME), 1986.
• KHAN, A.; ZHANG, H. Finite deformation of a polymer: experiments and modeling. International Journal of Plasticity, Elsevier, v. 17, n. 9, p. 1167-1188, 2001.
• KLETSCHKOWSKI, T.; SCHOMBURG, U.; BERTRAM, A. Endochronic viscoplastic material models for filled ptfe. Mechanics of Materials, Elsevier, v. 34, n. 12, p. 795-808, 2002.
• LEE, E. H. Elastic-plastic deformation at finite strains. Journal of Applied Mechanics, v. 36, p. 1-6, 1969.
• LEMAITRE, J.; DESMORAT, R. Engineering damage mechanics: ductile, creep, fatigue and brittle failures. [S.l.]: Springer Science & Business Media, 2006.
• MARIN, E. et al. On developing a viscoelastic-viscoplastic model for polymeric materials. CAVS Internal Report, p. 2007, 2007.
• MIEHE, C.; SCHAENZEL, L.-M.; ULMER, H. 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. Computer Methods in Applied Mechanics and Engineering, Elsevier, v. 294, p. 449-485, 2015.
• MOLNÁR, G. et al. Thermodynamically consistent linear-gradient damage model in Abaqus. Engineering Fracture Mechanics, Elsevier, v. 266, p. 108390, 2022.
• NARAYAN, S.; ANAND, L. Fracture of amorphous polymers: A gradient-damage theory. Journal of the Mechanics and Physics of Solids, Elsevier, v. 146, p. 104-164, 2021.
• NUNES, L.; DIAS, F.; MATTOS, H. da C. Mechanical behavior of polytetrafluoroethylene in tensile loading under different strain rates. Polymer Testing, Elsevier, v. 30, n. 7, p. 791-796, 2011.
• PARKS, D.; ARGON, A.; BAGEPALLI, B. Large elastic-plastic deformation of glassy polymers. Mechanics of Materials, v. 7, n. 1, p. 15-33, 1984.
• RAE, P.; BROWN, E. The properties of poly (tetrafluoroethylene)(ptfe) in tension. Polymer, Elsevier, v. 46, n. 19, p. 8128-8140, 2005.
• SCHÄNZEL, L.-M. Phase-field modeling of fracture in rubbery and glassy polymers at finite thermo-viscoelastic deformations Thesis (PhD) — Stuttgart University, 2015.
• SILVA, S. P. L. D. d. Phase-field models for tumor growth uin the avascular phase Thesis (PhD) — Coimbra University, 2009.
• SRIVASTAVA, V. et al. A thermo-mechanically-coupled large-deformation theory for amorphous polymers in a temperature range which spans their glass transition. International Journal of Plasticity, Elsevier, v. 26, n. 8, p. 1138-1182, 2010.
• SUN, P.; XU, J.; ZHANG, L. Full eulerian finite element method of a phase-field model for fluid-structure interaction problem. Computers & Fluids, Elsevier, v. 90, p. 1-8, 2014.
• TERVOORT, T. et al. A constitutive equation for the elasto-viscoplastic deformation of glassy polymers. Mechanics of Time-Dependent Materials, Springer, v. 1, n. 3, p. 269-291, 1997.
• TÓTH, L. F.; BAETS, P. D.; SZEBÉNYI, G. Thermal, viscoelastic, mechanical and´ wear behaviour of nanoparticle filled polytetrafluoroethylene: A comparison. Polymers, Multidisciplinary Digital Publishing Institute, v. 12, n. 9, p. 1940, 2020.
• TRELOAR, L. G. The physics of rubber elasticity. OUP Oxford, 1975.
• TRUESDELL, C.; NOLL, W. Encyclopedia of physics. The non-linear field theories of mechanics, v. 3, p. 3, 1965.
• WU, P.; GIESSEN, E. V. D. On improved network models for rubber elasticity and their applications to orientation hardening in glassy polymers. Journal of the Mechanics and Physics of Solids, Elsevier, v. 41, n. 3, p. 427-456, 1993.

Edited by

Editor: Marco L. Bittencourt and Josué Labaki

Publication Dates

• Publication in this collection
  18 Sept 2023
• Date of issue
  2023

History

• Received
  27 Feb 2023
• Reviewed
  24 June 2023
• Accepted
  13 July 2023
• Published
  25 July 2023