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. Graphical Abstract


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.

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 rubberlike 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 largestrain 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 −4 s −1 to 1×10 −1 s −1 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 andPopelar (1986) andAnderson (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.

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.
Latin American Journal of Solids and Structures, 2023, 20(6), e504 4/31 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.

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  is used, as in Anand and Ames (2006).

Kinematics
Consider a bar of length  0 in a given one-dimensional reference configuration at time  = 0, whose arbitrary material point is denoted by  and boundaries, denoted by , given by the two points at the bar ends, with normals with values -1 and 1.The motion of the bar, denoted by (, ), maps the undeformed to the current configuration with spatial points indicated by  and length , as illustrated in Figure 1.
The motion of each particle is described by a smooth one-to-one mapping  = (, ), in which (, ) =  + (, ), where  is the displacement field.The deformation gradient, which maps the material infinitesimal fibre  at point  to the spatial infinitesimal fibre  at point , velocity and velocity gradient in spatial description are given, respectively, by (, ) = (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   and plastic   parts of  as  =     . (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.
Latin American Journal of Solids and Structures, 2023, 20(6), e504 5/31 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: where  indicates a mechanism.Figure 2 illustrates the generalized Zener model with  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  should be considered, since finding the correct material parameters become an important aspect of these phenomenological continuum-level models.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   () and internal variables  () ) in parallel with an energy storage mechanism associated to plastic deformation (driven by  () ).Physically,   () represents an elastic deformation, such as stretching and rotation of intermolecular bonds and polymer chains.  () is a local inelastic deformation at , due to mechanisms such as slippage of polymer molecules.This last transformation takes material coordinates to a coherent structure in the relaxed space   , being   the range of   () , as shown in Figure 3.This relaxed space, also called sometimes intermediate or structural space, is defined as the output of the linear transformation   () and the input of the linear transformation   () .Thus,   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 with where is the inelastic velocity gradient in the intermediate configuration and   () 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 Analogously, polar decompositions of   () and   () provide with   () ,   () ,   () and   () being variables describing deformations.
The right and left Cauchy-Green deformations, for both elastic and plastic cases, are given, respectively, by ) , We define also which for 1D problems is equal to () in the case of compressibility and equal to one in case of incompressibility.

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 , a microscopic velocity () referring to the time rate of change of phase-field  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 .
Applying the PVP, consider that the external power comes from the macroscopic stresses  on the bar ends and the axial distributed loads () per unit of mass along the bar length.Moreover, it is considered a density of energy, , supplied by exterior to evolving structures (such as phase-field describing damage in material) and a surface density of energy given by flux  ℎ .Virtual macroscopic and microscopic velocities are denoted by  � and , respectively.Based on the previous definitions, the expression of the external virtual power is given by Considering the internal power, due to the existence of elastic and plastic stresses (  () and   () ) for each branch, the expression must be split in terms of virtual rates  �  and  �  .Analogously to the external power, microscopic terms   and ℎ are considered, representing, respectively, the volumetric density of energy exchanged by variation of a unit of Latin American Journal of Solids and Structures, 2023, 20(6), e504  7/31  and the flux of energy associated to the spatial variation of a unit of  in a unit of time.Based on previous definitions, the expression of the internal virtual power for a material with  mechanisms is given by The virtual power of the inertia (acceleration) loads is expressed as follows Given the previous expressions, the PVP is written as follows: Substituting Eqs.( 11), ( 12) and ( 13) in Eq.( 14), we have To obtain the dynamic equations, consider that , � is arbitrary and set  � and ̂ to zero.Therefore, Applying integration by parts in the first term on the right hand side and rearranging terms, the following equations are obtained: ) ) works as Cauchy stress for the problem.
Next, to find an expression for , set  � and ̂ to zero, and from Eq. ( 4) Equation ( 15) reduces to and therefore, ) . (20) Similarly to the cases above, let  �,  �  () and  �  () be zero.In this case, applying integration by parts and rearranging the terms, we obtain Latin American Journal of Solids and Structures, 2023, 20(6), e504  8/31 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 () and internal () energies.
Balance of energy states that the total energy rate is equal to the sum of the external power (  ) and the heat flux.Therefore, where the kinetic energy  does not consider microscopical velocities and is given by Moreover,  is the specific heat source density;  is the specific internal energy density, and  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  � =  and ̂= ̇, gives the reduced form of the balance of energy in the integral form, given by Due to the conservation of mass, we have ̇ , and, using also the divergence theorem, we obtain Next, the second law of thermodynamics is given in local form of the Clausius-Duhem inequality by with  the possible specific entropy density,   the entropy flux and   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 where  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, , which is given in Eulerian coordinates by Using the balance of energy and the Helmholtz free energy definition, dissipation inequality can be rewritten as Collecting the previous results, the basic governing equations in the spatial reference frame are ) . (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 Latin American Journal of Solids and Structures, 2023, 20(6), e504 9/31 with the Cauchy stress   () replaced by the first Piola-Kirchhoff stress  () .They are related as The microscopical balance of Eq. ( 21) is given now by with Rewriting the balance of energy of Eq. ( 24) in material reference leads to following equation: With For describing the dissipation inequality it is convenient to use the Clausius-Duhem and the Helmholtz free-energy definitions both in material reference.Therefore, and again the classical choices for the entropy flux and production terms are considered resulting in It is useful to rewrite the term   ()   ) of Eq.( 38) without variables in the spatial reference.For this purpose, partial right Cauchy-Green   () and   () and partial Green-Lagrange   () and   () strains are used and defined as The derivatives of the previous expressions are given by () . (40) Using these definitions and remembering that   () and are the symmetric parts of   () and   () , respectively, and using also the fact that   () is symmetric, results in the following identity: with where   () is the elastic second Piola-Kirchhoff stress.
Considering the Kroner-Lee decomposition of Eq.( 2), the first Piola-Kirchhoff stress is expressed by or, equivalently, Thus, the motion equation in Lagrangian coordinates is given by After simplifications, we obtain for the plastic term that Given that the skew-symmetric part of the gradient of velocity   is zero because there are no rotations, we have that with It is possible, then, to rewrite the dissipation inequality as Finally, from Eqs.( 42) and ( 48), the relation   () =   () obtained in the deformed configuration is then redefined in the Lagrangian frame of reference as follows Collecting the previous results, the basic governing equations in the material reference frame are ) . (51) Latin American Journal of Solids and Structures, 2023, 20(6), e504 11/31

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).
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  () and its evolution equation, which is the same as Eq.(52.v).This variable represents a dimensionless squared stretchlike quantity, whose evolution considers static (  () ) and dynamic ( () ) recovery terms.
Similarly to Frémond and Shitikova (2002),  �  ,   ,  0 , ℎ 0 and  0 are split in their reversible (non-dissipative) and irreversible (dissipative) parts, which are indicated, respectively, by the superscripts (. ) () and (. ) () as below ) , ) , ) , ) , Non-dissipative (reversible) parts may, in general, depend on , while the dissipative (irreversible) parts may depend both on  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 ̇ and  in Eq.( 49), and that the heat flux is purely dissipative (irreversible).Thus, we consider 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  () and using the evolution law of  () given in Eq.(52.v),we obtain  ()  ̇() with The entropy inequality (written in terms of the free-energy), after some manipulation and collecting similar terms is rewritten as 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, Since in Eq.( 59) we can give any values to  ̇,  ̇, ̇, ̇,  ė and  ̇ their coefficients can be set to zero.Thus, ()    () = 0, (61) By using the previous results, Eq.( 58) is reduced to The dissipation inequality above is, then, split into four inequalities, being the first related to the mechanical dissipation with representing the dissipative part of   () .Thus, it can be noticed that the plastic stress   () , conjugated to   () , is composed by a dissipative and a reversible (energetic) terms as Latin American Journal of Solids and Structures, 2023, 20(6), e504 13/31 () �  () ,  () � =   () �  () ,  () � + �2 −1  0  ()    ()  () �. (69) The energetic part is also called back-stress, and given by   () = 2 −1  0 �  () �  () , () ,� Proceeding with the analysis of the dissipation inequality, the last term in Eq.( 66) is the heat conduction inequality, whose positiveness is imposed by with ̃ a non-negative coefficient.
Analogously, the first term in Eq.( 66) corresponds to the microscopic energy inequality, whose positiveness is imposed by with  � a non-negative coefficient.For this work, the irreversible part of  �  to ensure positiveness is 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 , being the first one an elastic energy (isotropic with relation to   () ), the second a defect energy (isotropic with relation to  () ), associated to plastic flow, and the third related to damage.Therefore, For evaluation of stresses in plastic flow, considering that   () is the Mandel stress and defined as   () =   ()  �  ()() , the effective Mandel stress  () () is given by From the constitutive relation given in Eq.( 69) follows the definition: with  being a scalar flow strength of the material for each .
Given the relation for   () , using   () as the direction for plastic flow for each  (remembering that in 1D domain it is the sign of driving stress), the flow rule is given by = sign �  () � . (77)

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.
Latin American Journal of Solids and Structures, 2023, 20(6), e504 14/31 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   , 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.
Branches 2 and 3 represent responses under and above   , 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 crosslinks 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.

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  11 =  is considered in the model.Moreover, given the Kroner-Lee decomposition (Eq.( 3)) with simplification given in Eq.( 7), we have In the first branch, the elastic free energy is given by Latin American Journal of Solids and Structures, 2023, 20(6), e504 15/31 with   (1) being the elastic logarithmic strain measure,  is the Young modulus,  ℎ is the coefficient of thermal expansion and  0 is the reference temperature.
The dependence on temperature of  (following Dupaix and Boyce (2007)) is given by () = The defect free energy in the first branch is with   ( = 1,2,3) given under uniaxial conditions by  1 =  (1) , and , in accordance with Anand et al. (2009).Parameter () is the backstress modulus, assumed to be a linearly decreasing function of temperature with slope  (to be defined as material coefficient), when  <   .This backstress modulus is assumed zero when  ≥   .
The Mandel stress is defined as 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 Finally, the driving stress for plastic flow is given by the difference between Mandel stress and backstress as   = (1) −   .

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, with  () representing the ground state rubbery shear modulus of the material and   () is the maximum value of � 1 () − 3�, associated with the limited extensibility of the polymer chains.The term  1 () is given by where   () is given by Eq.( 78).
The stresses for branches 2 and 3 can be simplified to with =2 or 3.

Flow rules
The evolution equation in the first branch is given by Latin American Journal of Solids and Structures, 2023, 20(6), e504 16/31 with   (1) given by with  (1) the equivalent plastic shear strain rate;  (1) is the equivalent shear stress, given by In Eq.( 88), several material parameters must be calibrated, since their physical meaning are not easily extracted from experiments. 0 (1) is a pre-exponential factor,  is an activation energy,   is the Boltzmann constant,  is activation volume,  1 is a strain-rate sensivity parameter and   is the pressure sensivity of plastic flow.
The term   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 which has others material parameters (ℎ  ,   ,   ) to be defined and  given by Still in branch 1, the evolution of  is given by with  being a parameter driving the dynamic recovery of .
The flow rules for branches 2 and 3 are given by the followin simple power laws, respectively: with  0 () being the reference plastic shear strain rate and () the positive-valued strain-rate sensitivity parameter.

𝑆𝑆
(2) is a positive-valued shear resistance.In branch 3, which is linked to behavior in temperatures above   , the internal variable  (3) is related to a dissipative resistance during molecules movement, known as reptation movement.Its evolution is given by the following rule: with  given in Eq.( 91) and the initial conditions for  (3) and parameter ℎ 3 given by where  and  are parameters to be found in the calibration process.
The free energy related to elasticity corresponds to the first term on the right hand side of Eq. ( 97), expressed by a degradation term (()) and an energetic term ( 0 (1) ) as 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 where  > 0 is related to the damage rate for the virgin material that is, when  = 0.In Fig. 5, we show the graph of this function for several values of ; in this work, we take  = 0.001.
Further explanation about why we choose (99) as the degradation function can be found in Section 9.The free energy related to damage is where   is the Griffith's critical energy release rate,  a positive constant related to fracture layer width and  a (small) positive constant equal to 0.001 in this work.The expression for the potential ℋ is given by 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).

Approximation of damage equation
Using the constitutive relations for variables ℎ 0 ,  0 and  0 = 0 from Eq.(51.(iii)) and applying the free energy defined in Eq. ( 97), the damage equation is rewritten as Latin American Journal of Solids and Structures, 2023, 20(6), e504 18/31 The previous expression is multiplied by a test function  at both sides and integrated in the bar length  0 .After integrating by parts, the first term on the right-hand side, we have The functions with non-linear terms in , indicated generically as (), and the test functions  are approximated in a generic element  of the mesh, respectively, by where   () are the element shape functions,  is the polynomial order and  + 1 is the number of nodes of the element.
The approximation of the derivatives of  and  are given, respectively, by where   () () are the global derivatives of the element shape functions.
(106) Grouping terms and rearranging, the previous equation is rewritten as with 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: with  the time step and  the step number.
The starting value for  +1 (0) is given by converged values in last step   .The updated values in step  + 1 are given iteratively by

Approximation of motion equation
Based on the motion equation of Eq.(17.i) and given that the second Piola-Kirchhoff stress  is work conjugate to the virtual Green-Lagrange strain , the following expression for the virtual internal work is valid The expression for the virtual external work, for any virtual material axial displacement , is The virtual work   for the inertia load is given by The Principle of Virtual Work, given by   +   =   is, then written as The Newmark method employs the following time approximations for velocity and acceleration: with where   and   are the Newmark constants.
The expression for the global residue of the motion equation,   , is obtained as 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   [  ] in the direction of displacement increment, , considering that the forces are not dependent on the displacement field, is given by We have the following identities for the Green-Lagrange strain: is the constitutive function, also called tangent modulus, given by Latin American Journal of Solids and Structures, 2023, 20(6), e504 20/31 The residual expression, therefore, is rewritten as Analogously to the damage equations, the approximation of the material coordinate, in an element , using the local normalized coordinate system  is given by where   are the nodal coordinates of the initial mesh,   () () are the element shape functions,  is the polynomial order and  + 1 is the number of nodes of the element.
For an isoparametric element, the material axial displacement and strain are interpolated as where   () () 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   () , as follows Finally, the displacement increment  and its derivative are interpolated as where {  () } is the vector of nodal displacement increments in the material configuration.

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.
Latin American Journal of Solids and Structures, 2023, 20(6), e504 21/31 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 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.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 with  being the stresses evaluated at each time step of the simulation and from experiment.The starting parameters were taken from Marin et al. (2007) andSrivastava 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.

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 × 10 −3  −1 .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.

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.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).

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  ̃ 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 −3 to be compliant with experiments described by Rae and Brown (2005).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.

CONSIDERATIONS ABOUT THE DEGRADATION FUNCTION 𝑫𝑫(𝝋𝝋)
The degradation function () 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.
Latin American Journal of Solids and Structures, 2023, 20(6), e504 26/31 Since degradation involves changes in the microstructure, () 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 (), 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, () is chosen to fit to the observed material response.The importance of the choice of () 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 () = (1 − ) 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 () = (1 − ) 2 instead of (99).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 ℐ(, ) has two parts: the bulk damage energy Here  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  has a rather large influence on the velocity of the crack propagation.The value of the material parameter  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.

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.

Figure 3 :
Figure 3: Reference spaces including structural space given by the Kroner-Lee decomposition.
�� − � −   �, with   and   being the Young modulus in glassy and rubbery regions, respectively.Parameter  is related to size of transition close to   , and  represents slope of the shear modulus in terms of temperature, being divided in two regions (pre and post   ).

Figure 5 :
Figure 5: Degradation function for different values of parameter a.

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

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

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

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

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

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

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

Figure 14 :
Figure 14: Stress-strain curves of the 1D finite element model at -15 o C using D(ϕ) = (1− ) 2 .For each plot above, only one of the parameters of the damage equation (, , ) and the damage energy at the faces of the smeared cracks:

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