Prof K Ravi-Chandar

Prof P Le Tallec

J Dequiedt

me, during the early stages of my thesis, to formulate the dynamical

Résumé

La localisation de la déformation dans un milieu ductile déformé est le mécanisme d'instabilité qui provoque la défaillance finale. Ce phénomène se produit sous chargement statique ainsi que dynamique. Les exemples sont présentés pour le cas de chargement quasi-statique dans Fig. 1.1a, montrant la localisation de la déformation sous la forme de bandes de cisaillement dans une bande de métal en tension et pour le chargement dynamique en L'instabilité du matériau dans des conditions de chargement quasi-statique et la connexion entre la localisation de la déformation dans les solides microstructurés et la perte de leur ellipticité macroscopique est étudiée dans le chapitre 2. Un exemple d'un tel phénomène est donné sur la Fig. 1.1e, décrivant la localisation macroscopique de la déformation sous la forme de kink-bands originaires du flambement microstructural de la fibre en sein de bois de balsa. L'approche de modélisation en milieu continu du phénomène de localisation étudie les conditions des lois constitutives conduisant à la perte d'ellipticité des équations gouvernantes, propriété qui permet des solutions d'équilibre discontinues. Les modèles de micro-mécanique et les théories d'homogénéisation nonlinéaires nous aident à comprendre les origines de ce comportement et l'on pense qu'une perte d'ellipticité macroscopique (homogénéisée) s'aboutit aux motifs de déformation localisés. Bien que cela soit le cas dans de nombreuses applications d'ingénierie, il pose une question intéressante: existe-t-il toujours un motif de déformation localisé apparaissant dans les solides quand le chargement critique est dépassé et qu'il y a une perte d'ellipticité macroscopique?

Dans un souci de simplicité relative et de trac ¸abilité analytique, nous répondons à cette question dans le cadre restrictif d'un solide multicouche, nonlinear (hyperélastique) en déformation plane et plus spécifiquement sous compression axiale le long de la direction de laminage. La clé de la réponse se trouve dans la solution post-bifurquée homogénéisée du problème, qui pour certains matériaux est supercritique (augmentation de la force et du déplacement), conduisant dans ces composites à des chemins d'équilibre post-bifurqués qui ne présent pas de localisation de déformation pour la contrainte macroscopique bien au-dessus de celui correspondant à la perte d'ellipticité.

L'instabilité du matériau en relation avec la striction avec les conditions de chargement dynamiques dans des plaques rapidement déformées, lorsque l'inertie est prise en compte, est présentée dans le chapitre 3. Nous suivons l'évolution temporelle des perturbations spatialement localisées et leurs interactions sur une plaque infiniment large sous tension biaxiale dont les lois constitutives présentent une perte d'ellipticité à des niveaux de déformation adéquats. L'évolution temporelle non linéaire d'une perturbation spatialement localisée est étudiée analytiquement et numériquement. La méthode analytique, basée sur la linéarisation, est utilisée pour définir la taille de la zone d'influence d'une perturbation ponctuelle et nous étudions sa dépendance à l'égard des lois constitutives et des conditions de chargement. Les calculs numériques montrent comment la zone de déformation localisée se propage et explique l'augmentation apparente de la ductilité dans les plaques minces par le temps requis par la zone de striction pour atteindre les bornes de la plaque. Les interactions de défauts donnent également une idée des modes de défaillance observés expérimentalement.

L'instabilité structurale sous chargement dynamique est le sujet du chapitre 4 qui étudie la localisation des motifs de déformation apparaissant dans des anneaux élastiques et élastoplastiques élastiques compressés électromagnétiquement. Contrairement à l'approche largement utilisée dans la littérature pertinente, qui repose sur la méthode d'analyse modale pour déterminer le mode propre le plus rapide de la structure-éloquent seulement pour les cas où la vitesse de la structure parfaite est significativement inférieure à l'onde caractéristique associée. Dans cette thèse, nous analysons la réponse temps-dépendante d'un anneau rapidement et hydrostatiquement compressé à des perturbations spatialement localisées de sa solution principale (radialement symétrique), afin de comprendre l'initiation des mécanismes de défaillance correspondants. Il est montré que pour de petites valeurs de la vitesse de chargement appliquée, la structure échoue par un mode global, alors que pour de grandes valeurs de la vitesse de chargement appliquée, la structure échoue par un mode localisé de déformation. Nous constatons également que les lois constitutives dépendant à la vitesse ne sont pas nécessaires pour modéliser les phénomènes associés, conformément aux observations expérimentales qui montrent qu'il n'y a pas d'augmentation de la ductilité dans les anneaux minces lorsque les taux de charge augmentent et l'importance des défauts statistiquement distribués qui expliquent l'absence d'une longueur d'onde dominante dans les modes de défaillance. nonlinear homogenization theories help us understand the origins of this behavior and it is thought that a loss of macroscopic (homogenized) ellipticity results in localized deformation patterns. Although this is the case in many engineering applications, it raises an interesting question: is there always a localized deformation pattern appearing in solids loosing macroscopic ellipticity when loaded past their critical state?

List of Figures

In the interest of relative simplicity and analytical tractability, we answer this question here in the restrictive framework of a layered, nonlinear (hyperelastic) solid in plane strain and more specifically under axial compression along the lamination direction. The key to the answer is found in the homogenized post-bifurcated solution of the problem, which for certain materials is supercritical (increasing force and displacement), leading to post-bifurcated equilibrium paths in these composites that show no localization of deformation for macroscopic strain well above the one corresponding to loss of ellipticity.

Material instability under dynamic loading conditions in connection to necking in rapidly strained plates, when inertia is taken into account, is presented in Chapter 3, We follow the time evolution of spatially localized perturbations and their interactions in biaxially strained thin plates whose constitutive laws exhibit loss of ellipticity at adequate strain levels. The nonlinear time evolution of a spatially localized perturbation is studied analytically and numerically. The analytical method, based on linearization, is used to define the size of the influence zone of a point-wise perturbation and we study its dependence on constitutive laws and loading conditions. Numerical calculations show how the localized deformation zone propagates and explain the apparent increase in ductility in thin plates by the time needed by the necking zone to reach the boundaries of the plate. Defect interactions also give an idea of the failure patterns observed experimentally.

Structural instability under dynamic loading conditions is the object of Chapter 4 that studies the localization of deformation patterns appearing in electromagnetically compressed thin elastic and elastoplastic rings. In contrast to the widely used approach in the relevant literature, which is based on the method of modal analysis to determine the structure's fastest growing eigenmodemeaningful only for cases where the velocity of the perfect structure is significantly lower than the associated characteristic wave propagation speeds, in this thesis we analyze the time-dependent response of a rapidly, hydrostatically compressed ring to spatially localized perturbations of its principal (radially symmetric) solution, in order to understand the initiation of the corresponding failure mechanisms. It is shown that for small values of the applied loading rate, the structure fails through a global mode, while for large values of the applied loading rate the structure fails by a localized mode of deformation. We also find that rate-sensitive constitutive laws are not necessary to model the associated phenomena, in accordance with experimental observations that show no ductility increase in thin rings when loading rates increase and the importance of statistically distributed defects that explain the absence of a dominant wavelength in failure patterns.

Chapter 2

Localization of deformation and loss of macroscopic ellipticity in microstructured solids

Introduction

Localization of deformation in finitely strained ductile solids is the instability mechanism leading to failure by rupture. The general principles were introduced for the study of this fascinating and important phenomenon in the context of continuum mechanics by [START_REF] Hadamard | Lec ¸ons sur la propagation des ondes et les équations de l'hydrodynamique[END_REF] and subsequently advanced in his spirit by [START_REF] Hill | Acceleration waves in solids[END_REF], [START_REF] Mandel | Conditions de stabilité et postulat de drucker[END_REF] and [START_REF] Rice | The localization of plastic deformation[END_REF]. The underlying mathematical concept in the continuum model is the loss of ellipticity in the governing equations, which allows for discontinuous strain solutions. With the advent of homogenization theories since the 1960's, a vast amount of work has been dedicated to the bridging of scales and understanding how micromechanical features in solids lead to their macroscopic (homogenized) loss of ellipticity at adequate levels of strain or stress. A plethora of applications for a wide range of solids has appeared in the literature, covering rubber elasticity, various types of composites (porous, fiber-reinforced, particlereinforced, cellular solids etc.), metal plasticity, granular media, rocks, just to name a few. Since the review of such a large and diverse body of work is unfortunately not possible, only key references relevant to the points made in the present chapter will be cited.

To avoid difficulties related to microstructure geometry and the identification of associated scale and representative volume, our attention is restricted to solids with a well defined scale, i.e. to architectured materials with periodic microsctructures. The role played by buckling at the microscopic scale, as the onset of instability mechanism leading to macroscopic localization of deformation in these materials has been established and subsequently analyzed by a long series of investigations. For the case of fiber reinforced composites, the connection between local buckling and global localization started with the work of [START_REF] Rosen | Mechanics of composite strengthening[END_REF], who recognized microbuckling as the onset of instability mechanism. Subsequent investigations of [START_REF] Budiansky | [END_REF], [START_REF] Budiansky | Compressive failure of fibre composites[END_REF], [START_REF] Kyriakides | On the compressive failure of fiber reinforced composites[END_REF], [START_REF] Vogler | On the initiation and growth of kink bands in fiber composites. Part II: analysis[END_REF] and many others showed, with progressively more sophisticated experiments and detailed modeling, how the buckling instability evolves into a localized deformation pattern (kink band formation) and studied in detail the characteristics of these bands. The same basic mechanism, i.e. buckling initiated at the microstructural level, has been recognized in materials science as the cause for localization of deformation in cellular solids (crushing zones) and the interested reader is referred to the comprehensive monograph by [START_REF] Gibson | Cellular solids: structure and properties[END_REF]. Detailed experimental and theoretical investigations followed in mechanics with a particular interest in studying the initiation and evolution towards localization of the deformation pattern in cellular solids by [START_REF] Papka | In-plane compressive response of crushing of honeycomb[END_REF], [START_REF] Papka | Experiments and full-scale numerical simulations of in-plane crushing of a honeycomb[END_REF], Papka and Kyriakides (1999a), [START_REF] Papka | Biaxial crushing of honeycombs -Part II: Analysis[END_REF] for two-dimensional microstructures and Jang et al.

(2010), [START_REF] Wilbert | Buckling and progressive crushing of laterally loaded honeycomb[END_REF] for three-dimensional microstructures and in establishing conditions where local or global buckling is the critical mechanism at the onset of failure by [START_REF] Triantafyllidis | Onset of failure in aluminum honeycombs under general in-plane loading[END_REF], Gong et al. (2005), Lopez-Jimenez and[START_REF] Lopez-Jimenez | Buckling of rectangular and hexagonal honeycomb under combined axial compression and transverse shear[END_REF].

Progressing in parallel, the nonlinear homogenization theories that appeared in mechanics first addressed questions on macroscopic response in plasticity, viscoelasticity and nonlinear elasticity with various microstructures (e.g. see [START_REF] Suquet | Analyse limite et homogénéisation[END_REF], [START_REF] Talbot | Variational principles for inhomogeneous nonlinear media[END_REF], Ponte Castañeda (1991)) and subsequently explored localization of deformation issues (e.g. see [START_REF] Kailasam | A general constitutive theory for linear and nonlinear particulate media with microstructure evolution[END_REF]Ponte Castañeda (1998), Lopez-Pamies and[START_REF] Lopez-Pamies | Second-order estimates for the macroscopic response and loss of ellipticity in porous rubbers at large deformations[END_REF]). For periodic solids the question asked was the possibility of detecting instabilities at the microscopic level from their homogenized properties, thus formally connecting buckling at the microscopic level to localization of deformation. For these composites it has been shown, initially for layered solids by [START_REF] Triantafyllidis | On the comparison between microscopic and macroscopic instability mechanisms in a class of fiber-reinforced composites[END_REF] and subsequently for the general three-dimensional periodic case by [START_REF] Geymonat | Homogenization of nonlinearly elastic materials, microscopic bifurcation and macroscopic loss of rank-one convexity[END_REF], that microstructural bifurcation phenomena (micro-buckling) is the mechanism responsible for macroscopic loss of ellipticity and that a long wavelength critical mode (based on Bloch wave analysis of the perfect infinite composite) coincides with the loss of ellipticity in its homogenized incremental moduli. Further work for porous elastomers by [START_REF] Michel | Microscopic and macroscopic instabilities in finitely strained porous elastomers[END_REF] and for particle reinforced elastomers by [START_REF] Michel | Microscopic and macroscopic instabilities in finitely strained fiber-reinforced elastomers[END_REF] has been done to connect local buckling to the macroscopic loss of ellipticity and compare periodic to random isotropic media with the same volume fractions.

Since loss of ellipticity is the property allowing for discontinuous equilibrium solutions, it is 2.1. INTRODUCTION thought (and supported by micromechanical calculations in most of the known-to the best of our knowledge-engineering applications), that a loss of macroscopic (homogenized) ellipticity results in a localized deformation pattern in the post-bifurcated regime. However, two questions arise: is there always a localized deformation appearing in the post-bifurcation of solids loosing macroscopic ellipticity and what are the necessary conditions in the homogenized response leading to localization?

In the interest of relative simplicity and analytical tractability, the present chapter answers these questions in the restrictive framework of an infinite, layered, nonlinear (hyperelastic) solid under plane strain loading conditions and more specifically under axial compression along the lamination direction. For this problem, one can find macroscopic loads where the homogenized moduli of the principal solution loose ellipticity (and since the solid has an energy density, the corresponding homogenized energy looses rank-one convexity). Moreover one can also ensure that the critical (i.e. corresponding to the lowest applied load) bifurcation eigenmode of the infinite solid is global (infinite wavelength eigenmode), a property that for this problem allows us to find a homogenized solution for the post-bifurcated equilibrium path. The answer to the localization question posed lies in the homogenized, initial post-bifurcation response of the perfect layered solid, as seen in Fig. 2.1; it will be shown that for a composite with a monotonically increasing force (and displacement) postbifurcation response (Λ 2 > 0, λ 2 > 0), no localized deformation solution develops in spite of a loss of ellipticity found at the macroscopic critical strain λ c .

Increasing force & displacement bifurcated path: ! 2 > 0 , ! 2 > 0 ! (force or stress)

" (displacement or strain) ! c ! c " principal solution
Max. force bifurcated path: ! 2 > 0 , ! 2 < 0 Max. displacement bifurcated path: ! 2 < 0 , ! 2 < 0 Figure 2.1: Different cases for the homogenized, initial post-bifurcation behavior of a perfect, nonlinear (hyperelastic) layered composite under plane strain loading conditions which is subjected to axial compression along its lamination direction. Stable paths are marked by solid lines and unstable ones by broken lines. For a composite with a monotonically increasing force (and displacement) post-bifurcation response (λ 2 > 0, Λ 2 > 0), no localized deformation solution develops in spite of a loss of ellipticity found at the macroscopic critical strain λ c . The presentation is organized as follows: The model of the perfect, laminated, periodic compos-ite is presented in Section 2.2; more specifically the bifurcation load and corresponding eigenmode and their nature (local or global) are discussed in Subsection 2.2.1. The exact solution for the postbifurcation equilibrium path corresponding to a global eigenmode is given in Subsection 2.2.2 while the asymptotic analysis of the homogenized post-bifurcation equilibrium path and its connection to the homogenized moduli of the composite is given in Subsection 2.2.3.

The results are presented in Section 2.3, starting with the choice of constitutive laws in Subsection 2.3.1 and continuing with the general homogenized solution for the post-bifurcated equilibrium path for an infinite, perfect, incompressible hyperelastic layered solid in Subsection 2.3.2. This model allows us to investigate all possible scenarios: cases under which this bifurcation involves a maximum displacement, a maximum force or a bifurcation occurring under increasing force and displacement (see Fig. 2.1). The case of neo-Hookean composites is presented in Subsection 2.3.3, where it is shown that they always have a stable, homogenized post-bifurcation response under increasing force and displacement. Composites with decreasing homogenized post-bifurcation force or displacement are given next in Subsection 2.3.4. The important question of how an adequately large, but finite-size volume of such a composite will behave is addressed next in Subsection 2.3.5.

It is shown, by means of introducing a small geometric imperfection at the middle of a large sample that under these "soft" boundary conditions, the monotonically increasing force (and displacement) composites will evolve towards a uniform shearing solution away from the macroscopic critical load and show no localization of deformation pattern past the critical load, in spite of a macroscopic loss of ellipticity; as expected the composites with the snap-through (i.e. maximum displacement) macroscopic response will evolve into a solution with a single strong localized deformation zone.

Concluding remarks are presented in Section 2.4. Finally some complementary material of interest is presented in the appendices: Details of the general bifurcation analysis of the infinite, perfect, rate-independent, layered composite in 2.A, the post-bifurcation equilibrium of the compressible neo-Hookean composite in 2.B and the influence of constitutive model choice on the critical load (i.e. comparison of the hyperelastic model with its deformation theory counterpart using the same uniaxial response) in 2.C.

Modeling

This section pertains to the modeling of the onset of bifurcation and post-bifurcation response of the axially compressed hyperelastic layered solid. Finding the critical load (i.e. lowest macroscopic compressive strain or stress) at the onset of bifurcation and the corresponding eigenmode is presented in Subsection 2.2.1. The exact solution of the post-bifurcation problem for the case of a 2.2. MODELING global critical eigenmode is given in Subsection 2.2.2, followed by the corresponding asymptotic solution of this problem near the critical load in Subsection 2.2.3.

Bifurcation of a layered solid in plane strain: local vs. global critical mode

The presentation starts with the solution for the plane strain bifurcation problem for an infinite, perfectly periodic, layered, hyperelastic solid. The composite is subjected to uniaxial compression along the fiber direction characterized by a monotonically increasing "load parameter" λ ≥ 0, which designates the absolute value of the applied macroscopic strain under displacement control (or its corresponding work-conjugate stress Λ ≥ 0, when force is controlled). The goal is to find the lowest critical load λ c (or Λ c ) and corresponding eigenmode as the load parameter increases away from λ = 0 (or Λ = 0), which is the stress-free, reference configuration of the solid. Without loss of generality, it is assumed in this section that the composite is loaded under displacement control. The infinite, weightless, perfectly periodic solid is composed of a self-repeating sequence of two layers, as depicted in Fig. 2.2a. The two layers f (fiber) and m (matrix) have initial thickness H f and H m in the stress-free reference configuration. The solid is deformed under finite, plane strain conditions with perfect bonding between layers which guarantees traction and displacement continuity across each interface for all possible deformations. A full Lagrangian formulation of the problem is adopted with respect to a fixed Cartesian coordinate system X 1 -X 2 , where X 1 is the lamination direction. The procedure for determining the onset of the first bifurcation hinges on finding a non-trivial solution to the difference between principal and bifurcated equilibrium solutions at the onset of a bifurcation:

local mode λ c global mode λ H λ λ H f H m X 1 X 2 unit cell Reference configuration (stress-free) (a) (b)
∆Π ji,j = 0; X 1 ∈ R, X 2 ∈ [0, H], 1 (2.2.1)
and interface conditions:

[ [∆Π 2i ] ] = 0, [[∆u i ]] = 0; X 1 ∈ R, X 2 ∈ {0, H m , H}, (2.2.2)
where ∆Π is the difference in the first Piola-Kirchhoff stress tensors of the principal and bifurcated equilibrium solutions, ∆u is the corresponding displacement field difference and

H = H m + H f
is the initial thickness of the unit cell. Moreover, [[g]] denotes a difference in the values of any field quantity g when evaluated on both sides of an interface.

The constitutive response of the hyperelastic composite takes the form:

Π ji = ∂W F ij , (2.2.3)
where W (F) is the strain energy density of the fiber or matrix layer and F is the corresponding deformation gradient. Consequently, ∆Π from (2.2.1) can be expressed in terms of ∆u by:

∆Π ji = L ijkl ∆F kl ; L ijkl = ∂ 2 W ∂F ij ∂F kl , ∆F kl = ∆u k,l .
(2.2.4)

The fourth rank tensor L is termed the "incremental moduli tensor" and is a function of the position X and the load parameter λ. The above formulation pertains to the case of a compressible solid.

The slightly simpler formulation for the incompressible case has already been presented by [START_REF] Triantafyllidis | On the comparison between microscopic and macroscopic instability mechanisms in a class of fiber-reinforced composites[END_REF]. The compressible version for the onset of bifurcation in the axially compressed, hyperelastic layered solid under plane strain is given in [START_REF] Geymonat | Homogenization of nonlinearly elastic materials, microscopic bifurcation and macroscopic loss of rank-one convexity[END_REF], where the interested reader can find the complete derivations. However, for reasons of completeness of the presentation, the main results are outlined in this section and detailed derivations are given in 2.A.

The critical load λ c , corresponding to the onset of the first (as load parameter increases away from zero) bifurcation is found to be:

λ c := inf ω 1 H>0 λ(ω 1 H), (2.2.5)
where λ(ω 1 H) is a function of the dimensionless wavenumber ω 1 H of the eigenmode along the X 1 direction, as defined in (2.A.14), and depends on the two layer thicknesses and their incremental 2.2. MODELING moduli L m (λ) and L f (λ). When the wavelength L 1 of the critical mode is commensurate with the layer thickness H, i.e. (L 1 /H) c = 2π/(ω 1 H) c = 0, the critical mode is termed "local", and when this is not the case, i.e. (L 1 /H) c = 2π/(ω 1 H) c → ∞ ⇐⇒ (ω 1 H) c → 0, the critical mode is termed "global" (see Fig. 2.2b). The case (ω 1 H) c = 0 corresponds to an X 1 -independent solution, which is excluded by the local rank-one convexity condition for each layer, as shown in 2.A, thus explaining the use of infimum in (2.2.5).

For the present chapter it is important to ensure that the bifurcation occurring at λ c is global, i.e.

(ω 1 H) c → 0. In this case it can be shown that the critical load λ c = λ H , where λ H is the lowest load corresponding to the first loss of rank-one convexity of the homogenized moduli L H (λ) of the composite defined by:

λ H := min n =1 λ > 0 det L H ijkl (λ)n j n l = 0 , (2.2.6)
where the expressions for the components of the homogenized moduli L H (λ) tensor are given in (2.A.21).

For all the layered composites considered here it is shown that (2.2.6) is satisfied along the critical direction (n 1 , n 2 ) c = (1, 0) and hence:

L H 2121 (λ c ) = 0, (2.2.7)
which implies a vanishing shear stiffness of the composite at that load, perpendicular to the lamination direction. according to which all layers rotate by the same angle, while each layer experiences a uniform state of strain (and stress). All fiber layers share the same strain F f and all the matrix layers share

Post

F m , but F f = F m .
More specifically, the idealized post-buckling equilibrium is a periodic solution with a unit cell depicted in Figure 2.4. This post-bifurcated equilibrium path is found as follows: From the kinematics of deformation one has the following relations 2 : From equilibrium we obtain the following relations for Π ji , the average first Piola-Kirchhoff stresses which are work-conjugate to F ij :

f= -Λ δ = λ (a) (b) (c) δ = λ f= -Λ f= -Λ f= -Λ
F 11 = F f 11 = F m 11 = 1 -λ, F 21 = F f 21 = F m 21 = γ, F 12 = ζ f F f 12 + ζ m F m 12 = 0, F 22 = ζ f F f 22 + ζ m F m 22 , (
Π 11 = ζ f Π f 11 + ζ m Π m 11 = -Λ, Π 12 = ζ f Π f 12 + ζ m Π m 12 = 0, Π 21 = Π f 21 = Π m 21 , Π 22 = Π f 22 = Π m 22 = 0.
(2.2.9)

The above relations reflect the fact that the bifurcated configuration is subjected to an average compressive axial stress ( Π 11 = -Λ < 0) along the X 1 direction, zero average shear stress ( Π 12 = 0) and zero lateral normal stress ( Π 22 = 0). To complete the above system of equations (2.2.8) and (2.2.9), one has to add the constitutive response of the fiber and matrix, given by (2.2.3). A closed-form analytical solution is possible only for simple constitutive laws, as detailed in Section 2.3.

Post-bifurcation equilibrium for global critical mode: asymptotics near critical point

The bifurcated equilibrium path can be obtained from its homogenized energy density, defined by:

W H (λ, γ) = W := ζ f W f + ζ m W m .
(2.2.10)

By extremizing the homogenized potential energy W -Λλ with respect to the displacement parameters λ, γ introduced in Section 2.2.2, we show that the onset of bifurcation coincides with the vanishing of the homogenized moduli L H 2121 (λ c ) in (2.2.7). Indeed, a slightly stronger result is available by showing that along the principal equilibrium path (γ = 0):

W 0 ,γγ = L H 2121 (λ), (2.2.11)
where the 0 superscript or 0 subscript denotes evaluation on the principal equilibrium path. Indeed from (2.2.10) and the definitions of the first Piola-Kirchhoff stress and the incremental moduli in (2.2.3) and (2.2.4) one obtains:

W 0 ,γγ = L 0 ijkl ∂F ij ∂γ γ=0 ∂F kl ∂γ γ=0 + Π 0 ji ∂ 2 F ij ∂γ 2 γ=0 . (2.2.12) Notice first that in view of the kinematics (2.2.8) 1 , (2.2.8) 2 that ∂ 2 F 11 /∂γ 2 γ=0 = ∂ 2 F 21 /∂γ 2 γ=0 = 0. Moreover the absence of lateral stresses ∂W m /∂F m 22 = ∂W f /∂F f 22 = 0 in (2.
2.9) 4 and the vanishing of the average shear F 12 according to (2.2.8) 3 lead to the vanishing of the Π 0 ji ∂ 2 F ij /∂γ 2 γ=0 term of W 0 ,γγ in (2.2.12). Consequently and since from kinematics ∂F 11 /∂γ = 0, ∂F 21 /∂γ = 1 according to (2.2.8) 1 and(2.2.8) 2 respectively, (2.2.12) can be rewritten as:

W 0 ,γγ = L 0 2121 + L 0 2112 + L 0 1221 ∂F 12 ∂γ γ=0 + L 0 1212 ∂F 12 ∂γ 2 γ=0 .
(2.2.13)

In the derivation of (2.2.13) use is made of the fact that in view of the orthotropy of the principal solution L 0 1211 = L 0 1222 = L 0 2111 = L 0 2122 = 0 (plus all the principal symmetry related moduli) vanish, plus the fact (∂F 22 /∂γ) γ=0 = 0, which follows from the principal solution's orthotropy and the vanishing of the lateral stress Π m 22 = Π f 22 = 0. The last remaining ingredient to prove (2.2.11) is the calculation of (∂F 12 /∂γ) γ=0 for the matrix and fiber layers respectively. Indeed from taking the γ-derivatives of (2.2.8) 3 and (2.2.9) 3 one has:

ζ m ∂F m 12 ∂γ γ=0 + ζ f ∂F f 12 ∂γ γ=0 = 0, L 0m 1212 ∂F m 12 ∂γ γ=0 + L 0m 1221 = L 0f 1212 ∂F f 12 ∂γ γ=0 + L 0f 1221 .
(2.2.14)

From the above linear system one can calculate (∂F m 12 /∂γ) γ=0 and ∂F f 12 /∂γ γ=0 , which upon substitution into (2.2.13) and recalling the definition of L H 2121 (λ) in (2.A.21) 4 , gives (2.2.11). For the case of a hyperelastic layered composite of arbitrary energy density and volume fractions where a closed-form analytical solution is not possible, the initial stability of the bifurcated equilibrium path near the critical load can be found asymptotically. The bifurcation amplitude parameter of the system is the shear γ (principal solution γ = 0). Due to the symmetry of the problem, one has near the critical point the following asymptotic expressions for the applied average stress (Λ) or average strain (λ):

Λ = Λ c + Λ 2 γ 2 2 + O(γ 4 ), λ = λ c + λ 2 γ 2 2 + O(γ 4 ). (2.2.15)
The goal of the asymptotic analysis is to obtain Λ 2 and λ 2 as functions of geometry and material properties of the composite. Consequently, according to the general theory of stability of an elastic system with a simple eigenmode at criticality (the homogenized perfect composite has a finite energy given in (2.2.10)), for displacement control, the stability of the bifurcated path near the critical point requires λ 2 > 0 (or equivalently Λ 2 > 0) for force control.

RESULTS

Results

Following the presentation of the constitutive laws for the composite in Subsection 2.3.1, we proceed to the general formulation of the post-bifurcated equilibrium path for infinite, incompressible, hyperelastic layered composites in Subsection 2.3.2. It is subsequently shown in Subsection 2.3.3, that all perfect neo-Hookean composites irrespectively of fiber-to-matrix thickness ratio, have an initial post-bifurcation response with increasing load and displacement. We then present in Subsection 2.3.4 a more general layered composite, consisting of an equal thickness neo-Hookean layer and a softer nonlinear matrix, that can exhibit all possible post-critical responses. It is worth noticing that all these different post-bifurcation responses can be achieved with composites made of locally stable, i.e. strongly elliptic layers. Numerical (FEM) calculations for boundary value problems

showing the absence or presence of localized deformation zones in these composites, in accordance with the predicted initial post-bifurcation response of their perfect counterparts, are presented in Subsection 2.3.5.

Constitutive laws

Two different versions of a rank-one convex, plane strain, isotropic, hyperelastic constitutive law are employed. The incompressible version-used for its analytical tractability-of any isotropic hyperelastic solid subjected to plane strain can be written in terms of a one variable scalar function

g(z): 3 W (F) = µ 2 g(z); z := I 1 -2, I 2 = 1, g(0) = 0, g (0) = 1, (2.3.1)
where µ > 0 is the initial shear modulus of the material and I i (i = 1, 2) are the two invariants of the right Cauchy-Green tensor C, (related to the deformation gradient tensor F by C ij = F ki F kj ), namely:

I 1 = tr C, I 2 = 1 2 (tr C) 2 -tr C 2 = det(C). (2.3.2)
Rank-one convexity is guaranteed when g > 0 and g ≥ 0 (rank-one convexity requires a weaker condition:

g + 2zg > 0, since z ≥ 0).
A compressible version of (2.3.1) will also be used:

W (F) = µ 2 g(y) + κ 2 (I 1/2 2 -1) 2 ; y := I 1 -2I 1/2 2 , (2.3.3)
where µ and κ µ are the initial shear and bulk moduli of the solid. This material's rank-one convexity follows from its polyconvexity (see [START_REF] Ball | Convexity conditions and existence theorems in nonlinear elasticity[END_REF]) which is guaranteed, since y ≥ 0, by: g > 0 and g ≥ 0.

Post-bifurcation equilibrium: general setting

As previously mentioned, a closed-form analytical solution of the post-bifurcation equilibrium path described in Subsection 2.2.2. is not in general possible, thus motivating the need for the asymptotic analysis presented in Subsection 2.2.3. However, for the special case of incompressibility or of the case of compressible neo-Hookean composites one can find analytically tractable expressions for the post-bifurcation equilibrium paths.

For the incompressible composite in (2.3.1), the constitutive law of (2.2.3) a takes the form:

Π ji = ∂W F ij -p ∂ ∂F ij (det F -1) , det F = (I 2 ) 1/2 = F 11 F 22 -F 12 F 21 = 1, (2.3.4)
where p is the Lagrange multiplier associated to the incompressibility constraint (the undefined part of the hydrostatic pressure).

Recalling from the kinematics in (2.2.8) 1 , (2.2.8) 2 that F 11 = 1λ and F 21 = γ one obtains that:

Π 11 = µ g 1 -λ - (1 + γ F 12 ) 2 (1 -λ) 3 , Π 22 = 0 ⇒ p = µ g (1 + γ F 12 ) (1 -λ) 2 , Π 21 = µ g F 12 + γ (1 + γ F 12 ) (1 -λ) 2 , Π 12 = µ g γ + F 12 (1 + γ F 12 ) (1 -λ) 2 .
(2.3.5)

From shear traction continuity Π f 21 = Π m 21 from (2.2.9) 3 and the kinematic constraint that F 12 = 0 from (2.2.8) 3 one has the following system for the two unknowns F f 12 and F m 12 :

1 + γ 2 (1 -λ) 2 m f F f 12 -m m F m 12 + γ (1 -λ) 2 (m f -m m ) = 0, ζ f F f 12 + ζ m F m 12 = 0, m f := µ f g f , m m := µ m g m .
(2.3.6)

RESULTS

Solving the above linear system for F f 12 , F m 12 yields:

F f 12 = γ (1 -λ) 2 + γ 2 (m m -m f ) ζ m ζ f m m + ζ m m f , F m 12 = γ (1 -λ) 2 + γ 2 (m f -m m ) ζ f ζ f m m + ζ m m f .
(2.3.7)

The λ-γ relationship along the bifurcated equilibrium path can now be obtained by combining Π 12 = 0 from (2.2.9) 2 with the expressions for Π f, m 12 and F f, m 12 in (2.3.5) 4 and (2.3.7):

Π 12 = γ      m G - (∆m) 2 m H 1 (1 -λ) 2 + γ 2 2      , ∆m := m f -m m , m G := ζ f m f + ζ m m m , m H := m f ζ f + m m ζ m .
(2.3.8)

Thus the sought-after λ-γ relationship takes the form:

(1 -λ) 2 + γ 2 = ∆m (m G m H ) 1/2 , (2.3.9)
where without loss of generality we have tacitly assumed ∆m > 0 i.e. m f > m m .

One important remark about (2.3.9) is in order at this point, namely that this is an implicit equation for λ and γ, given the fact that m f and m m are functions of I f 1 and I m 1 which in turn depend on λ and γ by: (2.3.10) where the expressions for F 12 in the fiber and matrix layers are given in (2.3.7). The critical axial strain λ c is found from (2.3.9) for γ = 0.

I 1 = F ij F ij = (1 -λ) 2 + 1 + γ F 12 1 -λ 2 + (F 12 ) 2 + γ 2 ,

Neo-Hookean composites

As it turns out, the simplest case of a neo-Hookean composite is always stable under either force of displacement control with monotonically increasing forces and displacements as functions of the shear γ. Moreover, a closed-form solution is possible in this case, since m f = µ f and m m = µ m are now constants: (2.3.11) where it is tacitly assumed that the shear strain satisfies:

(1 -λ) 2 + γ 2 = ∆µ (µ G µ H ) 1/2 =⇒ λ = 1 - ∆µ (µ G µ H ) 1/2 -γ 2 1/2 , ∆µ := µ f -µ m , µ G := ζ f µ f + ζ m µ m , µ H := µ f ζ f + µ m ζ m ,
0 ≤ γ < (∆µ) 1/2 /(µ G µ H ) 1/4 .
Using (2.3.11), we calculate Λ, the absolute value of the average compressive stress Π 11 : The absolute value of the dimensionless axial stress Λ/µ G vs. its work-conjugate strain λ is depicted in Fig. 2.5, which shows the stable equilibrium paths of four different volume fraction neo-Hookean composites for two different values of fiber-to-matrix stiffness ratios µ f /µ m = 10, 100. As expected, increasing the fiber-to-matrix stiffness ratio, leads to a sharp decrease in the slope of the bifurcated equilibrium path, in view of an increasingly softening matrix response in the post-bifurcated path. From (2.3.11) one obtains, by setting γ = 0, the following simple expression for the critical load λ c under displacement control for the neo-Hookean composite:

Λ = -Π 11 = µ G (1 -λ) 3 1 -(1 -λ) 2 + γ 2 2 = µ G        1 - (∆µ) 2 (µ G µ H ) ∆µ (µ G µ H ) 1/2 -γ 2 3/2        . (2.3.12) µ f /µ m = 10 µ f /µ m = 100
λ c = 1 - (∆µ) 2 µ G µ H 1/4 . (2.3.13)
The influence of fiber volume fraction and stiffness contrast on the critical strain λ c is depicted in Fig. 2.6. The reference critical strain λ ref used for calculating the influence of fiber volume fraction under fixed stiffness contrast corresponds to ζ f = 0.5 (and changes according to the stiffness contrast), while the reference critical strain λ ref used for calculating the influence of stiffness contrast under fixed fiber volume fraction corresponds to µ f /µ m = 2 (and changes according to the 2.3. RESULTS fiber volume fraction). As expected the critical strain increases monotonically with increasing fiber volume fraction, while for a fixed fiber volume fraction, the critical strain decreases monotonically with increasing fiber-to-matrix stiffness ratio. As expected the critical stress increases monotonically with increasing fiber volume fraction, while for a fixed fiber volume fraction, the critical stress after a steep initial decrease with increasing fiber-to-matrix stiffness ratio, reaches a plateau after about µ f /µ m = 10 but not in a monotonic fashion, as seen in particular for ζ f = 0.95. Since an analytic form is available for the entire homogenized, post-bifurcated equilibrium path in a perfect neo-Hookean composite according to (2.3.11) and (2.3.12), one can find the stability of this equilibrium path by checking the positive definiteness of the homogenized energy. For algebraic simplicity, we calculate here the initial curvatures at criticality λ 2 and Λ 2 , whose positivity implies stability, according to the general theory, at least in a neighborhood of the critical point.

2 10 100 1000 0 1 Μ f Μ m Λ c Λ ref
At criticality, the strain curvature λ 2 of the neo-Hookean composite's bifurcated equilibrium path is found from (2.3.11):

λ 2 = 1 1 -λ c =⇒ λ 2 λ c = 1 λ c (1 -λ c ) (2.3.14)
The corresponding curvature Λ 2 is obtained by evaluating d 2 Λ/dγ 2 γ=0 from (2.3.12):

Λ 2 = 3 µ G 1 -λ c 1 (1 -λ c ) 4 -1 =⇒ Λ 2 Λ c = 3 (1 -λ c ) 2 .
(2.3.15)

We have thus shown that for the perfect, incompressible neo-Hookean composite, its homogenized post-bifurcation path has monotonically increasing displacements and forces; it is thus stable near the critical point and even beyond.

This result is not surprising; in contrast to classical fiber reinforced composites where the much stiffer fiber is linearly elastic while the soft matrix has a low hardening which favors shearing of the matrix and unloading of the fiber leading to a decreasing displacement (snap-back), here each layer of the composite stiffens at the same rate, resulting in an increasing force and displacement at the bifurcated equilibrium path.

It is worth checking if this strong post-bifurcation stability result found for arbitrary neo-Hookean composites is influenced by compressibility; it turns out that it is not and the corresponding calculations are in 2.B.

RESULTS

Rank-one convex composites exhibiting all possible cases of post-bifurcation response

We have proved that the homogenized, post-bifurcation equilibrium path of a perfect layered composite consisting of two neo-Hookean layers of arbitrary shear moduli and volume fractions is stable (near the critical load), under either force (Λ 2 > 0) or displacement (λ 2 > 0) control. An even stronger result holds: the homogenized, post-bifurcation equilibrium paths of the perfect neo-Hookean composite have a monotonically increasing macroscopic stress and strain as the bifurcation amplitude γ increases (dΛ/dγ > 0 and dλ/dγ > 0), as one can easily show from (2.3.11) and

(2.3.12).

The question that naturally arises is whether a different choice of a locally stable (i.e. rank-one convex) constitutive laws, i.e. a nonlinear g(z) in (2.3.1), can lead to an unstable post-bifurcation

response. The answer to this question is affirmative and we present below composites that exhibit unstable post-bifurcated equilibrium paths even under displacement control, i.e. λ 2 < 0.

For analytical tractability, we consider composites of two equal thickness layers (ζ f = ζ m = 0.5) but different energy densities, defined by (2.3.1):

g f (z) = 1 2 z, g m (z) = α 2 z n ; 1 ≥ α > 0, 1 ≥ n > 0.5. (2.3.16)
Since the material is isotropic, incompressible and under plane strain conditions, it is best described by its response in simple shear γ. Recalling that under simple shear (γ) 2 = z (:= I 1 -2) and that its shear stress τ = dW/dγ, we record in Fig. 2.8 the response of the fiber and matrix layers.

The response of the neo-Hookean fiber is linear (µ f = 1) and always stiffer than the response of the matrix, which exhibits softening under increasing applied strain. The singularity of the matrix material at the origin can be removed 4 and the shear strain γ m where the matrix response starts differing from the fiber's is adjusted by the parameter α as seen in Fig. 2.8. The equal layer thickness neo-Hookean composite is recovered for a matrix with n = 1 and α = µ m /µ f .

For the case ζ f = ζ m = 0.5 one obtains from (2.3.7):

F f 12 = -γ, F m 12 = γ, (2.3.17)
thus yielding with the help of (2.3.10) the following expressions for the strain invariant I 1 in each 4 We set g f (z) = gm(z) for 0 ≤ z ≤ zm and gm(z) = 0.5αz n for z ≥ zm, where zm = (γm) 2 with γm the shear strain that marks the onset of the matrix nonlinear response. Continuity of shear stress at γm dictates that:

nα = z (1-n) m
.

Note that the matrix material is always strongly rank-one convex since d 2 W/dγ 2 > 0.

layer:

I f 1 = (1 -λ) 2 + 1 (1 -λ) 2 + 2γ 2 1 - 1 (1 -λ) 2 + γ 4 (1 -λ) 2 = z f + 2, I m 1 = (1 -λ) 2 + 1 (1 -λ) 2 + 2γ 2 1 + 1 (1 -λ) 2 + γ 4 (1 -λ) 2 = z m + 2.
(2.3.18)

Consequently the λ-γ relationship along the bifurcated equilibrium path in (2.3.9) takes the form (after recalling also the definitions for m f , m m in (2.3.6) 3 ):

(1 -λ) 2 + γ 2 = 1 -nα(z m ) n-1 1 + nα(z m ) n-1 . (2.3.19)
The corresponding expression for the macroscopic stress along the bifurcated equilibrium path is found to be:

Λ = -Π 11 = 1 + nα(z m ) n-1 2(1 -λ) 3 1 -(1 -λ) 2 + γ 2 2 . (2.3.20)
Solving equation (2.3.19) for γ = 0 yields the following implicit equation for the critical strain λ c :

nα (1 -λ c ) 2 + 1 (1 -λ c ) 2 -2 (n-1) = 1 -(1 -λ c ) 2 1 + (1 -λ c ) 2 , (2.3.21)
where the bifurcation occurs at the lowest positive root λ c of the implicit equation (2.3.21).

Of interest is the initial post-bifurcation behavior of this composite, i.e. the curvatures λ 2 and Λ 2 of the bifurcated equilibrium path at critical load. To this end, by taking into account (2.3.18), we differentiate (2.3.19) with respect to γ at γ = 0 to find:

λ 2 = (1 -λ c ) 1 -(1 -λ c ) 2    1 -(1 -λ c ) 2 -(1 -n) 1 + (1 -λ c ) 2 2 (1 -λ c ) 2 + 1 -n 2 1 + (1 -λ c ) 2 2    . (2.3.22)
Initial post-bifurcation stability of this composite under displacement control requires λ 2 > 0, i.e.

that the numerator in the above expression for λ 2 be positive:

1 -(1 -λ c ) 2 -(1 -n) 1 + (1 -λ c ) 2 2 > 0. (2.3.23)
To find the initial stability under force control, by taking into account (2.3.21), we differentiate

(2.3.20) with respect to γ at γ = 0 to find Λ 2 :

Λ 2 = (1 -λ c ) -3 1 -(1 -λ c ) 2    3 1 -(1 -λ c ) 2 2 -2(1 -n) 2 -(1 -λ c ) 2 1 + (1 -λ c ) 2 2 (1 -λ c ) 2 + 1 -n 2 1 + (1 -λ c ) 2 2    .
(2.3.24) Initial post-bifurcation stability of this composite under force control requires Λ 2 > 0, i.e. that the numerator in the above expression for Λ 2 be positive:

3 1 -(1 -λ c ) 2 2 -2(1 -n) 2 -(1 -λ c ) 2 1 + (1 -λ c ) 2 2 > 0. (2.3.25)
The range of matrix constitutive parameters (α, n, see definition in (2.3.16)), for which the homogenized, post-bifurcated solution is initially stable under either force or displacement control (Λ 2 > 0, λ 2 > 0, red-shaded area), is initially stable only under displacement control (Λ 2 < 0, λ 2 > 0, green-shaded area) and initially unstable under either force or displacement control (Λ 2 < 0, λ 2 < 0, non-shaded area), is given in Fig. 2.8. The same figure also records contours of equal critical strain λ c as a function of the matrix constitutive parameters α, n. Notice that for large critical strains (approximately λ c > 0.5) the composite is initially stable, while for small strains (depending on the value of n) the composite exhibits a snap-back (λ 2 < 0). A word of caution: for the post-bifurcation results to be meaningful, one should select matrix parameters α, n for which the critical mode is global in nature, i.e. exclude composites with matrix parameters in the blue-shaded area of Fig. 2.8 (see discussion in Section 2.2). The results obtained in Fig. 2.8 are important in selecting composites for the numerical calculations of the boundary value problems reported in the next section.

Boundary value problem calculations

Thus far we have established the homogenized post-bifurcation behavior of an infinite, perfect, fiberreinforced, hyperelastic composite compressed axially along the direction of its fibers. We want to find out the response in the bulk of such a composite (i.e. away from boundaries) by solving an appropriate boundary value problem, in order to establish whether a localized pattern of deformation will actually emerge or not, once the critical load is approached, given that the energy density of the homogenized principal solution looses its rank-one convexity once a critical macroscopic strain (λ c )

or stress (Λ c ) has been reached. We seek the equilibrium path of a "realistic" such composite, i.e.

one with the inevitable small imperfections and finite boundaries; in the interest of avoiding edge effects, an adequately large sample should be considered and the perturbation should vanish on the sample edges. No analytical solution being possible (to the best of our knowledge) for a finite-size, imperfect composite, a numerical one is sought based on FEM discretization of the corresponding boundary value problem. It should be mentioned at this point that there is a vast amount of work dedicated to the problem of multiscale calculations, i.e. how to find the response of a microstructured solid which is prone to localization by selecting an appropriate representative volume element, consistent boundary conditions and an efficient computational strategy. Our main concern for the FEM calculations is the avoidance of boundary strain concentrations and this is achieved by the procedure described below. However, the reader interested in this aspect of the problem is referred to the recent article by [START_REF] Coenen | A multiscale approach to bridge microscale damage and macroscale failure: a nested computational homogenization-localization framework[END_REF], who propose an interesting computational strategy for arbitrary microstructures and also present a very comprehensive literature review.

To this end we seek the solution of a finite-size, initially square segment of the composite (of dimensions L × L), subjected to a macroscopic (average) deformation gradient with The reference configuration of the sample has a small geometric imperfection with respect to its perfect counterpart, which is necessary to trigger the bifurcated solution. Accordingly the following initial imperfection which is designed to have a maximum amplitude at the middle of the finite-size composite and vanish at its boundaries:

F 11 = 1 -λ,
∆X 1 = 0, ∆X 2 = ξ sin(πX 2 /L) arctan[β(X 1 /L -1/2)], (2.3.26)
where the parameter ξ = 10 -3 controls the maximum X 2 -amplitude of the perturbation and the parameter β = 8 controls its X 1 -range. The origin of the coordinate system is taken at the bottom left corner of the domain. The control parameter is the macroscopic axial strain λ and a Riks continuation method is employed to bypass the limit loads which are associated to snapbacks.

The FEM calculations presented correspond to three different composites, described in Subsection 2.3.4; the material parameters of their matrix layers are indicated by a dot in Fig. 2.8, which

gives the initial post-bifurcation of the homogenized response for the corresponding perfect composite (the n = 1 axis corresponds to the fully neo-Hookean composite with α = µ m /µ f , as described in Subsection 2.3.3).

The response of a neo-Hookean composite with ζ f = ζ m = 0.5 and α = µ m /µ f = 0.22, whose perfect configuration exhibits a stable post-bifurcation behavior (Λ 2 > 0, λ 2 > 0), is shown in Fig. 2.9. The composite has two equal thickness neo-Hookean (n = 1) layers (see point A in Fig. 2.8). On the left is the macroscopic stress-strain response of the homogenized perfect composite (in red) and of its finite-size imperfect counterpart, based on FEM computations using 40 unit cells (in black). On the middle and right are depicted the deformed configurations at points a and b,

showing also contours of the Lagrangian shear strain component E 12 . The macroscopic stress-strain curve of the imperfect composite is almost indistinguishable from its perfect counterpart. When the critical strain λ c is reached at point a, the imperfect composite shows a higher strain zone in the middle, as seen in the middle picture of Fig. 2.9. Upon further increase of the applied macroscopic strain, the composite reaches at point b a uniform solution corresponding to the bifurcated state of its perfect counterpart, as shown in the right picture of Fig. 2.9. In contrast to what is known in typical engineering composites, the present example shows a case where in spite of loosing its macroscopic rank-one convexity, the composite exhibits no localized zone of deformation in its post-critical response. Similarly to the neo-Hookean case presented in Fig. 2.9, the macroscopic stressstrain curve of the imperfect composite is almost indistinguishable from its perfect counterpart.

When the critical strain λ c is reached at point a, the imperfect composite is covered by alternating higher and lower strain zones (a mix of principal and bifurcated solutions of the corresponding infinite, perfect composite), as shown in the middle picture of Fig. 2.10. Upon further increase of All the above calculations pertain to locally stable (i.e. locally rank-one convex) hyperelastic composites that share the same feature: a critical (i.e. occurring at lowest applied macroscopic load)

2.4. CONCLUSION long wavelength bifurcation mode which corresponds to a loss of ellipticity in the homogenized principal solution. Our calculations show that a localization of deformation in these composites does not always appear in the neighborhood of critical load and beyond; the result depends on the post-bifurcation behavior of the homogenized, perfect composite. In contrast to the widely studied cases of elastoplastic composites, with high stiffness contrast between fiber and matrix that exhibit kink band solutions, we have shown here the existence of composites that have stable, homogenized post-bifurcated solutions with increasing macroscopic stresses and strains. These composites do not exhibit localized deformation post-bifurcated solutions, in spite of the fact that their homogenized energy looses its rank-one convexity as the applied loading increases.

Conclusion

Localization of deformation in solids is the instability mechanism leading to failure by rupture.

In the framework of continuum modeling, this phenomenon is captured by the loss of ellipticity in the governing equations, may lead to discontinuous strain solutions. To better understand the origins of continuum models that exhibit loss of ellipticity at adequate levels of strain or stress, a substantial amount of work has been dedicated to the nonlinear homogenization of microstructured solids to study how geometry and constitutive laws at the microscopic level lead to a macroscopic loss of ellipticity. Since a loss of ellipticity is the property allowing for discontinuous equilibrium solutions, it is thought that a loss of macroscopic (homogenized) ellipticity results in a localized deformation pattern in the post-bifurcated regime. Although this is the case in many engineering applications, it raises interesting questions: is there always a localized deformation appearing in the post-critical equilibrium path of solids loosing macroscopic ellipticity and what are the sufficient conditions in the homogenized response that lead to localization?

The present chapter answers these questions in the framework of a simple, but analytically tractable microstructure, namely an infinite, layered, locally stable (i.e. point-wise rank-one convex) nonlinear (hyperelastic) solid under plane strain loading conditions and more specifically under axial compression along the lamination direction. For this problem, one can find macroscopic loads where the moduli of the homogenized principal solution loose ellipticity (and since the solid has an energy density, the corresponding homogenized energy looses rank-one convexity). Moreover we can ensure that the critical (i.e. corresponding to the lowest applied load) buckling mode of the infinite, perfect solid is global in nature (infinite wavelength eigenmode), a property that allows us to find a homogenized solution for the bifurcated equilibrium path of the infinite, perfect structure. Within this homogenized framework, we prove that perfect, neo-Hookean composites (incompressible or compressible) of arbitrary fiber/matrix volume fraction, always have a stable initial post-bifurcation response with increasing force and displacement (Λ 2 > 0, λ 2 > 0). We also construct a more general composite, consisting of a neo-Hookean fiber and a softer matrix of equal thicknesses which, depending on the constitutive details of the matrix, can exhibit all possible initial post-bifurcation responses: a snap-back (Λ 2 < 0, λ 2 < 0), a maximum force (Λ 2 < 0, λ 2 > 0) or a stable bifurcation occurring under increasing force and displacement (Λ 2 > 0, λ 2 > 0).

The important question addressed next is how a finite-size volume of such a composite will behave in the bulk, i.e. away from its boundaries where strain concentration can easily appear. To this end, we consider a square sample (containing 40 unit cells and using a refined FEM mesh for the corresponding calculations) of the composite with periodic boundary conditions that can capture the homogeneous principal and bifurcated equilibrium paths of the perfect structure. We show, by means of introducing a small geometric imperfection at the middle of the sample, that under these boundary conditions composites with a monotonically increasing force (and displacement) homogenized initial post-bifurcation response (Λ 2 > 0, λ 2 > 0) will evolve towards a uniform shearing solution away from the macroscopic critical load and show no localization of deformation pattern past the critical load, in spite of a macroscopic loss of ellipticity of the principal solution. As expected, composites exhibiting a the snap-back (i.e. maximum displacement) in their macroscopic response (Λ 2 < 0, λ 2 < 0), will evolve into a solution with a single strong localized deformation zone.

The key concept for finding whether loss of macroscopic ellipticity leads to localization of deformation lies in the post-bifurcation behavior of the solid under investigation. Providing consistent criteria based on homogenization ideas, for the absence or presence of localized deformation zones in the post-critical regime of finitely strained solids, is possible, as the current layered composite model shows, for cases when the critical bifurcation mode is also global in nature and well separated from other eigenmodes. For more complex problems, such as solids with two-or three-dimensional periodic microstructures, the presence or absence of localized deformation patterns cannot be answered by using homogenization ideas, because of the local nature of the bifurcated solutions (finite number of unit cells involved). Instead efficient numerical calculations of their equilibrium paths that take into account the symmetry groups of these structures (both point and space) are needed in order to find their solutions well past the onset of a first instability (see [START_REF] Combescure | Post-bifurcation and stability of a finitely strained hexagonal honeycomb subjected to equi-biaxial in-plane loading[END_REF]).

2.A. CRITICAL LOAD OF AN AXIALLY COMPRESSED LAYERED SOLID

2.A Critical load of an axially compressed layered solid

2.A.1 Finite wavelength (local) bifurcation eigenmodes

To find conditions for a non-trivial solution to the system (2.2.1), (2.2.2) and (2.2.3), one can take advantage of its X 1 translational invariance by considering a Fourier transform with respect to X 1 .

This linear system of partial differential equations and interface conditions with piecewise constant coefficients that governs the onset of bifurcation of the layered solid is thus reduced to the following linear system of ordinary differential equations and interface conditions in X 2 :

ω 2 1 L i1k1 vk -iω 1 (L i2k1 + L i1k2 )v k,2 -L i2k2 vk,22 = 0, (2.A.1) [[iω 1 L i2k1 vk + L i2k2 vk,2 ]] = 0, [[v i ]] = 0, (2.A.2)
where v(ω 1 , X 2 ) is the Fourier transform of ∆u(X 1 , X 2 ) and the real number ω 1 is the Fourier transform variable corresponding to X 1 . It is assumed that the field ∆u is uniformly bounded and has adequate continuity, in which case its Fourier transform v exists in the sense of distributions.

To determine a non-trivial solution v (up to a multiplicative constant), for the periodic system of ordinary differential equations in (2.A.1) and (2.A.2), the system is solved on just one unit cell together with some additional boundary conditions at its ends X 2 = 0 + and X 2 = H + . These conditions are provided by Floquet theory, which applies to linear systems of ordinary differential equations in X 2 , with periodic coefficients (period is the unit cell thickness H), according to which:

vi (ω 1 , H + ) = exp (iω 2 H) vi (ω 1 , 0 + ), (2.A.3)
where the real number ω 2 (ω 2 H ∈ [0, 2π)) is the Floquet parameter of the solution.

The general solution to the system of ordinary differential equations with piecewise constant coefficients (2.A.1) is found in each layer to be the sum of four linearly independent partial solutions:

vk (ω 1 , X 2 ) = 4 j=1 m C (j) k exp iω 1 m Z (j) X 2 ; X 2 ∈ (0, H m ), vk (ω 1 , X 2 ) = 4 j=1 f C (j) k exp iω 1 f Z (j) X 2 ; X 2 ∈ (H m , H), vk (ω 1 , X 2 ) = 4 j=1 m C * (j) k exp iω 1 m Z (j) X 2 ; X 2 ∈ (H, H + H m ), (2.A.4)
where Z (j) (j = 1, 4) are the four complex roots of the following fourth order, biquadratic polynomial in Z:

det L i2k2 Z 2 + (L i2k1 + L i1k2 )Z + L i1k1 = 0, (2.A.5)
and C (j) is the eigenvector of the 2 × 2 matrix shown in brackets in (2.A.5) associated with the root Z (j) . The eigenvector components C (j)

1 and C

(j)

2 are related by:5

C (j) 2 = D (j) C (j) 1 ; D (j) := - L 1212 Z 2 (j) + L 1111 (L 1221 + L 1122 )Z (j) , (2.A.6)
Note that equations (2.A.5) and (2.A.6) are valid for each of the two layers and that superscripts m and f are omitted from these equations in the interest of notational simplicity. The requirement that the roots Z (j) are complex for the loading parameter λ of interest, stems from the assumed strong rank-one convexity of each layer, which implies the absence of any discontinuous deformation gradients in each layer for all loading paths considered here.

The Fourier transform of the interface conditions (2.A.2), after substituting equations (2.A.4) and (2.A.6), gives the following equations for the coefficients

m C 1 , f C 1 and m C * 1 in matrix form: m V exp iω 1 m ZH m m C 1 = f V exp iω 1 f ZH m f C 1 , f V exp iω 1 f ZH f C 1 = m V exp iω 1 m ZH m C * 1 , (2.A.7)
for the interfaces X 2 = H m and X 2 = H m + H f = H, respectively. The components of the 4 × 4 matrices V and Z are defined by:

V 1j = 1, V 2j = D (j) , V 3j = L 1212 Z (j) + L 1221 D (j) , V 4j = L 2211 + L 2222 Z (j) D (j) , Z ij = δ ij Z (j) .
(2.A.8)

The components of the vector C 1 are the four constants C

C * 1 = exp (iω 2 H) exp -iω 1 m ZH m C 1 .
(2.A.9) 2.A. CRITICAL LOAD OF AN AXIALLY COMPRESSED LAYERED SOLID Finally, after employing the above result (2.A.9) into equations (2.A.7), a non-trivial solution v(ω 1 , X 2 ) = 0 (or equivalently m C 1 = 0) exists if the matrix with constant coefficients K has unimodular eigenvalues: A.10) where I is the 4 × 4 identity matrix. It should be noted here that the 4 × 4 matrix K(λ, ω 1 H) satisfies: A.11) and that has real invariants I K i ∈ R, (i = 1, . . . , 4). Thus, the critical load parameter λ c , which represents the first occurrence of a bifurcation in the layered solid during a monotonically increasing loading history, corresponds to the first occurrence of a singular matrix in (2.A.10) 1 , as the loading parameter λ increases from zero, over all possible pairs of dimensionless wavenumbers ω 1 H and

det [K(λ, ω 1 H) -exp (iω 2 H) I] = 0; K := f K m K, l K := l V exp iω 1 l ZH l l V -1 ; l = m, f, (2.
det[K] = 1, K -1 (λ, ω 1 H) = K(λ, -ω 1 H), (2.
ω 2 H.
The calculation works as follows: At criticality y c := exp (i(ω 2 H) c ) is an eigenvalue of matrix K, and in view of the fact that |y c | = 1), it corresponds to the first occurrence of a unimodular root of the following fourth order equation:

y 4 -y 3 I K 1 (λ, ω 1 H) + y 2 I K 2 (λ, ω 1 H) -yI K 1 (λ, ω 1 H) + 1 = 0, (2.A.12)
where the two real invariants I K j of matrix K, satisfy in view of (2.A.11) 3 :

I K 1 := tr K, I K 2 := 1 2 (tr K) 2 -tr K 2 , I K j (λ, ω 1 H) = I K j (λ, -ω 1 H).
(2.A.13)

For a fixed dimensionless wavenumber ω 1 H, λ(ω 1 H) is the lowest value of λ for which the characteristic equation (2.A.12) admits a unimodular solution for y and corresponds to the first bifurcation load with a dimensionless wavelength ω 1 H along the X 1 direction. It is given by the lowest positive λ root of one of three equations:

λ(ω 1 H) = min                      λ > 0 2I K 1 -I K 2 -2 = 0, (ω 2 H) c = 0, 2I K 1 + I K 2 + 2 = 0, (ω 2 H) c = π, 1 4 (I K 1 ) 2 -I K 2 + 2 = 0, (ω 2 H) c = cos -1 (I K 1 /4).                      . (2.A.14)
The critical load parameter λ c is then found by a numerical search as the lowest value of λ(ω 1 H), when the infimum is taken over an adequately large interval ω 1 H ∈ R, in the process also giving the corresponding critical dimensionless wavenumber (ω 1 H) c , as defined in (2.2.5).

Notice that in view of (2.A.13) 3 only positive wavenumbers (ω 1 H > 0) need to be considered, thus explaining the definition of the critical load in (2.2.5). Some additional comments are in line at this point. Depending on whether the infimum in (2.A.14) occurs at the first, second or third equation, one can determine the nature of the critical eigenmode. More specifically, when the infimum occurs for (ω 2 H) c = 0, the corresponding eigenmode is H-periodic, as shown in Fig. 2.12d, when the infimum occurs for (ω 2 H) c = π, the corresponding eigenmode is H-antiperiodic, as shown in Fig. 2.12e, while for an infimum reached at (ω 2 H) c = cos -1 (I K 1 /4), the corresponding eigenmode is shown in Fig. 2.12f.

2.A.2 Long wavelength (global) bifurcation eigenmodes

It should be noted here, that in all the previous calculations it was tacitly assumed that the dimensionless wavenumber ω 1 H = 0. The function λ(ω 1 H) has a singular point at the origin, i.e. λ(0) = λ(0 + ) := lim ω 1 H→0 + λ(ω 1 H), since two physically different types of modes can exist in the neighborhood of ω 1 H = 0: an X 2 -periodic, X 1 -independent eigenmode v(0, X 2 ) for ω 1 H = 0 and a long wavelength eigenmode with dimensionless wavelength along the

X 1 direction 2π/(ω 1 H) = L 1 /H → ∞ when ω 1 H → 0 + .
For the latter case it has been shown by [START_REF] Triantafyllidis | On the comparison between microscopic and macroscopic instability mechanisms in a class of fiber-reinforced composites[END_REF], that the limit value of λ(ω 1 H), as ω 1 H → 0 + , is λ H which corresponds to the loss of rank-one convexity of the homogenized incremental moduli L H . Although the proof of 2.A. CRITICAL LOAD OF AN AXIALLY COMPRESSED LAYERED SOLID this assertion for the general case of periodic composites is detailed in [START_REF] Geymonat | Homogenization of nonlinearly elastic materials, microscopic bifurcation and macroscopic loss of rank-one convexity[END_REF], for reasons of completeness a brief outline for layered composites in plane strain is presented here.

Starting point for the proof of the above assertion is (2.A.10), where the asymptotic expansion of K(λ, ω 1 H) with respect to ω 1 H gives to first order:

det [G(λ) -(ω 2 /ω 1 ) I] = 0; G 11 (λ) = G 22 (λ) = G 33 (λ) = G 44 (λ) = 0, G 12 (λ) = -L 1221 (L 1212 ) -1 , G 21 (λ) = -L 1122 (L 2222 ) -1 , G 13 (λ) = (L 1212 ) -1 , G 31 (λ) = -L 1111 -(L 1122 ) 2 (L 2222 ) -1 , G 24 (λ) = -(L 2222 ) -1 , G 42 (λ) = L 2121 -(L 1221 ) 2 (L 1212 ) -1 , G 34 (λ) = L 1122 (L 2222 ) -1 , G 43 (λ) = L 1221 (L 1212 ) -1 , G 14 (λ) = G 41 (λ) = G 23 (λ) = G 32 (λ) = 0.
(2.A.15)

Hence (ω 2 /ω 1 ) is an eigenvalue of the matrix G, and must satisfy the following biquadratic equation:

(ω 2 /ω 1 ) 4 + (ω 2 /ω 1 ) 2 I G 2 (λ) + I G 4 (λ) = 0, (2.A.16)
where the invariants of G are real and given in terms of its components by:

I G 2 (λ) = -(G 12 G 21 + G 13 G 31 + G 24 G 42 + G 34 G 43 ), I G 4 (λ) = (G 12 G 43 -G 13 G 42 )(G 21 G 34 -G 31 G 24 ).
(2.A.17)

The critical load parameter corresponding to the long wavelength λ(0

+ ) := lim ω 1 H→0 + λ(ω 1 H)
is then found as the lowest value of λ for which the biquadratic (2.A.16) admits real solutions in (ω 2 /ω 1 ), namely the lowest λ root of one of two equations:

λ(0 + ) = min          λ > 0 I G 4 = 0, (ω 2 /ω 1 ) c = 0, 1 4 (I G 2 ) 2 -I G 4 = 0, (ω 2 /ω 1 ) c = -(I G 2 /2) 1/2 .          . (2.A.18)
The proof that λ(0 + ) coincides with λ H , the lowest load parameter corresponding to the loss of the rank-one convexity of the homogenized moduli L H (λ) follows from (2.A.15) and (2.A.16) and the expressions for the homogenized moduli, which are presented next.

For the orthotropic layered medium at hand subjected to axial loading the determination of L H (λ) starts with the following definition:

L H ijkl (λ) = 1 H H 0 L mnpq δ im δ jn + ij ϕ m,
L H 1111 (λ) = L 1111 -(L 1122 ) 2 (L 2222 ) -1 + L 1122 (L 2222 ) -1 2 (L 2222 ) -1 -1 , L H 1122 (λ) = L 1122 (L 2222 ) -1 (L 2222 ) -1 = L H 2211 (λ), L H 2222 (λ) = (L 2222 ) -1 -1 , L H 2121 (λ) = L 2121 -(L 1221 ) 2 (L 1212 ) -1 + L 1221 (L 1212 ) -1 2 (L 1212 ) -1 -1 , L H 1221 (λ) = L 1221 (L 1212 ) -1 (L 1212 ) -1 = L H 2112 (λ), L H 1212 (λ) = (L 1212 ) -1 -1 .
(2.A.21)

The first, as the load parameter increases, loss of ellipticity for the homogenized, layered solid, corresponds to the lowest load parameter λ H for which the homogenized incremental moduli L H (λ) loose rank-one convexity, i.e. λ H is the lowest λ-root of (2.2.4). From the orthotropy of the homogenized incremental moduli L H (λ), the determinant of the homogenized acoustic tensor L H ijkl (λ)n j n l in (2.2.6) is the following biquadratic equation:

(n 2 /n 1 ) 4 + (n 2 /n 1 ) 2 I H 2 (λ) + I H 4 (λ) = 0, (2.A.22)
where the coefficients I H 2 and I H 4 are expressed in terms of the components of the homogenized moduli tensor by:

I H 2 (λ) = (L H 1111 L H 2222 + L H 1212 L H 2121 ) -(L H 1122 ) + L H 2222 ) 2 (L H 1212 L H 2222 ) -1 , I H 4 (λ) = (L H 1111 L H 2121 )(L H 1212 L H 2222 ) -1 .
(2.A.23) However, from the expressions for the invariants of G in (2.A.15) and the expressions for the components of the homogenized moduli tensor in (2.A.21) one can show that:

I G 2 (λ) = I H 2 (λ), I G 4 (λ) = I H 4 (λ), (2.A.24)
thus proving our assertion that λ(0 + ) = λ H and moreover, since the biquadratics in (2.A.16) and

(2.A.22) coincide, that (ω 2 /ω 1 ) c = (n 2 /n 1 ) c .

2.A.3 Periodic (X 1 -independent) bifurcation eigenmodes

The last case remaining to be checked is ω 1 H = 0, i.e. when the eigenmode is independent on X 1 . The corresponding critical load parameter λ(0) is also found from the transformed governing equations (2.A.1) and (2.A.2). In this case it can be seen from ( 2.A.1), that v1 (0, X 2 ) and v2 (0, X 2 ) are piecewise linear functions in X 2 within each layer. A non-trivial solution vi (0, X 2 ) exists when L 1212 (λ) = 0 or L 2222 (λ) = 0 in either the fiber or the matrix layer, thus giving the critical stretch ratio as the lowest λ root of one of the four equations:

λ(0) = min λ > 0 L f 1212 (λ) = 0, L f 2222 (λ) = 0, L m 1212 (λ) = 0, L m 2222 (λ) = 0 . (2.A.25)
For the hyperelastic solids investigated here, their rank-one convexity guarantees that L 1212 (λ) > 0 and L 2222 (λ) > 0 for both fiber and matrix. Even for rate-independent solids that might loose ellipticity, one can see that λ(0) from (2.A.25) is strictly larger than λ(0 + ) from (2.A.18), itself being by definition larger than λ c , the lowest λ-root of (2.2.5), namely:

λ(0) > λ(0 + ) = λ H ≥ λ c (2.A.26)
Note that the existence of a singularity at λ at ωH 1 = 0 explains the use of the infimum in the definition of (2.2.5).

2.B. POST-BIFURCATION ASYMPTOTICS FOR A COMPRESSIBLE NEO-HOOKEAN COMPOSITE

one obtains the following expressions for p f and p m :

p f = c m (1 + r f ) + γ 2 [(1 + r f ) /ζ m + (1 + r m ) /ζ f ] /µ H c f c m + γ 2 [c m /ζ f + c f /ζ m ] /µ H , p m = c f (1 + r m ) + γ 2 [(1 + r m ) /ζ m + (1 + r m ) /ζ f ] /µ H c f c m + γ 2 [c m /ζ f + c f /ζ m ] /µ H , c f := (1 -λ) 2 + r f /µ f , r f := µ f /k f , c m := (1 -λ) 2 + r m /µ m , r m := µ m /k m .
(2.B.4)

The δ-γ relationship along the bifurcated equilibrium path is obtained from the requirement Π 12 = 0 in (2.2.9) 2 , and the help of (2.B.1), (2.B.3) giving:

Π 12 = γ µ G - (p f -p m ) 2 µ H = 0. (2.B.5)
Consequently, with the help of (2.B.4) the sought δ-γ relationship is:

c m (1 + r f ) -c f (1 + r m ) = (µ G µ H ) 1/2 c f c m + γ 2 µ H c m ζ f + c f ζ m , (2.B.6)
where without loss of generality it was tacitly assumed that p f > p m or equivalently from (2.B.4)

that c m (1 + r f ) > c f (1 + r m )
. Notice that at the incompressible limit r f = r m = 0, in which case (2.B.6) reduces to (2.3.11) 1 as expected.

From (2.B.6) at bifurcation, i.e. for γ = 0, recalling also the definitions of c f , c m , r f , r m in (2.B.4) 3 which give an additional relation between c f and c m :

c m µ m -c f µ f = r m -r f := ∆r, (2.B.7)
one obtains the following result for the critical load λ c under displacement control (λ = -δ):

λ c = 1 - 1 2 ∆M (µ G µ H ) 1/2 + D -(r f + r m ) 1/2 , ∆M := µ f (1 + r f ) -µ m (1 + r m ) , D := (∆M ) 2 µ G µ H + (∆r) 2 + 2 ∆r (µ G µ H ) 1/2 [µ f (1 + r f ) + µ m (1 + r m )] 1/2 . (2.B.8)
As expected, at the incompressible limit r f = r m = 0, the above expression reduces to (2.3.13).

The expression for the curvature λ 2 (λ 2 = -δ 2 ) of the bifurcated equilibrium path at the critical

2.C. INFLUENCE OF CONSTITUTIVE LAW CHOICE ON CRITICAL LOAD

2.C Influence of constitutive law choice on critical load

In all calculations presented in this chapter, the response of each layer of the composite is described by a hyperelastic constitutive law, which cannot take into account the presence of a possible unloading in the bifurcated solution. Since unloading will occur in the stiffer layer (matrix), the neo-Hookean model is adequate for this case, due to its linear τ -γ (shear stress-shear strain) response.

Hence, to make connections with elastoplasticity, one can compare the predictions of the hyperelastic composites to models where the softer (matrix) layer follows a deformation theory constitutive law, for this layer will continue loading in the homogenized post-bifurcated solution. However, only the onset of bifurcation can be compared for the two different models; post-bifurcation calculations for deformation theory of plasticity would require integrating the corresponding rate-independent (hypoelastic) model, unlike the hyperelastic case where no such integrations are needed.

For comparison purposes we revisit below the composite investigated in Subsection 2.3.4, consisting of two equal thickness layers, the stiffer (fiber) being neo-Hookean with energy density

W f = (1/2)(I 1 -2)
and the softer (matrix) being in turn a) hyperelastic, b) deformation theory elastoplastic, both sharing the same uniaxial stress-strain law derived from the matrix energy density. For incompressible plane strain conditions, the incremental form of a rate-independent, pressure-insensitive, initially orthotropic material, as first noted by [START_REF] Biot | Mechanics of incremental deformation[END_REF], takes the general form: respectively. It can be shown, (e.g. see Abeyaratne and Triantafyllidis (1981a)), that for the deformation theory model proposed by Stören and Rice (1975a), which has the same uniaxial stress-strain law as the hyperelastic model (energy density W = (α/2)g(I 1 -2)), only the incremental modulus µ is different between the hyperelastic and deformation theory under loading condition:

µ = α 2 (1 -λ) 2 + 1 (1 -λ) 2 g (z), hyperelastic model, µ = α 4 ln(1 -λ) (1 -λ) 2 - 1 (1 -λ) 2 g (z), deformation theory model, (2.C.2) where z := I 1 -2 = (1 -λ) 2 + (1 -λ) -2 ) -2.
Using the matrix energy density from (2.3.16) in conjunction with (2.C.2), one obtains the following expressions for the critical load under displacement control λ c :

nα (1 -λ c ) 2 + 1 (1 -λ c ) 2 -2 (n-1) = 1 -(1 -λ c ) 2 1 + (1 -λ c ) 2 , hyperelastic model, nα (1 -λ c ) 2 + 1 (1 -λ c ) 2 -2 (n-1) = -ln(1 -λ c ) 2 1 + (1 -λ c ) 4 [1 + (1 -λ c ) 2 ] 2 , deformation theory model. (2.C.3)
Comparing the critical strain λ c of the composite using the two different matrix models is presented in Fig. 2.13. As expected from (2.C.2), since for a given axial strain λ the shear moduli µ are larger for the hyperelastic matrix (while remaining unchanged for the fiber), the stiffness contrast between the two layers is greater for the deformation theory version of the matrix and hence result in a lower critical strain for the corresponding composite. An alternative way of comparing the critical strain λ c of the composite as a function of the initial matrix-to-fiber stiffness ratio α for the two different matrix models for a given hardening exponent (n = 0.55, 0.57, 1) is presented in Fig. 2.14. As explained above, for a given set of material parameters (α, n), the hyperelastic model has a critical load which is consistently higher than its hypoelastic counterpart. However, for small values of α, i.e. large differences between the initial stiffness of the two layers, there is practically no difference between the predictions of the two different matrix models, which means almost identical critical strains in the range 0 ≤ λ c < 0.2, as seen in Fig. 2.14.

Chapter 3

Dynamic stability of biaxially strained thin sheets under high strain rates

Introduction

The issue of dynamic stability of structures is an important engineering problem and as such has drawn considerable attention. The first investigation in this area appears to be the work of [START_REF] Koning | Impact buckling of thin bars in the elastic range hinged at both ends[END_REF], who investigated the influence of inertia in a simply supported imperfect column subjected to a sudden axial load. A substantial amount of work followed that investigated the response of, mainly elastic, structures to impulse or time-dependent loads. As a result, and due to the many possible definitions for the stability of time-dependent systems, the term dynamic stability encompasses many classes of problems and different physical phenomena and has many interpretations, with inertia being the only common denominator.

In the absence of inertia, the processes of failure by a bifurcation instability mode in elastic solids and structures is well understood (e.g. [START_REF] Brush | Buckling of bars, plates, and shells[END_REF]) and a general asymptotic analysis, termed Lyapunov-Schmidt-Koiter (LSK), has been developed for their study. The first effort to use the LSK general analysis for the dynamic stability problem of an elastic structure appears to be [START_REF] Budiansky | Dynamic buckling of imperfection sensitive structures[END_REF], where the authors proposed an asymptotic analysis of the time-dependent problem using the eigenmodes of the static problem.

Another idea, popular in fluid mechanics, has also been adopted for the dynamic stability analysis of solids with more general constitutive laws under high rates of loading, according to which one seeks the solid's fastest growing eigenmode, or the wavelength associated with lowest necking strain. This method has been repeatedly applied in the study of dynamic stability of various elastoplastic structures (bars, rings, plates, shells etc.) under high loading rates where the failure pattern and size of fragments is of interest (e.g. see [START_REF] Guduru | The dynamics of multiple neck formation and fragmentation in high rate extension of ductile materials[END_REF]; [START_REF] Jouve | Analytic study of the onset of plastic necking instabilities during biaxial tension tests on metallic plates[END_REF]; [START_REF] Mercier | Multiple necking duringt hedynamic expansion of hemispherical metallic shells, from experiments to modelling[END_REF]; [START_REF] Mercier | Predictions of bifurcation and instabilities during dynamic extension[END_REF]; [START_REF] Shenoy | Necking bifurcations during high strain rate extension[END_REF]; [START_REF] Sorensen | Unstable neck formation in a ductile ring subjected to impulsive radial loading[END_REF]; [START_REF] Xue | Material aspects of dynamic neck retardation[END_REF]; [START_REF] Zhou | An elastic-visco-plastic analysis of ductile expanding ring[END_REF]).

However, recent experimental evidence from rapidly expanding electromagnetically loaded metallic rings by Zhang andRavi-Chandar (2006, 2008) finds no evidence of a dominant wavelength at the necked pattern of the rings. Moreover, they find no experimental evidence of influence of strain rate on the necking strains, which are consistent with maximum force criterion of a rate-independent constitutive law (Considère criterion).

As explained by these authors, using the fastest growing eigenmode to predict the onset of failure is physically meaningful provided that some characteristic velocity of the principal solution-e.g.

ring/tube expansion rate-is much slower than the speed of propagation of perturbations in the solid or structure at hand. For high loading rates, commensurate with some characteristic wave propagation speed in the structure, a novel approach to the stability analysis is required, namely the study of evolution of localized perturbations.

We are motivated by the experimental studies of Zhang andRavi-Chandar (2006, 2010) on the high strain rate expansion of thin rings and tubes, that show no evidence of a dominant wavelength in their failure mode and no influence of strain rate sensitivity on the necking strains-the onset of failure of an electromagnetically loaded, dynamically expanding tube is shown in Figure 3.1, where one can observe a rather random failure mode. Recently, [START_REF] Ravi-Chandar | Dynamic stability of a bar under high loading rate: Response to local perturbations[END_REF] studied the dynamic stability of an incompressible, nonlinearly elastic bar at different strain rates by following the evolution of localized small perturbations introduced at different times. The same approach is followed here for the biaxial stretching of thin plates, where we follow the time evolution of spatially localized perturbations and their interactions. Following this introduction, in Section 3.1, the formulation of the problem (definition of influence zones and algorithm for the nonlinear FEM calculations) is presented in Section 3.2. Results are given in Section 3.3, starting with the description of the constitutive models and following with the numerical study of the nonlinear time evolution of a such a perturbation showing that these structures are stable until the time when the condition for the loss of ellipticity is reached. In the same section we present an analytical method, based on linearization, to define the size of the influence zone of a point-wise perturbation and we study its dependence on constitutive laws and loading conditions. The presentation is concluded in Section 3.4, while details on the constitutive models used are presented in the appendix.

Problem Formulation

This section starts with the presentation of the model for the propagation of perturbation about a point defect on a biaxially strained, rate-independent, flat plate of infinite extent by studying the evolution of influence zones (linearized approach). The setting of the corresponding nonlinear problem, which is treated numerically, is presented subsequently.

The influence zones of a biaxially strained elastoplastic plate

We consider a two-dimensional thin, flat plate (idealized as a membrane) of infinite extent and uniform initial thickness H subjected to a biaxial stretching as shown in Figure 3.2. To avoid inplane acceleration terms in the unperturbed solution of the perfect plate, the following stretch ratios are being imposed at infinity:

λ 1 = 1 + (c cos ψ)t, λ 2 = 1 + (c sin ψ)t. (3.2.1)
Hence a uniaxial straining corresponds to tan ψ = 0, balanced biaxial straining to tan ψ = 1, while a uniaxial stressing is approximated by tan ψ = -1/2 (assuming incompressibility and valid only initially for small strains. A more accurate value for finite strains depends on constitutive response).

In the absence of body forces the equations of motion of the thin plate can be put in the form1 :

∂N αβ ∂X β = ρ 0 ∂ 2 u α ∂t 2 , (3.2.2)
where N αβ is the nominal (force/reference thickness) stress resultant, ρ 0 the reference mass density, X α the reference geometric coordinates and u α (X, t) the corresponding displacement of a material point initially at X in the Lagrangian description. For simplicity the reference configuration is identified with the stress-free configuration of the plate.

The plate's constitutive equation is assumed to be rate-independent with the following relation between the time derivative of the stress measure Ṅαβ and rate of deformation gradient Ḟαβ2 : Of interest here is the propagation of the localized perturbation about X = 0. Using ∆f to denote the difference between the perturbed and unperturbed values respectively of a field quantity f and exploiting the fact that the principal solution is homogeneous (i.e. independent of X), one obtains the following systems governing the evolution of perturbation3 :

Ṅαβ = L αβγδ Ḟγδ ; F ζ = δ ζ + ∂u ∂X ζ , ( 3
L αβγδ ∂ 2 ∆u γ ∂X δ ∂X β = ρ 0 ∂ 2 ∆u α ∂t 2 . (3.2.4)
We follow the propagation of the perturbation in all directions n and for this purpose we consider solutions of the form:

∆u α (X, t) = U α f (V t -n α X α ) , (3.2.5)
where f (z) is an arbitrary function of z ∈ R, U α is the perturbation amplitude and V its speed of propagation.

Introducing (3.2.5) into (3.2.4) one finds that (V ) 2 is one of the two eigenvalues of L αβγδ n β n δ := A αγ (n) which is the acoustic tensor along n:

L αβγδ n β n δ -ρ 0 (V ) 2 δ αγ U γ = 0. (3.2.6)
The absence of any characteristic length ratio in the infinite flat plate, leads to the following nondimensionalization of perturbation velocity (ν), distance (χ), and time (τ ) for this problem

ν := V [G/ρ 0 ] -1/2 ; χ := Xc [G/ρ 0 ] -1/2 ; τ := ct, (3.2.7)
where [G/ρ 0 ] -1/2 is the shear wave propagation speed at zero strain (G being the corresponding shear modulus of the plate at zero strain) and c the straining speed introduced in (3.2.1).

We are now in a position to estimate the extent of the zone influenced by a point-wise perturbation introduced at the onset of deformation (τ = 0) in the plate until the onset of a necking localization, characterized by the loss of ellipticity in the perfect thin plate. We also call this perturbation an initial defect, since the introduction of the spatially localized perturbation at the beginning of the loading process is equivalent to the presence of a defect in the plate. To this end we define the lowest and highest dimensionless wave propagation speeds ν -and ν + which are the lowest and highest eigenvalues of the acoustic tensor according to (3.2.6) appropriately non-dimensionalized with the help of (3.2.7).

For a given time τ , one can thus define the influence zones χ -and χ + determined by: (3.2.8) which are the distances travelled at time τ by the slow and fast wave fronts respectively along direction φ (recall n 1 = cos φ, n 2 = sin φ).

χ -(φ, τ ) := τ 0 ν -φ, τ dτ , χ + (φ, τ ) := τ 0 ν + φ, τ dτ ,
The perturbation can no longer propagate along a direction φ once its lowest speed reaches ν -= 0, which occurs at time τ e (φ) (subscript e standing for loss of ellipticity in the incremental equilibrium equations):

ν -(φ, τ e (φ)) = 0 (ν -(φ, τ ) > 0 for 0 ≤ τ < τ e (φ)) . (3.2.9)
The locus of points reached by the wave propagating along a given direction φ until ν -= 0 is given by:

χ e (φ) := τe(φ) 0 ν -φ, τ dτ , ( 3
δ -= χ -(φ m , τ m ) = τm 0 ν -φ m , τ dτ . (3.2.12)
In a similar way we are interested in the maximum size disc, centered at X = 0, that the perturbation can reach. At time τ M the plate has already lost ellipticity for all possible directions of wave propagation φ and the zone influenced by the perturbation at X = 0 is χ + (φ, τ M ) as seen in Figure 3.4. In analogy to δ -, one can also define δ + the radius of the maximum disc covering entirely the range of influence of the perturbation at X = 0 when the plate has lost ellipticity along all possible directions φ, namely:

δ + = χ + (φ M , τ M ) = τ M 0 ν + φ M , τ dτ . (3.2.13)
It should be noted here that δ -exists as long as the model looses ellipticity for the loading considered. Moreover even if a δ -exists, a δ + might not. 51

PROBLEM FORMULATION

δ - δ + χ -(φ, τ m ) χ + (φ, τ M ) χ e (φ)

Proof of linearized stability

In order to get an analytically tractable way to investigate the stability of the plate's uniform stretching solution (3.2.1), we start by studying the linearized system of the evolution of a perturbation in (3.2.2). By defining the displacement and particle velocity perturbations as: (3.2.14) and introducing them to the time derivative of (3.2.2), one obtains upon linearization about the uniform stretching solution the following equation:

δu α (X, t) := u α (X, t) -u 0 α (t), δv α (X, t) := v α (X, t) -v 0 α (t),
∂L αβ ζ ∂F γδ 0 Ḟ 0 ζ δu γ,δβ + L 0 αβγδ δv γ,δβ = ρ δv α (3.2.15)
using the constitutive relation of the principal solution

L 0 αβγδ δu γ,δβ = ρ δ vα . (3.2.16)
We can write down the the Fourier transform of (3.2.15) and (3.2.16) using the mapping of notations of δu → ∆u and X -→ ω:

-ω δ ω β ∂L αβ ζ ∂F γδ 0 Ḟ 0 ζ ∆u γ -ω δ ω β L 0 αβγδ ∆v γ = ρ ∆v α (3.2.17) -ω δ ω β L 0 αβγδ ∆u γ = ρ ∆ vα (3.2.18)
By replacing ∆u terms in (3.2.17) using (3.2.18), we obtains an equation expressed in particle velocity ∆v

ρ ∂L αβ ζ ∂F γδ 0 Ḟ 0 ζ ω δ ω β [L γηξθ ω η ω θ ] -1 ∆ vξ -ω β ω δ L 0 αβγδ ∆v γ = ρ ∆v α (3.2.19)
We now define the wavenumber ω := ω 2 1 + ω 2 2 and the polarization unit vector

(n 1 , n 2 ) satisfying ω 1 = ω n 1 , ω 2 = ω n 2 , and n 2 1 + n 2 2 = 1. (3.2.19
) can thus be rewritten as:

-M (n) αγ ∆ vγ + ω 2 A (n) αγ ∆v γ + δ αγ ∆v γ = 0 (3.2.20)
with acoustic tensor A (n) and damping tensor M (n) defined by

A (n) αγ := (1/ρ) L αβγδ n β n δ (3.2.21) M (n) αξ := (1/ρ) ∂L αβ ζ ∂F γδ 0 Ḟ 0 ζ n β n δ A -1 γξ (3.2.22) 53 3.2. PROBLEM FORMULATION
We assume that the tangential moduli L αβγδ is derived from an internal energy, namely L αβγδ = ∂ 2 W/∂F αβ ∂F γδ . With this assumption, (3.2.22) can be considerably simplified in form:

M (n) αξ = (1/ρ) ∂ 3 W ∂F αβ ∂F γδ ∂F ξ 0 Ḟ 0 ζ n β n δ A -1 γξ = D Dt (1/ρ) L αβγδ n β n δ A -1 γξ = Ȧαγ A -1 γξ (3.2.23)
Using the method of frozen coefficients, i.e. assuming that the rate of growth/decay of the perturbation is much higher than the loading rate, (3.2.20) is considered as a constant coefficient equation, which admits solutions of form ∆v = I exp (s I t) V I where s I are the 4 roots of eigenvalue equation:

det -sM αγ + ω 2 A αγ + δ αγ s 2 = 0 (3.2.24)
Our next step is to prove that the real parts of all s I are negative. This relies on that fact that A is derived from real-world material behavior, and has special properties. Before the material loses ellipticity, A (n) is for all polarization unit vectors n a positive-definite matrix admitting two positive eigenvalues, noted as α

+ (n) = ν 2 + (n) and α -(n) = ν 2 -(n),
   - α+ α + s + ω 2 α + + s 2 -β α + -α - α - s -β α + -α - α + s - α- α - s + ω 2 α -+ s 2      V + V -   = 0 (3.2.28)
which is a fourth-order polynomial equation

s 4 + a 3 s 3 + a 2 s 2 + a 1 s + a 0 = 0 (3.2.29)
with coefficients (3.2.30)

a 3 = - α+ α + - α- α - , a 2 = det Ȧ α + α - + (α + + α -) ω 2 ,
a 1 = -α + α- α - + α - α+ α + ω 2 , a 0 = α + α -ω 4 . (3.2.31)
The weakening behavior of the metallic material in plastic regime grants two more properties:

1) det Ȧ > 0 and 2) α+ < 0, α-< 0, with which we can define five positive real values

A = -( α+ /α + + α-/α -) , B = det Ȧ, Σ = (α + + α -) , Γ = -( α-α + /α -+ α+ α -/α + ) , ∆ = α + α -. (3.2.32)
and express (3.2.31) with them:

a 3 = A, a 2 = B + Σ ω 2 , a 1 = Γω 2 , a 0 = ∆ ω 4 . (3.2.33)
We recall the special case of the Routh-Hurwitz theorem: the neccessary and sufficient conditions for the fourth order equation a 4 s 4 + a 3 s 3 + a 2 s 2 + a 1 s + a 0 = 0 to have all four roots with negative real parts are:

         a 1 > 0, a 2 > 0, a 3 > 0, a 4 > 0 a 3 a 2 -a 4 a 1 > 0 a 3 a 2 a 1 -a 4 a 2 1 -a 2 3 a 0 > 0 . (3.2.34)
It is evident that (3.2.34) 1 is satisfied. For (3.2.34) 2 we have

a 3 a 2 -a 4 a 1 = AB -(Γ -AΣ)ω 2 = AB -( α+ + α-) ω 2 > 0 (3.2.35)
and for (3.2.34) 3 we have (after some algebra):

a 3 a 2 a 1 -a 4 a 2 1 -a 2 3 a 0 = ABΓω 2 + AΓΣ -A 2 ∆ -Γ 2 ω 4 = ABΓω 2 + α+ α-(α + -α -) 2 / (α + α -) ω 4 > 0 (3.2.36)
The conditions are hence satisfied, and we can prove that all four s I have negative real parts, this is to say for all wavenumbers ω, the corresponding solution ∆v = I exp (s I t) V I decays in time. Now we want to characterize the dependence of s I on ω when ω tends to infinity and 0 by partially solving (3.2.29). Using Ferrai-Cardano formula, one can express the roots of the fourth order polynomial s 4 + a 3 s 3 + a 2 s 2 + a 1 s + a 0 = 0 as:

s = - a 3 2 (1) ± γ 2 (2) ± 1 2 a 2 3 2 - 4a 2 3 -ζ (1) ± η 4γ 1/2 (3.2.37) 55 3.2. PROBLEM FORMULATION with some new notations γ = a 2 3 /4 -2a 2 /3 + ζ 1/2 , η = -a 3 3 + 4a 3 a 2 -8a 1 , ζ = g(k/2) 1/3 + (2/k) 1/3 / 3, g = a 2 2 -3a 3 a 1 + 12a 0 , k = h + h 2 -4g 3 1/2 , h = 2a 3 2 -9a 2 (a 3 a 1 + 8a 0 ) + 27 a 2 1 + a 2 3 a 0 . (3.2.38)
By investigating high order terms of these parameters, we have for ω → ∞:

η ∼ (4AΣ -8Σ) ω 2 , g ∼ 12∆ + Σ 2 ω 4 , h ∼ 2Σ Σ 2 -36∆ ω 6 , γ ∼ ζ1/2 , ζ ∼ ζ, k ∼ 2 Σ + 2 √ 3i∆ 1/2 3 ω 6 (3.2.39)
where ζ := (1/3) 12∆ + Σ 2 Σ + 2 √ 3i∆ 1/2 ω 6 . We conclude that s have the same order as ζ1/2 = O ω 3 when ω tends to infinity. By investigating the constant terms of these parameters, we have for ω → 0:

η ∼ A 4B -A 2 , g ∼ B 2 , h ∼ 2B 3 , γ ∼ γ, ζ ∼ B 3 + 1/B / 3, k ∼ 2B 3 (3.2.40) with γ := A 2 /4 -2B/3 + B 3 /3 + 1/(3B) 1/2 and s I ∼ - A 2 (1)
± γ 2

(2)

± 1 2 A 2 2 - 4B 3 - B 3 3 - 1 3B (1)
± A 4B -A 2 4γ 1/2 (3.2.41)
We also need the initial conditions to complete our proof of linearized stability. Note V + I , V -I as β(α +α -)/α -, -s I α+ /α + + ω 2 α + + s 2 I , the solutons of (3.2.20) can be written as

∆v = I λ I exp (s I t) V + I e + + V - I e - (3.2.42)
where λ I is real valued function of ω and should be determined by initial conditions. Without loss of generality, let's consider a localized initial condition with the help of δ D (x) the Dirac function located at the origin

δv (X; 0) = δ D (x) e + , δ v (X; 0) = 0. (3.2.43)
With these initial conditions, λ I are solutons of the following system: s I have negative real parts and are of order ω 3 when ω tends to infinity and tends to constant when ω tends to zero, λ I and V I are all of polynomial order in ω on the open boundaries of R + , hence there is always a large Ω ∈ R + that the before-mentioned integral is well defined on R 2 /{ω < Ω}. It is also defined on the compact {X ω < Ω}, so the inverse of Fourier transform exists and decays to zero when the time t tends to infinite. Thus is proven the linearized stability before the material loses its ellipticity. Fig. 3.5 shows a verification of the linearized stability of perturbation before the material loses ellipticity using a set of interconnected numerical simulations. In these calculations, instead of introducing an imperfection in the material at the beginning of the simulation, a small axisymmetric perturbation in displacement is directly introduced into the perfect solution at a chosen time of the simulation. The axisymmetric perturbation in displacement takes the following form, which always gives finite deformation measure over its domain of definition:

I λ I V + I = 1/ √ 2π, I λ I V - I = I λ I s I V + I = I λ I s I V - I = 0 (3.2.
∆u r = A(r/R) exp (-r/R) , ∆u θ = 0; A < Re 2 (3.2.46)
The condition on A (used as amplitude parameter here, definition different from (3.2.32)) and R (size parameter) should be ensured in order to prevent inversion of the space orientation. In the Fig. 3.5, the perturbations introduced before τ m first rapidly decay, and then eventually localize after time passes τ m (lines in black to orange color). On the other hand, the perturbations introduced after τ m localize directly without decaying (lines in light orange or yellow color).

Numerical Calculations

Nonlinear dynamics calculations for the evolution of a spatially localized initial perturbation are done using the finite element method using the simplest constitutive law, namely a hyperelastic, finite (logarithmic) strain model fitted a uniaxial power law. The use of this model for elastoplastic materials is justified by the absence of unloading in the calculations (of interest is the response of the plate up to the loss of ellipticity). A brief description of the algorithm used is presented here.

The starting point of the calculations is the weak formulation of equilibrium equations (3.2.2)

in Lagrangian (reference) configuration:

A [N αβ δu α,β + ρ 0 ∂ 2 u α ∂t 2 δu α ]dA = 0, (3.2.47)
with δu as the test function. A spatial FEM discretization leads to the solution of the following system of equations:

M • A(t) + F(U(t)) = 0; A(t) := ∂V(t)/∂t, V(t) := ∂U(t)/∂t (3.2.48)
where F is the force vector, M the mass matrix and U(t), V(t) and A(t) respectively vector of nodal displacements, velocities and accelerations.

The time marching algorithm chosen for the solution of (3.2.48) is the HHT-α method (see [START_REF] Hilber | Improved numerical dissipation for time integration algorithms in structural dynamics[END_REF] which uses the following updating scheme for the displacement and velocity vectors:

U t+∆t = U t + ∆tV t + (∆t) 2 2 (1 -2β)A t + 2βA t+∆t , V t+∆t = V t + ∆t (1 -γ)A t + γA t+∆t , (3.2.49)
which are in turn used for the iterative solution of (3.2.48) by driving its residual vector R at each time step to zero according to:

0 = R t+∆t = M • A t+∆t + (1 + α)F(U t+∆t ) -αF(U t ).
(3.2.50)

In the above expressions, the constants α, β and γ govern the stability and numerical dissipation of the algorithm and are related by β = (1α) 2 /4 and γ = 1/2α (see [START_REF] Hilber | Improved numerical dissipation for time integration algorithms in structural dynamics[END_REF]). For the calculations reported here we choose α = -0.05. The remainder R t+∆t i at iteration i at time step t + ∆t is updated using the tangent stiffness matrix K t+∆t i of the algorithm:

R t+∆t i = K t+∆t i • (U t+∆t i+1 -U t+∆t i ), K t+∆t i := ∂R ∂U t+∆t i = 1 β∆t 2 M + (1 + α) ∂F ∂U (U t+∆t i ), (3.2.51) until convergence in the displacement is reached, i.e. U t+∆t i+1 -U t+∆t i ≤ U t+∆t i
, where is a conveniently chosen tolerance parameter.

The spatial discretization of the plate uses standard two-dimensional isoparametric quadrilateral elements. The boundary conditions imposed at any time are the displacements and velocities of the principal (perfect) solution given in (3.2.1). The reference configuration rectangular domain used in the calculations covers completely a disc of radius δ + , thus ensuring that no perturbation wave ever reaches any boundary. The initial conditions are the displacements and velocities of the principal solution. Instead of prescribing a slightly perturbed initial displacement or velocity field, we chose the equivalent approach of using an initial imperfection through a central element with a slightly lower shear modulus, as detailed in the results section. A special element incorporating the constitutive law and time solution algorithm described here is then introduced in a commercial FEM code (ABAQUS) to calculate the results for this chapter.

Results

This section starts with the constitutive models chosen. It continues with the study of the evolution of a single spatially localized perturbation with different amplitudes, followed by the study on interactions of such perturbations. These results show the stability of the biaxially strained plate, as long as none of its points has reached the loss of ellipticity condition and determine the actual (nonlinear) zone of influence of the perturbation. The section concludes by investigating the influence of constitutive law and load orientation on the size of the minimum and maximum influence zones.

Constitutive laws

The analysis presented in Section 3.2 is general; any rate-independent constitutive law (which can be put in the form of (3.2.3)) can be accommodated, provided that its membrane (plane stress)

version looses ellipticity at some strain level. Results presented here correspond to the three such models: a hyperelastic (deformation theory) type model of plasticity, the J 2 deformation theory model of Stören and Rice (1975b) and a finite strain generalization of the J 2 flow theory. All models are fitted to the same power law uniaxial stress-strain curve and share the same principal solution. Since no unloading occurs in the perturbed plate prior to reaching a loss of ellipticity, the use of deformation theory type constitutive models is adequate for analyzing its stability.

We start with the hyperelastic constitutive model, which is described by a strain energy W , a function of the equivalent logarithmic strain e as follows:

W = E y 2 1 1 + χ e y χ+1 + 1 2 χ -1 χ + 1 ,    χ = 1 for e ≤ y , χ = n for e > y ,
(3.3.1)

RESULTS

where the equivalent strain e is given in terms of the principal logarithmic strain components α :

e = 2 √ 3 2 1 + 2 2 + 1 2 1/2 ; α = ln λ α , (3.3.2)
with λ i the stretch ratios of the deformation (principal values of the stretch tensor U, the rotationless part of the deformation gradient

F = R • U; U = F T • F 1/2 ).
The above isotropic model is fitted with a piecewise power law uniaxial stress-strain curve4 :

σ e σ y = e y χ , (3.3.3)
where the exponent χ is given in (3.3.1) and the equivalent stress σ e is the Von-Mises stress given in terms of the principal Cauchy stress σ α by:

σ e = σ 2 1 + σ 2 2 -σ 1 σ 2 1/2 . (3.3.4)
Since the principal solution is biaxial straining, the principal stresses are related to the principal logarithmic strains by:

σ α = ∂W ∂ α ; σ 1 = 2 3 E s 2 1 + 2 , σ 2 = 2 3 E s 1 + 2 2 , (3.3.5)
where E s = σ e / e is the secant modulus, E t = dσ e /d e is the tangent modulus of the equivalent uniaxial stress-strain curve in which the equivalent stress and strain are related by: σ e = dW/d e .

When we are no longer along the principal axes of deformation (as is the case of numerical FEM calculations) the stress measures and incremental moduli of this model (see (3.2.3)) are found by:

N αβ = ∂W ∂F αβ , L αβγδ = ∂ 2 W ∂F αβ F γδ . (3.3.6)
In addition to the above-presented hyperelastic constitutive model used in numerical calculations, for comparison purposes two more constitutive models will be employed for the calculation of influence zones under different loading orientations: the J 2 deformation theory model by Stören and Rice (1975b) and the J 2 flow theory model, both in their finite strain version. Details on the incremental moduli derivations of these models from the initial three-dimensional formulation to the two-dimensional plane stress version are presented in the appendix.

As previously mentioned, all three constitutive models share the same uniaxial stress-strain curve and are so constructed as to have the same response when loaded with fixed principal axes of deformation. Since the calculation of δ -and δ + (the minimum and maximum) influence disc sizes for the different constitutive laws requires the principal solution, the evaluation of the corresponding incremental moduli are presented below along the fixed principal axes.

The non-zero components of the plane stress moduli in (3.2.3) are given below in two groups;

the normal moduli components are:

L 1111 = 1 λ 2 1 4 3 E * + E t -E * σ 1 σ e 2 -σ 1 , L 1122 = 1 λ 1 λ 2 2 3 E * + E t -E * σ 1 σ 2 σ 2 e , L 2222 = 1 λ 2 2 4 3 E * + E t -E * σ 2 σ e 2 -σ 2 , (3.3.7)
where for the J 2 deformation theory model as well as the hyperelastic model in (3.3.1), (3.3.2) the normal incremental moduli are the same with E * = E s = σ e / e while for the J 2 flow theory model

E * = E.
The shear moduli components are given by:

L 1212 = 1 λ 2 2 E * 3 + σ 2 -σ 1 2 , L 2121 = 1 λ 2 1 E * 3 + σ 1 -σ 2 2 , L 1221 = L 2112 = 1 λ 1 λ 2 E * 3 - σ 1 + σ 2 2 , (3.3.8)
where for the J 2 flow theory E * = E, for the J 2 deformation theory E * = E s while for the

hyperelasitc model E * = E s [(λ 2 1 + λ 2 2 )/(λ 2 1 -λ 2 2 )](ln λ 1 -ln λ 2 ).
The principal stresses σ α for all three models are identical and given by (3.3.5).

The hyperelastic and J 2 deformation theory models loose ellipicity at realistic strain levels for all load path orientations ψ (see definition (3.2.1)) while the J 2 flow theory gives unrealistic results for load orientations ψ > 0 (and hence the need for the deformation theory models used). Moreover, only values of ψ for which both principal stresses σ i are tensile (σ α = ∂W/∂ α > 0) will be investigated, since a compressible membrane stress is unsustainable (thin plate will immediately buckle).

Evolution at a spatially localized perturbation

We start by analyzing the influence of an initial imperfection, located at the origin, in the form of a square domain of size ∆χ = 6 × 10 -5 whose shear modulus is G(1 + ξ), where ξ = 4 is the imperfection amplitude. Equivalently, one could have taken a perturbation in the principal solution (dimensionless) displacement field u 1 (χ, 0) = χ 1 τ cos ψ, u 2 (χ, 0) = χ 2 τ sin ψ or in the principal solution (dimensionless) velocity field ν 1 (χ, 0) = χ 1 cos ψ, ν 2 (χ, 0) = χ 2 sin ψ. The shear modulus imperfection used here is equivalent to an isotropic displacement or velocity field perturbation that would have resulted from a sudden isotropic dilation/contraction of the perturbed domain at τ = 0. Hence we use the terms initial imperfection or perturbation indistinguishably.

RESULTS

To better visualize the influence of the initial perturbation on the dymamic behavior of the biaxially strained plate, we plot two different measures of the perturbations: the perturbation in the total Green-Lagrange strain:

∆E := E(χ, τ ) -E 0 (χ, τ ) = E αβ -E 0 αβ E αβ -E 0 αβ 1/2 > 0, E αβ = 1 2 F γα F γβ -δ αβ ; E 0 αβ = diag 1 2 λ γ 2 -1 , (3.3.9)
and the perturbation in the shear component of the Green-Lagrange strain, which in view of the absence of a shear strain component in the principal solution, is the shear strain component of the perturbed solution:

∆E 12 := E 12 (χ, τ ) -E 0 12 (χ, τ ) = E 12 > 0. (3.3.10)
The influence of the initial perturbation is shown in Figure 3.6 and Figure 3.7 which depict the evolution of contours of ∆E and ∆E 12 in χ 1 , χ 2 space (only the positive quadrant is shown here due to symmetry). More specifically, the evolution of the strain perturbation due to an initial imperfection of amplitude ξ = 4 (contours of perturbations with magnitudes below 10 -3 are not plotted, for visual enhancement of the evolution of the initial perturbation) for three different dimensionless times τ (0.17, 0.35 and 0.52 = τ m ) are depicted in these figures. Results correspond to a hyperelastic material with n = 0.22 and y = 10 -3 , strained along a loading path with ψ = -1/2, which corresponds to uniaxial stressing for small strains.

For the lowest value of τ = 0.17, only a small region near the origin χ = 0 is affected, while for τ = 0.35 the emergence of a localized band of deformation in the direction of the loss of ellipticity of the material (π/2φ m ) is obvious. At the time of onset of loss of ellipticity, τ m the localized deformation band in the direction π/2φ m is more pronounced. What is worth mentioning is that the localized deformation appears to propagate in three tongues. This phenomenon can be explained by the square shape of the initial perturbation domain, where each corner acts as a source.

A static analogue of this phenomenon has been found in Abeyaratne and Triantafyllidis (1981b).

Moreover the width of the localized deformation zone is considerably larger than the size of the initial perturbation due to propagation of the signal.

Notice that results in Figure 3 These results show that, due to wave propagation, the width of the localized deformation zones are considerable larger than the width of the initial imperfection, but also a fraction of the linearized estimate χ -(influence zone for the slower wave) for the same time, a phenomenon also observed for the growth of a localized perturbation in the nonlinear bar model of [START_REF] Ravi-Chandar | Dynamic stability of a bar under high loading rate: Response to local perturbations[END_REF].

To study the stability of the structure under a spatially localized perturbation, we follow the time (2015), who find stability of spatially localized perturbations at all times prior to reaching the bar's Considère point (maximum force).

The physical meaning of the minimum δ -and maximum δ + influence disc sizes is illustrated in

Size of influence zones for various constitutive laws and loading orientations

The following three figures give the minimum δ -and maximum δ + influence disc sizes as functions of the load orientation angle ψ for the three different plasticity models considered and for three range are terminated when one of the stresses becomes compressive (σ 2 < 0) (applicable for the calculation of δ -) or when a finite value of the influence disc cannot be found (applicable for the calculation of δ + ).

As expected, for a given material and load orientation ψ, both δ -and δ + are increasing functions of the hardening exponent n. There is practically no difference for the minimum influence disc size δ -between the J 2 deformation and hyperelastic theory models over the entire range of load orientations of interest. However the maximum influence disc size δ + predictions for the same two constitutive models coincide only for a certain range of ψ > 0. As the uniaxial strain is approached ψ = 0, the stiffer hyperelastic theory predicts no finite maximum influence disc size δ + , in contrast to the J 2 deformation theory that predicts finite δ + for a significant range of ψ < 0.

There is however a significant difference in the predictions of the much stiffer J 2 flow theory that considerably overestimates δ -over the other two constitutive models (the difference increasing with increasing hardening exponent n) for the range that a reasonable loss of ellipticity strain can

δ -[flow] δ -[hyperelastic] δ -[deformation] δ + [deformation] δ + [hyperelastic]
Figure 3.14: Minimum δ -and maximum δ + influence disc sizes, in dashed and solid lines respectively, as functions of the load orientation angle ψ for the three different constitutive models considered, all sharing the same uniaxial stress-strain curve with y = 0.002, n = 0.1 be found (essentially in the range ψ ≤ 0) and which does not have a δ -for strain paths with ψ ≤ 0,

given that the J 2 flow theory model does not predict loss of ellipticity for these loadings. Also notice that the J 2 flow theory has no finite δ + for any loading.

The difference in the minimum δ -and maximum δ + influence disc sizes predicted by the different constitutive models (in particular between deformation and flow theories) is indicative of the difficulty in predicting failure patterns in these structures and their extreme sensitivity to the constitutive model chosen.

Conclusion

This chapter pertains to the influence of loading rate on the stability of structures when inertia plays a dominant role. The currently established approach to study these stability problems is the method of modal analysis, which determines the structure's fastest growing eigenmode. This method supposes that all points in the structure can be perturbed simultaneously, an assumption that is not appropriate for cases when the velocity of material points in the structure are comparable to the associated wave propagation speeds.

The novel idea here is to analyze the evolution of spatially localized perturbations of the timedependent, high strain rates states of these structures, in order to understand the initiation of the corresponding failure mechanisms. Following the recent analysis by [START_REF] Ravi-Chandar | Dynamic stability of a bar under high loading rate: Response to local perturbations[END_REF] in one-dimensional bars, we study the high strain extension of a two-dimensional, incom-69 Using a finite strain deformation theory of plasticity (based on logarithmic strain), we follow the time evolution of spatially localized perturbations and their interactions. The nonlinear time evolution of such a perturbation is studied numerically using FEM and it is shown that these structures are stable until the time when the condition for the loss of ellipticity is reached. An analytical method, based on linearization, is used to define the size of the influence zone of a point-wise perturbation and we study its dependence on constitutive laws and loading conditions.

3.4. CONCLUSION δ -[flow] δ -[hyperelastic] δ -[deformation] δ + [deformation] δ + [hyperelastic]
The above approach is useful for the stability analysis of more realistic structures under high strain rates. As one such example we cite the recent work by [START_REF] Putelat | Dynamic stability of externally pressurized elastic rings subjected to high rates of loading[END_REF] on the stability of a pressurized thin ring at high rates, where it is shown that for small values of the applied loading rate, the structure fails through a global mode, while for large values of the applied loading rate the structure fails by a localized mode of deformation, as also found recently in the experiments of [START_REF] Mainy | Dynamic buckling of thin metallic rings under external pressure[END_REF]. Our study also shows the sensitivity of the size of minimum and maximum influence zones with respect to the constitutive model used, and hence the caution needed in using such calculations to predict failure patterns.

3.A. PLANE STRESS INCREMENTAL MODULI FOR THE CONSTITUTIVE MODELS USED

3.A Plane stress incremental moduli for the constitutive models used

Three different nonlinear constitutive models are used in this chapter; all adapted for finite strains:

J 2 deformation theory model, J 2 flow theory model and a hyperelastic model, all sharing the same uniaxial stress strain curve, which can be arbitrary and fit to experimental data. The J 2 deformation/flow theory models are incompressible rate-independent (hypoelastic) models that can be put in the from:

σ = C : D -ṗI, (3.A.1)
where ( ) denotes the Jaumann rate of the Cauchy stress tensor ( σ = σ -Ω • σ + σ • Ω with Ω the spin tensor), D is the strain rate tensor and ṗ the hydrostatic pressure rate.

For the case of finite strains, the above current configuration relation can be transformed into its reference configuration counterpart:

Ṡij = L ijkl Ėkl -ṗC -1 ij ; E ij = 1 2 C ij -δ ij , (3.A.2)
where S is the second Piola-Kirchhoff stress, E its work-conjugate Green-Lagrange strain and the reference configuration components of moduli tensor L are:

L ijkl := 2 3 E * 1 2 C -1 ik C -1 jl + C -1 il C -1 jk - 3 2 1 - E t E * S ij S kl σ 2 e - 1 2 C -1 ik S jl + C -1 jk S il + C -1 il S jk + C -1 jl S ik , (3.A.3)
where S is the deviatoric part of the stress tensor S and σ e the Mises equivalent stress, namely:

S ij = S ij - 1 3 C -1 ij C kl S kl , σ 2 e = 3 2 C ik C jl S ij S kl . (3.A.4)
In the above expressions E * = E for the J 2 flow theory while for the J 2 deformation theory of Stören and Rice (1975b) E * = E s = σ e / e is the secant modulus of the uniaxial stress-strain curve.

In both models E t is the tangent modulus of the uniaxial stress-strain curve E t = dσ e /d e .

The principal axes expressions in three dimensions for the equivalent stress σ e and the equivalent strain e , which are useful in (3.3.7), (3.3.8) in the sequel in view of the biaxial loading of the plate are:

σ e = σ 2 1 + σ 2 2 + σ 2 3 -σ 1 σ 2 -σ 2 σ 3 -σ 3 σ 1 1/2 , e = 2 3 2 1 + 2 2 + 2 3 -1 2 -2 3 -3 1 1/2 . (3.A.5)
Due to plane stress loading conditions:

Ṡ3i = 0, Ėα3 = 0, Ė33 = -C 33 C -1 γδ Ėγδ , (3.
Ṡαβ = M αβγδ Ėγδ , M αβγδ = L αβγδ -C 33 L αβ33 C -1 γδ + C -1 αβ L 33γδ + C -1 αβ C -1 γδ C 2 33 L 3333 .
(3.A.7)

The above result, combined with (3.A.3) gives the following expression for M αβγδ :

M αβγδ = 2 3 E * 1 2 C -1 αγ C -1 βδ + C -1 αδ C -1 βγ + C -1 αβ C -1 γδ - 3 2 1 - E t E * S αβ S γδ σ 2 e - 1 2 C -1 αγ S βδ + C -1 βγ S αδ + C -1 αδ S βγ + C -1 βδ S αγ . (3.A.8)
Recalling the relations between S & N and E & F , the moduli L in (3.2.3) are found to be:

L αβγδ = M βζδ F α F γζ + S δβ δ αγ . (3.A.9)
For the case of biaxial loading of interest here F αβ = diag λ γ and thus: A.11) and are based on successive application of chain rule of differentiation using W ( e ) where the equivalent strain e = e (I 1 , I 2 ) is expressed in terms of the invariants of C, which in turn depend on the principal stretch ratios λ α by:

C 11 = λ 2 1 , C 22 = λ 2 2 , C 12 = 0, S 11 = σ 1 /λ 2 1 , S 22 = σ 2 /λ 2 2 , S 12 = 0, (3.
M αβγδ = ∂ 2 W ∂E αβ ∂E γδ = 4 ∂ 2 W ∂C αβ ∂C γδ , (3.
I 1 = tr C = λ 2 1 + λ 2 2 , I 2 = det C = λ 1 λ 2 2 . (3.A.12)
After some lengthy aglebra we end in expression in the expressions given in section 3.3.1.

Chapter 4

Localization of deformation of metallic rings under high loading rate compression

Introduction

As expounded in the previous chapter (Chapter 3), for high loading rates that are commensurate with some characteristic wave propagation speed in the solid/structure, a new approach to analyzing the dynamic stability is needed, namely to study the time evolution of localized perturbation introduced onto the principal solution of the system. In this chapter, in contrast to a solid under rapid extension, of particular interest here is the influence of loading rate on the stability of structures under compression that exhibits an instability even under quasistatic loading. As a model structure to illustrate these ideas, we select an elastoplastic ring subjected to external hydrostatic pressure which is applied at different rates (appropriately non-dimensionalized with respect to elastic axial wave speed). Of course such a classical topic has been treated repeatedly in the mechanics literature; following the work of [START_REF] Carrier | On the vibrations of the rotating ring[END_REF], different linear and nonlinear versions of the ring dynamical equations of increasing complexity have been proposed (e.g. [START_REF] Boresi | Energy methods in parametric excitation of rings[END_REF][START_REF] Dempsey | Dynamic nonlinear forcing of elastic rings[END_REF][START_REF] Goodier | The elastic cylindrical shell under nearly uniform radial pressure[END_REF][START_REF] Graff | On dispersion of elastic waves in rings[END_REF][START_REF] Morley | Elastic wave in a naturally curved rod[END_REF][START_REF] Simmonds | Accurate nonlinear equations and a perturbation solution for the free vibrations of a circular elastic ring[END_REF][START_REF] Wah | Dynamic buckling of thin circular rings[END_REF] to study their vibrations. The stability of rings subjected to impulsive or step loadings has also been repeatedly studied (e.g. [START_REF] Amabili | Review of studies on geometrically nonlinear vibrations and dynamics of circular cylindrical shells and panels, with and without fluid-structure interaction[END_REF][START_REF] Anderson | Dynamic pulse buckling of cylindrical shells under transient lateral pressures[END_REF][START_REF] Florence | Buckling of viscoplastic cylindrical shells due to impulsive loading[END_REF][START_REF] Goodier | The elastic cylindrical shell under nearly uniform radial pressure[END_REF][START_REF] Lindberg | Buckling of a very thin cylindrical shell due to an impulsive pressure[END_REF][START_REF] Lindberg | Stress amplification in a ring caused by dynamic instability[END_REF][START_REF] Lindberg | Dynamic pulse buckling[END_REF][START_REF] Simmonds | Accurate nonlinear equations and a perturbation solution for the free vibrations of a circular elastic ring[END_REF]. These studies rely on modal analysis using Fourier series whose truncation leads nonlinear amplitude equations and showed that dynamic buckling is triggered by flexural modes.

Problem setting

We consider a homogeneous elastoplastic ring of rectangular section with thickness h, width a and cross-sectional area A = h × a. The ring has a mid-line radius r and follows small strain-moderate rotation Timoshenko kinematics described by ṽ(θ), w(θ), and ψ(θ) respectively the tangential and normal displacements of the ring's reference mid-line at point θ and the rotation of the section perpendicular to the mid-line, initially at θ.

Before deriving the governing equations for the system, we would want to introduce a set of useful dimensionless variables and physically relevant parameters. Length, time, and stress are nondimensionalized by r, r/c and G where r is the radius of the ring in initial reposed configuration, c = G/ρ is proportional to the ring's longitudinal wave speed. In addition, the slenderness parameter η and the dimensionless pressure λ are also proven to be expressive and listed below:

η := I/ Ar 2 = (h/r) 2 /12; λ := λra/ (GA) = λ/G (r/h) . (4.2.1)
To find the system's Lagrangian, We need to determine its potential and kinetic energies P and K respectively. The potential energy P consists of two parts: the stored strain energy P int plus P ext the work potential of the externally applied uniform pressure lambda. All the quantities are expressed with before-mentioned non-dimensionlization in length, time, and stress.

The stored strain energy P int is

P int = 2π 0 h/2 -h/2 W ( θθ ) + 1 2 Gχ γ 2 rθ dz ar dθ (4.2.2)
where the axial and shear strains θθ and γ rθ are given by

θθ = v + w + 1 2 v -w 2 + ψ z r , γ rθ = v -w -ψ (4.2.3)
with f denoting the θ-derivative of the corresponding function, W the axial strain energy, χ the shear correction factor (χ = 2/3 for a rectangular section) and G the material's shear modulus.

By reformulating the integration in thickness and collecting all the dimensional terms, the following expression is obtained for the internal energy

P int = Garh 2π 0 1/2 -1/2 W e + φ 2 2 + ζ h r ψ dζ + 1 2 χ(φ -ψ) 2 dθ (4.2.4)
where e := v + w, φ := vw , and W := W /G denote respectively the dimensionless axial strain, rotation of the ring's mid-line, and the dimensionless strain energy.

The work potential P ext of the external pressure loading λ applied on the ring equals λ ∆S

where ∆S is the change of area due to deformation (ṽ, w) enclosed by the ring's mid-line, which is given by (e.g. [START_REF] Brush | Buckling of bars, plates, and shells[END_REF])

P ext = Garh λ 2π 0 w + 1 2 v 2 -vw + w v + w 2 dθ = Garh λ 2π 0 w + 1 2 (vφ + we) dθ (4.2.5)
where λ is taken positive when acting inwards (resulting in compressive hoop stress σ θθ < 0) in the ring.

Using the same kinematic assumptions as for the derivation of (4.2.3), the kinetic energy of the ring is

K = Garh 2π 0 1 2 v2 + ẇ2 + η ψ2 dθ (4.2.6)
where η is defined in (4.2.1) and also served as slenderness parameter.

The system's Lagrangian is now determined by (4.2.4), (4.2.5), and (4.2.6)

L = K -P = K -(P int + P int ) . ( 4
          v = σ -σφ + χ(ψ -φ) -λφ, ẅ = -σ + χ ψ -φ -λ (1 + e) -(σφ) , η ψ = h r τ -χ(ψ -φ) (4.2.8) 4.2. THEORY
where σ and τ denote respectively the dimensionless axial force and the dimensionless bending moment defined by

σ := 1/2 -1/2 dW d dζ, τ := 1/2 -1/2 ζ dW d dζ, (4.2.9)
To the above equations initial conditions for v, w, ψ and v, ẇ, ψ must be added.

Principal solution

Of particular interest here now is the perfect structure's principal solution (v 0 = 0, w 0 , ψ 0 = 0), i.e.

the response of the perfect ring to a uniform pressure loading at constant rate (starting at t = 0): λ =t. Due to axisymmetry, there is zero tangential displacement and zero rotation of section, the governing equation of axial displacement w 0 (t) is simplified to:

ẅ = -σ -λ(1 + w), (4.2.10)
with inital conditions w 0 (0) = ẇ0 (0) = 0.

For the quasistatic case with elastic material, we have σ = kw and w 1 with k = 2(1 + ν).

By using the fact that dt 2 = -2 dλ 2 , the previous equation can be written as

2 d 2 w dλ 2 = kw -λ (4.2.11)
with corresponding initial conditions, which admits solution [START_REF] Putelat | Dynamic stability of externally pressurized elastic rings subjected to high rates of loading[END_REF])

w(λ) = √ k sin √ k λ -λ k (4.2.12)
For the case where inertia effect is taken into account, we have w ∼ 1, and the force σ could present substantial nonlinearity regarding w. Due to these high nonlinearities, a numerical approach is adopted to calculate its solution using given initial conditions. By virtue of the axisymmetry of the principal solution, this numerical simulation could be very computationally light-weighted by using axisymmetric elements to mesh only the ring's rectangular cross-section. The same set of numerical results, used for computation of the principal solutions, will also serve in numerical stability analysis reported in Section 4.3.

Linearized stability analysis

Assuming a perturbation (∆v, ∆w, ∆ψ) superposed on the principal solution (0, w 0 , 0), the linearized perturbation equations obtained from (4.2.8) are where σ 0 and E 0 t denote respectively the axial force and the tangential modulus, both computed with the principal solution:

         ∆v = E 0 t ∆v + (E 0 t + σ 0 + χ + λ)∆w -(λ + σ 0 + χ)∆v + χ∆ψ, ∆ ẅ = -(E 0 t + σ 0 + χ + λ)∆v + (χ + σ 0 )∆w -(E 0 t + λ)∆w + χ∆ψ , η ∆ ψ = ηE 0 t ∆ψ -χ(∆ψ -∆v + ∆w ) (4.2.13) (a) Elastic (b) Bilinear Et/E = 0.22 (
σ 0 := 1/2 -1/2 dW d 0 dζ, E 0 t := 1/2 -1/2 d 2 W d 2 0 dζ. (4.2.14)
The equation (4.2.13) is completed by the initial conditions for (∆v, ∆w, ∆ψ) and their time derivatives at the time of the onset of perturbation t = t 0 .

Using the method of frozen coefficients, i.e. assuming that the rate of growth/decay of the perturbation is much higher than the loading rate, (4.2.13) is considered as a constant coefficient equation, which admits a base of solutions of form

(∆v, ∆w, ∆ψ) = (v(k), ŵ(k), ψ(k)) exp i(ωt + kθ) ; ω ∈ R, k ∈ N, k ≥ 2. ( 4
          χ + λ + σ 0 + E 0 t k 2 -ω 2 -ik E 0 t + σ 0 + χ + λ -χ ik E 0 t + σ 0 + χ + λ E 0 t + λ + (χ + σ 0 )k 2 -ω 2 -ikχ -χ ikχ χ + η E 0 t k 2 -ω 2           = 0 (4.2.16)
which yields the following bi-cubic polynomial in ω 2 a 0 + a 2 ω 2 + a 4 ω 4 + a 6 ω 6 = 0, (4.2.17) 4.3. NUMERICAL SIMULATIONS with coefficients

a 0 = -χ(E 0 t + λ)(λ + σ 0 ) + σ 0 χ -η(E 0 t + λ)(λ + σ 0 + χ) E 0 t k 2 + η(σ 0 + χ) E 0 t 2 k 4 k 2 -1 a 2 = -χ(E 0 t + 2λ + σ 0 ) -η(E 0 t + λ)(λ + σ 0 + χ) + ηλ 2 -η(E 0 t + λ)(E 0 t -σ 0 -χ) -χ(E 0 t + σ 0 ) k 2 -ηE 0 t (E 0 t + 2σ 0 + 2χ) k 4 a 4 = χ + η(E 0 t + 2λ + σ 0 + χ) + η(2E 0 t + σ 0 + χ) k 2 a 6 = -η (4.2.18)
It should be noted that for k = 1 the dispersion equation has a zero root ω 2 = 0, reflecting the fact that k = 1 corresponds to a rigid body mode of the ring, easily verified since the corresponding strain measure vanishes.

We can find out the time corresponding to the onset of static buckling by calculating the lowest time t b required for a zero root ω 2 , i.e. solving a 0 = 0 in (4.2.18). Using the fact that up to static buckling the thin ring always has λ 1, σ 0 1, and λ ≈ -σ 0 , a 0 = 0 can be simplified into

λ = η E 0 t k 2 -1 (4.2.19)
which corresponds to the buckling pressure of a quasistatically loaded ring, achieved for the lowest integer value of k = 2 [START_REF] Brush | Buckling of bars, plates, and shells[END_REF].

For finite loading rates inertia effects are important and perturbations travel at finite speeds. For such cases, failure occurs by a localized deformation mode, which correspond to short wavelength k 1. The time corresponding to their onset of instability can be found by investigating the behavior of the dispersion relation (4.2.17) for different values of k, when the lowest root ω 2 = 0, which is equivalent to a 0 = 0. Under stability theory framework, all types of infinitesimal deviation from the ideal system are "equivalent" in the sense that they can equally trigger instability (if there is any) under the same loading condition intrinsic to the ideal system and then draw corresponding bifurcated solutions. In the real experiments of ring compression, they can be incarnated into various forms:

voids/inclusions inside aluminum, machining error, non-uniformity of the external fields, etc. In numerical simulations, these can be modeled by imperfections manually introduced on the level of elements ("imperfection"), or direct modification of the converged solution between any two consecutive time-integration iterations ("perturbation"). In all the following calculations, both the "imperfection" type and "perturbation" type of modification are used for the ring compression system. More precisely, for the "imperfection" type, small radial shifts ∆w (θ) is introduced onto the mid-plane of the ring as part of the reposed configuration, and boundary surfaces are defined as in/outward equidistant offset of the mid-plane by h/2 (Figure 4.3b), where h is the thickness of the ring used in (4.2.1). Such geometric imperfection imitates roughly the machining error of the ring. For the "perturbation" type, small radial velocity ∆ ẇ(θ, t 0 ) is added to the principal solution at a chosen time t = t 0 . Such fluctuation of the solution imitates roughly the inhomogeneity of the external fields. When no confusion is possible, we refer both types of "small modification" of the system as imperfections and specify the first type as "geometrical imperfection" when individual clarification is needed.

Imperfections of various shapes can be implemented by different amplitudes a (θ), of interest here is two specific kinds of imperfection, namely the spatially isolated imperfection and randomly distributed imperfections, as illustrated in As read from Figure 4.4b, the result of discrete Fourier transform shows that the imperfection amplitude a (θ) introduced by these randomly distributed imperfections blend in a wide range of periodic shapes without promoting a significant dominant wavelength, i.e. it cannot be directly approximated by a sinusoidal function. One could expect that this random profile is free from bias for any unstable periodic mode.

Results

In the case of linearly time-varying loading case, the external surface of the ring is subject to a uniform pressure field, with time-varying amplitude λ = t where is the dimensionless loading rate. By the dimensionless loading rate alone one can freely parameterize the system around the quasistatic regime and the highly dynamic regime. This freedom marks a great difference with the case of two-dimensional extension of metallic sheets, as in the latter case the "free" parameter-the extent of the sheet-can be treated as infinite in the context of spatially localized perturbation. is An interesting empirical practice is proposed by [START_REF] Putelat | Dynamic stability of externally pressurized elastic rings subjected to high rates of loading[END_REF]. In these numerical experiments, they observed that the axial component (as opposed to bending and shearing components) of internal energy stored in the elastic ring first increased with the load, then followed by a steep drop once the absolute maximum has been reached. And more importantly, this maximum point is proven to be a robust indicator to detect the onset of instability for both the local and global modes. Inspired by this observation, we propose a similar approach to serve as an indicator of the onset of instability for elastoplastic rings by averaging mid the axial strain along the mid-plane (i.e.

the red lines highlighted in Figure 4.3) and calculating its minimum within time.

is due to the fact that no collision detection is set up for these calculations, and some segments of the ring eventually intersect themselves which ends in an inversion of inner/outer surfaces. One can roughly divide the figures of six different loading rate into two groups according to the ratio between t l the time of onset of general instability and t b the time of the onset of quasistatic buckling instability, with 4.5a, 4.5b, and 4.5c in the first group for t l /t b 1, 4.5e and 4.5f in the second group for t l ∼ t b , and 4.5d entering into either group as a transitional element. As expected, the imperfect ring loaded at a high rate is initially stable for pressures well above the static buckling, which is calculated by ignoring the speed of propagation of the perturbation. The rings in 4.5a and 4.5b have their t l strictly less than than this lower bound, which implies that the influence generated by their isolated/distributed imperfection(s) never has enough time to affect the whole structure. Whereas for 4.5e and 4.5f, the influence would have totally disseminated over the whole structure.

For all the six figures, especially in the first four: 4.5a, 4.5b, 4.5c, and 4.5d where local effect overwhelms global effect, The t l and the corresponding strain mid (t l ) all have a good match between two extensively different imperfection shapes (isolated/distributed). This proves reciprocally the robustness of this indicator, i.e. using the axial strain limit to predict the onset of instability, for both quasistatic and high loading rate cases.

The influence of the time t 0 at the introduction of isolated perturbation in velocity and the meaning of the localization time t l introduced in (4.3.5) is illustrated in Fig. 4.6 for two case of high loading rate ( = 10 -2 and 10 -4 ) for a ring with slenderness η = 10 -4 . We plot the norm of the perturbation ∆ ẇ as a function of time for perturbations with different introduction time t 0 . Notice that t l is an upper bound for the real time of onset of localization. The perturbation amplitude first decreases provided that the time of perturbation introduction is less than t l , which is a direct consequence of the stability of the perturbed solution before onset of localization.

Once the t l is well defined through the results of calculations, it's possible to investigate the shape and patterns formed by the metallic ring in post-bifurcation phase. In We first observe that independently of the loading rate, the initial disturbance caused by the introduced imperfection splits in two. As the loading rate ranges from the smallest = 10 -6 to the In the regime where localization of deformation is in favor, we observe a stationary wave packet whose maximum amplitude of its envelope grows with time, and that compared to elastic ring case [START_REF] Putelat | Dynamic stability of externally pressurized elastic rings subjected to high rates of loading[END_REF], the localizations on the elastoplastic ring has a much faster spatial decay outside of the two principal localizations. For the latter one, an interesting analog can be made with the calculation of stretching of a two-dimensional elastoplastic sheet reported in the previous chapter.

In profile is traced according to its mid-plane, evolving from outmost (the pure blue line) inwardly (purple lines than the pure red line). Green lines correspond to the onset of localization t l . The figures are ordered from high loading rate to low loading rate.

t s . In Figure 4.8, it could imply that the strong short-range effect has somehow "locked" distributed imperfections' information and force them to express "in place".

Conclusion

As a model structure to study the influence of inertia and loading rates on the stability of a structure that becomes unstable even at static loads (structures with buckling modes), we study an elastoplastic ring subjected to external hydrostatic pressure applied at different rates (appropriately non-dimensionalized with respect to elastic axial wave speed). Unlike existing analyses of this the extremum of the average axial strain is reached. This chapter presents a new approach for investigating the dynamic stability of structures that exhibit instabilities even under static loadings. For a future perspective, the influence of different types of dynamic loading (i.e. step, pulse) needs to be considered as well as a more realistic modeling of the electromagnetic loading conditions.

Chapter 5

Concluding remarks

This thesis pertains to the study of localization of deformation in bulk of solids and in structures under quasistatic and dynamic loading conditions. The novel approach of this work consists of introducing geometrically localized perturbations and investigating their time evolution to study failure by localization of deformation in the corresponding problems. Analytical as well as numerical (FEM) tools are used and the results are interpreted in light of experiments available in the literature.

Following the brief introduction, the problem of localization of deformation in microstructured solids under static loading condition is examined in Chapter 2. The specific problem discussed is the stability of uniform compressive deformation in a unidirectionally oriented fiber reinforced composite. This thesis refutes the common belief that there is always a localized deformation appearing in the post-critical equilibrium path of solids loosing macroscopic ellipticity. It's proposed in this

Chapter that the key concept for finding whether loss of macroscopic ellipticity leads to localization of deformation lies in the post-bifurcation behavior of the solid under investigation. A phase diagram has been created that outlines the various kinds of instabilities that arise in this problem for cases when the critical bifurcation mode is also global in nature and well separated from other eigenmodes.

In Chapter 3 and Chapter 4, two problems involving dynamic instability are presented. The main objective is to explore the influence of high rate of background deformation on the growth of a perturbation arising from a point defect. This thesis refutes the prevalent approach of using the method of modal analysis, which determines the structure's fastest growing eigenmode, to study its stability. Modal analysis method supposes that all points in the structure can be perturbed simultaneously, an assumption that is not appropriate for cases when the velocity of material points in the structure are comparable to the associated wave propagation speeds. The novel idea proposed CHAPTER 5. CONCLUDING REMARKS 90 in these two Chapters is to analyze the evolution of spatially localized perturbations of the timedependent, high strain rates states of these structures, in order to understand the initiation of the corresponding localization of deformation.

In Chapter 3 where biaxial tension incorporating dynamics/inertial effects is considered, it is shown that the deformation under biaxial tension is dynamically stable until loss of ellipticity of material is reached. Moreover, we found that under high strain rates the plate has not shown any localized deformation mode for strains well above the ones corresponding to the static loss of ellipticity of material, thus provides explanation for ductility increase observed in experiments.

In Chapter 4 where an analysis of localization of deformation in thin metallic rings under high strain-rate compression is presented, it is demonstrated that for high strain rates the failure pattern involves highly localized deformation zones and a resulting "irregular flower pattern" indicates interactions of randomly existing imperfections on structures as observed in experiments.

  Fig. 1.1b, désignant l'apparition de la striction dans un tube en aluminium en expansion électromagnétiques. Elle peut se trouver au sein des matériaux et on parle alors d'instabilité matérielle, ou sur la structure entière et dans ce cas il s'agit d'une instabilité géométrique. Les exemples sont présentés pour le cas de chargement quasistatique dans Fig. 1.1c, décrivant le pliage localisées dans des tubes carré de métal sous l'écrasement axial et pour le chargement dynamique en Fig. 1.1d, montrant la localisation de la déformation dans un anneau d'aluminium électromagnétiquement compressé. Cette thèse étudie le phénomène de localisation de déformation dans des contextes matériaux ou géométriques et avec des conditions de chargement statique ou dynamique. Dans tous les cas, un outil unifié est utilisé : l'évolution de la perturbation à support localisé. La classification des sujets selon la condition de chargement mécanique (statique vs. dynamique) et la mécanisme de l'instabilité (matièrielle vs. structurelle) est présentée dans le tableau figurant dans la partie gauche de la Fig. 1.2, où chaque cellule contient un expriment représentative. Le tableau donnant le chapitre correspondant dans cette thèse est donné dans la partie droite de la Fig. 1.2.

1. 1 CHAPTER

 1 Fig. 2.8. On the left is the macroscopic stress-strain response of the homogenized perfect composite (in red) and of its finite-size imperfect counterpart computed by FEM (in black). On the middle and right are deformed configurations at points a and b, showing also contours of the Lagrangian shear strain component E 12 . Undeformed configuration is shown by its bounding square. . . . . . . . . . . . . . . . 26 2.10 Response of a composite with a stable homogenized post-bifurcation behavior Λ 2 > 0, λ 2 > 0. The composite consists of two equal thickness layers: a neo-Hookean (n = 1) and a softer layer (n = 0.75, α = 1.0, see point B in Fig. 2.8). On the left is the macroscopic stress-strain response of the homogenized perfect composite (in red) and of its finite-size imperfect counterpart computed by FEM (in black). On the middle and right are depicted the deformed configurations at points a and b, showing also contours of the Lagrangian shear strain component E 12 . Undeformed configuration is shown by its bounding square. . . . . . . . . . . . . . . . . . . . . 27

Figure 2 . 2 :

 22 Figure 2.2: Reference configuration with a unit cell in (a) and bifurcation eigenmode type (local or global) in (b), for axially compressed layered solid deformed under plane strain conditions.

  -bifurcation equilibrium for global critical mode: exact solution By ignoring the influence of boundary conditions, the onset and evolution of the long wavelength bifurcated solution emerging from λ c = λ H is depicted in Figure 2.3, according to which the initially parallel and straight layers in the principal solution in Figure 2.3a evolve in a long-wavelength wavy pattern shown in Figure 2.3b. The post-bifurcation equilibrium solution is idealized in Figure 2.3c,

Figure 2 . 3 :

 23 Figure 2.3: Post-buckling mechanism for the case of a global critical mode in axially compressed layered media subjected to a macroscopic compressive stress f = -Λ with a corresponding macroscopic strain δ = -λ.

Figure 2

 2 Figure 2.4: Unit cell for the post-buckling equilibrium solution of an axially compressed composite with a global critical mode.

2

  Henceforth g := ζ f g f + ζmg m denotes the weighted average of a function g with different values of the fiber and matrix part of the composite.

  where ζ f := H f /H and ζ m := H m /H are respectively the volume fractions of the fiber and matrix.

Figure 2 . 5 :

 25 Figure 2.5: Dimensionless macroscopic axial stress Λ/µ G vs. its work-conjugate strain λ for the principal and bifurcated equilibrium path of perfect neo-Hookean composites with µ f /µ m = 10 (left) and µ f /µ m = 100 (right) calculated for four different volume fractions.

Figure 2 . 6 :

 26 Figure 2.6: Influence of fiber volume fraction on the critical strain for two different fiber-to-matrix stiffness ratios (left) and influence of the fiber-to-matrix stiffness contrast on the critical strain for four different volume fractions (right).

  Figure 2.7: Influence of fiber volume fraction on the critical stress for two different fiber-to-matrix stiffness ratios (left) and influence of the fiber-to-matrix stiffness contrast on the critical stress for four different volume fractions (right).

Figure 2

 2 Figure 2.8: Left: Uniaxial response in simple shear (τ -γ) of the neo-Hookean fiber (linear, in black) and of the matrix (in color, for different values of the hardening exponent n) of a composite containing two equal thickness layers. Right: Initial post-bifurcation stability behavior of this composite at a given critical load (λ c , in solid black lines) as a function of matrix constitutive parameters α, n. Area shaded red: Λ 2 > 0, λ 2 > 0, area shaded green: Λ 2 < 0, λ 2 > 0, remaining non-shaded area: Λ 2 < 0, λ 2 < 0. The blue shaded area at the bottom of the graph indicates composites where a local bucking mode precedes the global one.

F

  12 = 0 and a macroscopic (average) first Piola-Kirchhoff stress with Π 12 = Π 22 = 0, in agreement with (2.2.8). Moreover, periodic boundary conditions are imposed on opposite faces of this square segment, so that the perfect, finite-size composite can exhibit both the perfect principal and perfect bifurcated solutions presented in Subsection 2.2.2, thus avoiding the development of large deviations from the perfect solution at the boundaries.Several such finite-size composites are investigated; the results presented here correspond to an initially square segment containing 40 unit cells, with each such cell consisting of a fiber layer (initial dimension L/80 × L and a matrix layer (initial dimension L/80 × L; only composites with 2.3. RESULTS equal matrix and fiber volume fractions (ζ m = ζ f = 0.5) are modeled. Each fiber and matrix layer are discretized using a Cartesian grid of 4 quadrilateral elements through the thickness and 320 elements along their length for a grid consisting of a total of 102,400 such elements. The FEM calculations are performed using ABAQUS by implementing the constitutive relations in(2.3.16) into the incompressible element CPE4H.

Figure 2

 2 Figure 2.9: Response of a neo-Hookean composite, exhibiting a stable homogenized postbifurcation behavior (Λ 2 > 0, λ 2 > 0). The composite consists of two equal thickness neo-Hookean (n = 1) layers with a stiffness ratio α = µ m /µ f = 0.22, see point A in Fig. 2.8. On the left is the macroscopic stress-strain response of the homogenized perfect composite (in red) and of its finitesize imperfect counterpart computed by FEM (in black). On the middle and right are deformed configurations at points a and b, showing also contours of the Lagrangian shear strain component E 12 . Undeformed configuration is shown by its bounding square.

  Figure 2.10: Response of a composite with a stable homogenized post-bifurcation behavior Λ 2 > 0, λ 2 > 0. The composite consists of two equal thickness layers: a neo-Hookean (n = 1) and a softer layer (n = 0.75, α = 1.0, see point B in Fig. 2.8). On the left is the macroscopic stress-strain response of the homogenized perfect composite (in red) and of its finite-size imperfect counterpart computed by FEM (in black). On the middle and right are depicted the deformed configurations at points a and b, showing also contours of the Lagrangian shear strain component E 12 . Undeformed configuration is shown by its bounding square.

  (2.A.4). Here again, the superscripts m and f are omitted from (2.A.8) in the interest of notational simplicity since the components of V, Z and C 1 are evaluated on the corresponding layer. Substituting equation (2.A.4) into the Floquet conditions (2.A.3) results in the additional relation:

  m

Figure 2

 2 Figure 2.12: Different types of bifurcation modes in axially compressed layered media.

σ

  denotes the Jaumann rate of the Cauchy stress, D the strain rate tensor, ṗ the hydrostatic pressure rate and the quantities µ and µ are the incremental moduli associated with an infinitesimal simple shear, superposed on a homogeneous deformation, parallel to the principal axes and at π/4

Figure 2 .

 2 Figure 2.13: Critical load (λ c , in solid black lines) as a function of matrix constitutive parameters for a composite containing two equal thickness layers: a stiff neo-Hookean fiber and a soft nonlinear matrix. On the left figure the matrix is hyperelastic while on the right figure the matrix is based on a deformation theory model of elastoplasticity, with both models sharing the same uniaxial stressstrain response. The blue shaded area in each graph indicates composites where a local bucking mode precedes the global one.

Figure 2 .

 2 Figure 2.14: Influence of constitutive law choice on the critical strain λ c as a function of the initial matrix-to-fiber stiffness ratio α for the hyperelastic (red line) and the hypoelastic (deformation theory, blue line) matrix models presented in Fig. 2.13, for three different matrix hardening exponents: n = 0.55, 0.57, 1. A solid line indicates a global critical (lowest strain) mode and a dotted line indicates a local one.

Figure 3 . 1 :

 31 Figure 3.1: Unfolded conical mirror image for an electromagnetically expanding Al 6061-O tube test (from Zhang and Ravi-Chandar (2010)), showing the onset and evolution of necks under high strain rate loading; notice the absence of a dominant wavelength in the failure pattern.

Figure 3

 3 Figure 3.2: A schematic diagram of a biaxially strained plate.

  .2.3) where L αβγδ are the plane stress incremental moduli of the plate, which in general depend on the current state of stress plus the deformation history represented by a set of internal variables. This moduli obtained from the three-dimensional version of the constitutive equation Ṅij = L ijkl Ḟkl (relating the rate of the norminal stress Ṅij to its work-conjugate quantity Ḟkl ) plus the plane stress condition Ṅi3 = 0 and the orthotropy of the plate with respect to the thickness direction. For the case of a hyperelastic material the stress measures are derivable from a potential i.e. N αβ = ∂W/∂F αβ in which case L αβγδ = ∂ 2 W/∂F αβ ∂F γδ where W (F) is the two-dimensional strain energy density of the plate.

  Figure 3.3: Typical influence cone of an intial (τ = 0) perturbation at X = 0 showing the evolution of the influence zones χ -(red), χ + (blue) of the slow (ν -) and fast (ν + ) wave speeds as a function of time τ .

Figure 3

 3 Figure 3.4: Influence zones of an initial (τ = 0) perturbation at X = 0 for biaxially stretched thin, flat plates: χ -(φ, τ m ), χ + (φ, τ M ) in solid lines. The minimum δ -and maximum δ + discs of influence are depicted in dotted line. Also plotted in dashed line is χ e (φ), the locus of points reached by the wave propagating along a given direction φ until loss of ellipticity occurs. Results correspond to uniaxial strain ψ = 0 of a power-law type material with hardening exponent n = 0.22.

  along with their orthogonal eigenvectors, noted as e + and e -: A = α + e + e + + α -e -e - (3.2.25) Since e + and e -are orthogonal, one can define a real value β such that ė+ = βe -, ė-= -βe + (3.2.26) and the time derivative of A can be given accordingly: Ȧ = α+ e + e + + αe -e -+ β (α +α -) (e + e -+ e -e + ) . (3.2.27) Using α + , α -, and β and their time derivatives, we can transform (3.2.24) into a equivalent yet more explicit form:

λ

  Figure 3.5: Influence of the time of perturbation introduction on the subsequent growth/decay of the perturbation. The vertical (blue) line indicates the localization time τ m , which for the parameter values used here is τ m = 0.52.

  Figure 3.6: Green-Lagrange strain perturbation contours ∆E(χ, τ ) at three different dimensionless times a) τ = 0.17, b) τ = 0.35 and c) τ = 0.52 = τ m , where only contours of ∆E ≥ 10 -3 are shown in color. The extent of the influence zones χ -(φ, τ ) and χ + (φ, τ ) for the slowest and fastest wave speeds ν -and ν + respectively, are also shown in these figures. Results calculated correspond to a hyperelastic constitutive law with a piecewise power law uniaxial curve ( y = 0.002, n = 0.22) and a loading angle tan ψ = -1/2.

Figure 3 . 7 :

 37 Figure 3.7: Shear strain perturbation contours ∆E 12 (χ, τ ) at three different dimensionless times a) τ = 0.17, b) τ = 0.35 and c) τ = 0.52 = τ m , where only contours of ∆E ≥ 10 -3 are shown in color. The extent of the influence zones χ -(φ, τ ) and χ + (φ, τ ) for the slowest and fastest wave speeds ν -and ν + respectively, are also shown in these figures. Results calculated correspond to a hyperelastic constitutive law with a piecewise power law uniaxial curve ( y = 0.002, n = 0.22) and a loading angle tan ψ = -1/2.

  Figure 3.8: Green-Lagrange strain perturbation ∆E plotted in polar coordinates (χ 1 , χ 2 ) = (ρ cos θ, ρ sin θ) for θ ∈ [0, π/2] and at different distances ρ from the origin (distance is increasing by constant ∆ρ = 1.5×10 -3 and ρ+∆E(ρ, θ) is plotted in the y-axis). Results are shown for three different dimensionless times a) τ = 0.17, b) τ = 0.35 and c) τ = 0.52 = τ m . The extent of the influence zones χ -(τ ) and χ + (τ ) for the slowest and fastest wave speeds ν -and ν + respectively are also shown in these figures. Result calculated correspond to a hyperelastic constitutive law with a piecewise power law uniaxial curve ( y = 0.002, n = 0.22) and a loading angle tan ψ = -1/2.

  Figure 3.9: Shear strain perturbation ∆E 1 2 plotted in polar coordinates (χ 1 , χ 2 ) = (ρ cos θ, ρ sin θ) for θ ∈ [0, π/2] and at different distances ρ from the origin (distance is increasing by constant ∆ρ = 1.5 × 10 -3 and ρ + ∆E(ρ, θ) is plotted in the y-axis). Results are shown for three different dimensionless times a) τ = 0.17, b) τ = 0.35 and c) τ = 0.52 = τ m . The extent of the influence zones χ -(τ ) and χ + (τ ) for the slowest and fastest wave speeds ν -and ν + respectively are also shown in these figures. Result calculated correspond to a hyperelastic constitutive law with a piecewise power law uniaxial curve ( y = 0.002, n = 0.22) and a loading angle tan ψ = -1/2.

Figure 3

 3 Figure 3.10: Profile of Green-Lagrange strain perturbation ∆E plotted at different dimensionless times 0 ≤ τ ≤ τ m (in increments of ∆τ = 0.02 ) at a distance ρ from the origin and for two different values of polar angle a) θ = φ m and b) θ = π/2φ m . The extent of the influence zones χ -(τ ) and χ + (τ ) for the slowest and fastest wave speeds ν -and ν + respectively are also shown in these figures. Result calculated correspond to a hyperelastic constitutive law with a piecewise power law uniaxial curve ( y = 0.002, n = 0.22) and a loading angle tan ψ = -1/2.

  Figure 3.12: Influence of the initial amplitude ξ of a localized perturbation at χ = 0. Notice that the magnitude of the perturbation (measured by its maximum, at a given time τ , over the entire plate, i.e. ∆E max (τ ) := max χ∈R 2 ∆E(χ, τ ) in (a) and ∆E max 12 (τ ) := max χ∈R 2 ∆E 12 (χ, τ ) in (b)) decreases, for each value of τ < τ m , with decreasing ξ, thus showing the stability of the structure as long as it stays in the elliptic domain (non-shaded area τ < τ m in the graphs).

Figure 3 .

 3 Figure3.13 that shows contours of strain perturbation ∆E ≥ 10 -3 at the time of loss of ellipticity τ m in a plate with two localized imperfections of the same size and initial amplitude spaced at a distance smaller than δ -in a) and at distance larger than δ -in b).It appears from Figure3.13a, that when the localized deformation zones of the two perturbations meet and interact, the width of the resulting localized deformation zone is bigger than the width of the single localized imperfection. The failure pattern for the stretched plate can be explained as resulting from interaction of statistically distributed such localized defects-inevitable in realityas observed experimentally in the tube expansion experiments of[START_REF] Zhang | On the dynamics of necking and fragmentation -iv. expansion of al 6061-o tubes[END_REF]

  different power-law hardening exponents, Figure 3.14 for n = 0.1, Figure 3.15 for n = 0.22 (typical of Al alloys) and Figure 3.16 for n = 0.40 (typical of steel alloys). Curves in the ψ < 0

Figure 3 .

 3 Figure 3.15: Minimum δ -and maximum δ + influence disc sizes, in dashed and solid lines respectively, as functions of the load orientation angle ψ for the three different constitutive models considered, all sharing the same uniaxial stress-strain curve with y = 0.002, n = 0.22

Figure 3 .

 3 Figure 3.16: Minimum δ -and maximum δ + influence disc sizes, in dashed and solid lines respectively, as functions of the load orientation angle ψ for the three different constitutive models considered, all sharing the same uniaxial stress-strain curve with y = 0.002, n = 0.40

  A.10) which upon substitution into (3.A.8), (3.A.9) gives the incremental moduli expressions in (3.3.7), (3.3.8) for J 2 flow & J 2 deformation theories.Calculations of the moduli for the hyperelastic model also use (3.A.9) & (3.A.10) but M is derived from the strain energy potential W through:

  Figure 4.1: Overlay of sequential images of a radially compressed ring electromagnetically compressed with a charge level of 3 kV. Images are captured between t = 0 and 91.8 µs (Mainy, 2012). Courtesy of Prof. K. Ravi-Chandar, University of Texas Austin.

  .2.7) Using now Hamilton's principle, i.e. by extremizing the action integral T 0 L dt over time paths with fixed initial and final time values of the independent variables v, w, and ψ, we deduce the following Euler-Lagrange equations governing respectively the axial, normal and rotational motion of the ring 

  Figure 4.2: Radial displacement versus external load for principal solutions with different loading rates. Color red, green, and blue represent respectively different loading rate = 10 -3 , 10 -2 , and 10 -1 .

  modelingA model of J 2 plasticity associated with power-law isotropic hardening is chosen to match the mechanical properties of annealed Al 6061-O and will be used in all following calculations. The plastic response of Al 6061-O was investigated by[START_REF] Zhang | On the dynamics of necking and fragmentation -I. realtime and post-mortem observations in Al 6061[END_REF] through uniaxial tensile tests and is formulated here as:σ := σ/G = 2(1 + ν) y 1 + β p n (4.3.1)expressed in measured quantities: Young's modulus E = 70 GPa, Poisson's ratio ν = 0.25, yield strain y = 3.57 × 10 -4 , β = 14165, and hardening power n = 0.22.

  Figure 4.3: A zoom-in of two exemplary meshes for thin rings (η = 10 -4 ). Red lines represent the mid-plane of rings and hold constant distance r to the origin in the unperturbed configuration.Blues lines represent the inner/outer surfaces of the rings which always retain a constant distance h between themselves. On the right-hand side, an additional radial shifts ∆w (θ) in the form of a Gaussian function is added onto the mid-plane, producing a highly localized bump on the original perfect ring.

  Figure 4.4: Schematic graphs showing the two kinds of imperfections used in the study, along with their wave components' renormalized amplitudes plotted versus wavenumbers. In both figures, dashed lines indicate the reference circle of unity radius; red lines indicate the perturbed reposed configurations with exaggeration. For the spatially isolated imperfection, the tuple of parameters A = 10 -2 , κ = 10, and θ = π is used, and the red line is 10 times exaggerated. For the randomly distributed imperfections, the discussed set of parameters is used, and the red line is 100 times exaggerated.

  Figure 4.7, the visual representations of the time evolution of the ring's profile are shown for rings with isolated imperfection and for different loading rate . The influence of loading rate is clear shown in the figures.

Figure 4 . 6 :

 46 Figure 4.6: Influence of the time of perturbation introduction on the subsequent growth/decay of the perturbation. The vertical dashed lines delimitated the localization time t l , which for the parameter values used here is a) t l ≈ 37 b) t l ≈ 8.8. Imperfections used here are velocity perturbation with isolated extent.

  Figure 4.8c, no dominant wavelength can be spotted. For the irregular spacing pattern showing in these three figures, we shall adopt a rationalization based on the local interaction between localized imperfection present in the structure. In Figures 4.8aand 4.8b, the post-bifurcation patterns have similar shapes up to amplitude

Figure 4 . 7 :

 47 Figures 4.8a and 4.8b have their time of onset of localization t l strictly less than the lower bound

  n δ kp δ lq +

		kl ϕ p,q dX 2 ,	(2.A.19)
	where δ ij is the Kronecker delta. The characteristic field	ij ϕ(X 2 ) is the unit cell's response to the
	ij-th component of the unit macroscopic deformation and is a periodic function in X 2 , with period
	the unit cell thickness H. It is calculated by solving the following boundary value problem given in
	its variational form:

H 0 L mnpq δ kp δ lq + kl ϕ p,q δϕ m,n dX 2 = 0.

(2.A.20)

From (2.A.19)-(2.A.20

), and recalling the orthotropy of the principal solution, one obtains the following expressions for the nonzero components of the homogenized moduli tensor:

  A.6) which substituted in (3.A.2) give the following relation between Ṡαβ and its work conjugate Ėαβ :

  To further illustrate the connections between finite propagation speed of signals and delay of onset of instability, one can calculate t s , a lower bound of time taken by the axial wave to propagation through the entire

	ring:
	t s := lim k→∞	2π/c g (k) = 2π/ 2(1 + ν) ≈ 3.16π.	(4.3.6)

Here and subsequently in this chapter, Latin indexes range from 1 to

2, unless indicated differently. Einstein's summation convention is implied over repeated indexes. Repeated indexes in parentheses are not summed, unless indicated explicitly.

Here and subsequently in this chapter g , g etc. denote first and second derivatives of g with respect to its argument.

Repeated indexes in parentheses are not summed, unless explicitly indicated by a summation symbol.

Note: Here and subsequently in this chapter Greek indexes range from 1 to

while Latin indexes range from 1 to

2 Note: Here and subsequently superimposed dot denotes time differentiation ḟ := ∂f (X, t)/∂t 3 For hyperelastic material, a linearization of perturbation yields ∆N αβ = L αβγδ ∆F γδ . To avoid algebraic complications we further assume that the same holds true for a rate-independent material.

Note: for a uniaxial stress state 2 = -1/2 and e = ; Moreover y and σy = E y are the yield strain and stress respectively in a uniaxial loading path.

Acknowledgements

Firstly, I would like to express my sincere gratitude to my advisor Prof. N. Triantafyllidis for the continuous support and encouragement during the course of my thesis. Under his guidance, I had the chance not only to work on my thesis but also to learn some other interesting problems. I could not have imagined having a better advisor and mentor for my Ph.D. studies.

Acknowledgements vii

corrections, and suggestions that ensued. My thanks also go to all the members of my thesis committee for their support and encouragement. This thesis would not have been possible without the generous support from the Monge Doctoral Fellowship and the Laboratoire de Mécanique

2.B Post-bifurcation asymptotics for a compressible neo-Hookean composite

It is now of interest to see if the stability result for the neo-Hookean composite is influenced by compressibility. To this end we now consider the compressible neo-Hookean composite in (2.3.3) a which also admits a closed form solution. Recalling again from the kinematics in (2.2.8) that F 11 = 1λ and F 21 = γ one obtains from (2.2.9) the following expressions for the first Piola-Kirchhoff stresses:

(2.B.1)

Again recalling shear traction continuity Π f 21 = Π m 21 from (2.2.9) 3 and the kinematic constraint F 12 = 0 from (2.2.8) 3 one obtains the following linear system for the two unknowns F f 12 and F m 12 :

(2.B.2)

The solution of the above linear system for F f 12 , F m 12 gives:

where point is found by taking the second γ-derivative of (2.B.6) at γ = 0, giving after some straightforward by lengthy algebra:

which reduces, as expected, to its incompressible limit (2.3.14) as r f = r m = 0. Notice that λ 2 > 0,

showing that the compressible neo-Hookean composite has like its incompressible counterpart a stable post-bifurcation response under displacement control.

At leading order, the dynamics of flexural modes are governed by Mathieu-Hill equations whose characteristic curves of associated Mathieu functions delineate boundaries of instability domains within the control parameter plane of load versus ring's slenderness.

All the above-mentioned works were concerned with the stability of ring vibrations and not

with their stability at high loading rates as is the case of interest here. Our investigation is further motivated by work involving rings with high strain rate using electromagnetic loading-since this method avoids propagating waves-under tension that study the influence of high loading rate on metal ductility [START_REF] Gourdin | Analysis and assessment of electromagnetic ring expansion at a high-strain-rate test[END_REF][START_REF] Triantafyllidis | Onset of necking in electro-magnetically formed rings[END_REF]Zhang andRavi-Chandar, 2006, 2008), and in particular by experiments in ring and cylinder under electromagnetic compression by [START_REF] Anderson | Dynamic pulse buckling of cylindrical shells under transient lateral pressures[END_REF] and [START_REF] Jones | Dynamic plastic buckling of rings and cylindrical shells[END_REF], since these experiments combine structural instability with rapid loading. It is the most recent experimental work of [START_REF] Mainy | Dynamic buckling of thin metallic rings under external pressure[END_REF] that serves as the starting point for this investigation, and in particular, the localized failure patterns observed (see Fig. 4.1), which are in remarkable contrast to global buckling modes of externally pressurized rings under quasistatic loading rates. In order to keep essential features such as buckling under static loading and finite wave speeds for all wavenumbers, we concentrate on the dynamics of an elastic ring following a von Karman-Timoshenko theory allowing for small strains, moderate rotations, transverse shear and rotational inertia. The ring's stability is studied by following the evolution of a localized small perturbation. It is shown that for small values of the applied loading rate the structure fails through a global mode, while for large values of the applied loading rate the structure fails by a localized mode of deformation. Following Section 4.1 the presentation of the work continues with Section 4.2, where we derive the equations of motion, outline the numerical scheme for the principal solution of these equations, and present the linearized analysis of the initial growth/decay of a perturbation. The results are given in Section 4.3 where we conduct numerical calculations of the evolution of different types of spatially localized imperfection/perturbation and a discussion in Section 4.4 concludes this chapter.

Theory

In the first subsection, we derive the equations of motion from Hamilton's variational principle, from which we deduce the structure's Euler-Lagrange equations. In the second subsection, we study the behavior of the principle solution when dynamic loading condition with different loading rates are applied onto the ring. In the third subsection, we carry out linearized analysis on the problem to study the initial growth/decay of a perturbation introduced at a finite moment t = t 0 to the system. Suppose there are N nodes on the mid-plane indexed by i, then the t l is calculated as:

The relations between time and mid are plotted in Figure 4.5 for a product combination of (isolated imperfection, randomly distributed imperfections) and ∈ 10 -1 , 10 -2 , 10 -3 , 10 -4 , 10 -5 , 10 -6

for geometrically imperfect rings. All time increments after mid > 1 are cut off on purpose, as this kind of axial tension in a real compression experiment is quite unphysical. Numerically, this Résumé : Localization of deformation in finitely strained ductile solids is the instability mechanism leading to their failure by rupture. This phenomenon occurs under static and dynamic loading conditions. It can appear in the bulk of solids, in which case it is referred to as a material instability phenomenon or in a structure, in which case one talks about a structural instability problem. The thesis at hand studies localization in the material and structural context, both under static and dynamic conditions, using a common tool: the evolution of a geometrically localized perturbation.

An introduction to the localization of deformation problem in solid mechanics is presented in Chapter 1. The material instability aspect of localization of deformation in microstructured solids under quasistatic loading and its connection to macroscopic ellipticity-the continuum criterion for the presence of a discontinuous strain field-are addressed in Chapter 2. In this part we show the connection of the homogenized post-bifurcation response to the presence or absence of a localized deformation field in an infinite, fiber-reinforced composite under plane strain compression. The material instability aspect of localization of deformation under dynamic loading, i.e. where inertia becomes important, is addressed in Chapter 3. In this part we study the influence cones for the wave propagation emerging from a point perturbation in an infinite, biaxially strained plate whose constitutive response loses ellipticity at finite level of strain. The structural instability aspect of localization of deformation is investigated under dynamic loading conditions in Chapter 4 by studying the dynamic compression on an electromagnetically loaded metallic ring. In contrast to the quasistatic case, where a global failure mode is observed, the failure pattern for the rapidly compressed ring shows highly localized deformation areas. Le premier chapitre sert d'introduction au problème de localisation de déformation en mécanique des solides. Le deuxième chapitre quant à lui, porte sur l'instabilité matérielle de la localisation de déformation dans des milieux à microstructure sous compression quasi-statique, ainsi que son lien à l'ellipticité macroscopique (qui n'est autre que le critère continue de la présence d'un champ de déformation discontinue). Dans ce chapitre nous démontrons la relation entre la solution de postbifurcation homogénéisée et la présence ou l'absence d'un champ de déformation localisé dans un domaine de composites renforcées de fibres infiniment large sous compression. Le troisième chapitre est consacré à l'aspect de l'instabilité matérielle de la localisation de déformation sous chargement dynamique, où l'effet d'inertie devient non négligeable. Dans ce chapitre nous étudions une perturbation singulière sur une plaque infiniment large sous tension biaxiale et ses cônes d'influence, avec une loi de comportement qui perd l'ellipticité. Le quatrième chapitre étudie l'instabilité géométrique de la localisation de déformation sous chargement dynamique et pour ceci nous nous intéressons à la compression dynamique d'un anneau métallique sous chargement électromagnétique. Contrairement au cas quasi-statique, des domaines de déformation localisée sont observés dans le mode de défaillance de l'anneau.