Georges Heyihin
email: gheyihin@gmail.com

Christophe Awanto
email: christophe.awanto@gmail.com

Franois Lanzetta
email: francois.lanzetta@univ-fcomte.fr

Thermodynamic analysis of the Stirling Duplex machine 1 st

Keywords: Stirling Duplex, adiabatic model, engine, cooler

The Stirling Duplex concept, in which a Stirling engine drives a Stirling heat pump, has many energy and environmental benefits. This machine is essentially composed of three movable elements, a working piston and two displacers, in the same enclosure. The machine has two circuits of nonpolluting working fluid, in both parts, engine and refrigerator. The compatibility of the Stirling machine with any type of thermal energy as an external combustion engine contributes to its industrial interests and scientific research. In this paper, a study is presented to estimate Stirling Duplex s parameters and efficiencies by considering the adiabatic model.

I. INTRODUCTION

Ever-growing energy needs and the resulting environmental problems have intensified research into renewable energies and new cleaner engines. The operation of Stirling machines is in line with the efficient use of energy and the protection of the environment, and as a result they have become the main focus of recent scientific and industrial research. Numerous scientific institutes or commercial companies have invested to contribute to the development of Stirling machines [START_REF] Walker | Future coal-burning Stirling engines[END_REF]- [START_REF] Thombare | Technological development in the Stirling cycle engines[END_REF]. The Stirling engine was invented by Robert Stirling in September 1876, and was used to train wagons in coal mines. Due to the lack of materials and the high cost above those of others machines, the Stirling engine has not had much industrial success. But the scarcity of fossil fuels and the use of renewable energies have led to renewed research around Stirling machines. The Stirling engine is a high efficiency, simple, quiet and highly reliable tool that is non-polluting and requires little maintenance [START_REF] Beale | Understanding Stirling engines[END_REF]. Its compatibility with any type of thermal energy as an external combustion engine contributes to its scientific and industrial interests. The reversibility of the Stirling cycle is used for cold generation. Driven by a motor, the machine becomes a refrigerating machine or a heat pump. The Duplex machine used as a heat pump or refrigerator is essentially two freepiston Stirling machines sharing a common piston and an under pressure enclosure. When the high temperature heat is supplied to the heater, the "motor" part produces mechanical energy that is used by the heat pump to extract heat from the cold source to the hot well [START_REF] Penswick | Duplex-Stirling-Machines[END_REF]. The domestic refrigeration and air-conditioning market is currently dominated by refrigerant compression machines containing chlorofluorocarbons (CFCs) which destroy the ozone layer. In fact, CFCs released into the atmosphere have an important power of chemical transformations of the ozone layer, protecting the earth against ultraviolet rays. The Stirling refrigerator would therefore be an environmental solution against the use of chlorinated refrigerant compression machines. The use of Stirling Duplex will solve two problems of the century simultaneously: -Energy solution: the use of renewable energies such as solar energy, biomass, etc., as thermal energy at the input of the engine, to the detriment of fossil energies; -Protection of the environment: the reduction of refrigerants destroying the ozone layer. In this document, the performances of the Stirling Duplex for refrigeration production are estimated.

II. DESCRIPTION OF THE STIRLING DUPLEX MACHINE A. Presentation of the Stirling Duplex machine

The Stirling Duplex machine is the type of combined cycle machine where Stirling cycle engine drives the reverse Stirling cycle. This device used as a refrigeration tool is essentially composed of two free-piston Stirling machines arranged (backto-back), linked together by a common piston in a same enclosure. A set of displacer and piston operating in motor, drives another set, which in turn works as a heat extractor. In both parts, the displacer is driven by the pressure differential between the compression and expansion spaces on either side of the displacer. Each section of the Stirling Duplex (engine or cooler) contains all the main elements of a Stirling machine such as the heater, the regenerator, the cooler, the displacer and the piston [START_REF] Penswick | Duplex-Stirling-Machines[END_REF]. The uniqueness of the design comes from the fact that the two Stirling machines have the same piston, resulting from the configuration of three moving components (Figure 1). When the high temperature heat is supplied to the heater, the "motor" portion generates mechanical energy that is used by the heat pump to extract heat from the cold source to the hot well. The system from these three components is a potentially efficient, simple and durable heat pump pump, with, of course, all the advantages associated with free-piston Stirling engines. The concept of the free-piston Stirling Duplex engine has its roots in the history of free-piston engines. Current Stirling machines such as refrigerators or heat pumps are built at SUNPOWER, and have been used all over the world since 1979 [START_REF] Penswick | Duplex-Stirling-Machines[END_REF].

B. Ideal Stirling Duplex cycle

The ideal Stirling cycle is a reversible thermodynamic cycle with two isothermal transformations and two isochoric transformations. As the two parts of the Stirling Duplex work with separate working fluids, the cycle described by the Stirling Duplex machine is dual: the motor cycle and the cooler cycle (Figure 2) [START_REF] Chen | Thermally Actuated Duplex Free-Piston Stirling System for domestic, cooling and heating application[END_REF]. in this phase, the displacer continues to move while the piston remains stationary. Like any relaxation, the pressure of the fluid drops in front of the increase of the volume of the expansion chamber. The temperature remains constant by the external heat input which compensates the cooling. At the engine, work is provided and at the heat sink is absorbed from the medium to be cooled. • Isochoric cooling: (Engine: 4-1, Refrigerator 4'-1 ').

During this phase, the two pistons move simultaneously transferring the working fluid from the expansion space to that of compression through the regenerator, without volume variation. During this transfer of fluid, heat is transferred to the regenerator thereby reducing the temperature of the gas. At the chiller, the amount of heat is transferred from the regenerator to the fluid, thereby increasing its fluid temperature.

III. THERMODYNAMIC EQUATIONS

The adiabatic cycle assumes that the compression and expansion spaces are perfectly insulated. This is more realistic and adiabatic conditions can be safely assumed since the machines operate at high frequencies allowing little or no time for heat transfer in the working spaces. All heat input to the cycle occurs in the heater, and all heat output occurs in the cooler. Gases leave the heater at the heat source temperature and are mixed perfectly as soon as they enter the expansion space. Similarly, gases leave the cooler at the heat sink temperature and are mixed perfectly as soon as they enter the compression space. In this case too, perfect regeneration is assumed. This section describes the adiabatic model proposed by Urielli [START_REF] Martini | Stirling Engine Design Manual[END_REF] , where each of the two machines parts is represented in five compartments :

• Compression space (c) • Cooler exchanger (k) • Regenerator (r) • Heater exchanger (h) • Expansion space (e)
This study is based on the following assumptions:

• The expansion and compression processes are adiabatic • The temperature of hot end heat exchanger, cold end heat exchanger and regenerator is constant. • The working fluid obeys ideal gas law.

• The specific heats for the working fluid are constant.

• The pressure is uniform throughout each of two machines parts at all instants. • Perfect regeneration occurs in the regenerator volume. 1) Pressure and volume variation: The evaluation of the pressure along the Stirling machine is calculated considering the mass balance applied to the overall machine.

m = m c + m k + m r + m h + m e (1)
The mass balance is re-written in a form that allows the calculation of the pressure.

p = m • r Vc Tc + V k T k + Vr Tr + V h T h + Ve Te (2)
The Duplex machine is similar to a combination of two (2) beta-type Stirling machines associated "back-to-back"; the geometrical parameters of each section of the Duplex are those of the Stirling beta. The volumes of the compression and expansion chambers as a function of the harmonic movement of the piston are given by the following formulas [START_REF] Awrejcewicz | Wydawnictwo, and Conference on Dynamical Systems -Theory and Applications[END_REF]:

V c = V swe 2 • (1 -cos(ϕ)) + V swc 2 • (1 -cos (ϕ -α)) (3)
V e = V swe 2 • (1 -cos(ϕ)) (4)
2) Mass drivatives and mass flow: The derivatives of the mass i respect to the angle position of the power piston for each heat exchanger, which are assumed isothermal, are calculated from equation ( 5)

dm i dϕ = m i p • ∂p ∂ϕ , i = k, r, h (5)
The derivatives for the compression and expansion spaces, which are not isothermal, are evaluated by the equations below:

dm c dϕ = p • ∂Vc ∂ϕ + Vc• ∂p ∂ϕ γ R • T ck (6
)

dm e dϕ = p • ∂Ve ∂ϕ + Ve• ∂p ∂ϕ γ R • T he (7)
The mass flows at the interfaces are evaluated with equations ( 8), ( 9) , [START_REF] Awrejcewicz | Wydawnictwo, and Conference on Dynamical Systems -Theory and Applications[END_REF], [START_REF] Araoz Ramos | Thermodynamic analysis of Stirling engine systems applications for combined heat and power[END_REF] applied to each control volume:

• Between the compression space and the cooler:

m ck = - ∂m c ∂ϕ (8)
• Between the cooler and the regenerator:

m kr = m ck - ∂m c ∂ϕ (9)
• Between the regenerator and the heater:

m rh = m he + ∂m h ∂ϕ (10)
• Between the expansion space and the heater:

m he = ∂m e ∂ϕ (11)
3) Pressions variation: The analytical formula for change in pressure based on energy balance is given by:

dp dϕ = - γ • p • ∂Vc ∂ϕ T ck + ∂Ve ∂ϕ T he Vc T ck + γ • ( V k T k + Vr Tr + V h T h ) + Ve T he (12)
4) Temperature variation: The variation of the temperature inside the Stirling machine is assumed to follow the scheme shown below [START_REF] Araoz Ramos | Thermodynamic analysis of Stirling engine systems applications for combined heat and power[END_REF]. This scheme shows that the heater and cooler temperatures are assumed to be constant and the regenerator temperature is assumed to vary linearly. The linear variation permits the evaluation of a mean effective temperature for the regenerator space.

T r = T h -T k ln( T h T k ) (13)
The temperatures in the compression and expansion spaces change over the cycle. These are calculated with the differential form of the equation of state for the working gas.

dT c dϕ = T c • ( ∂p ∂ϕ ) p + ( ∂Vc ∂ϕ ) V c - ( ∂mc ∂ϕ ) m c ( 14
)
dT e dϕ = T e • ( ∂p ∂ϕ ) p + ( ∂Ve ∂ϕ ) V e - ( ∂me ∂ϕ ) m e (15)
The boundary temperatures Tck and The change with the direction of the flow. Therefore the following conditional equations are used to calculate these temperatures.

If m ck > 0 then T ck = T c else T ck = T k (16) If m he > 0 then T he = T h else T he = T e (17)
The other pair of boundary temperatures (Tkr and Trh), are assumed constant and equal to the cooler and heater temperatures respectively for the ideal adiabatic model these

T kr = T k (18) T rh = T h (19)
5) Energy balances : The energy balance applied to each heat exchanger, taking into account that these are isothermal spaces, gives the following expressions:

• Cooler:

dQ k dϕ = V k • ( ∂p ∂ϕ ) • C v R -C p • (T ck • m ck -T kr • m kr ) (20)
• Regenerator:

dQ r dϕ = V r • ( ∂p ∂ϕ ) • C v R -C p • (T kr • m kr -T rh • m rh ) (21)
• Heater:

dQ h dϕ = V h • ( ∂p ∂ϕ ) • C v R -C p • (T rh • m rh -T he • m he ) (22)
The compression and expansion work are calculated from the energy balance applied to the corresponding working spaces, which are assumed adiabatic.

dW c dϕ = p • ( ∂V c ∂ϕ ) ( 23
)
dW e dϕ = p • ( ∂V e ∂ϕ ) ( 24
)
The work output of the Stirling engine or refrigerator is the sum of the compression and expansion work on the machine.

W = W c + W e ( 25
)
The mechanical powers absorbed in the engine section, and transmitted in the chiller section by the power piston are given by and [START_REF] Penswick | Duplex-Stirling-Machines[END_REF]: These powers are a function of the phase () of the movement of the displacer relative to to the movement of the power piston

P m = p m • A • X • π • f • sin(ϕ m ) (26)
P f = p f • A • X • π • f • sin(ϕ f ) (27)
Like all combined cycle machines, the work indicated at the output of the engine section is equal to that of the heat pump [START_REF] Kagawa | An experimental study of a 3-kW Stirling engine[END_REF].This makes it possible to write the following relationship between the average pressures in the Stirling engine and the Stirling cooler:

p f = p m • sin(ϕ m ) sin(ϕ f ) (28)
6) Coefficient of performance: Two performance factors can be considered for the Stirling Duplex machine: it is COP h when it is used as a heat pump, and the COP f as a refrigerating machine.

COP h = Q km + Q kf Q h (29) COP f = Q f Q h (30)

IV. SIMULATION RESULTS AND DISCUSSIONS

In what follows, the numerical solution algorithm is presented, see Fig. 4. The adiabatic Stirling Duplex model is analysed and simulated to evaluate its performance and operating parameters. The basic parameters and simulation conditions are presented in Table 1. The geometric dimensions are those of a Stirling engine used by A. Asnaghy et al in their work in 2012 [START_REF] Asnaghi | Thermodynamics Performance Analysis of Solar Stirling Engines[END_REF] These works, represented numerically by the curved areas in the P-V diagram, are called the indicated works of the machine. Table 2 shows the characteristics and performance of the Stirling machine for one operating cycle. 1) Volume variation: Fig. 6 and Fig. 7 show the volume variations respectively the engine-Stirling and cooler-Stirling sections. The variation of the total volume of each Stirling section depends almost on the variation of the volume of the compression space. This is explained by the shape of the Duplex machine composed of two Stirling beta where the instantaneous volume of the compression space is the sum of the volumes "scanned" by the displacer on the one hand and by the power piston on the other hand share.

At the Fig. 8 illustrating the volumes in the compression and expansion spaces, the variations of the expnsion volumes while those of the compression volumes are out of phase by π 4 rad; this is the phase's difference of the displacers with respect to the piston [START_REF] Penswick | Duplex-Stirling-Machines[END_REF].

2) Pressure variation: Fig. 9 shows the variation of the average pressures in the two sections of the machine. The pressures variation are similar in that they are linked by the periodic phase evolution of the gas in the two sections (28).

3) Thermal balance: The Stirling Duplex can be used as a refrigerating machine or a heat pump. Fig. 10 shows the variations of heat quantities, supplied at the input of the machine, produced as a heat pump and extracted as a refrigerating machine. It appears that the Stirling Duplex is more efficient as a heat pump than as a refrigerating machine (Table II). In this study, several parameters of the Stirling Duplex machine are discussed. The working fluid considered is helium for its thermal characteristics (high heat transmission capacity) and chemical characteristics (stability) [START_REF] Asnaghi | Thermodynamics Performance Analysis of Solar Stirling Engines[END_REF]. The adiabatic model used, based on the work of W. Martini [START_REF] Martini | Stirling Engine Design Manual[END_REF] and L. Penswick [START_REF] Penswick | Duplex-Stirling-Machines[END_REF], made it possible to highlight these studied parameters, such as the average pressures in the two sections of the machine, the different volumes depending on the position of the power piston, the temperatures the amounts of heat and work supplied and received, the coefficients of performance of the machine used as a refrigerating machine and as a heat pump.

Fig. 1 .

 1 Fig. 1. Representation of the Stirling Duplex machine.

Fig. 2 .

 2 Fig. 2. : Operation Cycle of the Duplex machine

Fig. 3 .

 3 Fig. 3. : Temperature variation inside the control volumes

Fig. 4 .

 4 Fig. 4. Flowchart of the simulation

Fig. 5 .

 5 Fig. 5. Pression-Volume diagram of the Duplex

Fig. 6 .Fig. 7 .Fig. 8 .Fig. 9 .Fig. 10 .

 678910 Fig. 6. Variation volume in the engine-section


T Barriere

X Gabrion

S Holopainen
email: sami.holopainen@tuni.fi

A compact constitutive model to describe the viscoelastic-plastic behaviour of glassy polymers: Comparison with monotonic and cyclic experiments and state-of-the-art models

Keywords: Cyclic plasticity, Thermodynamics, Implementation, Experimentation, Ratcheting

come

A compact constitutive model to describe the viscoelastic-plastic behaviour of glassy polymers : Comparison with monotonic and cyclic experiments and state-of-the-art models

Thierry Barriere, Xavier Gabrion, Sami

Introduction

Examples of engineering components that are manufactured from polymers and are subjected to cyclic loadings are from fields such as automotive and aerospace systems, marine structures, healthcare instruments, and sporting

τ τ L e a) b) c) τ β dev τ L e a) (1)
(2)

c) µ B τ B µ 1 τ A b) τ L e a) (1)
(2)

(N+1) c) µ B τ B µ 0 µ 1 µ N τ A b)
Figure 1: Representation of the three dimensional version of the [START_REF] Haward | The use of a mathematical model to describe isothermal stress-strain curves in glassy thermoplastics[END_REF]-model (according to [START_REF] Boyce | On the kinematics of finite strain plasticity[END_REF]), the proposed model, and the [START_REF] Anand | On modeling the micro-indentation response of an amorphous polymer[END_REF]-model. The models are governed by the following elements: a) an elastic spring, b) a nonlinear system (from a single dashpot to Kelvin-Voigt-like elements, i = 0, 1, ..., N ), and c) a nonlinear spring.

This paper studies the modeling of the mechanically motivated cyclic longterm deformation behavior without a remarkable temperature rise and which behavior occurs at relatively small stress frequencies; practically, the stress levels under such conditions remain below the peak yield stress. The paper aims at developing a model capable of predicting more accurately, the shape of the hysteresis loops and amount of energy-dissipation under cyclic loading processes. Another important objective is an exact prediction of the ratcheting behavior including cyclic hardening. Furthermore, the model is suitable for the simulation of very large deformations, i.e. it is based on the decomposition of the deformation gradient into viscoelastic-plastic and purely elastic parts, cf. [START_REF] Boyce | On the kinematics of finite strain plasticity[END_REF]; [START_REF] Anand | A theory of amorphous solids undergoing large deformations with application to polymer glasses[END_REF]. The formulation of the model follows both the micro-mechanics (anisotropic hardening, [START_REF] Boyce | On the kinematics of finite strain plasticity[END_REF], and yielding including free volume, [START_REF] Anand | On modeling the micro-indentation response of an amorphous polymer[END_REF]), and internal state variable methodologies (viscoelasticity, [START_REF] Anand | On modeling the micro-indentation response of an amorphous polymer[END_REF]), and stems from the rheological (spring-dashpot) representations displayed in Fig. 1.

In detail, the proposed approach relies on the pioneering model of [START_REF] Haward | The use of a mathematical model to describe isothermal stress-strain curves in glassy thermoplastics[END_REF] for glassy polymers (based on Eyring's theory) that is augmented by a single set of internal variables able to predict the missing viscous deformation behavior (area/shape of the hysteresis loops) and the nonlinear cyclic hardening (ratcheting behavior) in three dimensions. Against the state-of-the-art models of [START_REF] Anand | A theory of amorphous solids undergoing large deformations with application to polymer glasses[END_REF]; [START_REF] Anand | On modeling the micro-indentation response of an amorphous polymer[END_REF]; [START_REF] Jiang | A viscoelasticplastic constitutive model for uniaxial ratcheting behaviors of polycarbonate[END_REF] that can also predict the cyclic deformation behavior well, the elastic portion of the deformation is described by a single element a) clearly separated from viscoelasticplastic elements in accordance with the classical [START_REF] Haward | The use of a mathematical model to describe isothermal stress-strain curves in glassy thermoplastics[END_REF] -model. As a result, the nonlinear Langevin spring c) applied in the model to describe the anisotropic hardening of amorphous network structures is modeled solely using the viscous deformation (as originally proposed in [START_REF] Boyce | On the kinematics of finite strain plasticity[END_REF] for three dimension). Moreover, the specified stress (τ A in Fig. 1) defines the stress in each micromechanism. Compared with state-of-the-art models, [START_REF] Anand | On modeling the micro-indentation response of an amorphous polymer[END_REF]; [START_REF] Jiang | A viscoelasticplastic constitutive model for uniaxial ratcheting behaviors of polycarbonate[END_REF], the proposed model requires a reduced set of internal state variables and material parameters to be defined. The model predictions are evaluated by comparison with the experiments for uniaxial tension under different stress ratios, as well as with the consistent stateof-the-art models' predictions, [START_REF] Boyce | On the kinematics of finite strain plasticity[END_REF]; [START_REF] Anand | On modeling the micro-indentation response of an amorphous polymer[END_REF]; [START_REF] Jiang | A viscoelasticplastic constitutive model for uniaxial ratcheting behaviors of polycarbonate[END_REF].

This paper is arranged as follows: Kinematics and a thermodynamically consistent constitutive model are introduced subsequently. Section 3 starts with the description of the applied tests, and the calibration procedure for the material parameters is given in Section 3.2. Then, the effectiveness of the approach under creep/recovery and different cyclic loadings for the technologically imperative polymer is discussed in Section 4. Finally, the important conclusions and hotspots for future research are presented.

The approach

Deformation kinematics

At time t ∈ R + , deformation of a body P is determined by the map y : X → x in which X ∈ B and x ∈ b denote the material points determined in the undeformed reference and deformed current placements B and b, respectively1 . The motion of b is then given by the map x = y(X, t) and a deformation gradient F := ∇y, in which J := det(F ) > 0 is the Jacobian and the operator ∇ denotes the gradient with respect to the point X.

The proposed model is founded on the Kröner-Lee decomposition, i.e. the decomposition

F = F e F vep , (1)
where F e and F vep determine the local deformations for the intermediate placement N resulting from elastic and viscoelastic-viscoplastic mechanisms. Whereas the elastic part represents the reversible elastic mechanisms of the chain network, the viscoelastic-plastic part denotes the partially reversible mechanisms of molecular chains, [START_REF] Anand | A theory of amorphous solids undergoing large deformations with application to polymer glasses[END_REF]. The viscoelastic-plastic part manifests as long-term creep strain and stress relaxation, due to the inertia of the microstructure to reach its equilibrium. Purely viscoelastic effect appears as recovery of strain after a stress removal, whereas the plastic effect is due to irreversible, dissipative mechanisms, such as chain breakage and the slippage of chain entanglements, [START_REF] Bouvard | A general inelastic internal state variable model for amorphous glassy polymers[END_REF].

The intermediate placement N , discussed above, is merely locally determined: F e and F vep respectively, being not real deformation gradients. However, det(F e ) > 0 and det(F vep ) > 0 hold, i.e. their inverses exist. Because F e and F vep particularize the elastic and viscoelastic-plastic deformations, respectively, they are regarded as uncoupled, and N is stress-free.

Applying the polar decomposition,

F e = v e R e , (2)
where v e is the symmetric stretch tensor and R e is the rotation associated with the elastic deformation2 . The viscoelastic-viscoplastic deformation is defined by a symmetric stretch tensor (in the intermediate placement)

Cvep := F vep F vep,T .
(3)

In ( 3), the transpose of tensor is symbolized by the letter T. To represent purely elastic deformation, symmetric stretch tensor

b e := F e F e,T =: (v e ) 2 (4)
given in the spatial placement is also defined. Moreover, the symmetric stretchlike tensors Cvp , Cve

are introduced. These tensors are given in the intermediate placement N characterizing the internal state and its change due to the viscoplastic and viscoelastic deformations, elements (1) and (2) in Fig. 1. They can be viewed as the linear transformations that map vectors in the relaxed space.

Rate kinematics

To model the motion of a solid body, the tensorial velocity gradient

l := grad(v) = Ḟ F -1 , (6)
in which v is a velocity field, is defined; grad is the gradient with respect to a deformed point x defined above, and the marker dot means the time derivative. The decomposition (1) in ( 6) yields

l = l e + l vep (7)
where

l e := Ḟ e F e-1 , l vep := φ * ( Lvep ), Lvep := Ḟ vep F vep-1 . ( 8
)
The operator φ * in (8) denotes the push-forward operation of the quantity in the parentheses by using F e , cf. (Belytschko et al., 2000, Box 5.16). The specific form of φ * will be defined in Section 2.3. Using (3) and (8)(right), the rate of Cvep becomes

Ċvep = Ḟ vep F vep,T = Ḟ vep F vep,T + F vep Ḟ vep,T = Lvep Cvep + Cvep Lvep,T .
(9) Based on ( 9), ( 4) in ( 8) results in ḃe = l e b e + b e l e,T .

(10)

The complements (symmetric, sym, and skew-symmetric, skew) are given by: 12) respectively. It then follows from (7) as follows:

d := sym(l),
l := d + ω = d e + d vep + ω e + ω vep (13)
where d := d e + d vep and ω := ω e + ω vep . As a departure from the previous formulations for polymers, cf. e.g. (Anand and Ames, 2006, eq. (19) with the presumption that the spins vanish therein), it is further defined that l vep := l vp + l ve , l vp := d vp + ω vp , l ve := d ve + ω ve .

(14)

Thermodynamics of the model

Thermodynamic treatment as introduced in [START_REF] Anand | A theory of amorphous solids undergoing large deformations with application to polymer glasses[END_REF]; [START_REF] Holopainen | Modeling of mechanical behavior of amorphous solids undergoing fatigue loadings, with application to polymers[END_REF] is used here to describe the dissipative, irreversible and prolonged effects governing the plastic and viscoelastic deformations, respectively. The basis of the theory is the virtual power comprising of two prerequisites:

[1 ] Power balance: the external and internal virtual powers are equal, i.e. W ext (A, V) = W int (A, V) within a part of the body A(t) ⊆ b(t) for all virtual velocities V;

[2 ] Invariancy under all changes in frame: the internal power W int (A, V) must be frame-indifferent.

In order to treat the force and moment balances (symmetry of stress), which are the consequences of an expenditure of the external virtual power, we refer to previous works, see e.g. [START_REF] Anand | On modeling the micro-indentation response of an amorphous polymer[END_REF]; [START_REF] Holopainen | Modeling of mechanical behavior of amorphous solids undergoing fatigue loadings, with application to polymers[END_REF].

Frame-indifference

A theory is required to be invariant under changes in a frame of the space (Euclidean), that considers time-dependent transformations as follows:

y(X, t) → q(t) + Q(t)(y(X, t) -O) (15)
in which O is a fixed origin, q(t) a vector, and Q(t) is an orthogonal rotation tensor defined at time t for which

Q -1 = Q T holds.
Considering the reference placement to be invariant, whereas the current and intermediate placements depend on the choice of changes in frame,

F → QF , F vep → QF vep F e → QF e QT ( 16
)
where Q and Q represent the rotations of the current and intermediate placements. Then, the deformation measures in (3) and ( 4) are transformed by:

Cvep → Q Cvep QT , b e → Qb e Q T . (17)
The transformations ( 16) in ( 6) and ( 8) result in

l → QlQ T + QQ T and d → QdQ T , Lvep → Q Lvp QT + Q QT , Dvep → Q Dvp QT , l vep → Ql vep Q T -QF e-T QT QF e-1 Q T , l e → Ql e Q T + QQ T + QF e-T QT QF e-1 Q T . (18)
When leading the results, the property QT Q = -QT Q has been used. Then, based on (18) l = l e + l vep → Q(l e + l vep )Q T + QQ T = QlQ T + QQ T . Supposing that virtual velocity fields transform as their non-virtual counterparts do, the prerequisite [2] argues that in the current frame, W int (A, V) = W int * (A * , V * ), where a part of the solid body A converts to a region A * and V to a virtual velocity V * . Similarly, the stress field τ transforms to τ * conjugate to virtual velocities when

W int (A, V) := A τ : ldv = W int * (A * , V * ) := A (τ * : (Q lQ T + QQ T )dv (19)
where A : B := trace(AB). Because the body part A is arbitrary, (19) can be localized, i.e. τ : l = τ * : (Q lQ T + QQ T ). Since the change in frame, l, as well as QQ T = -Q QT (skew-symmetric) are arbitrary, the stress (symmetric) transforms by τ → Qτ Q T .

Benefiting (18) in ( 19) and considering the symmetry of τ yields

Constitutive theory

The proposed model includes a linear spring a) and a nonlinear spring c) resulting in an initial elastic response and anisotropic effect of material in large deformations, see Fig. 1. To predict increasing macro-yielding followed by a notable softening effect, the model also includes a nonlinear dashpot (1). To improve the predictions of the cyclic hysteresis loops, the model is augmented by a system for viscoelastic micro-mechanisms consisting of a KelvinVoigt-like element (2) arranged in series with the dashpot (1). Moreover, the element (2) governs creep and recovery, and its action together with the spring a) is used to model stress relaxation.

Using the equations ( 7) and ( 11), and noting the symmetry of stress τ , the localized dissipation in its spatial form becomes

D = τ : d e + τ : d vep -ψ ≥ 0, ( 21
)
where D is the power of the local dissipation and ψ stands for the rate of the Helmholtz free energy per unit volume. The Helmholtz free energy ψ is taken to be a sum of several independent potentials, [START_REF] Chaboche | Thermodynamic formulation of constitutive equations and application to the viscoplasticity and viscoelasticity of metals and polymers[END_REF]; [START_REF] Anand | On modeling the micro-indentation response of an amorphous polymer[END_REF], i.e.

wherein ψ e , ψ vep , and ψ ve are the contributions associated with the elastic spring a), the nonlinear Langevin spring c), and the nonlinear spring in the viscoelastic element (2), respectively. The viscoelastic potential ψ ve is also taken to be the function of the free volume, ϕ4 .

Remark 2.1. To guarantee the invariance of the free energy, finally invariants of the tensor arguments in the potential function (22) will be used.

The rate φ is given by ψ = ∂ψ e ∂b e : ḃe +

∂ψ vep ∂ Cvep : Ċvep + ∂ψ ve ∂ Cve : Ċve + ∂ψ ve ∂ϕ φ. ( 23
)
Taking advantage of ( 10 (24)

In accordance with (24) and noting the symmetry of Cvep by ( 9),

∂ψ vep ∂ Cvep : Ċvep = 2 ∂ψ vep ∂ Cvep Cvep : Dvep =: Bvep : Dvep = β vep : d vep (25)
in which Bvep represents an internal back stress quantity that controls the anisotropic hardening at large strains owing to the viscous work accumulated in a material body. The spatial complement β vep is derived by the push-forward operation determined in ( 20) and represents the stress in branch B (Fig. 1), i.e. τ B := β vep . Moreover, using the symmetry of Cve given by ( 5), one can define:

∂ψ ve ∂ Cve : Ċve =: 2 ∂ψ ve ∂ Cve Cve : Dve =: Bve : Dve = β ve : d ve ( 26
)
where Bve is an internal stress acting on the spring in the element (2). Its spatial counterpart in the deformed current placement is β ve . Consistent with (26), Ċve := Dve Cve + Cve Dve .

A complementary dissipation potential, associated with the dashpot (1) and that in the element (2) (cf. Fig. 1), is as

φ(τ (1) , τ (2) ) := φ vp (τ (1) ) + φ ve (τ (2) ) (28)
where τ (1) := τ A and τ (2) := τ Aβ ve stand for the driving stresses in the dashpots of the elements ( 1) and ( 2), respectively. Benefiting the potential φ, a part of the dissipation power becomes

D = ∂φ vp ∂τ (1) : τ (1) + ∂φ ve ∂τ (2) : τ (2) . (29)
Considering ( 29), substitution of ( 24), (25), and ( 26) into ( 21) results in τ -2 ∂ψ e ∂b e b e : d e = 0,

(τ -β vep ) : d vep -β ve : d ve - ∂ψ ve ∂ϕ φ - ∂φ vp ∂τ (1) : τ (1) - ∂φ ve ∂τ (2) : τ (2) = 0. (30)
Equation ( 30) must be valid under any possible thermodynamically admissible mechanisms, when the constitutive equation

τ = 2 ∂ψ e ∂b e b e (31)
is obtained. Defining ṡ(1) := -∂φ vp ∂τ (1) : τ (1) and ṡ(2) := -

∂φ ve ∂τ (2) : τ (2) , (32)
which represent the rates of the decreasing internal variables in the dashpots of the elements (1) and ( 2), (30) 2 can be written in the compact form as:

τ A : d vep = β ve : d ve + ∂ψ ve ∂ϕ φ -ṡ(1) -ṡ(2) . ( 33
)
Theorem 2.2. The model satisfies the dissipation (Clausius-Duhem) inequality D ≥ 0.

Proof. Considering ( 31), (32), and ( 33) in ( 30), the local dissipation power in (21) becomes:

D = -ṡ(1) -ṡ(2) . ( 34
)
Since ṡ(1) ≤ 0 and ṡ(2) ≤ 0 in (34), D always is greater than or equal to zero

Distinct constitutive equations 2.4.1. Elasticity

The approach is refined to address amorphous polymers. The chain network structure of the amorphous polymers is practically disorderly when a usual isotropic strain energy function,

ψ e = 1 2 κ e (I e 1 ) 2 + 2µ e J e 2 , (35)
where κ e and µ e denote the bulk and the shear moduli, is suitable. The invariants in ( 35) are given by:

I e 1 := trace(ln v e ) = ln J e , J e 2 := 1 2 (ln v e ) dev : (ln v e ) dev ,
wherein J e := det(v e ). The stress equation follows from ( 31) and ( 35):

τ = L e : ln v e , (36)
where L e is the standard fourth order elasticity tensor, see (Holopainen and Wallin, 2012, eqs. (11,12)), (Holopainen et al., 2017, eq. ( 47)).

Inelasticity Flow rule

According to ( 14), the rate of viscoelastic-plastic deformation is given by the sum

d vep = d vp + d ve . ( 37
)
The viscoelastic contribution d ve , associated with the mechanism (2) in Fig. 1, is determined by a rule as:

d ve := γve n (2) , n (2) := τ (2),dev √ 2τ (2) , τ (2) := 1 2 τ (2),dev : τ (2),dev (38)
where τ (2),dev = τ A,devβ ve,dev and γve represents a positive, effective viscoelastic strain rate defined subsequently. The deviatoric part dev for each tensor [•] is defined via the identity i as:

[•] dev := [•] -1/3trace([•])i.
Consistent with (38), the rate d vp , associated with the mechanism (1), is determined as

d vp := γvp n (1) , n (1) := τ (1),dev √ 2τ (1) , τ (1) := 1 2 τ (1),dev : τ (1),dev (39)
where τ (1),dev = τ A,dev = τ dev -β vep,dev and γvp represents a positive, effective viscoplastic strain rate. It then follows from ( 38) in ( 27 (40) or det( Cve ) = 1 (equals to the undeformed value). Similarly, using (9) and noting (37), det( Cvep ) = 1, i.e. the viscoelastic-plastic stretching is incompressible when J e := det(F e ) =: J.

Plastic deformation

Plastic deformation, which governs the irreversible mechanisms of the material, is considered the most dominant in amorphous solid polymers. Physically motivated [START_REF] Boyce | On the kinematics of finite strain plasticity[END_REF]-model and Anand and Ames (2006)-model represent two approaches capable of governing the development of the plastic deformation through the element (1) in Fig. 1. Because the Anand and Ames (2006)-model can predict the transient yield peak that is typical to amorphous polymers more accurately, this model is used in this study.

The evolving plastic deformation through the dashpot (1) is modeled by the following power-law type strain rate γvp = v0

τ (1) s (1) + αp 1/m0 (41)
where v0 , m 0 , and α denote the material parameters, and p = -trace(τ A /J e )/3 is the pressure. Strain softening is modeled by an internal variable as:

ṡ(1) = h 0 γvp (1 -s (1) /s (1) (ϕ)), s (1) (0) = s 0 (42) where s(1) := s sv (1 + b(ϕ sv -ϕ)) (43) with φ := g 0 γvp ( s (1) s sv -1), ϕ(0) = 0 (44)
for evolving free volume. Additional parameters needed are h 0 , s 0 , g 0 , s sv , b, and ϕ sv .

Anisotropic hardening

To model the anisotropic hardening behavior in large strains, several choices for the plastic portion of the free energy are available. The most prominent ones are based on the 8-chain and full network models, cf. [START_REF] Boyce | On the kinematics of finite strain plasticity[END_REF]; [START_REF] Wu | On improved network models for rubber elasticity and their applications to orientation hardening in glassy polymers[END_REF] and Fig. 2.

According to the Anand and Ames (2006)-model, the anisotropic deformation behavior in large strains is based on the entire deformation (spring c) in Fig. 1). However, as discussed in Introduction, this treatment can be questioned for amorphous solid polymers, and the idea that the anisotropic deformation behavior is solely given in terms of the viscoelastic-plastic response is applied in the proposed model. Based on the 8-chain modeling concept, cf. [START_REF] Boyce | On the kinematics of finite strain plasticity[END_REF], the viscoelastic-plastic portion of the free energy in the proposed model (the element c) in Fig. 1) is governed by the function:

ψ vep (λ vep ec ) = C R λ vep ec √ N ξ + ln ξ sinh(ξ) ≥ 0, ψ vep (1) = 0 (45)
wherein C R and N stand for parameters, and ξ := L -1 (λ vep ec / √ N ) means the inverse of the Langevin function, L. Applying a non-affine network stretch

λ vep ec := 1 √ 3 trace( Cvep ), ( 46
)
as has been done in [START_REF] Boyce | On the kinematics of finite strain plasticity[END_REF], the backstress Bvep representing kinematic hardening becomes

Bvep = C R 3λ vep ec √ N ξ Cvep . ( 47
)
When calculating ξ in (47), the Padé approximation is applied, cf. [START_REF] Cohen | A pade approximant to the inverse langevin function[END_REF].

Corollary 2.3. The spatial counterpart β vep of Bvep takes the form

β vep = C R 3λ vep ec √ N ξc vep , (48)
where c vep = J e-2/3 F e Cvep F e,T .

Proof. Considering the kinetic, stress measures contravariant, cf. (Belytschko et al., 2000, Box 5.16), Cvep in ( 47) is in contravariant form. Then the push-forward operation (determined in connection with (20)) of Cvep results in ĉvep := F e Cvep F e,T . Using ( 38) and ( 39) in (37) reveals that d vep is deviatoric, when

0 = trace(d vep ) = trace( Dvep ) = trace( Cvep-T Ċvep,T ) = 1 det( Cvep ) d dt det( Cvep ), i.e. det(F e-1 ĉvep F e-T ) = det(F e-1 )det(ĉ vep )det(F e-T ) = 1. It is then re- quired that ĉvep = J e+2/3 c vep

Viscoelastic deformation

Consistent with (41) for viscoplasticity, the viscoelastic deformation through the dashpot in element ( 2) is governed by

γve = v0 τ (2) s (2) + αp 1/m1 ≥ 0 (49)
where m 1 is a material parameter. Against the internal variable s (1) defined above, s (2) is regarded as a constant parameter.

The back stress influencing the viscoelastic element is defined by the function of the neo-Hookean form as follows:

ψ ve = 3 2 µ 1 (λ ve ec ) 2 -1 , (50)
where

λ ve ec := 1/ √ 3 trace( Cve ) (51)
is an effective viscoelastic stretch, and the back stress modulus µ 1 := μ1 (ϕ) is considered a function of the local free volume ϕ as

μ1 := c 1 1 - µ 1 µ 1,sat φ, µ 1 (0) = µ 0 1 (52)
involving the material parameters c 1 , µ 1,sat , and µ 0 1 . Expecting µ 1,sat ≤ µ 0 1 , µ 1 reduces to its saturation value µ 1,sat as the free-volume increases. Then, in view of (51),

∂ψ ve ∂ Cve = µ 1 2 Ī
where Ī is the identity given in the intermediate placement, gives

Bve := µ 1 Cve . (53)
Benefiting Corollary 2.3, the spatial form becomes β ve := µ 1 c ve , where c ve = J e-2/3 F e Cve F e,T . The treatment of the dissipation inequality (21) based on the constitutive equations is given in Appendix A. A numerical solution method is introduced in Appendix B. A brief summary of the approach is presented in Table 1.

Experimentation

An amorphous thermoplastic polycarbonate (PC) (Lexan 223R granulate) was used to calibrate and evaluate the proposed model. Dogbone-shaped specimens for both monotonic and cyclic loading tests were injected in accordance with the standard ASTM D638-03 with the type IV specimen, ASTM Committee (2003), [START_REF] Barriere | Modeling, simulation, and experimentation of fatigue behavior in amorphous solids[END_REF]. The geometry of the specimen is demonstrated in Fig. 3. 2. Kinetics: τ = L e : ln v e = 2µ(ln v e ) dev + κ ln

J e i, τ B = C R 3λ vep ec √ N L -1 ( λ vep ec √ N )c vep (anisotropic hardening), τ A = τ -τ B , β vep := τ B .
2.1 Micro-macro transition:

λ vep ec = 1 √ 3 ( Ī : Cvep ) ∈ [0, √ N ].
3. Viscoelastic and viscoelastic-plastic micromechanisms:

3.1 Flow rules:

d vep = d vp + d ve = γvp n (1) + γve n (2) , n (1) := τ dev A √ 2τ dev A , τ dev A := 1 2 τ dev A : τ dev A , n (2) := τ dev A -β ve √ 2τ (2) , τ (2) := 1 2 (τ dev A -β ve : τ dev A -β ve ).
3.2 Loading, unloading (viscoelasticity): γve := v0 ( Model parameters: v0 , α, s (2) , m 1 , c 1 , µ 1,sat , µ 0 1 (viscoelasticity), s 0 , m 0 , C R , N, h 0 , g 0 , b, scv, ϕcv, ϕ 0 (viscoplasticity).

τ (2) s (2) + αp ) 1 m 1 , s (2) = s (2) (0), β ve = µ 1 c ve , μ1 := c 1 (1 - µ 1 µ 1,sat ) φ, µ 1 (ϕ(0)) = µ 0 1 > µ 1,sat . 3.3 Yielding, softening (viscoplasticity): φ = g 0 ( s (1) scv -1) γvp , ϕ(0) = ϕ 0 , γvp := v0 ( τ dev A s (1) + αp ) 1 m 0 , p = - 1 3 trace(τ A ), ṡ(1) = h 0 (1 - s (1) s(1) ) γvp , s(1) (ϕ)) := ssv(1 + b(ϕsv -ϕ)), ϕ(0) < ϕcv.

Uniaxial tensile tests

The experiments were conducted under isothermal conditions (room temperature and low stress/strain rates) without a marked temperature rise (below 10 • C), i.e. the tests are considered quasi-static. An Instron Electropulse E10000 machine equipped with a 10 kN load sensor was used. An extensometer with the gauge length of 25 mm was also applied to measure the axial elongation u, see Fig. 3. The strain was defined by ǫ := u/L, where L is the gauge length of the extensometer used. The stress was defined by σ := F/A, where F is the force (measured by the machine) and the area of the initial cross section is given by A5 . The monotonic tensile tests were performed at different strain rates (displacement control) ǫ = 0.001 s -1 , ǫ = 0.01 s -1 , and ǫ = 0.1 s -1 . The acquisition frequency for these tests was 20 Hz. The force-controlled cyclic tests were conducted under sinusoidal wave form ranging between the maximum stress and minimum stress based on different stress ratios of R = 0.1 and R = 0.5 (the ratio between the smallest and largest stresses). The values of maximum stresses investigated were 97 %, 90 %, 75 %, and 50 % of the peak yield stress, 60 MPa. The data acquisition frequency for cyclic tests was 1000 Hz while the loading frequency was f = 5 Hz, i.e. 200 points per each cycle were recorded. At least two specimens were tested at each level, and the tests interrupted due to the failure of the specimens. For calibration and simulations, the proposed model and peer state-of-theart models were implemented in a computer program using the Intel Fortran application. The evolution equations used in the numerical solution method are presented in the Appendix B.

Model calibration

In the gauge section of test specimens, uniaxial loading conditions prevailed, and the model reduces to a simplified form facilitating the model calibration. The elastic model parameters were extracted from the strain-controlled, monotonic cold drawing tests. The Young's modulus E was initially calculated from the linear strain range (0 -0.001) of the σǫ curve, and then, during calibration, the range was progressively increased (E was reduced) to sharpen the predictions. The Poisson's ratio ν typical for solid PC-polymer was found from [START_REF] Holopainen | Modeling of long-term behavior of amorphous glassy polymers[END_REF]; [START_REF] Holopainen | Modeling of the mechanical behavior of amorphous glassy polymers under variable loadings and comparison with state-of-the-art model predictions[END_REF]. A detailed calibration procedure to find the remaining parameters is discussed next.

Stage 1. The stress-strain response around the peak yield stress and for subsequent softening and hardening at large deformations were first determined under a displacement-controlled monotonic tension by least squares fitting and the data is shown in Fig. 4. Initially, the viscoelastic mechanism (2) was considered constant, i.e. the initial values of both the variables and parameters were applied in ( 49) and ( 52). Moreover, s (1) was first considered constant in (41). Then, only the initial flow strength s 0 , the strain-rate-sensitivity parameter m 0 , strain rate parameter v0 , pressure sensitivity parameter α, the rubbery modulus C R , and the locking stretch N were needed to be defined. The initial parameters C R and N governing anisotropic hardening at large srains were taken from [START_REF] Boyce | On the kinematics of finite strain plasticity[END_REF] (for a PC). Other initial values were taken from (Anand andAmes, 2006, p. 1168) (for a PMMA-polymer).

Stage 2. Continuing with data given in Fig. 4, the deformation resistance s (1) was allowed to evolve along with the free volume ϕ when the parameters h 0 , g 0 , b, s cv , and ϕ cv were calibrated. Assuming that the change of ϕ is significantly small, the value ϕ cv = 0.001 proposed in [START_REF] Anand | On modeling the micro-indentation response of an amorphous polymer[END_REF] chosen. The parameter h 0 directly influences the pre-peak slope, whereas the resting parameters influence the peak yield stress and its post-peak slope, and were determined by the least squares fitting. The initial values of the parameters were extracted from (Anand andAmes, 2006, p. 1168). The parameters determined in this stage have only a minor effect to the ones that had already been determined in the previous stage for large deformations. The development of the free volume ϕ and the deformation resistance s (1) with strain is depicted in Fig. 4. The comparison of Fig. 4(right) with the macroscopic stress-strain response in Fig. 4(left) (at strain values of 0.03 -0.13) reveals that the deformation resistance influenced by the free volume is the one that strongly influences the shape of the yield peak.

Stage 3. The remaining parameters for viscoelasticity (m 1 , s (2) , c 1 , µ 1,sat , µ 0 1 ) were finally defined. Fig. 4 schematically shows the point, in terms of s (2) , as the viscoelastic micromechanism approximately starts taking effect. The point at which the viscoplastic deformations start governing the deformation behavior (in terms of s (1) ) has higher strain. Thus, an inspection of Fig. 4 indicates that the initial values of these deformation resistances should range between 10 -30 MPa, and satisfy s (1) > s (2) . The slope that connects these two points is governed by the modulus µ 1 and is, based on the experimentally observed σǫ curves, greater than the elastic modulus E for the initial response (strains 0 -0.012). The evolution of µ 1 is shown in Fig. 5. The parameters c 1 and µ 1,sat in (52) define the slope of softening and the final value of µ 1 , respectively.

To attain finer details of the experimentally observed stress-strain responses including unloadings, the remaining parameter m 1 was defined based on the assumption m 0 < m 1 . Only a few ad hoc iterations were needed to obtain its reasonable value.

The final set of model parameters is given in Tables 2 and3. For comparison, the [START_REF] Anand | On modeling the micro-indentation response of an amorphous polymer[END_REF]-model parameters for the tested PC have been determined, and the values are presented in Tables 4, 5, and6

Evaluation of the approach

Displacement-control

Based on the final set of parameters, the stress versus strain responses up to the large strain are shown in Fig. 4. For comparison, the Anand and Ames (2006)-model response is also shown. Despite a fairly reduced parameter set, the proposed model captures the experimental response well. A low post-peak slope observed in the experimental response right after the peak yield stress is due to an inhomogeneous three-axial deformation behavior exhibiting short-term necking (strain increases at almost a constant stress level).

The capability of the model in predicting the strain rate dependence is demonstrated in Fig. 6. The model can predict the increasing peak yield stresses with increasing strain rates well. At the highest strain rate, the test specimens showed a brittle fracture, immediately when the peak yield stress was reached. Further, other experimental responses up to this strain value are presented in Fig. 6. Because the model predictions represent a single material (integration) point without failure, the model responses are able to replicate significant softening followed by hardening at large strains.

Load-control

A creep-recovery test was first conducted because during cyclic loadings a part of the accumulated deformation comes from creep, and also the capacity to predict recovery implicates the capacity to predict cyclic response. The isochronous test comprises the elastic domain (the first load-unload), viscoelastic domain (second load-unload), and the viscoelastic-plastic domain (third loadunload), see Fig. 7, ASTM Committee (2001). As it has been shown, the proposed model captures the creep strains (at nonzero constant stress levels) and subsequent recovery at zero-stress fairly well. Although the model seems to exaggerate the irreversible plastic strain, the value of the residual strain at the end of the test is small.

Figure 8 shows the predicted and experimentally observed hysteresis loops, ǫ Figure 6: Experimental tensile stress-strain curves on PC at constant strain rates of 0.001 s -1 (black solid curve), 0.01 s -1 (black dashed curve), and 0.1 s -1 (black dashand-dot curve) (left). The blue colored curves implicate the model responses. Peak yield stress vs strain rate in tension (right). as the maximum stress is right below the peak yield stress. Despite the reduced set of material parameters, comparison of the responses reveals that the proposed model is even more accurate than the Anand and Ames (2006)-model. However, due to a large number of parameters needed in the Anand and Ames (2006)-model, the optimum set of calibrated parameters is difficult to find. It can also be observed that the incorporation of only one viscoelastic micromechanism in the proposed model makes possible the nonlinear loading and unloading responses with acceptable accuracy. This characteristic is due to the evolution of µ 1 shown in Fig. 5 that results in a suitable material stiffness which is primarily responsible for slightly increasing stiffness of the initial loading response as well as the nonlinear unloading behavior, [START_REF] Anand | On modeling the micro-indentation response of an amorphous polymer[END_REF].

For comparison, the prediction by the [START_REF] Boyce | On the kinematics of finite strain plasticity[END_REF]-model is also demonstrated in Fig. 8. Owing to the linear elastic constitutive description without viscoelasticity, [START_REF] Boyce | On the kinematics of finite strain plasticity[END_REF]-model is not able to predict the nonlinear hysteresis loops or ratcheting under the applied maximum stress.

The influence of the loading frequency, predicted by the proposed model, is demonstrated in Fig. 9(left). The material shows hardening, i.e. stretching reduces apparently with increasing loading frequency. This behavior is typical for many amorphous polymers, [START_REF] Boyce | On the kinematics of finite strain plasticity[END_REF]; [START_REF] Anand | On modeling the micro-indentation response of an amorphous polymer[END_REF].

The selected loops under lower maximum stress values are depicted in Figures 9(right) and 10. Hysteresis loops exhibit ratcheting (also termed cyclic creep), i.e. the mean strain grows as the applied load varies between its maximum and minimum values, cf. [START_REF] Jiang | A viscoelasticplastic constitutive model for uniaxial ratcheting behaviors of polycarbonate[END_REF]; [START_REF] Holopainen | Continuum approach for modeling fatigue in amorphous glassy polymers. applications to the investigation of damageratcheting interaction in polycarbonate[END_REF]; [START_REF] Kang | Cyclic Plasticity of Engineering Materials: Experiments and Models[END_REF]. A significant initial increase of ratcheting deformation is due to the great mean stress level used. Then the growth of mean strain gradually reduces during the fluctuating loading, which characteristic can be termed as cyclic hardening. In our test specimens, cyclic hardening was recognized as a slightly inhomogeneous deformation behavior (macro-scopically manifested as initial necking), i.e. stretch initially increased at fixed maximum/minimum tensional stress levels, but rapidly attained its stabilized final value.

To demonstrate the ratcheting behavior more precisely, the ratcheting (mean) strain is defined by:

ǫ r := 1 2 (ǫ min + ǫ max ) (54)
wherein ǫ min and ǫ max mean the minimum strain and the maximum strain in each cycle. Development of ǫ r is demonstrated in Fig. 11(left). Under the high mean stresses applied, the proposed model slightly overestimate the ratcheting strain as the number of cycles N < 2000. Considering the highest maximum stress, the model response also seems to differ from the experimental curve in the intermediate region as N =2000-6000. The difference in this region is due to an experimentally observed temperature peak resulting in slightly nonisothermal conditions and thus, increased ratcheting strain, Fig. 12. Otherwise, the temperature increase remains low, and the proposed model well captures the experimentally observed ratcheting deformation at different stress levels.

Using the parameters given in Tables 4,5 Table 7, no difference with the Anand and Ames (2006)-model response was observed, cf. (Jiang et al., 2015, Fig. 5).

The loops and the development of the ratcheting strain under different maximum stresses for R = 0.5 are demonstrated in Figs 12, 13, and 14, respectively. Because the temperature now does not influence the mechanical behavior (the rice remains under 5 • C), the proposed model is able to predict the ratcheting strain with reasonable accuracy, especially when compared to the Anand and Ames (2006)-model, Fig. 14.

Concluding remarks and further research avenues

The article introduces a compact viscoelastic-viscoplastic constitutive model for solid polymers capable of improving the predictions of the cyclic deformation behavior at stress levels below the peak yield stress. It is also evidenced that the model is accurate under very large strains, as monotonic loadings are applied. The proposed model is based on the celebrated Haward and Thack- -a reduced set of internal state variables and material parameters need to be defined; -elastic portion of the deformation is clearly described by a linear spring separate to viscoelastic-plastic elements according to the classical Haward and Thackray (1968) -model; -the nonlinear Langevin spring is defined solely using the viscous (viscoelasticplastic) deformation in accordance with the micromechanics-based Boyce et al.

(1989)-model; -purely plastic and viscoelastic deformation are defined by the specified stress controlling the corresponding micromechanisms.

From a number of uniaxial tension experiments, it is shown that this compact formulation, in combination with the calibration procedure can reasonably predict the shape of the hysteresis loops and long-term ratcheting behavior. The comparison with [START_REF] Anand | On modeling the micro-indentation response of an amorphous polymer[END_REF]; [START_REF] Jiang | A viscoelasticplastic constitutive model for uniaxial ratcheting behaviors of polycarbonate[END_REF] -models indicates that, despite the larger number of parameters needed in these distinguished models, the proposed model is at least as accurate. Creep-recovery tests evidenced also the capability of the model. Due to the compact formulation, the proposed model is easy to implement and can be used as a built-in feature in different finite-element packages.

Although the uniaxial experimentation has been comprehensive, the proposed model needs to be evaluated under cyclic shear/torsion and arbitrary multidimensional loadings. It should be noted that the research concerning the experimentation of multi-axial cyclic plasticity is practically nonexistent in lit-erature. Considering the high loading rates, the effect of temperature rises and needs to be encompassed in the model. Furthermore, an important problem is the modeling of fatigue damage (incl. microstructural changes) that supposedly develops under the applied load levels.

proposing a spin W vep which is algorithmically consistent, i.e. it becomes finally skew-symmetric at the end of the integration step, [START_REF] Holopainen | Modeling of long-term behavior of amorphous glassy polymers[END_REF]; [START_REF] Engqvist | Modelling multi-scale deformation of amorphous glassy polymers with experimentally motivated evolution of the microstructure[END_REF]; [START_REF] Holopainen | Continuum approach for modeling fatigue in amorphous glassy polymers. applications to the investigation of damageratcheting interaction in polycarbonate[END_REF]. Then, supposing the current deformation gradient F is known (as it does in standard finite-element implementations), the variables needed to be solve become: Y := [F e Cve W vp s (1) µ 1 ϕ].

(B.1)

The deformation measure Cve defines the backstress needed in the flow-type rule (38) and its evolution is defined in (27). Moreover, the viscoelastic spin is considered zero. Consistent with ( 27), but noting the evolution of viscoplastic deformation can include a spin, the evolution of Cvp can be defined by Ċvp := Lvp Cvp + Cvp Lvp,T . However, this strain measure is optional, and only used for post plotting purposes, see (39) for viscoplastic flow. Once F e is known, F vep is solved from the decomposition (1).

Based on an implicit Euler scheme, the exponential update of F vep in accordance with [START_REF] Miehe | Associative multiplicative elasto-plasticity formulation and aspects of the numerical implementation including stability analysis[END_REF] is used, i.e. where Lvep = Dvep + W vp . To simplify notation, n + 1 for indicating the updated state is omitted, whereas n is used to refer to the known state t n .

To evaluate W vp , its skew-symmetry-property as well as the symmetry of F e at the end of the integration interval are used. Finally, the following equations must to be solved: (B.4) The system (B.3) is well suitable for standard finite-element codes, [START_REF] Holopainen | Modeling of long-term behavior of amorphous glassy polymers[END_REF]; [START_REF] Engqvist | Modelling multi-scale deformation of amorphous glassy polymers with experimentally motivated evolution of the microstructure[END_REF]. Therein, the algorithmic tangent stiffness (ATS) tensor is of major importance. The procedure to obtain the ATS tensor can be found from (Holopainen and Wallin, 2012, p. 8). An initial simulation result based on the system (B.3) is shown in Fig. 3.

Reduction to uniaxial stress

The evolution equations for uniaxial loadings used in the calibration and simulations are defined below. The total deformation gradient under such conditions is expressed by: F = v 11 e 1 ⊗ e 1 + v 22 e 2 ⊗ e 2 + v 33 e 3 ⊗ e 3 , (B.5)

where {e i , i = 1, 2, 3} are the unit base vectors and v ii F ii denote the principal stretches having the relation J = det(F ) = v 11 v 2 22 > 0 with the rotation R being unity. The only nonzero stress component is given by:

σ 11 = 1 J E

  d e := sym(l e ), d vep := sym(l vep ), Dvep := sym( Lvep ) (11) and ω := skew(l), ω e := skew(l e ), ω vep := skew(l vep ), W vep := skew( Lvep ), (

τ

  : l e → τ : (d e + QQ T ) + T : QT Q = τ : d e , τ : l vep → τ : d vep + T : QT Q = τ : d vep , (20) i.e. τ : l = τ : (l e + l vep ) = τ : (d e + d vep ) = τ : d holds. The stress T := F e-1 τ F e-T in (20) is the counterpart of τ given in the intermediate placement 3 .

ψ

  = ψ(b e , Cvep , Cve , ϕ) = ψe (b e ) + ψvep ( Cvep ) + ψve ( Cve , ϕ)

  ), noting the symmetry of b e , and considering ψ e as an isotropic function of b e result in ∂ψ e ∂b e : ḃe = 2 ∂ψ e ∂b e b e : d e .

  ) that 0 = trace(d ve ) = trace( Dve ) = trace( Cve-T Ċve,T

Figure 2 :

 2 Figure2: Idealized chain structure according to the 8-chain model (left) and the full network model (right). The dimension of the 8-chain cube is given by r0 and φ0, and R0 denotes the radius of the micro-sphere.

  1. Kinematics: F = F e F vep , b e = F e F e,T , v e = √ b e , Cvep := F vep F vep,T , c vep = J e-2/3 φ * ( Cvep ), Cve , c ve = J e-2/3 φ * ( Cve ), Dve = φ * (d ve ), d vep = d vp + d ve , Dvep = φ * (d vep ), Lvep = Dvep + W vep .

3. 5

 5 Evolution of kinematical variables: Ḟ vep := Lvep F vep , Ċve := Dve Cve + Cve Dve .

Figure 3 :

 3 Figure 3: Flowchart for the test program comprising design of the test specimens, uniaxial tensile test set up, and simulations.

Figure 4 :

 4 Figure4: Stress-strain curves (using the final set of material parameters) in tension on PC at the constant strain rate of 0.001 s -1 (left). The markers indicate the points as the viscoelastic (black circle) and viscoplastic (blue circle) micro-mechanisms (2) and (1) start displaying a notable effect. Evolution of the free-volume ϕ for PC according to the proposed model (blue curve) and the[START_REF] Anand | On modeling the micro-indentation response of an amorphous polymer[END_REF]-model (red curve).

Figure 5 :

 5 Figure5: Evolution of the shear resistance s(1) for PC according to the proposed (blue curve) and the[START_REF] Anand | On modeling the micro-indentation response of an amorphous polymer[END_REF]-model (red curve) based on the final set of calibrated parameters. Corresponding evolutions of the stiffness µ1 (left). For the[START_REF] Anand | On modeling the micro-indentation response of an amorphous polymer[END_REF]-model, also the stiffnesses µ2 (dashed curve) and µ3 (dashedand-dotted curve) for additional viscoelastic micro-mechanisms are shown.

  . Moreover, a comparison is made with the Jiang et al. (2015)-model that is a straightforward extension of the Anand and Ames (2006)-model. Its additional parameters are presented in Table 7. Compared to the Anand and Ames (2006)-model and Jiang et al. (2015)-model, each requiring 28 and 34 material parameters, only 19 parameters are required in the proposed model.

Figure 7 :

 7 Figure 7: Loading for the creep-recovery test (left). The corresponding strain responses according to the experiment (black) and the proposed model (blue) (right).

Figure 8 :Figure 9 :

 89 Figure 8: First hysteresis loops of the Anand and Ames (2006)-model vs data with the stress ratio R = 0.1 and the frequency f = 5 Hz (left). Comparison of the loops by the proposed model, data, and Boyce et al. (1989)-model predictions (right).

Figure 10 :

 10 Figure 10: 30th and 7000th hysteresis loops. The black and blue colors imply the experimental data and predicted results according to the proposed model.

Figure 11 :

 11 Figure 11: Ratcheting strain responses with R = 0.1 and the frequency f = 5 Hz. The black, blue, and red colors imply the experimental data and predicted results according to the proposed model and the Anand and Ames (2006)-model, respectively (left). Corresponding hysteresis loops by the Anand and Ames (2006)-model (right).

Figure 12 :R

 12 Figure12: The temperature rise during the loading as the maximum stress is 54 MPa (left). 30th, 3000th, and 10,000th hysteresis loops (right). The black and blue colors imply the experimental data and predicted results according to the proposed model.

Figure 14 :

 14 Figure 14: Ratcheting strain responses with the stress ratio R = 0.5. The colors black and blue imply the experimental data and predicted results according to the proposed model (left). The Anand and Ames (2006)-model predictions are also shown (red color). Corresponding predictions of the hysteresis loops by the Anand and Ames (2006)-model (right). The loops up to 3000 cycles have been shown.

F

  vep = exp(∆t Lvep )F vep n , (B.2)

R

  1 : = F e -F F vep-1 n exp(-∆t Lvep ), R 2 : = Cve -Cve n -Dve Cve -Cve Dve , R 3 : = skew(F e ), R 4 : = sym( W vp ),R 5 : = (s (1)s (1) n -∆t ṡ(1) )/s 0 , R 6 : = (µ 1µ 1,n -∆t μ1 )/µ 0 1 , R 7 : = ϕϕ n -∆t φ, (B.3)where the first system R 1 encases (B.2) in the decomposition (1). The residuals R 2 and R 4 in (B.3) consist of six equations and R 3 of three equations (27 state variables in all).Solution of (B.3) is obtaioned by the Newton-Raphson method. The solution becomesY n+1 = Y n + ∆Y , where ∆Y = -J -1 R, J := ∂R ∂Y (Jacobian), R := [R 1 R 2 R 3 R 4 R 5 R 6 R 7 ].

Table 1 :

 1 Summary of the proposed model.

Table 2 :

 2 Elastic and viscoelastic model parameters for a PC adapted to a cold drawing test and a cyclic test using the proposed model.

	Parameter Unit .........	E MPa	ν	v0 s -1	α	s (2) MPa	m 1	c 1 • 10 -6 MPa	µ 1,sat MPa	µ 0 1 MPa
	Value ......... 2100	0.37	0.038	0.204	11	0.19	4.5	1100	3500

Table 3 :

 3 Viscoplastic parameters using the proposed model (ϕ(0) = 0).

	Parameter	s 0	m 0	C R	N	h 0	g 0	b	scv	ϕcv
	Unit .........	MPa		MPa		MPa			MPa
	Value .........	28.5	0.055	14.0	1.65	4500	0.017	755	27.5	0.001

Table 4 :

 4 Elastic and viscoelastic parameters for the PC using the Anand and Ames (2006)-model. Moreover, m2 = m3 = m1.

	Parameter	E	ν	v0	α	s (1)	s (2)	s (3)	m 1
	Unit .........	MPa		s -1		MPa	MPa	MPa
	Value ......... 2100	0.37	0.038	0.204	11	17	21	0.19

Table 5 :

 5 Viscoelastic model parameters using the[START_REF] Anand | On modeling the micro-indentation response of an amorphous polymer[END_REF]-model (continued).

	Parameter Unit .........	c 1 • 10 -6 MPa	c 2 • 10 -6 MPa	c 3 • 10 -6 MPa	µ 1,sat MPa	µ 2,sat MPa	µ 3,sat MPa	µ 0 1 MPa MPa µ 0 2	µ 0 3 MPa
	Value .........	4.5	1.8	1.3	1100	600	250	3500	3500	3900

Table 6 :

 6 Viscoplastic parameters for the PC using the Anand and Ames (2006)-model (ϕ(0) = 0).

	Parameter Unit .........	s 0 MPa	m 0	C R MPa	N	h 0 MPa	g 0	b	scv MPa	ϕcv	c 0 • 10 -6 MPa
	Value .........	28.0	0.055	14.0	1.65	16500	0.028	860	27.0	0.001	0

Table 7 :

 7 Additional viscoelastic constitutive model parameters for the PC using the Jiang et al. (2015)-model.

	Parameter	s 1,sat	s 2,sat	s 3,sat	d 1	d 2	d 3
	Unit .........	MPa	MPa	MPa	MPa	MPa MPa
	Value .........

  ln(v e 11 ) (B.6) where E is the Youngs modulus and v e 11 is the principal stretch of v e defined in (2). Owing to the uniaxial stress state, the elastic deformation is constrained by e 33 ) = ln(v e 22 ) = -ν ln(v e 11 ) (B.7) where ν is the Poissons ratio. Due to the symmetry, v 33 = v 22 in (B.5), when it follows from the decomposition (1) and (B.7) that F vep 33 = F vep 22 . Since F vep and F e are diagonal, Cvep 33 = Cvep 22 in (3). Then, noting volume-preserving viscoelastic-plastic deformation, the relation Cvep 22 = (F vep 22 ) 2 = 1/F vep Because the uniaxial stress state is inherently spin less and results in F e being symmetric, the tensor valued equations in the system (B.3) reduces to Considering displacement-control: once Cvep 11 and thus, F vep 11 are known, F e 11 can be solved from the decomposition (1) during iterations, i.e. F e 11 = F 11 /F vep 11 . When force-control is applied, F e 11 = v e 11 is directly available from the constitutive equation (B.6). Then, once Cvep 11 and thus, F vep 11 are known, F 11 is available from the decomposition (1).

	11	(B.8)
	holds. Similarly, based on (40) also Cvp 33 = Cvp 22 = 1/ Cvp 11 and Cve 33 = Cve 22 = 1/ Cve 11 in (5) are isochoric strains.
	R 1 : = Cvep 11 -Cvep n,11 -2 Dvep 11 R 2 : = Cve 11 -Cve n,11 -2 Dve 11 Cve Cvep 11 , 11 .	(B.9)

The term placement means a certain configuration of a continuum body P, i.e. the placement is given by S = k(P) under the configuration k : P → S, cf.(Zienkiewicz et al.

, 2014, p. 449). Examples are the reference and current placements, B and b.

Owing to path dependency of the viscoplastic deformation, observation of stretching and rotations is challenging and thus, the polar decomposition of F vep is not used.

It then follows that the operator φ * = F e-T (•)F e-1 in (8).

Evolution of the local free-volume is considered an important source for the transient peak associated with yielding of amorphous polymers,[START_REF] Anand | On modeling the micro-indentation response of an amorphous polymer[END_REF].

The measured stress is actually the 1st Piola Kirchoff stress,[START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF]; the Kirchhoff stress applied in the model is τ := σλ, where σ is the Cauchy/true stress and λ = 1 + ǫ is the axial stretch. However, the difference between the Kirchhoff, Cauchy, and 1st Piola Kirchoff stress is small under the strains ǫ < 0.3.

Evolution of the viscoelastic stretching Ċve (and thus, Cve ) is considered independent on the free volume ϕ which is meant for the viscoplastic stretching in accordance with its evolution law (44). Actually, change in the free-energy due to the free-volume change can be shown to be small, cf.(Anand and Ames, 2006, p. 1159).

Acknowledgements

The authors acknowledge the technician Ly Phanousith from FEMTO-ST institute for conducting all the experiments and sharing his expertise in the field of measurements. We want to show also thanks to Professor Jean-Luc Bouvard (Mines-Paris Tech., France) for his valuable comments.

Appendix A. Dissipation inequality vs specific constitutive equations Using (33) in ( 21) for the dissipation inequality, yields

which needs to be positive or zero. Using the function ( 50), ∂ψ ve /∂ϕ φ = 3 2 μ1 (λ ve ec ) 2 -1 6 . Since μ1 is negative or zero (see (52) and Fig. 5), (λ ve ec ) 2 -1 must be positive or zero in order to get -∂ψ ve /∂ϕ φ ≥ 0. Using (51) for λ ve ec , gives trace( Cve ) ≥ 3. It can be shown (using principal stretches and noting incompressibility, det( Cve ) = 1 by ( 40)) that this condition is always valid.

It remains to show that τ A : d vepβ ve : d ve ≥ 0. Using (37), (38), and (39), 2) tr β ve (τ A -β ve ) +tr τ A β ve (A.1) where the abbreviation tr means trace. Since γvp , γve , τ (1) , and τ (2) are positive, it is sufficient to show that the last term in (A.1) is less than the sum of other terms. Since trace

which is less than the sum of other terms on r.h.s of (A.1).

Appendix B. Numerical solution method

System of the equations To define the intermediate placement uniquely, F e is regarded as symmetric, see [START_REF] Boyce | On the kinematics of finite strain plasticity[END_REF]; [START_REF] Holopainen | Modeling of long-term behavior of amorphous glassy polymers[END_REF]. This implies that the visco-plastic spin W vp is generally nonzero. This restriction is adjusted by


Ting-Ting Wang

Yan-Feng Wang

Yue-Sheng Wang

Vincent Laude
email: vincent.laude@femto-st.fr

Complex Resolvent Band Structure of Phononic Crystals

We examine the extension of the definition of the resolvant band structure (RBS) to complex wavenumbers and complex frequencies. The RBS is obtained by applying a stochastic excitation to the unit-cell of a crystal and monitoring the response in dispersion space. It is found that the poles of the RBS appear in complex dispersion space around the complex band structure.

Obtaining the band structure (BS) of a particular phononic or sonic crystal, including all its peculiarities, has been a central problem in the field of phononics from the start [1][START_REF] Laude | Phononic crystals[END_REF] . Indeed, the band structure is a representation of the dispersion relation for Bloch waves, giving possible combinations of the wavevector k and of the angular frequency . For a closed and bounded domain, and in the absence of loss,  is a real function of real k and the (real) band structure is easily defined and obtained by solving an eigenvalue problem. When loss is added, complex-valued (k) [START_REF] Hussein | [END_REF] or k() 4 band structures (CBS) can be defined by extending the domain of definition of the wavevector and/or of the frequency to the complex plane. The question of which quantity should be made complex -or even if both should be made complex at the same time -remains largely open 5 .

Recently, the concept of the resolvent band structure (RBS) was introduced 6 . With the RBS, one abandons the idea of obtaining a functional relationship in the form of bands and instead proceeds to map the resolvent set over the full dual (k,) space. The resolvent set is by definition the complement of the spectrum of eigenvalues: it is composed of the whole complex plane with the isolated eigenvalues composing the spectrum removed. In practice, the RBS can be obtained easily by considering a spatially random source distributed in the phononic crystal. Indeed, the random source will excite all modes of vibration and it then suffices to map the response of the crystal as a function of w and k.

The RBS was previously obtained for real k and  and it was shown that it could adequately take into account material loss and radiation at infinity in open structures such as waveguides 6 . Nothing, however, prevents against letting k and w venture freely in the complex plane. This is in sharp contrast with the classical BS and the CBS. In this paper, this freedom is exploited to verify that the RBS gives dispersion information fully consistent with the complex k() and (k) band structures, in a case where they all can be directly compared.

PHONONICS-2019-0074

We consider a 2D sonic crystal composed of a periodic repetition of square resonators grafted along an air tube (as in chapter 2 of Ref. [2]). The unit-cell of the sonic crystal is closed and bounded: it thus admits well-defined bands. It also supports both Bragg and locally-resonant band gaps. Frequency-dependent material loss is added by making the bulk modulus of air complex. Figure 1 shows the RBS in case the wavenumber is allowed to live in the complex plane while the frequency is kept real, first in the absence of loss. For display reasons, a 3D representation of the RBS is shown on a cube. It is observed that the stochastic response concentrates exactly along the complex k() bands. Next, we consider real k and let  wander in the complex plane; such a choice is particularly desirable in case material loss is considered, since one is interested in the attenuation measured by the imaginary part of the frequency. It is found again in figure 2 that the stochastic response concentrates exactly along the complex (k) bands also presented in the figure. As a conclusion, it is found on the example of simple sonic crystal that the complex resolvant band structure identifies with both the k() and (k) complex band structures. The latter two CBS can be obtained exactly in the case considered, which is not necessarily true for all phononic and sonic crystals, especially if their unit cell is infinite and thus radiation at infinity is involved. The complex RBS is a useful tool to obtain the full complex dispersion relation of Bloch waves and it applies to a much larger range of periodic structures than band structure methods based on eigenvalue solvers.

Figure 1 :

 1 Figure 1: (a-b) Complex resolvant band structure shown for complex wavenumber and real frequency in the case of a lossless sonic crystal of rectangular resonators grafted along an air tube. The model is 2D. a is the lattice constant. (c) Stochastic excitation of the sonic crystal.

Figure 2 :

 2 Figure 2: (a-b) Complex resolvant band structure shown for real wavenumber and complex frequency in the case of a lossy sonic crystal of rectangular resonators grafted along an air tube.