# AN ARC-LENGTH APPROXIMATION FOR ELLIPTICAL ORBITS

Ashim B. Karki1 and Aayush Jha2
1Wolfram Research, Champaign, IL, USA
2Department of Physics, St. Xavier’s College, Maitighar, Kathmandu 44600,
Nepal
E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

In this paper, we overlay a continuum of analytical relations which
essentially serve to compute the arc-length described by a celestial body in
an elliptic orbit within a stipulated time interval. The formalism is based
upon a two-dimensional heliocentric coordinate frame, where both the
coordinates are parameterized as two infinitely differentiable functions in
time by using the Lagrange inversion theorem. The parameterization is firstly
endorsed to generate a dynamically consistent ephemerides for any celestial
object in an elliptic orbit, and thereafter manifested into a numerical
integration routine to approximate the arc-lengths delineated within an
arbitrary interval of time. As elucidated, the presented formalism can also be
orchestrated to quantify the perimeters of elliptic orbits of celestial bodies
solely based upon their orbital period and other intrinsic characteristics.

###### keywords:

Methods: Numerical – Celestial Mechanics – Ephemerides

††pubyear: 2019††pagerange: AN ARC-LENGTH APPROXIMATION FOR ELLIPTICAL
ORBITS–AN ARC-LENGTH APPROXIMATION FOR ELLIPTICAL ORBITS

##  1 Introduction

Over the past few decades, the advancement in astrodynamics has instigated a
new era of celestial ephemerides. With the advent of novel observational
routines - namely, Microwave ranging, lunar laser ranging, VLBI measurements,
etc. - astronomers have not only subsided ambiguity in estimations of orbital
elements of celestial objects but have also deduced results that are
consequently more viable and pragmatic. Throughout the course of its
refinement, data procured from theoretical approaches for determination of
ephemerides has advanced to a machine-level precision with the incorporation
of modern computational software. The refurbishment of classical analytical
theories proposed by Newcomb (Newcomb, 1882), rendered to make
reinvestigations into the Nautical Almanac, eventually came to a halt in the
late 1960s. Inevitably, digital computation was set on par with modern
astronomy when Eckert (1951) pioneered to numerically integrate the equations
of orbital motion of the planets and the Moon, thereafter followed by
Oesterwinter & Cohen (1972). The enhancement course was overtaken by American
groups at JPL, Pasadena (Standish Jr et al., 1976; Newhall et al., 1983) with
the commencement of a series of Development Ephemerides (which took into
account empirical discrepancies, such as perturbations, obliquity, nutations,
and librations of the planets) as the need for more accurate planetary
ephemerides became evident with the advancement in space exploration, and a
plethora of other rectifications have been made in this arena hitherto.
Notwithstanding the meticulously stringent observations made till date, the
contemporary ephemerides are mostly modelled upon the Earth-centred inertial
(commonly, the J2000) coordinate frames for numerical integration, whereas, a
miniscule progression has been made to ameliorate corresponding concepts
reliant upon heliocentric coordinate system. In this paper, we set forth
analytical relations to compute the positions, and ultimately the arc-lengths
of celestial bodies by initially parameterizing the elliptic (in this case,
heliocentric) coordinates as direct functions of time elapsed from the epoch.
The parameterization that follows can be thought of as a rudimentary
application of the open-form solution to the Inverse Kepler's problem as it
makes a predominant use of the Lagrange Inversion theorem to revert the
Kepler's equation and express the eccentric anomaly in terms of mean anomaly
and hence, secure two parametric functions in time. The first order time-
derivatives of these functions, in the light of the inverted series of the
eccentric anomaly, can thereafter communally yield an expression for arc-
length traced by the celestial body in between two arbitrary time periods
referenced from the epoch. As a corroboration to our work, values of the
orbital elements of Earth-Moon barycentre and other celestial bodies -
rectified to an eighth decimal-digit precision - have been deployed from
Murray & Dermott (1999) to numerically integrate the afore-mentioned
expression, and hence compute the perimeter of their elliptic orbits in the
form of an arc-length approximation.
In Section 2, the Kepler's equation is inverted using the Lagrange inversion
theorem, supplemented with functions of eccentricity anomaly and its first
order time-derivative. In Section 3, the heliocentric coordinates are
parametrized in time and ephemerides of Earth-Moon barycentre is tabulated
along with an illustration. Section 4 constructs a general expression to
numerically compute the arc-lengths described by celestial bodies advancing in
elliptic orbits. In Section 5, we demonstrate how this expression can be
effectuated to calculate the perimeters of elliptic orbits. Section 6
scrutinizes about the probable causes of ambiguities in approximations secured
by the equations deduced in this paper. The final verdicts are demarcated in
the last section.

##  2 DISAMBIGUATION OF THE INVERSE KEPLER'S EQUATION

On account of the fact that the Kepler's equation cannot be analytically
solved for eccentric anomaly E, a myriad of fixed-point iterative solutions,
which employ the Newton-Raphson method, have been proposed and revised so far.
However, it is worthwhile to note that the Lagrange inversion theorem is
undeniably a more pertinent approach in this case, as it facilitates the
reversion of the equation by expressing it into an open-form convergent
Maclaurin series in terms of the mean anomaly M, allowing it to be solved a
priori. The accuracy of the series can thereby be optimized by truncating it
as per convenience. The relation between mean anomaly M and the eccentric
anomaly E can be set up in accordance to the Kepler's equation as

$M=h(E)=E-e\sin E,$ (1)

where e is the eccentricity of the elliptic orbit and h(E) represents a
function of M in terms of E. The function is clearly analytic and
differentiable at E=0, and as $\textit{h}^{\prime}(0)=1-e\neq
0,\textit{h}(0)=0$, E can reasonably be enumerated in the form of M as

$E=\sum_{k=1}^{\infty}\frac{M^{k}}{k!}\lim_{E\rightarrow
0^{+}}\left(\frac{\mathrm{d}^{k-1}}{\mathrm{d}E^{k-1}}\left(\left(\frac{E}{E-e\sin
E}\right)^{k}\right)\right),\forall e\neq 1.$ (2)

On evaluation, the equation (2) yields

$E=\frac{M}{1-e}-\frac{e}{(1-e)^{4}}\frac{M^{3}}{3!}+\frac{\left(9e^{2}+e\right)}{(1-e)^{7}}\frac{M^{5}}{5!}\\\
-\frac{\left(225e^{3}+54e^{2}+e\right)}{(1-e)^{10}}\frac{M^{7}}{7!}+\mathcal{O}\left(M^{9}\right),\forall
e\neq 1.$ (3)

The mean anomaly $M$ can also be rendered as $M=nt$, where t is the time
elapsed since pericentre passage and the mean motion111The standard
gravitational parameter is denoted here by
$\mu=\mathcal{G}\left(m_{1}+m_{2}\right)$, where $m_{1}$ and $m_{2}$ are the
two corresponding masses of the bodies in the orbiting system and
$\mathcal{G}$ is the Newtonian Gravitational constant. The orbital period of
the celestial body is denoted by P, a is its semi-major axis.
$n=\sqrt{\mu/a^{3}}=2\pi/P.$
The eccentricity anomaly E can now be deduced in terms of t by reformulating
equation (3) as

$E=\frac{nt}{1-e}-\frac{e}{(1-e)^{4}}\frac{n^{3}t^{3}}{3!}+\frac{\left(9e^{2}+e\right)}{(1-e)^{7}}\frac{n^{5}t^{5}}{5!}\\\
-\frac{\left(225e^{3}+54e^{2}+e\right)}{(1-e)^{10}}\frac{n^{7}t^{7}}{7!}+\mathcal{O}\left(t^{9}\right).$
(4)

Taking the first order time-derivative of equation (4) on a term-by-term basis
leads to

$\dot{E}=\frac{n}{1-e}-\frac{e}{(1-e)^{4}}\frac{n^{3}t^{2}}{2!}+\frac{\left(9e^{2}+e\right)}{(1-e)^{7}}\frac{n^{5}t^{4}}{4!}\\\
-\frac{\left(225e^{3}+54e^{2}+e\right)}{(1-e)^{10}}\frac{n^{7}t^{6}}{6!}+\mathcal{O}\left(t^{8}\right).$
(5)

##  3 PROCESSION OF EARTH IN RESPECT OF HELIOCENTRIC COORDINATE FRAMES

Affixing the Sun at one of the two foci of the elliptic orbit, it can be
derived from general geometry that the heliocentric cartesian coordinates222As
a matter of convenience, the orbital inclination can be set to 0 for
computation and the ascending node can be placed in the reference direction
allowing the coordinate frames to be parsed in a cartesian plane.of the Earth-
Moon barycentre during the procession are

$x=a(\cos E-e),$ (6) $y=b\sin E,$ (7)

where a is the semi-major axis of the elliptical orbit and b is its semi-minor
axis. Substituting the expression for E in equations (6) - (7), we get

$x=a\cos\left[\sum_{k=1}^{\infty}\frac{M^{k}}{k!}\lim_{E\rightarrow
0^{+}}\left(\frac{\mathrm{d}^{k-1}}{\mathrm{d}E^{k-1}}\left(\left(\frac{E}{E-e\sin
E}\right)^{k}\right)\right)\right]-ae,$ (8)
$y=b\sin\left[\sum_{k=1}^{\infty}\frac{M^{k}}{k!}\lim_{E\rightarrow
0^{+}}\left(\frac{\mathrm{d}^{k-1}}{\mathrm{d}E^{k-1}}\left(\left(\frac{E}{E-e\sin
E}\right)^{k}\right)\right)\right],$ (9)

which on further simplification can be written respectively as

$x=a\cos\Bigg{[}\frac{M}{1-e}-\frac{e}{(1-e)^{4}}\frac{M^{3}}{3!}+\frac{\left(9e^{2}+e\right)}{(1-e)^{7}}\frac{M^{5}}{5!}\\\
-\frac{\left(225e^{3}+54e^{2}+e\right)}{(1-e)^{10}}\frac{M^{7}}{7!}+\mathcal{O}\left(M^{9}\right)\Bigg{]}-ae,$
(10)

$y=b\sin\Bigg{[}\frac{M}{1-e}-\frac{e}{(1-e)^{4}}\frac{M^{3}}{3!}+\frac{\left(9e^{2}+e\right)}{(1-e)^{7}}\frac{M^{5}}{5!}\\\
-\frac{\left(225e^{3}+54e^{2}+e\right)}{(1-e)^{10}}\frac{M^{7}}{7!}+\mathcal{O}\left(M^{9}\right)\Bigg{]}.$
(11)

In order to manifest the elliptical orbit as a parametric curve such that the
set of arbitrary points $(x,y)=(f(t),g(t))$ describe a continuous and
differentiable trajectory, $0\leq t\leq P$, one can plug in the series of $E$
secured from equation (4) into equations (6) - (7) to write333It is to be
noted that, from the manner in which the equations have been derived, both $E$
and $M$ should be expressed in circular measure.

$x=f(t)=a\cos\Bigg{[}\frac{nt}{1-e}-\frac{e}{(1-e)^{4}}\frac{n^{3}t^{3}}{3!}+\frac{\left(9e^{2}+e\right)}{(1-e)^{7}}\frac{n^{5}t^{5}}{5!}\\\
-\frac{\left(225e^{3}+54e^{2}+e\right)}{(1-e)^{10}}\frac{n^{7}t^{7}}{7!}+O\left(t^{9}\right)\Bigg{]}-ae,$
(12)
$y=g(t)=b\sin\Bigg{[}\frac{nt}{1-e}-\frac{e}{(1-e)^{4}}\frac{n^{3}t^{3}}{3!}+\frac{\left(9e^{2}+e\right)}{(1-e)^{7}}\frac{n^{5}t^{5}}{5!}\\\
-\frac{\left(225e^{3}+54e^{2}+e\right)}{(1-e)^{10}}\frac{n^{7}t^{7}}{7!}+\mathcal{O}\left(t^{9}\right)\Bigg{]}.$
(13)

Table 1: If one centres the origin (here, the centre of mass of the Sun) at the focus of the elliptical orbit and overlooks perturbative effects of other planets and apsidal procession of the orbit, then the course of maneuver of the Earth-Moon barycentre can be described using equations (12) - (13). The variations in the heliocentric coordinates are quantified below apropos of their corresponding arbitrary time periods $t$ (time elapsed since periastron passage). Note that $t$ has been expressed in Earth days (d), and the horizontal and vertical displacements from the origin have been computed in Astronomical Units (AU). $t$ | $x$ | $y$
---|---|---
30 | 0.849832888 | 0.499037433
60 | 0.485074866 | 0.864872995
90 | -0.013617045 | 0.999852941
120 | -0.513134358 | 0.867961230
150 | -0.880147059 | 0.504387701
180 | -1.016691176 | 0.006185495
210 | -0.886323529 | -0.493667781
240 | -0.523836898 | -0.861751337
270 | -0.025989706 | -0.999819519
300 | 0.474334225 | -0.871016043
330 | 0.843475936 | -0.509718583
360 | 0.983215241 | -0.012370722
Figure 1: The scatter represents the positions of the Earth-Moon barycentre at
varying time periods $t$. At this point, it can be inferred that the Earth-
Moon barycentre reaches back to the periastron after 360 Earth days
(approximately, one year) have passed.

##  4 ARC-LENGTHS OF ELLIPTICAL ORBITS: NUMERICAL INTEGRATION ROUTINE

Contingent upon the relations we have established so far, here we construct a
generalized equation that generates the arc-length covered by a celestial
object in between any two arbitrary times - both referenced from the
periastron passage. The first order time-derivatives of the heliocentric
coordinates can be maintained from equations (6) - (7) as

$\dot{x}=-a\cdot\dot{E}\sin E,$ (14) $\dot{y}=b\cdot\dot{E}\cos E.$ (15)

One can now presume the elliptical orbit to be a parametrically defined curve
continuously differentiable on the interval $0\leq t\leq P$ in consonance with
equations (12) - (13). Also, $f^{\prime}(t)$ and $g^{\prime}(t)$ are not
simultaneously zero, as every ellipse is devoid of cusps or corners. In
account of its geometry, the elliptical orbit is transversed exactly once and
does not double back on itself or reverse its direction of motion, since
$\left(f^{\prime}\right)^{2}+\left(g^{\prime}\right)^{2}>0$ throughout the
interval $0\leq t\leq P$. Therefore, one can satisfactorily define the arc-
length in this case to be

$\ell=\int_{t_{1}}^{t_{2}}\sqrt{\dot{x}^{2}+\dot{y}^{2}}{\mathrm{d}t},$ (16)

where $\ell$ represents the arc-length covered by the celestial object within
two arbitrary intervals of time $t_{1}$ and $t_{2}$ commencing from the
periastron passage. From equations (14) - (15) one can modify equation (16)
into

$\ell=\int_{t_{1}}^{t_{2}}\dot{E}\sqrt{a^{2}\sin^{2}E+b^{2}\cos^{2}E}{\mathrm{d}t}.$
(17)

From general geometry, the relationship between semi-major axis $a$, semi-
minor axis $b$, and the eccentricity $e$ of an ellipse is

$e=\sqrt{1-\frac{b^{2}}{a^{2}}}.$ (18)

One can modulate equation (17) purely in terms of eccentricity $e$ and semi-
major axis $a$ using equation (18) as

$\ell=a\int_{t_{1}}^{t_{2}}\dot{E}\sqrt{1-e^{2}\cos^{2}E}{\mathrm{d}t}.$ (19)

Using equation (2), one can write equation (19) as

$\ell=a\int_{t_{1}}^{t_{2}}\frac{\mathrm{d}}{\mathrm{d}t}\left(\sum_{k=1}^{\infty}\frac{M^{k}}{k!}\lim_{E\rightarrow
0^{+}}\left(\frac{\mathrm{d}^{k-1}}{\mathrm{d}E^{k-1}}\left(\left(\frac{E}{E-e\sin
E}\right)^{k}\right)\right)\right)\cdot\\\
\sqrt{1-e^{2}\cos^{2}\left[\sum_{k=1}^{\infty}\frac{M^{k}}{k!}\lim_{E\rightarrow
0^{+}}\left(\frac{\mathrm{d}^{k-1}}{\mathrm{d}E^{k-1}}\left(\left(\frac{E}{E-e\sin
E}\right)^{k}\right)\right)\right]}{\mathrm{d}t}.$ (20)

On expansion of equation (20) and successive substitution of definition of
$M=nt$ (see also equations 4 and 5), one can get

$\ell=a\int_{t_{1}}^{t_{2}}\bigg{(}\frac{n}{1-e}-\frac{e}{(1-e)^{4}}\frac{n^{3}t^{2}}{2!}+\mathcal{O}\left(t^{4}\right)\bigg{)}\cdot\\\
\sqrt{1-e^{2}\cos^{2}\left[\frac{nt}{1-e}-\frac{e}{(1-e)^{4}}\frac{n^{3}t^{3}}{3!}+\mathcal{O}\left(t^{5}\right)\right]}{\mathrm{d}t}.$
(21)

Hereafter, one can numerically integrate equation (21) to estimate the length
of the arc advanced by any celestial object in its elliptical orbit. In order
to comply with the Second law of Kepler, the equations have been preliminarily
framed in a way that numerically integrated arc-lengths are relatively higher
when the celestial body is closer to the periapsis than when it is approaching
the apoapsis for equal durations of time; the effect is ensued as a result of
tangential velocity variations throughout the course of its orbit. (see also
table 2).

Table 2: The arc-lengths delineated by the Earth-Moon barycentre at discrete intervals of time $\Delta t$ for approximately one half of the elliptic orbit (containing the apoapsis) are demonstrated. Within these eighteen different time intervals (instigating from the 90th day to the 270th day), one can notice that, even though the time intervals are evenly spaced, the arc-lengths first gradually diminish until the Earth-Moon barycentre has reached the apoapsis (around 180th day), and then moderately increase as the intervals proceed. Time-interval $\Delta t$ (Earth days) | Arc-length (AU)
---|---
90 to 100 | 0.173978721
100 to 110 | 0.173977334
110 to 120 | 0.173974685
120 to 130 | 0.173971091
130 to 140 | 0.173966982
140 to 150 | 0.173963285
150 to 160 | 0.173959195
160 to 170 | 0.173956450
170 to 180 | 0.173954945
180 to 190 | 0.173954862
190 to 200 | 0.173956209
200 to 210 | 0.173958827
210 to 220 | 0.173962399
220 to 230 | 0.173966500
230 to 240 | 0.173970637
240 to 250 | 0.173974314
250 to 260 | 0.173977090
260 to 270 | 0.173978633

##  5 RENDITION OF THE CIRCUMFERENCE OF THE ELLIPTIC ORBIT AS AN ARC-LENGTH

The parametrically defined arc-length obtained from equation (21) must result
in the total circumference of the elliptical orbit when numerically integrated
within the limits 0 to $P$, where $P$ is the orbital period of the celestial
body. If the origin is set at the elliptic centre, one can parametrically
define the elliptic orbit in polar form to be

$u=a\sin\theta,$ (22) $v=b\cos\theta,$ (23)

where $a$ is its semi-major axis, and $b$ is its semi-minor axis such that
$a>b$ and $0\leq\theta\leq 2\pi$. Then,

$\left(\frac{\mathrm{d}u}{\mathrm{d}\theta}\right)^{2}+\left(\frac{\mathrm{d}v}{\mathrm{d}\theta}\right)^{2}=a^{2}\cos^{2}\theta+b^{2}\sin^{2}\theta,$
(24)

which from equation (18) can be written as

$\left(\frac{\mathrm{d}u}{\mathrm{d}\theta}\right)^{2}+\left(\frac{\mathrm{d}v}{\mathrm{d}\theta}\right)^{2}=a^{2}\left(1-e^{2}\sin^{2}\theta\right),$
(25)

where $e$ is the orbit's eccentricity. Then, the complete elliptic integral of
second kind (perimeter of the ellipse), $\wp$, can be given as

$\wp=4a\int_{0}^{\pi/2}\sqrt{1-e^{2}\sin^{2}\theta}{\mathrm{d}\theta}.$ (26)

The integral in equation (26) can be evaluated after binomial expansion of
$\sqrt{1-e^{2}\sin^{2}\theta}$ to

$\wp=2\pi a\left[1-\left(\frac{1}{2}\right)^{2}e^{2}-\left(\frac{1\cdot
3}{2\cdot 4}\right)^{2}\frac{e^{4}}{3}-\left(\frac{1\cdot 3\cdot 5}{2\cdot
4\cdot 6}\right)^{2}\frac{e^{6}}{5}-\mathcal{O}\left(e^{8}\right)\right].$
(27)

Finally, the celestial body will have covered an arc-length equal to the
circumference of the elliptical orbit, so $\ell=\wp$, when numerically
integrated within the limits 0 to $P$. From equations (21)-(27), one can write

$\int_{0}^{P}\bigg{(}\frac{n}{1-e}-\frac{e}{(1-e)^{4}}\frac{n^{3}t^{2}}{2!}+\mathcal{O}\left(t^{4}\right)\bigg{)}\cdot\\\
\sqrt{1-e^{2}\cos^{2}\left[\frac{nt}{1-e}-\frac{e}{(1-e)^{4}}\frac{n^{3}t^{3}}{3!}+\mathcal{O}\left(t^{5}\right)\right]}{\mathrm{d}t}\\\
=2\pi\left[1-\left(\frac{1}{2}\right)^{2}e^{2}-\left(\frac{1\cdot 3}{2\cdot
4}\right)^{2}\frac{e^{4}}{3}-\mathcal{O}\left(e^{6}\right)\right].$ (28)

Table 3: The circumferences of the elliptic orbits of five celestial bodies (with the least orbital inclinations) are computed and enlisted below alongside their deviations from standard measurements.  Celestial body | Circumference (AU) | Uncertainty (%)
---|---|---
Earth | 6.354241 | 1.137%
Neptune | 190.728037 | 0.955%
Adrastea (satellite) | 0.005370 | 0.866%
Metis (satellite) | 0.005376 | 0.013%
Atlas (satellite) | 0.005787 | 0.113%

##  6 APSIDAL PROCESSION AND TRUNCATION ERRORS

It is customary that the argument of periapsis in every elliptical orbit must
experience a steady procession over time resulted by either general
relativistic effects, rotational quadrupole bulges on the planet, stellar
quadrupole moments, Lidov-Kozai mechanism, or perturbative effects of other
planets (Kipping, 2011). The work presented so far relies upon an isolated
orbiting system; however, perturbations can ensue numerous, yet insignificant,
discrepancies onto the existent conditions.
The inconsistencies induced in the computed estimates tend to vanish off when
orbits with low inclinations are designated. The series can be truncated to
any degree of accuracy, and the use of robust computational software should
suffice to curtail any obscurity in calculations.

##  7 Conclusions

Parameterization of heliocentric coordinates as functions of angular
quantities is unequivocally a straightforward task at hand, as they can be
reformulated using elementary geometrical relations. However, this paper
devises a new framework to revisit such conventionality by parameterizing the
coordinates fundamentally in terms of time. This approach, in assent with the
Lagrange inversion theorem, is then reinforced into a set of analytical
relations predominantly governed by time. Complementary to the interposition
of a new module for ephemerides calculations, this procedure also doubles as a
numerical integration routine to compute arc-lengths traced by celestial
bodies advancing in elliptical orbits.

## Acknowledgements

The author acknowledges the support conferred by Wolfram Research Inc.,
Champaign, Illinois, USA. Tabulated and graphical data presented throughout
this paper has been exclusively computed using Mathematica.

## References

  * Eckert (1951) Eckert W. J., 1951, [United States. Nautical Almanac Office. Astronomical paper; v. 12 (1951)],[Washington: US Nautical Almanac Office, 1951], xxxi 327 p. 29 cm., 12, 29
  * Kipping (2011) Kipping D., 2011, The Transits of Extrasolar Planets with Moons. Springer Theses, Springer Berlin Heidelberg
  * Murray & Dermott (1999) Murray C. D., Dermott S. F., 1999, Solar system dynamics. Cambridge university press
  * Newcomb (1882) Newcomb S., 1882, [United States. Nautical Almanac Office. Astronomical paper; v. 1, pt. 6 (1882)],[Washington: US Nautical Almanac Office, 1882], p. 363-487: charts; 29 cm., 1, 363
  * Newhall et al. (1983) Newhall X., Standish E. M., Williams J. G., 1983, Astronomy and astrophysics, 125, 150
  * Oesterwinter & Cohen (1972) Oesterwinter C., Cohen C. J., 1972, Celestial Mechanics, 5, 317
  * Standish Jr et al. (1976) Standish Jr E., Keesey M. S., Newhall X., 1976


# A non incremental variational principle for brittle fracture

Géry de Saxcé

Univ. Lille, CNRS, Centrale Lille, UMR 9013 – LaMcube -
Laboratoire de mécanique multiphysique multiéchelle
F-59000, Lille, France, e-mail<EMAIL_ADDRESS>
###### Abstract

The aim of the paper is to propose a paradigm shift for the variational
approach of brittle fracture. Both dynamics and the limit case of statics are
treated in a same framework. By contrast with the usual incremental approach,
we use a space-time principle covering the whole loading and crack evolution.
The emphasis is given on the modelling of the crack extension by the internal
variable formalism and a dissipation potential as in plasticity, rather than
Griffith’s original approach based on the surface area. The new formulation
appears to be more fruitful for generalization than the standard theory.

Keywords: Linear elastic fracture mechanics, brittle material, symplectic
mechanics, calculus of variation

##  1 Introduction

Many dynamical systems are subjected to energy loss resulting from
dissipation, for instance collisions, surface friction, viscosity, plasticity,
fracture and damage. On the other hand, deformations of solids and motions of
fluids are modeled through constitutive laws. Due to collisions, brittle
fracture and threshold effects, most dissipative laws are non smooth and
multivalued. Moreover, experimental testing suggests that convexity is a
keystone property of these phenomenological laws. Continuum mechanics with
internal variables provides a convenient foundation to develop constitutive
models that describe the inelastic behavior of various materials, despite the
large differences in their physical nature and the relevant scale. The
phenomenological law link the vector of internal variable rates
$\dot{\bm{\alpha}}$ to the dual force vector $\bm{A}$. The theory of
Generalized Standard Materials [41] is based on an hypothesis of normal
dissipation.

$\dot{\bm{\alpha}}\in\partial\varphi(\bm{A})$

where occurs the subdifferential of a convex and lower semicontinuous function
$\varphi$, not everywhere differentiable. Equivalently, the law reads

$\bm{A}=\mbox{argmax}\,(\left\langle\dot{\bm{\alpha}},\bm{A}\right\rangle-\varphi(\bm{A}))$

or, introducing the Fenchel polar $\varphi^{*}$ [33]

$(\dot{\bm{\alpha}},\bm{A})=\mbox{argmin}\,(\varphi(\bm{A})+\varphi^{*}(\dot{\bm{\alpha}})-\left\langle\dot{\bm{\alpha}},\bm{A}\right\rangle)$

that does not favour one of the two dual variables. In the present work, the
latest formulation is chosen as starting point, leading to variational methods
and unconstrained or constrained optimization problems. Among them, Brezis-
Ekeland-Nayroles principle [13] [63], or in short BEN principle, is based on
the time integration of the sum of dissipation potential $\phi$ and its
Fenchel polar (analogous of Legendre polar for convex functions). Although
used a few in the literature, this principle is noteworthy in the sense that
it allows covering the whole evolution of the dissipative system.

In a previous paper [15], Buliga proposed the formalism of Hamiltonian
inclusions, able to model dynamical systems with $1$-homogeneous dissipation
potential (for laws such as brittle damage using Ambrosio-Tortorelli
functional [3]). This formalism is a dynamical version of the quasi-static
theory of rate-independent systems of Mielke ([55], [56], [57]).

Latter on, Buliga and the author merged in a symplectic framework this
formalism with BEN principle to extend it to dynamics [16]. The key-idea is to
decompose additively the time rate $\dot{z}$ into reversible part
$\dot{z}_{R}$ (the symplectic gradient) and dissipative or irreversible one
$\dot{z}_{I}$, next to define the symplectic subdifferential
$\partial^{\omega}\phi(z)$ of the dissipation potential. To release the
restrictive hypothesis of $1$-homogeneity (in particular to address
viscoplasticity), we introduce in this work the symplectic Fenchel polar
$\phi^{*\omega}$, that allows to build theoretical methods to model and
analyse dynamical dissipative systems in a consistent geometrical framework
with the numerical approaches not very far in the background. Numerical
simulation with the BEN principle were performed for elastoplastic structures
in statics [17, 19] and in dynamics [18].

Closer to this approach, we can cite the contributions of Aubin [5], Aubin,
Cellina and Nohel [4], Rockafellar [66], Stolz [70], which considered various
extensions of Hamiltonian and Lagrangian mechanics. In the article [7] by
Bloch, Krishnaprasad, Marsden and Ratiu, Hamiltonian systems are explored with
an added Rayleigh dissipation. A theory of quasistatic rate-independent
systems is proposed by Mielke and Theil [55], Mielke [56], and developed
towards applications in many papers, among them Mielke and Roubíček [57], see
also Visintin [75]. In ([39], [64]), Grmela and Öttinger proposed the
framework GENERIC (General Equation for Non-Equilibrium Reversible-
Irreversible Coupling), a systematic method to derive thermodynamically
consistent evolution equations. A variational formulation of GENERIC is
proposed in [53]. In [58], Mielke proposed a GENERIC formulation for
Generalized Standard Materials quite similar to the one in the present paper.
It would be worth to make the present formalism thermodynamic. In [69],
Stefanelli used Brezis-Ekeland-Nayroles variational principle to represent the
quasistatic evolution of an elastoplastic material with hardening in order to
prove the convergence of time and space-time discretizations as well as to
provide some possible a posteriori error control. Finally, Ghoussoub and
MacCann characterized the path of steepest descent of a non-convex potential
as the global minimum of Brezis-Ekeland-Nayroles functional [37].

Moreover, another advantage of Brezis-Ekeland-Nayroles principle is the
easiness to be generalized. Indeed, it is worth to know that many realistic
dissipative laws, called non-associated, cannot be cast in the mould of the
standard ones deriving of a dissipation potential. To skirt this pitfall, the
author proposed in [24] a new theory based on a function called bipotential.
It represents physically the dissipation and generalizes the sum of the
dissipation potential and its Fenchel polar, reason for which extension of
Brezis-Ekeland-Nayroles principle is natural. The applications of the
bipotential approach to solid Mechanics are various: Coulomb’s friction law
[25], non-associated Drücker-Prager [26] and Cam-Clay models [78] in Soil
Mechanics, cyclic Plasticity ([25],[9], [10], [54], [11]) and Viscoplasticity
[44] of metals with non linear kinematical hardening rule, Lemaitre’s damage
law [8], the coaxial laws ([27],[73]). Such kind of materials are called
implicit standard materials. A synthetic review of these laws can be found in
the two later references. It is also worth to notice that monotone laws but
which does not admit a convex potential can be represented by Fitzpatrick’s
function [35] which is a bipotential.

There is an abundant literature on the crack propagation criteria and our
intention is not to draw here an exhaustive picture. For a deep survey, the
reader is referred to [61]. Let us cite only the most popular ones: the
maximum tensile hoop stress criterion formulated in 1963 by Cherepanov [21],
Erdogan and Sih [32], the criterion of minimum strain energy density by Sih
[68] in 1973, the principle of local symmetry formulated in 1974 by Goldstein
and Salganik [38], the maximum energy release rate criterion considered in
1974 by Hussain et al. [45] and in 1978 by Wu [76], and the generalized
maximum energy release rate criterion [42, 40, 47, 48]. These criteria give
predictions more or less closed to the experimental results and the possible
link with the variational approach is puzzling, at least apparently. This is
one of the challenge we have to address.

Our aim here is to study the propagation of already initiated cracks. Then we
exclude of this paper the problem of the crack nucleation or initiation that
would require extra information about stresses in the damaged region. Two
models are used in the literature, the Coupled Criterion proposed par
Leguillon [50] and the Cohesive Zone Model, originated in the pioneer works by
Barenblatt [6] and Dugdale [30], next developed by Tvergaard and al [72], Xu
[77] and, for the variational aspects, by Bourdin et al [12].

The paper is organized as follows.

  * •

In Section 2, the initial crack and the loading history being known, the
problem is to find the crack evolution. Instead of the usual modelling of the
crack extension as a time-parameterized family of surfaces, we propose to
introduce a flow of a vector field on the final crack manifold.

  * •

Section 3 is devoted to the definition of the variable of the problem in the
symplectic framework: the displacement and flow fields and the corresponding
dynamic momenta.

  * •

In Section 4, the symplectic BEN principle is applied to the fracture
mechanics, that leads to the definition of a driving force, dual of the crack
flow rate.

  * •

In Section 5, we propose a constitutive law for brittle fracture in the form
of a crack stability criterion and the normality law for the crack extension
rule.

  * •

We prove in Section 6 for structure with uniform toughness the equivalence
with the classical variational approach in the form of a Munford-Shah
functional.

  * •

In Section 7, we deduce from the space-time principle an incremental version
suitable for the usual step-by-step numerical approaches.

  * •

In Section 8, we show how a suitable interpretation of the experimental data
prove the relevancy of the normality law and the link with the principle of
local symmetry.

  * •

Section 9 is devoted to the calculation of the crack driving force using the
calculus of variation on the jet space of order one in the general case of
dynamics.

##  2 Modelling of the crack extension

###  2.1 Data

Figure 2.1: Crack Flow

Let $\displaystyle\Omega\subset\mathbb{R}^{n}$ be a bounded, open set, with
piecewise smooth boundary $\partial\Omega$ occupied by the uncracked body
(Figure 2.1). Our aim is to study a crack extension during a time interval
$\left[0,T\right]$

The data of the problem is:

  * •

the initial crack $\Gamma_{0}$

  * •

the imposed displacements $\bar{\bm{u}}$ on the part $\partial\Omega_{0}$ of
the boundary (the supports)

  * •

the surface forces $\bar{\bm{f}}$ on the remaining part $\partial\Omega_{1}$

  * •

the volume forces $\bm{f}$

All the external actions, $\bar{\bm{u}},\bar{\bm{f}}$ and $\bm{f}$, are
prescribed during the time interval. No stress is transmitted through the
crack.

###  2.2 Crack

A crack bifurcates when its shape suffers a change of topology. The most
common example is a crack which develops in time new branches. During this
phenomenon the number of crack fronts increases. In this paper we do not study
the bifurcation of an existing crack.

We denote $S$ the set of admissible surfaces, i.e. closed countably $2$
rectifiable subsets of $\Omega$ without change of topology. The crack
evolution is described by a time-parameterized family of surfaces (Figure 2.1)

$\Gamma:\left[0,T\right]\rightarrow S:t\mapsto\Gamma_{t}=\Gamma(t)$

such that $\Gamma(0)=\Gamma_{0}$ and the map $\Gamma$ is monotone increasing,
$S$ being equipped with the inclusion order. The crack front $c_{t}$ at time
$t$ is parameterized by the arc length $s$.

The cornerstone of the formulation is to model the crack extension by a flow
on the cracked surface $\Gamma_{T}\backslash\Gamma_{0}$ during the time
interval

$\left[0,T\right]\times
c_{0}\rightarrow\Gamma_{T}\backslash\Gamma_{0}:(t,\bm{x}_{0})\mapsto\bm{x}=\bm{\psi}(t,\bm{x}_{0})$
(1)

such that

$\bm{\psi}(\\{t\\}\times c_{0})=c_{t},\qquad\bm{\psi}(\left[0,t\right]\times
c_{0})=\Gamma_{t}\backslash\Gamma_{0}$

It is worth to remark that:

  * •

From a calculus viewpoint, the field $\bm{\psi}$ represents both the crack and
its evolution. The main advantage is that it is easier working with fields
living in a functional space than with surfaces.

  * •

In particular, our approach is well suited for numerical applications.
$s\mapsto\bm{x}_{0}=\bm{f}(s)$ being the parameterization of the initial crack
$c_{0}$ by the arc length $s$, the cracked surface can be parameterized
according to
$(t,s)\mapsto\bm{x}=\tilde{\bm{\psi}}(t,s)=\bm{\psi}(t,\bm{f}(s))$. The node
$\bm{x}_{jk}=\tilde{\bm{\psi}}(t_{j},s_{k})$ being the position at time
$t_{j}$ of the node $\bm{x}_{0k}=\bm{f}(s_{k})\in c_{0}$ and $(s,t)\mapsto
N_{jk}(t,s)$ the corresponding shape function, the cracked surface
$\Gamma_{T}$ is parameterized by

$\bm{x}=\tilde{\bm{\psi}}(t,s)=\sum_{j}\sum_{k}N_{jk}(t,s)\,\bm{x}_{jk}$ (2)

In the forthcoming variational formulation, the nodal values
$\bm{x}_{jk}=\tilde{\bm{\psi}}(t_{j},s_{k})$ will be unknowns of the
discretized problem.

  * •

A streamline $t\mapsto\bm{x}=\bm{\psi}(t,\bm{x}_{0})$ does not correspond to
the motion of a material particle initially at the position $\bm{x}_{0}\in
c_{0}$. The flow is not Lagrangian but Eulerian in the sense it represents
just the evolution of the crack front.

If the cracked surface is smooth enough, it is generated by a vector field
$\dot{\bm{\psi}}$, solution on $\Gamma_{T}\backslash\Gamma_{0}$ of the ODE
(Figure 2.1)

$\frac{d}{dt}(\bm{\psi}(t,\bm{x}_{0}))=\dot{\bm{\psi}}(\bm{\psi}(t,\bm{x}_{0}))\qquad\mbox{with}\quad\bm{\psi}(0,\bm{x}_{0})=\bm{x}_{0}$

As in the FEM method, it is easy to compute the vector field $\dot{\bm{\psi}}$
as the partial derivative of (2) with respect to $t$.

##  3 Variables of the problem

To do not enter into too complex formulations, we consider the case of an
elastic body in small strains, although the generalization to the
elastoplasticity and the finite strains is rather straightforward. The cracked
solid at time $t$ is denoted $\Omega_{t}=\Omega\backslash\Gamma_{t}$. Within
the body, the displacement field at time $t$ is

$\Omega_{t}\rightarrow\mathbb{R}^{3}:\bm{x}\mapsto\bm{u}(t,\bm{x})$

According to the symplectic BEN principle [16], the variables of the problem
of the elastodynamics are:

  * •

the field couple $\bm{\xi}=(\bm{u},\bm{\psi})$ of the displacement and the
flow

  * •

the momenta $\bm{\eta}=(\bm{p},\bm{\pi})$ where $\bm{p}$ is the classical
linear momentum and $\bm{\pi}$ is the momentum associated to the flow

The system evolution is given by a curve

$t\mapsto\bm{z}(t)=(\bm{\xi}(t),\bm{\eta}(t))$

##  4 The symplectic BEN variational principle

The duality between the space $X$ of degrees of freedom $\bm{\xi}$ and the
space $Y$ of momenta $\bm{\eta}$ has the form

$\langle\dot{\bm{\xi}},\dot{\bm{\eta}}\rangle=\langle\dot{\bm{u}},\dot{\bm{p}}\rangle+\langle\dot{\bm{\psi}},\dot{\bm{\pi}}\rangle=\int_{\Omega_{t}}\dot{\bm{u}}\cdot\dot{\bm{p}}\mbox{
d}^{3}x+\int_{c_{t}}\dot{\bm{\psi}}\cdot\dot{\bm{\pi}}\mbox{ d}s$

The symplectic form is the $2$-form field on $X\times Y$ defined by

$\omega(\dot{\bm{z}},\dot{\bm{z}}^{\prime})=\omega((\dot{\bm{\xi}},\dot{\bm{\eta}}),(\dot{\bm{\xi}}^{\prime},\dot{\bm{\eta}}^{\prime}))=\langle\dot{\bm{\xi}},\dot{\bm{\eta}}^{\prime}\rangle-\langle\dot{\bm{\xi}}^{\prime},\dot{\bm{\eta}}\rangle$

The total Hamiltonian of the structure is taken of the integral form

$H_{t}=H(t,\bm{z})=\int_{\Omega_{t}}\left\\{\dfrac{1}{2\rho}\parallel\bm{p}\parallel^{2}+w(\nabla_{s}\bm{u})-\bm{f}(t)\cdot\bm{u}\right\\}\mbox{
d}^{3}x-\int_{\partial\Omega_{1}}\bar{\bm{f}}(t)\cdot\bm{u}\mbox{
d}(\partial\Omega_{1})$ (3)

The first term is the kinetic energy. The second one is the elastic strain
energy $w$ depending on the strain tensor $\bm{\varepsilon}=\nabla_{s}\bm{u}$
and the two latter terms are the works of the external forces. As usual, the
stress tensor is

$\bm{\sigma}=\frac{\partial w}{\partial\bm{\varepsilon}}$

A crack being a material discontinuity, we must distinguish two material
surfaces $\Gamma^{+}_{t}$ and $\Gamma^{-}_{t}$ that occupy the same position
as $\Gamma_{t}$ but are the two sides of the crack and have opposite unit
normal vectors, exterior to $\Omega_{t}$: $\bm{n}^{+}=-\bm{n}^{-}$.

A curve is said admissible if it satisfies

  * •

the boundary conditions:

$\forall
t\in\left[0,T\right],\quad\bm{u}=\bar{\bm{u}}\;\;\mbox{on}\;\;\partial\Omega_{0},\quad\bm{\sigma}\cdot\bm{n}=\bar{\bm{f}}\;\;\mbox{on}\;\;\partial\Omega_{1},\quad\bm{\sigma}\cdot\bm{n}^{\pm}=\bm{0}\;\;\mbox{on}\;\;\Gamma^{\pm}_{t}$

  * •

the initial conditions: $\bm{z}(0)=\bm{z}_{0}$

The SBEN formalism of dissipative media is based on the additive decomposition
of the velocity into reversible and irreversible parts

$\dot{\bm{z}}=\dot{\bm{z}}_{R}+\dot{\bm{z}}_{I}$

As the variables of the problem are fields, we use functional derivatives

$\langle\delta\bm{\xi},D_{\bm{\xi}}H\rangle+\langle
D_{\bm{\eta}}H,\delta\bm{\eta}\rangle=\lim_{\epsilon\to
0}\frac{1}{\epsilon}(H(t,\bm{z}+\epsilon\,\delta\bm{z})-H(t,\bm{z}))$

The reversible part of the velocity is given by the symplectic gradient of the
Hamiltonian (or Hamiltonian vector field)

$\dot{\bm{z}}_{R}=\nabla^{\omega}H=(D_{\bm{\eta}}H,-D_{\bm{\xi}}H)=((D_{\bm{p}}H,D_{\bm{\pi}}H),(-D_{\bm{u}}H,-D_{\bm{\psi}}H))$

For the Hamiltonian (3), we obtain

$D_{\bm{p}}H=\frac{\bm{p}}{\rho},\qquad D_{\bm{\pi}}H=\bm{0}$
$-D_{\bm{u}}H=\nabla\cdot\bm{\sigma}+\bm{f},\qquad-D_{\bm{\psi}}H=\bm{G}$

The last equation is just a definition of the driving force $\bm{G}$ of which
the explicit expression will be discussed afterwards.

In [16], we formulated a general variational principle for the dynamical
dissipative systems which claims that

Symplectic BEN principle 1. The natural evolution of the system minimizes the
functional

$\Pi(\bm{z})=\int_{0}^{T}\left\\{\phi(\dot{\bm{z}})+\phi^{*\omega}(\dot{\bm{z}}_{I})-\omega(\dot{\bm{z}}_{I},\dot{\bm{z}})\right\\}\mbox{
dt}$ (4)

among the admissible curves $t\mapsto\bm{z}(t)$ and the minimum is zero.

In this principle occurs the symplectic form $\omega$, the Hamiltonian through
the irreversible part of the velocity $\dot{z}_{I}=\dot{z}-\nabla^{\omega}H$,
a convex dissipation potential $\phi$ and its symplectic polar
$\phi^{*\omega}$. We calculate now the detailed expression of this ingredients
for the brittle fracture. We start with a simplifying hypothesis similar to
the one introduced in [16] to recover the classical elastoplasticity. We claim
that in the dissipation potential, the variables other than $\dot{\bm{\pi}}$
are ignorable

$\phi(\dot{\bm{z}})=\varphi(\dot{\bm{\pi}})$

The symplectic polar is defined in [16] as

$\phi^{*\omega}(\dot{\bm{z}}^{\prime})=\sup\left[\omega(\dot{\bm{z}}^{\prime},\dot{\bm{z}})-\phi(\dot{\bm{z}})\;:\;\dot{\bm{z}}\right]$

Under the previous hypothesis, one has

$\phi^{*\omega}(\dot{\bm{z}}^{\prime})=\sup\left[\left\langle\dot{\bm{u}}^{\prime},\dot{\bm{p}}\right\rangle+\left\langle\dot{\bm{\psi}}^{\prime},\dot{\bm{\pi}}\right\rangle-\left\langle\dot{\bm{u}},\dot{\bm{p}}^{\prime}\right\rangle-\left\langle\dot{\bm{\psi}},\dot{\bm{\pi}}^{\prime}\right\rangle-\varphi(\dot{\bm{\pi}})\;:\;\dot{\bm{z}}\right]$

that gives

$\phi^{*\omega}(\dot{\bm{z}}^{\prime})=\chi_{\left\\{\bm{0}\right\\}}(\dot{\bm{u}}^{\prime})+\varphi^{*}(\dot{\bm{\psi}}^{\prime})+\chi_{\left\\{\bm{0}\right\\}}(\dot{\bm{p}}^{\prime})+\chi_{\left\\{\bm{0}\right\\}}(\dot{\bm{\pi}}^{\prime})$

where $\chi_{K}$ is the indicatory function of the set $K$, equal to $0$ in
$K$ and $+\infty$ otherwise, and $\varphi^{*}$ is the Fenchel polar of
$\varphi$. The second term in the functional (4) becomes:

$\phi^{*\omega}(\dot{\bm{z}}_{I})=\chi_{\left\\{\bm{0}\right\\}}\left(\dot{\bm{u}}-\frac{\bm{p}}{\rho}\right)+\varphi^{*}(\dot{\bm{\psi}})+\chi_{\left\\{\bm{0}\right\\}}(\dot{\bm{p}}-\nabla\cdot\bm{\sigma}-\bm{f})+\chi_{\left\\{\bm{0}\right\\}}(\dot{\bm{\pi}}-\bm{G})$

As the minimum of the functional is the finite value zero, it will be reached
when the arguments of the indicatory functions vanish. Then we may take the
zero value of the indicatory functions in the functional while we introduce
extra corresponding constraints

$\bm{p}=\rho\,\dot{\bm{u}},\qquad\nabla\cdot\bm{\sigma}+\bm{f}=\dot{\bm{p}},\qquad\dot{\bm{\pi}}=\bm{G}$
(5)

Next we can transform the last term in the functional, taking into account the
velocity decomposition, the linearity and the antisymmetry of the symplectic
form

$-\omega(\dot{\bm{z}}_{I},\dot{\bm{z}})=-\omega(\dot{\bm{z}}-\dot{\bm{z}}_{R},\dot{\bm{z}})=-\omega(\dot{\bm{z}},\dot{\bm{z}})+\omega(\dot{\bm{z}}_{R},\dot{\bm{z}})=\omega(\dot{\bm{z}}_{R},\dot{\bm{z}})$
(6)

or in detail

$-\omega(\dot{\bm{z}}_{I},\dot{\bm{z}})=\left\langle\frac{\bm{p}}{\rho},\dot{\bm{p}}\right\rangle-\left\langle\dot{\bm{u}},\nabla\cdot\bm{\sigma}+\bm{f}\right\rangle-\left\langle\dot{\bm{\psi}},\bm{G}\right\rangle$

Owing to the constraints (5), it holds

$-\omega(\dot{\bm{z}}_{I},\dot{\bm{z}})=\left\langle\dot{\bm{u}},\dot{\bm{p}}-\nabla\cdot\bm{\sigma}-\bm{f}\right\rangle-\left\langle\dot{\bm{\psi}},\bm{G}\right\rangle=-\left\langle\dot{\bm{\psi}},\bm{G}\right\rangle$

Moreover, thanks to the former and latter constraints in (5), the momenta can
be eliminated from the functional and the intermediate constraint becomes

$\nabla\cdot\bm{\sigma}+\bm{f}=\rho\,\ddot{\bm{u}}$

while the initial condition on the linear momentum can be transformed into an
initial condition on the velocity because $\bm{p}(0)=\rho\,\dot{\bm{u}}(0)$.
Taking into account these transformations, the minimum can be searched only on
the space of the degrees of freedom $\bm{u}$ and $\bm{\psi}$. Then we obtain a
second version of the variational principle

Symplectic BEN principle 2. The natural evolution of the system minimizes

$\Pi(\bm{\xi})=\int_{0}^{T}\left\\{\varphi(\bm{G})+\varphi^{*}(\dot{\bm{\psi}})-\left\langle\dot{\bm{\psi}},\bm{G}\right\rangle\right\\}\mbox{
d}t$ (7)

among the admissible curves $t\mapsto\bm{\xi}(t)$ such that
$\nabla\cdot\bm{\sigma}+\bm{f}=\rho\,\ddot{\bm{u}}$ and the minimum is zero.

It could seem puzzling that the displacement does not appear explicitly in the
expression of the functional but in fact $\bm{u}$ and $\bm{\psi}$ are coupled
in the minimization problem:

  * •

$\bm{\psi}$ is controlled by $\bm{G}$ which depends on $\bm{u}$ as derivative
of the Hamiltonian.

  * •

$\bm{u}$ must satisfy the constraint and the admissibility conditions that are
defined on $\Omega_{t}=\Omega\backslash\Gamma_{t}$ controlled by $\bm{\psi}$.

##  5 Constitutive laws for brittle fracture

As discussed in the Introduction, most of the crack stability criteria in the
literature, in particular the most popular, are of local nature and not
variational. Nevertheless, Strifors proposed in [71] to predict the onset of
brittle fracture and the extension direction thanks to the crack extension
force that can be identified to the driving force $\bm{G}$. Now we show how to
cast it into the mold of our variational approach. Based on the works by
Strifors [71] and Hellen et al. [43], we propose a simple constitutive law. At
least for plane problems, it seems reasonable to think that the influence on
the fracture extension of the projection of $\bm{G}$ onto the tangent to the
crack front can be neglected. Then we introduce the deviatoric force

$\bm{G}_{\perp}=\bm{G}-(\bm{G}\cdot\bm{\tau})\,\bm{\tau}$

where $\bm{\tau}$ is the unit tangent vector to the crack front in the
direction of increasing arc length (Figure 2.1). As in the FEM method, it is
easy to compute $\bm{\tau}$ from the partial derivative of (2) with respect to
the arc length $s$ of $c_{0}$. If we assume that the toughness of the material
is isotropic, the critical energy release rate $G_{c}$ being a parameter
measuring it, we introduce the convex crack stability domain

$K=\left\\{\bm{G}\;\;\mbox{such that}\;\;\parallel\bm{G}_{\perp}\parallel\leq
G_{c}\right\\}$

and the crack extension rule

$\dot{\bm{\psi}}\in\partial\varphi(\bm{G})=\partial\chi_{K}(\bm{G})$

that is the normality law

if |  $\bm{G}\in\mathring{K}$ then |
---|---|---
| $\dot{\bm{\psi}}=\bm{0}$ | ! crack stability
else |  $\bm{G}\in\partial K$ and $\exists\lambda\geq 0,\qquad\dot{\bm{\psi}}=\lambda\,\bm{G}_{\perp}$ | ! crack extension

and $\dot{\bm{\psi}}$ is perpendicular to $\bm{\tau}$ when the crack extends.
The Fenchel polar is the support function

$\varphi^{*}(\dot{\bm{\psi}})=\int_{c_{t}}G_{c}\parallel\dot{\bm{\psi}}\parallel\,\mbox{d}s$

where $s$ is the arc length of $c_{t}$. With this choice of constitutive law,
we particularize the variational principle in the form

Symplectic BEN principle 3. The natural evolution of the system minimizes

$\Pi(\bm{\xi})=\int_{0}^{T}\left\\{\int_{c_{t}}G_{c}\parallel\dot{\bm{\psi}}\parallel\,\mbox{d}s-\left\langle\dot{\bm{\psi}},\bm{G}\right\rangle\right\\}\mbox{
d}t$ (8)

among the admissible curves $t\mapsto\bm{\xi}(t)$ such that
$\nabla\cdot\bm{\sigma}+\bm{f}=\rho\,\ddot{\bm{u}}$,
$\parallel\bm{G}_{\perp}\parallel\leq G_{c}$ and the minimum is zero.

##  6 Link with the classical variational approach to fracture

First we recall that the last term in (8) is in fact the last one of the
integrand of (4). Owing to (6), we have

$-\omega(\dot{\bm{z}}_{I},\dot{\bm{z}})=\omega(\dot{\bm{z}}_{R},\dot{\bm{z}})=\langle
D_{\bm{\eta}}H,\dot{\bm{\eta}}\rangle-\langle\dot{\bm{\xi}},-D_{\bm{\xi}}H\rangle=\dot{H}-D_{t}H$

Besides the crack surface element comprised between the crack fronts $c_{t}$
and $c_{t+dt}$, and the streamlines of $\bm{x}_{0}$ and
$\bm{x}_{0}+d\bm{x}_{0}$ is an infinitesimal parallelogram with two sides
supported by the vectors $d_{1}\bm{x}=\bm{\tau}\,ds$ and
$d_{2}\bm{x}=\dot{\bm{\psi}}\,dt$. Its area is

$d\Gamma=\parallel\bm{\tau}\times\dot{\bm{\psi}}\parallel\,ds\,dt=\parallel\bm{\tau}\parallel\,\parallel\dot{\bm{\psi}}\parallel\,\mid\sin{(\bm{\tau},\dot{\bm{\psi}})}\mid\,ds\,dt=\parallel\dot{\bm{\psi}}\parallel\,ds\,dt$

because $\bm{\tau}$ is a unit vector and $\dot{\bm{\psi}}$ is perpendicular to
$\bm{\tau}$.

If we assume that the toughness properties of the material are homogeneous,
the critical energy release rate is uniform

$G_{c}(\bm{x})=G_{c}$

The area of a surface will be denoted by $\mid\bullet\mid$. The variational
principle can be recast as

Symplectic BEN principle 4. The natural evolution of the system minimizes

$\Pi(\Gamma,\bm{u})=G_{c}\,\mid\Gamma_{T}\backslash\Gamma_{0}\mid-\int_{0}^{T}D_{t}H\,\mbox{
d}t+H_{T}-H_{0}$ (9)

among the admissible curves $t\mapsto(\Gamma(t),\bm{u}(t))$ such that
$\nabla\cdot\bm{\sigma}+\bm{f}=\rho\,\ddot{\bm{u}}$,
$\parallel\bm{G}_{\perp}\parallel\leq G_{c}$ and the minimum is zero.

Two cases must be considered in the applications:

  * •

If the loading is controlled by the forces, the structural evolution is in
general dynamic and the functional to minimize is

$\Pi(\Gamma,\bm{u})=G_{c}\,\mid\Gamma_{T}\backslash\Gamma_{0}\mid-\int_{0}^{T}D_{t}H\,\mbox{
d}t+H_{T}-H_{0}$

  * •

If the loading is slow and controlled by the displacements (for instance,
$\bar{\bm{f}}$ is null, the support is divided into two disjoint parts, one of
them is fixed to the foundation, the other one is loaded), the evolution is
quasi-static and the functional to minimize is reduced to

$\Pi(\Gamma,\bm{u})=G_{c}\,\mid\Gamma_{T}\backslash\Gamma_{0}\mid+H_{T}-H_{0}$

Regardless of the constant $H_{0}$, we recover the Munford-Shah functional
[62, 36, 14, 12]

##  7 Incremental method

Our goal now is to find an incremental or step-by-step formulation as a by-
product of our general principle. As usual, the time interval
$\left[0,T\right]$ is divided into $N$ small sub-intervals
$\left[t_{n},t_{n+1}\right]$ with $t_{0}=0$ and $t_{N}=T$. the value of a
variable $\bm{a}$ at time $t_{n}$ is denoted $\bm{a}_{n}$. The time step is
$\Delta t=t_{n+1}-t_{n}$ and $\Delta\bm{a}=\bm{a}_{n+1}-\bm{a}_{n}$ is the
increment of $\bm{a}$.

$\bm{\xi}_{n}=(\bm{u}_{n},\bm{\psi}_{n})$ being known from the initial
conditions for $n=0$ or from the previous increment otherwise and the time
step being prescribed, the problem is to calculate $\bm{\xi}_{n+1}$ or
equivalently $\Delta\bm{\xi}$. The idea is to apply the symplectic BEN
principle 3 to the sub- interval $\left[t_{n},t_{n+1}\right]$. The time step
being small, the functional can be approximated by considering the integrand
is constant. As for the radial return algorithm in plasticity, we use the
implicit scheme by evaluating the crack extension rate $\dot{\bm{\psi}}$ and
the driving force $\bm{G}$ at the end of the step. The crack extension
increment is

$\Delta\bm{\psi}=\dot{\bm{\psi}}_{n+1}\,\Delta t=\lambda\,\bm{G}_{n+1}\Delta
t$ (10)

Hence the incremental principle reads

Symplectic BEN principle $3^{*}$. The natural evolution of the system
minimizes

$\Delta\Pi(\bm{\xi}_{n+1})=\int_{c_{t_{n+1}}}G_{c}\parallel\Delta\bm{\psi}\parallel\,\mbox{d}s-\left\langle\Delta\bm{\psi},\bm{G}_{n+1}\right\rangle$
(11)

among the $\bm{\xi}_{n+1}$ such that
$\nabla\cdot\bm{\sigma}_{n+1}+\bm{f}_{n+1}=\rho\,\ddot{\bm{u}}_{n+1}$,
$\parallel(\bm{G}_{n+1})_{\perp}\parallel\leq G_{c}$ and the minimum is zero.

Likewise, the last version of the principle in the previous section gives rise
to the incremental form

Symplectic BEN principle $4^{*}$. The natural evolution of the system
minimizes

$\Delta\Pi(\Gamma_{n+1},\bm{u}_{n+1})=G_{c}\mid\Delta\Gamma\mid-(D_{t}H)_{n+1}\,\Delta
t+\Delta H$ (12)

among the couples $(\Gamma_{n+1},\bm{u}_{n+1})$ such that
$\nabla\cdot\bm{\sigma}_{n+1}+\bm{f}_{n+1}=\rho\,\ddot{\bm{u}}_{n+1}$,
$\parallel(\bm{G}_{n+1})_{\perp}\parallel\leq G_{c}$ and the minimum is zero.

If the loading is slow and controlled by the displacements, the evolution is
quasi-static and the functional to minimize is reduced to

$\Delta\Pi(\Gamma_{n+1},\bm{u}_{n+1})=G_{c}\mid\Delta\Gamma\mid+\Delta H$

and we recover Buliga’s incremental formulation [14].

##  8 Comparison of the constitutive law with experiments

Figure 8.1: Kinked crack

The aim of this section is to check the validity of the normality law with
respect to the experimental results on PMMA specimens [65]. A good test is to
considered an initial straight crack in a plate which, under mixed mode
loading, may extend suddenly in a direction deviating of an angle $\beta$ from
the original one by kinking (Figure 8.1). The length of the kink crack is
taken very small. The idea is to use the incremental formulation of the
previous section with only one step. To simplify the notations, the initial
value $a_{0}$ of a given quantity is denoted simply $a$ while the value
$a_{1}$ at the end of the step is denoted $a^{*}$. In Hellen et al. [43],
Euler explicit scheme is used that can read with the simplified notations

$\Delta\bm{\psi}=\dot{\bm{\psi}}\,\Delta t=\lambda\,\bm{G}\Delta t$

In the frame of origin $O$ at the initial crack tip, the axis $Ox_{1}$ in the
direction ahead the initial crack and the axis $Ox_{2}$ perpendicular to
$Ox_{1}$ within the plane (Figure 8.1), the driven force is given in terms of
stress intensity factors (SIFs) by

$G_{1}=\frac{1}{\bar{E}}\,(K^{2}_{I}+K^{2}_{II}),\qquad
G_{2}=-\frac{1}{\bar{E}}\,2\,K_{I}K_{II}$

where $\kappa=3-4\,\nu$ in plane strain, $\kappa=(3-\nu)/(1+\nu)$ in plane
stress and

$\frac{1}{\bar{E}}=\frac{(1+\nu)\,(1+\kappa)}{4\,E}$

The deviation angle $\beta$ of the kink crack is given by the slope of the
driven force

$\beta=\tan^{-1}\left(-\frac{2\,K_{I}K_{II}}{K^{2}_{I}+K^{2}_{II}}\right)$

In mode $I$, the formula gives $\beta=0^{\circ}$, according to the experience.
For $K_{I}=K_{II}$, $\beta=-45^{\circ}$ but the experimental values are above
$-49^{\circ}$. In mode $II$, $\beta=0^{\circ}$ while the experimental values
are in the range from $-70^{\circ}$ to $-74^{\circ}$. Clearly, the predictions
are disastrous. Our opinion is that the constitutive law must not be a priori
rejected but it is the explicit scheme which is problematic.

The problem of the kinked crack were studied by many authors. It is an awkward
problem of Elasticity and only approximated formula were proposed to express
the SIFs at the kink crack tip in terms of the SIFs at the initial crack tip.
At the limit of vanishing length of the kink crack, we adopt in the sequel the
expression proposed in [23]

$K^{*}_{I}(\beta)=\cos^{3}\frac{\beta}{2}\,K_{I}-3\,\sin\frac{\beta}{2}\,\cos^{2}\frac{\beta}{2}\,K_{II}$
$K^{*}_{II}(\beta)=\sin\frac{\beta}{2}\,\cos^{2}\frac{\beta}{2}\,K_{I}+\cos\frac{\beta}{2}\,\left(1-3\,\sin^{2}\frac{\beta}{2}\right)\,K_{II}$

However, it is noticeable that a correction is proposed in [34] for the
largest values of the $\mid\beta\mid$ and the dependency with respect to the
length of the kink crack is taken into account in [49, 2] but these
improvements will be not considered here. With the simplified notations, the
implicit scheme (10) reads

$\Delta\bm{\psi}=\dot{\bm{\psi}}^{*}\,\Delta t=\lambda\,\bm{G}^{*}\Delta t$
(13)

In the frame of origin $K$ at the kink crack tip, the axis $Kx^{*}_{1}$ in the
direction ahead the kink crack and the axis $Kx^{*}_{2}$ perpendicular to
$Kx^{*}_{1}$ within the plane (Figure 8.1), the driven force is given in terms
of stress intensity factors (SIFs) by

$G^{*}_{1}=\frac{1}{\bar{E}}\,\left[(K^{*}_{I})^{2}+(K^{*}_{II})^{2}\right],\qquad
G^{*}_{2}=-\frac{1}{\bar{E}}\,2\,K^{*}_{I}K^{*}_{II}$

The normality law (13) entails

$\bar{E}\,G^{*}_{2}=-2\,K^{*}_{I}K^{*}_{II}=0,\qquad\bar{E}\,\parallel\bm{G}^{*}\parallel=(K^{*}_{I})^{2}+(K^{*}_{II})^{2},\qquad$
(14)

To satisfy the former condition, two scenarios may be considered for the crack
extension:

  1. 1.

Scenario 1. The kink crack inclination is $\beta_{1}$, solution of
$K^{*}_{II}(\beta_{1})=0$, and
$\bar{E}\,\parallel\bm{G}^{*}\parallel=(K^{*}_{I}(\beta_{1}))^{2}$

  2. 2.

Scenario 2. The kink crack inclination is $\beta_{2}$, solution of
$K^{*}_{I}(\beta_{2})=0$, and
$\bar{E}\,\parallel\bm{G}^{*}\parallel=(K^{*}_{II}(\beta_{2}))^{2}$

According to the stability criterion, the crack extends for the inclination
$\beta_{i}$ with the maximum driven force magnitude. Then scenario 1 is
realized if $(K^{*}_{I}(\beta_{1}))^{2}>(K^{*}_{II}(\beta_{2}))^{2}$ and
scenario 2 otherwise.

The admissible solutions must be searched in the interval
$\left]-\pi,\pi\right[$. The mode $I$ singularity disappearing in compression,
$K_{I}$ is positive. In mode $II$, the SIF $K_{II}$ may be positive or
negative, depending on the sign of the shear loading. Before tackling the
mixed mode, let us consider the limit case:

  * •

Mode I. For scenario 1, there is a unique admissible solution
$\beta_{1}=0^{\circ}$ and $(K^{*}_{I}(\beta_{1}))^{2}=K^{2}_{I}$. For scenario
2, there is no admissible solution. Then the scenario 1 occurs, according to
the experimental observations.

Let us examine now the general case ($K_{I}\,K_{II}\neq 0$). For scenario 1,
there are two admissible solutions

$\beta_{1\pm}=\sin^{-1}\left(\frac{\frac{K_{II}}{K_{I}}\pm
3\,\frac{K_{II}}{K_{I}}\,\sqrt{8\,\left(\frac{K_{II}}{K_{I}}\right)^{2}+1}}{9\,\left(\frac{K_{II}}{K_{I}}\right)^{2}+1}\right)$

For scenario 2, there is a unique admissible solution

$\beta_{2}=2\,\tan^{-1}\left(\frac{K_{I}}{3\,K_{II}}\right)$

In particular, let us discuss the case:

  * •

Mixed mode $\bm{K}_{\bm{I}}=\bm{K}_{\bm{II}}>0$. For scenario 1, there are two
admissible solutions: $\beta_{1-}=\sin^{-1}(-4/5)=-53.13^{\circ}$ then
$(K^{*}_{I}(\beta_{1-}))^{2}=(1.847\,K_{I})^{2}$ and
$\beta_{1+}=\sin^{-1}(1)=90^{\circ}$ then
$(K^{*}_{I}(\beta_{1+}))^{2}=(-1.14\,K_{I})^{2}$. For scenario 2, there is a
unique admissible solution $\beta_{2}=36.87^{\circ}$ and
$(K^{*}_{II}(\beta_{2}))^{2}=(0.948\,K_{I})^{2}$. Comparing the values of
$(K^{*}_{I}(\beta_{1-}))^{2}$, $(K^{*}_{I}(\beta_{1+}))^{2}$ and
$(K^{*}_{II}(\beta_{2}))^{2}$, we conclude that the scenario 1 occurs with the
angle $\beta=\beta_{1-}=-53.13^{\circ}$ closed to the experimental values
($-49^{\circ}$ for $K_{II}=0.874\,K_{I}$ and values in the range from
$-56^{\circ}$ to $-59^{\circ}$ for $K_{II}=1.142\,K_{I}$).

In this respect, it is worth to remark that, if the sign of $K_{II}$ is
reversed, the sign of the kink angles are reversed too :
$\beta_{1-}=53.13^{\circ}$, $\beta_{1+}=-90^{\circ}$ and
$\beta_{2}=-36.87^{\circ}$ but the corresponding values of
$(K^{*}_{I}(\beta_{1-}))^{2}$, $(K^{*}_{I}(\beta_{1+}))^{2}$ and
$(K^{*}_{II}(\beta_{2}))^{2}$ are the same. We conclude that the sign of the
kink angle is reversed: $\beta=53.13^{\circ}$.

Finally, let us deal with the other limit case:

  * •

Mode II. For scenario 1, there are two admissible solutions
$\beta_{1\pm}=2\,\sin^{-1}(\pm\,1/\sqrt{3})=\pm\,70.47^{\circ}$ and
$(K^{*}_{I}(\beta_{1\pm}))^{2}=4\,K^{2}_{II}/3$. For scenario 2, there is a
unique admissible solution $\beta_{2}=0^{\circ}$ and
$(K^{*}_{II}(\beta_{2}))^{2}=K^{2}_{I}$. Then the scenario 1 occurs and the
sign of the kink angle depends of the one of mode $II$. For $K_{II}>0$, the
prediction $\beta=-70.47^{\circ}$ is just above the experimental values.

Comment 1. On the basis of the experimental data, Richard proposed an
empirical criterion in terms of the SIFs at the original crack tip

$\frac{K_{I}}{K_{Ic}}+\left(\frac{K_{II}}{K_{IIc}}\right)^{2}=1$

According to (14), the crack extends when

$(K^{*}_{I})^{2}+(K^{*}_{II})^{2}=\bar{E}\,G_{c}$

that gives, for the scenario 1 of mode I, $K^{2}_{I}=\bar{E}\,G_{c}$ then

$K^{2}_{Ic}=\bar{E}\,G_{c}$

and, for the scenario 1 of mode II,
$4\,K^{2}_{II}/3=\bar{E}\,G_{c}=K^{2}_{Ic}$ then

$\frac{K_{IIc}}{K_{Ic}}=\frac{\sqrt{3}}{2}=0.86$

a value better than $\sqrt{2/3}=0.81$ proposed in [34], by comparison with the
experimental data covering the range from $0.88$ to $0.95$. Moreover, our
prediction for the mixed mode $K_{I}=K_{II}$ leads to values
$\frac{K_{I}}{K_{Ic}}=0.54$ and $\frac{K_{II}}{K_{IIc}}=0.62$ closed to the
curve of Richard criterion.

Comment 2. At least in this example, the relevant scenario corresponds always
to the Principle of Local Symmetry $K^{*}_{II}(\beta)=0$ proposed in 1974 by
Goldstein and Salganik [38]. This principle is valid at least for a kinked
straight crack with vanishing kink crack but there can be no assurance that it
is valid for arbitrary cracks and we recommend to replace it by the normality
law proposed in this paper.

Comment 3. The present criterion should not be confused with the Maximum
Energy Release Rate criterion, proposed in 1974 by Hussain et al. [45], for
which the maximum must be searched among all the directions of the crack
extension while with our criterion the maximum is found only among the
directions satisfying the normality law, according to the different scenarios.
As argued in [20], the questions of when and how a crack propagates should be
simultaneously investigated and the energy conservation is not sufficient for
such a task.

Comment 4. This example shows also that the implicit scheme must be preferred
to the explicit scheme in fracture mechanics.

##  9 Calculation of the crack driving force

As we are working at time $t$ in this Section, the dependency with respect to
the time will be omitted. In order to obtain the expression of $\bm{G}$, we
use a special form of the calculus of variation performed on the jet space of
order one ([1], [31], [52]). For more details about the jet spaces, the reader
is referred for instance to [67]. The first jet prolongation of the smooth
function $\bm{u}:\Omega_{t}\rightarrow\mathbb{R}^{3}$ is the function
$j^{1}\bm{u}$ from $\Omega_{t}$ into the jet bundle
$J^{1}(\Omega_{t},\mathbb{R}^{3})$ such that

$j^{1}\bm{u}\ (\bm{x})=\left(\bm{x},\bm{u}\
(\bm{x}),\nabla\bm{u}(\bm{x})\right)$

The Hamiltonian (3) at time $t$ has the form

$H=\int_{\Omega_{t}}h(\bm{x},\bm{u},\nabla\bm{u},\bm{p})\mbox{
d}^{3}x-\int_{\partial\Omega_{1}}\bar{\bm{f}}\cdot\bm{u}\mbox{
d}(\partial\Omega_{1})$ $H=\int_{\Omega_{t}}h(j^{1}\bm{u},\bm{p})\mbox{
d}^{3}x-\int_{\partial\Omega_{1}}\bar{\bm{f}}\cdot\bm{u}\mbox{
d}(\partial\Omega_{1})$

Figure 9.1: Cylindrical neighbourhood

The new viewpoint which consists in replacing the original field $\bm{u}$ by
its first jet prolongation $j^{1}\bm{u}$ leads to perform variations not only
on the field and its derivatives but also on the variable $\bm{x}$. We want to
calculate the functional derivative of $H$ in the direction $\delta\bm{\psi}$
defined on the crack front $c_{t}$. The idea is to prolong it locally by a
field $\delta\bm{x}$ on an open cylindrical neighbourhood $C_{\eta}$ of radius
$\eta$ around the crack front, the field vanishing at radius $\eta$ and beyond
(Figure 9.1). To avoid the singularity of the stress field, we replace
$C_{\eta}$ by $C_{\eta,r}=C_{\eta}\backslash\bar{C}_{r}$ where $C_{r}$ is an
open cylindrical neighbourhood of radius $r<\eta$ around the crack front.
Next, we pass to the limit $\eta\rightarrow r\rightarrow 0$.

We consider now a new parameterization given by a regular map
$\bm{x}=\bm{\theta}(\bm{y})$ of class $C^{1}$ on
$C^{\prime}_{\eta,r}=\bm{\theta}^{-1}(C_{\eta,r})$ and we perform the
variation of the function $\bm{\theta}$, the new variable being $\bm{y}$. The
variation of the action is

$\delta
H=\delta\int_{C^{\prime}_{\eta,r}}h(\bm{\theta}(\bm{y}),\bm{u},\nabla_{\bm{y}}\bm{u}\cdot\nabla_{\bm{x}}\bm{y},(\nabla_{\bm{x}}\bm{y})^{T}\cdot\bm{p}^{\prime})\,\det(\nabla_{\bm{y}}\bm{x})\mbox{
d}^{3}y$ (15)

where $\bm{p}^{\prime}$ is the linear momentum in $\bm{y}$ coordinates. We do
not give here the details of the calculations that can be found in Appendix A
but only indicate the sketch of the method and the result. As usual in
calculus of variation, we integrate by part. Next we consider the limit case
where the function $\bm{\theta}$ approaches the identity of $\Omega_{t}$,
$\bm{y}$ approaches $\bm{x}$, $\bm{p}^{\prime}$ approaches $\bm{p}$ and
$C^{\prime}_{\eta,r}$ approaches $C_{\eta,r}$. Finally, the variation of the
Hamiltonian has the form

$\delta H=\int_{\partial C_{\eta,r}}\bm{n}\cdot\bm{T}\cdot\delta\bm{x}\mbox{
d}(\partial
C_{\eta,\rho})-\int_{C_{\eta,r}}\left(\nabla\cdot\bm{T}+\bm{f}\right)\cdot\delta\bm{x}\mbox{
d}^{3}x$ (16)

where

$\bm{T}=h\,\bm{1}_{\mathbb{R}^{3}}-\bm{\sigma}\cdot\nabla\bm{u}-\dot{\bm{u}}\otimes\bm{p}$
(17)

The contribution on the cylindrical part $\partial C_{\eta,r}$ of radius
$\eta$ vanishes as $\delta\bm{x}$. When $\eta$ approaches $r$, the
contributions on the remaining parts of $\partial C_{\eta,r}$ excepted the
cylindrical part of radius $r$, vanish too. Let $\bm{n}^{\prime}$ the outward
unit vector to $\partial C_{r}$, opposite to the outward unit vector $\bm{n}$
to $\partial C_{\eta,r}$ (Figure 9.1)

$\bm{n}^{\prime}=-\bm{n}$

Then using the polar coordinates $r,\theta$, one has

$\delta H=-\lim_{r\rightarrow 0}\int_{\partial C_{r}\cap\partial
C_{\eta,r}}\bm{n}^{\prime}\cdot\bm{T}\cdot\delta\bm{x}\;r\mbox{ d}\theta\mbox{
d}s$

Finally, recalling that the functional derivative is a density on the crack
front, the value of the driving force corresponding to the arc length $s$ is

$\bm{G}=-D_{\delta\bm{x}}H=\lim_{r\rightarrow
0}\int_{c_{r,s}}\bm{n}^{\prime}\cdot\bm{T}\;r\mbox{ d}\theta$ (18)

where $c_{r,s}$ is the line at the intersection between the cylindrical-shaped
surface $\partial C_{r}\cap\partial C_{\eta,r}$ and the plane normal to the
crack front $c_{t}$ at the arc length $s$. This limit is called by Gurtin et
al. [40] a tip integral and denoted

$\bm{G}=\oint_{\mbox{\tiny{tip}}}\bm{n}^{\prime}\cdot\bm{T}$

but we do not recover the last term of (17) in Gurtin’s and Stolz’s [70]
expressions of $\bm{T}$. In this form, the driving force appears as a
generalization of Rice-Eshelby $J$ integral. In statics, this quantity has the
well known property of path-independence and it is relevant to wonder whether
this property remains true in dynamics. The elastic stress field being
singular at the crack front, the volume force field may be neglected when $r$
approaches zero. In Appendix B, it is shown that for the natural evolution of
the system

$\nabla\cdot\bm{T}=(\nabla\bm{u})\cdot\dot{\bm{p}}-(\nabla\cdot\dot{\bm{u}})\,\bm{p}-2\,(\nabla_{s}\bm{p})\cdot\dot{\bm{u}}$
(19)

Unfortunately, this expression does not vanish in general but only in the
limit case of statics where we recover the classical Rice-Eshelby identity

$\nabla\cdot\bm{T}=\nabla\cdot\left[w\,\bm{1}_{\mathbb{R}^{3}}-\bm{\sigma}\cdot\nabla\bm{u}\right]=\bm{0}$

Then except the particular case of statics, there is no path-independence of
the integral in dynamics and passing to the limit is required in (18).

##  10 Conclusions and perspectives

In this work, we presented a particular version of the non incremental
symplectic BEN principle. It can be generalized in various ways:

  * •

The Hamiltonian $H$ depends on $\bm{u},\bm{p}$ and $\bm{\psi}$ (through
$\Omega_{t}=\Omega\backslash\Gamma_{t}$) but not on $\bm{\pi}$. It was a bias
motivated by the wish to recover classical formulations of the fracture
mechanics but we could give up this restrictive assumption by introducing an
explicit dependence of the Hamiltonian on $\bm{\pi}$.

  * •

For the same motivation, the dissipation potential $\phi$ depends on
$\dot{\bm{\pi}}$ but not on $\dot{\bm{u}},\dot{\bm{p}}$ and $\dot{\bm{\psi}}$.
Nevertheless there is nothing preventing us from lifting this restriction, in
particular concerning $\dot{\bm{\psi}}$.

  * •

The potential $\varphi$ could fully depend on $\bm{G}$, not only on
$\bm{G}_{\perp}$. It could be also generalized to anisotropic behaviours with
symmetry group techniques.

  * •

For easiness, we considered only 1-homogeneous potential $\varphi^{*}$ in this
work to illustrate the method but we can add extra terms to represent for
instance the dependence of the critical energy release rate with respect to
the crack extension vector as in [22, 46, 47, 48]. The crack extension by
fatigue can be also modelled thanks to homogeneous potentials of degree linked
to the crack growth law ([51], [12]).

  * •

We could easily take into account more sources of dissipation such as damping,
plasticity, damage and so on, by considering in the symplectic formalism new
internal variables and the corresponding momenta, and by simply adding terms
in the dissipation potential $\phi$.

  * •

Sometimes there is no getting away from experimental facts and the kinetic law
is non associated. For such atypical events, the present principle can be
generalized, replacing the sum of the dissipation potential and its Fenchel
polar by a bipotential.

We are convinced of the interest of this non incremental principle in
computational structural mechanics because the error can be controlled
uniformly on the whole evolution, in contrast to incremental methods which
accumulate the errors and degrade the accuracy over time.

We hope that our application of the normality law to the kinked crack allows
to take forward the discussion in the literature on the choice of the crack
stability criterion.

In the future, we plan to develop numerical methods based on the present
approach. We already have given a hint about the discretization of the crack
flow at the end of Section 2.2. In order to avoid a cumbersome remeshing, an
efficient technique is the XFEM where the standard displacement-based
approximation is enriched by incorporating both discontinuous fields and the
near tip asymptotic fields ([59], [28], [29], [60]).

Acknowledgements

The author would like to thank Marius Buliga, Djimedo Kondo, Abdelbacet
Oueslati, Pierre Gosselet, Céline Bouby and Long Cheng for their questions and
comments during seminars that allowed to improve the paper.

Appendix A

To calculate the expression of the crack driven force, we start from the
variation of the Hamiltonian (15)

$\delta
H=\delta\int_{C^{\prime}_{\eta,r}}h(\bm{\theta}(\bm{y}),\bm{u},\nabla_{\bm{y}}\bm{u}\cdot\nabla_{\bm{x}}\bm{y},(\nabla_{\bm{x}}\bm{y})^{T}\cdot\bm{p}^{\prime})\,\det(\nabla_{\bm{y}}\bm{x})\mbox{
d}^{3}y$

and we take the variation with respect to the fields $\bm{u}$ and
$\bm{\theta}$

$\delta
H=\int_{C^{\prime}_{\eta,r}}\\{\left[Tr\left(\bm{\sigma}\,\delta(\nabla_{\bm{y}}\bm{u}\cdot\nabla_{\bm{x}}\bm{y})\right)+\dot{\bm{u}}\cdot\delta((\nabla_{\bm{x}}\bm{y})^{T}\cdot\bm{p}^{\prime})-\bm{f}\cdot\delta\bm{x}\right]\,\det(\nabla_{\bm{y}}\bm{x})$
$\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+h\,\delta(\det(\nabla_{\bm{y}}\bm{x}))\\}\,\mbox{
d}^{3}y$ (20)

Next we calculate the variation of the field derivative in terms of the
derivative of the variation of the field

$\delta(\nabla_{\bm{y}}\bm{u}\cdot\nabla_{\bm{x}}\bm{y})=-\nabla_{\bm{y}}\bm{u}\cdot\nabla_{\bm{x}}\bm{y}\cdot\nabla_{\bm{y}}(\delta\bm{x})\cdot\nabla_{\bm{x}}\bm{y}$
$\delta(\det(\nabla_{\bm{y}}\bm{x}))=Tr\left(\nabla_{\bm{y}}(\delta\bm{x})\cdot\mbox{adj}(\nabla_{\bm{y}}\bm{x})\right)$

where $adj(\bm{A})$ is the transposed of the matrix of cofactors of $\bm{A}$

$\mbox{adj}(\bm{A})=\det(\bm{A})\,\bm{A}^{-1}$

Likewise, one has

$\delta((\nabla_{\bm{x}}\bm{y})^{T}\cdot\bm{p}^{\prime})=-(\nabla_{\bm{x}}\bm{y}\cdot\nabla_{\bm{y}}(\delta\bm{x})\cdot\nabla_{\bm{x}}\bm{y})^{T}\cdot\bm{p}^{\prime}=-(\nabla_{\bm{x}}\bm{y})^{T}\cdot(\nabla_{\bm{y}}(\delta\bm{x}))^{T}\cdot(\nabla_{\bm{x}}\bm{y})^{T}\cdot\bm{p}^{\prime}$

Inserting the three previous expression into (20) and after simple
manipulations, we obtain

$\delta
H=\int_{C^{\prime}_{\eta,r}}\\{Tr\left(\mbox{adj}(\nabla_{\bm{y}}\bm{x})\cdot\left(h\,\bm{1}_{\mathbb{R}^{3}}-\bm{\sigma}\cdot\nabla_{\bm{y}}\bm{u}\cdot\nabla_{\bm{x}}\bm{y}-(\dot{\bm{u}}\otimes\bm{p}^{\prime})\cdot\nabla_{\bm{x}}\bm{y}\right)\cdot\nabla_{\bm{y}}(\delta\bm{x})\right)$
$\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad-\det(\nabla_{\bm{y}}\bm{x})\,\bm{f}\cdot\delta\bm{x}\\}\,\mbox{
d}^{3}y$

Integrating by parts leads to

$\delta H=\int_{\partial
C^{\prime}_{\eta,r}}\bm{n}\cdot\mbox{adj}(\nabla_{\bm{y}}\bm{x})\cdot\left(h\,\bm{1}_{\mathbb{R}^{3}}-\bm{\sigma}\cdot\nabla_{\bm{y}}\bm{u}\cdot\nabla_{\bm{x}}\bm{y}-(\dot{\bm{u}}\otimes\bm{p}^{\prime})\cdot\nabla_{\bm{x}}\bm{y}\right)\cdot\delta\bm{x}\,\mbox{d}(\partial
C^{\prime}_{\eta,r})$
$-\int_{C^{\prime}_{\eta,r}}\\{\nabla\cdot\left(\mbox{adj}(\nabla_{\bm{y}}\bm{x})\cdot\left(h\,\bm{1}_{\mathbb{R}^{3}}-\bm{\sigma}\cdot\nabla_{\bm{y}}\bm{u}\cdot\nabla_{\bm{x}}\bm{y}-(\dot{\bm{u}}\otimes\bm{p}^{\prime})\cdot\nabla_{\bm{x}}\bm{y}\right)\right)$
$\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\det(\nabla_{\bm{y}}\bm{x})\,\bm{f}\\}\cdot\delta\bm{x}\,\mbox{
d}^{3}y$

Considering the limit case where the function $\bm{\theta}$ approaches the
identity of $C_{\eta,r}$, $\bm{y}$ approaches $\bm{x}$, $\bm{p}^{\prime}$
approaches $\bm{p}$ and $C^{\prime}_{\eta,r}$ approaches $C_{\eta,r}$, we
obtain (16)

$\delta H=\int_{\partial C_{\eta,r}}\bm{n}\cdot\bm{T}\cdot\delta\bm{x}\mbox{
d}(\partial
C_{\eta,\rho})-\int_{C_{\eta,r}}\left(\nabla\cdot\bm{T}+\bm{f}\right)\cdot\delta\bm{x}\mbox{
d}^{3}x$

where occurs the tensor (17)

$\bm{T}=h\,\bm{1}_{\mathbb{R}^{3}}-\bm{\sigma}\cdot\nabla\bm{u}-\dot{\bm{u}}\otimes\bm{p}$

Appendix B

Because of the stress singularity, the volume force that are regular may be
consider as uniform near the crack front. In one hand, the variation of the
Hamitonian density $h$ in the direction $\delta\bm{x}$ is

$\delta h=\nabla
h\cdot\delta\bm{x}=\nabla\cdot(h\,\bm{1}_{\mathbb{R}^{3}})\cdot\delta\bm{x}$
(21)

In the other hand, according to (3), $h$ is a function of $\bm{x}$ through
$j^{1}\bm{u}$ and $\bm{p}$

$h=\dfrac{1}{2\rho}\parallel\bm{p}\parallel^{2}+w(\nabla_{s}\bm{u})-\bm{f}(t)\cdot\bm{u}$

Using the chain rule and assuming the volume force are uniform, one has

$\delta
h=Tr(\bm{\sigma}\cdot\nabla(\delta\bm{u}))+\dot{\bm{u}}\cdot\delta\bm{p}-\bm{f}\cdot\bm{u}$

Using the linear momentum balance (intermediate constraints in (5) which is
satisfied for the natural evolution of the system), the first term of the
right hand member becomes

$Tr(\bm{\sigma}\cdot\nabla(\delta\bm{u}))=\nabla\cdot(\bm{\sigma}\cdot\delta\bm{u})-(\nabla\cdot\bm{\sigma})\cdot\delta\bm{u}=\nabla\cdot(\bm{\sigma}\cdot\delta\bm{u})+(\bm{f}-\dot{\bm{p}})\cdot\delta\bm{u}$

Hence

$\delta
h=\nabla\cdot(\bm{\sigma}\cdot\delta\bm{u})-\dot{\bm{p}}\cdot\delta\bm{u}+\dot{\bm{u}}\cdot\delta\bm{p}$

The infinitesimal variations of the displacement and linear momentum fields
resulting from an infinitesimal arbitrary variation $\delta\bm{x}$ being

$\delta\bm{u}=\nabla\bm{u}\cdot\delta\bm{x},\qquad\delta\bm{p}=\nabla\bm{p}\cdot\delta\bm{x}$

and in particular if the field $\delta\bm{x}$ is uniform

$\delta
h=\left[\nabla\cdot(\bm{\sigma}\cdot\nabla\bm{u})-\dot{\bm{p}}\cdot\nabla\bm{u}+\dot{\bm{u}}\cdot\nabla\bm{p}\right]\cdot\delta\bm{x}$
(22)

The expressions (21) and (22) being equal for every $\delta\bm{x}$, we obtain

$\nabla\cdot(h\,\bm{1}_{\mathbb{R}^{3}}-\bm{\sigma}\cdot\nabla\bm{u})=(\nabla\bm{u})^{T}\cdot\dot{\bm{p}}-(\nabla\bm{p})^{T}\cdot\dot{\bm{u}}$

Taking into account the identity

$\nabla\cdot(\dot{\bm{u}}\otimes\bm{p})=(\nabla\cdot\dot{\bm{u}})\,\bm{p}+(\nabla\bm{p})\cdot\dot{\bm{u}}$

we obtain (19)

$\nabla\cdot\bm{T}=(\nabla\bm{u})\cdot\dot{\bm{p}}-(\nabla\cdot\dot{\bm{u}})\,\bm{p}-2\,(\nabla_{s}\bm{p})\cdot\dot{\bm{u}}$

where occurs the tensor (17)

$\bm{T}=h\,\bm{1}_{\mathbb{R}^{3}}-\bm{\sigma}\cdot\nabla\bm{u}-\dot{\bm{u}}\otimes\bm{p}$

## References

  * [1] Aldaya, V., Azcárraga, J. A., 1978. Variational principles on rth order jets of fibre bundles in field theory, J. Math. Phys. 19(9), 1869-1975.
  * [2] Amestoy, M., Leblond, J.B., 1992. Crack path in plane situation – II. Detailed form of the expansion of the stress intensity factors. Int J Solids Struct. 29, 465–501.
  * [3] Ambrosio, L., Tortorelli, V., 1992. On the approximation of free discontinuity problems. Bollett UMI. 7(6B), 105–123.
  * [4] Aubin, J.P., Cellina, A.,Nohel, J., 1997. Monotone trajectories of multivalued dynamical systems. Annali di Matematica Pura ed Appl. 115, 99-117.
  * [5] Aubin, J.P., 2002. Boundary-Value Problems for Systems of Hamilton-Jacobi-Bellman Inclusions with Constraints. SIAM J. Control. 41, 425-456.
  * [6] Barenblatt, G.I., 1959. The formation of equilibrium cracks during brittle fracture: general ideas and hypotheses, Axially-symmetric Cracks. Prikladnaya Matematika i Mekhanika. 23, 622-636.
  * [7] Bloch, A.M., Krishnaprasad, P.S., Marsden, J.E., Ratiu, T.S., 1994. Dissipation induced instabilities. Ann. de l’Institut Henri Poincaré. Analyse non linéaire. 11(1), 37-90.
  * [8] Bodovillé, G., 1999. On damage and implicit standard materials. C. R. Acad. Sci. Paris Série IIB. 327(8), 715-720.
  * [9] Bobovillé, G. de Saxcé, G., 2001. Plasticity with non linear kinematic hardening : modelling and shakedown analysis by the bipotential approach. Eur. J. Mech. A/Solids. 20, 99-112.
  * [10] Bouby, C., de Saxcé, G., Tritsch, J.B., 2009. On shakedown of structures under variable loads with a kinematic non linear and non associated hardening rule, in: Weichert, D., Ponter, A. (Eds.), Limit States of Materials and Structures: Direct Methods. Springer.
  * [11] Bouby, C., Kondo, D., de Saxcé, G., 2015. A comparative analysis of two formulations for non linear hardening plasticity models: Application to shakedown analysis. European Journal of Mechanics A/Solids. 53, 48-61.
  * [12] Bourdin, B., Francfort, G., Marigo, J.J., 2008. The variational approach to fracture. J of Elasticity. 91, 5-148.
  * [13] Brezis, H., Ekeland, I., 1976. Un principe variationnel associé à certaines équations paraboliques. I. Le cas indépendant du temps, II. Le cas dépendant du temps. C. R. Acad. Sci. Paris Série A-B. 282, 971-974, and 282, 1197-1198.
  * [14] Buliga, M., 1999. Energy minimizing brittle fracture propagation. J. of Elasticity. 52, 201-238.
  * [15] Buliga, M., 2009. Hamiltonian inclusions with convex dissipation with a view towards applications. Math Appl 1(2), 228–225.
  * [16] Buliga, M., de Saxcé, G., 2016. A symplectic Brezis-Ekeland-Nayroles principle. Mathematics and Mechanics of Solids. 22(6), 1-15.
  * [17] Cao, X., Oueslati, A., Nguyen, A.D., de Saxcé, G., 2020. Numerical simulation of elastoplastic problems by Brezis-Ekeland-Nayroles non-incremental variational principle. Computational Mechanics. 65(4), 1006-1018
  * [18] Cao, X., Oueslati, A., Shirafkan, N., Bamer, F., Markert, B., de Saxcé, G., 2021. A non-incremental numerical method for dynamic elastoplastic problems by the symplectic Brezis-Ekeland-Nayroles principle. Computer Methods in Applied Mechanics and Engineering. 384, 11908.
  * [19] Cao, X., Oueslati, A., de Saxcé, G., 2021. A non-incremental approach for elastoplastic plates basing on the Brezis-Ekeland-Nayroles principle. Applied Mathematical Modelling. 99, 359-379.
  * [20] Chambolle, A., Francfort, G.A., Marigo, J.J., 2009. When and how do cracks propagate? J of the Mechanics and Physics of Solids. 57, 1614-1622.
  * [21] Cherepanov, G.P., 1963. One problem of indentor testing with the formation of cracks. Prikladnaya Matematika i Mekhanika. 27(1) (in Russian).
  * [22] Chopin, J., Bhaskar, A., Jog, A., Ponson, L., 2018. Depinning dynamics of crack fronts. Phys. Rev. Lett. 121, 235501.
  * [23] Cotterell, B., Rice, J.R., 1980. Some remarks. on elastic crack-tip stress fields. Int. J. of Fracture. 16, 155-169
  * [24] de Saxcé, G., Feng, Z.Q., 1991. New inequation and functional for contact with friction : the implicit standard material approach. Int. J. Mech. of Struct. and Machines. 19(3), 301-325.
  * [25] de Saxcé, G., 1992. Une généralisation de l’inégalité de Fenchel et ses applications aux lois constitutives. C. R. Acad. Sci. Paris Série II. 314, 125-129.
  * [26] de Saxcé, G., Bousshine, L., 1998. Limit Analysis Theorems for the Implicit Standard Materials: Application to the Unilateral Contact with Dry Friction and the Non Associated Flow Rules in Soils and Rocks. Int. J. Mech. Sci. 40(4), 387-398.
  * [27] de Saxcé, G., 2002. Implicit standard materials, in: Weichert, D., Maier, G. (Eds.), Inelastic Behaviour of Structures under Variable Repeated Loads, CISM International Centre for Mechanical Sciences, Courses and Lectures N∘ 432\. Springer.
  * [28] Dolbow, J., Moës, N., Belytschko, T., 2000. Discontinuous enrichment in finite elements with a partition of unity method. Finite Elements in Analysis and Design. 36, 235-260.
  * [29] Dolbow, J., Moës, N., Belytschko, T., 2000. Modeling fracture in Mindlin-Reissner plates with the extended finite element method. Int. J. of Solids and Structures. 37, 7161-7183.
  * [30] Dugdale, D.S., 1960. Yielding of steel sheets containing slits. J. Mech. Phys. Sol. 8, 100-104.
  * [31] Edelen, D.G.B, 1985. Applied Exterior Calculus. John Wiley and sons.
  * [32] Erdogan, F., Sih, G.C., 1963. On the crack extension in plates under plane loading and transverse shear. Trans. ASME, ser. D. 85(4), 519-527.
  * [33] Fenchel, W., 1949. On conjugate convex functions. Canadian Journal of Mathematics. 1, 1105.
  * [34] Fett, T., Munz, D., 2002. Kinked cracks and Richard fracture criterion. Int. J. of Fracture. 115, L69-L73
  * [35] Fitzpatrick, S., 1988. Representing monotone operators by convex functions, in: Workshop/Miniconference on Functional Analysis and Optimization. Canberra, Proc. Centre Math. Anal. Austral. Nat. Univ., 20, Austral. Nat. Univ, pp. 59-65.
  * [36] Francfort, G.A., Marigo, J.J., 1998. Revisiting brittle fracture as an energy minimization problem. J. Mech. Phys. Solids. 46, 1319–1342.
  * [37] Ghoussoub, N., MacCann, R.J., 2004. A least action principle for steepest descent in non-convex landscape. Contemporary Mathematics. 362, 177-187.
  * [38] Goldstein, R.V., Salganik, R.L., 1974. Brittle fracture of solids with arbitrary cracks. Int. J. of Fracture. 10, 507–23.
  * [39] Grmela, M., Öttinger, H.C., 1997. Dynamics and thermodynamics of complex fluids. I. Development of a general formalism. Phys. Rev. E 56(6), 6620-6632.
  * [40] Gurtin, M., Podio-Guidugli, P., 1998. Configurational forces and a constitutive theory for crack propagation that allows for kinking and curving. J. Mech. Phys. Solids. 46, 1343–1378.
  * [41] Halphen, B., Nguyen Quoc, S., 1975. Sur les matériaux standard généralisés. Journal de Mécanique. 14, 39-63.
  * [42] He, M., Hutchinson, J., 1989. Crack deflection at an interface between dissimilar elastic materials. Int. J. Solids Struct. 25, 1053–1067.
  * [43] Hellen, T.K., Blackburn, W.S., 1975. The calculation of stress intensity factors for combined tensile and shear loading. Int J Fract. Mech. 11, 605–17.
  * [44] Hjiaj, M., Bobovillé, G. de Saxcé, G., 2000. Matériaux viscoplastiques et loi de normalité implicites. C. R. Acad. Sci. Paris Série IIb. 328, 519-524.
  * [45] Hussain, M.A., Pu, S.L., Underwood, J.H., 1974. Strain energy release rate for a crack under combined mode I and mode II. Fract Anal ASTM STP. 560, 2–28.
  * [46] Kolvin, I., Cohen, G., Fineberg, J., 2015. Crack front dynamics: the interplay of singular geometry and crack instabilities. Phys. Rev. Lett. 114, 175501.
  * [47] Lebihain, M., Leblond, J.B., Ponson, L., 2020. Effective toughness of periodic heterogeneous materials: the effect of out-of-plane excursions of cracks. J of the Mechanics and Physics of Solids. 137, 103876.
  * [48] Lebihain, M., Ponson, L., Leblond, J.B., Kondo, D., 2021. Effective toughness of disordered brittle solids: A homogenization framework. J of the Mechanics and Physics of Solids. 153, 104463
  * [49] Leblond, J.B., 1989. Crack paths in plan situation – I. General form of the expansion of the stress intensity factors. Int J Solids Struct. 25, 1311–25.
  * [50] Leguillon, D., 2002. Strength or toughness? A criterion for crack onset at a notch. Eur. J. Mech. A/Solids. 21, 61-72.
  * [51] Lemaitre, J., Chaboche, J.L., 2012. Mechanics of solid materials, Cambridge University Press.
  * [52] Mangiarotti, L., Modugno, M., 1983. Some results on the calculus of variation on jet spaces, Ann. Inst. H. Poincaré. 23(1), 29-43.
  * [53] Manh Hong D., Peletier, M.A., Zimmer, J., 2013. GENERIC formalism of a Vlasov-Fokker-Planck equation and connection to large-deviation principles. Nonlinearity. 26, 2951-2971.
  * [54] Magnier, V., Charkaluk, E., Bouby, C., de Saxcé, G., 2006. Bipotential Versus Return Mapping Algorithms: Implementation of Non-Associated Flow Rules, in: Topping, B.H.V., Montero, G., Montenegro, R. (Eds.), Proceedings of The Eighth International Conference on Computational Structures Technology (las Palmas de Gran Canaria, sept. 12-15, 2006). Civil-Comp Press, Stirlingshire, United Kingdom.
  * [55] Mielke, A., Theil, F., 1999. A mathematical model for rate-independent phase transformations with hysteresis, in: Alber, H.D., Balean, R., Farwig, R. (Eds.), Workshop on Models of Continuum Mechanics in Analysis and Engineering, Shaker-Verlag, pp.117-129.
  * [56] Mielke, A., 2005, in: Feireisl, E. (Ed.), Handbook of Differential Equations, Evolutionary Equations, vol. 2. Elsevier, pp.461-559.
  * [57] Mielke, A., Roubíček, T., 2006. Rate-independent damage processes in nonlinear elasticity. Mathematical Models and Methods in Applied Sciences (M3AS). 16(2), 177-209.
  * [58] Mielke, A., 2011. Formulation of thermo-elastic dissipative material behavior using GENERIC. Contin. Mech.Thermodyn. 23, 233-256.
  * [59] Moës, N., Dolbow, J., Belytschko, T., 1999. A finite element method for crack growth without remeshing. Int. J. for Numerical Methods in Engineering 46:131-150.
  * [60] Moës, N., Belytschko, T., 2002. Extended finite element method for cohesive crack growth. Engineering Fracture Mechanics. 69, 813-833.
  * [61] Mróz, K.P., Mróz, Z., 2010. On crack path evolution rules. Engineering Fracture Mechanics. 77, 1781-1807
  * [62] Mumford, D., Shah, J., 1989. Optimal approximation by piecewise smooth functions and associated variational problems. Comm. on Pure and Appl. Math. 42(5), 577–685.
  * [63] Nayroles, B., 1976. Deux théorèmes de minimum pour certains systèmes dissipatifs. C. R. Acad. Sci. Paris Série A-B. 282, A1035-A1038.
  * [64] Öttinger, H.C., Grmela, M., 1997. Dynamics and thermodynamics of complex fluids. II. Illustrations of a general formalism. Phys. Rev. E. 56(6), 6633-6655.
  * [65] Richard, H.A., 1984. Examination of brittle fracture criteria for overlapping mode I and II loading applied to cracks, in: Sih, G.C. et al. (Eds.), Applications of Fracture Mechanics to Materials and Structures. Nijhoff Publ., The Hague, pp. 309-316.
  * [66] Rockafellar, R.T., 1970. Generalized Hamiltonian equations for convex problems of Lagrange. Pacific J. of Math. 33(2), 411-427.
  * [67] Saunders, D.J., 1989. The Geometry of Jet Bundles, Cambridge University Press.
  * [68] Sih, G.C., 1973. Some basic problems in fracture mechanics and new concepts. Engineering Fracture Mechanics. 5, 365-377
  * [69] Stefanelli, U., 2008. A variational principle for hardening elasto-plasticity. SIAM J. Math. Anal. 40(2), 623-652.
  * [70] Stolz, C., 1995. Functional approach in nonlinear dynamics. Arch. Mech. 47, 421-435.
  * [71] Strifors, H.C., 1973. A generalized force measure of conditions at crack tips. Int. J. Solids Struct. 10,1389–404.
  * [72] Tvergaard, V., Hutchinson, J.W., 1992. The relation between crack growth resistance and fracture process parameters in elastic-plastic solids. J. Mech. Phys. Sol. 40, 1377-1397.
  * [73] Vallée, V., Lerintiu, C., Fortuné, D., Ban, M., de Saxcé, G., 2005. A bipotential expressing simultaneous ordered spectral decomposition between stress and strain rate tensor, in: International conference New Trends in Continuum Mechanics (ed Theta), Constanta (Romania), 8-12 September 2003. Published under the title ”Hill’s bipotential”, New Trends in Continuum Mechanics, pp.339-351.
  * [74] Visintin, A., 2008. Extension of the Brezis-Ekeland-Nayroles principle to monotone operators. Adv. math. Sci. Appl. 18, 633-650.
  * [75] Visintin, A., 2013. Structural stability of rate-independent nonpotential flows. Discrete and Continuous Dynamical Systems Series S. 6, 257-275.
  * [76] Wu, C.H., 1978. Fracture under combined loads by maximum-energy-release-rate-criterion. ASME J Appl Mech. 45, 553–8.
  * [77] Xu, X., Needleman, A., 1994. Numerical simulations of fast crack growth in brittle solids. J. Mech. Phys. Sol. 42, 1397-1434.
  * [78] Zouain, N., Pontes Filho, I., Vaunat, J., 2010. Potentials for the modified Cam-Clay model. European Journal of Mechanics A/Solids. 29, 327-336.