<figure><img src="image_1.png" /></figure>

HAL
open science

# Gradient damage analysis of a cylinder under torsion: Bifurcation and size effects

Tianyi Li, Radhi Abdelmoula

► To cite this version:

Tianyi Li, Radhi Abdelmoula. Gradient damage analysis of a cylinder under torsion: Bifurcation and size effects. Journal of Elasticity, 2021, 143, pp.209-237. 10.1007/s10659-021-09815-x. hal-03103019v2

HAL Id: hal-03103019

https://hal.science/hal-03103019v2

Submitted on 8 Feb 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Gradient Damage Analysis of a Cylinder Under Torsion: Bifurcation and Size Effects

Tianyi Li<fnref n="1" /> · Radhi Abdelmoula<fnref n="2" />,<fnref n="3" />

Received: 5 June 2020 / Accepted: 12 January 2021
© The Author(s) 2021

**Abstract** In this work, an elastic-damage evolution analysis is carried out for a cylinder under torsion made of a material obeying a gradient damage model with softening. Both semi-analytical and asymptotic approaches are developed to analyze the elastic, axisymmetric and bifurcation stages. We show the existence of a fundamental branch where the damage field is asymmetric and localized within a finite thickness from the boundary. By minimizing a generalized Rayleigh quotient, the bifurcation time and modes are obtained as a function of the length scale $\epsilon = l/R$ involving a material internal length and the cylinder radius. We will then focus on these size effects by assuming that $\epsilon$ is a small parameter in an asymptotic setting. After justification, specific spatial and temporal rescaled variables are introduced for the boundary layer problem. It is shown that the axisymmetric damage evolution and the bifurcation are governed by two universal functions independent of the length scale. The simulation results obtained by the semi-analytical approach are formally justified by the asymptotic methods.

**Keywords** Gradient damage models · Variational methods · Bifurcation · Asymptotic method · Size effects

**Mathematics Subject Classification (2010)** 47J30 · 74G10 · 74G60 · 74R10

<figure><img src="image_3.png" /></figure>

<footnote>R. Abdelmoula
radhi.abdelmoula@univ-paris13.fr; radhi.abdelmoula@upmc.fr</footnote>

<footnote>T. Li
tianyikillua@gmail.com</footnote>

<footnote n="1">1 Paris, France</footnote>

<footnote n="2">2 LSPM (CNRS-UPR 3407), Université Paris XIII, 99 avenue Jean-Baptiste Clement, 93430 Villetaneuse, France</footnote>

<footnote n="3">3 Institut Jean Le Rond d'Alembert (CNRS-UMR 7190) Sorbonne Universités, Univ. Paris 6, 4 Place Jussieu, 75252 Paris Cedex 05, France</footnote>
# 1 Introduction

Since the advent of gradient damage models [26], these phase-field methods are gaining popularity among the computational mechanics community to simulate fracture phenomena of various physical natures, see [22] for a quick overview. Considerable efforts have been made to better understand the behaviors of these models governed by variational principles. Initially, the energy functional being minimized is conceived as a regularized approximation to the variational approach to fracture [5, 8], where the vanishing parameter $\ell$ is purely numerical in nature [4]. Recently, it has been recognized that this parameter $\ell$ can also be interpreted as an internal length scale characteristic of the material, see [6, 9, 25, 32]. In [16, 35], a fundamental link has been established between the gradient damage models and Griffith's fracture model. When $\ell$ is small compared to the characteristic size of the body, the temporal propagation of an existing phase-field crack is asymptotically governed by a Griffith-type law involving a generalized energy release rate $G^\ell$ and a corresponding fracture toughness $G_c^\ell$. In [38], crack nucleation induced by stress concentrations is investigated numerically. The internal length $\ell$ introduces a size-effect in structural failure governed by a competition between the material strength and a Griffith toughness criterion.

In this paper, we will focus on more fundamental properties of gradient damage models and particularly on the possibility of bifurcation and the associated size-effects. In a unidimensional setting, a large body of literature has been devoted to the bifurcation and stability issues of the gradient damage models, see [3, 27–29] but also a quick summary in [22]. The uniqueness and stability properties of the initial elastic solution depend strongly on the damage constitutive functions in the energy functional. In general, the homogeneous damage field is more likely to become unstable and bifurcate to a localized solution for sufficiently long bars compared to the internal length. In [36], initiation of periodic crack array in a semi-infinite domain is studied following analytical and numerical approaches. The spatial distribution of cracks can be explained by bifurcation of a fundamental branch. Since the internal length is the only length scale in the problem, it is shown that the bifurcation time is proportional to $\ell^2$, while the wavelength scales with $\ell$. The objective of this paper is to further explore bifurcation properties of the gradient damage models and the corresponding the size-effects.

In [11, 13], a multiscale echelon structure of semi penny-shaped cracks is observed under mixed mode I and III loading. In [37], the development of lance-like fracture facets is investigated for a cylinder under combined tensile and antiplane shear stresses. A mechanism of the combined growth of the parent and daughter cracks under mixed loading is suggested in [18, 30] to capture the formation of the echelon crack patterns. The most intriguing phenomena observed in these experiments are the presence of repetitive fracture patterns, which we believe can be explained by the bifurcation mechanism in gradient damage models. In [31], such phase-field models are used in the numerical simulation for mixed mode I and III problems. It is found that the echelon cracks can be captured only if the material internal length parameter is small compared to the structure dimension. In this paper, we propose a theoretical analysis of such formation of repetitive damage patterns by performing a bifurcation analysis of the gradient damage models. Instead of considering the full complex three-dimensional mixed-mode problem, which is beyond the scope of this study, we propose to analyze a simplified pure antiplane shear problem with a cylindrical bar subject to torsional loading. It will be shown that a periodic damage pattern can be formed according to the variational principles of gradient damage models.
The loading prescribed on the cylindrical bar is controlled by an increasing angle $t$ of torsion on the upper section of cylinder. Three elastic-damage evolution stages will be analyzed: elastic, axisymmetric and bifurcation, as indicated by the following diagram

$$ t: 0 \xrightarrow{(u^F, 0)} t_c \xrightarrow{(u^F, \alpha^F)} t_* \xrightarrow{(u^F, \alpha^F) + (u^B, \alpha^B)} t_r $$

The elastic stage is characterized by an identically zero damage while the displacement solution $u^F$ satisfies static equilibrium. The existence of this stage is conditioned by the choice of damage constitutive functions, see [22]. The variational damage criterion will be reached at a finite critical time $t = t_c$. In the asymmetric stage, damage doesn't depend on the angular variation and can be expressed by $\alpha = \alpha(r)$. This is the fundamental branch from which bifurcation possibilities will be explored. Similarly to the thermal shock problem, the damage solution in the fundamental branch is not spatially homogeneous, since the damage criterion is not satisfied everywhere in the domain. For the bifurcated solution, we are looking for an out-of-plane displacement field $u_z(r, \theta)e_3$ independent of axial variable $z$. These approximations can be justified in the following settings

* We assume in this work that the length of cylinder $L$ is relatively small such that a two-dimensional analysis can be carried out. More precisely, as can be seen from the prescribed boundary conditions (1), the displacement is left free in the axial direction such that we only seek the displacement solutions of type $u + we_3$ for arbitrary $w \in \mathbb{R}$. In such setting, it is known that for a circular section there is no out-of-plane warping displacement of the section [21]. We will show further that the axisymmetric damage solution also does not induce warping: at first approximation, damage remains translationally invariant along the axial direction.

* The boundary conditions for the damage variable are of Neumann type so that the damage does not explicitly depend on the axial variable.

In this paper, both the semi-analytical and asymptotic approaches are developed. In Sect. 2, we recall the main ingredients of gradient damage models and define the problem settings for the torsion problem. Section 3 is devoted to the semi-analytical approach, where the elastic-damage evolution problem is solved analytically without further physical approximations. Numerical methods are used when closed-form solutions cannot be obtained, hence the solutions are both analytical and numerical in nature. Using the variational principles of the gradient damage models, we show the existence of a fundamental branch where the damage field is asymmetric, localized in a finite crown and decreases monotonically from a maximum value at the boundary to zero with a finite thickness. The axisymmetric damage is solved and is expressed using the Bessel and Struve functions. Numerical results are presented to illustrate the influence of the length scale $\ell/R$. Torque evolution is also computed to characterize structural softening behavior of the cylinder.

The analysis of bifurcation from this fundamental branch is then performed. Similar to [36], bifurcated solutions follow a periodic localized damage pattern and are sought using a partial Fourier series in the circumferential direction characterized by a wave number. We then arrive at two energy minimization problems involving a generalized Rayleigh quotient: one at a fixed wave number that governs the bifurcated displacement and damage solutions, the other with respect to the wave number. The first problem is equivalent to a generalized eigenvalue problem where the smallest eigenvalue is to be determined. Using the finite element method for spatial discretization and numerical algorithms, the Rayleigh quotient is obtained as a function of time, wave number and the length scale $\ell/R$. On the one hand, numerical simulations suggest that the axisymmetric solution will always lose uniqueness at a
time approaching $t_c$ as $\ell/R$ becomes smaller. The wave number, on the other hand, appears to be inversely proportional to $\ell/R$. These two numerical observations are the main results of this paper, and call for a more formal derivation using the asymptotic approach.

In Sect. 4, the elastic-damage evolution problem is solved under the assumption that the length scale $\epsilon = \ell/R$ is small following an asymptotic approach. Using the perturbation technique and the boundary layer method, we construct an asymptotic approximation of the semi-analytical solution by using two rescaled spatial and temporal coordinates and by performing an asymptotic expansion with respect to $\epsilon$. It is shown that the axisymmetric damage solution can be expressed with the help of a universal function that does not explicitly depend on $\epsilon$. The influence of the length scale $\ell/R$ can be successfully explained by the introduced rescaled spatial and temporal variables. Comparison with the previous semi-analytical results shows satisfying agreements.

We then continue the asymptotic analysis for the bifurcation from the fundamental branch. By using the definition of the Rayleigh quotient and thanks to the variational nature of the formulation, the orders of magnitude of the bifurcated solutions can be obtained with respect to the small parameter $\epsilon$. We obtain again a universal expression for Rayleigh quotient that does not explicitly depend on $\epsilon$. Its temporal evolution is governed by that of the previous universal function. The asymptotic Rayleigh quotient varies monotonically from $\infty$ at the onset of damage initiation to 0, which shows the existence of a universal bifurcation time independent of $\epsilon$. The size-effects previously obtained by the semi-analytical approach are formally justified by the asymptotic method.

# 2 Gradient Damage Modeling

## 2.1 Variational Formulation

We consider a three-dimensional homogeneous and isotropic body whose natural reference configuration $\Omega$ is a straight cylinder of length $L$ with a cross circular section $S$ of radius $R$, see Fig. 1. Formally, we thus have

$$ \Omega = S \times [0, L], \quad S = \{(x, y) \in \mathbb{R}^2 | x^2 + y^2 \le R^2\}. $$

A frame of cylindrical coordinate system $(r, \theta, z)$ along with its unit vectors $(e_r, e_\theta, e_\zeta)$ is attached to the body $\Omega$. The lateral surface of the cylinder is stress-free

$$ \sigma n = 0 \quad \text{on } \partial S \times (0, L), $$

and the body is subject to the following displacements and stress conditions on the lower and upper cross sections

$$ \left\{ \begin{array}{l} u_r = u_\theta = 0 \quad \text{and} \quad \sigma_{33} = 0 \quad \text{on } S \times \{0\}, \\ u_r = 0, \quad u_\theta = \text{tr} \quad \text{and} \quad \sigma_{33} = 0 \quad \text{on } S \times \{L\}, \quad \text{for } t \ge 0, \end{array} \right. \quad (1) $$

where $t$ denotes the increasing angle of torsion applied to the upper section and is regarded as a loading parameter. The geometric and boundary conditions of the cylinder under torsion is summarized in Fig. 1.

The quasi-static structural evolution of the cylinder is described by the equilibrium state $(u, \alpha)$ of the displacement and damage fields at time $t$. We assume that at $t = 0$ the body
<figure><img src="image_4.png" /><figcaption>**Fig. 1** Geometric and boundary conditions of the cylinder
$\Omega = S \times [0, L]$ under torsion</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 2** The discrete crack $\Gamma \subset \Omega$ is regularized by a continuous damage field $0 \le \alpha \le 1$ in (a). The cross-section of the damage perpendicular to the phase-field crack is characterized by an optimal damage profile $\alpha_*$ indicated in (b) obtained with the damage constitutive model (14)</figcaption></figure>

is at rest $u_0 = 0$ and the cylinder is undamaged $\alpha_0 = 0$. The evolution problem consists in finding $t \mapsto (u, \alpha)$ for all $0 < t < t_r$ where $t_r$ corresponds to the onset of fracture, i.e., when $\max \alpha(\mathbf{x}) = 1$.

The cylinder is composed of a gradient-damage material. Its general formulation under a quasi-static setting is now recalled as follows. For a more complete physical and numerical overview, the interested readers are referred to [22, 26]. Contrary to a sharp interface description of cracks, the gradient damage approach introduces a continuous phase field replacing strong displacement discontinuities by strain localizations within a finite band. The damage parameter $\alpha$ is a scalar which can only grow from 0 to 1. It ensures a smooth transition between the undamaged state $\alpha = 0$ and the crack $\alpha = 1$, see Fig. 2(a).

The evolution of such gradient-damage bodies is completely governed by variational arguments based on the definition of a total potential energy. Given the displacement and damage fields $(u, \alpha)$ belonging to appropriate functional spaces, the total potential energy is defined by

$$ \mathcal{P}(u, \alpha) = \int_{\Omega} W(\boldsymbol{\varepsilon}(u)(\mathbf{x}), \alpha(\mathbf{x}), \nabla\alpha(\mathbf{x})) \, d\mathbf{x}, \qquad (2) $$

where $W$ is the bulk energy density depending on the linearized strain tensor $\boldsymbol{\varepsilon} = \frac{1}{2}(\nabla u + \nabla^T u)$, damage $\alpha$ and its gradient. This bulk energy density is the sum of two components: the stored elastic energy density and the dissipated energy density which quantifies the amount of energy consumed up to a damage state

$$ W(\boldsymbol{\varepsilon}, \alpha, \nabla\alpha) = \frac{1}{2} A(\alpha) \boldsymbol{\varepsilon} \cdot \boldsymbol{\varepsilon} + w(\alpha) + w_1 \ell^2 \nabla\alpha \cdot \nabla\alpha. \qquad (3) $$

1. The first term in (3) is the elastic energy density characterized by a damage-dependent elasticity tensor $A(\alpha)$ which represents local stiffness degradation due to damage. The
corresponding damage-dependent stress tensor reads $\sigma = A(\alpha)\varepsilon$. As in previous work, we assume that the Poisson's ratio is unaffected by the damage variable, which leads to

$$A(\alpha) = a(\alpha)A_0,$$

where $\alpha \mapsto a(\alpha)$ denotes the non-dimensional stiffness degradation function with $a(1) = 0, a(0) = 1$.

2. The last two terms in (3) are the dissipated energy density, which introduces a length parameter $\ell$ that controls the damage band width from a geometric point of view (see Fig. 2). In (3), $\alpha \mapsto w(\alpha)$ represents another damage constitutive function representing the local damage dissipation during a homogeneous damage evolution and its maximal value $w(1) = w_1$ is the energy completely dissipated during such process when damage attains 1. Contrary to local strain-softening constitutive models, here the damage dissipation mechanism becomes non-local and localization is systematically accompanied by finite energy consumption, due to the presence of the gradient term, see [2, 3].

To properly define the total potential energy, the admissible functional spaces for $u$ and $\alpha$ must be specified. Before final failure $t < t_r$, the displacement field belongs to the classical Sobolev space $H^1(\Omega, \mathbb{R}^3)$ to ensure that the total potential energy that will be defined below is finite. The kinematically admissible space $C_t$ takes into account the essential boundary conditions mentioned in (1). It is an affine space and its associated linear space is denoted by $C_0$. Formally, they are respectively defined by

$$ \begin{aligned} C_t &= \{ \mathbf{v} \in H^1(\Omega, \mathbb{R}^3) | v_r = v_\theta = 0 \quad \text{on } S \times \{0\}, \quad v_r = 0, v_\theta = \text{tr} \quad \text{on } S \times \{L\} \}, \\ C_0 &= \{ \mathbf{v} \in H^1(\Omega, \mathbb{R}^3) | v_r = v_\theta = 0 \quad \text{on } S \times \{0, L\} \}. \end{aligned} \quad (4) $$

Since damage is modeled as an irreversible defect evolution, its admissible space will be built from a current damage state $0 \le \alpha \le 1$. Due to the presence of the damage gradient in the dissipated energy density (3), the above Sobolev space $H^1(\Omega, \mathbb{R})$ is also used. Formally it is defined by

$$ \mathcal{D}(\alpha) = \{ \beta \in H^1(\Omega, \mathbb{R}) | 0 \le \alpha \le \beta \le 1 \}. \quad (5) $$

It can be seen that a virtual damage field $\beta$ is admissible, if and only if it is accessible from the current damage state $\alpha$ verifying the irreversibility condition, i.e., the damage only grows. Consequently, the admissible space for the damage rate $\dot{\alpha} \in \dot{\mathcal{D}}$ is a convex cone and is given by

$$ \dot{\mathcal{D}} = \{ \dot{\beta} \in H^1(\Omega, \mathbb{R}) | \dot{\beta} \ge 0 \}. \quad (6) $$

With all variational ingredients set, the temporal evolution of the displacement and damage pair $t \mapsto (u, \alpha)$ is governed by the following three physical principles

1. **Irreversibility**: the damage $t \mapsto \alpha$ is a non-decreasing function of time.

2. **Meta-stability**: the current state $(u, \alpha)$ is always directionally stable in the sense that for all admissible perturbations $v \in C_0$ and $\beta \in \dot{\mathcal{D}}$, there exists a $\bar{h} > 0$ such that for all $h \in [0, \bar{h}]$ we have

$$ (u + h\mathbf{v}, \alpha + h\beta) \in C_t \times \mathcal{D}(\alpha), \quad \mathcal{P}(u, \alpha) \le \mathcal{P}(u + h\mathbf{v}, \alpha + h\beta). \quad (7) $$
3. **Energy balance:** the only energy dissipation is due to damage such that we have the following energy balance

$$ \mathcal{P}(\mathbf{u}, \alpha) = \mathcal{P}(\mathbf{u}_0, \alpha_0) + \int_0^t \left( \int_{\Omega} \boldsymbol{\sigma}_s \cdot \boldsymbol{\varepsilon}(\dot{\mathbf{U}}_s) \, \mathrm{d}\mathbf{x} \right) \mathrm{d}s. \quad (8) $$

In practice, the gradient-damage evolution problem is solved by invoking the first-order necessary conditions of these physical principles. Interested readers are referred to [22] and references therein for detailed derivations based on variational reasoning. If the involved fields are sufficiently regular in space and in time, it can be shown that $\mathbf{u}$ satisfies the following static equilibrium equation

$$ \mathrm{div} \, \boldsymbol{\sigma} = \mathrm{div}(\mathcal{A}(\alpha)\boldsymbol{\varepsilon}(\mathbf{u})) = 0 \text{ in } \Omega, \quad (9) $$

while $\alpha$ evolves according to the following damage criterion and consistency condition

$$ \frac{1}{2} \mathcal{A}'(\alpha) \boldsymbol{\varepsilon}(\mathbf{u}) \cdot \boldsymbol{\varepsilon}(\mathbf{u}) + w'(\alpha) - 2w_1 \ell^2 \Delta \alpha \geq 0 \text{ in } \Omega, \quad \nabla \alpha \cdot \mathbf{n} \geq 0 \text{ on } \partial \Omega, \quad (10) $$

$$ \left( \frac{1}{2} \mathcal{A}'(\alpha) \boldsymbol{\varepsilon}(\mathbf{u}) \cdot \boldsymbol{\varepsilon}(\mathbf{u}) + w'(\alpha) - 2w_1 \ell^2 \Delta \alpha \right) \dot{\alpha} = 0 \text{ in } \Omega, \quad (\nabla \alpha \cdot \mathbf{n}) \dot{\alpha} = 0 \text{ on } \partial \Omega. \quad (11) $$

These first order conditions are not sufficient for bifurcation and stability analysis for a given state $(\mathbf{u}, \alpha)$, which requires the second-order directional derivatives of $\mathcal{P}$, see for instance [3, 29]. In the most general case, it is given by

$$ \begin{aligned} \mathcal{P}''(\mathbf{u}, \alpha)(\mathbf{v}, \beta) = & \int_{\Omega} \mathcal{A}(\alpha)\boldsymbol{\varepsilon}(\mathbf{v}) \cdot \boldsymbol{\varepsilon}(\mathbf{v}) \, \mathrm{d}\mathbf{x} + \int_{\Omega} \left( \frac{1}{2}\mathcal{A}''(\alpha)\boldsymbol{\varepsilon}(\mathbf{u}) \cdot \boldsymbol{\varepsilon}(\mathbf{u}) + w''(\alpha) \right) \beta^2 \, \mathrm{d}\mathbf{x} \\ & + \int_{\Omega} 2\mathcal{A}'(\alpha)\boldsymbol{\varepsilon}(\mathbf{u}) \cdot \boldsymbol{\varepsilon}(\mathbf{v}) \beta \, \mathrm{d}\mathbf{x} + \int_{\Omega} 2w_1 \ell^2 \nabla\beta \cdot \nabla\beta \, \mathrm{d}\mathbf{x}. \end{aligned} \quad (12) $$

## 2.2 Basic Properties of Gradient Damage Materials

We assume that the stiffness degradation function $\alpha \mapsto \mathbf{a}(\alpha)$ and the local damage dissipation function $\alpha \mapsto w(\alpha)$ verify certain physical properties which characterize the behavior of a strongly brittle material, see [27, Hypothesis 1]. During a homogeneous uniaxial traction experiment, it leads to the definition of the critical stress $\sigma_c$ beyond which damage grows and the maximal stress that the material can sustain:

$$ \sigma_c = \sqrt{\frac{2Ew'(0)}{s'(0)}}, \quad \sigma_m = \max_{\beta \in [0,1]} \sqrt{\frac{2Ew'(\beta)}{s'(\beta)}}. \quad (13) $$

In (13), $E$ denotes the Young's modulus of the undamaged stiffness tensor $\mathcal{A}_0$.

The physical properties of the gradient damage model, in particular its softening or hardening character, directly depend on these two damage constitutive functions. An abundant literature is devoted to a theoretic or numerical analysis of these damage constitutive laws. The interested readers are referred to [17, 20, 22, 38] and references therein for a discussion on this point. In this contribution, the following damage constitutive model initially proposed in [26] is used

$$ \mathbf{a}(\alpha) = (1 - \alpha)^2, \quad w(\alpha) = w_1\alpha. \quad (14) $$
If can be shown that damage does not evolve as long as a non-zero critical stress is not reached, a rather appreciated property when modeling brittle fracture. Then a strain-softening behavior is observed as damage grows for $\alpha \in (0, 1)$, which implies that the critical stress coincides with the maximal stress. According to (13), during a unidimensional bar traction test we have

$$ \sigma_c = \sigma_m = \sqrt{w_1 E}. \qquad (15) $$

The damage dissipation energy density (last two terms) in (3) can be regarded as an equivalent Griffith crack surface functional density in phase-field models, see for example [23]. An effective fracture toughness $G_c$, i.e., the energy required to create a unit Griffith-like crack surface, can be identified as the energy dissipated during the optimal damage profile creation in a uniaxial traction experiment, see [26] for a detailed discussion on this point. The optimal damage profile $\alpha_*$ can be considered as the theoretic cross-section perpendicular to a gradient-damage crack, see Fig. 2. If $x$ refers to the transverse coordinate axis centered at the crack where $\alpha_*(x) = 1$, the optimal damage profile for the (14) model is given by

$$ \alpha_*(x) = \begin{cases} (1 - \frac{|x|}{2\ell})^2 & \text{if } |x| < 2\ell, \\ 0 & \text{else.} \end{cases} \qquad (16) $$

It is also indicated in Fig. 2(b). With such optimal damage profile, the corresponding dissipated energy can be consequently computed by injecting $\alpha_*$ into (3) and then integrated along the cross-section, which leads to

$$ G_c = 4\ell \int_0^1 \sqrt{w_1 w(\beta)} d\beta = \frac{8}{3} w_1 \ell. \qquad (17) $$

This equation prescribes a relation between the fracture toughness $G_c$, the maximal damage dissipation $w_1$ and the internal length $\ell$.

# 3 Semi-Analytical Solutions

## 3.1 Purely Elastic Solution

By adopting the damage constitutive laws (14), the cylinder will first behave elastically with an identically zero damage field $\alpha = 0$ when $t$ increases from 0. Due to static equilibrium (9) and the prescribed boundary conditions (1), the corresponding displacement field is unique up to an arbitrary constant axial translation in the $e_3$ direction and its only nonzero component is the tangential displacement in the polar coordinate system

$$ \mathbf{u}^{\mathrm{F}} = \frac{trz}{L} \mathbf{e}_{\theta}. \qquad (18) $$

In the sequel, $\mathbf{u}^{\mathrm{F}}$ denotes the fundamental displacement field from which perturbations in the sense of the meta-stability condition (7) will be then considered in Sect. 3.3. Its corresponding strain and stress fields in the cylinder are given by

$$ \boldsymbol{\epsilon}^{\mathrm{F}} = \frac{tr}{2L} (\boldsymbol{e}_{\theta} \otimes \boldsymbol{e}_z + \boldsymbol{e}_z \otimes \boldsymbol{e}_{\theta}), \quad \boldsymbol{\sigma}^{\mathrm{F}} = \frac{\mu_0 tr}{L} (\boldsymbol{e}_{\theta} \otimes \boldsymbol{e}_z + \boldsymbol{e}_z \otimes \boldsymbol{e}_{\theta}), \qquad (19) $$

where $\mu_0$ denotes the shear modulus of the sound material.
Due to the consistency condition (11), the elastic solutions (18)–(19) remain valid if the damage criterion (10) is satisfied with a strict inequality almost everywhere in the cylinder. Using the fundamental displacement field (18), with an identically zero damage field, (10) now reads

$$ \frac{\mu'(0)t^2r^2}{2L^2} + w'(0) \geq 0, \qquad (20) $$

where $\mu(\alpha) = a(\alpha)\mu_0$ denotes the shear modulus of a damaged material point. As $a'(0) \le 0$ and $w'(0) \ge 0$ (stiffness degradation and dissipation hypotheses, see [27]), the damage criterion (20) will be first attained (becomes an equality) on the lateral surface at $r = R$, i.e., on a subset of measure zero in $\mathbb{R}^3$. The corresponding time will be denoted by $t_c$ (critical time) and is given by

$$ t_c = \sqrt{-\frac{2w'(0)L^2}{\mu'(0)R^2}}. $$

For the damage constitutive law (14) used in this paper, with (17), the critical time thus reads

$$ t_c = \sqrt{\frac{w_1 L^2}{\mu_0 R^2}} = \sqrt{\frac{3G_c L^2}{8\mu_0 \ell R^2}} > 0. \qquad (21) $$

It can be verified that a non-zero elastic phase indeed exists in this case. The elastic phase is summarized in (22) for $0 \le t \le t_c$, with the fundamental displacement field given by (18) and a zero damage field. The subsequent gradient-damage evolution will be considered in Sect. 3.2.

$$ t: 0 \xrightarrow{(\text{Elastic}, 0)} t_c \longrightarrow \dots \qquad (22) $$

## 3.2 Axisymmetric Damage Solution

From $t > t_c$, the damage criterion (10)–(11) is gradually uniformly (in the $e_3$ direction) and axisymmetrically (with respect to the $\theta$ coordinate) attained from the outer lateral surface $\partial S \times (0, L)$ into the cylinder. The main objective of this section is thus to analytically investigate the induced axisymmetric damage field $\alpha^F$ with the property

$$ \alpha(r, \theta, z) = \alpha^F(r). \qquad (23) $$

It can be shown that the associated displacement field is necessarily the fundamental displacement $u^F$ given by (18) and the associated stress and strain fields read as

$$ \boldsymbol{\epsilon}^{\mathrm{F}} = \frac{tr}{2L}(\boldsymbol{e}_{\theta} \otimes \boldsymbol{e}_{z} + \boldsymbol{e}_{z} \otimes \boldsymbol{e}_{\theta}), \quad \boldsymbol{\sigma}^{\mathrm{F}} = \frac{\mu(\alpha)tr}{L}(\boldsymbol{e}_{\theta} \otimes \boldsymbol{e}_{z} + \boldsymbol{e}_{z} \otimes \boldsymbol{e}_{\theta}). $$

Due to the fact that $\alpha$ only depends on $r$, the stress continues to verify the static equilibrium equation (9) and the prescribed boundary conditions (1). The axisymmetric damage phase is summarized in (24) for $t_c < t \le t_*$, with the fundamental displacement field given by (18) and the fundamental damage field.

$$ t: 0 \xrightarrow{(\text{Elastic}, 0)} t_c \xrightarrow{(\text{Axisymmetric}, \alpha^F)} t_* \longrightarrow \dots \qquad (24) $$
<figure><img src="image_6.png" /><figcaption>Fig. 3 Axisymmetric damage solution: (a) an illustration of the damage contour in an arbitrary section of the cylinder; (b) temporal evolution of the axisymmetric damage</figcaption></figure>

In the sequel, we will make use of the derived damage evolution law (10)–(11) to solve the axisymmetric damage field for $t > t_c$. Inserting the fundamental displacement field (18) into (10), we obtain the following damage criterion

$$ \mu'(\alpha) \frac{t^2 r^2}{2L^2} + w'(\alpha) - 2w_1 \ell^2 \Delta\alpha \geq 0 \quad \text{for all } r \in [0, R). \quad (25) $$

We suppose that the criterion (25) is reached with an equality within a radial thickness of $e$ at time $t$, see Fig. 3(a). Formally, the following subset of the cross-section $S^d \subset S$ can be defined

$$ S^d = [R - e, R) \times [0, 2\pi) = \left\{ (r, \theta) \in S | \mu'(\alpha) \frac{t^2 r^2}{2L^2} + w'(\alpha) - 2w_1 \ell^2 \Delta\alpha = 0 \text{ with } z = \bar{z} \right\}, \quad (26) $$

where without loss of generality $\bar{z} \in (0, L)$ designates an arbitrary cross-section.

To solve this second-order elliptic differential equation in $S^d$, several additional conditions must be specified. On the one hand, according to (11), on the lateral surface $r = R$ where damage evolves the normal damage gradient must be zero. On the other hand, due to continuity of the damage field (implicitly contained in the variational formulation due to the presence of the damage gradient in the energy functional), it is equal to zero at the interface between the undamaged part and $S^d$, see Fig. 3(b). Hence, we have

$$ \frac{\partial \alpha^F}{\partial r}(R) = 0, \quad \alpha^F(R - e) = 0. \quad (27) $$

Contrary to classical boundary value problems, here the boundary itself $S^d$ is yet to be specified due to the unknown thickness parameter $e$. Moreover, in this case the boundary $S^d$ (and thus the thickness) evolves with time as can be seen in Fig. 3(b). The determination of the thickness parameter $e$ requires hence another condition. It can be shown (see [15, 36] that the damage gradient is also continuous at the interface

$$ \frac{\partial \alpha^F}{\partial r}(R - e) = 0. \quad (28) $$

This equation completes the previous boundary conditions (27) in order to solve the differential equation (26). In summary, with the chosen damage constitutive law (14), we arrive
at the following system

$$ \left\{ \begin{array}{l} -\frac{\mu_0(1-\alpha)t^2r^2}{L^2} + w_1 - 2w_1\ell^2\Delta\alpha = 0 \quad \text{in } [R-e, R), \\ \frac{\partial\alpha^F}{\partial r}(R=0) = 0, \quad \alpha^F(R-e) = 0, \quad \frac{\partial\alpha^F}{\partial r}(R-e) = 0. \end{array} \right. \qquad (29) $$

We now propose a semi-analytical method to solve the unique axisymmetric damage field governed by the system (29) for $t \ge t_c$. The Laplacian operator is now expressed in the cylindrical coordinate system, taking into account the axisymmetry of the damage field which is also constant in the $z$ direction. Written with the non-dimensional radial coordinate $\alpha^F(r) = d(r) = d(r/R)$, the differential equation (29) becomes

$$ r d''(r) + d'(r) - \frac{1}{2} \left(\frac{tR}{t_c \ell}\right)^2 r^3 d(r) = \frac{1}{2} \left(\frac{tR}{t_c \ell}\right)^2 \left(t_c^2 r/t^2 - r^3\right) \quad \text{for } \Sigma = \frac{R-e}{R} \le r < 1, \quad (30) $$

where the non-dimensional variable $\Sigma$ denotes the interface between the purely elastic region $r < \Sigma$ and the region where the damage criterion (25) is attained with an equality. The general solution to (30) reads

$$ d(y) = 1 + C_1 I_0(y) + C_2 K_0(y) + \frac{\pi t_c R}{4\sqrt{2t\ell}} L_0(y), \quad y = \frac{tR}{2\sqrt{2t_c\ell}} r^2, \qquad (31) $$

where $I_v$ and $K_v$ are respectively the modified $v$-order Bessel function of the first kind and the second kind and $L_v$ is the modified $v$-order Struve function. Using the following identifies concerning the derivative of these special functions

$$ I'_{0} = I_{1}, \quad K'_{0} = -K_{1}, \quad L'_{0} = \frac{1}{2} \left( \frac{2}{\pi} + L_{-1} + L_{1} \right) $$

and the boundary conditions (27), the two constants $C_1$ and $C_2$ can be readily expressed as an explicit function of the thickness parameter $e$. Hence, the last boundary condition of (29) provides a nonlinear equation on $e$. The standard Newton-Raphson algorithm is applied to solve numerically this equation. The obtained damage solution is thus exact and semi-analytical in nature.

In Fig. 4(a), the damage thickness evolution is presented for three values of $\ell$ with respect to the normalized time scale $t/t_c$. Due to irreversibility and in the absence of unloading, the thickness is strictly increasing. For smaller $\ell$, the damaged zone is also decreased at a fixed $t/t_c$. The evolution is non-linear in general, however its initial slope near $t = t_c$ seems to be independent of $\ell$. This property will be established in Sect. 4.

At a fixed $t/t_c = 1.1$, the damage field along the radial direction $0 \le r \le R$ for three internal lengths $\ell$ is shown in Fig. 4(b). The maximal damage on the lateral surface $r = R$ increases as the internal length $\ell$ decreases. It should be reminded that the critical time $t_c$ at which damage begins depends itself on the internal length $\ell$, see (21).

The torque $M$ generated on the surface $z = L$ characterizes the global structural response of the cylinder under torsion. Using the stress vector $F = \sigma^F e_3$ on the section $S$ during the axisymmetric stage, the torque $M$ can be calculated by

$$ \begin{aligned} M &= e_3 \cdot \int_S r e_r \times F \, ds = \int_S \frac{\mu(\alpha^F)tr^2}{L} \, ds \\ &= \frac{2\pi \mu_0 t}{L} \int_0^R a(\alpha^F)r^3 \, dr = \frac{2\pi \mu_0 t}{L} \left( \frac{(R-e)^4}{4} + \int_{R-e}^R a(\alpha^F)r^3 \, dr \right). \end{aligned} \qquad (32) $$
<figure><img src="image_7.png" /><figcaption>Fig. 4 Axisymmetric solutions: (a) damage thickness evolution; (b) damage solutions at $t/t_c = 1.1$; (c) torque evolution</figcaption></figure>

Its evolution is illustrated in Fig. 4(c). For each $ℓ$, it is normalized by its value $M_c$ at $t = t_c$. Similarly to damage thickness, its initial temporal rate of change appears to be a constant independently of $ℓ$. During the axisymmetric stage, the moment first increases and reaches a maximal value. Afterwards the cylinder displays a structural softening behavior, where the moment decreases with $t$.

## 3.3 Bifurcation Analysis of the Fundamental Branch

The axisymmetric state ($u^F, \alpha^F$) given respectively in (18) and (29) is obtained by considering only the first-order stability conditions (9), (10) and (11). This state is not guaranteed to remain always meta-stable and hence bifurcation may produce as the torsion angle $t$ increases. This section is thus mainly devoted to a bifurcation analysis of the axisymmetric state ($u^F, \alpha^F$). Several stages of the complete evolution considered in this paper are recalled in (33) for $t_* < t < t_r$, with $t_r$ indicating the instant of fracture.

$$ t: 0 \xrightarrow{(u^F, 0)} t_c \xrightarrow{(u^F, \alpha^F)} \text{Axisymmetric} \xrightarrow{t_*} \text{Bifurcation} \xrightarrow{(u^F, \alpha^F)+(u^B, \alpha^B)} t_r. \quad (33) $$

Although the axisymmetric solution in question is independent of the $z$ variable, the bifurcated state ($v, \beta$) may be not. However, such full three dimensional analyses appear mathematically less tractable. Hence, in the sequel we will only seek the displacement field $u$ in the following form

$$ \begin{aligned} u &= u^F(r, z) + u^B(r, \theta) \\ &= u^F(r, z) + u^B(r, \theta) e_3, \end{aligned} \quad (34) $$
where $\mathbf{u}^F \in C_t$ is the fundamental displacement given by (18) and $\mathbf{u}^B = \mathbf{u}^B(r, \theta) e_3 \in C_0$ denotes the axial component regarded as the new unknown kinematic field. Due to static equilibrium (9), it can be shown that the damage field $\alpha$ is also independent of $z$

$$ \alpha = \alpha^F(r) + \alpha^B(r, \theta). \quad (35) $$

By substituting (34) and (35) in the expression of the total potential energy (2), we obtain the new expression for our two-dimensional analysis

$$ \begin{aligned} \mathcal{P}(\mathbf{u}, \alpha) &= L \hat{\mathcal{P}}(\mathbf{u}^B, \alpha) \\ &= L \int_S \left( \frac{1}{2} \mu(\alpha) \left( \nabla \mathbf{u}^B \cdot \nabla \mathbf{u}^B + \left(\frac{tr}{L}\right)^2 + \frac{2t}{L} \frac{\partial \mathbf{u}^B}{\partial \theta} \right) + w(\alpha) + w_1 \ell^2 \nabla \alpha \cdot \nabla \alpha \right) \, ds, \end{aligned} \quad (36) $$

where we recall that in a Cartesian coordinate system the derivative with respect to the tangential direction can be evaluated by

$$ \frac{\partial \mathbf{u}^B}{\partial \theta}(x, y) = -y \cdot \frac{\partial \mathbf{u}^B}{\partial x}(x, y) + x \cdot \frac{\partial \mathbf{u}^B}{\partial y}(x, y). $$

The stability and bifurcation properties of the gradient damage model require the second derivative (12) of the potential energy. Within the two-dimensional setting (34)–(35), the second derivative of the total potential energy (12) evaluated at the axisymmetric state $(\mathbf{u}^F, \alpha^F)$ is a quadratic form in $(\mathbf{v}, \beta)$ given by

$$ \mathcal{P}''(\mathbf{u}^F, \alpha^F)(\mathbf{v}, \beta) = L \hat{\mathcal{P}}''(0, \alpha^F)(\mathbf{v}, \beta) = L(\mathbf{A}(\mathbf{v}, \beta) - \mathbf{B}(\beta)), \quad (37) $$

with A and B two quadratic forms defined respectively by

$$ \begin{aligned} \mathbf{A}(\mathbf{v}, \beta) &= \int_S \left( \mu_0 (1 - \alpha^F)^2 \left(\frac{\partial \mathbf{v}}{\partial r}\right)^2 + \mu_0 \left( \frac{1 - \alpha^F}{r} \frac{\partial \mathbf{v}}{\partial \theta} - \frac{2tr}{L} \beta \right)^2 + 2w_1 \ell^2 \nabla \beta \cdot \nabla \beta \right) \, ds, \\ \mathbf{B}(\beta) &= \int_S 3\mu_0 \frac{t^2 r^2}{L^2} \beta^2 \, ds, \end{aligned} \quad (38) $$

where the chosen damage constitutive law (14) is used. Note that $\mathbf{v} = \mathbf{v} e_3$ and $\beta$ represent a perturbation to the axisymmetric state $(\mathbf{u}^F, \alpha^F)$.

The uniqueness of the axisymmetric solution can be analyzed by formulating the associated rate problem (see [3, 36]) which involves the following Rayleigh quotient

$$ \lambda = \mathcal{R}(\mathbf{v}, \beta) = \begin{cases} \frac{\mathbf{A}(\mathbf{v}, \beta)}{\mathbf{B}(\beta)} & \text{if } \beta \neq 0, \\ +\infty & \text{otherwise.} \end{cases} \quad (39) $$

The current evolution $t \mapsto (\mathbf{u}, \alpha)$ is unique, if the second derivative remains positive-definite

$$ \inf \lambda > 1 \text{ for all } (\mathbf{v}, \beta) \in H^1(S) \times \mathcal{D}^d(S), \quad \mathcal{D}^d(S) = \{\beta \in H^1(S) | \beta = 0 \text{ in } S \setminus S^d\}, \quad (40) $$

where the set $\mathcal{D}^d$ denotes the set of damage rates that remain zero in the undamaged part. The objective of the bifurcation analysis is to determine $t$ such that $\lambda = 1$, i.e., when the current axisymmetric solution is no longer unique.
In this work, we will assume the following particular ($v$, $\beta$) perturbations with 3 parameters ($k$, $A_k$, $C_k$)

$$ 
\begin{aligned}
v(r, \theta) &= \sum_k A(r) \sin k\theta, \quad A \in H^1(0, 1) \equiv H^1, \\
\beta(r, \theta) &= \sum_k C(r) \cos k\theta, \quad C \in D^d(0, 1) \equiv D^d,
\end{aligned}
\quad (41) $$

where $k \in \mathbb{N}^*$. In (41), the functions $A(r)$ and $C(r)$ represent the radial variation of $v$ and $\beta$, while the wave number $k$ characterizes the frequency of their angular variation. Note that at fixed $r$, the functions $v(r, \cdot)$ and $\beta(r, \cdot)$ are square-integrable continuous periodic functions, and can be well represented by the Fourier basis. The phase $\pi/2$ between $v$ and $\beta$ results in a natural coupling between them in the Rayleigh quotient.

Substituting (41) in (38) and using some basic integration properties of the trigonometric functions, the Rayleigh quotient now reads

$$ \lambda = \mathbb{R}(k, A, C) = \frac{A(A, C)}{B(C)}, \quad (42) $$

where at a fixed wave number $k$ the two quadratic forms $A$ and $B$ read

$$ 
\begin{aligned}
A(A, C) &= \int_0^R \left( 2w_1 \ell^2 \left(C'^2 + \frac{C^2 k^2}{r^2}\right) + \mu_0 (1 - \alpha^F)^2 A'^2 \right. \\
&\qquad \left. + \mu_0 \left( \frac{1 - \alpha^F}{r} kA - \frac{2trC}{L} \right)^2 \right) r \, dr, \\
B(C) &= 3\mu_0 \frac{t^2}{L^2} \int_0^R r^3 C^2 \, dr.
\end{aligned}
\quad (43) $$

To minimize (42), we will perform first the minimization with respect to $(A, C)$, and then to the wave number $k$

1. Minimization respect to $(A, C)$. Suppose that $(A_k, C_k)$ minimizes the Rayleigh ratio (42) at fixed $k \in \mathbb{N}^*$ with a minimum value $\lambda_k = \mathbb{R}(k, A_k, C_k)$. The corresponding necessary condition reads thus

$$ A'(A_k, C_k)(\hat{A}, \hat{C}) = \lambda_k B'(C_k) \hat{C} \quad \text{for all } (\hat{A}, \hat{C}) \in H^1 \times D^d, \quad (44) $$

where $A'$ and $B'$ denotes the first directional derivative of $A$ and $B$ respectively. This equation corresponds to a generalized symmetric eigenvalue problem, where $\lambda_k$ is the smallest eigenvalue minimizing the Rayleigh quotient (42) and $(A_k, C_k)$ denotes the associated eigenvector.

2. Minimization with respect to $k \in \mathbb{N}^*$. Given $(A_k, C_k)$ the resulting minimizing eigenvector of the previous problem (44), the next step is to minimize the reduced Rayleigh quotient $k \mapsto \mathbb{R}(k, A_k, C_k)$, or equivalently the eigenvalue $\lambda_k$, for $k \in \mathbb{N}^*$. It consists of an integer optimization problem. The minimizing wave number is then referred to as the optimal wave number $k_*$. The minimum Rayleigh quotient evaluated at the optimal wave number $k_*$ will be denoted by

$$ \lambda_* = \lambda_{k_*} = \mathbb{R}(k_*, A_{k_*}, C_{k_*}), \quad k_* = \underset{k}{\operatorname{argmin}} \mathbb{R}(k, A_k, C_k). \quad (45) $$
<figure><img src="image_8.png" /><figcaption>**Fig. 5** Rayleigh quotient $\lambda_k$ as a function of the wave number $k$ for $\ell/R = 0.1$: (a) at $t/t_c = 1.06$; (b) at $t/t_c \approx 1.149$</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 6** The optimal eigenfunctions with $k_* = 6$ for $\ell/R = 0.1$ at $t/t_c \approx 1.149$: (a) normalized $A$; (b) normalized $C$</figcaption></figure>

The eigenvalue problem and the minimization are solved numerically in Appendix.

For a specific internal length $\ell/R = 0.1$, the k-dependence of the Rayleigh quotient $\lambda_k$ is illustrated in Fig. 5 for two values of torsion angles $t$. It can be seen that first at $t/t_c = 1.06$ the optimal wave number $k_* = 5$ achieves the minimum of the Rayleigh quotient and then it shifts to $k_* = 6$ when the minimum Rayleigh quotient $\lambda_* = 1$, i.e., the instant when the axisymmetric solution is no longer unique.

First let us take a look at the optimal ($v_*, \beta_*$) perturbations for $\ell/R = 0.1$ at $t/t_c \approx 1.149$. They are optimal in the sense that they are evaluated with the optimal wave number $k_*$. The corresponding eigenfunctions ($A_*, C_*$) are illustrated in Fig. 6.

The eigenfunctions ($A_*, C_*$) can not be normalized separately, we have thus normalized ($A_*, C_*$) such that the rescaled $\tilde{A}(R) = 1$. We observe that using the same normalization coefficient, the rescaled $\tilde{C}(R)$ is strictly larger than $\tilde{A}(R)$. This property will be established in Sect. 4. The function $A$ is non-zero almost everywhere in the whole radial direction $0 < r \le R$ while the function $C$ is non-zero only in the damaged part for $r > R - e$. This originates from the space definitions (41).

For a smaller internal length $\ell/R = 0.05$, the k-dependence of the Rayleigh quotient as well as the optimal perturbation function $\tilde{v}$ at the instant of bifurcation $\lambda_* = 1$ for $t/t_c \approx 1.073$ are presented in Fig. 7. Compared to $\ell/R = 0.1$, a larger optimal wave number $k_* = 11$ is found for $\ell/R = 0.05$. Also, the instant of bifurcation $t_*/t_c$ is also decreasing.

For $\ell/R = 0.05$ and $\ell/R = 0.1$, the time dependence of the optimal wave number $k_*$ and that of the minimum Rayleigh quotient $\lambda_*$ are presented in Fig. 8, where the observations above are clearly illustrated. In Fig. 8(a), while the optimal wave number $t \mapsto k_*$ is not
<figure><img src="image_10.png" /><figcaption>**Fig. 7** Bifurcation analysis for $\ell/R = 0.05$ at $t/t_c \approx 1.073$: (a) Rayleigh quotient $\lambda_k$ as a function of the wave number $k$; (b) the optimal normalized perturbation function $\tilde{v}$ with $k_* = 11$</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 8** Time dependence of the Rayleigh minimization problem for $\ell/R = 0.05$ and $\ell/R = 0.1$: (a) optimal wave number $k_*$; (b) minimum Rayleigh quotient $\lambda_*$</figcaption></figure>

constant after damage initiation, its variation is minimal ($\pm 1$ in the current time interval) and in any case we are only interested in its value at the instant of bifurcation when $\lambda_* = 1$. In Fig. 8(b), we can see the optimal Rayleigh quotient gradually decreases with time. The bifurcation time can be identified when $\lambda_* = 1$.

The instant of bifurcation $t_b/t_c$ and the optimal wave number $k_*$ are presented with respect to $\ell/R$ in Fig. 9. It illustrates size-dependance of bifurcation and will be analytically studied in Sect. 4. On the one hand, $t_*/t_c$ scales linearly with the length scale $\ell/R$. In particular, for vanishing non-zero internal length $\ell \to 0^+$ the bifurcation takes place just after damage initiation $t_* \to t_c$. On the other hand, $k_*$ scales inversely proportional with the length scale $\ell/R$ in particular for small $\ell/R$. These two properties illustrate the size-dependance of bifurcation and will be theoretically studied in Sect. 4.

# 4 Asymptotic Behaviors for Small Internal Lengths

In the asymptotic approach, the internal length $\ell$ is assumed to be a small parameter compared to the cylinder radius $R$ in an asymptotic setting. The perturbation analysis is justified for the two following reasons:

1. The differential equation of asymmetric damage (30) involves a small coefficient $\ell$ in front of the highest order derivative, which suggests a boundary layer analysis similar to [1, 34].
<figure><img src="image_12.png" /><figcaption>Fig. 9 Size effect of bifurcation: (a) bifurcation time $t_*/t_c$; (b) optimal wave number $k_*$</figcaption></figure>

2. In the bifurcation problem (or more precisely the eigenvalue problem) (44), the linear op- erators contain this small parameter $\ell$, which also calls for eigenvalue expansion analysis like [10, 12].

Therefore we investigate the elastic-damage problem in the neighborhood of the boundary $r = R$ by applying a zoom in space. A corresponding time rescaling is then carried out. The asymptotic analysis of the bifurcated solution from the fundamental solution is then performed for the boundary layer problem with respect to the small parameter. Finally, the bifurcation time and the wave number is then deducted from this analysis. For more refer- ences and applications of asymptotic method, we refer the reader to [39] for applications of this technique in fluid mechanics, [7, 24] for general boundary value problems and [1, 34] for solid mechanics applications.

Following the work of [35], we consider in this work that the radial thickness $e$ of the damaged region is of order $\ell$. With this hypothesis, the following non-dimensional variables are introduced

$$ \epsilon = \sqrt{2} \frac{\ell}{R} \ll 1, \quad \zeta = \frac{1-r}{\epsilon} = \mathcal{O}(1) \ge 0, \qquad (46) $$

where $\epsilon$ is the small parameter in our asymptotic analysis and $\zeta$ can be regarded as the new (stretched) spatial coordinate for the damage zone. Note that $\mathcal{O}(\epsilon^n)$ denotes a term of the order $\epsilon^n$. Using this change of variable (46), we propose to analyze the asymptotic properties of damage evolution for small values of $\epsilon$.

## 4.1 Axisymmetric Damage Solution

We will first introduce a new loading parameter $t = \mathcal{O}(1)$ compatible with the spatial coor- dinate $\zeta$ introduced in (46) by the following

**Proposition 1** The new loading parameter $t$ that governs damage evolution is defined by

$$ t = \frac{1 - t_c / t}{\epsilon} = \mathcal{O}(1) \ge 0. \qquad (47) $$

*Proof* To prove the above relation, we perform asymptotic expansion in the damage region near the boundary $r = R$. Rewriting the axisymmetric equation (30) using (46), we obtain

$$ -\frac{t^2}{t_c^2} \epsilon^2 (1 - d(\zeta))(1 - 3\epsilon\zeta + \dots) + 1 - \epsilon\zeta - ((1 - \epsilon\zeta)d''(\zeta) - \epsilon d'(\zeta)) + \mathcal{O}(\epsilon^2) = 0, $$
$$ \forall \zeta \in (0, \Sigma_{\zeta}), \qquad (48) $$

where $\Sigma_{\zeta}$ is the interface $\Sigma$ in (30) expressed in the stretched variable $\zeta$

$$ \Sigma_{\zeta} = \frac{1 - \Sigma}{\epsilon}. \qquad (49) $$

The boundary conditions in (29) are also rewritten using the new spatial variable $\zeta$, leading to

$$ d'(0) = 0, \quad d'(\Sigma_{\zeta}) = d(\Sigma_{\zeta}) = 0. \qquad (50) $$

We assume an asymptotic expansion of damage field $d(r)$ in the form

$$ d(r) = d_0(\zeta) + \epsilon \cdot d_1(\zeta) + \dots \qquad (51) $$

Inserting (51) into (48), we obtain at order $O(1)$

$$ d''_{0}(\zeta) - \frac{t^2}{t_c^2} d_{0}(\zeta) = 1 - \frac{t^2}{t_c^2} \implies d_{0}(\zeta) = A \cosh \zeta + B \sinh \zeta + 1 - \frac{t_c^2}{t^2}, $$

where A and B are two constants. Using the boundary conditions $d'_{0}(0) = d'_{0}(\Sigma_{\zeta}) = 0$ derived from (50), we obtain $A = B = 0$, which implies

$$ d_0(\zeta) = 1 - \frac{t_c^2}{t^2}. $$

But the above expression do not verify the another boundary condition $d_{0}(\Sigma_{\zeta}) = 0$. We conclude hence that $d_{0} = O(\epsilon)$ and we deduce that

$$ 1 - \frac{t_c}{t} = O(\epsilon). $$

The proposed loading parameter (47) is thus justified. $\square$

Using the new loading parameter (47) in (48), we obtain thus

$$ (1 - d(\zeta))(-1 + \epsilon(-2t + 3\zeta + \dots)) + 1 - d''(\zeta) + \epsilon(-\zeta + d''(\zeta) + d'(\zeta)) + O(\epsilon^2) = 0, \quad \forall \zeta \in (0, \Sigma_{\zeta}). \qquad (52) $$

With the asymptotic expansion (51)), the differential equations at order $\epsilon^0$ and $\epsilon$ are given by

$$ \text{At order } \epsilon^0: \quad d''_0(\zeta) - d_0(\zeta) = 0, \qquad (53) $$

$$ \text{At order } \epsilon^1: \quad d''_1(\zeta) - d_1(\zeta) + (2t - 3\zeta)(1 - d_0(\zeta)) + \zeta - d''_0 - d'_0 = 0, \qquad (54) $$

with $d_0(\zeta)$ and $d_1(\zeta)$ satisfying the boundary conditions (50)). This leads to

$$ \begin{aligned} & d_0(\zeta) = 0, \\ & d_1''(\zeta) - d_1(\zeta) - 2\zeta + 2t = 0. \end{aligned} \qquad (55) $$
The interface $\Sigma$ can be written with the new loading parameter $t$ as

$$ \Sigma = 1 - \epsilon \cdot f(t) + \mathcal{O}(\epsilon^2), \quad \Sigma_{\zeta} = f(t) = \mathcal{O}(1), \qquad (56) $$

where $f(t)$ is positive by construction and is a universal function in the sense that it does not depend on the internal length $\ell$. The normalized thickness is hence given by

$$ \mathbf{e} = \frac{e}{R} = 1 - \Sigma = \epsilon \cdot f(t) + \mathcal{O}(\epsilon^2). \qquad (57) $$

After calculation, we obtain the following expressions for damage

$$ d(\zeta) = 2\epsilon \left( (f(t) - t) \cosh(f(t) - \zeta) - \sinh(f(t) - \zeta) - \zeta + t \right) + \dots, \qquad (58) $$

$$ d'(\zeta) = 2 \left( (f(t) - t) \sinh(f(t) - \zeta) - \cosh(f(t) - \zeta) + 1 \right) + \dots, \qquad (59) $$

where the universal function $f(t)$ defined by (56) is governed by

$$ (t - f(t)) \sinh f(t) + \cosh f(t) = 1. \qquad (60) $$

**Remark 1** Using the rescaled spatial (46) and temporal (47) variables, the asymptotic axisymmetric damage solution (58) does not explicitly depend on the small parameter $\epsilon$ aside from the multiplication factor. The temporal evolution of the axisymmetric damage can be directly explained by $t$ and the universal function $f(t)$.

Since no-closed form solution $f(t)$ can be obtained from the nonlinear equation (60), it is solved numerically as a function of $t$ in Fig. 10. On the one hand, it can be easily shown that

$$ f(t) \approx 2t \quad \text{for } t \approx 0. \qquad (61) $$

On the other hand, using the asymptotic behaviors $\cosh(f) \sim \exp(f)/2$ and $\sinh(f) \sim \exp(f)/2$ valid for larger $t \gg 1$, we have

$$ f(t) \approx 1 + t - 2 \exp(-f(t)) \approx 1 + t \quad \text{for } t \gg 1. \qquad (62) $$

Using (57) and (61), it can then be easily deduced that

$$ \dot{e}(t) \equiv \frac{de}{dt}(t) = \frac{df}{dt}(t) \cdot \frac{t_c}{t^2} \Rightarrow \dot{e}(t_c) = 2, \qquad (63) $$

which formally proves that the normalized thickness evolves independently of the internal length at the onset of damage initiation. This result agrees well with our previous semi-analytical result in Fig. 4. The comparison with the asymptotic thickness (57) is presented in Fig. 11(a). The approximation is excellent even for moderate values of $\epsilon$ especially for $t/t_c \to 1$.

At fixed $t/t_c$, the asymptotic behavior of the damage thickness can be analyzed as follows. For small values of $\ell$, we have $t \gg 1$, hence the approximation (62) applies, which gives

$$ \mathbf{e} = \frac{e}{R} \approx \epsilon \cdot f(t) \approx \epsilon \left( 1 + \frac{1 - t_c/t}{\epsilon} \right) = 1 - \frac{t_c}{t} + \epsilon. \qquad (64) $$

This equation shows that at fixed $t/t_c \gg 1$, the thickness value decreases as $\ell$ becomes smaller, in agreement with the semi-analytical result in Fig. 4(a).
<figure><img src="image_13.png" /><figcaption>Fig. 11 Comparison between the semi-analytical (exact) axisymmetric solutions and the asymptotic approach: (a) damage thickness evolution; (b) damage solutions at $t/t_c = 1.1$; (c) torque evolution; (d) maximum values of the torque</figcaption></figure>

In Fig. 11(b), the asymptotic damage solution (58) is compared to the previous semi-analytical solution (31). The approximation is also remarkable especially for smaller values of $\ell$.

At fixed $t/t_c$, the damage value on the boundary $\alpha(R)$ is now analyzed with respect to $\ell$. Using (58) and the approximation (62) valid for small $\ell$, we have

$$ \alpha(R) \approx 2\epsilon \left( (f(t) - t) \cosh(f(t)) - \sinh(f(t)) + t \right) \approx 2\epsilon(t-1) = 2 \left( 1 - \frac{t_c}{t} - \epsilon \right), \quad (65) $$

which shows that $\alpha(R)$ increases for smaller values of $\ell$, at fixed $t/t_c$. This result is in agreement with the semi-analytical results in Fig. 4(b).

Using the asymptotic approximation (58), the torque initially defined by (32) becomes

$$ M = \frac{\pi \mu_0 t R^4}{2L} \left( 1 + 8\epsilon^2 \left( f(t)^2 - 2f(t) \cdot t \right) + \mathcal{O}(\epsilon^3) \right), \quad (66) $$
where we recall that $t$ is related to $t$ by (46). It can be deduced in particular that for $t \approx t_c$ ($t \approx 0$), $M$ evolves linearly with $t$ and is independent of $\ell$. This result confirms our previous result presented in Fig. 4(c). The comparison between the semi-analytical and the asymptotic solutions is given in Fig. 11(c). Similar to previous comparisons, the approximation by the asymptotic solution is outstanding for smaller $\ell$.

In Fig. 11(d), the maximum values of the torque achieved during the axisymmetric stage are presented for the semi-analytical and the asymptotic solutions. As expected, the difference becomes smaller as $\ell$ is decreased. It is also confirmed that the maximum torque is decreasing for smaller values of $\ell$.

## 4.2 Bifurcation Analysis of the Fundamental Branch

In this subsection, we perform an asymptotic analysis of the bifurcation of the axisymmetric solution for small values of $\ell$. Our analysis will be based on the minimization of the Rayleigh quotient (40) with the particular perturbations assumed in (41).

To begin with, we will assume that the bifurcation time is compatible with the new rescaled time parameter defined in (47). This implies in particular that the Rayleigh quotient remains of order $O(1)$ for smaller values of $\epsilon$ since at bifurcation we have necessarily $\lambda_* = 1$.

$$ \frac{1 - t_c/t_b}{\epsilon} = O(1), \quad \lambda_k = O(1). \qquad (67) $$

This assumption will be proven a posteriori.

The first step is to evaluate the relative orders of magnitude of $k$, $A_k$ and $C_k$ with respect to the small parameter $\epsilon$. Recall that they minimize the Rayleigh quotient, see (44) and (45). Since $A_k$ and $C_k$ are the solution to the eigenvalue problem (44), we can assume without the loss of generality

$$ A_k = O(\epsilon). \qquad (68) $$

The orders of magnitude of $k$ and $C_k$ will be determined from this choice.

**Proposition 2** *Assuming (68), we have the following estimates*

$$ k = O(1/\epsilon), \quad C_k = O(1). $$

*Proof* Let us rewrite A and B in (43) of the Rayleigh quotient (42) using the normalized radial coordinate $r = r/R$

$$ \begin{aligned} A(A_k, C_k) &= R^2 \int_0^1 \left( w_1 \epsilon^2 \left( C_k'^2 + \frac{C_k^2 k^2}{r^2} \right) + \frac{\mu_0}{R^2} (1-d(r))^2 A_k'^2 \right. \\ & \qquad \left. + \frac{\mu_0}{R^2} \left( \frac{1-d(r)}{r} k A_k - \frac{2tR^2}{L} r C_k \right)^2 \right) r dr, \\ B(C_k) &= 3\mu_0 \frac{t^2 R^4}{L^2} \int_0^1 r^3 C_k^2 dr. \end{aligned} \qquad (69) $$

Focus on the gradient damage term in A involving explicitly the small parameter $\epsilon$. For this term to be of the same order as $C_k$ due to (67), we deduce immediately

$$ k = O(1/\epsilon), \quad (\cdot)' = O(1/\epsilon), \qquad (70) $$
where $(\cdot)'$ denotes the derivation operator with respect to $r$. This second relation justifies the use of a stretched variable $\zeta$ introduced previously in (46).

We recall that $d(r) = O(\epsilon)$ according to (58), hence it can be neglected with respect to 1. Due to the definition (68) and with (70), the terms in $A(A_k, C_k)$ involving $A_k$ are thus of order $O(1)$. By comparison with $B(C_k)$, this means also $C_k = O(1)$. The proof is thus complete. □

**Differential Equation for $A_k$ at Fixed $k \in \mathbb{N}^*$** Integrating by parts (44) with $\hat{C} = 0$ rewritten in the normalized radial coordinate $r = r/R$, we obtain the following differential equation governing $A_k$

$$r(rA'_k(r))' - \frac{2d'(r)}{1-d(r)}r^2A'_k(r) - k^2A_k + \frac{2tR^2k}{L}r^2\frac{C_k}{1-d(r)} = 0 \quad \text{for } 0 < r < 1, \quad (71)$$

$$A'_k(0) = A'_k(1) = 0, \quad [[A_k]](\Sigma) = 0, \quad [[A'_k]](\Sigma) = 0, \quad (72)$$

where $[[\cdot]]$ denotes the usual jump operator. Note that in the undamaged region for $r < 1 - \epsilon \equiv \Sigma$, we have $d = \beta_k = 0$ which implies also $C_k = 0$. We have hence

$$r^2 A''_k(r) + r A'_k(r) - k^2 A_k(r) = 0 \quad \text{for } 0 < r < \Sigma. \quad (73)$$

Since $A_k$ is bounded, the solution to (73) is given by

$$A_k(r) = Cr^k \quad \text{for } 0 < r < \Sigma, \quad (74)$$

where $C$ is an arbitrary constant depending on $k$.

**Differential Equation for $C_k$ at Fixed $k \in \mathbb{N}^*$** After integrating by parts (44) with $\hat{A} = 0$ rewritten in the normalized radial coordinate $r = r/R$, we obtain the weak equation for $C_k$, we obtain the differential equation governing $C_k$

$$- \epsilon^2 r (r C'_k(r))' + \left( \epsilon^2 k^2 + (4 - 3\lambda_k) \frac{t^2}{t_c^2} r^4 \right) C_k(r) - 2 \frac{t}{t_c^2} \frac{L}{R^2} k r^2 (1 - d(r)) A_k(r) = 0 \\ \text{for } \Sigma \le r \le 1, \quad (75)$$

$$C_k(\Sigma) = 0, \quad C'_k(1) = 0. \quad (76)$$

We recall that the equation is reduced to the damaged interval $(\Sigma, 1)$ due to the definition of $D^d$ in (40). We can see that the small parameter $\epsilon$ appears in the highest derivative. This singular perturbation problem requires a zoom at $r = 1$.

**Asymptotic Analysis of Displacement and Damage Fields for Small Internal Lengths**

Using (67) and Proposition 2, we seek the leading-order solutions to (71) and (75) for $\Sigma \le r \le 1$ [10, 12]

$$\lambda_k = \Lambda + \dots, \quad \frac{t_c}{t_b} = 1 - \epsilon \cdot t_*, \quad k = \frac{k}{\epsilon} + \dots, \quad A_k = \epsilon \cdot a + \dots, \quad C_k = c + \dots \quad (77)$$

Using the stretched variables $\zeta$ and $\Sigma_\zeta$ defined respectively in (46) and (56) and the fact that $d(r) = O(\epsilon)$ and $d'(r) = O(1)$ according to (58) and (59), we obtain at leading order for
$$0 \le \zeta \le \Sigma_{\zeta}$$

$$\begin{aligned} a''(\zeta) - k^2 a(\zeta) + \frac{2t_c R^2 k}{L} c(\zeta) &= 0, \\ -c''(\zeta) + (k^2 + (4 - 3\Lambda)) c(\zeta) - \frac{2k L}{t_c R^2} a(\zeta) &= 0. \end{aligned} \quad (78)$$

We seek approximate solutions for $a(\zeta)$ and $c(\zeta)$ as a power series of $\zeta$

$$a(\zeta) = a_0 + a_1 \zeta + a_2 \zeta^2 + a_3 \zeta^3 + a_4 \zeta^4 + \dots, \quad c(\zeta) = b_0 + b_1 \zeta + b_2 \zeta^2 + b_3 \zeta^3 + b_4 \zeta^4 + \dots \quad (79)$$

The boundary conditions (72) and (76) are now expressed using the stretched variables

$$a'(0) = 0, \quad c'(0) = 0, \quad \llbracket a \rrbracket(\Sigma_{\zeta}) = 0, \quad \llbracket a \rrbracket'(\Sigma_{\zeta}) = 0, \quad c(\Sigma_{\zeta}) = 0 \quad \text{where } (\cdot)' = \frac{d}{d\zeta}. \quad (80)$$

Using the first two conditions of (80), we find immediately $a_1 = b_1 = 0$. Substituting (79) into (78), we obtain $a_3 = b_3 = 0$ and

$$\begin{aligned} 2a_2 - k^2 a_0 &= -\frac{2t_c k R^2}{L} b_0, & -2b_2 + (k^2 + (4 - 3\Lambda)) b_0 &= \frac{2k L}{t_c R^2} a_0, \\ 12a_4 - k^2 a_2 &= -\frac{2t_c k R^2}{L} b_2, & -12b_4 + (k^2 + (4 - 3\Lambda)) b_2 &= \frac{2k L}{t_c R^2} a_2. \end{aligned} \quad (81)$$

Using the remaining boundary conditions in (80) and the solution (74) in $r \in (0, \Gamma)$, we obtain

$$a_0 + a_2 \Sigma_{\zeta}^2 + a_4 \Sigma_{\zeta}^4 = C \Sigma^{k/\epsilon}, \quad 2a_2 \Sigma_{\zeta} + 4a_4 \Sigma_{\zeta}^3 = -kC \Sigma^{k/\epsilon-1}, \quad b_0 + b_2 \Sigma_{\zeta}^2 + b_4 \Sigma_{\zeta}^4 = 0, \quad (82)$$

where $C$ is the integration constant in (74) and $\Sigma_{\zeta}$ and $\Sigma$ are related by (49). Since $\epsilon$ is a small parameter, we have $k/\epsilon - 1 \approx k/\epsilon$ and (82) can be reduced to the following homogeneous system

$$\begin{aligned} k(a_0 + a_2 f(t)^2 + a_4 f(t)^4) + 2a_2 f(t) + 4a_4 f(t)^3 &= 0, \\ b_0 + b_2 f(t)^2 + b_4 f(t)^4 &= 0, \end{aligned} \quad (83)$$

where $\Sigma_{\zeta}$ is replaced by $f(t)$ due to (56).

Using (81), the variables ($a_2, b_2, a_4, b_4$) can be expressed using a linear combination of ($a_0, b_0$). Hence, the only unknowns in the homogeneous system (83) are ($a_0, b_0$). It admits a non-trivial solution if and only if the determinant $J$ is zero, which is given by the following quadratic equation in $\Lambda$

$$J = J(f(t), k, \Lambda) = P(f(t), k)\Lambda^2 + Q(f(t), k)\Lambda + R(f(t), k) = 0 \implies \Lambda = \Lambda(f(t), k), \quad (84)$$

with ($P, Q, R$) three polynomials involving different powers of $f$ and $k$ up to $f(t)^8 k^8$. The Rayleigh quotient $\Lambda = \Lambda(f(t), k)$ is thus the smaller positive solution to (84). Observe that the asymptotic Rayleigh quotient $\Lambda$ is not an explicit function of the time variable $t$ or $t$. Its time-dependence is through the value of the universal function $f(t)$. It is also a universal function in the sense that it also does not depend on the small parameter $\epsilon$. Its effective
<figure><img src="image_14.png" /><figcaption>**Fig. 12** Universal asymptotic optimal Rayleigh quotient $\Lambda_*$ as well as the optimal asymptotic wave number $k_*$</figcaption></figure>

$\epsilon$-dependence is explained by the rescaled time variable $t$. The optimal asymptotic wave number $k_*$ can be obtained by minimizing $\Lambda$ at fixed $f(t)$.

$$ \frac{\partial \Lambda}{\partial k}(f(t), k_*) = 0. \qquad (85) $$

Using $k_*$, we obtain thus the optimal Rayleigh quotient $\Lambda_*(t) = \Lambda(f(t), k_*)$. We can show that for $t \approx 0$, the Rayleigh quotient is given by

$$ \Lambda = \frac{2 - 2\sqrt{3}/3}{f(t)^2} + \frac{4}{3} + \frac{k^2}{3} + \mathcal{O}(f(t) \cdot k), \qquad (86) $$

which shows that for $t \to 0$, we have $k_* \to 0$ while $\Lambda_* \to \infty$. For $t \to \infty$, we have

$$ \Lambda = \frac{k^2}{3} + \mathcal{O}\left(\frac{1}{f(t)^2}\right), \qquad (87) $$

which implies $k_* \to 0$ and $\Lambda_* \to 0$ for $t \to \infty$. The optimal wave number and the Rayleigh quotient are presented in Fig. 12, which confirms (86) and (87).

**Remark 2** The limiting behaviors of the asymptotic Rayleigh quotient (86) and (87) actually confirms the previous assumption (67): the bifurcation time is actually a constant using the rescaled loading parameter (47), which implies $\lambda_k = \mathcal{O}(1)$.

The universal bifurcation time written using the variable $t$ can be obtained by solving $\Lambda_*(t_*) = 1$, which leads to

$$ f(t_*) \approx 1.685 \implies t_* \approx 0.9982, \quad k_* \approx 0.7946. \qquad (88) $$

Using (77), in our asymptotic analysis the optimal wave number as well as the bifurcation time are hence given by

$$ \frac{t_b}{t_c} = \frac{1}{1 - \epsilon \cdot t_*} \approx 1 + \epsilon \cdot t_*, \quad k_* = \frac{k_*}{\epsilon}. \qquad (89) $$

In Fig. 13(a), we compare the $k$-dependence of the Rayleigh quotient at fixed time, obtained respectively by the semi-analytical approach (see Fig. 7) and the asymptotic solution (84). Although the approximate Rayleigh quotient is slightly larger than the exact value, their dependence with respect to $k$ is in perfect agreement. In fact the optimal wave number in this case is 11 using the exact approach, and approximately 11.3 with the asymptotic solution.
<figure><img src="image_15.png" /><figcaption>**Fig. 13** Comparison between the semi-analytical (exact) bifurcation solutions and the asymptotic approach: (a) Rayleigh quotient $\lambda_k$ as a function of the wave number $k$ for $\ell/R = 0.05$ at $t/t_c \approx 1.073$; (b) time-dependence of the optimal Rayleigh quotient $\lambda_*$</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 14** Size effect of bifurcation with the semi-analytical approach and the asymptotic solution: (a) bifurcation time $t_*/t_c$; (b) optimal wave number $k_*$</figcaption></figure>

The time-dependence of the optimal Rayleigh quotient $\lambda_*$, evaluated with the optimal wave number $k_*$, is presented in Fig. 13(b). The previous semi-analytical result in Fig. 8 and the asymptotic solution are compared. It can be seen that the approximation is better with smaller values $\epsilon$. Due to the universal nature, the implicit $\epsilon$-dependence of $\lambda_*$ is introduced by the rescaled time parameter $t$ given by (47). Using the fact that $t \mapsto \lambda_*$ is decreasing (see Fig. 12), we deduce immediately that $\lambda_*$ decreases when $\epsilon$ becomes smaller. This result is confirmed in Fig. 13(b).

Finally, in Fig. 14 we review the size effect computed by the semi-analytical approach and the asymptotic solution. The asymptotic solution (89) predicts relatively well the $\epsilon$-dependence of the bifurcation time and the optimal wave number.

# 5 Conclusion and Perspectives

In this work the elastic-gradient damage evolution problem of a cylinder under torsion is successfully analyzed using two approaches. In the first approach, the fundamental branch and its subsequent bifurcation are solved using analytical and numerical methods, where the results obtained are valid for arbitrary values of the internal length. In the second asymptotic approach, we assume that the length scale $\epsilon = \ell/R$ between the material internal length $\ell$ and the cylinder radius $R$ is small. It is found that the simulation results obtained by the
<figure><img src="image_17.png" /><figcaption>**Fig. 15** (a) Bifurcated damage solution for $t_b < t < t_r$; (b) final fracture pattern at $t_r$ for $\ell/R = 0.05$</figcaption></figure>

semi-analytical (or analytical-numerical) approach can be formally justified by the asymptotic methods.

One of the main contribution of this work is the application of such boundary layer method for the study of gradient-damage models. The compatibility between the rescaled spatial and temporal variables is formally proved in Proposition 1. The orders of magnitude of the bifurcated solutions with respect to the small parameter $\epsilon$ are also carefully analyzed in Proposition 2. It is shown that the axisymmetric damage evolution and the bifurcation are governed by two universal functions independent of the small parameter. The $\epsilon$-dependence of the solutions can be explained largely by the introduced spatial and temporal rescaling.

1. The first universal function as illustrated in Fig. 10 governs the axisymmetric damage, and in particular its thickness evolution.

2. The second as shown in Fig. 12 gives a universal temporal evolution of the optimal Rayleigh quotient, which governs bifurcation from the fundamental branch.

The size effects of bifurcation formally derived in (89) and illustrated in Fig. 14 constitute the main results of this paper. Similar to the unidimensional problem as well as the 2d thermal shock problem, bifurcation is more likely to produce earlier for smaller values of the internal length.

In this paper, we focused on the loss of uniqueness of the fundamental branch characterized by the asymmetric damage solution. However, a more complete analysis should also involve a stability study of the bifurcated solutions. Furthermore, subsequent bifurcations may produce since it is not clear that the bifurcated *damage* pattern as presented in Fig. 7(b) may eventually evolve smoothly into fully-established *cracks*. Our numerical simulations not reported in this paper actually suggest that the bifurcated solution analyzed in this paper is stable for a finite time, see Fig. 15(a), during which the damage field is indeed periodic and agrees well with our semi-analytical and asymptotic analyses. However afterwards a second bifurcation, more abrupt than the previous one, takes place and leads to the final fracture of the cylinder, see Fig. 15(b). The numerical results will be reported in another paper in which this second bifurcation will be analyzed.

The current study is carried out using a particular set of damage constitutive functions (14). This is in general a preferred choice compared to other constitutive laws in which the elastic domain is absent. Future work could be devoted to other strongly-brittle gradient-damage materials since bifurcation behaviors could be different.

The bifurcation analysis and the boundary layer method could also be applied for three-dimensional gradient-damage problems. The physical role played by the internal length as well as the size-effects in damage evolution can thus be investigated analytically beyond the usual uni- or two-dimensional settings. As mentioned in the introduction, mixed mode I and III loading can also be considered in the future. We may hope that the experimentally ob-
served echelon cracks may correspond in fact to a bifurcation mode similar to that observed in Fig. 7(b).

**Conflict of interest** The authors declare that they have no conflict of interest.

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## Appendix. Numerical Minimization of the Rayleigh Quotient

The generalized symmetric eigenvalue problem defined by (44) and the subsequent minimization problem of the Rayleigh quotient (42) are solved numerically. Standard $\mathbb{P}_1$ finite elements are used to discretize the functions $(A, C)$ in the interval $[0, R]$. After discretization the eigenvalue problem (44) becomes

$$ K_k x_k = \lambda_k M x_k, \quad (90) $$

where

- The matrices $K_k$ and $M$ correspond to the discrete quadratic forms A and B;

- The eigenvector $x_k \neq 0$ regroups the nodal values of the functions $(A_k, C_k)$;

- The smallest eigenvalue $\lambda_k$ is a finite-element approximation of the continuous one.

In this paper spatial discretization is performed with the FEniCS framework [19]. A specific conjugate gradient procedure described in [14] implemented in the library SLEPc [33] is applied to solve the discrete eigenvalue problem (90).

The subsequent minimization problem (45) with respect to the wave number $k$ is performed through direct evaluations of the Rayleigh quotient $\lambda_k$ for $k \in \mathbb{K}$ where $\mathbb{K}$ is a subset of $\mathbb{N}^*$ which is supposed to contain the minimum wave number. The minimization procedure of the Rayleigh quotient is summarized in Algorithm 1.

<caption>**Algorithm 1** Minimizing the Rayleigh quotient (42) at a given instant $t > t_c$</caption>

1: Given a heuristic set $\mathbb{K} \subset \mathbb{N}^*$ containing the minimizing wave number $k_*$
2: **for** the wave number $k \in \mathbb{K}$ **do**
3: Solve the eigenvalue problem (90) for $(A_k, C_k)$
4: Evaluate the corresponding Rayleigh quotient $\lambda_k = \mathbb{R}(k, A_k, C_k)$
5: **end for**
6: Find the optimal Rayleigh quotient $\lambda_* = \lambda_{k_*}$ and the minimizing $k_*$ through $\lambda_* = \min_{k \in \mathbb{K}} \lambda_k$
The next objective would be to find the bifurcation time $t_*$ at which the optimal Rayleigh quotient reaches 1. Remark that the above Algorithm 1 materializes a scalar function $t \mapsto \lambda_*$ defined for $t > t_c$. A classical root-finding method is thus applied to find the corresponding instant $t_*$ when $\lambda_* = 1$.

## References

1. Abdelmoula, R., Leger, A.: Singular perturbation analysis of the buckling of circular cylindrical shells. Eur. J. Mech. A, Solids **27**(1), 706–729 (2008)

2. Benallal, A., Billardon, R., Geymonat, G.: Bifurcation and localization in rate-independent materials. Some general considerations. In: Nguyen, Q. (ed.) Bifurcation and Stability of Dissipative Systems, International Centre for Mechanical Sciences, vol. 327, pp. 1–44. Springer, Vienna (1993)

3. Benallal, A., Marigo, J.J.: Bifurcation and stability issues in gradient theories with softening. Model. Simul. Mater. Sci. Eng. **15**(1), S283 (2007)

4. Bourdin, B., Francfort, G.A., Marigo, J.J.: Numerical experiments in revisited brittle fracture. J. Mech. Phys. Solids **48**(4), 797–826 (2000)

5. Bourdin, B., Francfort, G.A., Marigo, J.J.: The variational approach to fracture. J. Elast. **91**(1–3), 5–148 (2008)

6. Bourdin, B., Marigo, J.J., Maurini, C., Sicsic, P.: Morphogenesis and propagation of complex cracks induced by thermal shocks. Phys. Rev. Lett. **112**(1), 014,301 (2014)

7. Cole, J., Kevorkian, J.: Perturbation Methods in Applied Mathematics. Springer, New York (1980)

8. Francfort, G.A., Marigo, J.J.: Revisiting brittle fracture as an energy minimization problem. J. Mech. Phys. Solids **46**(8), 1319–1342 (1998)

9. Freddi, F., Royer-Carfagni, G.: Regularized variational theories of fracture: a unified approach. J. Mech. Phys. Solids **58**(8), 1154–1174 (2010). https://doi.org/10.1016/j.jmps.2010.02.010

10. Freidrichs, K.: Perturbation of Spectra in Hilbert Spaces. Lectures in Applied Mathematics, vol. III. Am. Math. Soc., Providence (1965)

11. Goldstein, R.V., Osipenko, N.M.: Fracture structure near a longitudinal shear macrorupture. Mech. Solids **47**(5), 505–516 (2012). https://doi.org/10.3103/s0025654412050032

12. Kato, T.: Perturbation Theory for Linear Operators. Die Grundlehren der mathematischen Wissenschaften in Einzeldarstellungen, vol. 132. Springer, Berlin (1966)

13. Knauss, W.: An observation of crack propagation in anti-plane shear. Int. J. Fract. Mech. **6**(2), 183–187 (1970). https://doi.org/10.1007/bf00189825

14. Knyazev, A.V.: Toward the optimal preconditioned eigensolver: locally optimal block preconditioned conjugate gradient method. SIAM J. Sci. Comput. **23**(2), 517–541 (2001)

15. Li, T.: Gradient-Damage Modeling of Dynamic Brittle Fracture: Variational Principles and Numerical Simulations. Ph.D. thesis, Université Paris-Saclay (2016)

16. Li, T., Marigo, J.J.: Crack tip equation of motion in dynamic gradient damage models. J. Elast. **127**(1), 25–57 (2017). https://doi.org/10.1007/s10659-016-9595-0

17. Li, T., Marigo, J.J., Guilbaud, D., Potapov, S.: Gradient damage modeling of brittle fracture in an explicit dynamics context. Int. J. Numer. Methods Biomed. Eng. **108**(11), 1381–1405 (2016). https://doi.org/10.1002/nme.5262

18. Lin, B., Mear, M.E., Ravi-Chandar, K.: Criterion for initiation of cracks under mixed-mode I + III loading. Int. J. Fract. **165**(2), 175–188 (2010). https://doi.org/10.1007/s10704-010-9476-7

19. Logg, A., Mardal, K.A., Wells, G. (eds.): Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book, vol. 84. Springer, Berlin Heidelberg, Berlin (2012). https://doi.org/10.1007/978-3-642-23099-8

20. Lorentz, E., Cuvilliez, S., Kazymyrenko, K.: Modelling large crack propagation: from gradient damage to cohesive zone models. Int. J. Fract. **178**(1–2), 85–95 (2012)

21. Love, A.: A Treatise on the Mathematical Theory of Elasticity, 4th edn. Cambridge University Press, Cambridge (1927)

22. Marigo, J.J., Maurini, C., Pham, K.: An overview of the modelling of fracture by gradient damage models. Meccanica **51**(12), 3107–3128 (2016)

23. Miehe, C., Welschinger, F., Hofacker, M.: Thermodynamically consistent phase-field models of fracture: variational principles and multi-field FE implementations. Int. J. Numer. Methods Biomed. Eng. **83**(10), 1273–1311 (2010). https://doi.org/10.1002/nme.2861

24. Nayfeh, A.H.: Introduction to Perturbations Techniques. Wiley-Interscience Publication, New York (1981)
25. Nguyen, T.T., Yvonnet, J., Bornert, M., Chateau, C., Sab, K., Romani, R., Roy, R.L.: On the choice of parameters in the phase field method for simulating crack initiation with experimental validation. Int. J. Fract. **197**(2), 213–226 (2016). https://doi.org/10.1007/s10704-016-0082-1

26. Pham, K., Amor, H., Marigo, J.J., Maurini, C.: Gradient damage models and their use to approximate brittle fracture. Int. J. Damage Mech. **20**(4), 618–652 (2011)

27. Pham, K., Marigo, J.J.: From the onset of damage to rupture: construction of responses with damage localization for a general class of gradient damage models. Contin. Mech. Thermodyn. **25**(2–4), 147–171 (2013)

28. Pham, K., Marigo, J.J.: Stability of homogeneous states with gradient damage models: size effects and shape effects in the three-dimensional setting. J. Elast. **110**(1), 63–93 (2013)

29. Pham, K., Marigo, J.J., Maurini, C.: The issues of the uniqueness and the stability of the homogeneous response in uniaxial tests with gradient damage models. J. Mech. Phys. Solids **59**(6), 1163–1190 (2011)

30. Pham, K.H., Ravi-Chandar, K.: On the growth of cracks under mixed-mode I + III loading. Int. J. Fract. **199**(1), 105–134 (2016). https://doi.org/10.1007/s10704-016-0098-6

31. Pham, K.H., Ravi-Chandar, K.: The formation and growth of echelon cracks in brittle materials. Int. J. Fract. **206**(2), 229–244 (2017). https://doi.org/10.1007/s10704-017-0212-4

32. Pham, K.H., Ravi-Chandar, K., Landis, C.M.: Experimental validation of a phase-field model for fracture. Int. J. Fract. **205**(1), 83–101 (2017). https://doi.org/10.1007/s10704-017-0185-3

33. Roman, J.E., Campos, C., Romero, E., Tomas, A.: SLEPc Users Manual. Tech. Rep. DSIC-II/24/02 – Revision 3.7, D. Sistemes Informàtics I Computació, Universitat Politècnica de València, (2016)

34. Sanchez-Palencia, E.: Vibration and Coupling of Continuous Systems. Asymptotic Methods. Springer, Berlin Heidelberg New York, London Paris Tokyo (1989)

35. Sicsic, P., Marigo, J.J.: From gradient damage laws to Griffith’s theory of crack propagation. J. Elast. **113**(1), 55–74 (2013)

36. Sicsic, P., Marigo, J.J., Maurini, C.: Initiation of a periodic array of cracks in the thermal shock problem: a gradient damage modeling. J. Mech. Phys. Solids **63**, 256–284 (2014)

37. Sommer, E.: Formation of fracture ‘lances’ in glass. Eng. Fract. Mech. **1**(3), 539–546 (1969)

38. Tanné, E., Li, T., Bourdin, B., Marigo, J.J., Maurini, C.: Crack nucleation in variational phase-field models of brittle fracture. J. Mech. Phys. Solids **110**, 80–99 (2018). https://doi.org/10.1016/j.jmps.2017.09.006

39. Van Dyke, M.: Perturbation Methods in Fluid Mechanics. Parabolic Press, Stanford (1975)