Stability and Bifurcation in Nonlinear Mechanics

Abstract

Analysis of stability and bifurcation is studied in nonlinear mechanics with mechanisms of dissipation: plasticity, damage, fracture. With introduction of a set of internal variables, this framework allows a systematic description of the material behavior via two potentials: the free energy and the potential of dissipation. For standard generalized materials, internal state evolution is governed by a variational inequality depending on the mechanism of dissipation. This inequality is obtained through energetic considerations in an unified description based upon energy and driving forces associated with internal variable evolution. This formulation provides criterion for existence and uniqueness of the system evolution. Examples are presented for plasticity, fracture and damaged materials.

1. Introduction

Behavior of material based on an energetic approach is providing a large framework for the description of anelastic structures. Various approaches have been developed. The introduction of the internal variables allows a systematic description of the material behavior via two potentials: the free energy and the potential of dissipation.

The development of such description is due to the works of several authors [1, 2, 3, 4]. The purpose of this chapter is to study the quasistatic evolution of a anelastic structure. The system evolution is analyzed using the definition of functionals presented here in the case of nonlinear dynamics, firstly for internal variables associated with volume dissipation in nonlinear mechanics (plasticity and damage), secondly for dissipation due to singularities and discontinuity propagation (fracture, phase transformation). Quasistatic evolution is studied for dissipative materials. Stability and uniqueness of the response of the system under prescribed loading are discussed, due to the formulation of the rate boundary value problem in terms of velocity and evolution of internal parameters.

2. Preliminaries and general features

Let a body $\mathrm{\Omega }$ submitted to external forces described by vector fields $\underset{¯}{f}$ over $\mathrm{\Omega }$ and vector fields $\underset{¯}{T}$ along the boundary $\mathrm{\partial \Omega }$. The external forces are generally functions of time. Under loading the body is deformed. The actual position $\underset{¯}{x}$ of a material point is a function $\underset{¯}{\mathrm{\Phi }}$ of its initial position $\underset{¯}{X}$ and of the time. The displacement $\underset{¯}{u}$ is then defined by:

$$\underset{¯}{x}\left(\underset{¯}{X},t\right)=\underset{¯}{\mathrm{\Phi }}\left(\underset{¯}{X},t\right)=\underset{¯}{X}+\underset{¯}{u}\left(\underset{¯}{X},t\right)$$E1

Hence, a material element $\mathrm{d}\underset{¯}{X}$ is transported by the motion to the material element $\mathrm{d}\underset{¯}{x}$. The corresponding transformation is the linear application associated with the gradient of transformation $\mathbf{F}$:

$$\mathrm{d}\underset{¯}{x}=\frac{\partial \underset{¯}{x}}{\partial \underset{¯}{X}}.\mathrm{d}\underset{¯}{X}=\mathbf{F}.\mathrm{d}\underset{¯}{X}$$E2

The actual length of the material element is given by:

$$\mathrm{d}\underset{¯}{x}.\mathrm{d}\underset{¯}{x}=\mathrm{d}\underset{¯}{X}.{\mathbf{F}}^T.\mathbf{F}.\mathrm{d}\underset{¯}{X}=\mathrm{d}\underset{¯}{X}.\mathbf{C}.\mathrm{d}\underset{¯}{X}$$E3

The changes of the local geometry, the stretching, and the shearing of material fibers are determined by the Cauchy-Green tensor $\mathbf{C}={\mathbf{F}}^T.\mathbf{F}$. In small perturbations, the gradient of the displacement is small and the deformation is reduced to its linear contribution $\epsilon \left(\underset{¯}{u}\right)$: $2\epsilon \left(\underset{¯}{u}\right)=\nabla \underset{¯}{u}+{\nabla }^T\underset{¯}{u}$.

2.1 Notion of stability

For conservative system, when the loading $\underset{¯}{T}$ depends on one parameter $\lambda$, the dynamical system associated with the evolution of the body $\mathrm{\Omega }$ is defined by a functional

$$\underset{¯}{\dot{x}}=F\left(\underset{¯}{x},\lambda \right)$$E4

Then, positions of equilibrium are given by $F\left(\underset{¯}{x},\lambda \right)=0$. At this stage without any particular conditions, the uniqueness $\underset{¯}{x}\left(\lambda \right)$ is not ensured. But for a known position $\underset{¯}{x}\left(\lambda \right)$ under small perturbation $\mathit{d\lambda }$ it is possible to determine the corresponding variation $d\underset{¯}{x}$ of the position $\underset{¯}{x}\left(\lambda \right)$. Secondly, given some perturbation of equilibrium at fixed $\lambda$, if the response remains closed to that position, the equilibrium is say stable.

The stability of the position of equilibrium ${\underset{¯}{x}}_o\left(\lambda \right)$ is then determined with respect to Lyapounov definition:

$$\forall \epsilon \exists \alpha \left(\Vert \underset{¯}{x}\left(0,\lambda \right)-{\underset{¯}{x}}_o\left(\lambda \right)\Vert +\Vert \underset{¯}{\dot{x}}\left(0,\lambda \right)\Vert \right)\le \alpha \to \Vert \underset{¯}{x}\left(t,\lambda \right)-{\underset{¯}{x}}_o\Vert \le \epsilon$$E5

where $\underset{¯}{x}\left(t,\lambda \right)$ is solution of (Eq. 4) with initial conditions near the equilibrium state. It is clear that the notion of stability of an equilibrium position is a dynamical notion.

3. Study of conservative system

The evolution of the system is governed by the total potential energy, which is composed by the free energy of the material, and by the potential energy of the loading. A generic point of view is given by the study of the metronome (Figure 1) and introduction of asymptotic expansion to characterize the equilibrium solution [5, 6].

Figure 1.

The metronome.

A vertical rigid bar is clamped by a spiral spring with free energy $W\left(\theta \right)$, and we applied a vertical loading $\lambda$. A mass $M$ is attached at the top of the bar. The total potential energy $\mathcal{E}$ and the kinetic energy satisfy:

$$\mathcal{E}\left(\theta ,\lambda \right)=W\left(\theta \right)-\mathit{\lambda L}\left(1-cos\theta \right),\mathcal{K}=\frac{1}{2}M{\dot{\theta }}^2{L}^2$$E6

Near the position $\theta =0$ the energy is developed as

$$W=\frac{1}{2}{C}_1{\theta }^2+\frac{1}{3}{C}_2{\theta }^3+\frac{1}{4}{C}_3{\theta }^4+\dots$$E7

The dynamical system to study becomes

$${\mathit{ML}}^2\ddot{\theta }+\frac{\partial \mathcal{E}}{\partial \theta }=0$$E8

First, we characterize equilibrium position (${\theta }_o,{\lambda }_o$) and the research of equilibrium position near this point is determined by asymptotic expansion as proposed [5, 6], linking loading $\lambda$ to a position $\theta$.

3.1 Static equilibrium path

An equilibrium state ($\lambda ,\theta$) satisfies

$$\frac{\partial W}{\partial \theta }-\mathit{\lambda L}sin\theta =\theta \left(C_1+C_2\theta +C_3{\theta }^2+\dots -\mathit{\lambda L}\left(1-\frac{{\theta }^2}{6}..\right)\right)=0$$E9

then two equilibrium paths exists

$$\theta \left(\lambda \right)=0,\forall \lambda ,\mathit{\lambda L}=\frac{1}{sin\left(\theta \right)}\frac{\partial W}{\partial \theta }=C_1+\left(C_2+\frac{1}{6}{C}_1\right)\theta +\dots$$E10

The common point of the paths is the bifurcation point:

$${\lambda }_c=\frac{C_1}{L},\theta =0$$E11

3.2 Stability analysis

The dynamical behavior around this position is a weakly nonlinear dynamical system, taking account of a new asymptotic expansion

$$\begin{array}{l}\lambda ={\lambda }_o+{\lambda }_1\xi +{\lambda }_2{\xi }^2+\dots \\ \theta ={\theta }_o+{\theta }_1\left(\tau \right)\xi +{\theta }_2\left(\tau \right){\xi }^2+\dots \\ \tau ={\xi }^{m}t\left({\mathrm{\Omega }}_o+\xi {\mathrm{\Omega }}_1+{\xi }^2{\mathrm{\Omega }}_2+\dots \right)\end{array}$$E12

The characteristic time $\tau$ is chosen to satisfy the dependency of the pulsation of the system with respect to the loading.

$${\mathit{ML}}^2{\xi }^{2m}{\left(\mathrm{\Omega }\left(\xi \right)\right)}^2\ddot{\theta }+\frac{\partial W}{\partial \theta }-\mathit{\lambda L}sin\theta =0$$E13

The motion is then governed by

$$\begin{array}{l}{\mathit{ML}}^2{\xi }^{2m}\left({\mathrm{\Omega }}_o+\xi {\mathrm{\Omega }}_1+{\xi }^2{\mathrm{\Omega }}_2+\dots \right)\left({\ddot{\theta }}_1\xi +{\ddot{\theta }}_2{\xi }^2+\mathrm{..}\right)=\\ \left(L{\lambda }_o-C_1\right){\mathit{\xi \theta }}_1+{\xi }^2\left(L\left({\lambda }_o{\theta }_2-{\lambda }_1{\theta }_1\right)-\left(C_1{\theta }_2+C_2{\theta }_1^2\right)\right)\\ +\dots \end{array}$$E14

3.2.1 Discussion

• If ${\lambda }_o\ne {\lambda }_c$ then $m=0$ and we have

  $${\mathit{ML}}^2{\mathrm{\Omega }}_o^2{\ddot{\theta }}_1=\left({\lambda }_o-{\lambda }_c\right)L{\theta }_1$$E15

  then ${\lambda }_o\le {\lambda }_c$ the position $\theta =0$ is stable, ${\mathrm{\Omega }}_o^2=\frac{{\lambda }_c-{\lambda }_o}{\mathit{ML}}>0$; and unstable for $\lambda >{\lambda }_c$.

• If $\lambda ={\lambda }_c$ and ${\lambda }_1\ne 0$, then $m=\frac{1}{2}$ this implies that $\xi \ge 0$ and

  $${\mathit{ML}}^2{\mathrm{\Omega }}_o^2{\ddot{\theta }}_1=\left(L{\lambda }_1-C_2{\theta }_1\right){\theta }_1$$E16

  If $C_2\ne 0$, two positions of equilibrium exist: $\theta =0$ and ${\theta }_e={\lambda }_1\frac{L}{C_2}$. A position of the fundamental path $\left(\lambda ={\lambda }_c+{\lambda }_1\xi ,0\right)$ with ${\lambda }_1<0$ is stable, unstable otherwise.

The position ${\theta }_e$ is stable if ${\lambda }_1>0$. Finally, the position $\left({\lambda }_c,0\right)$ is unstable.

• $C_2=0$. It is necessary to consider ${\lambda }_1=0$, $m=1$, and the motion is governed by

  $${\mathit{ML}}^2{\mathrm{\Omega }}_o^2{\ddot{\theta }}_1=\left({\lambda }_{2}L-\left(C_3+\frac{L{\lambda }_c}{6}\right){\theta }_1^2\right){\theta }_1$$E17

Then, we have three positions of equilibrium, one along the fundamental path $\left(\lambda ={\lambda }_c+{\lambda }_2{\xi }^2,0\right)$, and two other

$$\lambda ={\lambda }_c+{\lambda }_2{\xi }^2,\theta =\pm \xi \sqrt{\frac{{\lambda }_{2}L}{C_3+\frac{{\lambda }_c}{6}}}$$E18

The fundamental path $\theta =0$ is stable if $\lambda <{\lambda }_c$ and stability of position along the secondary path if and only if ${\lambda }_2>0$, in this case the bifurcation point is a stable point of equilibrium.

The results are resumed on the following picture, with fundamental path ($\theta =0,\lambda$), and particular phase diagram (Figures 2–4).

Figure 2.

Case ${\lambda }_1<0\ne 0$, Stable (s) and unstable (u) paths. Phase diagram for $\lambda <{\lambda }_c$.

Figure 3.

Case ${\lambda }_1=0$, ${\lambda }_2>0$, Stable (s) and unstable (u) paths. Phase diagram for $\lambda >{\lambda }_c$.

Figure 4.

Case ${\lambda }_1=0$, ${\lambda }_2<0$, Stable (s) and unstable (u) paths. Phase diagram for $\lambda <{\lambda }_c$.

This description of conservative system is well known, the systematic proposed expansion can be used for study stability of beams, plates,…, the displacement $\theta$ is replaced by a vector displacement. The second derivative of the potential energy plays a fundamental rule, when this quadratic form is positive definite, then uniqueness is ensured, that is Lejeune-Dirichlet theorem. For non-conservative system, the proposed asymptotic expansion should be used, static-uniqueness does not ensure Lyapounov stability, as illustrated with a bi-pendulum under following load [7].

4. Mechanical behavior with time independent processes

Let us consider a local free energy $\psi \left(\epsilon ,\alpha \right)$ depending on internal parameters $\alpha$, the total potential energy becomes

$$\mathcal{E}\left(\underset{¯}{u},\tilde{\alpha },{\underset{¯}{T}}^d\right)={\int }_{\mathrm{\Omega }}\psi \left(\epsilon \left(\underset{¯}{u}\right),\alpha \right)\mathrm{d\Omega }-{\int }_{\partial {\mathrm{\Omega }}_T}{\underset{¯}{T}}^d.\underset{¯}{u}\mathrm{d}a.$$E19

The admissible fields $\underset{¯}{u}$ satisfy $\underset{¯}{u}={\underset{¯}{u}}^d$ along $\partial {\mathrm{\Omega }}_u$.

The evolution of internal parameters $\alpha$ is given by additional constitutive law.

Let us consider now time-independent processes, hence there is no viscosity. This framework permits description of dry friction, plasticity, damage, and fracture.

Let us consider the support function of a convex $\mathcal{C}$ defined by a regular convex function $f$ of the driving force $\mathrm{A}$

$$\mathrm{A}\in \mathcal{C}=\left\{\mathrm{B}/f\left(\mathrm{B}\right)\le 0\right\}$$E20

The evolution of internal variables verifies the normality rule:

$$\dot{\alpha }=\lambda \frac{\partial f}{\mathrm{\partial A}}=\lambda \mathcal{N},\lambda \ge 0,f\left(\mathrm{A}\right)\lambda =0$$E21

The internal parameter $\alpha$ evolves if the driving force $\mathrm{A}$ satisfies $f\left(\mathrm{A}\right)=0$, otherwise the internal parameter cannot evolve. The rate of $\alpha$ is normal to the equi-potential surface $f=0$, the notation $\mathcal{N}=\frac{\partial f}{\mathrm{\partial A}}$ is then adopted.

The loading history is described by a increasing parameter $\tau$. At each state $\tau$, the driving force $\mathrm{A}\left(x,\tau \right)$ satisfies the inequality $f\left(\mathrm{A}\left(x,\tau \right)\right)\le 0$ and the state equations $\mathrm{A}=-\frac{\partial \psi }{\partial \alpha }$.

The system is in equilibrium during the time, so that at the current state, the potential energy is stationary among the set of admissible displacements $\delta \underset{¯}{u}$ which satisfy $\delta \underset{¯}{u}=0$ over $\partial {\mathrm{\Omega }}_u$:

$$\frac{\partial \mathcal{E}}{\partial \underset{¯}{u}}.\delta \underset{¯}{u}=0$$E22

The variations of the potential energy are equivalent to

$$div\mathit{\sigma }=0,\mathit{\sigma }=\frac{\partial \psi }{\partial \epsilon },\mathit{\sigma }.\underset{¯}{n}={\underset{¯}{T}}^d\text{over}\partial {\mathrm{\Omega }}_T$$E23

These equations are true at each state of applied loading, so the evolution of equilibrium is given by ($\dot{f}$ is the derivative with respect to the fictitious time $\tau$):

$$\mathit{div}\dot{\sigma }=0\text{and}\dot{\sigma }=\frac{{\partial }^2\psi }{\partial \epsilon \partial \epsilon }:\dot{\epsilon }+\frac{{\partial }^2\psi }{\partial \epsilon \partial \alpha }:\dot{\alpha }\text{in}\mathrm{\Omega }\text{and}\dot{\sigma }.\underset{¯}{n}={\underset{¯}{T}}^d\text{over}\partial {\mathrm{\Omega }}_T$$E24

which is equivalent to

$$0=\frac{d}{\mathit{d\tau }}\left(\frac{\partial \mathcal{E}}{\partial \underset{¯}{u}}\right).\delta \underset{¯}{u}={\int }_{\mathrm{\Omega }}\epsilon \left(\delta \underset{¯}{u}\right):\left(\frac{{\partial }^2\psi }{\partial \epsilon \partial \epsilon }:\dot{\epsilon }+\frac{{\partial }^2\psi }{\partial \epsilon \partial \alpha }:\dot{\alpha }\right)\mathrm{d\Omega }-{\int }_{\partial {\mathrm{\Omega }}_T}{\underset{¯}{T}}^d.\delta \underset{¯}{u}\mathrm{d}a$$E25

The current state is determined by the evolution of the internal state $\dot{\alpha }$. To determine existence and uniqueness of the evolution of the system, the rate boundary value problem must be studied as pointed out in Refs. [8, 9].

4.1 Evolution of $\alpha$

Considering the normality rule, we can conclude that

$$\lambda \ge 0,f\le 0,\mathit{\lambda f}=0$$E26

For an internal state such that $f=0$, the evolution of $f$ satisfies $\dot{f}\le 0$, and simultaneously the time derivative of the condition $\mathit{\lambda f}=0$ implies that

$$\stackrel{.}{\overline{\mathit{\lambda f}}}=\dot{\lambda }f+\lambda \dot{f}=0$$E27

When $f=0$, $\lambda \dot{f}=0$ then $\lambda >0$ if and only if $\dot{f}=0$, that is the classical consistency condition. This provides the definition of the set $\mathcal{P}$ of the admissible fields $\lambda \left(x\right)$.

At each state $\tau$, the domain $\mathrm{\Omega }$ is decomposed into two complementary sub-domains ${\mathrm{\Omega }}_r$ and $I^{\tau }$ such that:

$$x\in {\mathrm{\Omega }}_r=\left\{x\in \mathrm{\Omega }|f\left(\mathrm{A}\left(x,\tau \right)\right)<0\right\},x\in I^{\tau }=\left\{x\in \mathrm{\Omega }|f\left(\mathrm{A}\left(x,\tau \right)\right)=0\right\}$$E28

Then $\mathcal{P}$ is defined as:

$$\mathcal{P}=\left\{\beta \left(x\right)|\beta \left(x\right)=0,\forall x\in {\mathrm{\Omega }}_r\text{and}\beta \left(x\right)\ge 0,\forall x\in I^{\tau }\right\}.$$E29

It is obvious that the field $\lambda \left(x\right)$ is an element of $\mathcal{P}$.

Considering now a point $x\in I^{\tau }$ then $\lambda \left(x\right)\ge 0$ and consequently $\dot{f}\lambda =0$. As $\dot{f}\le 0$, we deduce:

$$\lambda \left(x\right)\in \mathcal{P},\left(\lambda \left(x\right)-\beta \left(x\right)\right)\dot{f}\ge 0,\forall \beta \left(x\right)\in \mathcal{P}$$E30

and

$${\int }_{\mathrm{\Omega }}\left(\lambda \left(x\right)-\beta \left(x\right)\right)\dot{f}\mathrm{d\Omega }\ge 0$$E31

among the set $\mathcal{P}$ of admissible fields $\beta$. This is a variational inequality.

By using now the definition of $f$, considering the equations of state for $\mathrm{A}$ and the normality rule for $\dot{\alpha }$, the inequality (Eq. (31)) is rewritten as ($\mathcal{N}=\frac{\partial f}{\mathrm{\partial A}}$):

$${\int }_{\mathrm{\Omega }}\left(\lambda \left(x\right)-\beta \left(x\right)\right)\left(\mathcal{N}:\frac{{\partial }^2\psi }{\partial \alpha \partial \epsilon }:\dot{\epsilon }+\mathcal{N}:\frac{{\partial }^2\psi }{\partial \alpha \partial \epsilon }:\mathcal{N}\lambda \right)\mathrm{d\Omega }\le 0.$$E32

This inequality is a formulation similar to that of Ref. [8].

4.2 The rate boundary value problem

Let us consider the functional $\mathcal{F}$ based on the velocities:

$$\mathcal{F}\left(\underset{¯}{\tilde{v}},\tilde{\lambda },{\underset{¯}{\tilde{\dot{T}}}}^d\right)={\int }_{\mathrm{\Omega }}\frac{1}{2}\epsilon \left(\underset{¯}{v}\right):\mathbb{C}:\epsilon \left(\underset{¯}{v}\right)+\epsilon \left(\underset{¯}{v}\right):M\lambda +\frac{1}{2}\mathit{\lambda H\lambda }\mathrm{d\Omega }-{\int }_{\partial {\mathrm{\Omega }}_T}{\underset{¯}{\dot{T}}}^d.\underset{¯}{v}\mathrm{d}a$$E33

with the notations: $\mathbb{C}=\frac{{\partial }^2\psi }{\partial \epsilon \partial \epsilon },M=\frac{{\partial }^2\psi }{\partial \epsilon \partial \alpha }:\mathcal{N},H=\mathcal{N}:\frac{{\partial }^2\psi }{\partial \alpha \partial \alpha }:\mathcal{N}$.

The solution of the rate boundary value problem satisfies the variational inequality

$$\frac{\partial \mathcal{F}}{\partial \underset{¯}{\tilde{v}}}.\left(\underset{¯}{\tilde{v}}-{\underset{¯}{\tilde{v}}}^{\ast }\right)+\frac{\partial \mathcal{F}}{\partial \tilde{\lambda }}.\left(\tilde{\lambda }-{\tilde{\lambda }}^{\ast }\right)\le 0$$E34

among the set of admissible fields $\left({\underset{¯}{\tilde{v}}}^{\ast },{\tilde{\lambda }}^{\ast }\right)\in \mathcal{K}x\mathcal{P}$ with $\mathcal{K}=\left\{\underset{¯}{v},|,\underset{¯}{v},=,{\underset{¯}{v}}^d,\text{over},\partial ,{\mathrm{\Omega }}_u\right\}$.

Generally, the modulus of elasticity $\mathbb{C}$ is a quadratic positive-definite operator, then the field $\underset{¯}{v}$ is unique for given $\dot{\alpha }$. So the velocity $\underset{¯}{v}$ can be eliminated: ${\underset{¯}{v}}^{\mathit{sol}}=\underset{¯}{v}\left(\tilde{\lambda },{\underset{¯}{\tilde{\dot{T}}}}^d,{\underset{¯}{\tilde{v}}}^d\right)$ is linear of each argument and so a new functional, that is defined only on the internal variables, which is quadratic in $\tilde{\lambda }$ is defined:

$${\mathcal{F}}^{\ast }\left(\tilde{\lambda }\right)=\mathcal{F}\left({\underset{¯}{\tilde{v}}}^{\mathit{sol}}\left(\tilde{\lambda },{\underset{¯}{\tilde{\dot{T}}}}^d,{\underset{¯}{v}}^d\right),\tilde{\lambda },{\underset{¯}{\tilde{\dot{T}}}}^d\right)=\frac{1}{2}\tilde{\lambda }.Q.\tilde{\lambda }-\tilde{\lambda }.T\left({\underset{¯}{v}}^d,{\underset{¯}{\dot{T}}}^d\right)$$E35

Stability condition. It is known that a solution exist if

$$\forall \tilde{\beta }\in \mathcal{P},\tilde{\beta }.Q.\tilde{\beta }\ge 0,$$E36

where $\mathcal{P}$ is the set of admissible fields (Eq. 29). This condition of existence ensures that the current state is stable.

Uniqueness and no-bifurcation. The solution of the boundary value problem is also unique if

$$\forall \tilde{\beta }\in {\mathcal{P}}^{\ast },\tilde{\beta }.Q.\tilde{\beta }\ge 0,$$E37

where ${\mathcal{P}}^{\ast }$ is the set

$${\mathcal{P}}^{\ast }=\left\{\tilde{\beta }|\beta \left(x\right)=0,\forall x\in {\mathrm{\Omega }}_r,\beta \left(x\right)\ne 0,x\in I^{\tau }\right\}$$E38

This condition ensures that there is no bifurcation.

4.3 Property of the functional

Let us consider that a solution is determined, the domain $I^{\tau }$ is decomposed in three different domains depending on $\lambda >0$ or not:

• the loading zone $I_{+}^{\tau }=\left\{x\in \mathrm{\Omega }/x\in I^{\tau },{\mu }^{\tau }\left(x\right)>0,{\dot{f}}^{\tau }\left(x\right)=0\right\}$

• the unloading zone $I_{-}^{\tau }=\left\{x\in \mathrm{\Omega }/x\in I^{\tau },{\mu }^{\tau }\left(x\right)=0,{\dot{f}}^{\tau }\left(x\right)<0\right\}$

• the neutral zone $I_o^{\tau }=\left\{x\in \mathrm{\Omega }/x\in I^{\tau },{\mu }^{\tau }\left(x\right)=0,{\dot{f}}^{\tau }\left(x\right)=0\right\}$

Introducing asymptotic expansion to define a loading path with parameter $\tau$

$$\begin{array}{l}{\underset{¯}{T}}^d={\underset{¯}{T}}_0^d+\tau {\underset{¯}{T}}_1^d+{\tau }^2{\underset{¯}{T}}_2^d+\dots \\ {\underset{¯}{u}}^d={\underset{¯}{u}}_0^d+\tau {\underset{¯}{u}}_1^d+{\tau }^2{\underset{¯}{u}}_2^d+\dots \end{array}$$E39

A local response in terms of displacement and internal variable fields is assumed to be developed also as an asymptotic expansion:

$$\alpha ={\alpha }_0+{\mathit{\tau \alpha }}_1+{\tau }^2{\alpha }_2+\dots \underset{¯}{u}={\underset{¯}{u}}_0+\tau {\underset{¯}{u}}_1+{\tau }^2{\underset{¯}{u}}_2+\dots$$E40

The term of order one corresponds to the solution of the boundary value problem in velocities. Similar asymptotic expansions are deduced from the yielding function $f$ satisfying the normality rule at each order and constraints are then obtained on the successive orders of the internal state. Hence, the characterization of order two shows that

$${\lambda }_2{f}_1+{\lambda }_1{f}_2=0$$E41

The properties of ${\lambda }_2$ are given related to the decomposition of $I^{\tau }$ and the field ${\lambda }_2$ is an element of the set ${\mathcal{P}}_2$

$${\mathcal{P}}_2=\left\{\tilde{\mu }/\left\{\begin{array}{ll}\mu \left(x\right)& =0,\text{if}x\cancel{\in }{I}^{\tau }\cup I_{+}^{\tau }\cup I_o^{\tau },\\ \mu \left(x\right)& \ge 0,\text{if}x\in I_o,\\ \mu \left(x\right)& \in \Re ,\text{if}x\in I_{+}^{\tau }\end{array}\right.\right\}$$E42

The boundary value problem for the order two has the same form that for order one, except that the linear term contains terms due to order one [10]

$${\mathcal{F}}_2\left({\underset{¯}{\tilde{u}}}_2,{\tilde{\lambda }}_2\right)=\frac{1}{2}Q\left(\left({\underset{¯}{\tilde{u}}}_2,{\tilde{\lambda }}_2\right),\left({\underset{¯}{\tilde{u}}}_2,{\tilde{\lambda }}_2\right)\right)+F_2\left(\left({\underset{¯}{\tilde{u}}}_1,{\tilde{\lambda }}_1\right),\left({\underset{¯}{\tilde{u}}}_2,{\tilde{\lambda }}_2\right)\right)$$E43

and the solution of the rate boundary value problem of order $2$ (satisfies)

$$\frac{\partial {\mathcal{F}}_2}{\partial {\underset{¯}{\tilde{u}}}_2}\left({\underset{¯}{\tilde{u}}}_2^{\ast }-{\underset{¯}{\tilde{u}}}_2\right)+\frac{\partial {\mathcal{F}}_2}{\partial {\tilde{\lambda }}_2}\left({\tilde{\lambda }}_2^{\ast }-{\tilde{\lambda }}_2\right)\ge 0$$E44

among the set of admissible field ${\underset{¯}{\tilde{u}}}_2^{\ast }$ which may satisfied the boundary conditions at order two and ${\tilde{\lambda }}_2^{\ast }$ is an element of ${\mathcal{P}}_2$.

The condition of stability on order two is quite different than the condition of order one, due to the presence of unloading zone. The condition of no-bifurcation is also changed taking account of $\lambda =0$ on $I_o^{\tau }$. The loss of positivity of $Q$ on these new spaces changes the critical value ${\underset{¯}{T}}_c^d,{\underset{¯}{u}}_c^d$.

5. The Shanley column

This model has been used by many authors [6], especially to study plastic buckling as discussed in Ref. [11]. The rigid rod model has two degrees of freedom: the downward vertical displacement $u$ and the rotation $\theta$ (Figure 5). The column is supported by the a uniformly distributed springs along the segment $\left[-l,l\right]$. The behavior of the spring is elasto-plastic with linear hardening.

Figure 5.

The Shanley column.

$$\psi \left(\epsilon ,\alpha \right)=\frac{1}{2}E{\left(\epsilon -\alpha \right)}^2+\frac{1}{2}H{\alpha }^2$$E45

Let us consider a state for which the plastic domain $I_{+}^{\tau }$ is $\left[d,l\right]$. The value of $d$ is determined by the condition of neutral loading $\left[\dot{\alpha }\right]\left(d,t\right)=0$. The equations of equilibrium are deduced from the potential energy

$$\mathcal{E}\left(u\left(x\right),\alpha \left(x\right),T^d\right)={\int }_{-l}^l\psi \left(\epsilon ,\alpha \right)\mathit{dx}+T^d\left(u+L\left(1-{\theta }^2/2\right)\right),\epsilon \left(x\right)=u-\mathit{x\theta }$$E46

then the equilibrium state obeys to

$$0=T^d+{\int }_{-l}^{+l}E\left(x\right)\left(\epsilon -\alpha \right)\mathit{dx},0=-T^d\mathit{L\theta }+{\int }_{-l}^{+l}E\left(x\right)x\left(\epsilon -\alpha \right)\mathit{dx}$$E47

These equations are valid during the loading process, taking account of the determination of $d\left(t\right)$. Then, we obtain

$$0={\dot{T}}^d+{\int }_{-l}^{d}E\dot{\epsilon }\mathit{dx}+{\int }_d^l{E}_T\dot{\epsilon }\mathit{dx},0=-L\stackrel{\cdot }{\overline{T^d\theta }}+{\int }_{-l}^{+d}\mathit{Ex}\dot{\epsilon }\mathit{dx}+{\int }_d^l{E}_{T}x\dot{\epsilon }\mathit{dx}$$E48

A non-trivial solution in $\theta$ is obtained by introducing the time-scale $\tau$ such that the velocity $x_1=\dot{d}$ of propagation of the unloading domain is finite. The domain $I_{+}^{\tau }=\left[x_{\tau },l\right]$ is defined with a asymptotic expansion

$$d=x_{\tau }=\sum _i{x}_i{\tau }^i$$E49

At point $x_{\tau }$, the condition $\left[\dot{\alpha }\left(x_{\tau },\tau \right)\right]=0$ where $\alpha \left(x,t\right)={\sum }_i{\alpha }_i\left(x\right){\tau }_i$ gives conditions on the asymptotic expansion:

$$0=\left[{\alpha }_1\left(x_o\right)\right]0=\left[{\alpha }_2\left(x_o\right)\right]+x_1\left[{\alpha }_1^{\prime }\left(x_o\right)\right]0=\left[{\alpha }_3\left(x_o\right)\right]+2x_1\left[{\alpha }_2^{\prime }\left(x_o\right)\right]+\left(2x_1^2+x_2\right)\left[{\alpha }_1^{\prime }\left(x_o\right)\right]$$E50

A non-trivial solution is then obtained as

$${\alpha }^{\ast }\left(x,\tau \right)=\alpha \left(\tau \right)+m\left(\tau \right)x$$E51

We can take the time derivative of the equilibrium equations taking account of the position of $x_{\tau }$ and of discontinuities (Eq. (50)) of the mechanical quantities on this boundary. It is obvious that we have:

$$\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{-l}^{l}f\left(x,\tau \right)\mathit{dx}={\int }_{-l}^l\dot{f}\left(x,\tau \right)\mathit{dx}+{\left[f\left(x,\tau \right)\right]}_{x_{\tau }^{-}}^{x_{\tau }^{+}}{\dot{x}}_{\tau }$$E52

We find $T_c=T_T=\frac{2l^{3}E}{3L}=\frac{H}{E+H}{T}_E$, $m_1=0$ and $m_2=-T_2/2{\mathit{Hl}}^2$,

$$\begin{array}{l}{x}_1^2=\frac{4l^2}{3}\frac{T^2}{T_E-T_T},T=T_T+T_2\frac{{\tau }^2}{2},\theta =\frac{{\mathit{ElT}}_2}{3\mathit{LH}\left(T_E-T_T\right)}\frac{{\tau }^2}{2}+\dots \end{array}$$E53

This is a bifurcated path. The condition of stability of the fundamental path $\left(\theta =0,T\right)$ is preserved for loading near $T=T_c$ but for $T\ge T_c$ another path exists which is also a stable path. This is quite different of conservative system, for which the bifurcation point corresponds to a loss of stability of the fundamental path.

More applications can be found in many papers for elasto-plasticity [12] with implications on the constitutive laws [9, 13]. Influence of pre-bifurcation conditions have been also analyzed [14, 15].

5.1 A simple model of fracture

Let us consider a straight beam under bending with fixed extremities $l_1,l_2$, where the beam is clamped. The length of the beam is $l_1+l_2$, the vertical displacement is $v\left(x\right)$ defined on segment $\left[-l_1,l_2\right]$ and the strain is given by: $v$: $\epsilon =v^{″}\left(x\right)y$. We apply a load at the origin, and we study the possibility of decohesion at points $l_1,l_2$. We study two cases, first the applied load is a vertical displacement $v\left(o\right)=V$ and second the load is controlled at the origin $T\left(o\right)=F$. For the first case, the potential energy at the equilibrium is

$$W\left(l_1,l_2,V\right)=\frac{3}{2}{\mathit{EIV}}^2\frac{{\left(l_1+l_2\right)}^3}{l_1^3{l}_2^3}=\frac{1}{2}{\mathit{kV}}^2$$E54

We define $J_i=-\frac{\partial W}{\partial l_i}$

$$J_1=3\frac{{\left(l_1+l_2\right)}^2}{l_1^3{l}_2^3}\frac{l_2}{l_1}$$E55

and $J_2$ is obtained permuting indices. The evolution of the delamination is given by the normality law

$${\dot{l}}_i\ge 0,J_i\le G_c,\left(J_i-G_c\right){\dot{l}}_i=0$$E56

and existence and uniqueness are given with respect to positivity or not of $Q$ such that

$$Q_{\mathit{ij}}=-J_{\mathit{ij}}=\frac{{\partial }^{2}W}{\partial l_i\partial l_j}=18{\mathit{EIV}}^2\frac{l_1+l_2}{l_1^5{l}_2^5}\left[\begin{array}{cc}{l}_2^3\left(l_1+2l_2\right)& l_1^2{l}_2^2\\ l_1^2{l}_2^2& l_1^3\left(2l_1+l_2\right)\end{array}\right]$$E57

For $l_1=l_2$, $Q$ is always positive definite, the position is then always stable and we have no bifurcation.

When the force is controlled

$$W\left(l_1,l_2,F\right)=-\frac{F^2}{k},k=3\mathit{EI}\frac{{\left(l_1+l_2\right)}^3}{l_1^3{l}_2^3}$$E58

The associated $Q$ matrix becomes

$$Q=-\frac{2F^2}{\mathit{EI}}\frac{l_1{l}_2}{{\left(l_1+l_2\right)}^2}\left[\begin{array}{cc}{l}_2^3\left(l_2-l_1\right)& 2l_1^2{l}_2^2\\ 2l_1^2{l}_2^2& l_1^3\left(l_1-l_2\right)\end{array}\right]$$E59

is always negative definite. The symmetric equilibrium is always an unstable state with possible bifurcation, the eigenvalue of $Q$ having opposite signs.

For multi-cracking of a body, the rate boundary value problem has been formulated and condition of existence and uniqueness have been deduced [16, 17].

6. Stability of moving surfaces

We study now a moving surface associated with a change of mechanical properties. This framework is used to describe damage or phase transformation. Variational formulations were performed to describe the evolution of the surface between the sound and the damaged material [18, 19, 20]. Connection with the notion of configurational forces [21] can be investigated.

6.1 Some general features

The domain $\mathrm{\Omega }$ is composed of two distinct volumes ${\mathrm{\Omega }}_1,{\mathrm{\Omega }}_2$ of materials with different mechanical characteristics. The bounding between the materials is perfect and the interface is denoted by $\mathrm{\Gamma }$, ($\mathrm{\Gamma }=\partial {\mathrm{\Omega }}_1\cap \partial {\mathrm{\Omega }}_2$). The external surface $\mathrm{\partial \Omega }$ is decomposed in two parts $\partial {\mathrm{\Omega }}_u$ and $\partial {\mathrm{\Omega }}_T$ on which the displacement ${\underset{¯}{u}}^d$ and the loading ${\underset{¯}{T}}^d$ are prescribed, respectively. We consider isotherm processes. The material $1$ changes into material 2 as the motion of the interface $\mathrm{\Gamma }$ by an irreversible process. Hence, $\mathrm{\Gamma }$ moves with the normal velocity $\underset{¯}{c}=\varphi \underset{¯}{\nu }$ in the reference state, $\underset{¯}{\nu }$ is the outward ${\mathrm{\Omega }}_2$ normal, then $\varphi$ is positive.

Along $\mathrm{\Gamma }$, the mechanical quantities $f$ can have a jump denoted by ${\left[f\right]}_{\mathrm{\Gamma }}=f_1-f_2$, and any volume average has a rate defined by

$$\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{\mathrm{\Omega }\left(\mathrm{\Gamma }\right)}f\mathrm{d\Omega }={\int }_{\mathrm{\Omega }\left(\mathrm{\Gamma }\right)}\dot{f}\mathrm{d\Omega }-{\int }_{\mathrm{\Gamma }}{\left[f\left({\underset{¯}{x}}_{\mathrm{\Gamma }},t\right)\right]}_{\mathrm{\Gamma }}\varphi \mathrm{d}a$$E60

where $\varphi$ is the normal propagation of the interface.

The state of the system is characterized by the displacement field $\underset{¯}{u}$, from which the strain field $\epsilon$ is derived. The main internal parameter is the spatial distribution of the two phases given by the position of the interface boundary $\mathrm{\Gamma }$. We analyze quasi-static motion of $\mathrm{\Gamma }$ under given loading prescribed on the boundary $\mathrm{\partial \Omega }$.

Introducing the total potential energy of the system

$$\mathcal{E}\left(\underset{¯}{u},\mathrm{\Gamma },{\underset{¯}{T}}^d\right)={\int }_{{\mathrm{\Omega }}_1}{\psi }_1\left(\epsilon \right)\mathrm{d\Omega }+{\int }_{{\mathrm{\Omega }}_2}{\psi }_2\left(\epsilon \right)\mathrm{d\Omega }-{\int }_{\mathrm{\partial \Omega }}{\underset{¯}{T}}^d.\underset{¯}{u}\mathrm{d}a$$E61

The behavior of the phase $i$ is assumed linear elastic. The state equations are reduced to

$${\psi }_i=\frac{1}{2}\epsilon :{\mathbb{C}}_i:\epsilon ,\mathit{\sigma }=\frac{\partial {\psi }_i}{\partial \epsilon }$$E62

We can notice that the position of the interface $\mathrm{\Gamma }$ becomes an internal parameter for the global system. The characterization of any equilibrium state is given by the stationary point of the potential energy ($\frac{\partial \mathcal{E}}{\partial \underset{¯}{u}}\cdot \delta \underset{¯}{u}=0$) among the set of the admissible field $\delta \underset{¯}{u}$ satisfying $\delta \underset{¯}{u}=0$ over $\partial {\mathrm{\Omega }}_u$. This formulation is equivalent to the set of local equations:

• local constitutive relations: $\mathit{\sigma }=\rho \frac{\partial {\psi }_i}{\partial \epsilon }={\mathbb{C}}_i:\epsilon ,\mathrm{on}{\mathrm{\Omega }}_i,$

• momentum equations: $div\mathit{\sigma }=0\mathrm{on}\mathrm{\Omega },{\left[\sigma \right]}_{\mathrm{\Gamma }}.\underset{¯}{\nu }=0\text{over}\mathrm{\Gamma },\mathit{\sigma }.\underset{¯}{n}={\underset{¯}{T}}^d\text{over}\partial {\mathrm{\Omega }}_T,$

• compatibility relations: $2\epsilon =\nabla \underset{¯}{u}+{\nabla }^t\underset{¯}{u},{\left[\underset{¯}{u}\right]}_{\mathrm{\Gamma }}=0\text{over}\mathrm{\Gamma },\underset{¯}{u}={\underset{¯}{u}}^d\text{over}\partial {\mathrm{\Omega }}_u.$

This equation emphasized the fact that the position of the interface $\mathrm{\Gamma }$ plays the role of internal parameters (Figure 6).

Figure 6.

Propagation of the interface.

The driving force associated with the motion of the interface $\mathrm{\Gamma }$ is obtained as

$$-\frac{\partial \mathcal{E}}{\mathrm{\partial \Gamma }}.\delta \mathrm{\Gamma }={\int }_{\mathrm{\Gamma }}G\left({\underset{¯}{x}}_{\mathrm{\Gamma }},t\right)\mathit{\delta \varphi }\left(s\right)\mathrm{d}a,G\left(s\right)={\left[\psi \right]}_{\mathrm{\Gamma }}-\underset{¯}{n}.\mathit{\sigma }.{\left[\nabla \underset{¯}{u}\right]}_{\mathrm{\Gamma }}.\underset{¯}{n}={\left[\psi \right]}_{\mathrm{\Gamma }}-\mathit{\sigma }:{\left[\epsilon \right]}_{\mathrm{\Gamma }}$$E63

An energy criterion is chosen as a generalized form of the well-known theory of Griffith. Then, we assume

$$\varphi \ge 0,\mathcal{G}\left({\underset{¯}{x}}_{\mathrm{\Gamma }},,,t\right)-G_c\le 0,\left(\mathcal{G}\left({\underset{¯}{x}}_{\mathrm{\Gamma }},t\right)-G_c\right)\varphi =0,$$E64

This decomposed the interface into two part ${\mathrm{\Gamma }}^{+}$ where $\mathcal{G}=G_c$ and the complementary part. At a point ${\underset{¯}{x}}_{\mathrm{\Gamma }}$ in ${\mathrm{\Gamma }}^{+}$, where the propagation occurs

$$\frac{\mathrm{d}}{\mathrm{d}t}\left(\mathcal{G}\left({\underset{¯}{x}}_{\mathrm{\Gamma }}\left(t\right),t\right)-G_c\right)=0,\text{and}\varphi \ge 0$$E65

The critical value is conserved following the moving interface: ${\mathrm{D}}_{\varphi }\mathcal{G}=0$. This leads to the consistency solution, which determines $\varphi \left({\underset{¯}{x}}_{\mathrm{\Gamma }}\right)$

$$\left(\varphi -{\varphi }^{\ast }\right){\mathrm{D}}_{\varphi }\mathcal{G}\left({\underset{¯}{x}}_{\mathrm{\Gamma }}\left(t\right),t\right)\ge 0,\forall {\varphi }^{\ast }\ge 0,\text{over}{\mathrm{\Gamma }}^{+}$$E66

For a given loading ${\underset{¯}{v}}^d,{\underset{¯}{T}}^d$ and a propagation $\varphi \left(s\right)$ of the interface, the evolution of the internal state satisfies

• local constitutive relations: $\dot{\mathit{\sigma }}={\mathbb{C}}_i:\epsilon \left(\underset{¯}{v}\right),\mathrm{on}{\mathrm{\Omega }}_i,$

• momentum equations: $div\dot{\mathit{\sigma }}=0\mathrm{on}\mathrm{\Omega },{\mathrm{D}}_{\varphi }\left({\left[\mathit{\sigma }\right]}_{\mathrm{\Gamma }}.\underset{¯}{\nu }\right)=0\text{over}\mathrm{\Gamma },\dot{\mathit{\sigma }}.\underset{¯}{n}={\underset{¯}{T}}^d\text{over}\partial {\mathrm{\Omega }}_T,$

• compatibility relations: $2\epsilon \left(\underset{¯}{v}\right)=\nabla \underset{¯}{v}+{\nabla }^t\underset{¯}{v},{\mathrm{D}}_{\varphi }\left({\left[\underset{¯}{u}\right]}_{\mathrm{\Gamma }}\right)=0\text{over}\mathrm{\Gamma },\underset{¯}{v}={\underset{¯}{v}}^d\text{over}\partial {\mathrm{\Omega }}_u.$

where

$$\begin{array}{c}{\mathrm{D}}_{\varphi }{\left[\underset{¯}{u}\right]}_{\mathrm{\Gamma }}={\left[\underset{¯}{v}\right]}_{\mathrm{\Gamma }}+\varphi {\left[\nabla \underset{¯}{u}\right]}_{\mathrm{\Gamma }}.\underset{¯}{\nu }=0,\\ {\mathrm{D}}_{\varphi }{\left[\mathit{\sigma }\right]}_{\mathrm{\Gamma }}.\underset{¯}{\nu }={\left[\dot{\mathit{\sigma }}\right]}_{\mathrm{\Gamma }}.\underset{¯}{\nu }-{\mathit{div}}_{\mathrm{\Gamma }}\left({\left[\mathit{\sigma }\right]}_{\mathrm{\Gamma }}\varphi \right)=0\end{array}$$E67

with ${div}_{\mathrm{\Gamma }}F=divF-\underset{¯}{\nu }.\nabla F.\underset{¯}{\nu }$. The velocity $\underset{¯}{v}$ is the solution of a problem of heterogeneous elasticity with boundary conditions linear with respect to the propagation $\varphi$: ${\underset{¯}{v}}^{\mathit{sol}}=\underset{¯}{v}\left(\varphi ,{\underset{¯}{v}}^d,{\underset{¯}{T}}^d\right)$. And we obtain

$$\begin{array}{c}{\mathrm{D}}_{\varphi }\mathcal{G}={\left[\mathit{\sigma }\right]}_{\mathrm{\Gamma }}:\nabla {\underset{¯}{v}}_1-{\dot{\mathit{\sigma }}}_2:{\left[\nabla \underset{¯}{u}\right]}_{\mathrm{\Gamma }}-\varphi G_n\\ G_n=-{\left[\mathit{\sigma }\right]}_{\mathrm{\Gamma }}:\left(\nabla \nabla {\underset{¯}{u}}_1.\underset{¯}{\nu }\right)+\nabla {\mathit{\sigma }}_2.\underset{¯}{\nu }:{\left[\underset{¯}{u}\right]}_{\mathrm{\Gamma }}\end{array}$$E68

Finally, the evolution of the system is determined by the functional

$$\mathcal{F}\left(\underset{¯}{v},\varphi \right)={\int }_{\mathrm{\Omega }}\frac{1}{2}\epsilon \left(\underset{¯}{v}\right):\mathbb{C}:\epsilon \left(\underset{¯}{v}\right)\mathrm{d\Omega }-{\int }_{\partial {\mathrm{\Omega }}_T}{\underset{¯}{\dot{T}}}^d.\underset{¯}{v}\mathrm{d}a-{\int }_{\mathrm{\Gamma }}\varphi {\left[\mathit{\sigma }\right]}_{\mathrm{\Gamma }}:\nabla {\underset{¯}{v}}_1\mathrm{d}a+\frac{1}{2}{\varphi }^2{G}_n\mathrm{d}a$$E69

and the variational inequality

$$\frac{\partial \mathcal{F}}{\partial \underset{¯}{v}}\left(\underset{¯}{v}-{\underset{¯}{v}}^{\ast }\right)+\frac{\partial \mathcal{F}}{\partial \varphi }\left(\varphi -{\varphi }^{\ast }\right)\ge 0$$E70

The stability of the actual state is determined by the condition of the existence of a solution $\left(\text{with},W,=,\mathcal{F},\left({\underset{¯}{v}}^{\mathit{sol}},\varphi ,{\underset{¯}{T}}^d\right)\right)$

$$\mathit{\delta \varphi }\frac{{\partial }^{2}W}{\partial \varphi \partial \varphi }\mathit{\delta \varphi }\ge 0,\mathit{\delta \varphi }\ge 0\mathrm{on}{\mathrm{\Gamma }}^{+},\mathit{\delta \varphi }\ne 0,$$E71

and the uniqueness and non-bifurcation is characterized by

$$\mathit{\delta \varphi }\frac{{\partial }^{2}W}{\partial \varphi \partial \varphi }\mathit{\delta \varphi }\ge 0,\mathit{\delta \varphi }\ne 0\mathrm{on}{\mathrm{\Gamma }}^{+}.$$E72

6.2 Delamination of a thin membrane under pressure

The strain energy $\psi$ of the membrane is given as $\psi \left(u\right)=\frac{1}{2}K{\left(\nabla u\right)}^2$ where $u$ is the transverse displacement as depicted on (Figure 7). The potential energy of the whole system is:

Figure 7.

Delamination of a thin membrane.

$$\mathcal{E}\left(u\left(x,y\right),p\right)={\int }_{\mathrm{\Omega }}\frac{1}{2}K{\left(\nabla u\right)}^2\mathrm{d}a-{\int }_{\mathrm{\Omega }}\mathit{pu}\mathrm{d}a$$E73

The displacement $u=0$ over $\mathrm{\partial \Omega }$. When the boundary $\mathrm{\partial \Omega }$ is moving with normal velocity $\varphi$ the variation of energy determines the associated driving force

$$\partial \mathcal{E}={\int }_{\mathrm{\Omega }}\frac{\partial }{\partial u}\left(\psi -\mathit{pu}\right)\mathit{\delta u}\mathrm{d\Omega }-{\int }_{\mathrm{\partial \Omega }}\left(\psi -pu\right)\mathit{\delta \varphi }\left(s\right)\mathrm{d}a$$E74

where the displacement $\mathit{\delta u}$ is related to the boundary $\mathrm{\partial \Omega }$ which is moving with the velocity $\mathit{\delta \varphi }$. Along the front $u=0$ at each instant, then the variations are linked as:

$$\mathit{\delta u}+\nabla u.\underset{¯}{n}\mathit{\delta \varphi }=0$$E75

In the domain, the variations of the solution satisfies

$$K\mathrm{\Delta }\mathit{\delta u}=0\underset{¯}{x}\in \mathrm{\Omega }$$E76

The driving force $\mathcal{G}$ (satisfies)

$$D_m={\int }_{\mathrm{\partial \Omega }}\mathit{\psi \delta \varphi }\mathrm{d}a={\int }_{\mathrm{\partial \Omega }}\mathcal{G}\left(s\right)\mathit{\delta \varphi }\left(s\right)\mathrm{d}a,\mathcal{G}={\psi }_{\mid \mathrm{\partial \Omega }}$$E77

The variational inequality takes the form

$${\int }_{\mathrm{\partial \Omega }}\frac{\mathrm{d}}{\mathrm{d}t}\left(\mathcal{G}-G_c\right)\left(\mathit{\delta \varphi }-\varphi \right)\mathrm{d}a\ge 0,\forall \mathit{\delta \varphi }\ge 0$$E78

The boundary value problem is given by the functional

$$\mathcal{F}\left(\tilde{v},\tilde{\varphi }\right)={\int }_{\mathrm{\Omega }}\frac{1}{2}K{\left(\nabla v\right)}^2\mathrm{d\Omega }-\sqrt{2{\mathit{KG}}_c}{\int }_{\mathrm{\partial \Omega }}\underset{¯}{n}.\nabla \nabla u.\underset{¯}{n}{\varphi }^2\mathrm{d}a$$E79

$v$ and $\varphi$ are linked by the constrain $v+\varphi \nabla u.\underset{¯}{n}=0$ over $\mathrm{\partial \Omega }$. The evolution of $\mathcal{G}$ is given by

$$\delta \mathcal{G}=K\nabla u.\nabla \mathit{\delta u}+\left(K\nabla u.\nabla \nabla u.\underset{¯}{n}\right)\mathit{\delta \varphi }$$E80

The set of the admissible velocities $v$ is $\mathcal{K}$:

$$\mathcal{K}=\left\{\left(\tilde{v},\tilde{\varphi }\right)|v\left(s\right)+\varphi \left(s\right)\nabla u.\underset{¯}{n}=0,\varphi \ge 0,\mathcal{G}\le G_c,\varphi \left(\mathcal{G}-G_c\right)=0\right\}$$E81

For circular geometry, the displacement solution is $u=\frac{p}{4K}\left(R^2-r^2\right)$ and the propagation is possible when $\mathcal{G}=G_c$ that defines the critical pressure $p_c=\frac{2}{R}\sqrt{2{\mathit{KG}}_c}$. Consider a change of shape by a Fourier expansion

$$\mathit{\delta \varphi }=a_o+\sum _i{a}_{i}cos\left(\mathit{i\theta }\right)+b_{i}sin\left(\mathit{i\theta }\right)$$E82

the associated velocity solution of the rate boundary value problem is

$$v^{\mathit{sol}}=\frac{\mathit{pR}}{2K}\left(a_o+\sum _i\left(a_{i}cos\left(\mathit{i\theta }\right)+b_{i}sin\left(\mathit{i\theta }\right)\right){\left(\frac{r}{R}\right)}^i\right)$$E83

Evaluating the functional $W\left(\varphi \right)=\mathcal{F}\left(v^{\mathit{sol}},\varphi \right)$, the condition of stability is deduced as

$$2\pi G_c\left(-2a_o^2+\sum _i\left(i-1\right)\left(a_i^2+b_i^2\right)\right)\ge 0$$E84

hence the circular shape is unstable for pressure controlled system.

If now the volume is controlled, the pressure becomes the Lagrange multiplier associated with the condition ${\int }_{\mathrm{\Omega }}u\mathrm{d\Omega }=V^d.$ The condition of stability under this loading, becomes

$$2\pi G_c\left(6a_o^2+\sum _i\left(i-1\right)\left(a_i^2+b_i^2\right)\right)\ge 0.$$E85

The stability is ensured, but uniqueness is not, $a_1$ and $b_1$ can be defined such that $\mathit{\delta \varphi }=a_o+a_{1}cos\theta +b_{1}sin\theta \ge 0$. Many other examples are founded in literature for more complicated situations.

7. Conclusions

We have presented an introduction to the analysis of bifurcation and stability during the evolution of nonlinear system governed by potential energy, potential of dissipation and normality rule. This framework is used in elasto-plasticity, in fracture and for moving interfaces.

The rate boundary value problem has a formal identical structure and leads to variational inequalities that the evolution of internal state must satisfy. These inequalities are based on the second derivative of the energy of the system, and are quadratic operators. The properties of these operators give the condition of existence and uniqueness of the system evolution.

Some applications have been presented. Many other situations can be investigated as in phase transformation [19]. This last example shows how the analysis of stability-bifurcation has strong implications in homogenization for the definition of an homogenized constitutive behavior.

The conditions of stability and no-bifurcation can also be used to determine criterion of initiation of defect as pointed out in Refs. [22, 23].