Yves Rémond

M H Kaboli

S Ghobadi

A Riasi

M Zamanifar

Y Taghipour

J Pourhosseini

H Kermani

R Y Khaksar

R Naj

M Vaghari

H Jafari

H Hesni

R S Saloot

S Bagherzadeh

A Hadi

M Mehdipour

S Tamimi

A Rahimi

P Pooyan

M Safdari

A Rajabpour

E Zamanpour

M Siavashi

M Moghadami

M Doostmohammadi

M M Heyhat

A Parvizi

N Bahlouli

S M'guil

J P M Correia

Ch D George

N Wagner

M Barth

S Nierenberger

Ch Lhadi

Dh Bernard

S Kumar

K Umale

W Wang

A K Wen

R Mossi Idrissa

R Matadi

D Elanwar

M Martin

S Essa

A S Isel

Abstract

Interphase layer in semi-crystalline polyethylene has been the least known constituent of this widely used polymer, in terms of the mechanical properties. Because of the metastable nature and nanometric size of the interphase and its intimate mechanical coupling to the neighboring crystal and amorphous domains, detailed characterization of its mechanical properties have eluded any experimental evaluation. Mechanical characterization of the interphase layer in polyethylene is the focus of two major technical chapters of this dissertation. The characterization scenarios are deployed through applying the micromechanics and continuum mechanics relationships to the relevant atomistic simulation data. The third technical chapter deals with the large, viscoplastic deformation simulation of an aggregate of polyethylene using a multiscale, homogenization analysis.

Elastic characterization of the interphase layer is implemented by applying the relationships of two distinct micromechanical homogenization techniques to the Monte Carlo molecular simulation data available for the interlamellar domain. The micromechanical approaches consist of "Extended Composite Inclusion Model" and "Double-Inclusion Method". The atomistic data, on the other hand, includes the variation of the interlamellar stiffness as well as the amorphous and interphase thicknesses with temperature for 350-450 K. To implement this characterization, the temperature dependence of the amorphous stiffness is also required, which is established using the relevant findings from the literature. The interphase stiffness is successfully dissociated form that of the interlamellar domain using the abovementioned micromechanical techniques, whose results match perfectly. Interestingly and contrary to conventional materials, the interphase stiffness lacks the common feature of positive definiteness, which indicates its mechanical instability. An ad hoc sensitivity analysis is worked out to assess the impact of the existing uncertainties on the dissociated results. Moreover, the effective Young's modulus of the interphase is evaluated using its dissociated stiffness, which compares well with the effective interlamellar Young's modulus for highly crystalline polyethylene, reported in an experimental study. This satisfactory agreement along with the identical results produced by the two micromechanical approaches confirms the new information about the interphase elastic properties and endorses the proposed dissociation methodology to be applied to similar problems.

Hyperelastic characterization of the interlamellar domain and the interphase layer in polyethylene is undertaken in the next chapter. When polyethylene undergoes large deformations, its interphase layer together with the amorphous phase behaves hyperelastically. The proposed hybrid algorithm consists in applying the constitutive relations of an isotropic, compressible, hyperelastic continuum to the molecular dynamics simulation results of a polyethylene stack. Assuming a neo-Hookean model, the governing equations are derived, using which the hyperelastic parameters for the central amorphous phase, the interphase layer, and the interlamellar domain are identified with the help of the optimization notion and a set of nonnegative objective functions. The identified hyperelastic parameters for the interlamellar domain are in good agreement with the ones that have been estimated experimentally and are frequently used in the literature for the noncrystalline phase. The specifically developed uncertainty analysis indicates that the shear moduli are identified with a higher degree of certainty, in contrast to the bulk moduli. It is also revealed that the presented continuum mechanics analysis is able to capture the melting, recrystallization and rotation of crystalline chains that take place during the deformation. The by-product of this identification analysis is the estimation of the evolution of the interphase boundaries that might not be estimated otherwise.

The last chapter is devoted to reexamining the large, viscoplastic deformation of an aggregate of polyethylene. The novelties of this re-examination lie in the model adopted for the back stress of the noncrystalline phase, in correcting the projection tensor proposed for the modified Taylor approach, and also in enforcing the inequality constraints that used to be neglected in the previous studies. Properly speaking, in this re-examination, the Gent model is used for the back stress of the noncrystalline domain in lieu of the 3-chain model or 8-chain model that has previously been used. Additionally, a projection tensor is proposed for the modified Taylor approach which is more complete that the one available in the literature. Furthermore, in the newly developed computer code, the constraint that the resolved shear stresses in each slip system must not exceed the corresponding critical shear resistance is respected. And above all, the outstanding feature of this chapter is the proposed "optimization-based" methodology which is realized through defining a set of non-negative objective functions and finding their minima in the presence of a set of associated constraints.

v Figure 2. A) An electron micrograph of a chlorosulphonated section of a medium-density polyethylene with an insert cut through a computer-built spherulite [15], B) Schematic model of secondary crystalline branches (twigs) grown on a mother crystalline branch [16]. .......................... 7 Figure 3. A) Simulated fully grown spherulite [16], B) Spherulites, as they start to appear and grow in polyethylene [17] Figure 5. A) Unit cell model introduced in [START_REF] Tomita | Computational characterization and evaluation of deformation behavior of spherulite of high density polyethylene in mesoscale domain[END_REF] for multiscale modeling of polyethylene; B) Schematic of multiscale modeling of PE introduced in [START_REF] Van Dommelen | Micromechanical modeling of intraspherulitic deformation of semicrystalline polymers[END_REF]; C) The building block of an aggregate of polyethylene along with the representation of local crystalline and convected axes [START_REF] Van Dommelen | Micromechanical modeling of the elasto-viscoplastic behavior of semi-crystalline polymers[END_REF]; D) The relative positioning of the major constituents of polyethylene microstructure on which the authors in [START_REF] Hütter | Monte Carlo Simulations of Semicrystalline Polyethylene: Interlamellar Domain and Crystal-Melt Interface[END_REF] based their atomistic model; E) schematic of a spherulite + the composite inclusion model often adopted in the literature [START_REF] Lee | Simulation of large strain plastic deformation and texture evolution in high density polyethylene[END_REF] κ into the first relationship of Eq.( 30) is superposed.

List of Figures

The noisiness of the lateral stress computed by Eq.( 30) can be traced to the simulation data for 4. The slip systems of the crystallite in high density polyethylene along with their respective slip resistances adopted here [START_REF] Van Dommelen | Micromechanical modeling of the elasto-viscoplastic behavior of semi-crystalline polymers[END_REF][START_REF] Lee | Simulation of large strain plastic deformation and texture evolution in high density polyethylene[END_REF][START_REF] Lee | Micromechanical modeling of large plastic deformation and texture evolution in semi-crystalline polymers[END_REF]

General Introduction

POLYETHYLENE is the most applied polyolefin with a globally annual production of more than 80 million metric tons [1]. Polyethylene belongs to the family of polyolefins, which constitute the most widely used commodity polymers in the world with a worldwide production in excess of 120 million tons per year [2]. As compared to, say polypropylene, which is the next important member of the family of polyolefins, the volume of production of polyethylene is about twice larger. In everyday life, plastic is the most commonly used term for polyethylene. Shopping bags, shampoo bottles, children toys and different parts of irrigation systems are well known examples, whereas bullet proof vests [3], total knee replacement [4] and total hip arthroplasty [5] are less known examples of polyethylene applications. Such ubiquitousness is primarily due to its favorable resistance to physical loads as well as the inertness to most chemicals, which is offered at an economic price.

It was the outbreak of World War II that unveiled the significance of polyethylene. Owing to the huge material demand created in the wake of the war, much of the supply of polyethylene has been engineered for applications in radar technology, due to its excellent dielectric properties, and in the second place, in submarine communication cables [6]. In the aftermath of the war, much of the output of polyethylene in Britain was mostly directed to injection molded household items, while on the other side of the ocean it opened its way through packaging industry and film market [6]. Interesting properties of polyethylene led to the production of new grades and consequently new range of applications including Tetra Pak® containers for liquid foodstuff, impermeable sheets in agriculture and construction, thin-walled bottles, water piping, corrosion resistive tanks, and so forth, thanks to its numerous salient advantages. Nevertheless, low tensile strength, high flexibility, low softening temperature, tendency to creep, and environmental stress cracking are some of the basic drawbacks of polyethylene that have hindered its wider applications.

Potential applications of various grades of polyethylene, though, depend on the underlying microstructure. Polyethylene of any grade is partly crystalline and partly amorphous, hence the epithet "semi-crystalline". Furthermore, since 1950s there have been arguments that there exists in this biphase, polymer composite, so to speak, a third phase with intermediate properties that separates the fully amorphous phase from the crystalline one. This third phase, called "interphase layer" here and later on, has been known to be stiffer than the central amorphous phase but had eluded any characterization attempt. In other words, the mechanical properties of the crystalline and central amorphous phase have been identified by means of a variety of techniques during recent decades but the mechanical properties of the interphase layer had remained unknown due chiefly to its nanometric thickness as well as its intimate anchorage to the adjacent phases. Neither the state-of-the-art experimental techniques nor the atomistic simulation approaches, which are known as virtual laboratories, have been able to provide us with numerical values of the mechanical properties of the interphase layer. This task has been undertaken for the first time in this thesis where the elastic and hyperelastic properties of the interphase layer are successfully identified by employing the relationships from the micromechanics of heterogeneous materials and continuum mechanics in conjunction with the atomistic simulation data. The significance of this characterization is better appreciated when one wishes to estimate the effective properties of any grade of polyethylene, as a heterogeneous solid, which is a direct function of the constitutive properties of its constituents as well as their volume fraction.

The focus of the present dissertation is the mechanical characterization of the interphase layer in polyethylene together with the predication of the effective response of an aggregate of polyethylene subjected to large deformations, from a multiscale homogenization point of view. To this end, the microstructure of polyethylene, or strictly speaking, high density polyethylene is surveyed in the first chapter. A quick introduction to the chemistry of polyethylene and to the various production processes is also given, for completeness purposes, followed by a presentation of the most wellknown grades of polyethylene along with some of their features. After a brief discussion about the crystalline lamella of polyethylene, the non-crystalline phase and in particular the interphase layer is elaborated in details. A general review of the relevant literature is also provided in this chapter.

In the second chapter, a combinatory methodology is presented for the elastic characterization of the interphase layer. The methodology consists in applying two micromechanical homogenization techniques to Monte Carlo molecular simulation results available for the noncrystalline domain of polyethylene. The molecular simulation results include the variation of the tensorial components of the interlamellar domain together with the variation of the interphase thickness with temperature for the temperature range 350-450 K. The two homogenization techniques employed here provide analytical relationships in tensorial form. For the implemented dissociation analysis, the stiffness tensor of the amorphous phase is also required which is established based on the findings of different experimental and theoretical studies. The presented analysis leads to the useful by-results that are as valuable as the main findings of this chapter. Finally, the dissociated results are indirectly verified using a different homogenization technique and by comparison with the experimental data from the literature.

In the third chapter, the hyperelastic characterization analysis of the interphase, amorphous and interlamellar domains in semi-crystalline polyethylene is presented in detail. The deployed methodology is again a hybrid one: the governing constitutive equations of an isotropic, homogeneous, compressible, hyperelastic continuum are applied to molecular dynamics simulation data of a polyethylene stack. The identified hyperelastic parameters compare very well with the corresponding values determined experimentally for the interlamellar domain, indicating successful implementation of the proposed identification scenario. It is revealed that the presented continuum analysis can successfully capture the physical phenomena such as melting/re-crystallization that take place at the crystalline/noncrystalline interface, or the reorientation of the chain stems in the crystallites that are observed during the molecular simulation. Similar to the preceding chapter, the hyperelastic characterization analysis is accompanied by interesting by-results, including the evolution of the amorphous and interlamellar domains, which are equally valuable.

As a complementary part, the fourth chapter is dedicated to the analytical and numerical re-examination of the large deformation simulation of an aggregate of polyethylene. One of the objectives of this chapter is to compare different back stress models available in the literature that have been proposed for capturing the effect of molecular alignment of polymer chains in the noncrystalline domain under large deformations. As another objective, a complete, modified Taylor homogenization approach, which is specifically adapted for the large deformation simulation of polyethylene crystallites, is established based on the physical constraint of the inextensibility of the chain stems in the crystallites. This new approach together with the Sachs approach is used for the prediction of diagrams of equivalent stress vs. equivalent strain as well as the pole figures of polyethylene samples under different modes of deformation. This chapter contains some preliminary results and further results will be published subsequently.

Finally, chapter five recapitulates the novelties of this study, the contributions made, and also the conclusions drawn. Several practical suggestions are put forward as well for the follow-up work. Appended to the end of this dissertation are nine useful appendices that contain materials relevant to the discussions of the manuscript whose direct inclusions in the main body could obscure its readability. The most interesting appendices are Appendix A, which covers the mathematical details of calculating the orientational (volume) average of a fourth-order tensor, and Appendix E where an innovative algorithm for calculating the inverse of any fourth-order tensor that possesses minor symmetry is elaborated.

Semi-crystalline polyethylene

Chemistry of semi-crystalline polyethylene 1.1 SEMI-CRYSTALLINE POLYETHYLENE is chemically known as the simplest organic linear macromolecule, which is produced from ethylene (CH2=CH2), a colorless, flammable gas, with the IUPAC name "ethene". Upon polymerization, the double bonds between carboncarbon atoms are broken and (-CH2-)n polymer chains are formed. Therefore, in polyethylene, the base monomer is methylene,

H H | | C
--, hence the IUPAC designation "poly(methylene)" for the commercially named polyethylene. In a molecular chain the non-extreme carbon atoms are covalently bonded to each other with a pair of hydrogen atoms linked to each carbon. Chain ends, however, are terminated by methyl groups.

Different configurations of molecular chains in polyethylene give rise to different grades of polyethylene. Of various types of polyethylene, the principal varieties are high density polyethylene (HDPE), low density polyethylene (LDPE) and ultra-high molecular weight polyethylene (UHMWPE). Should some of the hydrogen atoms in the generic molecular chain of polyethylene be replaced with polyethylene chains, the resulting polyethylene is called branched polyethylene or low density polyethylene. The molecular weight in LDPE is around 25000. In case there are no branches on the backbone chain, the resulting polyethylene will be linear polyethylene or high density polyethylene. Low level of branching in HDPE is achievable by an appropriate choice of catalyst. The molecular weight in HDPEs ranges from 2-5×10 5 . When the molecular weight varies from 2×10 6 to 6×10 6 , a new grade of polyethylene known as ultra-high molecular weight polyethylene is obtained. The very long chains in UHMWPE are responsible for its comparatively more efficient load bearing and more favorable mechanical properties.

Production processes of semicrystalline polyethylene 1.2

Although there were sporadic reports of the production of polymerized ethylene [7][8][9], its commercial importance was not appreciated until 1935 when Perrin developed a reproducible process for the polymerization of ethylene whose patent was taken out by the British company of Imperial Chemical Industries [10]. This type of polyethylene was what we know today as low density polyethylene. The significance of high density polyethylene, on the other hand, was highlighted by the work of Ziegler and his research team who worked on a class of catalysts for producing unbranched polyethylene [11]. Phillips technique, developed by Phillips Petroleum researchers who were contemporary with Ziegler team, is another well-known process for producing high density polyethylene that was commercialized in 1961 [12]. This last process required higher pressure and temperature and the resulting product had a slightly higher density. In the late 50's, Standard Oil company developed its own process for producing high density polyethylene but, contrary to two other processes, has not been widely acknowledged by the market [6].

The primary step of nearly all polyethylene production lines consists of converting ethane to ethylene. Ethane is separated from natural gas after its exploitation and in a thermal cracking unit is converted to ethylene under heat and in the presence of catalyst. Of the numerous polymerization techniques, in one of the most commonly used ones, the ethylene gas is pumped under pressure and heat into a fluidized reactor bed containing metallic catalyst particles. Due to the high stability of ethylene, its polymerization is driven in the presence of appropriate catalysts. The polyethylene granules deposit at the bottom of the reactor and ethylene gas is recycled from the top of the reactor. Granulated polyethylene obtained in the form of granular powder is then ready to be sold as is or melted and extruded as required. It is worth mentioning that polyethylene is classified as a thermoplastic polymer, meaning that it can be melted and reshaped a lot of times.

Microstructure of high density polyethylene 1.3

High density polyethylene fills the first rank among other types of polyethylene by the volume of production per year. Its nearly linear molecular structure contains a low level of defects which allows for the organization of chains in a regular order. Indeed, like most other solid, linear macromolecules, polyethylene is neither amorphous nor fully crystalline, but semi-crystalline; it exhibits a glass transition as well as a melting point. Upon solidification/cooling from melt, a high fraction of linear chains start getting packed into regular structures at countless points to form base-centered orthorhombic crystals [13]. Two schematic illustrations of Figure 1 show the crystalline order observed in semi-crystalline polyethylene from two different angles. The fine crystalline parts of polyethylene, as illustrated in Figure 2, look like lamellar branches that have grown radially at many nuclei on which secondary lamellar branches (twigs) have developed. On the other hand, the noncrystalline/interlamellar part of polyethylene fills the spaces between the crystalline lamellae. Each group of crystalline lamellae that appear as radially nucleated branches from a nucleus plus the noncrystalline parts housed in-between forms a microstructural superstructure called spherulite (Figure 3A). In fact, as the polyethylene is cooled down, the spherulites start to appear and grow (Figure 3B) and when they reach each other, they form more or less flat boundaries. The concept of spherulite in semi-crystalline polymers is the analogue of grains in polycrystalline materials. The sphrulitic microstructure of polyethylene is well evident from the micrograph of Figure 3B. Therefore, should polyethylene be viewed as a heterogeneous solid, its major constituents are crystalline lamellae and noncrystalline/interlamellar region.

A B

Crystalline lamellae in semicrystalline polyethylene

Crystallinity of semicrystalline polyethylene depends mainly on the level of branching in molecules and cooling conditions. In general, chains with lesser degree of defects are potentially more appropriate for higher degree of crystallinity. Additionally, if cooling from melt is slow then the chains have enough time for recrystallization while in rapid quenching, the time required for most chains to arrange in crystalline order is taken away. The degree of crystallinity is directly proportional to the density of polyethylene because the packing of noncrystalline zones is looser and less efficient.

The prevailing lattice structure seen in semicrystalline polyethylene is orthorhombic. Bunn was the first who reported the crystalline structure of polyethylene and measured its lattice parameters [18]. Figure 1 illustrates the orthorhombic lattice structure of polyethylene crystallite along with its lattice parameters. The less stable monoclinic structure in mechanically stressed samples [19] as well as the hexagonal structure in samples subjected to high pressure and temperature [20] has also been observed. Numerous studies have been devoted to investigating the elastic properties of the crystalline phase in polyethylene [14,[21][22][23][24][25]. The literature on the crystalline phase in polyethylene is quite extensive but since the elasticity of the crystalline phase is not of main concern in this work, it is not further examined here.

Noncrystalline domain in semicrystalline polyethylene

The noncrystalline phase in semicrystalline polyethylene is composed of the loops, cilia and bridges, as illustrated in a cartoon in Figure 4. Segments of chains that exit a crystallite and enter the same crystallite constitute loops or foldings, segments of chains connecting two crystallites are called bridges or tie chains, and loose ends of chains projecting out of crystallites are called cilia or tails. As there is no predefined order within the noncrystalline phase, sometimes it is referred to as the amorphous phase in the literature. When thought of as a continuum, however, the noncrystalline domain is not homogeneous and is in turn composed of two constituents.

Interphase layer in semicrystalline polyethylene

From a theoretical standpoint, Flory argued that the continuity of polymeric chains implies the development of a new order intermediate between the crystalline lamellae and the central amorphous phase [26]. He added that the existing order within the crystalline lamellae cannot abruptly disappear at the interface of the crystalline/noncrystalline interface and the sharp transition from the ordered region with orthorhombic symmetry to a disordered region with isotropic symmetry is unlikely. He used the term "interphase" to call the layer intermediate between the crystallite and the central amorphous phase. Mandelkern et al. [27] presented another theoretical argument about the formation of the interphase layer based on the severe conformational differences between the ordered and disordered states.

The presence of the interphase layer as the third constituent in semicrystalline polyethylene has experimentally been demonstrated too. A wide range of experimental techniques, including broad-line proton NMR [28], high resolution 13 C NMR [29,30], electron microscopy [31], small angle neutron scattering [START_REF] Russell | Neutron and X-ray scattering studies on semicrystalline polymer blends[END_REF], dielectric relaxation [START_REF] Hahn | Dielectric relaxation of the crystal-amorphous interphase in poly(vinylidene fluoride) and its blends with poly(methyl methacrylate)[END_REF][START_REF] Hahn | Evidence for a crystal-amorphous interphase in PVDF and PVDF/PMMA blends[END_REF], and Raman spectroscopy [START_REF] Strobl | RAMAN SPECTROSCOPIC METHOD FOR DETERMINING THE CRYSTALLINITY OF POLYETHYLENE[END_REF][START_REF] Glotin | A Raman spectroscopic study of the morphological structure of the polyethylenes[END_REF] confirm the significant volumetric contribution of the interphase region into the conventional semicrystalline polymers. Mandelkern [START_REF] Mandelkern | RELATION BETWEEN STRUCTURE AND PROPERTIES OF CRYSTALLINE POLYMERS[END_REF] asserts that the difference between the crystallinities measured in density and heat of fusion experiments is directly related to the interfacial content. The theoretical and experimental studies cited above in concert with many other similar ones led to the acknowledgment of the presence of the partially ordered interphase layer as the third constituent in semicrystalline polyethylene.

Most of the relevant studies dedicated to the interphase properties are limited to measuring the thickness and volume fraction. Yoon and Flory estimated an average thickness of 6-8 Å for the interphase layer [START_REF] Yoon | Chain packing at polymer interfaces[END_REF]. In another theoretical study, Flory et al. suggest a thickness of 10-12 Å which they claim is in accord with experimental measurements [START_REF] Flory | The interphase in lamellar semicrystalline polymers[END_REF]. Mandelkern et al. [START_REF] Mandelkern | The interphase thickness of linear polyethylene[END_REF] made use of the technique of Raman internal and longitudinal acoustical modes to estimate the thickness of the interphase layer for a variety of polyethylene samples. In their simulation study, Kumar and Yoon [START_REF] Kumar | Lattice model for crystal-amorphous interphases in lamellar semicrystalline polymers: Effects of tight-fold energy and chain incidence density[END_REF] estimated a thickness of 10-30 Å for the interphase layer in polyethylene. In't Veld et al. employed the concept of Gibbs dividing surface in their molecular simulation study and reported a thickness of around 15 Å for the interphase layer at 350 K [START_REF] Veld | Temperature-dependent thermal and elastic properties of the interlamellar phase of semicrystalline polyethylene by molecular simulation[END_REF].

Although the existence of the interphase layer has been evidenced and techniques have been proposed in the literature for measuring its dimensions, there is a marked lack of knowledge about its mechanical properties. Experimental characterization of this layer has remained elusive due to the interference of signals from different phases [START_REF] Strobl | RAMAN SPECTROSCOPIC METHOD FOR DETERMINING THE CRYSTALLINITY OF POLYETHYLENE[END_REF][START_REF] Naylor | Raman spectroscopy employed for the determination of the intermediate phase in polyethylene[END_REF]. The attempts to isolate the noncrystalline phase in polyethylene by preferential degradation of the crystallite were unsuccessful due to the reconstruction of the material at the interface [START_REF] Olley | A permanganic etchant for polyolefines[END_REF]. The nanometric thickness of the interphase layer as well as its thermodynamical metastability act as essential impediments to any experimental characterization [START_REF] Veld | Temperature-dependent thermal and elastic properties of the interlamellar phase of semicrystalline polyethylene by molecular simulation[END_REF].

It should be emphasized that, although the role of the interfacial region is important as a connector between the two primary phases without which the polyethylene would be a very weak substance, there is no distinct boundary between the interphase layer and the central amorphous region and they cannot be demarcated by considering a sharp line separating them. Sometimes in the literature the designation "amorphous" is used to indicate the noncrystalline region housed between crystallites. This way of nomination can be misleading as the interlamellar domain is composed of partially ordered and fully disordered regions and the ensemble of the two regions is not necessarily isotropic. To avoid any confusion, in this dissertation, the term amorphous indicates the totally disordered, central region situated between partially ordered interphase layers and the ensemble of the two regions is called as "noncrystalline" or "interlamellar" domain.

Composite-inclusion model 1.4

In most publications where the effective response of high density polyethylene from a macroscopic point of view and using a multiscale homogenization approach has been investigated, a double-inclusion model has been proposed as the basic representative unit for the reconstruction of the real spherulitic microstructure [START_REF] Lin | Structure and plastic deformation of polyethylene[END_REF][START_REF] Tomita | Computational characterization and evaluation of deformation behavior of spherulite of high density polyethylene in mesoscale domain[END_REF][START_REF] Van Dommelen | Micromechanical modeling of intraspherulitic deformation of semicrystalline polymers[END_REF][START_REF] Van Dommelen | Micromechanical modeling of the elasto-viscoplastic behavior of semi-crystalline polymers[END_REF][START_REF] Nikolov | A micro/macro constitutive model for the small-deformation behavior of polyethylene[END_REF][START_REF] Lee | Simulation of large strain plastic deformation and texture evolution in high density polyethylene[END_REF][START_REF] Ahzi | Modeling of plastic deformation evolution of anisotropy in Semi-Crystalline Polymers[END_REF]. This basic structural element was first proposed by Ahzi et al. [START_REF] Ahzi | Modeling of plastic deformation evolution of anisotropy in Semi-Crystalline Polymers[END_REF] to investigate the large deformation analysis of semi-crystalline polyethylene. Indeed, this way of representation can be justified based on the spherulitic morphology of the undeformed semi-crystalline polymers which is observed when crystallized from melt. The radially oriented "fins" of each spherulite constitute the building blocks of this mesostructural unit in most multiscale studies of semi-crystalline polymers. The schematic of this unit element, which is often encountered in the literature, is given in Figure 5: a two-layer composite model made of bonding a crystalline layer to a noncrystalline one. As depicted in Figure 5, the orientation of each composite inclusion is denoted by n, the normal to the hypothetical interface plane separating the two layers. The crystallographic c-axis of the crystallite denotes the molecular chain alignment in the lamella. It has been demonstrated experimentally and theoretically that the misalignment between n and c is around 30° in the undeformed samples of polyethylene [START_REF] Mandelkern | The interphase thickness of linear polyethylene[END_REF][START_REF] Hütter | Polyethylene {201} crystal surface: interface stresses and thermodynamics[END_REF].

Uniformly random orientational distribution of a sufficiently large number of such composite inclusions in space leads to a configuration that is virtually equivalent to the macroscopically isotropic polyethylene sample. The higher the number of such building blocks is, the closer we get to the real morphology. In practice, however, the CPU time demand for such simulations grows dramatically as the number of grains exceeds several hundreds. This restriction comes from the solution of a system of couplednonlinear equations with a number of unknowns of (10) O for each inclusion, due to the constitutive law of each phase + their mechanical coupling requirements, and finally fulfilling the global compatibility and equilibrium conditions. The associated mathematical and numerical details are provided in Chapter 4.

Finally, it should be noted that the term "multiscale" analysis that is used for the numerical modeling and simulation of spherulitic microstructure in semi-crystalline polymers is due to the different size scales that are involved. Microscopic scale is reserved for the lamellar constituents due to their planar dimensions that reach a few micrometers. Any material point within the spherulites is treated on the mesoscopic scale, whereas an ensemble of multiple spherulites grouped together imparts to the macroscopic scale of the simulations. In a multiscale analysis, the macroscopic behavior of an aggregate in linked to the microscopic behavior of individual inclusions through a multiscale micromechanical homogenization approach. For want of a universal homogenization approach that captures all aspects of interactions between the inclusions, numerous approaches have so far been developed and are still being proposed [START_REF] Van Dommelen | Micromechanical modeling of the elasto-viscoplastic behavior of semi-crystalline polymers[END_REF][START_REF] Lee | Simulation of large strain plastic deformation and texture evolution in high density polyethylene[END_REF][START_REF] Ahzi | Modeling of plastic deformation evolution of anisotropy in Semi-Crystalline Polymers[END_REF][START_REF] Parks | Polycrystalline plastic deformation and texture evolution for crystals lacking five independent slip systems[END_REF][START_REF] Lee | Micromechanical modeling of large plastic deformation and texture evolution in semi-crystalline polymers[END_REF]. Their predictions may differ simply due to different underlying assumptions. The differences between some of the well-known available approaches are reviewed in Chapter IV where large deformation of an aggregate of polyethylene is re-investigated.

An expert is a person who has made all the mistakes that can be made in a very narrow field. (Niels Bohr-7 October1885 -18 November 1962)

Elastic characterization of the interphase layer in semi-crystalline polyethylene

INTERPHASE LAYER in semi-crystalline polyethylene is the least known constituent of this well-known polymer in terms of the mechanical properties. The nanometric size of this layer together with the lack of a criterion to precisely define its boundaries, and practical impossibility to prepare pure samples of interphase films in the laboratory, due to its intimate coupling to the coexisting adjacent phases, are major factors that impede its mechanical characterization, experimentally. The scale of the problem suggests resorting to atomistic/molecular simulation tools for investigating the mechanical properties of this constituent in question.

At the molecular simulation level, Hütter et al. [START_REF] Hütter | Polyethylene {201} crystal surface: interface stresses and thermodynamics[END_REF] employed the concept of a sharp Gibbs dividing surface in order to define a set of interfacial properties corresponding to the interphase. They obtained interfacial stresses and interfacial internal energies, but were not able to extract a value for the interfacial tension due to the significant contributions from its dependence on interfacial strain. For the interlamellar domain, however, the tool has proved effective in estimating the thermoelastic properties of the domain.

To date, the best attempts to take advantage of the molecular simulation tools with the aim of the mechanical characterization of the interlamellar domain in semicrystalline polyethylene have been made by Rutledge and co-workers who have developed a molecular level simulation capable of capturing the essential features of the interlamellar domain [START_REF] Veld | Temperature-dependent thermal and elastic properties of the interlamellar phase of semicrystalline polyethylene by molecular simulation[END_REF][START_REF] Balijepalli | Molecular simulation of the intercrystalline phase of chain molecules[END_REF][START_REF] Balijepalli | Conformational statistics of polymer chains in the interphase of semi-crystalline polymers[END_REF][START_REF] Gautam | Molecular simulations of the interlamellar phase in polymers: Effect of chain tilt[END_REF]. They employed the united atom model for simulating the interlamellar domain of polyethylene as a representative volume element surrounded by rigid crystalline lamellae [START_REF] Veld | Temperature-dependent thermal and elastic properties of the interlamellar phase of semicrystalline polyethylene by molecular simulation[END_REF]. Then, they utilized their model to compute the elastic properties of this domain over a range of temperatures.

Data provided by the molecular simulations 2.1

Part of the output of the Monte Carlo molecular simulations that is useful for our characterization analysis is the part that relates to the interlamellar stiffness components, il C , the interlamellar thickness, il t , and the interphase thickness, ip t , as a function of temperature for the range of 350-450 K. The diagrams of the normal stiffness components of il C vs. temperature are given in Figure 6. They suggest that all stiffness components do not follow the same trend since il 11 C does not show a monotonically increasing or decreasing behavior with temperature. The parameters of the quadratic curves fitted to the relevant atomistic data are given below. It is underlined that these dependences are only valid for 350-450K.

( ) ( ) ( ) ( ) (
)

il
C T T C T T C T T C T T C T T = - + - = - + = - + = - + = - + ( ) il 2 13
.1705 MPa 0.0521523 45.990307 10969.641

C T T = - + ( 1)
A B The other components of il C were not as important as the normal components, in addition to the fact that their calculation were computationally expensive. These components have been evaluated at the typical temperature 435 K and are taken to be constant across the temperature range of interest for the lack of any better estimate. The interlamellar stiffness evaluated at the typical temperature 435 K is given below

-     -     =   -     - -   -   ∼ C ( 2)
As suggested by the 6×6 matrix form of il C , the interlamellar domain has a monoclinic (or probably higher) symmetry. The likely physical grounds of the observed symmetry are discussed in Subsection 2.4. The uncertainty intervals reported for the interlamellar stiffness components are 30 MPa for the normal components, 100 MPa for il 44 C and 60 MPa for the two other shearing stiffnesses [START_REF] Veld | Temperature-dependent thermal and elastic properties of the interlamellar phase of semicrystalline polyethylene by molecular simulation[END_REF]. For the non-orthotropic terms, no uncertainty interval is reported.

In addition to the interlamellar stiffness, the interlamellar and interphase thicknesses have also been measured at different temperatures (Figure 7A,B). According to in 't Veld et al. [START_REF] Veld | Temperature-dependent thermal and elastic properties of the interlamellar phase of semicrystalline polyethylene by molecular simulation[END_REF], for calculating the interphase thickness, the concept of "sharp Gibbs dividing surface" has been employed. The interlamellar thickness, il t , varies almost linearly with temperature whereas the interphase thickness, ip t , can be well described using a quadratic estimate.

tip(nm) til(nm)

For the elastic characterization analysis, the volume fraction of the interphase and amorphous regions are also required. These quantities are easily calculable from the thickness of both phases as follows

il ip ip am ip il il 2 2 = , = t t t t t η η - ( 5)
Indeed, due to the periodic boundary conditions applied on the lateral planes of the interlamellar domain, the above relationships between the thicknesses and volume fractions are valid. But for implementing the dissociation analysis, the variation of the amorphous stiffness, am C , with temperature is also required. This last quantity is established based on the relevant findings from the literature, as explained in the following Subsection.

Temperature dependence of the central amorphous stiffness 2.2

Elastic properties of the central amorphous phase in polyethylene have been the topic of quite few experimental and theoretical studies (e.g. [START_REF] Crist | Mechanical properties of model polyethylenes: Tensile elastic modulus and yield stress[END_REF][START_REF] Hellwege | Die isotherme Kompressibilität einiger amorpher und teilkristalliner Hochpolymerer im Temperaturbereich von 20-250 °C und bei Drucken bis zu 2000 kp/cm2[END_REF][START_REF] Krigas | Model copolymers of ethylene with butene-1 made by hydrogenation of polybutadiene: chemical composition and selected physical properties[END_REF][START_REF] Janzen | Elastic moduli of semicrystalline polyethylenes compared with theoretical micromechanical models for composites[END_REF][START_REF] Humbert | A re-examination of the elastic modulus dependence on crystallinity in semi-crystalline polymers[END_REF]). As for the symmetry of this region, it has an isotropic symmetry [START_REF] Mandelkern | The interphase thickness of linear polyethylene[END_REF] simply due to its disordered structure. In practice, purely crystalline or purely amorphous polyethylene samples cannot be prepared in the laboratory [START_REF] Bédoui | Micromechanical modeling of isotropic elastic behavior of semicrystalline polymers[END_REF]. Therefore, the reported values for the elastic coefficients of the amorphous phase are based either on the theoretical arguments or on the modest extrapolation of corresponding coefficients of different polyethylene samples to zero crystallinity. Hellwege et al. [START_REF] Hellwege | Die isotherme Kompressibilität einiger amorpher und teilkristalliner Hochpolymerer im Temperaturbereich von 20-250 °C und bei Drucken bis zu 2000 kp/cm2[END_REF] measured the isothermal compressibility, which is the inverse of bulk modulus, for various polyethylene samples. Compared to the elastic modulus, bulk modulus shows a weaker sensitivity to the crystallinity. Therefore, extrapolation of this quantity provides an acceptable estimate for the central amorphous bulk modulus, am κ .

The extrapolated value reported for this quantity is am 1800MPa κ = (op. cit.). Krigas et al. [START_REF] Krigas | Model copolymers of ethylene with butene-1 made by hydrogenation of polybutadiene: chemical composition and selected physical properties[END_REF] conducted numerous tensile tests on different polyethylene specimens at several testing machine crosshead speeds less than 1 in/min. They concluded that except for low crystallinity samples, the measured Young's modulus remains independent of the strain rate and total strain for small strains. They also observed that the curve-fit of the measured Young's modulus and the one estimated from the plateau modulus intersect at the crystallinity of 0.03 ξ = . From the intercepts of the two diagrams, they concluded that am 3.5 0.5MPa

E = ±
, for the amorphous Young's modulus, where in conjunction with the bulk modulus reported in [START_REF] Hellwege | Die isotherme Kompressibilität einiger amorpher und teilkristalliner Hochpolymerer im Temperaturbereich von 20-250 °C und bei Drucken bis zu 2000 kp/cm2[END_REF], the amorphous Poisson's ratio will be am 0.4994 υ .

Using the same extrapolation technique, Crist et al. [START_REF] Crist | Mechanical properties of model polyethylenes: Tensile elastic modulus and yield stress[END_REF] From the kinetic theory of rubber elasticity [START_REF] Treloar | The Physics of Rubber Elasticity[END_REF], the plateau shear modulus of the liquid-like polyethylene is expressed as follows

0 N e = RT G M ρ ( 6
)
where ρ is the amorphous phase density, T the absolute temperature, R the ideal gas constant, and Due to its isotropic symmetry, the amorphous phase needs two independent elastic constants for the description of its stiffness tensor. The advantage of working with the amorphous Young's (or shear) modulus and Poisson's ratio lies in the understanding that we have from their temperature dependence behavior. As for am υ , the average of the above values is ~ 0.4996 which is already too close to the limiting value of 0.5. This closeness to 0.5 is due to the rubbery state of the amorphous phase at room temperature which leads to its further approaching this critical value at higher temperatures. It is therefore quite rational to assume that am υ remains nearly constant for the temperature range of interest with possibly negligible fluctuations around its mean value.

For am E , the average of the values reported in the previous studies for room temperature is ~ 5 MPa, which matches the mean value adopted by Humbert et al. [START_REF] Hellwege | Die isotherme Kompressibilität einiger amorpher und teilkristalliner Hochpolymerer im Temperaturbereich von 20-250 °C und bei Drucken bis zu 2000 kp/cm2[END_REF] for the amorphous phase of polyethylene. According to Eq.( 6) which has been introduced in the context of the kinetic theory of rubber elasticity [START_REF] Treloar | The Physics of Rubber Elasticity[END_REF], the elastic modulus of the amorphous phase is a linear function of temperature in the rubbery region. Moreover, using the first and second laws of thermodynamics and based on probabilistic discussions, it is demonstrated that the elastic modulus of a single chain in an amorphous polymer in the rubbery state is proportional to 2 3kT Nl ( 7) where k is the Boltzmann constant and N is the number of links in the chain each having an average length l [START_REF] Anil Kumar | Fundamentals Of Polymer Engineering (power Engineering[END_REF]. It is therefore quite justifiable to assume that am E is a linear function of temperature over the temperature range of interest where the amorphous phase of polyethylene is rubbery [START_REF] Nikolov | A micro/macro constitutive model for the small-deformation behavior of polyethylene[END_REF]. Accordingly, the following linear function is adopted for the temperature dependence of am

E for 350 K 400 K T ≤ ≤ ( ) am 5 MPa 293 T E = ( 8)
Using the above relationship and the constant value am 0.4996 υ = , the temperature dependent stiffness tensor of the central amorphous phase is built by substituting am E from Eq.( 8) and the temperature independent am 0.4996 υ = into the following relationship ( )

am am am am am 2 2 1 1 2 ijkl il jk ik jl ij kl E υ υ υ   = + +   + -   δ δ δ δ δ δ C ( 9
)
It should be noted that am E is a weak function of temperature (for instance ( )

am 400 K 6.8MPa E ∼
) and one may take am E to remain almost constant through the entire rubbery zone; this assumption would be consistent with the diagrams of the storage modulus vs. temperature for typical amorphous polymers in the rubbery regime.

Additionally, the two adopted amorphous elastic constants, namely am υ and am E , do not exhibit substantial changes within the temperature range considered here, to demand more precise estimations. Nevertheless, as demonstrated in Subsection 2.5, using an ad hoc sensitivity analysis the impact of the possible existing uncertainties in the adopted amorphous elastic constants on the generality of the conclusions is evaluated. Based on the careful examinations carried out, we believe that the basic conclusions remain essentially unaltered if more accurate forms of temperature dependence for am υ and am E were available. At any rate, the proposed methodology and dissociating tools presented herein remain applicable even if other forms of dependence are employed.

Micromechanical approaches of DIM and ECIM 2.3

Thus far, the temperature-dependent forms for the amorphous and interlamellar elastic tensors, namely am C and il C , respectively, together with the interphase and amorphous volume fractions are known and the only remaining unknown is the interphase stiffness. The schematic illustration of Figure 8 gives a brief graphical abstract of the problem. If there exists a relationship between all these quantities such as the following functional form ( )

il am ip am ip , , , , 0 F η η = C C C ( 10
)
then the interphase stiffness can be calculated or, as is alternatively referred to in this work, dissociated from the interlamellar stiffness. Such a relationship is provided by two distinct micromechanical approaches known as Double Inclusion Method (DIM) and Extended Composite Inclusion Model (ECIM).

Double Inclusion Method (DIM)

Developed by Hori and Nemat-Nasser [START_REF] Hori | Double-inclusion model and overall moduli of multi-phase composites[END_REF], DIM proposes an Eshelby-based formulation for evaluating the homogenized stiffness of an ellipsoidal inclusion encapsulating another ellipsoid with the entire double-inclusion being embedded in a reference medium. In the Monte Carlo molecular simulations, the periodic boundary conditions are imposed in a way that the interlamellar region can be treated as an inner, thin, disk-like ellipsoid, namely the core amorphous phase, wrapped by another hollow, thin, disk-like ellipsoid, namely the side interphase layers. Therefore, the problem under discussion fits the double-inclusion model if the Eshelby tensor of a disk-like ellipsoid (i.e. an ellipsoid with a very small aspect ratio) is used. After some mathematical manipulation, the DIM relationship for the interphase stiffness is rendered into the following form

( ) ( ) ( ) ( ) 1 1 1 1 1 1 1 1 ip ref ref ref il ref am am ip ip 1 η η η - - - - - - - - ∞ ∞ ∞ = - - - - - - +                             I C C S I C C S S C C C ( 11)
Here, I represents the fourth order identity tensor and ∞ S stands for the Eshelby tensor of a disk-like inclusion inserted in a reference medium whose stiffness is ref C . A concise summary of the governing equations of this method is brought in Appendix C. Further details about the assumptions and derivation procedure of DIM formulations are available in the work of Hori and Nemat-Nasser [START_REF] Hori | Double-inclusion model and overall moduli of multi-phase composites[END_REF].

Extended Composite Inclusion Model (ECIM)

Formulated first by Ahzi et al. [START_REF] Ahzi | Estimates of the overall elastic properties in semicrystalline polymers[END_REF][START_REF] Ahzi | Composite modeling for the effective elastic properties of semicrystalline polymers[END_REF] the composite inclusion model (CIM) is an attempt to find a compromise between the Voigt and Reuss mixture formulae for a layered composite inclusion by introducing strain and stress concentration tensors which serve also as weight functions. In this approach, the composite inclusion is made by stacking together two layers whose thicknesses are much less than the two other dimensions. The weight functions are established through the simultaneous enforcement of the continuity of deformation and equilibrium at the interface of the two composing layers. Here, the idea is extended to a three-layer composite inclusion, hence the designation "extended composite inclusion model" (ECIM). Again, owing to the nature of the periodic boundary conditions imposed in molecular simulations, the interlamellar region can be thought of as two thin interphase layers with one thin amorphous layer inserted in between. With reference to the notation defined in Appendix D, the dissociative analogue of Eq.( 10) solved for the unknown ip C is given by either of the following two equivalent equations:

( ) ( ) ( ) ( ) ( ) 1
η η η η - - - - = - - = - -         C C R C C Q Q C R C ( 12)
where ip S is a function of ref C as well as the aspect ratios (or geometry) of the ellipsoidal inclusion. Apart from the geometry of the problem, which has been taken into account during the derivation of the ECIM formulae, there is no need to resort to the concept of some "reference medium" in the ECIM formulation. It is, therefore, anticipated that for this specific case where the ellipsoidal inclusion is disk-like, the DIM results will be independent of the choice of C . Therefore, the equivalence between the independent variables of the two dissociative approaches is taken for granted.

Results and discussion

2.4

For reminding purposes, the general form of the interlamellar stiffness, which is adopted from the Monte Carlo molecular simulation study [START_REF] Veld | Temperature-dependent thermal and elastic properties of the interlamellar phase of semicrystalline polyethylene by molecular simulation[END_REF], is re-produced below.

C C C C C C C C C   -   -     =   ± -     - - ±   - ±     C ( 13
)
The entries of the upper left 3×3 submatrix are substituted from the relationships ( 1) keeping in mind that their uncertainty interval is 30 MPa ± . For the other nonzero components, the uncertainty interval of each component is given beside its respective mean value. Before explaining the details of the results, it is emphasized that since the melting temperature of polyethylene is around 407 K, only the temperature range 350-400 K is considered for the dissociation analysis.

The interphase stiffness, ip C , output by the two dissociation approaches, i.e. DIM and ECIM, for the typical temperature 370 K are given in Table 1. As explained earlier, DIM involves the Eshelby tensor for a disk-shaped inclusion, which must be evaluated numerically as it has no closed-form solution in the general case when the reference medium is anisotropic. Therefore, in the numerical code developed, a very small positive number, say 10 -6 , has been assigned to the aspect ratio of the disk-like ellipsoid for the computation of the Eshelby tensor. It was also observed that due to the recursive nature of the ECIM in the dissociation mode (see Appendix D), the method has shown numerical divergence despite deploying several stabilizing strategies. Therefore, the following numerical alternative was invoked. First, il C was symbolically calculated using an unknown ip C . As a result, a system of 13 coupled equations with 13 unknowns is thus obtained for the solution of which a hybrid optimization algorithm has been employed. The two-step, hybrid optimization algorithm consists of combining the Genetic Algorithm with another non-linear optimization technique called the Nelder-Mead (or simplex search) method [START_REF] Nelder | A Simplex Method for Function Minimization[END_REF]. In the first step, a ballpark estimate for the solution is found using the Genetic Algorithm which is used as the initial guess for the Nelder-Mead method in the second step. The dissociated ip C attributed to ECIM in Table 1 is the result of this combinatory numerical method.

A quick comparison reveals that the results of the two micromechanical methods agree perfectly. From a practical point of view, however, the ECIM formulation is fairly straightforward and simpler than the DIM formulation but is less efficient in the dissociation mode. Interestingly, and as discussed in Subsection 2.5, when the nonorthotropic elements of the interlamellar stiffness are neglected, the combination of the dual ECIM formulae ( 12) converges to the solution using the recursive method, which is much faster than the hybrid optimization technique. On the contrary, the DIM is very fast in both dissociation and homogenization modes although its major drawback is the development of the rather complicated numerical code for calculating the Eshelby tensor. In the direct/homogenization mode, however, when the dissociated ip C is used in combination with am C to produce the initial il C , the ECIM and the DIM produce the correct solution quite fast. In view of the details provided in Appendix D, the ECIM formulation in the homogenization mode takes the following explicit form

( ) ( ) ( )( ) ( ) il ip ip am am am am 1 1 1 il ip ip am am am am 1 or 1 η η η η - - - = - + = - + C C C C C Q C Q R R ( 14)
and the DIM in the direct mode of two-component homogenization takes the following form of mathematical representation

( ) ( ) ( ) ( ) 1 1 1 1 1 1 il ref am ref ref ip ref ref am ip η η - - - - - - ∞ ∞ ∞ = + + - + + - -                     I S C C C S C C C S C C ( 15)
Therefore, once ip C , am C and their volume fractions are known, the ECIM and the DIM directly return the solution, namely il C , in a single step without requiring any sophisticated numerical technique. Besides, it should be underlined that the interphase stiffness components shown in boldface in Table 1 indicate that ip C is not positive definite. This finding is discussed in depth in the following Subsections.

Finally, it can be seen that il C and ip C have the same form of symmetry. To justify, one may reason that since the central amorphous phase has the highest symmetry and the interlamellar region has the second lowest symmetry, it is quite anticipated that the interphase layer possesses the same form of symmetry as the interlamellar domain does. However, the orthorhombic crystalline structure of crystallite takes its known form of orthotropic representation in the material reference frame of the crystallite coinciding the abc-frame of the lattice. In the model of in 't Veld et al. [START_REF] Veld | Temperature-dependent thermal and elastic properties of the interlamellar phase of semicrystalline polyethylene by molecular simulation[END_REF] the ab-plane of the lattice makes an angle of ~30° with the crystallinenoncrystalline interface plane, which coincides the xy-plane [START_REF] Hütter | Monte Carlo Simulations of Semicrystalline Polyethylene: Interlamellar Domain and Crystal-Melt Interface[END_REF] (see Figure 5). When the crystalline stiffness with its familiar form of representation in the abc-frame is rotated into the xyz-frame, it takes a new representation that looks exactly like il C in terms of zeros and non-zero elements, as described previously by in 't Veld et al. [START_REF] Veld | Temperature-dependent thermal and elastic properties of the interlamellar phase of semicrystalline polyethylene by molecular simulation[END_REF]. For a better comprehension, this fact is illustrated in Figure 9. Therefore, it is very likely that the symmetry of the crystalline stiffness is responsible for the induced symmetry in ip C . However, because of the uncertainties available in the components of ip C , it is impossible to determine the number of independent components of the interphase stiffness and as a result impossible to determine the true symmetry of the interphase layer.

C C C C C C C C C C C C                     from to
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 C C C C C C C C C C C C C C C C C C C C ′ ′ ′ ′     ′ ′ ′ ′   ′ ′ ′ ′     ′ ′     ′ ′ ′ ′   ′ ′    

A probe into the shearing components of il C

Due to their critical role in the dissociation analysis, the shearing stiffnesses in il C are examined more carefully. C for several temperatures are plotted in Figure 10, Figure 11 and Figure 12, respectively. As is qualitatively evident from the diagrams and as demarcated by vertical dashed lines, there are intervals of il 44 C for which the dependent variables become unbounded, which is unacceptable. Therefore, these intervals must be excluded from

il 44 100MPa 100MPa C - ≤ ≤
. Specifically, Figure 10, Figure 11 and Figure 12 suggest that the imprecise intervals ( ) C takes incomparably large value. C , respectively. As demonstrated in section 2.5, this property is insensitive to the uncertainties of the adopted amorphous elastic constants and the uncertainties of il C .

Deviation of the interphase stiffness from positive definiteness

Positive definiteness of the stiffness tensor for stable materials found in nature is demonstrated based on the first law of thermodynamics and the positivity of the elastic strain energy. If the stiffness tensor is represented in a 6 6

× matrix form, positive definiteness requires positivity of the diagonal elements. Here, this requirement is violated at least for the interphase stiffness since ip 44 C and ip 55 C , although too close to zero, robustly take negative values at least in the temperature range 350-400 K. It is worth noting that unlike either the crystalline or amorphous phases, the interphase and interlamellar domains are not necessarily thermodynamically stable phases that can ever exist in the absence of the stabilizing influence of the adjoining crystalline lamellae. Thus there is no compelling reason to require their mechanical stability in isolation either. Here, in our example, the negative shear stiffnesses are only observed in the transversal plane of the interphase layer, whose thickness is ~ 1 nm and plays the role of the transition region between the crystallites and the amorphous phase. Moreover, negativity of the shear modulus has been observed earlier for nanoscale domains within an amorphous matrix [START_REF] Yoshimoto | Mechanical heterogeneities in model polymer glasses at small length scales[END_REF]. Other examples of the studies available in the open literature on the heterogeneous materials containing at least one component with non-positive definite stiffness include [START_REF] Lakes | Dramatically stiffer elastic composite materials due to a negative stiffness phase[END_REF][START_REF] Kochmann | Dynamic stability analysis of an elastic composite material having a negative-stiffness phase[END_REF][START_REF] Huang | Theoretical investigation of the behavior of an acoustic metamaterial with extreme Young's modulus[END_REF][START_REF] Lakes | Extreme damping in composite materials with negative-stiffness inclusions[END_REF][START_REF] Jaglinski | Negative Stiffness and Negative Poisson's Ratio in Materials which Undergo a Phase Transformation[END_REF]. Therefore, while one or more constituents in a mechanically stable, heterogeneous solid might have non-positive definite stiffness, the global stiffness of the entire heterogeneous material is positive definite.

For the crystalline/noncrystalline interface in semi-crystalline polymers, Cammarate and Eby [START_REF] Cammarata | Effects and measurement of internal surface stresses in materials with ultrafine microstructures[END_REF], Fisher et al. [START_REF] Fisher | Surface stresses in paraffin and polyethylene[END_REF] and Cammarata and Sieradzki [START_REF] Cammarata | Surface and interface stresses[END_REF] report negative interface stresses of the order of -0.1 to -0.3 J/m 2 for a { } 001 surface. In reference to these negative stress values, Hütter et al. [START_REF] Hütter | Polyethylene {201} crystal surface: interface stresses and thermodynamics[END_REF] argue that the fact that the interface is under pressure can be attributed to the way the polymer chains exit (or enter) the crystal. Therefore, the interphase layer is physically pre-strained. On the other hand, pre-strained objects, which contain stored energy, are well-known examples in which non-positive definite stiffness may take place [START_REF] Kochmann | Analytical stability conditions for elastic composite materials with a non-positive-definite phase[END_REF][START_REF] Kochmann | Stability criteria for continuous and discrete elastic composites and the influence of geometry on the stability of a negative-stiffness phase[END_REF]. As a result, non-positive definite stiffness of the interphase layer is in accord with one's anticipation in light of the earlier studies. Non-positive definite stiffness of the interphase layer means that under particular strain fields the associated strain energy of one or more interphase layer(s) becomes negative (equivalent to the release of some stored energy). However, due to the stabilizing effects of the amorphous and crystalline layers whose stiffness tensors are positive definite, the global strain energy of any polyethylene sample that exists in isolation becomes positive under any kinematically admissible strain field. In other words, since polyethylene samples usually have a uniformly random distribution of constituents, they are macroscopically isotropic satisfying the sufficient conditions of stability. According to Kochmann [START_REF] Kochmann | Stability criteria for continuous and discrete elastic composites and the influence of geometry on the stability of a negative-stiffness phase[END_REF] if the shear and bulk moduli of an isotropic solid are positive the sufficient conditions of stability are met and its strain energy is always positive under any nonzero strain field.

It is also worth noting that although ip C without displaying any asymptotic behavior. It should be reminded that negativity of the shear stiffnesses imply that upon imposition of positive corresponding shear strains, negative stresses will be produced. Then, one may reason that the negligible negativity of ip 44 C or ip 55 C produces negligible negative shear stresses, which may be tolerated by the surrounding media.

However, the negativity of ip 66 C is comparatively so large that, even with relatively small positive shear strains, it produces such large negative shear stresses that are not deemed to be balanced by the surrounding media. Additionally, when il

Verification of dissociation analysis results

Crist et al. [START_REF] Crist | Mechanical properties of model polyethylenes: Tensile elastic modulus and yield stress[END_REF] reported their measurements for different samples of PE, linear HPB, 3S HPB and 4S HPB spanning a crystallinity range of 0.35 0.68 ξ < < , which they were able to describe with a linear fit on a semilogarithmic scale with the correlation coefficient of 0.974. Assuming essentially only a two-component model comprising crystalline and noncrystalline (which they call amorphous) phases, they argue that in semi-crystalline PE, the average Young's modulus of the noncrystalline phase is nearly constant (~ 300 MPa) for 0.7 ξ > , but drops appreciably as the crystallinity is lowered below 70%. Therein, they proposed the following double-argument dependence for the average Young's modulus of the noncrystalline phase at room temperature

( ) ( ) am 2 exp 7.158 0.7 = MPa 300 0.7 E ξ ξ ξ ξ  ≤  ≥  ( 16
)
Given that Crist and co-workers did not consider the presence of a third component, i.e. the transitional interphase separating crystallites and the central amorphous phase, it is

reasonable to equate am E in their two-component model with il E . Indeed, il E is calculable form il C
which is, according to Eq.( 10), a function of the properties and volume fractions of its constituents. Therefore, mathematically speaking ( )

ip am il ip = , , E f η C C ( 17)
which is in agreement with the suggestion made by Crist et al. [START_REF] Crist | Mechanical properties of model polyethylenes: Tensile elastic modulus and yield stress[END_REF] since ip η , which denotes the interphase volume fraction within the interlamellar domain, is in its turn a function of the crystallinity. In other words, for an interlamellar region of arbitrary thickness one can write

ip ip ip ip am ip am 2 2 2 + 2 + V t V V t t η = ≃ ( 18
)
where ip t / ip V denotes the thickness/total volume of each interphase layer in a layered, "sandwich" model of the interlamellar domain. There is evidence that the interphase thickness is invariant with crystallinity, whereas the amorphous thickness varies to accommodate changes in interlamellar separation [START_REF] Yoon | Chain packing at polymer interfaces[END_REF]. Therefore am t and consequently ip η are functions of crystallinity:

( )

ip =h η ξ ( 19)
For the hypothetical state of no crystallinity, polyethylene is composed of the pure amorphous phase and there is no interphase which means:

ip am ip PE il am 0 1 1 if 0 then = = E E E η η η ξ → ⇒ = -→ →    ( 20
)
where at room temperature am E takes values between 2-11.4 MPa, as elaborated in Subsection 2.2. On the other hand, at high crystallinities, the amorphous phase disappears and the interlamellar domain will be dominated by the interphase layers, meaning that

ip am ip il ip 1 1 0 if 1 then = E E η η η ξ → ⇒ = -→ →      ( 21)
In summary, by increasing the crystallinity from zero, the interphase layers start to appear and the amorphous phase shrinks, implying that ip η is a positive and monotonically increasing function of crystallinity up to some critical crystallinity, cr ξ .

Since ip η cannot exceed unity in the interval of 0 1 ξ < < and eventually it has to go to unity as 1 ξ → , the most likely dependence form of ip η , which is in accord with the observation by Crist et al. [START_REF] Crist | Mechanical properties of model polyethylenes: Tensile elastic modulus and yield stress[END_REF], is the following: ) ( )

ip cr ip cr > η ξ ξ η ξ ≃ .
Consequently and in light of Eq.( 17), the average Young's modulus of the interlamellar domain, il E , becomes a function of crystallinity similar to that suggested by Crist et al. [START_REF] Crist | Mechanical properties of model polyethylenes: Tensile elastic modulus and yield stress[END_REF], while the constitutive properties of the constituents, namely ip C and am C , remain essentially independent of crystallinity. In other words, using a two-phase sandwich model to represent the interlamellar region in which the constitutive properties of the phases are independent of crystallinity, the form of dependence in Eq. ( 16) proposed by Crist et al. [START_REF] Crist | Mechanical properties of model polyethylenes: Tensile elastic modulus and yield stress[END_REF] can be justified.

Finally, to use the empirical Eq.( 16) for verification purposes, one can reason that according to Eq.( 21), at high crystallinity il ip E E ≃ , and therefore the average Young's modulus of the interphase layer at room temperature must be comparable to that of the interlamellar domain at high crystallinity, which is offered by the empirical relationship. Following the idea presented by Counts et al. [START_REF] Counts | A comparison of polycrystalline elastic properties computed by analytic homogenization schemes and FEM[END_REF], Hill's estimate [START_REF] Hill | The elastic behaviour of a crystalline aggregate[END_REF] is used to find an estimate for ip E from ip C . In the work of Counts et al. [START_REF] Counts | A comparison of polycrystalline elastic properties computed by analytic homogenization schemes and FEM[END_REF] the problem of estimating the overall shear and Young's moduli of a polycrystalline BCC Mg-Li, which takes nonpositive definite stiffnesses for some compositions of Mg-Li, is treated similarly. In fact, the average elastic moduli from FE and self-consistent analyses that are extrapolated to unstable regions agree very well with Hill's estimate for any composition. Therefore, for estimating the interphase average bulk and shear moduli using Hill's method, the bulk and shear moduli of the Voigt and Reuss approaches, calculable from ip C and

1 1 ip - - C
, where indicates the orientational (volume) averaging, are required. Thus, Hill's estimates of the interphase bulk and shear moduli are obtained as follows:

( )

( ) ip-Hill ip-V ip-R ip-Hill ip-V ip-R 1 = + 2 1 = + 2 G G G κ κ κ ( 22)
from which the average Young's modulus of the interphase layer obtained from Hill's estimate reads

ip-Hill ip-Hill ip-Hill ip-Hill ip-Hill 9 = 3 + G E G κ κ ( 23)
In the temperature range studied here, the closest to the room temperature is 350 K at which ip-Hill E is calculated to be 347 MPa, which compares well to the plateau value of 300 MPa proposed by Crist et al. [START_REF] Crist | Mechanical properties of model polyethylenes: Tensile elastic modulus and yield stress[END_REF], noting that the Young's modulus of amorphous polymers increases with temperature if the polymer is in the rubbery state. Furthermore, Ding et al. [START_REF] Ding | A molecular dynamics study of Young's modulus change of semi-crystalline polymers during degradation by chain scissions[END_REF] conducted a molecular simulation study on the Young's modulus change in a semi-crystalline polymer and observed that the Young's modulus of the interlamellar region increases with temperature in the rubbery state. It is therefore expected that the analogue of the empirical relation ( 16)

This result is in accord with our expectation that values higher than 300 MPa fall within the uncertainty interval of ip-Hill

E

. It is reiterated that for calculating the mean value of . For the entire temperature range, the evolution of ip-Hill

E

and the bounds of its uncertainty interval are similarly calculated. The plots of these evolutions are given in Figure 14. Interestingly, the three diagrams show a monotonically increasing trend with temperature in agreement with our anticipation. explained in most introductory textbooks on the behavior of polymers (see e.g. [START_REF] Anil Kumar | Fundamentals Of Polymer Engineering (power Engineering[END_REF] ), this observation is in accord with the common understanding that the elastic modulus of amorphous polymers in the rubbery region increases with temperature as the Young's modulus is controlled by the entropy change.

Sensitivity/uncertainty analysis 2.5

To examine the sensitivity of the non-positive definiteness of the calculated ip

Then in a total of 10 9 Monte Carlo cycles, at random temperatures, random am E and am υ are sampled from their assigned uncertainty intervals. Accordingly, a random am C is picked from its uncertainty space. Afterwards, at the same random temperature, a random il C is picked from its uncertainty space based on the uncertainty intervals obtained from the MC molecular simulations, except for il

44

C which is picked from ( ] 82,100 . Finally, using the DIM dissociation formulation, the new ip C and its eigenvalues are calculated. By carrying out this sensitivity analysis, none of the calculated interphase stiffnesses fulfilled the condition of positive definiteness. Given the extremely large number of the Monte Carlo cycles, it is very unlikely that one can find some temperature from 350K 400 K T < < and some random am C and il C , as described above, such that a positive definite ip C is obtained. Therefore, one can conclude, with a high degree of certainty, that the non-positive definiteness of

ip C within 350 K 400 K T < <
is an established fact and insensitive to the uncertainties of am C and il C . It should be noted that for the volume fractions involved in the calculations, no uncertainties were reported. However, if there were any, they would be treated similar to other uncertainties in the same way as described in the above.

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 320.5 0 0 0 0 0 109.3
As a last sensitivity check, only the orthotropic part of il C is preserved and the dissociation analysis at the same temperature of 370 K is carried out.

Hyperelastic characterization of the interlamellar domain and interphase layer in polyethylene

HYPERELASTIC characterization of the noncrystalline phase in polyethylene is the main objective of this chapter. Inspired by the characterization idea presented in the previous chapter, a similar methodology for the hyperelastic characterization of the interlamellar domain and interphase layer in high density polyethylene is presented here. Hyperelastic behavior of the interlamellar domain and its components, namely the interphase and central amorphous phase, is attributed to the rubbery state of this region

( PE g 173K T >
). This understanding is in accordance with what is adopted in the literature for the large deformation behavior of amorphous polymers in the rubbery regime [START_REF] Van Dommelen | Micromechanical modeling of the elasto-viscoplastic behavior of semi-crystalline polymers[END_REF][START_REF] Sedighiamiri | Micromechanical modeling of the deformation kinetics of semicrystalline polymers[END_REF][START_REF] Sedighiamiri | Micromechanics of semicrystalline polymers: Yield kinetics and long-term failure[END_REF]. Indeed, the need for the hyperelastic constitutive properties of the interlamellar domain and its constituents is better appreciated in numerical simulation of large deformation of semi-crystalline polyethylene [START_REF] Van Dommelen | Micromechanical modeling of the elasto-viscoplastic behavior of semi-crystalline polymers[END_REF][START_REF] Parks | Polycrystalline plastic deformation and texture evolution for crystals lacking five independent slip systems[END_REF][START_REF] Lee | Micromechanical modeling of large plastic deformation and texture evolution in semi-crystalline polymers[END_REF].

The basics of the methodology consist in applying the governing equations of the hyperelastic continuum to stress-strain data from molecular dynamics simulations of a polyethylene stack [START_REF] Lee | Plastic deformation of semicrystalline polyethylene by molecular simulation[END_REF]. Due to its simplicity, efficiency and widespread use, the neo-Hookean model of rubber elasticity is adopted for the interlamellar domain and its constituents [START_REF] Van Dommelen | Micromechanical modeling of the elasto-viscoplastic behavior of semi-crystalline polymers[END_REF][START_REF] Sedighiamiri | Micromechanical modeling of the deformation kinetics of semicrystalline polymers[END_REF][START_REF] Sedighiamiri | Micromechanics of semicrystalline polymers: Yield kinetics and long-term failure[END_REF]. It should be emphasized that although the presented methodology is founded on the application of the neo-Hookean model, the identification notion does not lose its generality and can be deployed in a similar way with any other hyperelastic model such as Yeoh model, Moony-Rivlin model, generalized Gent model, Arruda-Boyce model and so forth. [START_REF] Steinmann | Hyperelastic models for rubber-like materials: Consistent tangent operators and suitability for Treloar's data[END_REF][START_REF] Hossain | More hyperelastic models for rubber-like materials: consistent tangent operators and comparative study[END_REF].

Simulation stack of the molecular dynamics analysis and the

imposed boundary conditions

The MD simulation results that are employed here relate to a stack of PE that is comprised of two crystalline lamellae plus the interlamellar domain housed between them (see Figure 19A). The schematic of the simulated PE stack in its unloaded state along with the diagram of the density distribution in the longitudinal direction is reproduced from the basic MD study [START_REF] Lee | Plastic deformation of semicrystalline polyethylene by molecular simulation[END_REF] in Figure 19B. It is worth noting that although the existence of the interphase layer that separates the core amorphous zone from the crystallites has extensively been evidenced, the transition in the cup-shaped density diagram again confirms the presence of this third, transitional constituent. Further details vis-à-vis the atomistic simulation are provided in the MD study by Lee and Rutledge [START_REF] Lee | Plastic deformation of semicrystalline polyethylene by molecular simulation[END_REF]. deployed in this work. The main objective of this identification practice is the hyperelastic characterization of the amorphous, interphase and interlamellar regions. The adopted approach is the minimization of a set of nonnegative objective functions that are defined using the assumption that the three domains of interest behave as isotropic, compressible, hyperelastic continua.

Methodology

3.2

The polyethylene element for which the MD simulations have been conducted by Lee and Rutledge [START_REF] Lee | Plastic deformation of semicrystalline polyethylene by molecular simulation[END_REF] is subject to a global deformation gradient tensor, stack F . The incrementally imposed constraints of the stack during the simulation consists of longitudinal stretch of the stack along the z-axis, 1 z z λ ε = + , and equal lateral stretches, x λ and y λ . Whereas the interlamellar domain is assumed to be compressible, the lateral stretches are adjusted such that the entire stack deforms isochorically, meaning that ( )

1 stack 2 0 0 det det 0 0 1 0 0 x y x y z z λ λ λ λ λ λ -     = = ⇒ = =       F ( 26
)
Note that according to the simulation constraints, the lateral stretches are identical everywhere along the stack whereas the longitudinal stretch z λ represents the average stretch of the entire stack. In other words, at each increment

1 2 x y z λ λ λ - = =
for all the constituents of the stack while their respective longitudinal stretches may differ, i.e.

il am ip z z z z λ λ λ λ - - - ≠ ≠ ≠ .
The superscripts or subscripts "il", "am" and "ip", used here and later on, indicate the interlamellar, amorphous and interphase regions, respectively. Accordingly, the deformation gradient tensor of the interlamellar, amorphous and interphase domains takes the following generic form

1 2 1 β 2 -β 0 0 0 0 0 0 z z z λ λ λ - -       =         F ( 27
)
where β replaces "il", "am" and "ip". It is worth mentioning that, contrary to z λ , the evolution of -β z λ is a priori unknown and has not been recorded during the simulation.

Indeed,

il z λ -,

Hyperelastic constitutive law 3.3

The constitutive law of a compressible, isotropic, hyperelastic solid following the neo-Hookean model reads ( )

5 3 1 J J µ κ ′ = + - σ B 1 ( 28
)
where σ is the Cauchy stress, J represents the determinant of the deformation gradient tensor F, ′ B stands for the deviatoric part of the left Cauchy-Green deformation tensor

T =
B FF , 1 is the second-order identity tensor, and material parameters µ and κ denote shear and bulk moduli, respectively. On the other hand, however, only the interlamellar region and its constituents behave hyperelastically and not the crystallites. Therefore, the deformation gradient tensor of the interlamellar domain and its constituents, β F where β replaces "il" or "am" or "ip", must be substituted in Eq.( 28) for the identification of their hyperelastic parameters. In the following, the procedure of hyperelastic characterization for the central amorphous phase and the interlamellar domain is elaborated first. Then the governing relationships of the interphase layer are incorporated into the final objective function where the hyperelastic parameters of the three domains are identified concurrently.

Substitution of am/il

F into Eq.( 28) yields ( ) ( )

5
µ λ µ λ λ λ λ λ λ κ λ σ µ λ σ µ λ λ λ λ λ λ κ λ σ - - - -   - + + -          = - + + -       ( ) 5 5 5 1 2 1 3 3 3 3 3 3 am/il -am/il am/il -am/il -am/il -am/il am/il 2 1 3 z z z z z z z z µ λ µ λ λ λ λ λ λ κ λ -                        - + + -         ( 29)
Obviously, under the current simulation conditions, the shearing strains remain zero throughout the loading due to the simulation constraints. Additionally, it is understood from Eq.( 29) that the lateral stresses with the resolution of 1 Å during each increment. The two distributions at twelve sample strain increments are shown in Figure 20. Comparison of the two distributions for all 800 increments, not shown here, reveals that the two distributions show noticeable disagreements for the most part of the loading. For the crystalline domain, the disagreement would have quite been expected because of the anisotropy of the crystallites. Although the disagreement between the two lateral stresses for the amorphous phase is not as strong as it is for the crystalline phase, it does not agree with our expectations given its initially isotropic symmetry. In the latter case, the source of amorphous phase or interlamellar region, the mean value of am lat σ and il lat σ , respectively, averaged over their respective length is intended. In other words, the evolution of the amorphous/interlamellar length (or boundaries) in the z-direction, which is initially unknown, is required for calculating the evolution of the mean value of am/il lat σ . Strictly speaking, at any increment, the current length of each domain in the z-direction is required for calculating the average of lateral stress of the domains involved.

Given the equality of the right hand sides of the first two equations of the system of equations ( 29), one of the relationships is obviously trivial. Because of the redundancy of one of these relationships, the governing stress-stretch relationships for the amorphous/interlamellar region under the given mode of deformation are reduced to the following system of two independent equations ( )

( ) 5
z z z z z z z zz z z z z z µ λ σ λ λ λ λ κ λ µ λ σ λ λ λ λ κ λ - -    = - + -           = - + -       ( 30

Identification approach 3.4

Attempts were made to take advantage of the distribution of the "order parameter", P2, for extracting the evolution of the amorphous/interlamellar boundaries. Calculated by means of the second Legendre polynomial coefficient, P2 is often used to quantify the orientational ordering in nematic systems [START_REF] Lee | Plastic deformation of semicrystalline polyethylene by molecular simulation[END_REF]. P2 takes a maximum value of 1 for a perfectly uniaxially aligned system, 0 for a randomly oriented system, and -0.5 for a system with random orientational order in a plane perpendicular to the axis of symmetry. This parameter is also shown in Figure 20 at twelve sample increments. The lowest part of the P2 distribution belongs to the central amorphous region, the highest values of the distribution signify the crystalline lamellae on either side of the simulation domain, and the transitional parts are characteristic of the interphase layers. Nonetheless, low resolution of data recording (4 Å) used for P2, lack of a universally established benchmark founded on P2 for separating different regions, and finally melting and recrystallization phenomena that occur at medium strains impede any attempt to rely on this parameter for tracing the evolution of the initial amorphous/interlamellar elements.

After trying several identification algorithms (e.g. taking either xx σ or yy σ or their arithmetic mean as the true lateral stress for each domain under study with boundary evolutions roughly approximated from P2 distributions and lateral stress distributions, …), the best scenario that suits the problem at hand is the one in which the hyperelastic coefficients of all domains involved along with their respective boundaries are identified simultaneously. The details of the implemented identification algorithm are as follows.

The evolutions of the interlamellar/amorphous right and left boundaries, i.e. il-r z , am-r z , il-l z and am-l z , are evidently functions of the longitudinal strain, z ε . To clarify further, it is worth noting that, as discussed in the foregoing Subsection and in reference to the MD simulation conditions, the loading constraints are of displacement type, meaning that the longitudinal strain, z ε , is imposed incrementally and simultaneously with lateral stretches such that:

1 1 1 x y z ε ε ε = = + -. With
zz 1 z f σ ε =
, the second relationship of the system of equations ( 30) suggests that ( )

am/il 2 z z f λ ε - =
where in conjunction with the relationships of Eq.( 32), it is deduced that for this particular problem and simulation conditions, the interlamellar/amorphous right and left boundaries are functions of z ε . It should be noted that the right and left boundaries of the amorphous and interlamellar domains are not necessarily symmetric with respect to the origin, z=0. Therefore, these boundaries are specified with "r" and "l", designating right and left sides, respectively. Various functional forms, including polynomial, Fourier expansion, Gaussian expansion, etc., with unknown parameters can be assigned to each boundary. Here, a three-argument, polynomial form of dependence is adopted for each boundary. For instance, the right boundary of the interlamellar domain is described using the following form

z z z z z z z z z z z z z α α ε α ε α ε ε α α ε α ε α ε ε α α ε α ε α ε ε  + + + ≤ ≤  = + + + < ≤   + + + < ≤  ( 31)
where i α s are unknown parameters that remain to be identified. Forms of dependence similar to relation ( 31) are adopted for am-r z , il-l z and am-l z too, leading to a total of 48 design parameters, i α . The above form of dependence is easily deduced from the second relationship of Eq. ( 30) σ , therefore the amorphous and interlamellar boundaries will naturally exhibit similar dependence. This is the reason why a three-argument functional form for these boundaries is a more elegant option in lieu of a smooth, single argument one. Having adopted the functional forms of il-r z , am-r z , il-l z and am-l z with 48 unknown parameters, the stretches of the amorphous, interlamellar, and right and left interphase domains read as follows

am-r am-l il-r am-r am-l il-l il-r il-l z-am z-il z-ip-r z-ip-l am-r am-l il-r il-l il-r am-r am-l il-l 0 0 0 0 0 0 0 0 , , , z z z z z z z z z z z z z z z z λ λ λ λ - - - - = = = = - - - - ( 32)
where 0 z indicates the position at the unloaded state. With the above stretches, the hyperelastic parameters of the amorphous and interlamellar domains can be calculated using the second relationship of Eq.( 30) and a simple least squares technique. Expansion of the second relationship of Eq.( 30) for the amorphous and interlamellar domains for the 800 increments of loading yields -am -am -am 800 800

2 1 3 2 1 3 2 1 3 zz z z z z z z z z z z z z z z z z z z z σ σ µ κ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ - - - Λ     - -               - -       =                 - -           ⋮ ⋮ 2 800 z zz Y σ               ⋮ ( 33) ( ) ( ) ( ) 5 5 1 2 3 3 3 3 1 5 5 1 2 3 3 3 3 il il 5 5 1 2 3 3 3 3 il -il -il -il 1 1 -il -il -il 2 2 -il -il -il 800 800 2 1 3 2 1 3 2 1 3 zz z z z z z z z z z z z z z z z z z z z σ σ µ κ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ - - - Λ     - -               - -       =                 - -           ⋮ ⋮ 2 800 z zz Y σ               ⋮ ( 34)
Obviously, the above systems of equations do not have solutions for the amorphous and interlamellar hyperelastic parameters since the number of independent equations is way greater than the number of unknowns. However, to find the best estimate for these coefficients, the least squares technique is used. According to the least squares technique [START_REF] Nelles | Nonlinear System Identification: From Classical Approaches to Neural Networks and Fuzzy Models[END_REF], the best set of hyperelastic parameters for the amorphous and interlamellar domains are obtained as follows

Y Err Y Y Err Y µ µ κ κ µ µ κ κ - -     = Λ Λ Λ = -Λ             = Λ Λ Λ = -Λ         ( 35
)
where am/il Λ and Y are introduced in relationships ( 33) and [START_REF] Hahn | Evidence for a crystal-amorphous interphase in PVDF and PVDF/PMMA blends[END_REF] with am Err and il Err being the errors of the least squares estimates. In the ideal case, these error expressions go to zero. In general, however, they must be brought as close to zero as possible. Therefore, a multi-objective optimization problem is formulated in which the first and second objective functions are defined as the least squares errors of the above estimates:

( ) ( ) 1 1 T T T T 1 am am am am 2 il il il il , OF Y Y OF Y Y - - = -Λ Λ Λ Λ = -Λ Λ Λ Λ ( 36)
Moreover, the system of equations ( 30) can be similarly developed for both interphase layers.

( ) ( ) ( )

5 5 5 5 2 1 2 1 3 3 3 3 3 3 3 3 5 5 1 2 3 3 3 3 ip -ip-r ip -ip-l ip-r ip-l lat -ip-r -ip-r ip lat -ip-l -ip-l ip ip -ip-r ip -ip-r -ip-r ip 1 1 3 3 , 2 2 1 3 3 z z z z z z z z z z z z z zz z z z z zz z µ λ µ λ σ λ λ λ λ κ σ λ λ λ λ κ λ λ µ λ µ σ λ λ λ λ κ σ λ λ - - -      = - + - = - + -               = - + - =       ( ) 5 5 1 2 3 3 3 3 -ip-l -ip-l -ip-l ip 1 z z z z z z λ λ λ λ κ λ -        - + -       ( 37)
Expecting identical constitutive properties for both interphase layers, the second relationships of the systems of equations ( 37

In addition, the hyperelastic parameters of the interphase layer can be identified like the interlamellar and amorphous hyperelastic parameters:

( ) ( ) ( ) 5 5 1 2 3 3 3 3 5 5 1 2 3 3 3 3 5 5 1 2 3 3 3 3 ip ip ip -ip -ip -ip 1 1 -ip -ip -ip 2 2 -ip -ip -ip 800 800 2 1 3 2 1 3 2 1 3 z z z z z z z z z z z z z z z z z z z µ κ λ λ λ λ λ λ σ λ λ λ λ λ λ λ λ λ λ λ λ - - - Λ     - -               - -         =                   - -           ⋮ ⋮ ( ) ip-ls ip-ls 1 T T ip ip ip 1 2 LS 800 z zz zz Y Y µ κ σ σ - Λ Λ Λ           ⇒ =             ⋮ ( 39
)
Similar to the first and second objective functions, the fourth objective function is defined as the error of the least squares estimate of Eq.( 39).

( )

1 T T 4 ip ip ip ip OF Y Y - = -Λ Λ Λ Λ ( 40
)
On the other hand, the force equilibrium in the lateral direction requires that

( ) ( ) ( ) ( ) ip-r am ip-l il il-r am-r lat am-r am-l lat am-l il-l lat il-r il-l lat z z z z z z z z σ σ σ σ - + - + - = - ( 41)
The average lateral stresses appearing in the above equation are calculable by substituting the corresponding hyperelastic parameters into the first relationship of Eq. ( 30). Therefore, the fifth objective function is constructed as follows

( ) ( ) ( ) ( ) ( ) ( ) 5 5 2 1 3 3 3 3 5 5 2 1 3 3 3 3 5 5 2 1 3 3 3 3 800 ip z-ip 5 am-l am-r il-r il-l z-ip z-ip ip 1 am z-am am-r am-l z-am z-am am il z-i il-r il-l z-il z-il il 1 3 1 3 3 z z i z z z z z z OF z z z z z z z z µ λ λ λ λ λ κ λ µ λ λ λ λ λ κ λ µ λ λ λ λ λ κ - = - -     = - + - - + - +             - - + - -         - - + ∑ l 1 z i λ     -           ( 42)
The sixth objective function is built by requiring that the average lateral stresses of both interphase layers, which are calculated by integrating the lateral stress distributions from MD simulations, are as close as possible. Mathematically speaking, ( ) ( )

800 ip-r ip-r ip-l ip-l 6 lat-xx lat-yy lat-xx lat-yy i i OF σ σ σ σ = + - + ∑ ( 43
)
And finally, the last objective function is defined as follows ( )

800 am am am 7 lat xx / 2 yy i OF σ σ σ = - + ∑ ( 44)
where am lat σ denotes the lateral stress of the amorphous phase calculated by substituting am µ and am κ into the first relationship of the system of equations ( 30), and

( )

am am xx / 2 yy σ σ +
is the arithmetic mean of the amorphous lateral stresses averaged over its identified domain. The explanation for this particular selection is provided in Subsection 3.6.

A multi-objective optimization tool is the most suitable approach for finding the global minimum of the above problem. In the absence of such a tool, we opted to fusing all the objective functions into a single nonnegative objective function as follows

OF OF OF OF OF OF OF OF

= + + + + + + ( 45)
and the combination of "interior point" and "SQP" algorithms of the constrained nonlinear minimization tool of MATLAB ® software is used for finding the optimal point. The objective functions, however, are usually accompanied by a set of linear and / or nonlinear constraints in the form of equalities and / or inequalities. For the problem considered here, the following logical constraints are applicable:

( ) ( ) ( ) ( ) ( ) ( ) am-r am-l il-r il-l 0 0 0 0 am/il
z z z z z z z z z z z z z z z z µ κ < < - < < - < < - < < - > > = = = = = = ( ) ( ) l il-r il-r il-l il-l 800 800 800 Max Min , , z z z z = = ( 46)
The bounds of the amorphous and interlamellar domains at the no loading state are determined from the density distribution at the no loading state, Figure 19, as well as the distributions of the lateral stress at the very first loading increment. The bounds of the extremes of the boundaries come from the diagrams of Subsection 3.6. Minimization of tot OF and finding the optimal point is equivalent to the hyperelastic characterization of the domains involved as well as approximating the evolutions of the boundaries in question.

Results and discussion

3.5

Enforcing the continuity between the arguments of il-r z , am-r z , il-l z and am-l z , as introduced in relation ( 31), reduces the number of design parameters, i α s, from 48 to 40, which is yet rather high and demanding in terms of CPU time for the adopted combinatory optimization tool. As another drawback, the optimization tool that is employed here shows a strong dependence on the starting point. To circumvent this problem, the approximate solution obtained from the Genetic Algorithm optimization module after several generations can be used as the starting point.

α α α α α α α α α α α α = = = - = = = = = = = = - = i α , 13
α α α α α α α α α α α α = = = - = = = = - = = = = - = i α , 25
α α α α α α α α α α α α = - = - = = - = - = - = = - = - = - = = - i α ,
α α α α α α α α α α α α = - = - = = - = - = - = = - = - = - = = -
The initial length of the hyperelastic elements are calculable form the intercepts of the amorphous and interlamellar boundaries. Accordingly, the initial length of the interlamellar, amorphous and interphase domain identified in this study are 106.7 Å, 69.1 Å and 18.8 Å, respectively.

Having applied the above explained methodology together with the assumed neo-Hookean behavior, the hyperelastic parameters of the three domains of interest are identified as follows:

µ µ µ κ κ κ =  = =       = = =     ( 47)
Given that in the neo-Hookean model, the shear parameter is an index of the rigidity of the medium, the expectation that am ip µ µ

< is obviously satisfied because the interphase layer is known to be stiffer than the central amorphous phase. Furthermore, in the existing studies where the hyperelasticity of the noncrystalline domain is incorporated into the large deformation formulation of PE, it is assumed that at room temperature il 35MPa µ = [START_REF] Van Dommelen | Micromechanical modeling of the elasto-viscoplastic behavior of semi-crystalline polymers[END_REF][START_REF] Sedighiamiri | Micromechanical modeling of the deformation kinetics of semicrystalline polymers[END_REF][START_REF] Sedighiamiri | Micromechanics of semicrystalline polymers: Yield kinetics and long-term failure[END_REF]. Indeed, Sedighiamiri et al. [START_REF] Sedighiamiri | Micromechanical modeling of the deformation kinetics of semicrystalline polymers[END_REF] assert that the best fit to experimental results is seen when il 35MPa µ =

. Keeping in mind that the stiffness of amorphous polymers in the rubbery regime weakly increases with temperature, one would have expected values slightly higher than 35MPa for the interlamellar domain and much higher for the interphase layers. Rerunning MD simulation for a larger PE stack might help to fulfill this expectation. Although the bulk modulus is not as significant as the shear parameter, the identified bulk moduli are also in agreement with the mean values of 2 GPa and 3 GPa used elsewhere [START_REF] Van Dommelen | Micromechanical modeling of the elasto-viscoplastic behavior of semi-crystalline polymers[END_REF][START_REF] Bédoui | Micromechanical modeling of isotropic elastic behavior of semicrystalline polymers[END_REF][START_REF] Sedighiamiri | Micromechanical modeling of the deformation kinetics of semicrystalline polymers[END_REF][START_REF] Sedighiamiri | Micromechanics of semicrystalline polymers: Yield kinetics and long-term failure[END_REF]. As discussed in Subsection 3.6, this last parameter is identified with a lower degree of certainty as compared to the shear modulus.

To evaluate the success of the characterization algorithm in attaining the defined goal, that is finding the minimum of the positive objective functions, the boundaries and stretches of the left and right interphase layers are plotted in Figure 22 and Figure 23, respectively. Evidently, the evolutions of the four boundaries look qualitatively like those suggested by the plots of P2 parameter from Figure 12(b) of the underlying MD study [START_REF] Lee | Plastic deformation of semicrystalline polyethylene by molecular simulation[END_REF]. Additionally, if their respective intercepts are subtracted, the pair boundaries of the amorphous and interlamellar domains are almost each other's mirror, as one might have anticipated. As another sign of the successful optimization, the diagrams of ip-l λ and ip-r λ in Figure 23 of ip λ too, plotted by substituting the identified ip µ and ip κ into the second relationship of Eq.( 30) agrees very well with the two other ones. And finally, the evolutions of the average lateral stresses integrated over the two interphase layers are plotted in Figure 24. The two diagrams are in good agreement with slight, acceptable differences at medium and large strains that are probably due to the limited size of the simulation stack. Moreover, since smooth functional forms have been adopted for probably noisy boundaries, naturally a higher level of noise is pushed into the diagrams of ip-l lat σ and ip-r lat σ

. In other words, the disagreement between the two diagrams of Figure 24 may be further attenuated by allowing the boundaries to behave less smoothly. is not reliable enough for identification. Whether this disagreement is due to the small size of the simulation stack or has other reasons behind, warrants further simulations. Moreover, the melting phenomenon that has been observed during MD simulations at the interface of the crystalline-noncrystalline region is successfully captured by the presented continuum mechanics analysis. To elaborate, it is worth recalling that the plasticity of the crystallites produced by slip as well as the deformation mode of the entire PE stack is isochoric, while the noncrystalline domain deforms compressibly, simply because am/il/ip κ ≠ ∞ . Since the entire stack deforms isochorically, the volume change due to the interlamellar domain is compensated either by void creation or by phase transformation. The former has been clearly observed in a separate MD simulation carried out under a faster mode of deformation and the latter has taken place in the deformation mode under survey [START_REF] Lee | Plastic deformation of semicrystalline polyethylene by molecular simulation[END_REF]. Therefore, melting and/or recrystallization at the crystalline and noncrystalline interface should take place to accommodate the imposed volume conservation constraints.

The occurrence of such phenomena changes the current length of the hyperelastic element in addition to the elongation produced by the mechanical loading. In the adopted constitutive law, however, the elongation due to the mechanical loading is concerned and there is no need to take into account the change in the size of the hyperelastic element due to melting. Properly speaking, the logic of the problem and governing equations necessitates that am/il/ip z λ - appearing in the constitutive law relates to the same element that is initially marked out as the hyperelastic element without any mass exchange through their boundaries. Keeping this in mind, if melting occurs during the loading and part of the crystallites near the interface transforms into the noncrystalline phase, then the real, new interlamellar boundaries will be due to the stretching of the interlamellar domain engendered by loading plus the change in size due to the phase transformation. For characterization purposes and for the governing equations to remain valid, the change in length due to the phase transformation, if any, must be excluded from the evolution of the boundaries. For the problem examined in this study, the dominant phenomenon is melting; therefore there is no need to revise the initial length of the hyperelastic elements. Additionally, through a clear mathematical expression of the physics of the problem and defining relevant objective functions along with logical constraints, the need for reading the current position of the amorphous/interlamellar boundaries at each loading increment is avoided and these boundaries are identified along with the hyperelastic coefficients.

Another phenomenon captured by the presented continuum mechanics analysis is the rotation of the crystalline chains. In the basic MD simulation, the crystalline stems, which are initially oriented at ~ 30° with respect to the interface normal, start to rotate with loading to get aligned with the applied stress. This rotation continues up to 0.21 The density distribution at no loading state and lateral stress distributions at the very first loading increment can be utilized for locating the centers of the interphase layers. Two straight lines representing the average density values of the crystallites and the central amorphous phase are drawn (Figure 19) whose median intersects the density distribution at -45 Å and 43 Å. Moreover, the lateral stress distributions at the very first increment of loading have two picks at -45 Å and 43 Å (Figure 28A). It can therefore be concluded that the centers of the left and right interphase layers are located at -45 Å and 43 Å. Given that the interphase thickness varies from 10-20 Å [START_REF] Flory | The interphase in lamellar semicrystalline polymers[END_REF][START_REF] Veld | Temperature-dependent thermal and elastic properties of the interlamellar phase of semicrystalline polyethylene by molecular simulation[END_REF] distribution. On the other hand, it is evident that the amorphous lateral stresses show a weak dependence on the integration interval as compared to the interlamellar lateral stresses. This is the reason why the seventh objective function,

7
OF , is defined accordingly.

To evaluate the certainty of the identified hyperelastic parameters, the analytical diagrams of the lateral stress distribution for typical values of hyperelastic coefficients are plotted in Figure 31. In other words, typical shear and bulk moduli are substituted in the second relationship of the system of equations ( 30) µ < , the lateral stress distribution diagrams show a weak sensitivity and as such, the identified coefficient has a lower degree of certainty.

The identified shear moduli for the three hyperelastic domains of interest are greater than 30 MPa whereas the identified bulk moduli are around 2.7 GPa. Therefore, it can be asserted that the shear moduli are identified with a high degree of certainty whereas the identified bulk moduli have a lower degree of certainty. PLASTICITY of semicrystalline polyethylene has been extensively investigated both experimentally [START_REF] Frank | Deformation processes in polyethylene interpreted in terms of crystal plasticity[END_REF][START_REF] Seto | Phase transformation and deformation processes in oriented polyethylene[END_REF][START_REF] Young | Deformation mechanisms in oriented high-density polyethylene[END_REF][START_REF] Bowden | Deformation mechanisms in crystalline polymers[END_REF][START_REF] Haudin | Plastic deformation of amorphous and semi-crystalline materials[END_REF][START_REF] Bartczak | Deformation mechanisms and plastic resistance in single-crystal-textured high-density polyethylene[END_REF][100] and theoretically [START_REF] Van Dommelen | Micromechanical modeling of the elasto-viscoplastic behavior of semi-crystalline polymers[END_REF][START_REF] Lee | Simulation of large strain plastic deformation and texture evolution in high density polyethylene[END_REF][START_REF] Parks | Polycrystalline plastic deformation and texture evolution for crystals lacking five independent slip systems[END_REF][101][102][103]. The interest in the subject has been prompted by the daily increase in the structural applications of semicrystalline polymers, notably semi-crystalline polyethylene. Large plastic deformations in semi-crystalline polymers usually lead to the development of high anisotropy and preferential texturing simply due to their semi-crystalline microstructure. This new texturing and anisotropy may particularly take place during their forming processes where large strains are involved. The new anisotropies achieved in this way are sometimes quite desirable. Semi-crystalline polymers that are processed into fibers and ribbons with excellent mechanical performance in preferred directions are examples of desirable anisotropies.

The significance of the concept of semi-crystallinity is better appreciated when one notices that a purely amorphous polyethylene would be a soft, viscous substance and on the other side, a purely crystalline polyethylene sample would be extremely fragile. The polyethylene samples are, virtually, composites of both phases. The proportion of the constituents, their relative arrangements, and the level of connectivity between them decide the effective mechanical properties of the polyethylene sample. The end-use property of the semi-crystalline polyethylene depends on the interplay between the crystalline and noncrystalline domains. Properly speaking, the mechanical performance of semi-crystalline polyethylene is a function of the underlying microstructure. The most important microstructural features include crystallinity, initial crystallographic texture and molecular orientation. Experimental studies on the plasticity of semi-crystalline polymers have by and large explored the topic from the perspective of the morphological changes, the identification of deformation mechanisms, and the characterization of deformation resistances. The plastic deformation of semicrystalline polymers is in principle of crystallographic nature although being very complicated when considered on the scale of crystallite dimensions and from the viewpoint of interplay between the crystallinenoncrystalline phases. The mechanisms involved in the plasticity of semicrystalline polymers have been comprehensively reviewed by Bowden and Young [START_REF] Bowden | Deformation mechanisms in crystalline polymers[END_REF], Haudin [START_REF] Haudin | Plastic deformation of amorphous and semi-crystalline materials[END_REF] and Galeski [108], among many others. Of the existing semicrsytalline polymers, the data for the high density polyethylene is quite abundant basically owing to its simple molecular structure and ease of preparing samples with high degree of crystallinity. Single crystal state of high density polyethylene has been approximated by Bartczak et al. [START_REF] Bartczak | Deformation mechanisms and plastic resistance in single-crystal-textured high-density polyethylene[END_REF] with bulk samples compressed under plane strain state in a channel die. Large deformation of high density polyethylene has also been experimentally investigated under a variety of deformation modes including uniaxial tension and compression [100, 109] as well as simple shear [110,111]. The distinguishing feature of the plasticity in semi-crystalline polymers is the inextensibility of chain segments in the crystallites. Indeed, the folded chain structure of the crystalline phase constrains the arbitrary plastic deformation of crystallites. This means that there are less than five independent slip systems to accommodate arbitrary plastic deformation. This understanding needs to be taken into account when analyzing the plasticity of semi-crystalline polymers numerically.

In the available numerical studies that have dealt with the simulation of the plastic deformation and texture evolution in semi-crystalline polyethylene, the principal differences lie in the localization law and the constitutive laws adopted for each component. In this study, for the localization law, use is made of the complete version of modified Taylor approach proposed by Ahzi et al. [START_REF] Ahzi | Modeling of plastic deformation evolution of anisotropy in Semi-Crystalline Polymers[END_REF] together with Sachs approach. On the other hand, the viscoplastic constitutive law adopted for the crystalline phase is no different from what is chosen for this phase in most of the relevant studies. For the noncrystalline phase, however, three different statistical mechanics constitutive laws, viz. 3-chain model, 8-chain model and Gent model, are adopted to model the back stress and final results in terms of the number of rigid links required for each model are compared. As an objective, this chapter aims to determine which model requires a more realistic value for the number of rigid links.

In this reexamination, the stress-strain response and texture evolution of an aggregate of polyethylene subjected to finite strains is revisited. For simplification purposes, the elasticity and pressure sensitivity of the phases are neglected and it is assumed that each phase exhibits a viscoplastic behavior. As another simplifying assumption, the interphase layer and central amorphous domain are treated as a single unified domain. This assumption reduces the complexity level of developing the associated code without any appreciable compromise in the accuracy of the results. However, in case the hyperelastic properties of each of these domains are provided separately, the homogenized hyperelastic properties of the entire noncrystalline domain can be identified using the ideas of Chapter III, and then plugged into the governing equations.

Constitutive laws 4.2

The constitutive laws adopted for the crystalline and noncrystalline phases are identical with those presented in [START_REF] Ahzi | Modeling of plastic deformation evolution of anisotropy in Semi-Crystalline Polymers[END_REF]. The linear elasticity and pressure sensitivity of the components is neglected and their finite deformation is assumed to be strain ratedependent.

Viscoplasticity of the crystalline phase

The plastic mechanisms of deformation in polyethylene crystallites are crystallographic slip, twinning and stress-induced martensitic transformations. Of these mechanisms, the deformation due to slip is dominant, which is considered to be the only driving mechanism in this work. For the large deformation of crystalline phase in HDPE, the viscoplastic power law relationship correlates the microscopic strain rate and the resolved shear stress. Following the previous similar studies [START_REF] Ahzi | Modeling of plastic deformation evolution of anisotropy in Semi-Crystalline Polymers[END_REF][START_REF] Parks | Polycrystalline plastic deformation and texture evolution for crystals lacking five independent slip systems[END_REF][112][113][114], the following rate-dependent constitutive law is adopted for the crystallites

1 0 c n g g α α α α α τ τ γ γ - = ɺ ɺ ( 49)
where 0 γɺ is a reference strain rate, α τ ( )

g α ≤
is the resolved shear stress on system α , c n is a non-linear rate exponent and g α is the shear resistance of slip system α , which is given in Table 4 in front of each slip system. As discussed in chapter I, the dominant crystalline structure of high density polyethylene is orthorhombic with the lattice parameters a=7.4 Å, b=4.93 Å and c=2.54 Å, where c is the crystallographic axis coinciding the chain direction. Experimental results indicate that there exist four linearly independent slip systems in the orthorhombic crystal structure of HDPE, summarized in Table 4. The slip system strength, on the other hand, may evolve with deformation as follows

1 K g h α αβ β β γ = = ∑ ɺ ɺ ( 50
)
where h αβ are the components of the matrix of hardening moduli. Since the thin crystalline lamellae cannot retain dislocations, the strain hardening is neglected.

Additionally, the shear resistances, g α , are assumed to remain constant with the ongoing deformation.

        1.66 2.5 2.2

Viscoplasticity of the noncrystalline phase

Similar to the crystalline phase, a power law relationship is assumed between the local shear rate and its corresponding shear stress. This power law is an approximation of the exponential law proposed by Argon [115] for the glassy state of amorphous polymers. The simple viscoplastic power law relationship correlating the plastic shear rate, il γɺ , and the effective shear stress of the rubbery noncrystalline phase, il τ , in semicrystalline polyethylene reads

il il il 0 il 0 n a τ γ γ τ   =     ɺ ɺ ( 51)
In the above relationship, il n and 0 γɺ denotes the rate exponent and the reference strain rate,respectively, which are set equal to their counterparts in the crystalline phase, for simplicity and without loss of any generality. As such, the symbol n is chosen to represent the rate exponent of both phases, that is il c n n n = = . If 0 τ is set to represent the shear resistance of the easiest slip system in the crystalline phase, the dimensionless coefficient il a can be adjusted such that the product il 0 a τ equals the reference shear strength of the amorphous domain.

Application of the mechanical loading to the molecular chains in the noncrystalline region make them start getting aligned in the direction of the maximum stretch leading to anisotropy in the resistance of the noncrystalline domain to plastic deformation. To take this effect into account when deriving the tensorial constitutive law, a back stress tensor, il H , is introduced into the flow rule, following Boyce et al. [116]. This term is responsible for including the changes that the alignment of the molecular chains brings about in the plastic resistance of the interlamellar zone [START_REF] Ahzi | Modeling of plastic deformation evolution of anisotropy in Semi-Crystalline Polymers[END_REF]. Letting il D and il S denote the strain rate tensor and deviatoric Cauchy stress tensor of the noncrystalline domain, respectively, the driving stress created in the noncrystalline domain is defined as il il -S H . Therefore, the resolved shear stress, il τ , is defined as the norm of the driving stress:

( )( ) il il il il il 1 2 ij ij ij ij S H S H τ = - - ( 52)
Substitution into Eq.( 51), the tensorial form of the associated constitutive law reads as follows:

1 il il il il il 0 0 0 n a a γ τ τ - - - = S H S H D ɺ ( 53)
The form of the back stress model, il H , is defined according to the chain network model adopted for the interlamellar domain.

Back stress models

In modeling rubber elastic materials, there exist two perspectives: statistical mechanics (or micromechanically motivated) treatment and invariant-based or stretchbased continuum mechanics treatment. The phenomenological (or continuum-based) models which are based on the invariants of the left Cauchy-Green deformation tensor, lack a direct physical linkage to the underlying mechanisms of deformation. Of such models, one can mention Mooney, Rivlin, Valanis and Landel, Yeoh, and so forth [START_REF] Treloar | The Physics of Rubber Elasticity[END_REF]117]. The most sophisticated one belongs to Ogden [118] who proposes the following strain energy density function ( )

1 2 3 3 n n n n n n G W α α α λ λ λ α = + + - ∑ ( 54)
where n G and n α are fitting parameters and may take on any value including non- integers. The sum in Eq.( 54) is expanded as needed to achieve a reasonable fit to the data. According to Twizell and Ogden [119], stability of Eq.( 54) requires 0

n n G α > .
Additionally, Arruda and Boyce [120] argue that the Ogden model is an empirical relationship that requires more than one experiment to obtain the number of coefficients required for capturing the state of deformation dependence.

Although the phenomenological models have merits and are frequently used in investigating rubber deformation problems, many of them fail in successful description of the rubber response under deformation modes other than uniaxial tension, without changing the model parameters. In other cases, the models are either mathematically complicated or the number of required parameters is unreasonably high. Therefore, models of molecular chain network sound to be a better alternative. In the Gaussian statistics approach, as the first statistical mechanics attempt developed for describing the force-deformation relationship in a polymeric network, it is assumed that the chains never near their fully stretched length, lN , where N is the number of statistical rigid links of equal length l in the chain between chemical crosslinks [START_REF] Treloar | The Physics of Rubber Elasticity[END_REF]. The "rigid link" is considered to be that segment of the real chain that undergoes rigid body motion when the rubber material is deformed. The rigid link may consist of one or more monomers, depending on the rubber material under study. The strain energy density function proposed by the Gaussian statistics approach is of the following form ( )

2 2 2 r 1 2 3 r 1 3 , 2 W C C qkT λ λ λ = + + - = ( 55)
where i λ s are the applied stretches and the rubbery modulus r C is a function of the chain density, q , Boltzmann's constant, k , and absolute temperature, T . The stress- stretch relationship is obtained by differentiating the strain energy with respect to the stretch. In an incompressible rubber, the principal stresses can be obtained from the work of deformation as follows

d d i i i W p σ λ λ = - ( 56)
where i λ are principal stretches and the pressure term may be determined from the boundary conditions. The energy function ( 55) is derived with the assumption that the current length of the chain, r, remains well below the fully stretched length of the chain, Nl . At large deformations, i.e. / 0.4 r Nl > , the depart from the Gaussian model becomes more evident and the non-Gaussian nature of the chain stretch prevails [121].

The subsequent modified versions of chain statistics were an effort to allow for larger stretches than are afforded by the Gaussian statistics. A good review of these models are provided in [START_REF] Treloar | The Physics of Rubber Elasticity[END_REF]. The parameters that all these models have in common are the rubbery modulus, r C , and the chain locking stretch, L λ . The chain locking stretch is the chain stretch at which its full extensibility limit is reached. From the statistical discussions, the initial chain length is found from random walk statistics as Nl where in conjunction with the fully extended length, Nl , yields

L Nl N N l λ = = ( 57)
The mathematical models of molecular chain network in polymers are based on Langevin chain statistics. Kuhn and Grün [122] used this theory and derived the following formulation for the work of deformation due to the stretching of unstretched chains:

1 1 1 chain chain chain chain ln ln sinh r r r r W kTNq Tc Nl Nl Nl Nl - - -             ′ = + - -                         L L L ( 58
)
where q is the chain density, c′ is a combination of constants, and

1 -
L denotes the inverse Langenvin function. For further details on the inverse Langevin function, cf. Appendix I.

Wang and Guth made use of the Langevin statistics to present their three-chain model [123]. In this model, the representative cubic unit cell is made of three chains, as is schematically illustrated in Figure 32. This model assumes that the faces of the cube element are aligned with the principal stretch space during deformation. The stressstretch relationship in the three-chain model is given by

1 r L L 1 3 i i i C p λ σ λ λ λ -  = -     L ( 59
)
Although this model captures a small state of deformation dependence of the behavior at small stretches, it yields nearly identical predictions for all deformation modes at large stretches [121]. The results of this model are controlled by the contribution of chains that lie along the maximum principal stretches. Therefore it cannot be a good representative of the network response because it samples only the chains that lie along the principal stretch directions. The failure of this model to capture the state of deformation dependence lies in the absence of a cooperative nature of network chain deformation apart from the imposition of incompressibility conditions. This drawback was the major motivation for the introduction of the eight-chain model where the chains deform more cooperatively. In the eight-chain model proposed by Arruda and Boyce [120], it is assumed that the molecular network of any amorphous polymer can be represented by a model of eight chains, as illustrated in Figure 33. These eight molecular chains are attached to each corner of a cube at one end and join up at the center of the cube at the other end. This model requires an initial modulus and a "limiting chain extensibility" as the only parameters for describing the network. These two parameters are linked to the physics of molecular chains that contribute to the deformation. Due to the symmetry of the chain structure, the interior junction point remains in the center of the cube throughout the deformation and the stretch of each chain becomes the root mean square of the principal stretches The strain energy function of this model is given by

1 1 1 chain chain chain chain r ln ln sinh W C N N N N N λ λ λ λ - - -             = + -                           L L L ( 61)
where in conjunction with relationship ( 57) results in the following relation for principal Cauchy stresses

2 8ch 1 chain r chain 1 3 i i N C p N λ λ σ λ -  = -     L ( 62
)
Later on, Gent proposed the following empirical and less sophisticated relation for

W [124]: r 1 ln 1 2 m m C J J W J   = - -     ( 63)
Here, m J denotes the maximum value of

Gent 2 m r m 1 i i J C p J J σ λ = - - ( 64)
As compared to its predecessors, Gent model is known to be more successful. The results of the three models are compared in Subsection 4.6.

Composite inclusion 4.3

Following the leading work by Ahzi and co-workers [START_REF] Ahzi | Modeling of plastic deformation evolution of anisotropy in Semi-Crystalline Polymers[END_REF][START_REF] Ahzi | Estimates of the overall elastic properties in semicrystalline polymers[END_REF] and many others later, the composite inclusion model is used for modeling the crystalline and noncrystalline phase as a single inclusion. For modeling purposes, the crystalline and noncrystalline phases are thought as plate-like components perfectly bonded to one another. In the undeformed state, the crystalline-noncrystalline composite inclusions in HDPE are arranged, as described in Chapter I, in a spherulitic morphology. This initial configuration of polyethylene composite inclusions can be well approximated by the random distribution of composite inclusions in an aggregate of HDPE. This underlying assumption defines the basis of the initial texture considered in the subsequent simulations for investigating the finite deformation and texture evolution of HDPE aggregates under different modes of deformation.

The generalized form of the composite inclusion idea, i.e. the Extended Composite Inclusion Model, is described in Appendix E, where the effective stiffness of a three-layer inclusion is derived as a function of the stiffness and volume fraction of its constituents. As explained therein, in composite inclusion model, the compatibility and equilibrium at the interface ought to be fulfilled. Enforcement of compatibility requires the continuity of velocity across the crystalline-amorphous interface. This continuity condition means that if 0

x denotes an arbitrary vector lying in the interface plane such that

I 0 0 = i n x
, the interface compatibility condition can be expressed as ( )

Cr il 0 - = L L 0 x ( 65
)
Let I i e denote the orthonormal basis fixed to inclusion I, with I 3 e being co-directional with the normal to the interface of the two phases, I n . In this basis, the compatibility conditions in conjunction with the incompressibility conditions for both phases, necessitates that the following continuity equalities for the strain rate and spin components hold

W W W = = ( 68)
In Eq.( 66) the indices , 1, 2 α β = . On the other hand, for the tensors of strain rate and spin, the average theorem requires that ( )

I il Cr il il 1 η η = + - D D D ( 69)
( )

I il Cr il il 1 η η = + - W W W ( 70)
where in conjunction with relation ( 65) lead to the following relationships for the remaining components of strain rate and spin tensors that do not appear in Eqs.( 66)-

( )

I il Cr 3 il 3 il 3 1 D D D α α α η η = + - ( 71)
( )

Cr I Cr il 3 3 il 3 3 W W D D α α α α η = - - ( 72) ( )(
)

il I Cr il 3 3 il 3 3 1 W W D D α α α α η = + - - ( 73)
The force equilibrium at the interface necessitates the following equalities between the two corresponding components of deviatoric stresses in the I j e basis

I Cr il 3 3 3 S S S α α α = = ( 74)
The other components of the deviatoric stresses are correlated through the average theorem relation ( )

I il Cr il il 1 η η = + - S S S ( 75)
It is worth noting that for the normal component of the interface traction, it can be equilibrated by assuming that any jump in the normal component of the deviatoric stress is balanced by a corresponding jump in the pressure, that is il il cr cr [START_REF] Hahn | Dielectric relaxation of the crystal-amorphous interphase in poly(vinylidene fluoride) and its blends with poly(methyl methacrylate)[END_REF] 33

S p S p -= - ( 76)
where il p and cr p are the pressure available in the amorphous layer and crystalline lamella, respectively. Fortunately, the above pressure components do not appear in the solution algorithm because of the pressure insensitivity of the constitutive laws adopted for the phases involved.

Localization/interaction law 4.4

In any inhomogeneous aggregate, the local stress and strain of each inclusion is correlated to the global stress and strain of the aggregate through establishing a localization (or interaction) relationship. Proposing such a relationship has been the subject of numerous studies in crystal plasticity and mechanics of heterogeneous materials. Since finding the exact solution with all conditions of local compatibility and equilibrium for all grains being satisfied is impossible, approximate approaches have been proposed. The simplest and most well-known ones are Sachs model [125] and Taylor model [126]. In Sachs model, the local and global equilibrium are satisfied simply by assuming uniform stress everywhere. Global compatibility is enforced as a global volume average while local compatibility may be violated. Sachs approach in plasticity of inhomogeneous materials is the analogue of Reuss approach in elasticity of inhomogeneous composites.

The analogue of Voigt scheme in large deformation plasticity is Taylor approach which assumes that the plastic strain is identical everywhere in the polycrystal and equal to the macroscopic strain [126]. According to the assumptions of this approach, the conditions of local and global compatibility are fulfilled whereas the equilibrium at the grain boundaries is not generally satisfied. Molinari et al. [127] argue that in predicting the texture evolution, Taylor approach has often proved successful for metals with FCC and BCC Bravais lattices. This success is attributed to the close deformation resistances of the available slip systems [START_REF] Parks | Polycrystalline plastic deformation and texture evolution for crystals lacking five independent slip systems[END_REF]. Modified versions of Taylor approach have been proposed in the context of crystalline slip to investigate different slip resistances in different slip systems and to assess the strain rate sensitivity of slip systems [112,128], or to include the effects of elastic deformations as well as the hardening in the slip [114]. As the principal disadvantage, such models fail to capture the effects of grain geometry and orientation.

Due to its relatively better success in predicting the texture evolution in semicrystalline polyethylene, a modified version of Taylor approach is employed here as the localization law. Introduced first by Ahzi et al. [START_REF] Ahzi | Modeling of plastic deformation evolution of anisotropy in Semi-Crystalline Polymers[END_REF], their modified Taylor approach addresses the case when ψ , the angle between crystalline chain segments and normal to the crystalline-noncrystalline interface, falls below a threshold value, 0 ψ . In other words, in normally prepared polyethylene samples, ψ takes values around 30°. This configuration is supported by the molecular simulation study by Gautam et al. [START_REF] Gautam | Molecular simulations of the interlamellar phase in polymers: Effect of chain tilt[END_REF] which provides a thermodynamic explanation for the frequent observation that the crystalline-noncrystalline interface is parallel to { } 201 plane. However, during the deformation and as the texture evolves, the angle ψ in some grains might approach 0°.

Since the deformation of the crystallites in the chain direction remains inextensible, original Taylor approach yields unrealistic predictions for grains in which ψ is close to zero. To address such grains, Ahzi et al. [START_REF] Ahzi | Modeling of plastic deformation evolution of anisotropy in Semi-Crystalline Polymers[END_REF] proposed a modification to relax the strain rate component imposed in I n direction while conserving the global compatibility conditions. However, it is mathematically demonstrated that this modification is not complete and further modifications are required to include the case when ψ approaches the right angle.

Mathematical formulation of Modified Taylor

To formulate the modified Taylor approach of this work, let L denote the global macroscopic velocity gradient of the aggregate. The traceless, second-order tensor L can be uniquely and additively decomposed into two symmetric and skew-symmetric tensors as follows ( ) ( )

T T 1 1 2 2 = + + - = + L L L L L D W ( 77)
where symmetric D is the macroscopic deformation rate tensor and the skewsymmetric W denotes the macroscopic spin tensor. The assumption of Taylor approach requires that the deformation rate tensor for all inclusions are identical and equal to the macroscopic one, D. But we know that depending on the relative orientation of the chain segments in each crystalline lamella and its associated normal vector to the crystalline-noncrystalline interface, the imposed D might be impossible. Therefore, one needs to modify the portion of D that is allotted to such inclusions, hence the designation "Modified Taylor approach". Let us denote this share of D for inclusion "i" by I i D ɶ which is the closest projection of the initial D such that the violation of the kinematic constraints of inclusion "i" is avoided. To further clarify, let us assume, without loss of generality, that the convected coordinates system of inclusion "i", which is attached to the interface of this inclusion and indicates its orientation, is obtainable by simply rotating the abc-axes of the crystallite in the counterclockwise direction about baxis by angle i ψ . In this way, the c-axis of the crystallite, which indicates the direction of the chain stems, take the following representation in the convected coordinates system: [ ]

T i i sin 0 cos ψ ψ = c ( 78
)
As explained by Parks and Ahzi [START_REF] Parks | Polycrystalline plastic deformation and texture evolution for crystals lacking five independent slip systems[END_REF], the inextensibility of the chain stems in the crystallites is equivalent to

Cr i 0 D C = i ( 79
)
where Cr i D denotes the rate deformation tensor of the crystallite in inclusion "i" and C is a second-order tensor defined as the dyadic product of c, that is ⊗ C = c c . Expansion of the left-hand side of Eq.( 79) in inclusion axes by substitution from [START_REF] Jaglinski | Negative Stiffness and Negative Poisson's Ratio in Materials which Undergo a Phase Transformation[END_REF]

D D D D D D D D D D D D ψ ψ ψ ψ ψ ψ ψ ψ ψ ψ         = = ⇒             + + = C D D C = ɶ ɶ ɶ ɶ ɶ ɶ ɶ i ( 80)
According to the above relationship, when i 0 ψ → , i I 33 D ɶ must go to zero otherwise a mathematical error will be produced. Similarly, when i 90 ψ → , i I 11 D ɶ must go to zero otherwise it would not be possible to impose the inextensibility condition. Accordingly, the projection tensor i I P that projects D for inclusion "i" to produce I i D ɶ is a function of i ψ and should be constructed such that the two aforementioned limit conditions are satisfied. Thus the following three-argument mathematical representation is proposed for i I P :

i i i i I i i i 0 90 otherwise ψ ψ  -⊗ →  = -⊗ →    N N X X I P I I ( 81)
In the above projection rule, i i i

1 3 = ⊗ - N I n n
represents the deviatoric part of the dyadic product of the normal to inclusion "i", and similarly,

i i i 1 1 1 3 = ⊗ - X I e e
represents the deviatoric part of the dyadic product of i 1 e which is the unit vector of the inclusion axes, co-planar with i a , i c and i n . For numerical programming, however, a definite numerical value is needed to replace the concept of mathematical limit appearing in the previous relationship. To this end, the proposed projection is re-written in the following form

i i i i 0 I i i i 0 90 otherwise ψ ψ ψ ψ  ⊗ <  = ⊗ > -    N N X X I - P I - I ( 82)
where 0 ψ is a small angle close to 0°. Ahzi et al. [START_REF] Ahzi | Modeling of plastic deformation evolution of anisotropy in Semi-Crystalline Polymers[END_REF] proposed 15° as the threshold value, after running a large number of simulations. It should be mentioned that, as explained in [START_REF] Ahzi | Modeling of plastic deformation evolution of anisotropy in Semi-Crystalline Polymers[END_REF][START_REF] Parks | Polycrystalline plastic deformation and texture evolution for crystals lacking five independent slip systems[END_REF], the global compatibility requirement leads to the following normalization , is satisfied. Additionally, with some mathematical manipulation, the three-argument fourth-order tensor proposed in [START_REF] Kochmann | Analytical stability conditions for elastic composite materials with a non-positive-definite phase[END_REF] for calculating the local deformation rate tensor for inclusion "i" can be reduced to a single argument projection as follows

i I i i i i i i 0 0 90 3 3 exp exp 2 2 m m ψ ψ ω ω ψ ψ         -     = - - ⊗ - - ⊗                 N N X X P I ( 84
)
where for computational purposes 1

ω ≫ and m is of the order of 50.

Tips and remarks on implementation issues 4.5

As was discussed earlier, a macroscopically isotropic aggregate of polyethylene is assumed to be composed of a large number of composite inclusions that are uniformly randomly scattered in space. Such a scattering is equivalent to a uniformly random distribution of the inclusions' orientation. For such a distribution, one needs a uniformly random distribution of Euler angles for either the initial texture or the initial morphology. Here, we opted to create an initial random texture and let the morphology distribution be a function of the texture. The initial distribution of Euler angles for the texture is created as follows: The notation used here for Euler angles is similar to the one adopted in Appendix A and the same convention is respected for setting up the transformation tensor. If for inclusion "i" the angle straddling its c and n is denoted by i ψ , then the transformation tensor from the I j e basis of the inclusion to aggregate axes, i b , is obtained as follows

[ ] ( ) [ ] [ ]
( ) ( ) ( ) ( ) cos 0 sin : with 0 1 0 sin 0 cos i i i i i i i i ψ ψ ψ ψ     = =     -   b a r r ( 86
)
where i a is the transformation tensor from the crystalline axes of inclusion "i" to macroscopic axes. It should be noted that i ψ is measured from c to n in the clockwise direction.

Moreover, the numerical codes for the implemented schemes have been developed in MATLAB ® . This programming language is selected because of its powerful optimization tool. In other words, the codes are engineered from a system identification perspective at different levels. The constitutive equations for each phase in each inclusion together with the equations of local compatibility and equilibrium are solved by defining a set of relevant nonnegative objective functions and then finding their minima using the optimization toolbox of MATLAB ® . The requirement for the global compatibility and equilibrium, on the other hand, is also satisfied by finding the minima of a new set of appropriately defined nonnegative objective functions. The advantage of using the optimization tools lies in its capability to find the sought for minima in the presence of constraints, if any, no matter what their nature is. In this example, the most important constraint is imposed on the resolved shear stress of every slip system which is required to remain less than or equal to its corresponding shear resistance. Finally, the issues concerning the texture and morphology updating are treated as explained in earlier studies such as ref. [START_REF] Lee | Micromechanical modeling of large plastic deformation and texture evolution in semi-crystalline polymers[END_REF].

Results and discussion

4.6

The diagrams of The diagrams of Figure 34 suggest that in comparison with two other models and with identical parameters, Gent model shows a better agreement with experimental data, albeit using a fairly large number of rigid links. It is reiterated that a model that requires a lower N to agree with experimental data is rated as an acceptable and successful model. In the remaining diagrams, Gent model is preferred to other models for calculating the back stress because of its finer predictions. In Figure 35, the back stress follows Gent model and the number of rigid links is fixed at 50 N = . These diagrams indicate that the locking strain is a strong function of N and almost insensitive to other parameters. On the other hand, the intercept of the curves or, strictly speaking, the effective yield stress of the aggregate, shows a strong dependence on the viscoplastic resistance of the interlamellar domain, il a , the rate exponent, n , as well as the imposed dimensionless strain rate, The pole figures for { } 002 and { } 200 planes after 100% straining are given in Figure 37. Pole figures are informative, rather qualitatively than quantitatively, vis-à-vis the trend of the texture evolution. These particular planes have preferentially been selected to trace the evolution of chain axis (c axis) and a plane with normal perpendicular to c axis. In this regard, the pole figures suggest that under tensile loading, { } It is worth noticing that the above diagrams are all output of simple tensile straining and similar diagrams associated with other modes of deformations, such as shear and compression, are required to enable a more accurate calibration of the model parameters and also a more comprehensive evaluation of the presented formulation and the proposed methodology. For a complete parameter study, as a valuable objective, one needs to run a large number of simulations with various model parameters including n ,

Conclusions and Suggestions

MECHANICAL CHARACTERIZATION of the interphase layer has been the primary objective of this dissertation, which has successfully been implemented for the first time in this work. This accomplishment has been realized by way of a mechanistic look at the atomistic data available for a polyethylene stack.

Chapter II clearly demonstrates that when experimental techniques or even atomistic simulation tools are unable to provide us with information about the elastic properties of a constituent in a heterogeneous solid, a hybrid approach can effectively be the answer to the question. The hybrid approach consisted of applying the relationships from micromechanical homogenization techniques of DIM and ECIM to Monte Carlo molecular simulation data for the interlamellar domain of polyethylene. Elementary knowledge about the simulation conditions, viz. the boundary conditions of the simulation element, has been a must for adopting the appropriate micromechanical technique. Perfect agreement between the results of DIM and ECIM indicates that the problem can be viewed from different perspectives depending on the physics of the problem at hand. To elaborate, due to the periodic boundary conditions applied to lateral planes of the simulation element, the interlamellar domain was thought of either as an ellipsoidal inclusion embedded into another hollow ellipsoidal inclusion (tailored for DIM formulation), or as a rectangular layer housed between two other rectangular layers at both sides (tailored for ECIM). Successful application of the homogenization relationships for the elastic characterization of the interphase layer shows that, in addition to the primary function of homogenization approaches, they can be used in the reverse mode for dissociation purposes too, when needed.

As an essential step to carry out the dissociation analysis, the stiffness tensor of the central amorphous phase as a function of temperature was required. Since pure amorphous samples of polyethylene cannot be prepared in today laboratories to evaluate their mechanical properties, the stiffness tensor in question was established based on the theoretical arguments for amorphous polymers as well as the relevant experimental findings. The experimental measurements for different grades of polyethylene with different crystallinities were exploited by extrapolation to zero crystallinity for estimating the amorphous elastic properties. Of the amorphous elastic constants, Young's modulus and Poisson's ratio were selected for setting up the temperature dependent stiffness tensor because their response to temperature for 350-400 K was easily deducible in light of the theoretical arguments.

In the eventual dissociated interphase stiffness, shearing components in planes normal to the interface, i.e. C , robustly took small but negative values, leading to the non-positive definiteness of the interphase stiffness tensor, at least for the temperature range of interest. We believe that this non-positive definiteness was a valid outcome whose origin lay in the fact that the interphase is a transitional domain whose existence is always accompanied by neighboring crystalline and amorphous phases that mechanically stabilize the interphase. After running a specific sensitivity analysis, we ascertained that the non-positive definiteness of the interphase stiffness was insensitive to the uncertainties of the adopted amorphous elastic constants and to the uncertainties of the interlamellar domain, for the temperature range 350-400 K. This finding was consistent with negative stresses observed at the crystalline-noncrystalline interface of semi-crystalline polymers that were reported elsewhere. There were arguments that this negativity was an indicator of the pre-strained state of the interphase layer, attributed to the way the polymer chains exit the crystals. It was also observed that, contrary to the two other shearing stiffnesses, which was shrunk to 10 MPa ± in the course of the analysis. As another finding, it has been revealed that for dissociation purposes, DIM works perfectly without posing any numerical problems while ECIM is prone to numerical divergence problems if the recursive method is employed. However, for both techniques and regardless of their level of nonlinearity, re-definition of the solution as the global minimum of one (or a set of) nonnegative objective function(s) is a reliable numerical method that will definitely lead to the solution.

Finally, using the adopted two-component sandwich model, a plausible explanation has been suggested for an empirical relationship that describes the interlamellar average Young's modulus as a function of crystallinity. In the explanation provided, the constitutive properties of the composing phases are invariant with crystallinity while only the volume fractions vary with crystallinity. On the other hand, since the interphase layer is the dominant phase in the interlamellar region at high crystallinities, the average Young's modulus of the interphase should be comparable to that of the interlamellar domain at high crystallinities. Without taking the effect of uncertainties into account, Hill's estimate of the interphase average Young's modulus at 350 K is 347 MPa. This mean value compares well with the experimentally determined 300 MPa in addition to being consistent with the established fact that the elastic modulus of a rubbery amorphous polymer increases with temperature. This good agreement serves as a yardstick for the verification of the tailored methodology and the dissociation results.

Chapter III successfully dealt with the hyperelastic characterization of the interlamellar domain and its constituents in polyethylene, again using a hybrid approach: the governing equations of a hyperelastic, compressible, isotropic, neo-Hookean solid were applied to the molecular dynamics simulation data of a polyethylene stack run at 350 K and under volume conserving mode of deformation. For carrying out the characterization scenario, we had to introduce the evolution of the amorphous and interlamellar boundaries using a series of unknown coefficients as auxiliary variables, which are identified along with the hyperelastic parameters. Solution finding was rendered into finding the global minimum of a set of nonnegative objective functions in the presence of a set of constraints.

The identified shear moduli of three domains of interest agreed well with the experimentally determined shear modulus for the noncrystalline domain. Using an ad hoc sensitivity analysis it was demonstrated that the bulk moduli were identified with a higher degree of uncertainty in comparison with the shear moduli. It was also revealed that although the lateral stress distributions disagree over a wide range of strains, their arithmetic mean was a more reliable source for identification purposes. The small size of the initial simulation stack was speculated to be responsible for the disagreement observed.

Although the deployed characterization analysis concerned the noncrystalline domain, it successfully justified the melting and recrystallization at the crystallinenoncrystalline interface as well as the rotation of crystalline chains that took place within the crystallites during the deformation. Additionally, with the implemented identification idea, there was not only no need to have any information about the evolution of the initial hyperelastic elements beforehand, but the identification of the boundaries involved was the by-product of the presented methodology. Moreover, it was also concluded that the initial part of the lateral stress diagrams were not accurate enough for the identification purposes and, as a general conclusion, the larger the straining interval is the more reliable the identification scenario is carried out. And above all, the presented analysis can be viewed as a useful complementary study to the original molecular dynamics work that can be followed for similar problems.

Chapter IV was devoted to re-considering the viscoplastic response of an aggregate of polyethylene with the emphasis on the modified Taylor approach. It was demonstrated that mathematical expression of the inextensibility constraint for chain stems in crystallites can help to propose a complete projection tensor for the modified Taylor approach, which can be followed for cases of higher constraints. Besides, among the available back stress models, although 8-chain model and Gent model are mathematically close, Gent model is favored because of its simpler mathematical representation and lower computational cost. The substantial difficulty with 8-chain model and its predecessor, 3-chain model, is calculating the inverse Langevin function, for which there is no closed-form solution to date, and approximate solutions are used instead. From a more global point of view, the constitutive law for the viscoplastic response of the noncrystalline phase needs modifications on the basis of either phenomenological aspects or theoretical arguments. As compared to the modified Taylor approach, Sachs approach was less successful. However, although Sachs and Taylor approaches are known as lower bound and upper bound approaches, they can serve as benchmarks for evaluating the success of new interaction laws. Moreover, their estimates for the yield stress can be exploited for finding Hill's estimate for the yield stress, which is expected to be close to the true one.

The common feature of presented methodologies in the foregoing technical chapters, which is more or less highlighted where appropriate, is the "system identification" perspective that was employed for implementing the methodologies: after selecting the appropriate mathematical model that suits the physics of the problem in question, the "parameter identification" or strictly speaking, the "characterization", is carried out by finding the global minimum of a single nonnegative objective function or a series of nonnegative objective functions in the presence of a set of rational constraints. This global minimum finding is realized using the appropriate optimization tool. On the other hand, the nonnegative objective functions are defined in accordance with the physics of the problem and its mathematical representation. From a more general point of view, this system identification perspective is so powerful that can be applied to a wide range of problems from scientific research sphere to engineering applications.

Suggestions for the follow-up work:

• Monte Carlo molecular simulations should be re-run at lower temperatures and with a larger number of united atoms to lower the uncertainty intervals. For the thermo-elastic analysis, all interlamellar stiffness components should be computed at various temperatures and finally a more accurate interphase stiffness will be obtained by re-carrying out the dissociation algorithm presented in Chapter II.

• A framework should be developed for following the response of heterogeneous materials containing a non-positive definite stiffness in particular after the release of the stored energy in pre-strained phases.

• Molecular dynamics simulations should be re-run for a larger initial stack, under a variety of deformation modes and up to larger strains to enable applying the governing equations of an anisotropic, compressible, hyperelastic continuum and subsequently conducting a similar hyperelastic characterization scheme.

• The viscoplastic constitutive law for the noncrystalline phase should be revised based on the phenomenological observations or theoretical arguments available for amorphous polymers.

• Appropriate molecular dynamics ensembles should be designed for identifying the parameters of new viscoplastic laws and back stress models for the noncrystalline phase, in a reliable manner.

• Necessary modifications should be made to the crystalline constitutive law for higher strain rates. Additionally, since at higher strain rates, cavities are created inside polyethylene, the interaction law should be modified accordingly to take the effects of cavitation into account.

angles between two sets of coordinate axes are uniquely determined. When the z-and Zaxis are co-directional, 0 θ = and only the sum ϕ ψ + can be determined uniquely.

Similarly, when the z-and Z-axis are in the opposite directions, θ π = and only the difference ϕ ψ -can be determined uniquely. In these two last cases, also known as gimbal lock, none of the angles , ϕ ψ are individually identifiable.

In light of the above definitions and concepts, the new representation of any vector in the rotated frame is obtained simply by pre-multiplying its vectorial representation in the fixed frame by the three following rotation matrices cos sin 0 1 0 0 cos sin 0 a1 = sin cos 0 , a2 = 0 cos sin , a3 = sin cos 0 0 0 1 0 sin cos 0 0 1

ϕ ϕ ψ ψ ϕ ϕ θ θ ψ ψ θ θ             - -             -       (A1)
The three above rotation matrices may be collected in a single rotation matrix as follows

            = - - =                   - + - - - + a sin sin sin cos cos θ θ ϕ θ ϕ θ         -   (A2)
Therefore, for the transformation of any vector from the fixed frame into the rotated one, the vector needs to be pre-multiplied by the above matrix. On the other side, for the reverse transformation, the rotated vector needs to be pre-multiplied by the inverse of the above matrix to find its representation in the fixed frame 1 cos cos sin sin cos cos cos sin sin cos sin sin cos sin cos sin cos sin sin cos cos cos sin cos sin sin cos sin cos

ψ ϕ ϕ ψ θ ψ θ ϕ ψ ϕ ϕ θ θ ψ ϕ ϕ ψ ψ ϕ ψ θ ϕ θ ϕ θ ψ ψ θ θ - - - -     = + - + -       a (A3)
Thus in the continuum mechanics notation, we have

i ij j ′ = a v v (A4)
where v represents a vector in the fixed frame and ′ v is the same vector represented in the rotated frame. Since a second order tensor is the dyadic product of two vectors, the rotation of any second order tensor, B , from one frame into another is expressed as following

ij ip jq pq ′ = B a a B (A5)
The transformation of an arbitrary fourth-order tensor, C , form one frame into another is defined likewise

ijkl ip jq km ln pqmn = ′ a a a a C C (A6)
One of the major applications of the last two relationships is in the orientational averaging of tensorial quantities frequently encountered in the area of heterogeneous materials. In many practical applications in the field of heterogeneous materials, the composing phases have a uniform size distribution and are randomly distributed in space. In such a case, instead of transforming the property tensor of individual particles to the reference frame followed by weighting them with their respective volume fraction and then adding them up, one just needs to perform an orientational averaging on the transformational expression of the property tensor. Given that the volume element in Euler angles space is

= = C C C C C C C C C C C C C C C C ′ ′ ′ ′ ′ ′ ′ ′         ′ ⋯ ⋯ ⋯ ⋯ C C (
× × × = + + + = + + +     =          =       ′ ′ C C C C A A (A9)
By carrying out the integration operation (A7) on rs A , the averaging matrix, [ ] 81 81

Avg × , is obtained. [ ] [ ] ( ) ( ) ( ) ( ) ( ) ( )

Avg or

Avg rs s r

× × × = + + + = + + +     =          =       C C C C (A10)
Performing the integration reveals that only 21 columns and 21 rows of [ ] 81 81

Avg × are nonzero which shows that when the orientational distribution is uniformly random, only 21 elements of the initial fourth-order tensor C contribute to the final averaging.

C C C C C C C C C C C C C C C C C C C C C     -                             =                                   1 1 1 1 1 1 4 2 4 1 1 1 1 1 1 4 4 2 2 4 4 1 1 1 1 1 1 4 2 4 1 1 1 1 1 1 4 4 2 2 1 1 4 1 4 1 1 4 1 2 1 4 1 1 4 1 4 1 1 2 2 1 1 1 4 1 4 4 1 1 2 1 1 4 4 1 4 1 1 1 2 2 1 1 4 1 4 1 1 4 1 2 1 4 1 1 4 1 4 1 1 2 2 1 1 1 4 1 4 4 1 1 2 1 1 4 4 1 4 1 1 1 2 2 1 1 ----- ------ ------ ------ -- - -- - - -- - -- --- - -- -- - --- -- - -- - - -- - -- --- - -- -- - --- ---1 4 1 4 4 1 1 2 1 1 4 4 1 4 1 1 1 2 2 1 1 4 1 4 1 1 4 1 2 1 4 1 1 4 1 4 1 1 2 2 4 4 1 1 1 1 1 1 4 2 4 1 1 1 1 1 1 4 4 2 6 2 2 2 2 2 2 2 2 2 6 2 2 2 2 2 2 2 2 2 6 2 4 4 1 1 1 1 1 1 4 2 4 1 1 1 1 1 1 4 4 2 2 1 1 4 1 4 1 1 4 1 2 1 4 1 1 4 1 4 1 1 2 2 1 1 1 4 1 4 4 1 1 2 1 - -- -- - --- -- - -- - - -- - -- ------ ------ ------ ------ -- - -- - - -- - -- --- - -- --1 4 4 1 4 1 1 1 2 2 1 1 1 4 1 4 4 1 1 2 1 1 4 4 1 4 1 1 1 2 2 1 1 4 1 4 1 1 4 1 2 1 4 1 1 4 1 4 1 1 2 2 1 1 1 4 1 4 4 1 1 2 1 1 4 4 1 4 1 1 1 2 2 1 1 4 1 4 1 1 4 1 2 1 4 1 1 4 1 4 1 1 2 2 4 4 1 1 1 1 1 1 4 2 4 1 1 1 1 1 1 4 4 2 2 4 4 1 1 1 1 1 1 4 2 4 1 1 - --- --- - -- -- - --- -- - -- - - -- - -- --- - -- -- - --- -- - -- - - -- - -- ------ ------ ------ -- 1111
C C C C C C C C C C C C C C C C C C C C                                                             ----       333                                                                   (A11)
Thus the rest of the elements of ijkl C missing on the left hand side of (A11) are 0.

However, the majority of fourth-order tensors encountered in practice have, at least, minor symmetry as the lowest symmetry. Mathematically speaking, the fourth-order tensor C has minor symmetry if

ijkl jikl ijlk jilk = = = C C C C (A12)
Therefore, assuming that the fourth-order tensor C in Eq.(A11) has minor symmetry, the following equalities hold

= = = = = = = = = C C C C C C C C C C C C (A13)
The immediate conclusion of the above equalities is the reduction of (A11) to the following form

C C C C C C C C C C C C     -- -     -- -   -- - - --     -- - - --   -- -   =                   
C C C C C C C C C C C C                                   -- -       -- - - --       -- -    -- -                            (A14)
Careful examination of the coefficient matrix appearing on the right hand side of the above relationship indicates that the first, seventh and twelfth rows are identical, the second, third, sixth, eighth, tenth and eleventh rows are identical and finally the fourth, fifth and ninth rows are identical. This means that the corresponding components of the abovementioned sets in the left hand side vector are equal too. In other words,

= = 2 1 1 6 6 1 2 1 6 1 1 2 C C C C C C C C C C C C C C C C C C C C C C C C          = -- -         -- - - --                                          (A15)
Upon a more careful examination of (A15) it can be seen that the rows of the coefficient matrix are not independent and the second row can be reproduced simply by subtracting twice the third row from the first row.

= = = C C C C C C C C C C C C C C C C C C C C C C                       =       -- - - --                       2233 3311 3322 1111 1212 = = = -2 C C C C (A16)
Consequently, assuming that C is the stiffness tensor of a material, then Eq.(A16) confirms that the final averaged tensor has isotropic symmetry with solely two independent coefficients ; a fact that agrees with our initial intuition.

The foregoing arguments are valid if the orientational distribution of heterogeneities is uniform. In case the distribution of the fourth-order constitutive property is not uniformly random, the integrand of Eq. (A7) is multiplied by some orientational distribution function ( ) , ,

2π π 2π im jn kp lq mnpq 0 0 0 ijkl ijkl 2π π 2π 0 0 0 , , sin = , , sin

d d d d d d φ θ ψ θ ψ θ φ φ θ ψ θ ψ θ φ Φ = Φ ∫ ∫ ∫ ∫ ∫ ∫ a a a a C C C (A17)
Equation Section (Next)

Appendix B-Eshelby problem and Eshelby tensor

J.D Eshelby is known as the first micromechanician who shaped and analytically solved the problem of the stress and strain induced in a defect confined in an infinite homogeneous medium when the defect is in the form of an ellipsoidal inhomogeneity and undergoes a free stress strain [129]. In Eshelby problem, the substantial concept of eigen strain also known as transformation strain or polarized strain is often encountered. Indeed, eigen strain is a generic designation for such free stress strains engendered by thermal expansions, initially induced deformations, phase transformations, misfittings etc. In the elasticity of eigen strains, there is no need to look for the specific origin of the eigen strain; it can be due to any fictitious source which suits the equivalent inclusion problem [130].

Let us assume that the ellipsoidal domain Λ in an elastic solid is subjected to the arbitrary deformation * ε that would be a stressless deformation if were not bounded by the surrounding body. To start with, it is assumed that Λ and the host medium have identical constitutive properties. Because of the confining effect of the surrounding medium, stress fields are created inside and outside Λ . The emphasis on the ellipsoidal shape of Λ comes from the mathematical proof established by Eshelby that demonstrates that the induced stress and strain in Λ are uniform if Λ is ellipsoidal. Additionally, the ellipsoidal morphology is comprehensive enough to cover various geometries of inclusions and reinforcements often seen in heterogeneous materials. Assuming infinitesimal deformations, the resulting strain reads

ij ij ij + or * * = + = e e ε ε ε ε (B1)
where e and ε stands for the created elastic strain and perturbed strain, respectively. The compatibility conditions for the final strain necessitates

( ) ij i,j j,i 1 + 2 = ε u u (B2)
where u is the associated displacement field vector. The created elastic strain, however, is linked to the induced stress by virtue of Hooke's law:

( ) ( ) * * ij ijkl kl ijkl kl kl ijkl k,l kl - - = = = σ e ε ε ε C C C u (B3)
where C represents the stiffness tensor of the medium. It should be noted that in the last relationship, the symmetry of the stiffness tensor is exploited. In the absence of any kind of body forces, the static equilibrium requires

ij,j 0 = σ (B4)
The free boundary conditions are mathematically expressed as

ij j 0 = σ n (B5)
where n is the unit normal to the outer boundary of the medium. Substitution from (B3) into (B4) and (B5) results in

* ijkl k,lj ijkl kl,j = ε C C u (B6) * ijkl k,l j ijkl kl j = ε C C u n n (B7)
It can be seen from Eq.(B6) that The best and most famous solution method used for the differential equation (B6) to find u is fundamental solution or the Green function technique. The Green function of the above differential equation, G, is defined as follows

( ) ( ) mnpq pk,qn mk =-δ ′ ′ G δ C x -x x -x (B8) where ( ) δ

( ) (

)

* i jlmn mn ij,l d ∞ -∞ ′ ′ ′ = -∫ ε G C u x x -x x (B9)
Re-substitution of u into (B2) gives the perturbation strain

( ) ( ) { } * ij klmn mn ik,lj jk,li 1 2 d ∞ -∞ ′ ′ ′ = - + ∫ ε ε G G C x -x x -x x (B10)
Closed-form solution for the components of the Green function tensor in the general case where the medium is anisotropic does not exist. But a relevant tensor that is more important than the above Green tensor is the Eshelby tensor and is defined as follows

e * ij ijkl kl = ε ε S (B11)
It can be shown that for the general case where the medium is anisotropic, the components of the Eshelby tensor are calculable as follows [131]

( ) ( ) { } 1 2 e ijkl mnkl 3 imjn jmin 1 0 1 8π d d π ζ ω + - = + ∫ ∫ S C H H ξ ξ ( B12) where ( ) ( ) ( ) ( ) ( )
ij ijkl k l 2 2 i i 1 3 2 3 i mnl m1 n2 l3 ij ikl jmn km ln ik ijkl j l = ; = 1-cos ; = 1-sin ; 1 = ; ξ ; = 2 D a D ζ ζ ω ζ ζ ω = = N K K K N K K K H C ξ ξ ξ ξ ξ ζ ξ ξ ξ ξ e e e (B13)
In the above relationship, ijk e is the third-order alternating tensor also known as Levi-Civita symbol or permutation symbol and i a are half-diameters of the ellipsoid.

Like Green function tensor, the components of the Eshelby tensor have no closed-form solution and in practice, computer codes are developed to calculate their numerical values [131]. For specific cases where either the symmetry of the reference medium or the geometry of the ellipsoidal inclusion is simplified, closed-form solutions for the nonzero components of the Eshelby tensor are available.

Case B-1: If the reference medium is isotropic and

1 2 3 a a a ≠ ≠ : ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )
1-2 3 = + , 8π 1- 8π 1- 1-2 1 = - , 8π 1- 8π 1- 1-2 1 = - , 8π 1-υ 8π 1- 1-2 + = + + 16π 1- 16π 1- S a I I S a I I S a I I a a S I I I υ υ υ υ υ υ υ υ υ υ υ (B14)
The remaining nonzero elements of the Eshelby tensor are calculable similarly by cyclic permutation of indices and using the minor symmetry property of the Eshelby tensor. In the above relationships, 0 υ is the Poisson's ratio of the reference medium and the integrals 1 I and 3 I are expressed as follows ( )

( ) ( ) { } ( ) ( ) 2 2 2 1 3 1 2 3 1 2 3 1 3 2 2 2 2 2 2 2 2 1 3 1 2 1 3 2 3 1 3 - 4π 4π = , - , , = - , - - - - a a a a a a a a a I F k E k I E k a a a a a a a a a a θ θ θ           (B15)
where

( ) ( ) 2 2 2 2 0 0 2 2 2 -1 3 1 2 1 2 3 2 2 2 1 1 3 , = , , = 1-sin 1-sin - =sin 1-, = , > > - dw F k E k k wdw k w a a a k a a a a a a θ θ θ θ θ ∫ ∫ (B16)
To calculate the remaining integral values, the following identities are utilized.

1 3 -1 3 = 1-2 + -1-2 + 2 1- -1 -1 1 1 9 = = + 1-2 - 8 1- -1 4 1- 4 -1 1 3 = = -1-2 + 4 1- 2 -1 4 -1 1 1 = =- - 2 1- -1 4 1- S S S S S S S τ τ υ υ ϕ υ τ τ τ υ ϕ υ τ υ τ τ υ ϕ υ τ τ τ υ τ υ                                     ( ) ( ) ( ) ( ) ( ) ( ) ( )
1 3 = =- 1-2 + -1-2 + 2 1- -1 2 -1 1 3 = = + 1-2 - 4 1- 2 -1 4 -1 3 +1 1 +1 1 = = 1-2 - -1-2 - 4 1- -1 2 -1 S S S S S S τ υ ϕ τ τ υ υ ϕ υ τ τ τ υ ϕ υ τ τ τ τ υ υ υ τ τ                                               ϕ         (B18) 1
where

3 1 = a a τ and ( ) ( ) ( ) ( ) 1 2 1 2 3 2 2 1 1 2 2 3 2 2 -1 cosh if 1 -1 cosh -1- if 1 1- τ τ τ τ τ τ ϕ τ τ τ τ τ τ - -    - >        =     <        (B19)
In two previous relationships, 3 a is assumed to lie along the 3-axis.

Case B-3: If the reference medium is isotropic and τ → ∞ (long, narrow, cylindrical

Appendix C-Double inclusion method tailored for dissociating the interphase stiffness

Double-inclusion model (as the authors in [START_REF] Hori | Double-inclusion model and overall moduli of multi-phase composites[END_REF] named it) is an attempt to present a mean-field methodology to calculate the effective stiffness of an ellipsoidal inclusion enclosing another ellipsoid, with the entire inclusion being contained in a reference medium, as illustrated in Fig C 1. The underlying assumptions together with the technical details and derivation procedure are available op. cit. Here, the final derived relationship regarding the overall stiffness of the two-phase inclusion is adopted and adapted for our specific application. According to Hori and Nemat-Nasser [START_REF] Hori | Double-inclusion model and overall moduli of multi-phase composites[END_REF] the effective stiffness of the doubleinclusion, DI C , is obtained as follows

( ) ( )( ) 1 DI ref DI DI DI DI - = + - + C C I S I H I S H (C1) where ( ) ( ) ( ) ( ) (
)

DI 1 1 1 2 1 1 1 1 1 2 1 1 1 ref ref 1 2 2 ref ref DI 1 1 1 1 1 DI η η η η η η η η η η η Ω Γ - - Ω Ω Ω Ω - - Γ Ω Ω Ω - - Ω = + -       = - + + ∆ - ∆ + - ∆ +       - -               = -∆ + + - ∆ + - ∆ +       - -         = - = - ∆ = - Ω = H H H H S E S S S E S S E H S S E S S E S S E E C C C E C C C S S S (C2) Ω Γ ref DI
In the above relations, I is the fourth order identity tensor, Ω S and R S stand for the Eshelby tensor of the inner inclusion and the entire double inclusion, respectively, ref C denotes the stiffness of the reference medium, the superscripts DI , Ω and Γ represent, respectively, the entire double inclusion, the inner inclusion, and the cover of the inner inclusion, which equals D I -Ω . 1 C and 2 C stand for the stiffness of the inner inclusion and the outer cover, respectively, and η is the volume ratio of the inner inclusion to the entire double inclusion, which is equivalent to the interlamellar region in our problem. Since in the MC simulations of the interlamellar region, periodic boundary conditions are imposed on the lateral faces, one can assume that the double-inclusion and the inner inclusion are thin, penny-shaped disks. By analogy, the interphase layers are represented by the outer cover of the double-inclusion problem, i.e.

1 1 1 1 1 1 il ref am ref ref ip ref ref am ip η η - - - - - - ∞ ∞ ∞ = + + - + + - -                     I S C C C S C C C S C C (C4) ( ) ( ) ( ) ( )
When the above relationship is solved for ip C , we obtain

( )

1 1 1 1 1 1 1 1 ip ref ref ref il ref am am ip ip 1 η η η - - - - - - - - ∞ ∞ ∞ = - - - - - - +                             I C C S I C C S S C C C (C5)
The most serious practical difficulty faced in utilizing Eqs.(C4) and (C5) is the evaluation of Eshelby tensor, ∞ S , whose computation is demanding since the reference medium is not isotropic and as a result there is no closed-form solution for its non-zero components. An efficient numerical algorithm for computing the Esheby tensor of any ellipsoidal inclusion embedded in a reference medium of the general anisotropic symmetry can be found in [131].

Before wrapping up this appendix, it should be noted that in the last two relationships, ip C and il C depend on the stiffness of the selected reference medium. Interestingly, we observed that in the case of the aligned, disk-like inclusions, the effective stiffness of the double-inclusion is independent of the choice of the reference stiffness. Mathematical demonstration of this observation looks very difficult due to the nonlinear dependence of Ω S on ref C , on one hand, and the nonlinear dependence of il C on both, on the other hand. Therefore, any arbitrary ref C , can be substituted in the above relationships and the result will be independent of it. Although, the mathematical demonstration of the above statement for the general case of anisotropic ref C looks intimidating, for the simpler case of isotropic reference medium the demonstration is easier and rather straightforward. For the isotropic reference medium, ref C and ∞ S are represented as follows in the 6×6 matrix format:

( )( ) ( ) ( ) ( ) ref ref ref ref ref ref ref ref ref ref ref ref ref ref ref ref 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 1 2 2 0 0 1 1 2 0 0 0 0 1 2 2 0 0 0 0 0 0 1 2 2 E υ υ υ υ υ υ υ υ υ υ υ υ υ υ -     -     - =   - + -     -   -     C (C6) ref ref ref ref 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 υ υ υ υ ∞           - - =             S (C7)
where ,

C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C         =         C ip ip ip
C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C      =           C          (C8)
To show the independence of il C from the choice of isotropic ref C , one needs to substitute the above values into relationship (C4) and simplify the result. In the final simplified relationship, there must be no trace of either ref E or ref υ . The mathematical operations involved cannot be worked out manually and symbolic calculation software are indispensable. Here, Maple® has been chosen as the symbolic calculation aid. However, the abovementioned substitutions and the subsequent symbolic calculations are still too heavy for Maple®. Therefore, prior to substitution, relationship (C4) is reformulated as follows:

( ) ( ) ( ) ( ) ( ) ( ) 1 1
1 1

1 1 il ref 1 ref ref ref 1 1 1 1 ref am ref ref ip ref am ip where η η - - - - - - - ∞ ∞ - - - - ∞ ∞ = + - + -     =         = + - + + - C U S C C U I S C U C U S C C C S C C C (C9)
This new form of representation is simpler to handle for symbolic calculation software.

The following short code has been written in Maple® environment that takes care the above operations.

Due to the volume of calculations involved in simplifying il C , each component of il C is simplified separately. The simplified expressions for each component of il C are so long that take several pages, and this is the reason they are not given here. However, in the simplified expressions for each component of Equation Section (Next)

Appendix D -Extended Composite Inclusion Model

The original version of the composite inclusion model (CIM) was an effort to attain a more realistic estimation of the effective stiffness of a two-layer composite inclusion than those suggested by the Voigt and Reuss models [START_REF] Ahzi | Estimates of the overall elastic properties in semicrystalline polymers[END_REF][START_REF] Ahzi | Composite modeling for the effective elastic properties of semicrystalline polymers[END_REF]. In light of the notion presented therein, this approach is extended to calculate the effective stiffness/compliance of a three-layer composite inclusion. The assumptions made here are exactly similar to those discussed in the above-cited references: average theorems are implied; linear elasticity for each phase as well as the entire composite inclusion is assumed; the stress/strain in each phase is correlated with the stress/strain of the entire inclusion through initially unknown weight functions / R Q ; and finally, the uniformity of stress and strain in each phase is assumed.

The schematic of a three-layer composite inclusion is depicted in Fig D 1. From the average theorems, the average stress and strain of this composite inclusion, I σ and I ε , are expressed as

I 1 2 3 1 2 3 η η η = + + σ σ σ σ (D1) I 1 2 3 1 2 3 η η η = + + ε ε ε ε (D2)
where i η , i σ and i ε stand for the volume fraction, average stress and average strain of the ith layer/phase, respectively, with i 1,2,3 = . Let us assume that the governing linear elastic constitutive law for each phase follows as ( )

1 i i i i i i i i
or with

- = = = ε σ σ ε C S C S (D3)
where i S and i C are, respectively, the compliance and stiffness of the ith phase. If the effective compliance and stiffness of the three-layer inclusion are defined as coefficients correlating I σ and I ε as follows:

( )

1 I I I I I I I I
or with

- = = = ε σ σ ε C S C S (D4)
then substitution of (D3) and (D4) into (D1) and (D2) yields

I 1 1 2 2 3 3 1 2 3 η η η = + + C C Q C Q C Q (D5) I 1 1 2 2 3 3 1 2 3 η η η = + + S S R S R S R (D6)
where the weight functions i Q and i R , also called strain and stress concentrations, respectively, are defined as Therefore, once the stress or strain concentrations are determined, the effective stiffness/compliance of the composite inclusion is calculable. An important auxiliary assumption of CIM that is not explicitly stated in [START_REF] Ahzi | Estimates of the overall elastic properties in semicrystalline polymers[END_REF][START_REF] Ahzi | Composite modeling for the effective elastic properties of semicrystalline polymers[END_REF] but is invoked implicitly is that the stress and strain in each phase are assumed to be uniform. This assumption, in conjunction with the enforcement of the equilibrium conditions at the two interfaces, gives rise to

i i I = ε ε Q (D7) i i I = σ σ R (D8)
1 2 3 I β β β β = = = σ σ σ σ (D9)
where β 3, 4,5 = . In Eq.(D9) the Voigt contracted notation is used. Additionally, in agreement with the convention adopted in [START_REF] Veld | Temperature-dependent thermal and elastic properties of the interlamellar phase of semicrystalline polyethylene by molecular simulation[END_REF][START_REF] Ahzi | Estimates of the overall elastic properties in semicrystalline polymers[END_REF][START_REF] Ahzi | Composite modeling for the effective elastic properties of semicrystalline polymers[END_REF], the interfaces are parallel to the 12-plane, consequently the 3-axis is normal to the interfaces (see Fig D 1). Substitution of (D9) into (D8) leads to

i i I I i β βj j β βj βj = = → = σ σ σ δ R R (D10)
where δ denotes the extended Kronecker delta. Compatibility conditions at the interfaces along with the assumption of uniformity of strains in each phase require that if α 1, 2, 6

= then 1 2 3 I α α α α = = = ε ε ε ε (D11)
or equally Further decomposition of (D13) gives

( ) ( ) ( ) (
) ( ) ( )

1 1 1 1 2 1 2 1 2 3 2 3 2 I 2 αα α αβ β 1 αα α αβ β 3 αα α αβ β αj j 1 2 1 1 2 1 2 3 2 3 2 I 2 αα 1 αα α 2 αβ 1 αβ β 3 αα α 3 αβ β αj j η η η η η η η η η ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ + + + + + = + + + + + = σ σ
1 1 1 1 3 3 3 3 αα α αβ β αα α αβ β 1 3 3 1 1 1 3 1 α αα αα α αβ αβ β ′ ′ ′ ′ ′ ′ ′ ′ - ′ ′ ′ ′ ′ ′ ′ + = + = + - σ σ σ σ σ σ σ S S S S S S S S (D15)
Simultaneous use of Eqs.(D14) and (D15) yields

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1
η η η η η η η η η η η η η η ′ ′ ′ ′ ′ ′ ′ - - ′ ′ ′ ′ ′ ′ ′ ′ ′ - - ′ ′ ′ ′ ′ ′ ′ ′ ′
( )

1 1 1 2 2 3 1 1 1 2 2 2 3 1 3 1 2 I 2 αα 1 αα 3 αα αα αα α 2 αβ 1 αβ 3 αβ 3 αα αα αβ αβ β αj j 1 1 1 1 1 2 2 3 1 2 1 2 2 2 3 1 3 2 αα 1 αα 3 αα αα αα αj 2 αβ 1 αβ 3 αβ 3 αα αα αβ αβ β α η η η η η η η η η η η η η η - - ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ - - - ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ + + + + + + - = = + + - + + + - σ σ σ σ δ S S
η η η η η η η - - - ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ = + + - + + + -       δ S S S S S S S S S S S S S R (D20)
The elements of the six-by-six stress concentrations are now determined and can be substituted in Eq.(D6) for the calculation of the effective compliance. A similar procedure can be followed for the derivation of the strain concentrations, leading to

i αj αj = δ Q (D21) ( ) ( ) ( ) ( ) 1 1 1 1 2 1 2 3 1 2 1 2 2 3 1 3 βj 1 ββ 2 ββ 3 ββ ββ ββ ββ βj 2 βα βα 3 ββ ββ βα βα αj η η η η η - - - ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ = + + - - + -               δ δ C C C C C C C C C C C C Q (D22) ( ) ( ) ( ) ( ) 1
η η η η η - - - ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ = + + - - + -               δ δ C C C C C C C C C C C C Q (D23) ( ) ( ) ( ) ( ) 1 1 1 3 2 3 2 1 3 2 3 2 2 1 3 1 βj 3 ββ 2 ββ 1 ββ ββ ββ ββ βj 2 βα βα 1 ββ ββ βα βα αj η η η η η - - - ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ ′ = + + - - + -               δ δ C C C C C C C C C C C C Q (D24)
It is worth mentioning that Gueguen et al.

[132] made an attempt to derive similar relationships for the stress and strain concentrations but made errors and ended up with erroneous relationships.

In our example of the interlamellar region, the properties of phase 1 and phase 3 are identical as they represent the side interphase layers. Therefore, Eq. (D6) is rewritten as follows:

( ) ( ) ( )

1 1 1 I I 1 1 2 2 1 2 2η η - - - = = + C C C S R R (D25)
In the above equation, the unknown is Rearrangement of (D25) gives

( ) ( ) ( ) ( ) ( ) 1 1 1 1 1 I 2 2 1 2 1 1 1 ip ip il am am am am 2 or 1 η η η η - - - - - - = - = - -         C C R C C R C R C R (D26)
Similar treatment of Eq.(D5) leads to the following relationship ( ) ( )

1 ip il am am ip am am 1 η η -   = - -   C C C Q Q (D27)
Obviously, Eqs. (D26) and (D27) have the familiar form of

( ) = X X f (D28)
where f is a non-linear tensorial function of X . Apart from the arguments around the existence and uniqueness of the solution for Eq.(D26) or (D27), which are beyond the scope of this survey, the very first solution method which looks to suit the equation at hand is the numerical recursive method. In other words, some initial ip initial C is inserted into the right hand side of Eq. (D26) or (D27) whose result is fed into itself as many times as required until a certain convergence criterion is satisfied. For the problem examined in this study, it was observed that when il C is of orthotropic symmetry, Eq.(D28) converges to the solution using the recursive method while convergence problems are very likely when il C has monoclinic symmetry. In the latter case, optimization techniques are more promising alternatives. To this end, a nonnegative cost function is defined such that its optimization parameters correspond to the elements of X and its objective is minimizing the norm of ( )

S                             (E1)
It is underlined that the above matrix representation is quite subjective. The above representation is easy to reproduce and serves as a useful means to our objective which is finding the inverse of a fourth-order tensor. Another advantage of the above representation is that instead of tensorial product of two fourth-order tensors, their matrix forms can be multiplied like two ordinary matrices.

By definition, if there exists some fourth-order tensor ijkl J such that ijkl klmn ijkl klmn ijmn

= = = = SJ S J JS J S I (E2)
then J is the inverse of S and vice versa. In the above relationship, ijmn I stands for the fourth-order identity tensor and is defined as

( ) ijmn im jn jm in 1 = δ δ +δ δ 2 I
. Matrix expansion of the above relationship is as follows

J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J S ×                           
J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J
J J J J J J J                            
1 0 0 0 0 0 0 0 0 0 0.5 0 0.5 0 0 0 0 0 0 0 0.5 0 0 0 0.5 0 0 0 0.5 0 0.5 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0.5 0 0.5 0 0 0 0.5 0 0 0 0.5 0 0 0 0 0 0 0 0.5 0 0.5 0

0 0 0 0 0 0 0 0 1               =               (E3)
If the initial fourth-order tensor S has no symmetry and its rows and columns are independent then pre-multiplication of the right-hand side of (E3) by the inverse of the matrix form of S yields the matrix form of 1 -S . However, most fourth-order tensors in micromechanics have at least minor symmetry which means that the second row/column is identical with the fourth row/column, the third row/column is identical with the seventh row/column, and finally the sixth row/column is identical with the eighth row/column; consequently the matrix form of S is not invertible. Nevertheless, the property of minor symmetry of S can be exploited to eliminate the redundant elements and reduce the dimension of matrix form of S from 9×9 to 6×6.

Minor symmetry property of S implies that ijkl jikl ijlk jilk = = = S S S S . But the fourth- order identity tensor I has the same property of minor symmetry which entails that J has minor symmetry as well, i.e. ijkl jikl ijlk jilk = = = J J J J . Therefore, eliminating the duplicate rows of the matrix forms of S and duplicate columns of the matrix form of J , simplifies the relationship (E3) to the following simpler form

J J J J J J J J J J J J J J J J J J J J J S S S S S S S S S S S S S S                           
J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J                            
1 0 0 0 0 0 0 0.5 0 0 0 0 0 0 1 0 0 0 0 0 0 0.5 0 0 0 0 0 0 0.5 0

0 0 0 0 0 1               =               (
J J J J J J J J J J J J J J J J J J J J J J S S S S S S S S S S S S S S         ×             3
J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J                            
1 0 0 0 0 0 0 0.5 0 0 0 0 0 0 1 0 0 0 0 0 0 0.5 0 0 0 0 0 0 0.5 0

0 0 0 0 0 1         =           (E5)
But the duplicate columns in the matrix form of S and duplicate rows in the matrix form of J have not been eliminated. To have a clearer image of the way these two features are exploitable, the product of the second row of the matrix form of S and the first column of the matrix form of J is expanded as follows

J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J         ×           
1 0 0 0 0 0 0 0.5 0 0 0 0 0 0 1 0 0 0 0 0 0 0.5 0 0 0 0 0 0 0.5 0 0 0 0 0 0 1

                           =           (
J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J         ×           
1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0.5 0 0 0 0 0 0 0.5 0 0 0 0 0 0 0.5

                           =           (E8)
For easier recall, the previous relationship can be rewritten in the following form 4

J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J S J J J J J         ×             213 1212 J         =            
1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1

        =           (
J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J J        ×            1223 1213 2 2 J J J         =            
1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1

        =           (E10)
Given the modification made on the 9×9 matrix forms of S and J which have minor symmetry, the reduced forms of 6×6 matrices have independent rows and columns meaning that they are invertible. Therefore, once the inverse of any fourthorder tensor that has minor symmetry is sought, one needs to form either of the square matrices on the left hand side of (E10) as demonstrated and calculate its inverse. Using the calculated inverse and taking the coefficient 2 of some elements into account, 36 independent elements of the looked for inverse tensor become readily available.

Equation Section (Next)

Appendix F-Voigt notation and the 6×6 matrix form of stiffness and compliance tensors

As demonstrated in Appendix C, the fourth-order tensors having minor symmetry can be represented as 6×6 square matrices after eliminating their repeated elements. According to Voigt notation convention, the following indexing system can be used as a more convenient alternative for representing the entries of a second-order tensor, such as stress or strain tensors, or fourth-order tensors, such as stiffness or compliance tensors.

→ → → → → → (F1)
The following example of the Hooke's law explains how to use this convention. If the components of the stiffness tensor are collected in a 9×9 square matrix as demonstrated in Appendix C then the stress and strain components take the following sorting order in a column vector.

C C C C C C C C C C C C C C C C C C C C C C C C C C C C C σ σ σ σ σ σ σ σ σ               =              
C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C
C C C C C C C C C C C C ε ε ε ε ε ε ε ε ε                                                         (F2)
Similar to the arguments given in Appendix C, the symmetry of the stress tensor together with the symmetry of the stiffness allows for simplifying the above equation to the following form

C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C σ σ σ σ σ σ         =             2223
C C C C C C C C C C C C C C C C C C C C C ε ε ε ε ε ε ε ε ε                                                 (F3)
The following expansion helps to better understand the following simplification of Eq.(F3). Expansion of the second row in (F3) yields

C C C C C C C C C C C C C C C σ ε ε ε ε ε ε ε ε ε ε ε ε ε ε ε = + + + + + + + + = + + + + + (F4)
The coefficients "2" in the above expansion appear due to the symmetry of the strain tensor as well as the minor symmetry of the stiffness tensor. Similar expansions lead to the following simpler form of representation

C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C σ σ σ σ σ σ         =            
ε ε ε ε ε ε                                         (F5)
Although the coefficients "2" in the strain vector can be assigned to the corresponding elements of the stiffness matrix, in the continuum mechanics literature they appear in the strain vector. Performing the Voigt notation for the relationship (F5), we have

C C C C C C C C C C C σ ε σ ε σ ε σ ε σ ε σ ε σ ε σ ε σ ε σ ε σ ε σ ε                                 = =                                                 11
C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C         =             41
C C C C C C C C C C C C C C C C C C                     (F6)
In other words, Hooke's law in Voigt notation takes the following widely known representation

C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C σ ε σ ε σ ε σ ε σ ε σ ε                         =                                     (F7)
Similarly, it can be demonstrated that the other form of Hooke's law, which is expressed using the compliance tensor, takes the following representation

σ ε σ ε σ ε σ ε σ ε σ ε σ ε σ ε σ ε σ ε σ ε σ ε                                 = =                                                 11
                   (F8)
This last equivalence is evident in light of the discussions of Appendix C and in particular from relationship (E9) because the stiffness and compliance tensors are each other's inverse in the matrix form as they are in the tensorial form.

In the literature, the stiffness and compliance tensors are often reported in 6×6 matrix format instead of reporting all of their 81 elements. Besides, when 6×6 matrix forms of compliance and stiffness tensors of a material are expressed in their material reference frame, the type of available symmetry can be visually recognized. To establish a direct link between the stiffness/compliance of some material and its elastic constants the following simple algorithm can be used. For a fully anisotropic material, the 6×6 matrix form of the compliance tensor takes the following form as a function of the elastic constants

1 - - 1 - - 1 - - 1 1 g g g E E E E E E g g g E E E E E E g g g E E E E E E g g g E E E G G G g g g E E E G G G g υ υ υ υ υ υ κ κ κ κ = S 1 12,2 12
1 g g E E E G G G κ κ                                       (F9) E α , αβ υ , G αβ , , g αβ γ ,
, αβ γζ κ stand for elastic muduli, Poisson's ratios, shear moduli, shearextension coupling coefficients, and Chentsov coefficients, respectively [133]. Therefore, once such coefficients are available for a material, its compliance matrix is easily constructed from which its stiffness tensor can be readily calculated. On the contrary, if the stiffness matrix of a solid in its material reference frame is known, from its inverse and using the definition (F9), the associated elastic constants are determined.

It is worth noting that the major symmetry property of stiffness and compliance tensors is manifested in the symmetry of their 6×6 matrix formats. This means that the stiffness and compliance matrix of any material is necessarily symmetric irrespective of the degree of its material symmetry and this is equivalent to saying that every stiffness or compliance tensor possesses at most 21 independent components.

Equation Section (Next)

Appendix G-Crystalline structures and material symmetries

As discussed is Appendix D, the stiffness tensor of any material, as its constitutive property, has at most 21 nonzero independent components. Depending on the existing material symmetry, the number of independent stiffness elements may be less than 21. In linear elasticity of anisotropic solids, the arguments about material symmetry have been shaped based on the crystallographic considerations. In crystallography there exist 14 Bravis lattices which fall into seven lattice systems. On the other hand, it has been demonstrated that linear elastic symmetries do not exceed eight symmetries [134]. Of the available symmetries, if the isotropic one is excluded, each lattice system corresponds a unique linear elastic symmetry. Strictly speaking, there is a one-to-one correspondence between the set of lattice systems and the existing material symmetries in linear elasticity. In triclinic lattice system, the three primitive vectors are of unequal length and they are not orthogonal to one another. This system has the lowest symmetry among 3D crystal systems. In linear elasticity, the lowest symmetry takes place when there is no plane of symmetry which means that independent of the selected reference frame, the stiffness tensor has 21 independent components. Since there is no preferred reference frame, it can undergo arbitrary orthogonal rotations. Therefore, one can find three appropriate Euler angles such that in the rotated reference frame, three out of 21 independent stiffness components becomes null. Subtracting the three Euler angles equivalent to the three nulled components, the stiffness tensor of a triclinic material will have 18 nonzero, independent components. The next crystal system is the monoclinic lattice system where the basis vectors are of difference length and the crystal unit looks like a rectangular prism with parallelogram base. Corresponding to the monoclinic lattice system, the next lowest symmetry is monoclinic symmetry where there only exists a single material symmetry plane. In this case, it can be shown that not more than 13 stiffness components are nonzero if one of the planes of the reference frame coincide the symmetry plane. But again, the reference frame is indifferent with respect to the rotation about the normal to the symmetry plane. Therefore, one can find an appropriate rotation angle for the new reference frame where one of the 13 independent components becomes null and the independent, nonzero components of the stiffness tensor reduce to 12. It can be shown that if 12-plane of the reference frame coincide the symmetry plane, then the stiffness tensor take the following general form of matrix representation

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 C C C C C C C C C C C C C C C C C C C C                     (G1)
On the other hand, if the symmetry plane is 23-plane, then the stiffness matrix will assume the following configuration

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 C C C C C C C C C C C C C C C C C C C C                     (G2)
And finally, if the symmetry plane is 13-plane of the reference frame, then the stiffness matrix takes the following representation

C C C C C C C C C C C C C C C C C C C C                     (G3)
When the three basis vectors have nonidentical lengths but are mutually perpendicular, the crystal lattice is called orthorhombic and the associated linear elastic symmetry is orthotropic. In such a lattice system, there exist two symmetry planes that are perpendicular to one another which results in a third symmetry plane that is perpendicular to the first two planes. It can be demonstrated that the stiffness tensor expressed in the reference frame made by the normals of the symmetry planes, has the following general matrix form

C C C C C C C C C C C C                     (G4)
The next lattice system is tetragonal or square cuboid where two of the basis vectors are equi-length but the third vector has a different length. Similar to the orthorhombic lattice system, the vectors are mutually orthogonal. Tin, zirconium and magnesium chloride are examples of crystalline materials with tetragonal crystal system. The material symmetry in such solids is tetragonal where the normals to four symmetry planes lie in the fifth symmetry plane and this last symmetry plane is divided into eight equal sectors by the first four symmetry planes. If the 1-and 2-axis of the reference frame are each co-directional with one of the normals to the first four symmetry planes, then the stiffness tensor takes the following matrix form

C C C C C C C C C                     (G5)
The next lattice system is trigonal where the three primitive vectors have identical lengths and make equal, non-orthogonal angles with one another. The associated material symmetry is trigonal: the normals to the three symmetry planes are coplanar and lie at equal angles of 120° relative to each other. If the 3-axis of the reference frame is normal to the plane containing the normals to the symmetry planes and one of the other axes of the reference frame is codirectional with one of the normals of the symmetry planes, then the stiffness tensor takes the following matrix form ( )

0 0 0 - 0 0 0 0 0 - Sym. 0 1 - 2 C C C C C C C C C C C C C                       (G6)
Hexagonal crystal system is the six crystal system whose crystalline unit looks like a prism with regular hexagonal base. In this crystal system, two of the basis vectors have the same length making an angle of 120° with each other and the third vector is normal to the first two. Zinc, magnesium, cobalt, cadmium and alumina are examples of solids with hexagonal crystal system. The corresponding material symmetry is hexagonal which is sometimes referred to as plane isotropy. In this symmetry, there exist an infinite number of symmetry planes such that one of them contains the normals to the rest of the symmetry planes. If the 12-plane of the reference frame coincide the unique symmetry plane that is orthogonal to all other planes, then the stiffness tensor takes the following matrix form

0 0 0 0 Sym. 0 1 - 2 C C C C C C C C C C                       (G7)
Finally, the last crystal system is cubic which is the simplest one. In the cubic crystal system, all the basis vectors are of equal length and are orthogonal to each other. Aluminum, copper, silver and iron are examples of solids with cubic crystal system. The material symmetry associated with such materials is cubic symmetry with nine symmetry planes. If three of the normals to these symmetry planes are co-directional with 1-, 2-and 3-axis of the reference frame, then the other normals are aligned with the diagonals of the planes orthogonal to the first three normals. In this case the stiffness tensor has the following matrix form

C C C C C C C C C                     (G8)
For isotropic symmetry which is the simplest material symmetry, there is no corresponding crystal system. Every plane in the space is the plane of symmetry and the stiffness tensor in any reference frame has two independent constants and takes the following general matrix form of representation

0 0 0 0 0 0 0 0 0 1 - 0 0 2 1 Sym. - 0 2 1 - 2 C C C C C C C C C C C C                           (G9)
For isotropic symmetry, there is no preferred reference frame and the representation of the stiffness matrix is independent of the adopted reference frame. This is equivalent to saying that the stiffness tensor is insensitive to the rotation of the reference frame and to the orientational averaging operation.

Equation Section (Next)

Appendix H-Positive definiteness of stiffness and the theorems of upper bound and lower bound

Using the first law of thermodynamics, it can be shown that in an infinitesimal, linear elastic, adiabatic deformation, the strain energy density function, 0 W , is expressed as follows 0 1 2

p pq q W = ε ε C (H1)
Strain vector ε and stiffness tensor C are defined according to the Voigt notation convention discussed in Appendix D and summation over dummy indices is implied. Mechanical stability necessitates that when the material is unloaded from any linearly elastic loading, it retrieves its initially unloaded configuration. This understanding is equivalent to saying that the minimum state of energy of any material is its unloaded state and the strain energy function of any material is always positive. Therefore, if the strain vector is nonzero, the following inequality always holds

q 0 p pq > ε ε C (H2)
Additionally, from linear algebra we know that the real, square matrix n n D × is positive definite if and only if for any nonzero vector 1 n× a the following inequality is always true

T 0 D > a a (H3)
where the superscript T denotes the transpose. Comparison of the two recent inequalities suggest that the stiffness matrix of any material that can be found in isolation in nature is positive definite. Again from algebra, the inverse of any positive definite matrix is positive definite which leads to the positive definiteness of the compliance matrix. Additionally, in linear algebra it is demonstrated that the eigen values of any symmetric, positive definite matrix are all positive. For compliance and stiffness matrices, this last proposition is equivalent to the positivity of all diagonal elements in any reference frame. In the following paragraphs, the positive definiteness of the stiffness and compliance matrices are exploited to demonstrate that the overall stiffnesses estimated by the Voigt and Reuss techniques produce upper bound and lower bound of the effective stiffness of any nonhomogeneous material, respectively.

Assuming that the stress and strain induced in a nonhomogeneous body subjected to external loading X are σ and ε , respectively, if n denotes the outward normal to the external surface of the medium, ϒ , then

j ij i = on ϒ σ n Χ (H4)
Equation Section (Next)

Appendix I-Constitutive law of an isotropic, homogeneous, hyperelastic solid

A hyperelastic material is, by definition, a solid whose stress-strain relationship is derivable from an elastic potential energy that is necessarily a strain energy density function. Hyperelastic or Green elastic materials, which are typically characterized by low stiffness and large yield strain, are a subgroup of Cauchy elastic materials. The stress-strain response in hyperelastic materials is usually nonlinear and isotropic.

The Clausius-Duhem inequality for an isothermal process necessitates that in a hyperelastic material, the second Piola-Kirchhoff stress tensor, Y , has the following relationship with the strain energy density function of the material, W , and the deformation gradient tensor, F :

1 W -∂ = ∂ Y F F (I1)
In terms of the Cauchy stress tensor, the constitutive law of a homogeneous, compressible, hyperelastic material reads

1 t W J ∂ = ∂ σ F F (I2)
where J is the Jacobian of deformation. On the other hand, it can be shown that in an isotropically hyperelastic material, the strain energy density function can be uniquely expressed as a symmetric function of the principal stretches or equivalently as a function of the principal invariants of the left Cauchy-Green deformation tensor, t B = FF [137]. Mathematically speaking, ( ) ( ) ( ) By using the chain rule, the constitutive relationship for an isotropic, homogeneous, compressible, hyperelastic material takes the following mathematical representation

1 2 3 1 2 3 2 2 2 2 2 2 2 2 2 2 2 2 2 1 1 2 3 2 1 2 2 3 3 1 3 1 2 3 ˆ, , , , with: ,
2 1 3 1 2 2 3 3 2 2 W W W W I I I I I I I     ∂ ∂ ∂ ∂ = + - +     ∂ ∂ ∂ ∂     σ B B 1 ɶ ɶ ɶ ɶ (I4)
A variety of empirical and phenomenological strain energy functions have to date been proposed. Valuable, comparative reviews on the commonly used models are available in [121,[138][139][140].

For pressure insensitive, hyperelastic materials, a constant pressure term, p, is added to the right hand side of the relationship for the Cauchy stress to make it deviatoric. In other words, the normal Cauchy stress components in the principal axes is expressed as follows

i i i W p σ λ λ ∂ = - ∂ (I5)
Since the viscoplastic deformation of the noncrystalline domain in polyethylene is assumed to be pressure insensitive, the back stress term that appears in its constitutive law must be pressure insensitive, or equally deviatoric, as well. According to the relationships of Subsection 4.2.3 from Chapter IV, for 3-chain back stress model, we will have

3ch 1 1 1 1 3 1 2 1 2 3 1 1 1 3 1 2 1 2 3 1 or 3 1 1 1 0 3 3 3 1 9 i i r i r r r r C N p N C N p C N p C N p N N N p C N N N N λ σ λ λ λ λ λ λ λ λ λ λ λ λ λ - - - - - - -   = -           -+ -+ -=                     ⇒ = + +                 L L L L L L L (I6)
Therefore, the normal stress of deviatoric, Cauchy stress tensor in principal axes as a function of principal stretches takes the following representation according to 3-chain model:

3ch 1 1 1 1 3 1 2 1 2 3 1 3 3 i i r i C N N N N N λ λ λ λ σ λ λ λ λ - - - -             = - + +                         L L L L (I7)
Similarly for 8-chain back stress model one can write

2 chain 2 8ch 1 chain r chain 2 2 2 1 1 1 chain chain 3 chain 1 2 r r r chain chain chain 1 2 2 2 chain r 1 2 3 chain 3 1 or 3 1 1 1 0 3 3 3 1 9 i i N C p N N N N C p C p C p N N N N p C N λ λ λ σ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ - - - - -   = -           -+ -+ -=                 = + + ⇒       L L L L L 1 2 chain r chain chain 1 3 N p C N λ λ λ -  =     L (I8)
Therefore, the diagonal components of the deviatoric, Cauchy stress in principal axes and from the 8-chain model have the following form of representation ( )

2 2 chain 8ch 1 chain chain 1 3 i i r N C N λ λ λ σ λ - -   =     L (I9)
Similarly for Gent back stress model, we can write

( ) ( ) ( ) Gent 2 2 2 m r r r 2 2 m 1 chain chain 2 2 2 r 1 r 2 r 3 2 2 2 chain chain chain 2 r chain 2 chain 3 1 1 3 1 3 1 1 1 1 0 1 i i i i N J N C p C p C p J J N N N N N C p C p C p N N N N p C N σ λ λ λ λ λ λ λ λ λ λ λ λ λ - - = -= -= - - - -- - - - - -+ -+ -= ⇒ - - - - = - (I10)
Therefore, for the normal components of deviatoric, Cauchy stress tensor from Gent back stress model, we obtain ( )

Gent 2 2 r chain 2 chain 1 i i N C N σ λ λ λ - = - - (
- - = - - = - - =     = -             = -             = -         ∑ ∑ ∑ L L L L L L 9800 1 1 k 0 Exact This work ~74 10000 10000 i i i - - =     = -         ∑ L L (I16)
The interested reader may consult reference [143] for discussions about how to build the Taylor series expansion for ( )

1 y - L
although the use of Taylor series expansion is discouraged because it is not only numerically expensive but also prone to creating instability problems.

And finally, the arguments about the equivalence of 8-chain model and Gent model can be easily examined by equating the right hand sides of Eqs.(I9) and (I10). This assumption leads to The above relationship suggests that ( )

1 chain 2 chain chain 1 1 3 N N N N λ λ λ - -   =   -   L (I17)
1 y - L
is also a function of N , which does not appear as its argument! Because such a conclusion is not logically admissible, therefore the initial assumption is false and 8-chain model is not equivalent to Gent model. However, the disagreement between them is not considerable. This argument can be verified by replacing

i i i i N C C N N C N N N λ λ λ σ λ λ λ λ λ λ λ λ λ -  - -     = - ≈ - = - - - L (I19)
It can be seen that except for the numerator of the fractional expressions in Eqs.(I10) and (I19), the two relationships are identical. Therefore, the disagreement between the results of the two models is expected to be very small.

1-Introduction

Le polyéthylène est le polymère thermoplastique le plus répandu, le plus connu et le plus utilisé dans le monde. Les sacs plastiques, les bouteilles de shampoing, les jouets d'enfants, divers composants des systèmes d'irrigation sont des exemples bien connus tandis que les gilets pareballes [1], les arthroplasties complètes de la hanche et du genou sont des exemples moins connus des applications du polyéthylène [2,3]. Les Propriétés intéressantes du polyéthylène ont conduit à une large gamme d'applications. La bonne résistance à des chargements ainsi que la grande inertie chimique offertes, le tout à un prix économique font parties de ses meilleurs avantages.

Les applications potentielles de diverses nuances de polyéthylène dépendent de sa microstructure sous-jacente. Tous les types de polyéthylène sont, d'un point de vue microstructural, partiellement cristallins et partiellement non cristallins, d'où le qualificatif de « semi-cristallin ». Cependant, la partie non cristalline n'est pas monolithique et est constituée de deux composants : la phase amorphe centrale et la couche d'interphase qui rattache cette dernière à la lamelle cristalline. D'un point de vue théorique, Flory [4] a montré que la continuité des chaînes polymériques suggère le développement d'un nouvel ordre intermédiaire entre celui des cristallites et celui de la phase amorphe centrale. Il suggéra que l'ordre existant dans les cristallites ne pouvait pas disparaître subitement à l'interface des phases cristalline -non cristalline car la transition brutale d'une région ordonnée à une région désordonnée est impossible.

La présence de la couche d'interphase en tant que composant du polyéthylène a aussi été expérimentalement démontrée. Une large gamme de techniques expérimentales, y compris la RMN protonique [5], la RMN 13 C de haute résolution [6,7], la microscopie électronique [8], la diffusion de neutrons à petit angle [9], la relaxation diélectrique [10,11] et la diffusion Raman [12,13] confirment la contribution volumique importante de la région d'interphase dans les polymères semi cristallins. Mandelkern [14] maintient que la différence entre les cristallinités mesurées dans l'expérience de la fusion de chaleur et celle de densité est directement liée au contenu de l'interphase.

La plupart des études dédiées aux propriétés de l'interphase se limitent uniquement aux caractéristiques géométriques et à la fraction volumique. Bien que l'interphase soit plus raide que la phase amorphe centrale, elle a bien échappé à toutes les tentatives de caractérisation mécanique. Concrètement, les propriétés mécaniques des phases cristalline et amorphe ont été identifiées par divers moyens au cours des dernières décennies, mais celles de l'interphase sont restées inconnues essentiellement à cause de son épaisseur nanométrique ainsi qu'en raison de son ancrage étroit avec les phases adjacentes. La tâche de la caractérisation mécanique de l'interphase a été entreprise pour la première fois dans ce travail. Dans ce but, nous avons développé deux approches micromécaniques adaptées à la problématique de la caractérisation élastique de cette interphase. Ensuite, ces relations sont appliquées aux données de la simulation Monte Carlo moléculaire sur le tenseur de rigidité pour la zone interlamellaire. Le tenseur de rigidité de la phase amorphe centrale est alors requis. Il est établi à partir des résultats expérimentaux et fondés sur des arguments théoriques. Dans la partie suivante, la caractérisation hyperélastique de la zone interlamellaire et ses composants, c.à.d. la phase amorphe centrale et l'interphase, est mise en oeuvre en appliquant les équations d'un milieu continu hyperélastique, isotrope et compressible aux données de la dynamique moléculaire pour un élément représentatif du polyéthylène. Dans la dernière section, le problème de la simulation numérique de la grande déformation viscoplastique d'un agrégat de polyéthylène est revisité.

2-Caractérisation élastique de l'interphase dans le polyéthylène semi cristallin

Les positions relatives des cristallites, des interphases et de la phase amorphe centrale sont illustrées sur le schéma d'un élément représentatif du polyéthylène dans la Figure 1. L'objectif de cette section est la caractérisation élastique de l'interphase. L'analyse de l'identification élastique de l'interphase dans le polyéthylène est fondée, d'une part, sur les données de la simulation Monte Carlo moléculaire qui était destinée au calcul du tenseur de rigidité de la zone interlamellaire il In 't Veld et al [15]. Les diagrammes de la variation des composants de il C avec la température sont donnés sur la ( ) C montre que les résultats des deux approches s'accordent parfaitement. Par ailleurs, la symétrie matérielle de l'interphase ressemble à celle de la zone interlamellaire. Puisque la symétrie matérielle de la phase amorphe centrale est isotrope et celle de la phase cristalline est orthotrope, on pourrait attendre que la symétrie de la phase cristalline contribue à la symétrie de l'interphase. Pour vérifier cette hypothèse, la forme générique du tenseur de rigidité de la phase cristalline, est transférée de son référentiel matériel au référentiel adopté pour l'analyse de la dissociation. Les deux référentiels en question sont schématiquement représentés sur la Figure 3. Il est constaté que la forme monoclinique obtenue après la transformation du tenseur de rigidité de la phase cristalline est identique avec celle de la zone interlamellaire et de l'interphase.

-     -     =   -     - -   -   ∼ C (FR.
1 1 1 1 1 1 1 1 ip ref ref ref il ref am am ip ip 1 η η η - - - - - - - - ∞ ∞ ∞ = - - - - - - +                            
Tableau 1 : Résultats des approches de dissociation dans la température type de 370 K.

A T=370 K: am 6.31MPa Ce résultat confirme la méthodologie proposée pour l'analyse de dissociation.

C C C C C C C C C C C C C C C C C C C C ′ ′ ′ ′     ′ ′ ′ ′   ′ ′ ′ ′     ′ ′     ′ ′ ′ ′   ′ ′    

3-Caractérisation hyperélastique du domaine interlamellaire et pour l'interphase dans le polyéthylène

L'objectif de cette section consiste en la caractérisation hyperélastique du domaine interlamellaire et de ses composants. L'approche adoptée consiste à appliquer les équations d'un milieu hyperélastique, isotrope et compressible aux données de la simulation de dynamique moléculaire disponible pour un élément unitaire de polyéthylène. Les détails de cette simulation atomistique, faite pour la température de 350 K, sont disponibles dans le travail de Lee et Rutledge [21]. Le schéma de ce dernier ainsi que les conditions aux limites sont illustrés sur la Figure 4. La distribution de la densité le long de l'élément avant d'imposer le chargement montre qu'il y a une transition entre la phase amorphe centrale et les cristallites. Cette transition est encore une preuve de l'existence de l'interphase dans la zone non cristalline. Pour développer les équations caractéristiques des zones hyperélastiques, données cidessous, leur tenseur de déformation respectif en liaison avec le modèle néo-Hookéen est remplacé dans la loi de comportement des milieux hyperélastiques. L'approche d'identification des paramètres hyperélastiques des domaines concernés est fondée sur la définition d'un ensemble de fonctions coûts non négatives qui viennent de la physique du problème. Pour la mise en oeuvre de cette idée, les frontières des deux interphases sont introduites comme les variables auxiliaires. Finalement, en employant un outil d'optimisation globale, le minimum global des fonctions coûts est trouvé et permet de déterminer les meilleurs paramètres hyperélastiques. µ = [22][23][24]. Bien que le module de compréssibilié ne soit pas aussi important que le module de cisaillement, les modules de compressibilié identifiés sont en accord avec les valeurs moyennes de 2 GPa et 3 GPa reportés ailleurs [22,[24][25][26].

Pour évaluer le succès de l'algorithme proposé, les diagrammes des frontières de deux interphases ainsi que leur allongement sont représentés sur la figure 5. Par ailleurs, les diagrammes des contraintes normales et latérales de la zone interlamellaire et celles de la zone amorphe centrale sont tracés sur la figure 6. La similarité entre la courbe de la moyenne arithmétique des deux contraintes normales latérales et la courbe analytique, plaide pour la validité de la méthodologie présentée. Finalement, une analyse de l'incertitude, destinée à déterminer le niveau d'incertitude des paramètres identifiés, est effectuée. Cette analyse a mis en lumière que l'incertitude des modules de cisaillement identifiés est moins importante que celle des modules de compressibilité.

4-Grande déformation viscoplastique du polyéthylène

Dans cette section, la grande déformation viscoplastique du polyéthylène est revisitée. Pour réaliser ce réexamen, l'élément unitaire de l'inclusion composite à deux couches, illustré sur la Figure 7, est utilisé. L'une des couches représente la cristallite et l'autre correspond à la zone interlamellaire, aussi appelé non cristalline. La continuité du champ de déplacement et l'équilibre de force à l'interface des deux couches sont imposés. L'imposition de ceux dernières conditions en liaison avec le théorème de la moyenne amènent aux formulations suivantes : La loi de comportement adoptée pour chaque phase est identique avec celle choisie dans les publications similaires [22,27]. L'élasticité des phases et leur sensibilité à la pression sont négligées. D'après des articles similaires [28], la loi de comportement des cristallites est supposée dépendante de la vitesse de déformation, α γɺ , comme suit : Dans la relation ci-dessus il H et il S sont respectivement le tenseur de contrainte de rappel et le tenseur de la contrainte déviatorique de la zone interlamellaire. Les trois modèles qui s'utilisent pour la contrainte de rappel sont les modèles 3-chaînes, 8-chaînes, et de Gent [29]. Pour la loi d'interaction, la loi modifiée de Taylor est exploitée. L'hypothèse sous-jacente de cette loi consiste en l'uniformité complète de la déformation du matériau hétérogène. L'état inextensible des chaînes polymériques dans les cristallites exige la modification suivante La simplification de l'inextensibilité susmentionnée peut s'exprimer comme un tenseur de projection d'ordre quatre, i I P , qui est multiplié à la vitesse de déformation de l'agrégat. Ce tenseur-là qui est plus complet que son prédécesseur [28], peut s'exprimer comme suit après des manipulations mathématiques.

3ch 1 1 1 1 3 1 2 1 2 3 1 3 3 i i r i C N N N N N λ λ λ λ σ λ λ λ λ - - - -             = - + +                         L L L L (FR.
i I i i i i i i 0 0 90 3 3 exp exp 2 2 m m ψ ψ ω ω ψ ψ         -     = - - ⊗ - - ⊗                

Figure 1 .

 1 Figure 1. Two schematics of relative positioning of polyethylene molecular chains in crystalline lamellae. Carbon and hydrogen atoms are represented by large dark and small bright bullets [14]. .......................................................................................................................................................................................... 7

8 Figure 4 .

 84 Figure 3. A) Simulated fully grown spherulite [16], B) Spherulites, as they start to appear and grow in polyethylene [17]. ............................................................................................................................................. 8 Figure 4. Schematic illustration of the configuration of chains within the noncrystalline domain. Segments corresponding to loops (both ends of segment terminate at the same crystal lamella; short dashed lines), bridges (the ends of the segment terminate at different lamellae; dotted lines) and tails (one end terminates at a crystal lamella and the other terminates in the noncrystalline domain; long dashed lines). ............................................................................................................ 9

12 Figure 6 . 11 C ≡ , il 22 C 12 C ≡ , il 13 C ≡▲ and il 23 C 15 Figure 7 . 16 Figure 8 . 20 Figure 9 . 24 Figure 10 . Diagrams of ip 44 C vs. il 44 C 44 C 26 Figure 11 . Diagrams of ip 46 C vs. il 44 C 46 C 26 Figure 12 . Diagrams of ip 66 C vs. il 44 C 66 C 27 Figure 13 . 66 C vs. il 44 CFigure 14 . 33 Figure 15 . 33 Figure 16 . 35 Figure 17 . 12 C , ip 13 C and ip 23 C 35 Figure 18 . 15 C , ip 25 C , ip 35 C and ip 46 C 36 Figure 19 . 41 Figure 20 . 45 Figure 21 . 48 Figure 22 . 54 Figure 23 . 54 Figure 24 .. 55 Figure 25 .

 126112212132315716820924104444442611464446261266446627136644143315331635171213233518152535463619412045214822542354245525 Figure5. A) Unit cell model introduced in[START_REF] Tomita | Computational characterization and evaluation of deformation behavior of spherulite of high density polyethylene in mesoscale domain[END_REF] for multiscale modeling of polyethylene; B) Schematic of multiscale modeling of PE introduced in[START_REF] Van Dommelen | Micromechanical modeling of intraspherulitic deformation of semicrystalline polymers[END_REF]; C) The building block of an aggregate of polyethylene along with the representation of local crystalline and convected axes[START_REF] Van Dommelen | Micromechanical modeling of the elasto-viscoplastic behavior of semi-crystalline polymers[END_REF]; D) The relative positioning of the major constituents of polyethylene microstructure on which the authors in[START_REF] Hütter | Monte Carlo Simulations of Semicrystalline Polyethylene: Interlamellar Domain and Crystal-Melt Interface[END_REF] based their atomistic model; E) schematic of a spherulite + the composite inclusion model often adopted in the literature[START_REF] Lee | Simulation of large strain plastic deformation and texture evolution in high density polyethylene[END_REF]. ........................................................................................................................................ 12

zzσ 56 Figure 26 .+ 56 viiFigure 27 . 58 Figure 28 . 60 Figure 29 . 60 Figure 30 . 61 Figure 31 . 61 Figure 32 . 69 Figure 33 . 70 Figure 34 . 78 Figure 35 . 79 Figure 36 . 79 Figure 37 . 80 viiiFigure A 1 .Figure C 1 . 99 Figure D 1 .Figure I 1 . 1 -

 562656275828602960306131613269337034783579367937801199111 Figure 26. Evolution of il xx σ , il yy σ and ( ) il il 0.5 xx

bCErrError of the least squares estimate i 1 eFGF 1 JlN"qQRSSλΛσττAngle between n and c xv 0 ψ

 110 ...................................................................................... xi Second lattice vector in the orthorhombic crystallites of polyethylene b= b Second lattice parameter of the polyethylene crystallite B Left Cauchy-Green deformation tensor c Shortest lattice vector in the orthorhombic crystallites in polyethylene c = c Lowest lattice parameter of the polyethylene crystallite r Rubbery modulus in 3-chain and 8-chain models λ C Fourth-order stiffness tensor of material or phase λ eff C Effective/homogenized stiffness tensor of some heterogeneous medium ref CThe stiffness tensor of the reference medium in double-inclusion method cr Subscript for "critical" Projection of D for inclusion "i" in the modified Taylor approach E λ Young's modulus of phase λ First unit vector of the orthonormal basis attached to the interface plane of inclusion "i" λ Deformation gradient tensor of phase λ g α Shear resistance of slip system α αβ,γ g Mutual effect coefficient G λ Elastic shear modulus of phase λ xii 0 N Plateau shear modulus of a rubbery polymer h αβ Components of the tensor of hardening moduli H Second-order back stress tensor of the noncrystalline domain The invariants of the left/right Cauchy-Green deformation tensor I Fourth-order identity tensor il / il Subscript / superscript for interlamellar (noncrystalline) domain ip / ip Subscript / superscript for interphase layer Jacobian or determinant of the deformation gradient tensor The main variable of Gent model defined as 2 The average length of each rigid link in the molecular chain network model of amorphous polymers Average number of rigid links in each chain in the molecular chain network model of amorphous polymers Tensorial representation of the lamella normal direction for inclusion "I" OF Nonnegative objective function (also known as cost function) Order parameter" that quantifies the orientational ordering in nematic systems I P Fourth-order projection tensor for inclusion "I" used in modified Taylor method Volume chain density i Strain concentration tensor of phase "i" in the composite inclusion method r Current length of a polymer chain in 3-chain and 8-chain models i Stress concentration tensor of phase "i" in the composite inclusion method R Ideal gas constant i Deviatoric stress tensor of phase "i" i S Fourth-order compliance tensor of phase "i" Eshelby tensor of a penny-shaped (or disk-like) inclusion (i.e. aspect ratio ∞ ) Root mean square of the principal stretches in a molecular chain network model for amorphous polymers i Matrix of coefficients of phase "i" in least squares method µ Shearing modulus of the hyperelastic constitutive law υ Mean lateral stress of the amorphous/interlamellar/interphase region averaged over the respective lateral dimension α The resolved shear stress on slip system α 0 Shear resistance of the easiest slip system ψ Threshold value for ψ ω Coefficient of exponential argument of I P

Figure 1 .

 1 Figure 1. Two schematics of relative positioning of polyethylene molecular chains in crystalline lamellae. Carbon and hydrogen atoms are represented by large dark and small bright bullets [14].

Figure 2 .Figure 3 .

 23 Figure 2. A) An electron micrograph of a chlorosulphonated section of a medium-density polyethylene with an insert cut through a computer-built spherulite [15], B) Schematic model of secondary crystalline branches (twigs) grown on a mother crystalline branch [16].

Figure 4 .

 4 Figure 4. Schematic illustration of the configuration of chains within the noncrystalline domain. Segments corresponding to loops (both ends of segment terminate at the same crystal lamella; short dashed lines), bridges (the ends of the segment terminate at different lamellae; dotted lines) and tails (one end terminates at a crystal lamella and the other terminates in the noncrystalline domain; long dashed lines).

Figure 5 .

 5 Figure 5. A) Unit cell model introduced in [46] for multiscale modeling of polyethylene; B) Schematic of multiscale modeling of PE introduced in [47]; C) The building block of an aggregate of polyethylene along with the representation of local crystalline and convected axes [48]; D) The relative positioning of the major constituents of polyethylene microstructure on which the authors in [53] based their atomistic model; E) schematic of a spherulite + the composite inclusion model often adopted in the literature [50].

Figure 6 . 11 C ≡ , il 22 C ≡▲ and il 33 C ≡ B) il 12 C ≡ , il 13 C ≡▲ and il 23 C

 6112233121323 Figure 6. Variation of A) il 11 C ≡ , il 22 C ≡▲ and il 33 C ≡

Figure 7 .

 7 Figure 7. Variation of the thickness of A) interlamellar domain B) interphase layer in semicrystalline polyethylene with temperature [42].

eM

  the molecular mass between entanglements. Bédoui et al.[START_REF] Bédoui | Micromechanical modeling of isotropic elastic behavior of semicrystalline polymers[END_REF] and Sedighiamiri et al.[START_REF] Sedighiamiri | Micromechanical modeling of the elastic properties of semicrystalline polymers: A threephase approach[END_REF] used this relationship and estimated am 1.5MPa G = for room temperature. They both assumed a mean value of am 3000MPa

Figure 8 .

 8 Figure 8. Schematic illustration of the relative positioning of a polyethylene stack along with the known and unknown parameters of the presented dissociation analysis.

Figure 9 .

 9 Figure 9. Visual demonstration of the representation of the crystalline stiffness in its lattice abcframe and when rotated to the xyz-frame.

66 CThe shearing stiffness ip 55 C as a function of il 55 C has a behavior similar to ip 44 C as a function of il 44 C 55 C shows a weak dependence on the uncertainty of il 55 C 44 C and ip 55 C

 665555444455554455 Figure 12. Diagrams of ip 66 C vs. il 44 C . Within the approximate interval ( ) 15.5,16 delineated by

44 C and ip 55 C 66 C 66 C 44 C or ip 55

 5566664455 are negative, they are very close to zero in magnitude, compared to the other stiffnesses. The closeness to zero is such that they can be assumed independent of both temperature or the corresponding component in il C . But the situation for ip is totally different (see Figure12); since within the interval of takes incomparably negative values that are at least three orders of magnitude larger than ip

44 C 66 C 44 C and ip 55 C , negative values are not allowed for ip 66 C 44 C and ip 55 C 66 CThe sensitivity analysis of ip 66 C 66 C 66 C 66 C , the admissible interval of il 44 C will be further reduced to il 44 . 44 C 44 C

 44664455664455666666666644444444 approaches the right extreme of its allowable interval, the plateau-like behavior of ip is observed, supporting the speculation that, contrary to ip ; this comparison suggests that the allowable interval for il 44 C should shrink to ( ] 73,100 . It might also be argued that since ip are shearing resistances in the planes normal to the interface but ip is the shearing resistance in the plane parallel to the interface, the dissimilarity between their behaviors could be expected. has revealed that this shearing stiffness is robust to the uncertainties of all components of il C except for the uncertainties of il 44 C and il 66C . Accordingly, the diagrams of Figure13are plotted as follows: at different temperatures and across the interval of is varied within its uncertainty interval to find the minimum values of ip at each temperature. In light of the reasoning made in the preceding paragraph vis-à-vis the inadmissibility of negative values for ip This is the reason why in Table1the value of 90 MPa is assigned to il in lieu of its mean value calculated from MC molecular simulations. Consequently, for the rest of the calculations, the mean value and the uncertainty interval of il are, respectively, taken to be 90 MPa and (82,100) MPa. It is worth noting that the dissociation analysis has the unintended but useful by-product of confining the most uncertain component of il C .

Figure 13 . 66 C vs. il 44 C

 136644 Figure 13. Diagrams of the minimum values of ip 66 C vs. il 44 C at different temperatures. The arrow

  quickly but smoothly such that

(

Figure 14 .

 14 Figure 14. Average Young's modulus of the interphase layer calculated using Hill's estimate, ip-Hill E

Figure 15 .

 15 Figure 15. Average Young's modulus of the interlamellar domain estimated using Hill's estimate, il-Hill E , and the bounds of its uncertainty interval vs. temperature.

  uncertainties available in il C and am C , the following Monte Carlo sensitivity analysis has been carried out. First, the following uncertainty intervals were considered for the amorphous Poisson's ratio, am υ , and the amorphous Young's modulus, am E , for the temperature range 350K-400K.

Furthermore, the uncertainty intervals of the ip C components originated from the uncertainties of the il C, ip 33 C is the less sensitive component and ip 11 C

 il3311 components are calculated via the same Monte Carlo procedure: at a given temperature, the uncertainty space of il C is sampled 10 9 times, ip C s associated with each temperature are calculated, and finally the bounds of each component of ip C at each temperature are obtained (see Figure16-Figure17). Of the tensile components of ip C is the most sensitive component. Additionally, the non-orthotropic elements of ip C take small values close zero and exhibit a weak dependence on temperature.

Figure 16 .Figure 17 . Uncertainty intervals of ip 12 C , ip 13 C and ip 23 CFigure 18 .shearing stiffnesses ip 44 C and ip 55 CC have no impact on il 44 C , a weak impact on il 55 C and a strong impact on il 66 C , which may be viewed as more evidence for the dissimilarity between ip 66 C 44 C and ip 55 C 44 C and il 55 C

 16171213231844554455666644554455 Figure 16. Uncertainty intervals of normal components of ip C vs. temperature.

  A) Schematic illustration of a PE stack and the relative positioning of the constituents. The arrows indicate the direction of the applied displacement constraints in the MD simulations. B) Edited schematic of the simulated PE stack in its unloaded state together with the average density distribution along the stack [89]. The initial, unloaded length of the simulation PE stack (i.e. crystalline+noncrystalline phases) is 21.63 nm and it is loaded by imposing the longitudinal tensile stretch along the z-axis, and simultaneous compressive lateral stretches x λ and y λ . The xyz- frame is defined as described in Figure 19. The deformation of the stack is recorded in 800 increments up to the longitudinal strain 0.4917 z ε = during which the lateral and longitudinal stresses are recorded. This information, i.e. longitudinal stretch and stress of the entire stack, z λ and zz σ , respectively, together with the lateral stress distributions, xx σ and yy σ , constitutes the database for the identification procedure

  am z λ -and ip z λ -are treated as unknown functions of z λ (or equally z ε ) and are identified along with the hyperelastic coefficients.

  am/il xx σ and am/il yy σ are analytically identical, although the MD simulation results do not show such an identity. The lateral stresses xx σ and yy σ along the simulation stack have been recorded

Figure 20 .

 20 Figure 20. Distributions of the lateral stresses am xx σ and am yy σ along with the dimensionless order parameter P2 of the entire stack at twelve sample increments from MD simulations.

Figure 21 .

 21 Figure 21. Evolution of zz σ with z ε recorded during the MD simulations [89]. The slope of the curve changes sign at the extrema points

  ) imply z-ip-l z-ip-r λ λ = , where in conjunction with the first relationships results in ip-l ip-r lat lat σ σ = . Therefore, the third objective function is defined as follows

Figure 22 .

 22 Figure 22. Evolution of the amorphous and interlamellar boundaries with longitudinal strain. The boundaries are identified by minimizing

Figure 23 .

 23 Figure 23. The evolution of the right and left interphase stretches as defined in Eq.( 32) . The analytical curve of ip λ calculated by substituting the identified ip µ and ip κ into the second relationship of Eq.( 30) agrees well with the two others.

Figure 24 ..

 24 Figure 24. Diagrams of the average lateral stresses of the left and right interphase layers, ( ) ip-r/l ip-r/l 0.5 xx

Figure 25 .

 25 Figure 25. Evolution of am xx σ , am yy σ and ( ) am am 0.5 xx

Figure 26 .

 26 Figure 26. Evolution of il xx σ , il yy σ and ( ) il il 0.5 xx

.

  point the chains are fully aligned with the deformation direction.The fully aligned state continues up to 0.3 z ε ∼ , at which point the chain segments revert to their initial tilted state, albeit imperfectly. The stretch diagrams of the three domains along with the affine stretch line which are superposed in Figure27clearly show this phenomenon. Keeping in mind that in the absence of any phase transformation, the following relationship holds il zvolume fraction of the associated phase, one expects that the stretch of the least rigid phase will be higher than z λ . In comparison with the crystalline phase, the amorphous/interlamellar domain is more compliant, but the diagrams of Figure27suggests that z-cr λ , z-am λ and z-ip λ have almost the same contribution up to ~0.21 z ε Without resorting to the details of the MD observations andgiven the higher rigidity of the crystallites, the likeliest mechanism responsible for such a level of contribution from the crystalline phase is the rotation of the crystalline chains and their alignment in the direction of the applied stress. For 0still lower. The associated part of the stretch diagrams agrees well with the MD observation where melting takes place at the crystallinenoncrystalline interface within this straining interval. In other words, with the rotation mechanism being exhausted at ~0.21 z ε , the likeliest mechanism that contributes substantially to the overall stretch of the stack and lowers the contribution of il/am/ip z λ - is the melting of the crystalline phase. One may reason that looser packing and consequently lower density of the newly added, noncrystalline portion relative to the previously crystalline phase makes the melting phenomenon equivalent to a local dilatation which is responsible for lowering the contribution of il/am/ip z indicative of a rise in the contribution of the stretch from the corresponding phases. This part of the diagrams is consistent with the reversion of the chain stems from the aligned state to the slant state, due possibly to the interlamellar shear that develops at this stage[START_REF] Ward | An Introduction to the Mechanical Properties of Solid Polymers[END_REF].

Figure 27 .

 27 Figure 27. Evolution of am λ , il λ and ip λ with z ε . The three stretches evolve very closely. z λ of the

  section, first the authenticity of the arithmetic mean of the two lateral stresses within the amorphous interval, i.e. characterization is investigated. Then an analysis is conducted to examine the level of certainty of the identified hyperelastic coefficients.

  , the left and right ends of the interlamellar and amorphous domains are bounded to [-55..-50] Å, [48..53] Å, [-40..-35] Å and [33..38] Å, respectively. Therefore, as illustrated in Figure 28B, the shortest and longest amorphous intervals are confined to [-35..33] Å and [-40..38] Å, respectively. Similarly for the interlamellar domain, the shortest and longest intervals are confined to [-50..48] Å and [-55..53] Å, respectively. Assuming affine stretch evolutions for both amorphous and interlamellar regions, i.e. the shortest/longest amorphous/interlamellar interval is obtainable.The diagrams of am/il xx σ , shortest and longest intervals are given in Figure29and Figure30, respectively. In conjunction with the diagrams of Figure25and Figure26, it is concluded that the diagrams of () the typical analytical diagrams of the lateral stress

  stress distribution along the PE stack at the first increment of loading. The picks of the distribution denote the centers of the interphase layers. B) The shortest and longest amorphous and interlamellar intervals at the first loading increment which is assumed to be valid for the no-loading state.

Figure 29 .Figure 30 .Figure 31 .

 293031 Figure 29. Evolution of the average lateral stress distributions am xx σ , am yy σ and

  There are numerous experimental techniques available in the literature on the quantification of the crystallinity in polyethylene samples [104, 105]. Concerning the orientation of the chain axis of each inclusion, I c , and its lamella normal, I n , they can be experimentally traced by means of WAXS pole figures and 2D SAXS, respectively [106]. The molecular orientation of the noncrystalline domain, on the other hand, can be found by virtue of WAXS measurements in conjunction with pole figure technique and separation of the scattering produced by the crystalline and amorphous components [107].

Figure 32 .

 32 Figure 32. A cubic unit cell of three-chain model in the (a) undeformed, (b) uniaxial tension, and (c) biaxial tension configuration.

Figure 33 .

 33 Figure 33. A cubic unit cell of eight-chain model in the (a) undeformed, (b) uniaxial tension, and (c) biaxial tension configuration.

  -, at which the extension capacity of the rubber is exhausted. The associated Cauchy stress formulation is rendered into the following simpler form after substitution into relation[START_REF] Balijepalli | Molecular simulation of the intercrystalline phase of chain molecules[END_REF]

P

  denotes the reciprocal of the orientational averaging carried out over all I P . Using the modification introduced in ( 83), the global compatibility, i.e. I = D D ɶ

  . eq ε for various sets of parameters are plotted in the following figures. The material parameters associated with each set of diagrams are given below the figures. All variables and parameters of stress dimension are normalized to 0 τ . Similarly all variables and parameters of strain rate dimension are normalized to 0 γɺ to simplify numerical computations.

  is a more appropriate choice and helps to adjust the convexity of the curves. Moreover, with 200 N ∼ the locking strain is more effectively accommodated.

Figure 34 .

 34 Figure 34. The stress-strain curves plotted using different back stress models. These diagrams are stiffer than G'sell's experimental data with the adopted model parameters:

Figure 35 .Figure 36 .Figure 37 .

 353637 Figure 35. Stress-strain diagrams with Gent model at fixed 50 N = and eq 0 1 D γ = ɺ . Light green curve:

il a , N , r C and 0 ψ

 0 for different rates of straining and under different modes of deformation. This job is the subject of a future study. After climbing a great hill, one only finds that there are many more hills to climb. (Nelson Mandela 18 July 1918 -5 December 2013)

ip 66 C

 66 shows a different behavior, due possibly to its resistance in the plane parallel to the interface. As a valuable by-result, the dissociation analysis has had the favorable advantage of constraining the most uncertain component of the initial interlamellar stiffness, il 44 C . This particular component was the most uncertain one with an uncertainty interval 100 MPa ±

  averaging of an arbitrary fourth-order tensor, C , is calculated as follows relation, the overscore and the angle brackets denote the average of the quantity in question. The above integration can be expanded as follows:Because a fourth-order tensor in the general case has 81 nonzero components, each component in the transformed tensor, ijkl ′ C , consists of the sum of 81 terms which is the expanded form of im jn kp lq mnpq a a a a C . Therefore, if the transformed components are all collected in a column vector, column vector containing the elements of the initial tensor, a simpler representation of relationship (A6).

C

  acts as a body force because the equilibrium equation is rendered to ijkl k,lj i -= C u χ in the presence of the body force χ and zero eigen strain. Similarly, * ijkl kl j ε C n acts as a surface traction on the outer boundary. Therefore, one can say that the elastic displacement filed created by * ε in a free body is identical with the one created by the body force

′

  x -x is the Dirac delta function and mk δ is the Kronecker delta. Solution of the above differential equation returns G whose substitution into the following integral yields u .

4 :

 4 If the reference medium is isotropic and 0 τ → (thin, disk-like inclusions):

Fig C 1 .

 1 Fig C 1. Schematic illustration of the double-inclusion problem. There are a very inner ellipsoidal inclusion, indicated by "Ω", and a covering ellipsoidal inclusion, indicated by "Γ". The entire double-inclusion, indicated by "DI", is inserted in a reference medium, indicated by "ref".

  the amorphous phase by the inner encapsulated inclusion, i.e. to the following relationship

  il C there is no trace of ref E and ref υ which shows the independence of il C from ref C , when am C , ip C and am η are known. Similarly, it can be concluded that when the inclusions of the double-inclusion problem are thin and disk-like, ip C is independent of ref C if il C , am C and am η are known.

Fig D 1

 1 Fig D 1 Schematic of a three-layer composite inclusion along with the relative orientation of the selected reference frame.

3 α′σ

 3 from Eq.(D14), Eqs.(D2) and (D3) are exploited to obtain ( ) ( )

  as well as the associated volume fractions.

Figure I 1 . Comparison of different estimates proposed for 1 -

 11 Figure I 1. Comparison of different estimates proposed for 1 -L with the exact numerical solution

  effet de l'interphase.

  C . Ce dernier ainsi que l'épaisseur de la zone interlamellaire, notée il t , et celle de l'interphase, notée ip t , ont été calculés pour l'intervalle de température de 350 K à 450 K. Les détails de cette étude atomistique ainsi que ceux des calculs des il C , il t et ip t sont fournis dans l'article publié en 2006 par

Figure 1 :

 1 Figure 1 : Illustration schématique du positionnement relatif d'un élément de polyéthylène avec les paramètres connus et inconnus de l'analyse de dissociation.

Figure 2 .Figure 2 : 11 C ≡ , il 22 C ≡▲ and il 33 C ≡ B) il 12 C ≡ , il 13 C ≡▲ and il 23 C

 22112233121323 Figure 2. Pour les éléments restants, les valeurs calculées à la température de 435 K sont adoptées, faute d'avoir l'accès à de meilleures approximations. Le tenseur de rigidité interlamellaire évalué à la température de 435 K est donné ci-dessous. La forme de ce tenseur suggère l'existence d'une symétrie monoclinique, ou peut-être d'une symétrie plus haute. Les causes éventuelles de cette observation sont expliquées dans la suite.

6 )

 6 Dans la première relation, i.e. celle de ECIM , ip R et am R s'appellent les coefficients de concentration de contrainte. Dans la deuxième relation, i.e. celle de DIM, I représente le tenseur d'identité d'ordre quatre, ∞ S représente le tenseur d'Eshelby d'une inclusion discoïde et fine, et ref C indique le tenseur de rigidité du milieu de référence contenant la double inclusion. Les hypothèses détaillées de ces deux approches sont fournies dans les annexes C et D du mémoire. Le tenseur de l'interphase, noté ip C , produit par les approches DIM et ECIM pour la température type de 370 K est présenté dans le Tableau 1. Les valeurs négatives diagonales de ip C indiquent son état non défini positif. Il a été mis en évidence que ces deux valeurs négatives sont robustes aux incertitudes de il C et à celles de am C et donc l'interphase est mécaniquement instable. En outre, la comparaison rapide entre les deux ip

5 *

 5 Cette rigidité de cisaillement a été réglée à 90 MPa, au lieu de sa valeur moyenne rapportée dans l'étude de la simulation moléculaire, pour des raisons expliquées ultérieurement.

Figure 3 :

 3 Figure 3 : La forme du tenseur de rigidité de la phase cristalline dans son référentiel matériel, abc, ainsi que dans le référentiel de l'analyse, xyz.

7 ) 9 )=.

 79 Dans le polyéthylène à haute cristallinité, le domaine interlamellaire est dominé par l'interphase et donc ip il E E ≃ . Le module élastique de l'interphase peut être estimé en utilisant l'approximation de Hill comme suit.En calculant la moyenne volumique de ip C à 350 K, qui est la température la plus proche de la température ambiante pour ensuite la remplacer dans les relations ci-dessus, on obtient ip-Hill 347MPa E Cette estimation est en bon accord avec la valeur asymptotique de 300 MPa proposée dans la référence[18]. De plus, si les incertitudes de il C sont prises en compte, l'intervalle de ip-Hill E est déterminé comme suit:

Figure 4 A

 4 Figure 4 A) : Illustration schématique de l'élément unitaire de polyéthylène et du positionnement relatif des constituants. Les flèches indiquent la direction des déplacements appliqués dans la simulation atomistique. B) Schéma de l'élément unitaire de polyéthylène simulé dans son état non chargé, avec la distribution de la densité moyenne le long de l'élément.

  contrainte normale, moyenne et latérale, et la contrainte normale et longitudinale de la phase amorphe / du domaine interlamellaire. am/il µ et am/il κ sont respectivement le module de cisaillement et le module de compressibilité.

  longitudinal de l'ensemble de l'élément et -am/il z λ signifie l'allongement latéral de la phase amorphe / du domaine interlamellaire.

Figure 5 :Figure 6 :,

 56 Figure 5 : Evolution de la frontière droite et gauche de la zone amorphe et interlamellaire avec la déformation longitudinale.

Figure 7 : 13 ).

 713 Figure 7 : Élément unitaire de l'inclusion composite à deux couches employé pour étudier la déformation viscoplastique du polyéthylène.

  la vitesse de déformation de référence, α τ la contrainte de cisaillement résolue du système de glissement α , n exposant de vitesse de déformation, et g α la résistance de cisaillement du système de glissement α . Une loi simple de puissance viscoplastique pour le tenseur de la vitesse de cisaillement plastique, il D , s'écrit comme suit :

ψ

  représente l'angle entre la direction des chaînes et la direction perpendiculaire à l'interface cristalline-non cristalline dans l'inclusion ''i'', et 0 ψ est sa valeur de seuil. déviatorique du produit dyadique de la normale de l'inclusion ''i''déviatorique du produit dyadique du vecteur unitaire dans la direction ''1''.

1 ω

 1 ≫ et m font parties des paramètres du modèles. Le positionnement schématique des chaînes polymériques ainsi que les référentiels employés sont illustrés sur la Figure 8.

Figure 8 :

 8 Figure 8 : Schéma du positionnement des chaînes polymériques de la cristallite et référentiels employés dans l'analyse de grande déformation d'un agrégat de polyéthylène.

Figure 9 :Figure 10 :Figure 11 :

 91011 Figure 9 : Les courbes de contrainte-déformation tracées pour divers modèles de la contrainte de rappel. Ces courbes avec les paramètres

Table 1 .

 1 Output of the dissociation approaches at the typical temperature 370 K. ............................

Table 2 .

 2 The effect of negligible terms of ip C on the homogenized il C .................................................

Table 3 .

 3 Impact of nonorthotropic components of il C on the dissociated ip C .................................

	Table

Table 1 . Output of the dissociation approaches at the typical temperature 370 K.

 1

					at T=370 K: am 6.31MPa E =	, am 2.10 MPa G =	, am 0.66 η =	, ip 0.34 η =
	C	( am MPa	)						C	( il MPa	)	(taken from [42])
	3097.9 3093.7 3093.7 0	0	0	1749.9 1613.6 1092.9	0	-180	0
	3093.7 3097.9 3093.7 0	0	0	1613.6 2569.3	1150	0	-240	0
	3093.7 3093.7 3097.9 0	0	0	1092.9 1150	1249.6	0	50	0
		0			0		0	2.1	0	0		0	0	0	90 *	0	-200
		0			0		0	0	2.1	0	-180	-240	50	0	220	0
		0			0		0	0	0	2.1		0	0	0	-200	0	570
	Then:
											1972.2	1307.4	427.7	0	3.05	0
											1307.4	3824.3	524.2	0	3.92	0
			C	( ip M Pa	)	output by the DIM		427.7 0	524.2 0	567.3 0	0 -1.11	-0.33 0	0 7.31
											3.05	3.92	-0.33	0	-1.08	0
											0		0	0	7.31	0	320.5
											1972.2	1307.4	427.7	0	3.05	0
		C	( ip MPa hybrid optimization algorithm) ) output by the ECIM (using the	1307.4 427.7 0 3.05	3824.3 524.2 0 3.92	524.2 567.3 0 -0.33	0 0 -1.11 0	3.92 -0.33 0 -1.08	0 0 7.31 0
											0		0	0	7.31	0	320.5

* This

particular shearing stiffness was taken to be 90 MPa, as a safe value instead of its mean value reported in the molecular simulation study, for the reasons explained in Subsection 2.4.2.

  The most critical elements of il C are

					il 44 C , with an
	uncertainty of 100 MPa ±	, and	il 55 C , with an uncertainty of 60 MPa ±	. According to the
	sensitivity analyses carried out, the following observations were made:
	• Variation of	il 44 C within its uncertainty interval brings about the variation of ip 44 C ,
	ip 46 C and ip 66 C , while the other components of ip C are robust to this variation.
	• All components of il C , except for	il 44 C , were allowed to vary within their interval
	of uncertainty and were observed to affect the value of	ip 44 C only in the tenth
	decimal place. Similarly, holding	il 55 C fixed and varying the other components of
	il C within their uncertainty intervals was observed to affect the value of	ip 55 C only
	to the sixth decimal place. It can therefore be concluded that	ip 44 C / ip 55 C is only
	affected by the uncertainty in the corresponding	il 44 C / il 55 C and is very robust to
	the uncertainty of the other components of il C .
	• Variation of	il 66 C only varies ip 66 C , with the other components of ip C remaining
	robust to the fluctuations in	il 66 C .
	Accordingly, the diagrams of	ip 44 C , ip 46 C and ip 66 C vs. the uncertainty interval of il 44

Figure 10. Diagrams of ip 44 C vs. il 44 C . Within the approximate interval ( )

											2.2, 4
									2.2, 4 , (	) 2.5, 4 and (	-	15.5,16	)	must be excluded from the initial
	interval of il 44 C . Therefore, the rough interval	-	15.5 M Pa	≤	C	il 44	≤	16 M Pa	, within which
	ip 66 C >	10000MPa	, is excluded from the initial uncertainty interval of il 44 C . On the other hand,
	within [	] 100, 15.5 --	, ip 66 C >	4000 MPa	, which might be considered incomparable with the
	shearing components of	il C . Additionally, it looks rather unusual to assume that the
	admissible interval of il 44 C consists of two separate intervals, i.e.	[ C ∈ -il 44 100, 15.5 ] [ 16,100 ] -∪	.
	Therefore, if one sets the criterion for the admissibility of il 44 C to ip 66 C <	4000 MPa	, then the
	allowable interval of	C	il 44	shrinks to	26.5 MPa	≤	il 44 C	≤	100 MPa	since for
	16 M Pa <	il 44 C <	26.5 MPa	, ip 66 C takes values less than -4000 MPa.

delineated by dashed lines, ip 44 C takes incomparably large value. Figure 11. Diagrams of ip 46 C vs. il 44 C . Within the approximate interval ( )

	2.5, 4

delineated by dashed lines, ip 46

Table 2 . The effect of negligible terms of

 2

			at T=370 K: am 6.31MPa E =	, am 2.1MPa G =	, am 0.66 η =	, ip 0.34 η =
	( ip MPa terms beyond orthotropic symmetry) ) C (Taken from Table 1 after eliminating nonzero					C	( il MPa	)
	1972.2 1307.4 1307.4 3824.3 427.7 524.2	427.7 524.2 567.3	0 0 0	0 0 0		0 0 0	⇒	1600.5 1414.3 1134.4 1414.3 2304.6 1205.3 1134.4 1205.3 1238.1	0 0 0	0 0 0	0 0 0
	0	0	0	-1.11	0		0			0	0	0	90	0	0
	0	0	0	0	-1.08		0			0	0	0	0	220.7	0
	0	0	0	0	0	320.5			0	0	0	0	0	109.3
	( ip MPa negligible elements) ) C (Taken from Table 1 after eliminating small						C	( il MPa	)
	1972.2 1307.4	427.7	0	0		0		1600.5 1414.3 1134.4	0	0	0
	1307.4 3824.3 427.7 524.2	524.2 567.3	0 0	0 0		0 0	⇒	1414.3 2304.6 1205.3 1134.4 1205.3 1238.1	0 0	0	0

ip

C on the homogenized il C

Table 3

 3

confirms that ip C s calculated using the two approaches match perfectly. The tensile elements of the newly calculated ip C are still close to their corresponding components of ip C , given in Table 1 and calculated using a il C of monoclinic symmetry. Contrary to the situation in Table 1, here the ECIM converges to the same solution output by the DIM using the numerical recursive method. Indeed, a combination of dual formulae ( 12) is used to achieve the convergence. More interestingly, it is observed that the controversial shearing stiffnesses ip 44 C and ip 55 C appearing in Table 3 are equal to those appearing in Table 1, suggesting that these two shearing stiffnesses are almost independent of the non-orthotropic components of il C . The last shearing stiffness, ip 66 C , however, exhibits a strong dependence on the non-orthotropic components; another dissimilarity which might have been expected in advance.

Table 3 . Impact of nonorthotropic components of il C on the dissociated ip C

 3

					at T=370 K: am 6.31MPa E =	, am 2.1MPa G =	, am 0.66 η =	, ip 0.34 η =
	C	( am MPa	)						il C (taken from [42] after eliminating nonorthotropic elements)
	3097.9 3093.7 3093.7	0	0	0		1749.9	1613.6 1092.9	0	0
	3093.7 3097.9 3093.7	0	0	0		1613.6	2569.3	1150	0	0
	3093.7 3093.7 3097.9	0	0	0		1092.9	1150	1249.6	0	0
		0	0	0	2.1	0	0			0	0	0	90	0
		0	0	0	0	2.1	0			0	0	0	0	220
		0	0	0	0	0	2.1			0	0	0	0	0	570
	Then:
							2569.8	2074.1	362.03	0	0	0
							2074.1	4807.8	439.96	0	0	0
	C	( ip MPa	)	output by the DIM	362.03 0	439.96 0	574.47 0	0 -1.11	0 0	0 0
							0		0	0		0	-1.08	0
							0		0	0		0	0	1688.6
							2569.8	2074.1	362.03	0	0	0
	C		( ip MPa using numerical recursive ) output by the ECIM method	2074.1 362.03 0 0	4807.8 439.96 0 0	439.96 574.47 0 0	0 0 -1.11 0	0 0 0 -1.08	0 0 0 0
							0		0	0		0	0	1688.6

)

  In summary, am/il

		am/il lat σ is the
	average, lateral, normal stress of the amorphous/interlamellar phase. -am/il z λ , however, is
	a dependent variable that is an unknown function of z λ (or equally z ε ) and is required
	for evaluating	am/il lat σ .

µ and am/il κ are hyperelastic (neo-Hookean) parameters that are to be identified, am/il zz σ and z λ are the normal stress and stretch of the stack in the z- direction, respectively, which are outputs of the MD simulations, and

  . According to this relation, am/il/ip

			λ		is a direct function of
	zz σ which shows two sharp, angle-like corners at	z ε =	0.2165	and	z ε =	0.3046	, indicating
	that it can be split into three separate parts, each describable with a simple estimate
	such as polynomial. For a better comprehension of this selection, the evolution of zz σ
	obtained from the molecular simulation is given in Figure 21. Since am/il z λ -is a direct
	function of zz

  37, 38,..., 48

	i =	for the left interlamellar boundary:
	37	54.9 Å,	38	59.3Å,	39	44.5 Å,	40		97.9 Å,
	41	56.9 Å,	42	51.8 Å,	43	31.5 Å,	44		5.3 Å,
	45	23.9 Å,	46	277.6 Å,	47	556.7 Å,	48	463.4 Å

  and

						am/il λ -is calculated z
	numerically. Then	am/il λ -, am/il z µ and am/il κ are fed into the first relationship of the system
	of equations ( 30) to calculate the associated lateral stress distribution. It can be seen
	that when	κ >	0.3 GPa	or	5 MPa

Table 4 . Slip systems of the crystallite in high density polyethylene along with their respective slip resistances adopted from [48, 50, 55].

 4

	Type of slip	Slip system	Normalized shear
	Chain slip Transverse slip	( ) 100 001     ( ) 010 001     { } 110 001     100 010 ( ) 010 100 ( ) { } 110 1 10	strength 1.0 2.5 2.5

  

								1111
								1122
								1133
								1212
								1313
								2211
					2 6 2 8	2 2 6	2222
	2 4 4	4	4 4 2 4	4 4 4 2	2233
	2	1 1 6	6	1 2	1 6	1 1 2	2323
	2 4 4	4	4 4 2 4	4 4 4 2	3311
	2 4 4	4	4 4 2 4	4 4 4 2	3322
	6 2 2 8	8 2 6 2 8	2 2 6	3333

  For the Eshelby tensor of a disk-like inclusion, the relationships of Eq.(B21) have been used and for the 6×6 matrix representation, either of the matrices of the left hand side of Eq.(E10) can be employed. Interphase layer and amorphous phase are assumed to have full anisotropic symmetry, that is

		am	am	am	am	am	am
		11	12	13	14	15	16
		am	am	am	am	am	am
		12	22	23	24	25	26
		am	am	am	am	am	am
	am	13	23	33	34	35	36
		am	am	am	am	am	am
		14	24	34	44	45	46
		am	am	am	am	am	am
		15	25	35	45	55	56
		am	am	am	am	am	am
		16	26	36	46	56	66

ref E and ref υ are, respectively, Young's modulus and Poisson's ratio of the reference medium.

  

							2131	2132	2133
	2211	2212	2213	2221	2222	2223	2231	2232	2233
	2311	2312	2313	2321	2322	2323	2331	2332	2333
	3111	3112	3113	3121	3122	3123	3131	3132	3133
	3211	3212	3213	3221	3222	3223	3231	3232	3233
	3311

  

		122	2123	2131	2133
	2211	2212	2223	2231	2233
	2311	2312	2323	2331	2333
	3111	3112	3123	3131	3133
	3211	3212	3223	3231	3233
	3311	3312	3323	3331	3333

Résumé

  Dans ce travail de thèse, la caractérisation élastique et hyperélastique de l'interphase dans le polyéthylène a été abordée, suivi par le réexamen de la simulation en grande déformation viscoplastique d'un agrégat de polyéthylène. La caractérisation élastique de l'interphase est mise en oeuvre en appliquant deux approches micromécaniques distinctes à partir des données disponibles de la simulation moléculaire de type Monte Carlo pour la zone interlamellaire. Ces approches micromécaniques sont d'une part le modèle étendu d'un composite avec inclusion, et d'autre part la méthode de double inclusion. Après avoir établi le tenseur de rigidité de la phase amorphe en fonction de la température, le tenseur de rigidité de l'interphase est dissocié de celui de la zone interlamellaire en utilisant les deux approches mentionnées ci-dessus. Les résultats de celles-ci s'accordent parfaitement. Il a été mis en évidence qu'à l'inverse des matériaux conventionnels, le tenseur de rigidité de l'interphase n'est pas défini positif. Cette caractéristique indique l'instabilité mécanique de l'interphase. Une analyse de sensibilité ad hoc est faite afin d'évaluer l'impact des incertitudes des éléments du tenseur de rigidité interlamellaire sur celui de l'interphase. Par ailleurs, le module d'élasticité effectif de l'interphase estimé à partir de son tenseur de rigidité est comparable à celui de la phase interlamellaire dans le polyéthylène à haute cristallinité, module évalué expérimentalement. Cette similarité, ainsi que les résultats identiques des deux approches micromécaniques, confirment les nouvelles informations obtenues sur les propriétés élastiques de l'interphase et incitent à appliquer la méthodologie proposée à des problèmes similaires.La caractérisation hyperélastique du domaine interlamellaire et de la couche d'interphase est ensuite élaborée. Quand le polyéthylène est soumis à de grandes déformations, l'interphase et la phase amorphe centrale se comportent de façon hyperélastique. L'algorithme hybride proposé consiste à appliquer la loi de comportement d'un milieu continu isotrope, compressible et hyperélastique aux résultats de la simulation de la dynamique moléculaire d'un élément unitaire de polyéthylène. En supposant un modèle néo-Hookéen, les équations caractéristiques sont obtenues et permettent l'identification des paramètres hyperélastiques pour la phase amorphe centrale, l'interphase et la zone interlamellaire, à l'aide de la notion d'optimisation d'un ensemble de fonctions coûts non négatives. Les paramètres hyperélastiques identifiés pour la zone interlamellaire sont en bon accord avec ceux qui ont été estimés expérimentalement. L'analyse d'incertitude a mis en évidence que l'incertitude des modules de cisaillement est moins importante que celle des modules de compressibilité. Il a été constaté que les phénomènes de fusion, recristallisation, et rotation des chaînes cristallines qui ont lieu durant la déformation peuvent être identifiées. L'évolution des frontières de l'interphase avec la déformation est le résultat secondaire de cette analyse d'identification.La fin du travail est dédiée à la simulation numérique de la grande déformation viscoplastique d'un agrégat de polyéthylène avec ces nouveaux apports. Le modèle adopté pour la contrainte de rappel, le tenseur de projection proposé pour l'approche modifiée de Taylor, et l'imposition des inégalités liées à la loi de comportement des lamelles cristallines qui auparavant avaient été laissées de côté, font parties des contributions apportées. Concrètement, dans les relations développées, le modèle de Gent est employé pour la contrainte de rappel du domaine non cristallin, à la place du modèle 3-chaînes ou 8-chaînes. En outre, un tenseur de projection pour l'approche modifié de Taylor est proposé, qui est plus complet que son prédécesseur. De plus, dans le code numérique relatif à cette méthodologie, la borne supérieure de la contrainte de cisaillement de chaque système de glissement est respectée. Cette idée est mise en oeuvre grâce à une optimisation multiniveau réalisée via la définition d'un ensemble de fonctions coûts non négatives pour chaque inclusion ainsi que pour l'agrégat entier.

  1)Pour l'analyse de dissociation en question, la fraction volumique de chaque composant de la région interlamellaire est requise. Elle est calculée par le biais des relations suivantes : En outre, le tenseur de rigidité de la phase amorphe, noté am C , est aussi nécessaire pour l'analyse de dissociation. Ce dernier est construit avec les deux constantes élastiques indépendantes de cette phase, en raison de sa symétrie isotrope. Dans ce but, le module élastique et le coefficient de Poisson de la phase amorphe, am E et am υ respectivement, sont déterminés sur la base de diverses études expérimentales [16-19]. A la température ambiante, les valeurs moyennes choisies pour am E et am υ sont 5 MPa et 0.4996, respectivement. Dans le régime caoutchoutique, en augmentant la température, am υ demeure quasiment constant puisqu'il est déjà très proche de sa valeur critique, i.e. 0.5. Pourtant, am E augmente légèrement avec la température, comme il est généralement admis pour les polymères caoutchoutique [20]. Donc, la variation suivante de et ip η sont connus, ils peuvent être remplacés dans les relations des approches du modèle étendu d'inclusion composite (ECIM) et dans la méthode de double inclusions (DIM), qui se trouvent ci-dessous.

	En remplaçant am υ et am E dans l'équation suivante, am C est construit :
	C	am ijkl	=	( 2 1 E am υ +	am	)	  	il jk δ δ	+	ik jl δ δ	+	am am 1 2 2 υ υ -	ij kl δ δ	  	(FR.4)
	Maintenant que il C , am C , am
				η	am	=	t	il	il 2 t -t	ip	,	η	ip	=	ip il t 2 t	(FR.2)
																am E avec la
	température est adoptée.
						am E =	5	293 T	(	MPa	)	(FR.3)

η

  Etant donné que dans le modèle néo-Hookéen, le module de cisaillement représente la rigidité du milieu, l'inégalité am ip µ µ < est évidemment satisfaite parce que la couche d'interphase est reconnue plus rigide que la phase amorphe centrale. De plus, dans certaines études incorporant l'hyperélasticité du domaine interlamellaire dans la formulation de grande déformation du polyéthylène, il est supposé que il 35MPa

	am am µ κ   	= =	32.3MPa 2.73GPa	,	il il µ κ   	= =	33.4 MPa 2.73GPa	,	ip ip µ κ     	= =	35.5MPa 2.74 GPa	(FR.12)

Acknowledgement

encouragement, inspiration and guidance throughout the research have been an invaluable asset for me. Thanks to the financial support from the French embassy in Tehran and the French government within the "Bourse de cotutelle" program, the co-direction of my thesis was made possible, which is gratefully acknowledged.

My sincere gratitude goes to Prof. Gregory C. Rutledge,

Appendix A -Euler angles and orientational averaging of fourth-order tensors

Euler angles are a classical means whereby the relative spatial orientation of a Cartesian coordinate system with respect to another Cartesian coordinate system is defined by way of a series of rotations. As illustrated in Figure A 1, the rotational transformation between two sets of axes is described by means of their relative Eulerian angles. In other words, every directional configuration in space is attainable by performing three successive rotations, each around only one axis; this is equivalent to saying that every rotation matrix can be represented as the product of three simpler rotation matrices. If the axes of the fixed reference frame is represented in small letters x,y,z and the rotated frame in capital letters X,Y,Z, then the line of nodes, N, is defined by the intersection of the xy-and XY-plane which is perpendicular to the z-and Z-axis. Among the quite a few possible conventions for defining the Euler angles, the one demonstrated in Figure A 1 is quite popular. According to this convention ϕ (or α ) is the angle between the x-axis and the node line N, θ (or β ) is the angle spanning the z-and Z-axis, and finally ψ or γ is the angle measured from the node line to the X-axis. Because different conventions may be defined and employed by different authors, any use of Euler angles must be accompanied by their definition. It should be underlined that Euler angles between two sets of coordinate axes are defined when the handedness of the two frames is identical.

Positive direction of Euler angles in the adopted convention is counterclockwise. To measure the Euler angles, one needs to apply counterclockwise rotations , , ϕ θ ψ on a hypothetical frame that is initially coincident with the fixed frame to obtain the new orientation. Angle θ ranges from 0 to π and the two other angles vary from 0 to 2π .

Accordingly, except for the singular states where the xy-and XY-plane coincide, Euler Equation Section (Next)

Appendix E-Finding the inverse of a fourth-order tensor with minor symmetry

In many micromechanical relationships, there appears the inverse of fourthorder tensors such as the Eshelby tensor or the stiffness tensor. In this appendix, a numerical technique is presented that helps to find the inverse of fourth-order tensors that possess minor symmetry is presented. This algorithm is founded on a specific representation of fourth-order tensors in a 9×9 matrix format as well as the definition of the fourth-order identity tensor.

The elements of an arbitrary fourth-order tensor ijkl S may be represented as follows in a 9×9 matrix format

It is underlined that although the notion of material symmetry in linear elasticity is traditionally linked to the crystal systems, the concept of elastic symmetries is generalizable to noncrystalline materials, such as biological tissues, geological samples or synthetic materials, which may behave anisotropically. Strictly speaking, in linear elasticity it has been demonstrated that the available symmetries for both crystalline and noncrystalline media are restricted to the eight aforementioned material symmetries [134,135]. On the other hand, save for the isotropic symmetry, visual recognition of the existing symmetry from the stiffness matrix is impossible unless it is expressed in its proper material reference frame. Several algorithms have been proposed in the literature whose objective is to determine the type of symmetry of any given stiffness matrix, regardless of the reference frame in which they have been expressed (e.g. [135,136]).

125

If the strains are confined to the linear elastic range, we will also have

where C and S stand for the effective stiffness and compliance of the heterogeneous medium. If the true stress and strain fields are respectively denoted by ij σ and ij ε , the following equalities are valid

In the above relationship, V is the integration volume, dΩ is the volume element, 0 i u is the displacement vector over ϒ . Similar equalities can be developed for ij σ and ij ε

Given that ij σ and ij ε are constant and that the rightmost sides of the two last relationships are identical, therefore

The interpretation of this last equality is that the volume average of the strain energy equals the inner product of the average stress and average strain tensors. Now the auxiliary variables 0 ij σ and 0 ij ε are defined as follows

Moreover, for the stress and strain fields of the heterogeneous medium we have

The following equalities can be derived from the combination of the two above relationships

and consequently

The second terms on the left hand side of both relationships of (H12) are positive as they can be decomposed into ( )(

)

)

σ σ which are always positive due to the positive definiteness of the stiffness and compliance matrices. Moreover, volume integration over the second terms of the right hand sides of the two last equalities vanishes since

Therefore, the following inequalities are readily evident

By substitution from (H5) and (H9) it is concluded that

The first inequality suggests that every individual element of the overall stiffness tensor estimated by Voigt method is greater than or equal to the corresponding element in the true effective stiffness tensor. Similarly, the second inequality suggests that every individual element of the true effective compliance is less than or equal to the corresponding element in the effective compliance estimated by Reuss method. Given that the compliance and stiffness tensors are mutually inverse, a better interpretation of these last two inequalities is that Voigt and Reuss estimates of the effective stiffness of any heterogeneous material yield the upper and lower bounds of the real effective stiffness.

It should be noted that in deriving the relationships for 8-chain and Gent back stress models, a new variable called "current chain stretch",

, is defined and exploited. Additionally, in deriving Gent back stress relationship, m J denotes the limit value of 3 J λ λ λ = + + -. This material property is correlated to N according to the following equation

which is easily established by equating the locking chain stretch with N .

Remarks on the inverse Langevin function

1 -

L

The distinguishing feature of 3-chain and 8-chain models is the inverse Langevin function which is indispensable to them whereas Gent model is free from this expression. Mathematically, Langevin function is defined as ( ) ( )

Langevin inverse, however, has no closed form solution and rational functions by Cohen [141], Puso [142] and Treloar [START_REF] Treloar | The Physics of Rubber Elasticity[END_REF] have been proposed, based on Padé approximation.

( Before comparing the abovementioned estimates, one should notice that since Langevin function is an odd function, its inverse is an odd function too, i.e. ( ) ( )

. Therefore, if an acceptable estimate for the positive half domain is found, the negative half can be determined simply by multiplying the positive solution by -1. Keeping this in mind, the following rational function is proposed for estimating the inverse Langevin function for positive y :

( ) ( )

To compare the four estimates considered here, their diagrams together with that produced by exact numerical solution of ( )

Abstract

Elastic characterization of the interphase layer in polyethylene is implemented by applying the relationships of two micromechanical approaches, "Extended Composite Inclusion Model" and "Double-Inclusion Method", to the Monte Carlo molecular simulation data for the interlamellar domain. The results of the two approaches match perfectly. The interphase stiffness lacks the common feature of positive definiteness, which indicates its mechanical instability. Comparison with experimental results endorses the proposed methodology. For the hyperelastic characterization of the interlamellar domain and the interphase layer, the proposed hybrid algorithm consists in applying the constitutive equations of an isotropic, compressible, hyperelastic continuum to the molecular dynamics simulation results of a polyethylene stack. Evolution of the interphase boundaries are introduced as auxiliary variables and the notion of minimizing a set of nonnegative objective functions is employed for parameter identification. The identified hyperelastic parameters for the interlamellar domain are in good agreement with the ones that have been estimated experimentally. Finally, the large, viscoplastic deformation of an aggregate of polyethylene is reexamined. The Gent model adopted for the back stress of the noncrystalline phase, correcting the projection tensor for the modified Taylor approach, and the idea of multilevel optimization are among the contributions made.

Keywords: semicrystalline polyethylene, interphase layer, micromechanics, mechanical characterization, uncertainty analysis, optimization of cost functions, Viscoplasticity, back stress, localization law.