# Instabilities in mechanics

Romain Lagrange

► To cite this version:

Romain Lagrange. Instabilities in mechanics. Mechanics [physics]. Institut Polytechnique de Paris, 2023. tel-04535388

HAL Id: tel-04535388

https://hal.science/tel-04535388

Submitted on 6 Apr 2024

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.

Public Domain
Institut Polytechnique de Paris

HABILITATION À DIRIGER DES RECHERCHES

# Analyse de quelques problèmes d'instabilité en mécanique

par Romain LAGRANGE

le 12 décembre 2023

## Commission d'examen :

<table><tr><td>M. Jean-Marc CHOMAZ</td><td>Président du jury</td></tr><tr><td>Mme Lydia BOUROUIBA</td><td>Rapporteur</td></tr><tr><td>M. Bruno COCHELIN</td><td>Rapporteur</td></tr><tr><td>M. Stéphane ETIENNE</td><td>Rapporteur</td></tr><tr><td>M. Xavier AMANDOLESE</td><td>Examinateur</td></tr><tr><td>Mme Marianna BRAZA</td><td>Examinatrice</td></tr><tr><td>M. Arnaud LAZARUS</td><td>Examinateur</td></tr></table>
<figure><img src="image_1.png" /></figure>
REMERCIEMENTS

Mes premiers remerciements sont adressés à Patrick Le Tallec, référent en mécanique pour l'Habilitation à Diriger des Recherches (HDR) à l'Institut Polytechnique de Paris, pour ses conseils avisés et son accompagnement précieux tout au long du processus menant à la soutenance de l'HDR.

Je tiens ensuite à exprimer ma profonde gratitude envers les membres du jury de soutenance, en particulier les rapporteurs internationaux de ce travail, Lydia Bourouiba et Stéphane Etienne. Leur engagement et leur expertise ont été cruciaux, malgré leur découverte du processus spécifique de ce diplôme français. Parmi les rapporteurs, je tiens également à remercier chaleureusement Bruno Cochelin, dont l'influence dans ma formation à la mécanique des milieux continus, il y a maintenant une vingtaine d'années à l'École Centrale de Marseille, reste inestimable. Les rapports détaillés et circ constanciés témoignent de l'efficacité et du dévouement de Lydia Bourouiba, Stéphane Etienne et Bruno Cochelin à évaluer mon travail.

Mes remerciements s'étendent ensuite à Xavier Amandolèse, Arnaud Lazarus et Marianna Braza pour avoir accepté d'examiner mon travail d'HDR. Je tiens également à exprimer ma gratitude envers Jean-Marc Chomaz pour avoir accepté de présider le jury de soutenance. Présenter mes travaux a été un réel plaisir, et au-delà des enjeux liés au diplôme, j'ai particulièrement apprécié la bienveillance des membres du jury ainsi que la fluidité de nos échanges.

Ce travail d'HDR n'aurait pas été possible sans les interactions que j'ai eues avec de nombreux chercheurs rencontrés lors de mon parcours professionnel, notamment à l'IRPHE, à l'IFPEN, au MIT et au CEA. Je leur adresse donc un remerciement général, tout en soulignant le rôle essentiel de Maria-Adela Puscas dans la genèse de cette HDR. Je souhaite également remercier Natalie Lemaintec, secrétaire ultra efficace du laboratoire DYN, dont le soutien précieux a été déterminant pour l'organisation de la soutenance d'HDR.

"Last but not least", je remercie mes parents. Sans leur soutien indéfectible et leurs efforts constants pour me fournir les meilleures conditions de travail, rien (École d'Ingénieur, doctorat, HDR...) n'aurait été possible.
Table des matières

<table><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td>1.</td><td>Dossier personnel</td><td>1</td></tr><tr><td>1.1</td><td>Research Interests</td><td>1</td></tr><tr><td>1.2</td><td>Education</td><td>1</td></tr><tr><td>1.3</td><td>Professional Experience</td><td>1</td></tr><tr><td>1.3.1</td><td>Teaching</td><td>1</td></tr><tr><td>1.3.2</td><td>Research</td><td>2</td></tr><tr><td>1.4</td><td>Students guidance</td><td>3</td></tr><tr><td>1.4.1</td><td>Ph.D. Research Scholars</td><td>3</td></tr><tr><td>1.4.2</td><td>M.Sc./M.Eng. Students</td><td>4</td></tr><tr><td>1.5</td><td>Ongoing professional activities</td><td>5</td></tr><tr><td>1.6</td><td>Institutional/Administrative responsibilities</td><td>5</td></tr><tr><td>1.7</td><td>Projects</td><td>6</td></tr><tr><td>1.8</td><td>Awards</td><td>6</td></tr><tr><td>1.9</td><td>Languages</td><td>6</td></tr><tr><td>1.10</td><td>Computer skills</td><td>6</td></tr><tr><td>1.11</td><td>Publications</td><td>6</td></tr><tr><td>1.12</td><td>Communications</td><td>10</td></tr><tr><td>2.</td><td>Sypnosis des travaux de recherche</td><td>15</td></tr><tr><td>2.1</td><td>Instabilité d'un fluide dans un cylindre en précession</td><td>16</td></tr><tr><td>2.2</td><td>Instabilités de flambage/plissement</td><td>17</td></tr><tr><td>2.2.1</td><td>Flambage d'une poutre imparfaite sur fondation non-linéaire</td><td>17</td></tr><tr><td>2.2.2</td><td>Flambage/plissement d'un film mince sur substrat circulaire</td><td>18</td></tr><tr><td>2.3</td><td>Instabilité fluide-élastique d'un faisceau de tubes</td><td>19</td></tr><tr><td>3.</td><td>Instabilité d'un fluide dans un cylindre en précession</td><td>21</td></tr><tr><td>3.1</td><td>Introduction</td><td>25</td></tr><tr><td>3.2</td><td>Formulation of the problem</td><td>26</td></tr><tr><td>3.2.1</td><td>Base flow</td><td>27</td></tr><tr><td>3.2.2</td><td>Symmetry properties</td><td>29</td></tr><tr><td>3.3</td><td>Linear stability analysis</td><td>31</td></tr><tr><td>3.3.1</td><td>Order 1: free Kelvin modes</td><td>32</td></tr><tr><td>3.3.2</td><td>Triadic resonance</td><td>32</td></tr><tr><td>3.3.3</td><td>Order Ro: slow time equations</td><td>35</td></tr><tr><td>3.3.4</td><td>Amplitude equations using the horizontal shear de-composition</td><td>36</td></tr><tr><td>3.3.5</td><td>Case of a resonant cylinder</td><td>36</td></tr></tbody></table>
<table>
  <tbody>
    <tr>
      <td>3.3.6</td>
      <td>Viscous effects</td>
      <td>37</td>
    </tr>
    <tr>
      <td>3.4</td>
      <td>Numerical application far from a resonance</td>
      <td>38</td>
    </tr>
    <tr>
      <td>3.5</td>
      <td>Conclusion</td>
      <td>41</td>
    </tr>
    <tr>
      <td>3.6</td>
      <td>Appendix</td>
      <td>42</td>
    </tr>
    <tr>
      <td>3.6.1</td>
      <td>Operators of the Euler equations, forced mode vector<br/>v<sub>j</sub>, amplitudes a<sub>j</sub><sup>V</sup> and a<sub>j</sub><sup>H</sup>, free Kelvin modes vectors<br/>v<sub>l</sub><sup>+</sup> and v<sub>l</sub><sup>-</sup></td>
      <td>42</td>
    </tr>
    <tr>
      <td>3.6.2</td>
      <td>Derivation of amplitude equations</td>
      <td>43</td>
    </tr>
    <tr>
      <td>3.6.3</td>
      <td>Computation of<br/>&lt;v<sub>1</sub>e<sup>i(ω<sub>1</sub>t+m<sub>1</sub>φ)</sup>, (D<sup>i</sup>(ωt+φ) + c.c.) (v<sub>2</sub>e<sup>i(ω<sub>2</sub>t+m<sub>2</sub>φ)</sup>)&gt;</td>
      <td>45</td>
    </tr>
    <tr>
      <td>3.6.4</td>
      <td>Computation of &lt;v<sub>1</sub>e<sup>i(ω<sub>1</sub>t+m<sub>1</sub>φ)</sup>, N(v<sub>base</sub>, v<sub>2</sub>e<sup>i(ω<sub>2</sub>t+m<sub>2</sub>φ)</sup>)&gt;</td>
      <td>45</td>
    </tr>
    <tr>
      <td>3.6.5</td>
      <td>Computation of &lt;v<sub>1</sub>e<sup>i(ω<sub>1</sub>t+m<sub>1</sub>φ)</sup>, N(v<sub>2</sub>e<sup>i(ω<sub>2</sub>t+m<sub>2</sub>φ)</sup>, v<sub>base</sub>)&gt;</td>
      <td>46</td>
    </tr>
    <tr>
      <td>3.6.6</td>
      <td>Conditions of resonance and amplitude equations</td>
      <td>46</td>
    </tr>
    <tr>
      <td>4.</td>
      <td>Instabilités de flambage/plissement</td>
      <td>51</td>
    </tr>
    <tr>
      <td>4.1</td>
      <td>Introduction</td>
      <td>55</td>
    </tr>
    <tr>
      <td>4.2</td>
      <td>Definition of the problem</td>
      <td>57</td>
    </tr>
    <tr>
      <td>4.3</td>
      <td>Analytical approach</td>
      <td>60</td>
    </tr>
    <tr>
      <td>4.3.1</td>
      <td>Kinematics, energy formulation and equations of equi-librium</td>
      <td>60</td>
    </tr>
    <tr>
      <td>4.3.2</td>
      <td>Asymptotic expansion and reactive force of the substrate</td>
      <td>62</td>
    </tr>
    <tr>
      <td>4.3.3</td>
      <td>Critical stress, wavelength and comparison to wrinkling on a planar substrate</td>
      <td>64</td>
    </tr>
    <tr>
      <td>4.4</td>
      <td>Numerical simulations</td>
      <td>65</td>
    </tr>
    <tr>
      <td>4.5</td>
      <td>Results</td>
      <td>66</td>
    </tr>
    <tr>
      <td>4.5.1</td>
      <td>Hoop stress prior to wrinkling</td>
      <td>67</td>
    </tr>
    <tr>
      <td>4.5.2</td>
      <td>Critical pressure</td>
      <td>67</td>
    </tr>
    <tr>
      <td>4.5.3</td>
      <td>Phase diagram</td>
      <td>70</td>
    </tr>
    <tr>
      <td>4.5.4</td>
      <td>Comparison of wrinkling in our curved system with that on infinite planar substrate</td>
      <td>70</td>
    </tr>
    <tr>
      <td>4.6</td>
      <td>Conclusion</td>
      <td>74</td>
    </tr>
    <tr>
      <td>4.7</td>
      <td>Appendix</td>
      <td>75</td>
    </tr>
    <tr>
      <td>4.7.1</td>
      <td>Response of the substrate</td>
      <td>75</td>
    </tr>
    <tr>
      <td>4.7.2</td>
      <td>Terms of the linear analysis of stability</td>
      <td>78</td>
    </tr>
    <tr>
      <td>4.7.3</td>
      <td>Influence of the cavity size</td>
      <td>78</td>
    </tr>
    <tr>
      <td>5.</td>
      <td>Instabilité fluide-élastique d'un faisceau de tubes</td>
      <td>85</td>
    </tr>
    <tr>
      <td>5.1</td>
      <td>Introduction</td>
      <td>91</td>
    </tr>
    <tr>
      <td>5.2</td>
      <td>Position of the problem and dimensional analysis</td>
      <td>92</td>
    </tr>
    <tr>
      <td>5.2.1</td>
      <td>Single-phase flow instability threshold</td>
      <td>94</td>
    </tr>
    <tr>
      <td>5.2.2</td>
      <td>Two-phase flow instability threshold</td>
      <td>95</td>
    </tr>
    <tr>
      <td>5.3</td>
      <td>Single-phase flow experiments</td>
      <td>96</td>
    </tr>
    <tr>
      <td>5.3.1</td>
      <td>Experimental setup and instrumentation</td>
      <td>96</td>
    </tr>
    <tr>
      <td>5.3.2</td>
      <td>Results and discussion</td>
      <td>97</td>
    </tr>
    <tr>
      <td>5.4</td>
      <td>Two-phase flow experiments</td>
      <td>101</td>
    </tr>
    <tr>
      <td>5.4.1</td>
      <td>New criterion for the instability threshold</td>
      <td>101</td>
    </tr>
    <tr>
      <td>5.4.2</td>
      <td>Experimental setup and methodology</td>
      <td>102</td>
    </tr>
    <tr>
      <td>5.4.3</td>
      <td>Results and discussion</td>
      <td>103</td>
    </tr>
  </tbody>
</table>
<table><tr><td>5.5 Conclusion</td><td>106</td></tr><tr><td>6. Perspectives de recherche</td><td>113</td></tr><tr><td>Annexes</td><td>121</td></tr><tr><td>A. Flambage d'une poutre imparfaite sur fondation non-linéaire</td><td>123</td></tr><tr><td>B. Plissement d'un film mince attaché à un substrat non-plan</td><td>137</td></tr><tr><td>C. Interaction fluide-structure entre deux cylindres concentriques</td><td>157</td></tr><tr><td>D. Interaction fluide-structure entre deux cylindres parallèles</td><td>171</td></tr></table>
# 1. DOSSIER PERSONNEL

## 1.1 Research Interests

Fluid-Structure Interactions, Rotating Flows, Elasticity, Mechanical Instabilities, Nonlinear Dynamics and Chaos

## 1.2 Education

<table><tr><td>2009</td><td><b>Ph.D.</b>, Fluid Mechanics, <i>Aix-Marseille Université</i></td></tr><tr><td>2006</td><td><b>M.Sc.</b>, Fluid Mechanics, <i>Aix-Marseille Université</i>. Magna Cum Laude</td></tr><tr><td>2006</td><td><b>M.Eng.</b>, Solid Mechanics, <i>Centrale Marseille</i>, School of Engineering</td></tr></table>

## 1.3 Professional Experience

### 1.3.1 Teaching

<table><tr><td>2022-...</td><td><b>Lecturer in Mechanics</b>, <i>E.N.P.C.</i>, Paris, France</td></tr><tr><td>Course</td><td>Introduction to continuum mechanics (Fall 2022)</td></tr><tr><td>2015-2018</td><td><b>Lecturer in Mathematics and Mechanics</b>, <i>U.P.M.C.</i>, Paris, France</td></tr><tr><td>Courses</td><td>Instabilities (Fall 2018, 2017, 2016, 2015), Mathematics for engineers (Fall 2017, 2016), Numerical analysis (Spring 2018)</td></tr><tr><td>Recitations</td><td>Continuum Mechanics (Fall 2017), Mathematics for engineers (Fall 2017, 2016), Numerical analysis (Spring 2018)</td></tr><tr><td>2017-2018</td><td><b>Lecturer in Mathematics</b>, <i>Université Paris-Saclay</i>, Orsay, France</td></tr><tr><td>Recitations</td><td>Mathematical tools (Fall 2017)</td></tr><tr><td>2017-2018</td><td><b>Lecturer in Mathematics</b>, <i>C.N.A.M.</i>, Versailles, France</td></tr><tr><td>Courses</td><td>Linear algebra and geometry (Fall 2017)</td></tr><tr><td>2012-2015</td><td><b>Instructor in Mathematics</b>, <i>M.I.T.</i>, Cambridge, U.S.A</td></tr><tr><td>Courses</td><td>Nonlinear Dynamics and Turbulence (Spring 2015, with Prof. L. Bourouiba), Nonlinear Dynamics: Chaos (Fall 2014, Fall 2013), Tensor Calculus (Jan. 2014)</td></tr></table>
<table><tr><td>Recitations</td><td>Differential Equations (Spring 2015, Spring 2014), Multivariable Calculus (Spring 2013, Fall 2012)</td></tr><tr><td>2007-2009</td><td>Lecturer in Mechanics, Université Paul Cézanne, Marseille, France</td></tr><tr><td>Hands-on works</td><td>Solid and Fluid Mechanics for Undergraduates (102 hours in two years)</td></tr></table>

1.3.2 Research

<table><tr><td>2015-<br>In progress</td><td>Research Engineer, C.E.A., Gif-sur-Yvette, France<br>Vibrations of tube bundles subject to two-phase cross-flows / Fluid-structure interactions problems<br>Main collaborators: Framatome/E.D.F, Y. Fraigneau (L.I.M.S.I.), C. Josserand (LadHyX)</td></tr><tr><td>2012-2015<br>36 months</td><td>Instructor, M.I.T., Cambridge, U.S.A.<br>Wrinkling on curved surface</td></tr><tr><td>Details</td><td>Rationalized the effect of curvature on the wrinkling problem: stability analysis, Swift-Hohenberg model for pattern selection, crystallography properties of wrinkles.<br>Collaborators: P. M. Reis, J. Dunkel, D. Terwagne, M. Brojan, F. L. Jiménez, N. Stoop</td></tr><tr><td>2010-2011<br>15 months</td><td>Post-doctoral Researcher, I.F.P.E.N., Solaize, France<br>Buckling of thin structures</td></tr><tr><td>Details</td><td>Developed mathematical models and closed-form solutions for the buckling of slender structures used in the Offshore industry<br>Advisors: D. Averbuch, M. Martinez</td></tr><tr><td>2006-2009<br>36 months</td><td>Ph.D. Research Scholar, C.N.R.S., I.R.P.H.E., Marseille, France<br>Dynamics inside a precessing fluid cylinder: theory and PIV experiments</td></tr><tr><td>Details</td><td>Analyzed the Navier-Stokes equation in rotating frame: fixed points, stability, nature of bifurcations, model for transition to turbulence. Performed PIV measurements<br>Advisors: P. Meunier, C. Eloy, F. Nadal</td></tr><tr><td>2006<br>6 months</td><td>M.Sc. Student, C.N.R.S., I.R.P.H.E., Marseille, France<br>Aeroelastic instability of a flag in uniform flow</td></tr><tr><td>Details</td><td>Stability analysis and wind tunnel experiments<br>Advisors: C. Eloy, L. Schouveiler</td></tr></table>
2005-2006 **M.Eng. Student**, Centrale Marseille, Marseille, France

90 hours *Heat transfer through a damaged material*

Details Set-up heat transfer models and performed Abaqus simulations

Advisors: S. Bourgeois, T. Desoyer

1.4 Students guidance

1.4.1 Ph.D. Research Scholars

2023-2026 Recruitment in progress, funding secured

36 months Modeling and numerical simulation of the aeroelasticity of a solar panel mounted on flexible cables

Supervised with B. Prabel and M. A. Puscas (C.E.A.)

Details The objective of this thesis is to conduct theoretical and numerical studies to determine the vibration characteristics of photovoltaic panels mounted on flexible cables. The results of these studies should ultimately allow optimal sizing of the photovoltaic panels and the system of attachment to flexible cables.

2021-2024 Maud Kocher, I.P.P., Paris, France

36 months Vibrations of the primary circuit of a nuclear power plant

Supervised with P. Moussou (E.D.F.) and D. Panunzio
(C.E.A.)

Details The objective of this thesis is to build, with the usual tools (CFD simulations, tests), an acoustic-mechanical model for the reactor vessel and its internals. The model shall describe the transmission of acoustic waves between the inlets/outlets of the vessel, taking into account the turbulence in the annular space and the oscillations of the core envelope.

2020-2023 Clément Bazin, I.P.P., Paris, France

36 months Numerical and experimental studies of two-phase flows in interaction with a tube bundle

Supervised with C. Josserand (LadHyX) and M. G. Rodio
(C.E.A.)

Details The objective of this thesis is to numerically simulate the fluid flow regimes and fluid-structure interactions in square tube bundles subject to two-phase cross-flows. The numerical results shall be corroborated experimentally and shall feed theoretical works on the modeling of random excitation forces.
1.4.2 M.Sc./M.Eng. Students

2022 Lazar Lorand, M.Eng. student at E.S.T.A.C.A., Paris, France

4 months Theoretical study of the interaction between a flexible cylindrical beam and a fluid at rest

2021 Erwan Beaugendre, M.Eng. student at Polytech, Orléans, France

6 months Numerical study of the interaction of a transverse flow and a tube bundle

2021 Erika Diaferia, M.Sc. student at E.N.S.A.M., Paris, France

6 months Numerical simulation of the interaction between an axial flow and a slender cylinder

2020 Alexandre Epalle, M.Eng. student at I.N.S.A., Lyon, France

6 months Study of numerical dry friction modelling (modeling) in structural dynamics. Application to steam generator tubes in nuclear power plants.

2020 Matheus Melo Monteverde, M.Eng. student at E.N.S.T.A., Paris, France

3 months Interaction of two cylinders immersed in a viscous fluid

2020 Daniel Zegarra, M.Eng. student at Sup. Galilée, Paris, France

3 months Numerical simulation with TrioCFD of the hydrodynamic interaction of two cylinders subjected to small oscillations

2020 Domenico Panunzio, M.Sc. student at E.N.S.A.M., Paris, France

6 months Numerical simulation of a single-phase flow around a vibrating cylindrical tube

2019 Salvatore Ambriola, M.Sc. student at E.N.S.A.M., Paris, France

6 months Numerical simulation of the turbulent force acting on a cylinder subjected to a two-phase flow

2013-2014 Stanislav Y. Tsitkov, undergraduate student (UROP) at M.I.T., Cambridge, U.S.A.

12 months Transition to chaos in a precessing fluid cylinder

2012 Alice Nasto, graduate student at M.I.T., Cambridge, U.S.A.

2 months Fluid-Structure Interactions: from Fluttering Flags to Flapping Wings
## 1.5 Ongoing professional activities

**Examiner** Since 2020. CPGE maths program, Ecole Polytechnique Féminine, France

**Examiner** Since 2017. Bachelor program, Ecole Polytechnique, France

**Grader** 2019-2020. Centrale Supélec national exam of mathematics, *French Grandes Ecoles*, France

**Peer Review** Journal of Fluids and Structures, Physical Review Letters, Physics of Fluids, Chaos, Physical Review Fluids, European Journal of Mechanics - B/Fluids, Part C: Journal of Mechanical Engineering Science, Journal of Pressure Vessel Technology

## 1.6 Institutional/Administrative responsibilities

**2023-2026** Responsible for the 'Fiche d'action' titled 'FLOW'. This intern document, produced annually, defines the orientations of specific research works conducted at C.E.A., anticipating human, material needs, and associated funding

**2020-2023** Responsible for the 'Fiche d'action' titled 'Thèse diphasique', related to the Ph.D. work of C. Bazin

**Since 2020** Responsible for the 'Fiche d'action' titled 'F.F.A.T.T.', related to the theoretical work on vibrations of tube bundles, conducted for the institute C.E.A./Framatome

**2022** Member of the organizing committee, Flow-Induced Vibration conference, Paris-Saclay, 2022

**2018-2020** Responsible for the 'Fiche d'action' titled 'F.S.A./F.F.E.', related to the theoretical work on vibrations of tube bundles, conducted for the institute C.E.A./Framatome/E.D.F.

**2020** Representative of the Department of Systems and Structures Modeling (french acronym: D.M.2S.). Presentation of research activities on vibrations of tube bundles, for the evaluation of D.M.2S. by the scientific council of the Energy Division at Saclay (french acronym: D.E.S.)

**2019** Representative of D.M.2S. for its evaluation by the H.C.E.R.E.S. committee. Interview and discussion with the committee

**2019** Representative of the Laboratory of Dynamics at the general assembly of the Mechanical and Thermal Studies Service (french acronym: S.E.M.T.). Presentation of research activities on vibrations of tube bundles
**2016-2019** Member of the Unit Council (french acronym: C.U.) of S.E.M.T. This advisory council aims to regulate the life of the research unit

## 1.7 Projects

**Partner** GO-VIKING European project. Modeling and simulation of vibrations in steam generators and fuel assemblies

## 1.8 Awards

**2017** Nominated Specialist in Mechanics and Dynamics, C.E.A.

**2009** Elected Main Presentation of the Vortex Session, *Congrès Français de Mécanique*

**2008** Best Orator Award, C.E.A./C.E.S.T.A.

**2006** Best Poster Award, C.E.A./C.E.S.T.A.

## 1.9 Languages

**French** Native

**English** Fluent, B2 level, T.O.E.I.C.

**Spanish** Fluent, B2 level, D.E.L.E.

## 1.10 Computer skills

**Programming** Matlab, Mathematica, Maple, basic knowledge in C/C++ language and Abaqus

**Publishing** Latex, Microsoft Office

## 1.11 Publications

Articles available here

**Journal Publications** C. Bazin, M. G. Rodio, R. Lagrange, C. Josserand. Turbulence model effects on the two-phase all-regimes models prediction. *Submitted to Int. J. Multiph. Flow*, 2023

R. Lagrange, D. Panunzio, P. Piteau, X. Delaune, J. Antunes. A new criterion for the instability threshold of a square tube bundle subject to an air-water cross-flow. *J. Fluids Struct.*, vol. 122, 103980, 2023
P. Piteau, D. Panunzio, R. Lagrange, X. Delaune, J. Antunes. Experimental investigation of in-flow fluidelastic instability for rotated triangular tube bundles subjected to single-phase and two-phase transverse flows. *Submitted to J. Fluids Struct., special issue*, 2023

M. A. Puscas, R. Lagrange. Interaction of two cylinders immersed in a viscous fluid. On the effect of moderate Keulegan-Carpenter numbers on the fluid forces. *Eur. J. Mech. B Fluids*, vol. 101, 106-117, 2023

R. Lagrange, L. Lorand, M. A. Puscas. Forced beam vibrations of coaxial cylinders separated by a fluid gap of arbitrary size. Inviscid theory and numerical assessment of the fluid forces. *J. Fluids Struct.*, vol. 120, 103899, 2023

R. Lagrange, M. A. Puscas. Viscous theory for the vibrations of coaxial cylinders. Analytical formulas for the fluid forces and the modal added coefficients. *J. Appl. Mech.*, JAM-22-1411, 2023

R. Lagrange, M. A. Puscas, P. Piteau, X. Delaune, J. Antunes. Modal added-mass matrix of an elongated flexible cylinder immersed in a narrow annular fluid, considering various boundary conditions. New theoretical results and numerical validation. *J. Fluids Struct.*, vol. 114, 103754, 2022

D. Panunzio, M. A. Puscas, R. Lagrange. FSI - Vibrations of immersed cylinders. Simulations with the engineering open source code TrioCFD. Test cases and experimental comparisons. *C. R. Mecanique*, vol. 350, 451-476, 2022

R. Lagrange, M. A. Puscas. Hydrodynamic interaction between two flexible finite length coaxial cylinders: new theoretical formulation and numerical validation. *J. Appl. Mech.*, 89(8): 081006, 2022

R. Lagrange, Y. Fraigneau. New estimations of the added mass and damping of two cylinders vibrating in a viscous fluid, from theoretical and numerical approaches. *J. Fluids Struct.*, vol. 92, 102818, 2020

R. Lagrange, X. Delaune, P. Piteau, L. Borsoi, J. Antunes. A new analytical approach for modeling the added mass and hydrodynamic interaction of two cylinders subjected to large motions in a potential stagnant fluid. *J. Fluids Struct.*, vol. 77, 102-114, 2018

F.L. Jiménez, N. Stoop, R. Lagrange, J. Dunkel, P. M. Reis. Curvature-controlled defect localization in elastic surface crystals. *Phys. Rev. Lett.*, vol. 116, 104301, 2016
R. Lagrange, F.L. Jiménez, D. Terwagne, M. Brojan, P. M. Reis. From wrinkling to global buckling of a ring on a curved substrate. *J. Mech. Phys. Solids*, vol. 89, 77-95, 2016

R. Lagrange, P. Meunier, C. Eloy. Triadic instability of a non-resonant precessing fluid cylinder. *C. R. Mecanique*, vol. 344, 418-433, 2016

N. Stoop, R. Lagrange, D. Terwagne, P. M. Reis, J. Dunkel. Curvature-induced symmetry breaking of elastic surface patterns. *Nature Mater.*, vol. 14, 337-342, 2015

M. Brojan, D. Terwagne, R. Lagrange, P.M. Reis. Wrinkling crystallography on spherical surfaces. *Proc. Natl. Acad. Sci. U.S.A.*, vol. 112(1), 14-19, 2015

R. Lagrange. Limit point buckling of a finite beam on a nonlinear foundation. *Theor. Appl. Mech. Lett.*, vol. 4, 031002, 2014

R. Lagrange. Compression-induced stiffness in the buckling of a one fiber composite. *Theor. Appl. Mech. Lett.*, vol. 3, 061001, 2013

R. Lagrange, D. Averbuch. Solution methods for the growth of a repeating imperfection in the line of a strut on a nonlinear foundation. *Int. J. Mech. Sci.*, vol. 63(1), 48-58, 2012

R. Lagrange, P. Meunier, C. Eloy, F. Nadal. Precessional instability of a fluid cylinder. *J. Fluid Mech.*, vol. 666, 104-145, 2011

J.P. Lambelin, F. Nadal, R. Lagrange, A. Sarthou. Non-resonant viscous theory for the stability of a fluid-filled gyroscope. *J. Fluid Mech.*, vol. 639, 167-194, 2009

R. Lagrange, P. Meunier, C. Eloy, F. Nadal. Dynamics of a fluid inside a precessing cylinder. *Mec. Ind.*, EDP Sciences, vol. 10, 187-194, 2009

R. Lagrange, C. Eloy, F. Nadal, P. Meunier. Instability of a fluid inside a precessing cylinder. *Phys. Fluids*, vol. 20, 081701, 2008

P. Meunier, C. Eloy, R. Lagrange, F. Nadal. A rotating fluid subject to weak precession. *J. Fluid Mech.*, vol. 599, 405-440, 2008

C. Eloy, R. Lagrange, C. Souillez, L. Schouveiler. Aeroelastic instability of cantilevered flexible plates in uniform flow. *J. Fluid Mech.*, vol. 611, 97-106, 2008

M. Kocher, S. Benhamadouche, P. Moussou, A. Joly, D. Panunzio, R. Lagrange. Fluid structure interaction in the annular gap of a nuclear reactor vessel: design of a test rig. *ERCOFTAC 2023 Symposium on Multiphysics critical flow dynamics involving moving/deformable structures with design applications*, Toulouse, 2023

Conference publications
J. Antunes, P. Piteau, X. Delaune, R. Lagrange. Spectral extrapolation of frequency-dependent fluidelastic coupling coefficients from causality enforcing. *Proceedings of the 26th International Conference on Structural Mechanics in Reactor Technology*, Berlin, 2022

R. Lagrange, D. Panunzio, P. Piteau, X. Delaune, J. Antunes. A new criterion for the instability threshold of a square tube bundle subject to an air-water cross-flow. *Proceedings of the Flow Induced Vibration conference*, Paris-Saclay, 2022

D. Panunzio, R. Lagrange, P. Piteau, X. Delaune, J. Antunes. Experimental investigation of out-of-flow fluidelastic instability for rotated triangular bundles subjected to single-phase and two-phase transverse flow. *Proceedings of the Flow Induced Vibration conference*, Paris-Saclay, 2022

J. Antunes, P. Piteau, X. Delaune, D. Panunzio, R. Lagrange. Experimental investigation of in-flow fluidelastic instability for rotated triangular bundles subjected to single-phase and two-phase transverse flow. *Proceedings of the Flow Induced Vibration conference*, Paris-Saclay, 2022

P. Moussou, M. Kocher, D. Panunzio, R. Lagrange, A. Joly. Fluid structure interaction in a pressure vessel: a multipole approach for acoustic analysis. *Proceedings of the Flow Induced Vibration conference*, Paris-Saclay, 2022

M. Kocher, P. Moussou, D. Panunzio, R. Lagrange, A. Joly. Fluid structure interaction in a pressure vessel: turbulent forcing. *Proceedings of the Flow Induced Vibration conference*, Paris-Saclay, 2022

J. Antunes, P. Piteau, X. Delaune, R. Lagrange. Spectral extrapolation of frequency-dependent fluidelastic coupling coefficients from causality enforcing. *Proceedings of the SMIRT-26 conference*, Berlin-Potsdam, 2022

R. Lagrange, P. Piteau, X. Delaune, J. Antunes. Fluid-elastic coefficients in single phase cross flow: dimensional analysis, direct and indirect experimental methods. *Proceedings of the ASME PVP conference*, San-Antonio, 93984, 2019

R. Lagrange, P. Meunier, C. Eloy, F. Nadal. Instability of a fluid inside a precessing cylinder. 12e Rencontre du non-linéaire (Eds. C. Josserand, M. Lefranc, C. Letellier), *Non-Linéaire Publications*, Paris, 103-108, 2009

R. Lagrange, P. Meunier, C. Eloy, F. Nadal. Instability of a fluid inside a precessing cylinder. 11e Rencontre du non-linéaire, *Non-Linéaire Publications*, Paris, 103-108, 2008
Technical
notes

R. Lagrange. Analyse dimensionnelle des forces de couplage. Ecoulement diphasique eau-air, portance, pas carré. *NT/2021-67624/A*, C.E.A., 2021

R. Lagrange, D. Panunzio. Bilan sur l'adimensionnement des forces fluides diphasiques en pas carré. *NT/2021-69276/A*, C.E.A., 2021

R. Lagrange. Analyse dimensionnelle des forces aléatoires en portance, en écoulement diphasique eau-air et en faisceau à pas carré. *NT/2020-66059/A*, C.E.A., 2020

R. Lagrange, P. Piteau, M. Bellanger. Instabilité fluide-élastique en faisceau à pas carré. Proposition d'une nouvelle loi pour le seuil de stabilité en écoulement diphasique. *NT/2019-65374/A*, C.E.A., 2019

R. Lagrange. Analyse dimensionnelle des forces de couplage en écoulement diphasique eau-air. *NT/2018-63940/A*, C.E.A., 2018

R. Lagrange. Potential theory for the fluid-elastic instability: application to the two cylinders problem. *NT/2018-62987/A*, C.E.A., 2018

R. Lagrange. Analyse comparative du terme d'amortissement de la force de couplage mesuré par les méthodes directes et indirectes, en monophasique. *NT/2018-62652/A*, C.E.A., 2018

R. Lagrange. Dimensional analysis of the fluid-elastic coefficients in monophasic fluid flow. *Technical note, NT/2017-62195/A*, C.E.A., 2017

R. Lagrange. Potential-flow theory for the added mass and fluid force of two circular cylinders. *NT/16-018/A*, C.E.A., 2016

R. Lagrange. Flambage d'armure : modèles analytiques. *Technical note, I.F.P.E.N.*, 70 pages, 2011

Course hand-out
R. Lagrange. Introduction to continuum mechanics. *E.N.P.C.*, 2022. Lecture notes available here

Book chapter
R. Lagrange, P. Piteau, X. Delaune, D. Panunzio, B. Prabel, M. A. Puscas, M. G. Rodio, G. Ricciardi, V. Faucher. *Chapitre 2, section 2. Vibrations sous écoulement. Monographie CEA, Mécanique des structures*, Publication in progress, 2023

## 1.12 Communications

2023
M. A. Puscas, R. Lagrange. Vibrations de deux cylindres flexibles coaxiaux dans un fluide visqueux. Formules analytiques pour les forces fluides et les coefficients modaux ajoutés. *SMAI 2023*, Le Gosier, Guadeloupe, France
M. A. Puscas, R. Lagrange. Vibrations of two coaxial flexible cylinders in a viscous fluid. 22nd Computational Fluids Conference, Cannes, France

M. Kocher, P. Moussou, A. Joly, S. Benhamadouche, P. Panunzio, R. Lagrange. Fluid-structure interaction in the annular gap of a nuclear reactor vessel. ERCOFTAC 2023 Symposium on Multiphysics critical flow dynamics involving moving/deformable structures with design applications, Toulouse, France

C. Bazin, M. G. Rodio, R. Lagrange, C. Josserand. On the simulation of two-phase cross flow in a rod bundle guided by experimental forces spectra. 11th International Conference on Multiphase Flow, Kobe, Japan

C. Bazin, M. G. Rodio, R. Lagrange, C. Josserand. On the effect of turbulence models in the numerical simulation of bubbly flows in rectangular cross channels. 14th European Fluid Mechanics Conference (EFMC14), Athens, Greece

R. Lagrange, D. Panunzio, P. Piteau, X. Delaune, J. Antunes. A new criterion for the instability threshold of a square tube bundle subject to an air-water cross-flow. Flow Induced Vibration conference, Paris-Saclay, France

D. Panunzio, R. Lagrange, P. Piteau, X. Delaune, J. Antunes. Experimental investigation of out-of-flow fluidelastic instability for rotated triangular bundles subjected to single-phase and two-phase transverse flow. Flow Induced Vibration conference, Paris-Saclay, 2022

J. Antunes, P. Piteau, X. Delaune, D. Panunzio, R. Lagrange. Experimental investigation of in-flow fluidelastic instability for rotated triangular bundles subjected to single-phase and two-phase transverse flow. Flow Induced Vibration conference, Paris-Saclay, 2022

P. Moussou, M. Kocher, D. Panunzio, R. Lagrange, A. Joly. Fluid structure interaction in a pressure vessel: a multipole approach for acoustic analysis. Flow Induced Vibration conference, Paris-Saclay, 2022

M. Kocher, P. Moussou, D. Panunzio, R. Lagrange, A. Joly. Fluid structure interaction in a pressure vessel: turbulent forcing. Flow Induced Vibration conference, Paris-Saclay, 2022

J. Antunes, P. Piteau, X. Delaune, R. Lagrange. Spectral extrapolation of frequency-dependent fluidelastic coupling coefficients from causality enforcing. SMIRT-26 conference, Berlin-Potsdam, 2022

M. A. Puscas, R. Lagrange. Fluid-structure interaction problem of a flexible cylinder in a narrow coaxial cylindrical duct. WCCM-APCOM conference, Yokohama, Japan
M. A. Puscas, R. Lagrange. Fluid-structure interaction problem of two coaxial vibrating flexible cylinders separated by a thin layer of fluid. *8th European Congress on Computational Methods in Applied Sciences and Engineering*, Oslo, Norway

M. A. Puscas, R. Lagrange. Scaling law for large motions of two cylinders vibrating in a viscous Fluid. *International Conference on Wave Mechanics and Vibrations (10th WMVC)*, Lisbon, Portugal

2021 M. A. Puscas, R. Lagrange. Numerical simulation of the hydrodynamic interaction of two cylinders vibrating in a viscous fluid. *Pan American Congress on Computational Mechanics 2021*, Rio, Brazil

2019 R. Lagrange, P. Piteau, X. Delaune, J. Antunes. Fluid-elastic coefficients in single phase cross flow: dimensional analysis, direct and indirect experimental methods. *P.V.P. Meeting*, San Antonio, U.S.A.

P. Piteau, X. Delaune, J. Antunes, L. Borsoi, R. Lagrange. Experiments on the flow-induced streamwise dynamics of a rotated triangular tube bundle subjected to single-phase and two-phase transverse flows. *P.V.P. Meeting*, San Antonio, U.S.A.

2017 R. Lagrange. Quelques exemples d'instabilité en mécanique. *L.I.M.S.I. Seminar*, Orsay, France

2015 R. Lagrange, P. Meunier, C. Eloy, F. Nadal. Dynamique d'un fluide dans un cylindre en précession. *I.M.F.T. Seminar*, Toulouse, France

N. Stoop, R. Lagrange, D. Terwagne, P. M. Reis, J. Dunkel. Curvature-induced symmetry breaking selects elastic wrinkling patterns. *Annual Meeting of the A.P.S. March Meeting*, San Antonio, U.S.A.

R. Lagrange, P. Meunier, C. Eloy, F. Nadal. Dynamique d'un fluide dans un cylindre en précession. *L.M.F.A. Seminar*, Lyon, France

R. Lagrange, P. Meunier, C. Eloy, F. Nadal. Dynamique d'un fluide dans un cylindre en précession. *Experimental Physics and Modelling Seminar*, E.N.S. Lyon, France

2014 F. L. Jimenez, R. Lagrange, P. M. Reis. Buckling of a thin film on a curved soft substrate. *ASME 2014 International Mechanical Engineering Congress and Exposition*, Montreal, Canada

P. M. Reis, D. Terwagne, M. Brojan, R. Lagrange. Smart Morphable Surfaces for Aerodynamic Drag Control. *17th U. S. National Congress on Theoretical and Applied Mechanics*, East Lansing, U.S.A.
P. M. Reis, M. Brojan, D. Terwagne, R. Lagrange. Wrinkling Crystallography on Curved Surfaces. *Annual Meeting of the A.P.S. March Meeting*, Denver, U.S.A.

2013 R. Lagrange. Wrinkling on a curved surface. *Fluid lab group meeting*, Department of Applied Mathematics, Cambridge, U.S.A.

R. Lagrange. Wrinkling on a curved surface. *E.G.S. lab group meeting*, Department of Civil Engineering, Cambridge, U.S.A.

R. Lagrange, D. Averbuch. Buckling of an imperfect beam on a nonlinear foundation. *3th Annual Meeting of the A.P.S. March Meeting*, Baltimore, U.S.A.

2012 R. Lagrange. Imperfect beam on foundation. *New England Workshop on the Mechanics of Materials and Structures, Brown University*, Providence, U.S.A.

R. Lagrange, P. Meunier, C. Eloy, F. Nadal. Fluid flow inside a precessing cylinder. *Fluid lab group meeting*, Department of Applied Mathematics, Cambridge, U.S.A.

R. Lagrange, D. Averbuch. Buckling of an imperfect beam on a nonlinear foundation. *E.G.S. lab group meeting*, Department of Civil Engineering, Cambridge, U.S.A.

2011 R. Lagrange, P. Meunier, C. Eloy, F. Nadal. Ecoulement d'un fluide dans un cylindre en précession. *Séminaire I.F.P.E.N.*, Solaize, France

R. Lagrange, D. Averbuch, M. Martinez. Flambage d'armure : modèles analytiques. *I.F.P.E.N.-Technip meeting*, Solaize, France

R. Lagrange, D. Averbuch. Explicit solutions for the buckling of an imperfect strut on a nonlinear elastic foundation. *Engineering Mechanics Institute Conference*, Boston, U.S.A.

R. Lagrange, P. Meunier, C. Eloy, F. Nadal. Fluid flow inside a precessing cylinder. *Engineering Mechanics Institute Conference*, Boston, U.S.A.

R. Lagrange, D. Averbuch, M. Martinez. Modèle de cinématique exacte pour le flambage des fils d'armure. *I.F.P.E.N.-Technip meeting*, Paris, France

R. Lagrange, D. Averbuch, M. Martinez. Comportement non-linéaire d'un fil d'armure sous compression puis flexion. *I.F.P.E.N.-Technip meeting*, Paris, France

2010 R. Lagrange, P. Meunier, C. Eloy, F. Nadal. Nonlinear dynamics of a fluid inside a precessing cylinder. *63th Annual Meeting of the A.P.S. Division of Fluid Dynamics*, Long Beach, U.S.A.

R. Lagrange, D. Averbuch, M. Martinez. Post-flambage d'une poutre avec défaut. *I.F.P.E.N.-Technip meeting*, Paris, France
R. Lagrange, P. Meunier, C. Eloy, F. Nadal. Dynamique d'un fluide dans un cylindre en précession. *Invited Presentation at O.N.E.R.A.*, Meudon, France

2009

R. Lagrange, P. Meunier, C. Eloy, F. Nadal. Dynamics of a fluid flow inside a precessing cylinder. *62th Annual Meeting of the A.P.S. Division of Fluid Dynamics*, Minneapolis, U.S.A.

R. Lagrange, P. Meunier, C. Eloy, F. Nadal. From laminar to turbulent flow in a precessing cylinder. *Workshop "Waves and Instabilities in Geophysical and Astrophysical Flows"*, Porquerolles, France

R. Lagrange, P. Meunier, C. Eloy, F. Nadal. Dynamics of a fluid inside a precessing cylinder. *19e Congrès Français de Mécanique*, Marseille, France

2008

P. Meunier, R. Lagrange, C. Eloy, F. Nadal. Chaotic dynamics of a fluid inside a precessing cylinder. *12e Rencontre du non-linéaire*, Paris, France

R. Lagrange, P. Meunier, C. Eloy, F. Nadal. Dynamics of a fluid inside a precessing cylinder. *IUTAM symposium "150 Years of Vortex Dynamics"*, Copenhagen, Denmark

P. Meunier, R. Lagrange, C. Eloy, F. Nadal. Dynamics of a fluid inside a precessing cylinder. *7th Euromech Fluid Mechanics Conference*, Manchester, U.K.

R. Lagrange, P. Meunier, C. Eloy, F. Nadal. Dynamics of a fluid inside a precessing cylinder. *22nd International Congress of Theoretical and Applied Mechanics (I.C.T.A.M.)*, Adelaide, Australia

R. Lagrange, P. Meunier, C. Eloy, F. Nadal. Dynamique d'un fluide dans un cylindre en précession. *I.R.P.H.E. Seminar*, Marseille, France

2007

R. Lagrange, P. Meunier, C. Eloy, F. Nadal. Instabilité d'un fluide dans un cylindre en précession. *11e Rencontre du non-linéaire*, Paris, France

R. Lagrange, P. Meunier, C. Eloy, F. Nadal. Stabilité d'un fluide dans un cylindre en précession. *18e Congrès Français de Mécanique*, Grenoble, France

C. Eloy, P. Meunier, R. Lagrange, F. Nadal. Rotating fluid cylinder submitted to a weak precession. *60th Annual Meeting of the A.P.S. Division of Fluid Dynamics*, Salt Lake City, U.S.A.

R. Lagrange, P. Meunier, C. Eloy, F. Nadal. Stability of a fluid flow inside a precessing cylinder. *3rd European S.C.A.T. Workshop and Summer School*, Porquerolles, France

C. Eloy, L. Schouveiler, R. Lagrange. Flag flutter. *Chaos complexity, and transport C.C.T.07*, Marseille, France
## 2. SYPNOSIS DES TRAVAUX DE RECHERCHE

Mon travail de recherche est consacré à l'analyse de phénomènes non linéaires (en particulier les instabilités, que j'ai commencé à étudier dès mon Master, avec l'instabilité de flottement d'un drapeau soumis à un écoulement parallèle à sa surface) dans le domaine de la mécanique des fluides et des structures. Mon objectif est de répondre à des problèmes pratiques tels que la dynamique des tubes en U dans les générateurs de vapeur, ou à des problèmes plus fondamentaux tels que la relation entre la géométrie et les instabilités dans les systèmes mécaniques (par exemple les structures élastiques minces). Que les problèmes soient pratiques ou fondamentaux, le fil conducteur de mon travail reste l'utilisation de techniques de modélisation pour comprendre des phénomènes complexes.

Je mène principalement des études théoriques en étroite collaboration avec des expérimentateurs, des experts en simulation numérique et des étudiants (stagiaires et 2 doctorants) que j'encadre. Un objectif qui me tient particulièrement à cœur est de déduire des formulations analytiques prédictives, en considérant des modèles simplifiés de problèmes de mécanique intrinsèquement complexes. En général, ces modèles simplifiés se réduisent à un ensemble d'équations aux dérivées partielles comportant des termes non linéaires et plusieurs paramètres physiques. Les non-linéarités peuvent provenir de diverses sources, notamment des charges mécaniques externes, des changements de géométrie ou de position du système modélisé (non-linéarités géométriques) ou, plus rarement dans mes travaux, de comportements spécifiques des matériaux (non-linéarités matérielles).

En général, les termes non linéaires produisent un ensemble riche et complexe de solutions dont l'existence et la stabilité dépendent de nombres sans dimension. Pour découvrir ces solutions, déchiffrer leurs propriétés dynamiques, élucider les mécanismes de base et les paramètres de gouvernance, j'utilise des techniques classiques de mathématiques appliquées et de mécanique : analyse dimensionnelle, minimisation de fonctionnelles, analyses linéaires de stabilité et faiblement non-linéaires, expansions asymptotiques, théorie des écoulements potentiels, décomposition de Helmholtz des écoulements...

Pour illustrer mes contributions à l'analyse des phénomènes non linéaires en mécanique, j'ai choisi de présenter dans ce mémoire les résultats que j'ai obtenus dans les travaux de recherche qui suivent.
## 2.1 Instabilité d'un fluide dans un cylindre en précession

Les écoulements forcés par un mouvement de précession deviennent instables lorsque les vitesses de rotation et/ou l'angle de précession dépassent un seuil critique<fnref n="1" />. Pendant ma thèse (I.R.P.H.E., 2006-2009), j'ai effectué une analyse linéaire de stabilité de ces écoulements, en me concentrant sur le cas particulier d'un cylindre résonnant (cylindre avec un rapport d'aspect spécifique pour lequel l'écoulement de base se réduit à un unique mode hydrodynamique forcé). J'ai continué ce travail lors de mon séjour postdoctoal au M.I.T. (2012-2015), afin de développer une nouvelle théorie, valable quel que soit le rapport d'aspect du cylindre. J'ai montré que l'écoulement de base, induit par la précession du cylindre, est une superposition d'un écoulement de cisaillement (vertical ou horizontal) et d'une somme infinie de modes forcéés. J'ai ensuite effectué une analyse linéaire de stabilité de cet écoulement de base en considérant sa résonance triadique avec deux modes de Kelvin libres (solutions homogènes des équations de Navier-Stokes linéarisées). J'ai dérivé les équations d'amplitude des modes de Kelvin libres et obtenu une expression analytique du seuil d'instabilité et du taux de croissance. Ensuite, j'ai étudié la transition vers la turbulence en supervisant pendant un an un étudiant de premier cycle au M.I.T. Ensemble, nous avons effectué des simulations MATLAB des équations faiblement non linéaires que j'avais établies théoriquement. Les simulations ont montré que l'écoulement instable stationnaire devient périodique par le biais d'une bifurcation de Hopf supercritique. En augmentant davantage le nombre de Reynolds, nous avons montré que l'écoulement périodique subit de nombreuses bifurcations de doublement de période ainsi qu'un phénomène d'intermittence (coexistence de fenêtres périodiques et de chaos avec un attracteur étrange). Certains de ces résultats sont maintenant utilisés par P. Meunier (ancien directeur de thèse à l'I.R.P.H.E.) et R. Manasseh (Swinburne University of Technology) dans le cadre du développement d'un mélangeur sans pales, utilisant la précession pour assurer un mélange puissant pour un faible cisaillement (brevets FR1651703, US11311848). Mes travaux de recherche sur l'instabilité de précession ont été publiés dans les journaux suivants :

* R. Lagrange, P. Meunier, C. Eloy. Triadic instability of a non-resonant precessing fluid cylinder. C. R. Mecanique, 344, 418-433, 2016.

* R. Lagrange, P. Meunier, C. Eloy, F. Nadal. Precessional instability of a fluid cylinder. J. Fluid Mech., vol. 666, 104-145, 2011.

* J.P. Lambelin, F. Nadal, R. Lagrange, A. Sarthou. Non resonant viscous theory for the stability of a fluid-filled gyroscope. J. Fluid Mech., vol. 639, 167-194, 2009. PDF

* R. Lagrange, P. Meunier, C. Eloy, F. Nadal. Dynamics of a fluid inside a precessing cylinder. Mec. Ind., EDP Sciences, 10, 187-194, 2009.

<footnote n="1">¹ Vidéo I.R.P.H.E. de l'instabilité d'un fluide dans un cylindre en précession, visualisation avec des particules (kalliroscope) réfléchissantes : ici.</footnote>
* R. Lagrange, C. Eloy, F. Nadal, P. Meunier. Instability of a fluid inside a precessing cylinder. Phys. Fluids, vol. 20, 081701, 2008.

* P. Meunier, C. Eloy, R. Lagrange, F. Nadal. A rotating fluid subject to weak precession. J. Fluid Mech., vol.599, 405-440, 2008.

## 2.2 Instabilités de flambage/plissement

Les deux instabilités, flambage et plissement, appartiennent à la catégorie générique des instabilités de structures élastiques, élancées et/ou minces, soumises à des contraintes de compression ou de traction (i.e. compression transverse via l'effet Poisson). Lorsque la contrainte devient trop importante, la structure peut subir une instabilité se traduisant par une déformation hors-plan. Bien que la dénomination ne soit pas universelle, on parle généralement de flambage (buckling en Anglais) dans le cas d'une déformation globale de la structure. Le plissement (wrinkling en Anglais) est en général associé à un mode de flambage localisé (déformation exhibant une petite longueur d'onde), se traduisant par l'apparition de rides ou de plis au niveau de la surface de la structure.

### 2.2.1 Flambage d'une poutre imparfaite sur fondation non-linéaire

Dans le cadre d'une collaboration avec l'industrie parapétrolière (Technip), je me suis intéressé, avec D. Averbuch (I.F.P.E.N., 2010-2011), à la déformation des fils d'armure (tensile armours/wires en Anglais) qui assurent la tenue mécanique des risers flexibles aux efforts axiaux. Nous avons considéré le problème modèle d'une poutre imparfaite (défaut de courbure), soumise à un effort de compression axial et à une force de rappel transverse. L'imperfection est introduite pour prendre en compte la courbure géodésique du fil d'armure. La force de rappel (type bilinéaire ou régularisée) modélise quant à elle l'interaction non-linéaire d'un fil d'armure avec son support (type Winkler). La minimisation d'une énergie potentielle (incluant une énergie de déformation de flexion et les travaux des efforts extérieurs) conduit à une équation différentielle d'ordre 4 pour le déplacement transverse, non-linéaire, avec un terme de forçage directement lié à l'amplitude de l'imperfection. Pour résoudre cette équation différentielle (avec des conditions aux limites de type appui-appui) et déterminer les états d'équilibre statiques du système, nous avons proposé deux méthodes théoriques. La première méthode, applicable uniquement dans le cas d'une force de rappel de type bilinéaire, est basée sur la recherche d'une solution par morceaux de l'équation d'équilibre. Cette méthode, bien que contraignante car nécessitant des conditions de raccordement (continuité du déplacement, de la tangeante, de la courbure et de l'effort tranchant) en deux points de la structure, conduit à une solution exacte du problème d'équilibre. La deuxième méthode, valable quel que soit le type de force de rappel (bilinéaire ou régularisée), est basée sur une procédure de Galerkin (recherche d'une solution ayant la même forme que l'imperfection) et la minimisation d'un résidu. Cette approche, simple à mettre en place, produit une relation explicite entre la charge de compression et l'amplitude du déplacement, mais conduit à une solution approximative du problème
d'équilibre. Quelle que soit la méthode utilisée, nous avons montré que le système est sensible à l'imperfection et que la déformée est une amplification du défaut. De plus, nous avons montré l'existence d'une taille critique d'imperfection gouvernant l'existence d'une charge limite (i.e. un point limite) dans les courbes d'équilibre du système. Une expression analytique de cette imperfection critique et des lois d'évolution pour la charge limite ont été proposées et vérifiées numériquement (outil O.D.E. 45 de Matlab). Les travaux de recherche sur l'instabilité de flambage ont été publiés dans les journaux suivants :

* R. Lagrange. Limit point buckling of a finite beam on a nonlinear foundation. Theor. Appl. Mech. Lett., vol. 4, 031001, 2014.

* R. Lagrange. Compression-induced stiffness in the buckling of a one fiber composite. Theor. Appl. Mech. Lett., vol. 3, 061001, 2013.

* R. Lagrange, D. Averbuch. Solution methods for the growth of a repeating imperfection in the line of a strut on a nonlinear foundation. Int. J. Mech. Sci., vol. 63(1), 48-58, 2012.

### 2.2.2 Flambage/plissement d'un film mince sur substrat circulaire

Pour mieux comprendre l'instabilité de plissement, je me suis intéressé avec mes collègues du M.I.T. (P. Reis, J. Dunkel et al., 2012-2015), au cas de la déformation d'un film mince attaché à un substrat élastique circulaire (2D) contenant une cavité interne dépressurisée<fnref n="2" />. L'effet du substrat sur le film a été modélisé par une force de rappel, déterminée à partir d'une fonction d'Airy (résolution du problème statique d'élasticité linéaire 2D pour le substrat). A partir des équations d'Euler-Lagrange (obtenues via la minimisation d'une énergie potentielle) et d'un développement asymptotique du déplacement, j'ai déterminé les états d'équilibre du film (ordre 0, déformation purement radiale), la pression critique et les modes de l'instabilité (analyse linéaire de stabilité, ordre 1). Les modes de l'instabilité correspondent à un plissement local ou à un flambage global du film. J'ai établi un diagramme de phase montrant que la sélection d'un mode d'instabilité dépend du rayon effectif (rapport entre le rayon du substrat et l'épaisseur du film), du rapport des modules d'Young et de la taille de la cavité interne. Ce diagramme de phase a été corroboré par des simulations ABAQUS. Les résultats pour un substrat circulaire (2D) ont ensuite été étendus au cas d'un substrat sphérique (3D). Une théorie de coque et une minimisation d'une énergie de déformation ont permis d'obtenir une équation de Swift-Hohenberg généralisée pour les états d'équilibre du système. Nous avons montré, via une résolution numérique et une étude paramétrique de cette équation, l'existence d'une variété d'états d'équilibre caractérisés par différents patterns (hexagonaux, labyrinthiques). Nous avons montré que la sélection d'un pattern dépend du rayon effectif, du rapport des modules d'Young et de la contrainte excédentaire (contrainte de compression au-delà du seuil critique) appliquée au film.

<footnote n="2">² Vidéo M.I.T. de l'instabilité de plissement, cas d'un substrat semi-sphérique possédant une cavité intérieure dépressurisée : ici.</footnote>
Finalement, nous avons réduit l'équation de Swift-Hohenberg généralisée à une forme standard à partir de laquelle nous avons établi un diagramme de bifurcation et de phase des motifs de plissement. Les travaux de recherche sur l'instabilité de plissement ont été publiés dans les journaux suivants :

* F. Lopez Jimenez, N. Stoop, R. Lagrange, J. Dunkel, P.M. Reis. Curvature controlled defect localization in elastic surface crystals. Phys. Rev. Lett., 116, 104301, 2016.

* R. Lagrange, F. Lopez Jimenez, D. Terwagne, M. Brojan, P. Reis, From wrinkling to global buckling of a ring on a curved substrate. J. Mech. Phys. Solids, 89, 77-95, 2016.

* N. Stoop, R. Lagrange, D. Terwagne, P.M. Reis, J. Dunkel. Curvature induced symmetry breaking selects elastic surface patterns. Nature Materials, 14, 337-342, 2015.

* M. Brojan, D. Terwagne, R. Lagrange, P.M. Reis. Wrinkling crystallography on spherical surfaces. Proc. Natl. Acad. Sci. U.S.A., 112(1), 14-19, 2015.

## 2.3 Instabilité fluide-élastique d'un faisceau de tubes

L'instabilité fluide-élastique se produit lorsque des structures élastiques sont soumises à l'écoulement d'un fluide. Dans certains cas, la dynamique du fluide peut entraîner des vibrations de la structure, qui à leur tour modifient l'écoulement, créant ainsi un couplage fluide-structure. En tant que structures minces et flexibles, les tubes en U des générateurs de vapeur sont sensibles aux mécanismes d'interaction fluide-structure, en particulier dans leur partie supérieure, soumise à un écoulement transverse diphasique eau-vapeur. Dans ce contexte, je dirige depuis mon arrivée au C.E.A. en 2015 une activité de recherche sur la modélisation des forces fluide-élastiques (et des forces turbulentes) dans des faisceaux de tubes soumis à un écoulement transverse. Dans le cas d'un tube flexible, situé au centre d'un faisceau rigeide à pas carré, j'ai montré que ce dernier devient instable dans la direction de portance, lorsqu'il est soumis à un écoulement monophasique caractérisé par un nombre de Reynolds élevé<fnref n="3" />. À partir d'une analyse dimensionnelle, alimentée par des résultats expérimentaux obtenus sur des maquettes à échelle réduite, j'ai fourni une expression analytique du nombre de Reynolds critique. J'ai montré que le seuil d'instabilité dépend du nombre de Scruton (c'est-à-dire de l'effet masse-amortissement) et du nombre de Stokes (c'est-à-dire de l'effet fréquence/rigidité) du tube : i) pour un nombre de Scruton donné, le tube est d'autant plus stable que le nombre de Stokes est grand, ii) pour un nombre de Stokes donné, un tube est d'autant plus stable que le nombre de Scruton est grand. À la suite de ces résultats obtenus en écoulement monophasique, j'ai dérivé un nouveau critère pour le seuil d'instabilité en écoulement diphasique. Ce critère établit un lien direct avec le

<footnote n="3">³ Vidéo C.E.A. de l'instabilité fluide-élastique en portance, faisceau à pas carré, écoulement monophasique : ici.</footnote>
seuil d'instabilité en écoulement monophasique, ne dépend pas des régimes d'écoulement diphasique, et est corroboré par des expériences en écoulement air-eau. Ce nouveau critère doit être testé pour différentes géométries de faisceaux de tubes (par exemple en pas triangulaire) et pour des fluides diphasiques différents (eau-freon, vapeur-eau), avant sa mise en œuvre dans le code numérique thermo-hydraulique développé par Framatome et le C.E.A. Ce nouveau critère de stabilité a été présenté lors de la conférence internationale "Fluid-Induced Vibrations" que j'ai coorganisée en 2022 à Paris et a été retenu pour faire l'objet d'une publication dans un numéro spécial du journal "J. Fluids and Struct". Les travaux de recherche sur l'instabilité fluide-élastique ont été publiés/soumis dans les journaux suivants :

* R. Lagrange, D. Panunzio, P. Piteau, X. Delaune, J. Antunes. A new criterion for the instability threshold of a square tube bundle subject to two-phase cross-flow. Submitted to J. Fluids Struct., special issue, 2023.

* P. Piteau, D. Panunzio, R. Lagrange, X. Delaune, J. Antunes. Experimental investigation of in-flow fluidelastic instability for rotated triangular tube bundles subjected to single-phase and two-phase transverse flows. Submitted to J. Fluids Struct., special issue, 2023.

* R. Lagrange, P. Piteau, X. Delaune, J. Antunes. Fluid-elastic coefficients in single phase cross flow : dimensional analysis, direct and indirect experimental methods. Proceedings of the ASME, PVP2019, 93984, 2019.

Je développe dans la suite de ce mémoire les thèmes de recherche associés à ces trois types d'instabilité. Chaque chapitre commence par une introduction générale sur le phénomène d'instabilité considéré, suivie d'une publication (post-docotorale, identifiée par une couleur bleue ci-dessus) représentative du travail réalisé et de la méthodologie mise en place pour comprendre, modéliser et prédire les caractéristiques de cette instabilité. Les publications identifiées par une couleur rouge sont reportées en annexes A et B.
### 3. INSTABILITÉ D'UN FLUIDE DANS UN CYLINDRE EN PRÉCESSION

La compréhension de l'écoulement engendré par un mouvement de précession revêt une importance critique dans différents domaines. En aéronautique, certains engins volants contiennent un carburant liquide qui peut être mis en rotation par frottement sur les parois du réservoir. En l'absence de perturbation externe, ce liquide est en rotation solide. Toutefois, en raison de facteurs tels qu'une asymétrie dans la géométrie de l'engin volant, un petit angle de précession peut apparaître et entraîner un changement de direction de l'axe de rotation propre $\Omega_1$. L'axe peut se mettre à tourner autour d'un second axe, noté $\Omega_2$, donnant ainsi à l'engin volant un mouvement de précession par rapport à un référentiel galiléen. Dans le cas où l'écoulement du fluide carburant est en résonance, l'angle initial de précession peut augmenter, ce qui peut avoir une incidence considérable sur la stabilité de l'engin volant. Les problèmes de mécanique des fluides tournants ont également un impact significatif sur de nombreux phénomènes atmosphériques tels que les cyclones et les tornades. La dynamique de ces phénomènes est étroitement liée à la force de Coriolis, qui est omniprésente dans les fluides en rotation et qui agit perpendiculairement à la vitesse d'un corps en déplacement dans un référentiel en rotation. Un cyclone (i.e. une dépression) correspond à une structure fluide en forme de spirale, localisée à la surface de la Terre, en rotation autour de la verticale locale à la vitesse angulaire $\Omega_1 = \Omega_2 \sin \alpha$; $\Omega_2$ étant la vitesse de rotation de la Terre et $\alpha$ la latitude du cyclone. L'axe $\Omega_1$ du cyclone est lui-même en rotation autour de l'axe polaire et précesse donc à la vitesse angulaire $\Omega_2$. Ce mouvement de précession s'effectue alors avec un angle $\theta = 90-\alpha$ degrés. Dans ce contexte, une question importante est de savoir comment la précession via la force de Coriolis influence la dynamique de ce cyclone. Un autre grand domaine d'application des fluides tournants concerne la géophysique interne liée à la précession terrestre (rotation de la Terre autour de l'axe polaire, lui même en rotation autour d'un axe perpendiculaire au plan de l'écliptique). Le mouvement de précession de la Terre joue un rôle clé dans la dynamique de son noyau liquide (composé de métaux en fusion, principalement de fer et de nickel). Ce mouvement crée un forçage sur le noyau liquide qui génère des écoulements de fluides. Ces écoulements peuvent produire des courants électriques, donnant naissance à des champs magnétiques. Ces derniers renforcent alors à leur tour les courants, créant ainsi un effet dynamo. La maîtrise de la mécanique des fluides tournants est donc essentielle pour comprendre le processus de génération du champ magnétique terrestre.

La dynamique d'un fluide Newtonien (masse volumique $\rho$, viscosité cinématique
tique $\nu$) dans un cylindre (hauteur $H$, rayon $R$) en précession est gouvernée par l'équation de Navier-Stokes

$$ \frac{\partial \mathbf{U}}{\partial T} + (\mathbf{U} \cdot \nabla) \mathbf{U} = -\nabla \left(\frac{P}{\rho}\right) - 2\Omega \times \mathbf{U} - \Omega \times (\Omega \times \mathbf{R}) - \left(\frac{d\Omega}{dT}\right)_a \times \mathbf{R} + \nu \Delta \mathbf{U}, \quad (3.1) $$

avec $\mathbf{R}$ le vecteur position, $\mathbf{U}$ la vitesse du fluide dans le repère en rotation (représenté par le vecteur rotation $\Omega = \Omega_1 + \Omega_2$ et par sa variation temporelle $(d\Omega/dT)_a$ par rapport à un repère absolu). Dans cette équation le terme $\Omega \times (\Omega \times \mathbf{R})$ correspond à la force centrifuge. La composante $2\Omega \times \mathbf{U}$ correspond à la force de Coriolis. Le terme $(d\Omega/dT)_a \times \mathbf{R}$ correspond à l'accélération du repère tournant par rapport au repère absolu et est appélé accélération de Poincaré. Ce terme est essentiel dans un mouvement de précession car il représente un terme de forçage spécifique sur le fluide. Une analyse dimensionnelle de l'équation de Navier-Stokes montre que la dynamique du fluide est gouvernée par quatre nombres sans dimension

* le rapport d'aspect du cylindre, $h = H/R$,

* le rapport de fréquence, $\omega = \Omega_1/\Omega$,

* le nombre de Reynolds, $Re = \Omega R^2/\nu$, représentant le rapport entre les ordres de grandeur du terme de Coriolis et du terme visqueux de l'équation de Navier-Stokes,

* le nombre de Rossby, $Ro = U/(\Omega R)$, représentant le rapport entre les ordres de grandeurs du terme de transport convectif de l'équation de Navier-Stokes et du terme de force de Coriolis.

Pour des nombres de Reynolds et de Rossby modérés, nous avons montré que l'écoulement d'un fluide dans un cylindre en précession est une superposition d'un cisaillement et d'une infinité d'ondes hydrodynamiques, appelées modes de Kelvin. L'amplitude d'un mode de Kelvin diverge lorsque sa fréquence propre (fonction du rapport d'aspect $h$) est égale à la fréquence de précession $\Omega$. On dit alors que le cylindre de rapport d'aspect $h$ est résonnant. Dans un tel cylindre, le mode divergeant devient prépondérant devant le cisaillement et les autres modes de Kelvin. Il constitue alors l'écoulement de base à l'intérieur du cylindre en précession, voir Fig. 3.1. Lorsque le nombre de Reynolds augmente, cet écoulement de base devient instable. Grâce à des mesures PIV (Particle Image Velocimetry), nous avons montré que cette instabilité est due à une interaction entre le mode résonnant et deux modes de Kelvin libres (solutions de l'équation de Navier-Stokes linéarisée et homogène), voir Fig. 3.2. Une analyse linéaire de stabilité basée sur un mécanisme de résonance triadique entre modes de Kelvin (i.e. interaction entre trois modes de Kelvin via le terme non-linéaire de l'équation de Navier-Stokes) permet de prédire correctement le taux de croissance et le seuil de l'instabilité, voir Fig. 3.3. Une analyse faiblement non-linéaire permet de plus de prédire la saturation des modes instables qui forcent un mode axisymétrique (mode géostrophique) au travers des couches limites. Ce mode limite alors l'amplitude de l'instabilité en modifiant la fréquence des modes instables.
<figure><img src="image_2.png" /><figcaption>FIGURE 3.1 : Visualisations avec des particules réfléchissantes des modes de Kelvin dans un cylindre résonnant, obtenues pour différents rapports d'aspect *h*, par McEwan (1970).</figcaption></figure>

En plus de prédire la saturation de l'instabilité, la théorie faiblement non-linéaire prédit également des régimes d'intermittence et une transition vers le chaos. La prédiction de l'instabilité et de sa saturation sont en excellent accord avec les résultats expérimentaux, même pour des nombres de Reynolds une décade au dessus du seuil. Les résultats de cette analyse, obtenus dans le cadre de ma thèse, sont cependant uniquement valables dans le cas d'un cylindre résonnant.

L'extension de l'analyse de stabilité au cas d'un cylindre quelconque a été réalisée lors de mon séjour postdocoral au M.I.T., et est l'objet de la publication qui suit, parue dans la revue *Comptes Rendus Mécanique*, sous le titre "Triadic instability of a non-resonant precessing fluid cylinder".
<figure><img src="image_3.png" /><figcaption>FIGURE 3.2 : Champs de vorticité axial instantané (PIV) des modes de Kelvin dans un cylindre en précession. Mode de base résonnant forcé par la précession (a) et modes de l'instabilité (b et c). Les modes (a), (b) et (c) interagissent via un mécanisme de résonance triadique. Vidéo I.R.P.H.E. de l'instabilité de précession : ici.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIGURE 3.3 : Diagramme de stabilité d'un fluide dans un cylindre en précession, cas d'un cylindre résonnant. La droite verte correspond à la prédiction théorique du seuil d'instabilité.</figcaption></figure>
### Abstract

Flows forced by a precessional motion can exhibit instabilities of crucial importance, whether they concern the fuel of a flying object or the liquid core of a telluric planet. So far, stability analyses of these flows have focused on the special case of a resonant forcing. Here, we address the instability of the flow inside a precessing cylinder in the general case. We first show that the base flow forced by the cylinder precession is a superposition of a vertical or horizontal shear flow and an infinite sum of forced modes. We then perform a linear stability analysis of this base flow based on the mechanism of triadic resonance between a forced mode and two free Kelvin modes. Finally, we derive the amplitude equations of the free Kelvin modes and obtain an expression of the instability threshold and growth rate.

## 3.1 Introduction

Knowing the flow forced by precessional motion is of critical importance in several domains. In aeronautics, liquid propellants contained in flying objects can become resonant for specific geometries of their tank. The resulting flow can then create a destabilizing torque on the objects and dangerously modify their trajectories Stewartson (1958); Karpov (1965); Gans (1984); Garg et al. (1986); Herbert (1986); Pocha (1987); Agrawal (1993); Bao & Pascal (1997); Lambelin et al. (2009). Modeling the flow inside a precessing cylinder is thus a necessary step to design tank geometries that avoid these unwanted resonances.

In geophysics, most planets have a motion of slow precession, which is mainly governed by the planet aspect ratio. In the presence of a liquid core, this precessional motion creates a weak forcing that can drastically modify the flow inside the core due to the presence of resonances and critical layers. Flows inside liquid planet cores are of primordial interest to understand the generation of magnetic field by dynamo effect. For the present-day Earth, the magnetic field is likely due to the convection between the hot solid inner core and the colder mantle Malkus (1989); Grote et al. (1999, 2000); Grote & Busse (2001). However, a magnetic field was present on the early Earth although a solid inner core was not yet present. At that time, other mechanisms have generated and sustained the Earth's magnetic field. Tides (leading to elliptic streamlines) have often thought to be a source of energy sufficient for geodynamo Malkus (1989); Cébron et al. (2011), but recently it has also been shown numerically that precession could generate a magnetic field Tilgner (2005); Nore et al. (2011) (although this was not clearly proven for the case of the Earth). Besides, there is still some debate as to whether the production of kinetic energy due to precession is sufficient to balance the Ohmic energy loss induced by the magnetic field Rochester et al. (1975); Loper (1975); Roberts & Gubbins (1987); Vanyo et al. (1995); Kerswell (1996); Gans (1970); Tilgner (2000); Gailitis et al. (2002). In any case, even if precession is not the cause of magnetic field production on Earth, it might be different on other telluric planets.

To study the flow driven by a precessional motion, the cylindrical geom-
etry offers a good alternative to a planet-like spheroidal geometry because of its simplicity. In a precessing cylinder, the base flow is a sum of a shear flow and an infinite set of forced modes Kelvin (1880); Greenspan (1968). For particular precessional frequencies, a forced mode may become resonant when the height of the cylinder equals an odd number of half wavelengths Fultz (1959); McEwan (1970); Kerswell & Barenghi (1995); Meunier et al. (2008). In the framework of an inviscid and linear theory, this resonance leads to a divergence of the forced mode amplitude. Viscous effects however may saturate this amplitude to a value scaling as the inverse square root of the Ekman number (due to Ekman layers) (Gans, 1970). Nonlinear effects can also saturate the amplitude at a value scaling as the cubic root of the forcing (Meunier et al., 2008). This nonlinear saturation is due to the presence of a strong axisymmetric zonal flow (also called geostrophic flow), which tends to decrease the solid body rotation and thus detune the resonance of the forced mode Kong et al. (2015).

When the Ekman number is decreased or the precessing angle is increased above a critical value, a resonant forced mode can become unstable McEwan (1970); Thompson (1970); Manasseh (1992, 1994); Kobine (1995); Manasseh (1996). For small tilt angles, we have shown that this instability is due to a triadic resonance between the resonant forced mode and two free Kelvin modes Lagrange et al. (2008, 2009, 2011). However, outside of resonances, the forced modes have a small amplitude and the base flow is not made of a single mode anymore. Here, our principal objective is to perform a stability analysis of the complete base flow (made of a shear flow and a sum of forced modes) in the case of a non-resonant precessing fluid cylinder. We will consider the triadic interaction of the base flow with two free Kelvin modes to determine the conditions of instability and derive an expression for the growth rate.

This paper is organized as follows. Section 3.2 presents the problem of a precessing cylinder by introducing the governing equations. In this section we determine the off-resonance base flow and discuss the symmetry properties of the problem. In § 3.3, we develop a linear stability analysis of the complete base flow, based on the mechanism of triadic resonance. We discuss the conditions of resonance, derive the amplitude equations of the instability modes and provide an analytical expression of the instability growth rate and threshold. In § 3.4 we present an application for a particular case, and determine its stability diagram. Finally, some conclusions are drawn and discussed in the context of the transition to turbulence in precessing flows.

## 3.2 Formulation of the problem

Consider a cylinder of radius $R$, height $H$, axis of revolution along $\hat{k}$, entirely filled with a Newtonian fluid of kinematic viscosity $\nu$. The cylinder rotates at the angular speed $\Omega_1$ about $\hat{k}$, which also rotates at the angular speed $\Omega_2$ about the vertical axis and we denote by $\theta$ the precession angle, i.e. the angle between these two axes of rotation (Fig. 3.4).

To make the problem dimensionless, we introduce four numbers: the aspect ratio $h = H/R$, the frequency ratio $\omega = \Omega_1/\Omega$, with $\Omega = \Omega_1 + \Omega_2 \cos \theta$,
the Ekman number $Ek = \nu/(\Omega R^2)$ and the Rossby number $Ro = \Omega_2 \sin\theta/\Omega$, which will be assumed asymptotically small, i.e. $Ro \ll 1$ (weak precession). The dimensionless flow velocity in the cylinder's frame of reference $(O, \hat{\mathbf{i}}, \hat{\mathbf{j}}, \hat{\mathbf{k}})$ is denoted by $\mathbf{u} = \mathbf{U}/(R\Omega)$. The dimensionless cylindrical coordinates are $(r, \varphi, z)$, where $z = 0$ corresponds to the mid-height section of the cylinder and we note $\mathbf{r}$ the position vector of a fluid particle. In the cylinder's frame of reference, the dimensionless Euler equations (assuming an inviscid fluid) are Meunier et al. (2008); Lagrange et al. (2011)

$$ \frac{\partial \mathbf{u}}{\partial t} + 2(\hat{\mathbf{k}} + R\Omega \hat{\mathbf{\delta}}) \times \mathbf{u} + \nabla p = -2R\Omega \omega r \cos(\omega t + \varphi) \hat{\mathbf{k}} + \mathbf{u} \times (\nabla \times \mathbf{u}), \quad (1a) $$

$$ \nabla \cdot \mathbf{u} = 0, \qquad (1b) $$

with $\hat{\mathbf{\delta}} = \cos(\omega t) \hat{\mathbf{i}} - \sin(\omega t) \hat{\mathbf{j}}$. On the left hand side (LHS) of (1a), the first term is inertia, the second term is the Coriolis force and $p$ is the dimensionless pressure field defined as

$$ \begin{aligned} p = & \frac{P}{\rho\Omega^2 R^2} - \frac{1}{2}r^2 + Ro|1 - \omega|rz \cos(\omega t + \varphi) + \gamma_O \cdot \mathbf{r} \\ & - \frac{1}{2}Ro^2[z^2 + r^2 \sin^2(\omega t + \varphi)] \\ & + \frac{1}{2}\mathbf{u}^2, \end{aligned} \qquad (2) $$

where $\gamma_O = \Gamma_O/R\Omega^2$ is the dimensionless acceleration of the cylinder centroid $O$. On the right hand side (RHS) of (1a), the first term is the forcing due to precession, the second term is the convective nonlinear term. At this point, it is convenient to introduce the four components vector $\mathbf{v} = (\mathbf{u}, p)^T$ and recast equations (1) into a matrix formulation

$$ \left( \frac{\partial}{\partial t} \mathcal{I} + \mathcal{M} \right) \mathbf{v} = 2R\Omega \mathbf{F}_0 \cos(\omega t + \varphi) + \mathbf{N}(\mathbf{v}, \mathbf{v}) + Ro (\mathcal{D}^{i}(\omega t + \varphi) + \text{c.c.})\mathbf{v}, \quad (3) $$

where operators $\mathcal{I}$, $\mathcal{M}$, $\mathcal{D}$, the forcing vector $\mathbf{F}_0$ and the bilinear function $\mathbf{N}$ are reported in Appendix 3.6.1. The symbol c.c. stands for the complex conjugate.

### 3.2.1 Base flow

In the limit of small Rossby numbers, the base flow forced by the precessional motion can be found by solving the Euler equations (3) at first order. We thus have to solve the following inhomogeneous linear differential equation for $\mathbf{v}$

$$ \left( \frac{\partial}{\partial t} \mathcal{I} + \mathcal{M} \right) \mathbf{v} = 2R\Omega \mathbf{F}_0 \cos(\omega t + \varphi). \qquad (4) $$

Projecting this equation onto $\hat{\mathbf{k}}$ yields

$$ \frac{\partial v_z}{\partial t} + \frac{\partial v_p}{\partial z} = -2R\Omega r\omega \cos(\omega t + \varphi). \qquad (5) $$
There are two particular solutions of this equation, which can be found by assuming either $p = 0$ or $u_z = 0$. The first assumption leads to a vertical shear given by

$$ \mathbf{v} = Ro \mathbf{v}_{\text{shear}}^{\text{V}} e^{i(\omega t + \varphi)} + \text{c.c.} \quad \text{with} \quad \mathbf{v}_{\text{shear}}^{\text{V}} = \begin{pmatrix} 0 \\ 0 \\ ir \\ 0 \end{pmatrix}. \quad (6) $$

This vertical shear is schematically shown in Fig. 3.4 (middle) and corresponds to the flow that would be found in a cylinder with infinite height. However, there is an alternative particular solution given by the second assumption ($u_z = 0$) which leads to a horizontal shear flow

$$ \mathbf{v} = Ro \mathbf{v}_{\text{shear}}^{\text{H}} e^{i(\omega t + \varphi)} + \text{c.c.} \quad \text{with} \quad \mathbf{v}_{\text{shear}}^{\text{H}} = \frac{\omega z}{2 - \omega} \begin{pmatrix} i \\ -1 \\ 0 \\ r(\omega - 2) \end{pmatrix}. \quad (7) $$

This horizontal shear is schematically shown in Fig. 3.4 (right). None of these particular solutions satisfy the boundary conditions on the cylinder walls, and they have to be completed with homogeneous solutions of (4). In the case of the vertical shear $\mathbf{v}_{\text{shear}}^{\text{V}}$, the full solution with proper boundary conditions has been found to be (e.g. Meunier et al. (2008))

$$ \mathbf{v} = Ro \mathbf{v}_{\text{base}}^{\text{V}} = Ro \left( \mathbf{v}_{\text{shear}}^{\text{V}} + \sum_{j=1}^{\infty} a_j^{\text{V}} \mathbf{v}_j^{\text{V}} \right) e^{i(\omega t + \varphi)} + \text{c.c.}, \quad (8) $$

where $\mathbf{v}_j^{\text{V}}$ are forced modes of azimuthal wavenumber $m = 1$, frequency $\omega$, and axial wavenumbers $k_j^{\text{V}}$. The amplitudes $a_j^{\text{V}}$ and the structure of the forced modes are given in Appendix 3.6.1. Since here, these forced modes compensate the normal flow on the top and bottom walls of the cylinder, they have zero radial velocity $u_r = 0$ at $r = 1$. This boundary condition forces the axial wavenumbers $k_j^{\text{V}}$ through the dispersion relation $D(1, \omega, k_j^{\text{V}}) = 0$ given in Appendix 3.6.1 (it is stressed that the wavenumbers $k_j^{\text{V}}$ are not a multiple of $\pi/h$ since the forced modes do not satisfy the boundary condition $u_z = 0$ on the top and bottom wall).

In the case of a horizontal shear $\mathbf{v}_{\text{shear}}^{\text{H}}$, a similar procedure yields the following base flow solution

$$ \mathbf{v} = Ro \mathbf{v}_{\text{base}}^{\text{H}} = Ro \left( \mathbf{v}_{\text{shear}}^{\text{H}} + \sum_{j=1}^{\infty} a_j^{\text{H}} \mathbf{v}_j^{\text{H}} \right) e^{i(\omega t + \varphi)} + \text{c.c.}, \quad (9) $$

The amplitudes $a_j^{\text{H}}$ and the structure of the forced modes $\mathbf{v}_j^{\text{H}}$ are given in Appendix 3.6.1. In this case, note that the axial wavenumbers $k_j^{\text{H}}$ are odd multiples of $\pi/h$ ($k_j^{\text{H}} = (2j - 1)\pi/h$). Here, the forced modes $\mathbf{v}_j^{\text{H}}$ satisfy the boundary conditions on the top and bottom walls, but not on the lateral walls ($u_r \neq 0$ at $r = 1$), because these forced modes have been added to compensate an horizontal shear flow.
In both cases ($v_{base}^V$ or $v_{base}^H$), the base flow can be written as a superposition of a shear flow (vertical or horizontal respectively) and a sum of forced modes. These forced modes are similar to Kelvin modes Kelvin (1880); Greenspan (1968). However, we do not consider them as “real” Kelvin modes because they do not satisfy the free-slip boundary condition on all the cylinder walls. Figure 3.5 shows the two base flow solutions with an increasing number of forced modes. It is clear that both solutions tend to be equal $v_{base}^V = v_{base}^H = v_{base}$ when the number of forced modes is increased.

### 3.2.2 Symmetry properties

The stability analysis of the base flow described above will be presented below in §3. Before that, it is worth recalling some symmetry properties of the problem. It can first be noted that the vertical shear flow $v_{shear}^V$ has only a vertical velocity component $v_z$, which is an even function of $z$. In contrast, the horizontal shear flow $v_{shear}^H$ has no vertical component but its radial, azimuthal, and pressure components are odd functions of $z$. Both quadrivectors are thus of the type:

$$ \mathbf{v}^{-} = \begin{pmatrix} f^{-}(z, r, \varphi, t) \\ f^{-}(z, r, \varphi, t) \\ f^{+}(z, r, \varphi, t) \\ f^{-}(z, r, \varphi, t) \end{pmatrix}, \qquad (10) $$

where $f^+$ (resp. $f^-$) denotes even (resp. odd) functions of $z$. Equation (42) shows that the forced modes $v_j^V$ and $v_j^H$ also have the same $z$-parity. As a consequence, the base flow calculated above is of the type $\mathbf{v}^-$. This could have been anticipated since the first order terms of the Navier-Stokes equations only force this symmetry. However, at higher orders, the operators can generate a flow with the opposite symmetry of the type

$$ \mathbf{v}^{+} = \begin{pmatrix} f^{+}(z, r, \varphi, t) \\ f^{+}(z, r, \varphi, t) \\ f^{-}(z, r, \varphi, t) \\ f^{+}(z, r, \varphi, t) \end{pmatrix}. \qquad (11) $$

It is easy to show that the operators $\mathcal{D}$, $\overline{\mathcal{D}}$, and $\mathbf{N}$ have the following properties

$$ \mathcal{D}\mathbf{v}^{-} \sim \overline{\mathcal{D}}\mathbf{v}^{-} \sim \mathbf{v}^{+}, \qquad (12a) $$

$$ \mathbf{N}(\mathbf{v}_{\text{base}}, \mathbf{v}^{-}) \sim \mathbf{N}(\overline{\mathbf{v}}_{\text{base}}, \mathbf{v}^{-}) \sim \mathbf{v}^{+}, \qquad (12b) $$

$$ \mathbf{N}(\mathbf{v}^{-}, \mathbf{v}_{\text{base}}) \sim \mathbf{N}(\mathbf{v}^{-}, \overline{\mathbf{v}}_{\text{base}}) \sim \mathbf{v}^{+}, \qquad (12c) $$

where the symbol $\sim$ means “same parity as” and $v_{base}$ corresponds to either the vertical base flow $v_{base}^V$ or the horizontal base flow $v_{base}^H$. This means that, even if the base flow is of the type $\mathbf{v}^-$ at first order, the nonlinear terms may introduce a different symmetry in the problem. It is actually easy to show that the previous equations remain valid under the permutation of signs $(+, -) \rightarrow (-, +)$.
<figure><img src="image_5.png" /><figcaption>Figure 3.4: Sketch of a precessing cylinder of radius $R$ and height $H$ (left). The cylinder is rotating around its own axis with angular velocity $\Omega_1$, and this axis itself rotates around a second axis with angular velocity $\Omega_2$. The precession forces a shear flow which can be taken vertical (middle) or horizontal (right), but which does not respect the boundary conditions.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 3.5: Comparisons of the base flows computed from eqns. (8) and (9) for $h = 2.3$, $\omega = 1.34$, $t = 0$. (a-c) Flow $v_{\text{base}}^{\text{H}}$ obtained from eqn. (9) when the sum of forced modes is truncated to 0, 1, and 5 modes respectively. (d-f) Same plots for $v_{\text{base}}^{\text{V}}$ from eqn. (8). In both cases, the vectors show the projection of the velocity field in the plane $z = h/3$ and the color-coded map shows $v_z$ in the same plane.</figcaption></figure>
A direct consequence of these symmetry properties is that a triadic resonance cannot develop if the perturbation has a single parity ($\mathbf{v}^+$ or $\mathbf{v}^-$). Indeed, in a mechanism of triadic resonance, the base flow interacts with two free Kelvin modes $\mathbf{v}_1$ and $\mathbf{v}_2$ through the nonlinear operator $\mathbf{N}$. Let us assume that the free Kelvin modes have the same symmetry $\mathbf{v}^-$. The growth of the first free Kelvin mode is due to the nonlinear interaction of the second free Kelvin mode with the base flow, via the terms $\mathbf{N}(\mathbf{v}_{\text{base}}, \mathbf{v}_2^-)$, $\mathbf{N}(\mathbf{v}_2^-, \mathbf{v}_{\text{base}})$ and $\mathcal{D}\mathbf{v}_2^-$ which have the opposite symmetry $\mathbf{v}^+$. These forcing terms are thus perpendicular to the first free Kelvin mode and such a triadic resonance is non constructive. This can be properly shown by defining the dot product

$$ \langle \mathbf{X}, \mathbf{Y} \rangle = \int_V (\overline{X_r} Y_r + \overline{X_\varphi} Y_\varphi + \overline{X_z} Y_z + \overline{X_p} Y_p) d^3V, \quad (13) $$

where $\overline{X}$ refers to the conjugate of $X$ and $V$ is the volume of the cylinder. It is then straightforward to show, using (12a-c) that the dot products $\langle \mathbf{v}_1^-, \mathbf{N}(\mathbf{v}_2^-, \mathbf{v}_{\text{base}}) \rangle$, $\langle \mathbf{v}_1^-, \mathbf{N}(\mathbf{v}_{\text{base}}, \mathbf{v}_2^-) \rangle$ and $\langle \mathbf{v}_1^-, \mathcal{D}\mathbf{v}_2^- \rangle$ vanish because they only contain terms of the form

$$ \langle \mathbf{v}^+, \mathbf{v}^- \rangle = \int_{-h/2}^{h/2} f^+(z)f^-(z)dz = 0. \quad (14) $$

The same reasoning applies to free Kelvin modes with a symmetry $\mathbf{v}^+$. The general conclusion is that the constructive triadic resonances must couple an even free Kelvin mode $\mathbf{v}^+$ with an odd free Kelvin mode $\mathbf{v}^-$. We will now use this property to restrict the number of possible instabilities that may arise in the linear stability analysis.

## 3.3 Linear stability analysis

To study the stability of the base flow, we introduce a small perturbation in form of a four-components vector $\tilde{\mathbf{v}} = (\tilde{\mathbf{u}}, \tilde{\mathbf{p}})^T$, so that the total flow is

$$ \mathbf{v} = Ro \mathbf{v}_{\text{base}} + \tilde{\mathbf{v}} + o(Ro), \quad (15) $$

where $\mathbf{v}_{\text{base}}$ is either $\mathbf{v}_{\text{base}}^V$ or $\mathbf{v}_{\text{base}}^H$. Inserting this expansion into (3) yields an equation for the perturbation vector

$$ \left( \frac{\partial}{\partial t} \mathcal{I} + \mathcal{M} \right) \tilde{\mathbf{v}} = Ro \left[ \mathbf{N}(\mathbf{v}_{\text{base}}, \tilde{\mathbf{v}}) + \mathbf{N}(\tilde{\mathbf{v}}, \mathbf{v}_{\text{base}}) + (\mathcal{D}^{i(\omega t + \varphi)} + \text{c.c.}) \tilde{\mathbf{v}} \right] \\ + o(Ro) + o(|\tilde{\mathbf{v}}|), \quad (16) $$

where $|\tilde{\mathbf{v}}| = \sqrt{\langle \tilde{\mathbf{v}}, \tilde{\mathbf{v}} \rangle}$ is the magnitude of $\tilde{\mathbf{v}}$. The first two terms on the RHS of (16) represent the nonlinear interactions between the base flow and the perturbation. The third term represents the interaction between the forcing due to precession and the perturbation. The perturbation vector satisfies the inviscid boundary condition

$$ \tilde{\mathbf{u}} \cdot \mathbf{n} = 0 \text{ at the walls } (r = 1 \text{ or } z = \pm h/2). \quad (17) $$
To solve (16) and (17), we use a multiscale expansion in time with t a rapid time scale and $\tau = Rot$ a slow time scale. We then expand $\tilde{\mathbf{v}}$ as

$$ \tilde{\mathbf{v}} = \tilde{\mathbf{v}}_0(\mathbf{r}, \tau, t) + Ro \tilde{\mathbf{v}}_1(\mathbf{r}, \tau, t) + o(Ro). \quad (18) $$

Inserting (18) into (16) yields two equations: one of order 1 and one of order $Ro$ that shall be studied now. The equation at order one gives the form of the free Kelvin modes, and the equation at order $Ro$ gives their slow time dynamics, hence their stability properties.

### 3.3.1 Order 1: free Kelvin modes

At first order, the equation (16) and the inviscid boundary condition $\tilde{\mathbf{u}} \cdot \mathbf{n}$ are

$$ \left(\frac{\partial \mathcal{I}}{\partial t} + \mathcal{M}\right) \tilde{\mathbf{v}}_0 = \mathbf{0}, \quad (19a) $$

$$ \tilde{\mathbf{u}}_0 \cdot \mathbf{n} = 0 \quad \text{at the walls} \ (r = 1 \text{ or } z = \pm h/2). \quad (19b) $$

The solution to this homogenous problem is a linear combination of free Kelvin modes with different $z$-parities Lagrange et al. (2011)

$$ \tilde{\mathbf{v}}_0 = \sum_{l=1}^{\infty} A_l^+ \mathbf{v}_l^+ e^{i(\omega_l t + m_l \varphi)} + \sum_{l=1}^{\infty} A_l^- \mathbf{v}_l^- e^{i(\omega_l t + m_l \varphi)} + \text{c.c.} \quad (20) $$

Vectors $\mathbf{v}_l^+$ (resp. $\mathbf{v}_l^-$) have axial wavenumbers $k_l^+$ (resp. $k_l^-$) which are even (resp. odd) multiple of $\pi/h$ in order to respect the condition of no normal flow at the top and bottom ($z = \pm h/2$). This property is interesting because the wavenumbers are separated into two families, which will restrict the number of possible triadic resonances. The components of the free Kelvin modes are given in Appendix 3.6.1. In (20), $A_l^\pm$, $m_l$, and $\omega_l$ are the amplitude, azimuthal wavenumber, and angular frequency of the free Kelvin mode $\mathbf{v}_l^\pm e^{i(\omega_l t + m_l \varphi)}$. The wavenumbers are connected through the dispersion relation $D(m_l, \omega_l, k_l^\pm)$ such that the radial velocity of the mode vanishes at the cylinder wall $r = 1$.

To examine the mechanism of triadic resonance, the perturbation $\tilde{\mathbf{v}}_0$ is reduced to a combination of two free Kelvin modes $\mathbf{v}_1$ and $\mathbf{v}_2$ with unknown amplitudes $A_1(\tau)$ and $A_2(\tau)$

$$ \tilde{\mathbf{v}}_0 = A_1 \mathbf{v}_1 e^{i(\omega_1 t + m_1 \varphi)} + A_2 \mathbf{v}_2 e^{i(\omega_2 t + m_2 \varphi)} + \text{c.c.} \quad (21) $$

From now on, we attribute index 2 to the mode with the highest azimuthal wavenumber (i.e. $m_2 > m_1$).

### 3.3.2 Triadic resonance

We know from the symmetry properties presented in §3.2.2 that a triadic resonance between the base flow and the two free Kelvin modes can be constructive only if it involves free Kelvin modes with different $z$-parities. Therefore, the wavenumbers $k_1$ and $k_2$ must be integer multiple of $\pi/h$ with
different parities. It follows that the difference between the two wavenumbers must be an odd multiple of $\pi/h$:

$$k_2 - k_1 = (2p - 1)\pi/h, \qquad (22)$$

with $p$ an integer.

In addition, the base flow will interact with two free Kelvin modes if the operator $\mathbf{N}$ appropriately couples their time and azimuthal Fourier components. Thus the coupling term $\mathbf{N}(\mathbf{v}_{\text{base}}, \tilde{\mathbf{v}}_0) + \mathbf{N}(\tilde{\mathbf{v}}_0, \mathbf{v}_{\text{base}})$ in (16) has to contain the same Fourier components as $\mathbf{v}_l e^{i(\omega_l t + m_l \varphi)}$, $l = 1, 2$. Since these terms have the following time and azimuthal Fourier components

$$\mathbf{v}_l e^{i(\omega_l t + m_l \varphi)} : \text{Fourier components } (m_l, \omega_l), \qquad (23a)$$

$$\mathbf{N}(\mathbf{v}_{\text{base}}, \tilde{\mathbf{v}}_0) + \mathbf{N}(\tilde{\mathbf{v}}_0, \mathbf{v}_{\text{base}}) : \text{Fourier components}$$

$$(m_l + 1, \omega_l + \omega), (m_l - 1, \omega_l - \omega), \qquad (23b)$$

the base flow will interact with the two free Kelvin modes if

$$m_2 - m_1 = 1, \qquad (24a)$$

$$\omega_2 - \omega_1 = \omega. \qquad (24b)$$

We recognize on the RHS the azimuthal wavenumber $m_{\text{base}} = 1$ and the angular frequency $\omega_{\text{base}} = \omega$ of the base flow. The conditions (24a-b) are characteristic of triadic resonances occurring in various domains (surface waves, plate vibrations, etc.), and are the key ingredient of weak (or wave) turbulence theory.

To find a pair of free Kelvin modes that fulfill the conditions of resonance (24a-b), we proceed as shown in Fig 3.6. In the plane $(k_2, \omega_2)$ we plot the dispersion relation of the free Kelvin modes $(m_2, \omega_2, k_2)$ and the dispersion relation of the free Kelvin modes $(m_1, \omega_1, k_1)$ translated horizontally by $(2p - 1)\pi/h$ (with $p$ an arbitrary integer) and vertically by $\omega$. The intersection points correspond to free Kelvin modes satisfying the conditions of resonance (24a-b) and the condition induced by parity (22). In addition, the free Kelvin modes have to fit in the cylinder, so their axial wavenumbers need to be multiple of $\pi/h$. It means that the intersection point has to lie on the vertical dotted line of Fig 3.6. Such a tuned triadic resonance only occurs for particular aspect ratios (for a given forcing frequency $\omega$). Otherwise, in most cases, the three curves do not intersect at the same point. Figure 3.6 shows an example of a tuned triadic resonance for $h = 2.3$, $\omega = 1.34$ with $m_1 = 2$ and $m_2 = 3$. The label for these points is $(m_2, l_1, l_2)$, where $l_{1,2}$ is the branch number of the dispersion relations. There are an infinity of possible triadic resonances and they can all be studied in the inviscid case. However, when viscous effects are taken into account, the highest wavenumbers will be damped, such that in practice only the lowest axial, azimuthal and radial wavenumbers may be treated.

The procedure outline above shows how the axial, azimuthal and radial wavenumbers of two free Kelvin modes can be found to permit triadic resonance. In the next section, we consider the Euler equations at the next order in $Ro$, in order to calculate the slow temporal evolution of these free Kelvin modes and determine their stability.
<figure><img src="image_7.png" /><figcaption>Figure 3.6: Dispersion relations of the free Kelvin modes with azimuthal wavenumbers $m_1 = 2$ (dashed lines) and $m_2 = 3$ (solid lines). The dispersion relation for $m_1 = 2$ is translated by $\pi/h$ ($h = 2.3$) along the abscissae and by $\omega = 1.34$ along the ordinate. Vertical dotted lines indicate the discretisation of the axial wavenumber as a multiple of $\pi/h$ imposed by the inviscid boundary condition at the top and bottom of the cylinder. A point lying at the intersection of the branches of the dispersion relations and a vertical line corresponds to a pair of resonant free Kelvin modes. The combination ($m_2 = 3, l_1 = 1, l_2 = 1$) (marked with a black circle) is an example of resonant free Kelvin modes. Coordinates for this point are $k_2 = 2\pi/h$ and $\omega_2 = 0.874$.</figcaption></figure>
### 3.3.3 Order Ro: slow time equations

At order $Ro$, equation (16) becomes

$$ \left( \frac{\partial \mathcal{I}}{\partial t} + \mathcal{M} \right) \tilde{\mathbf{v}}_1 = \mathbf{N}(\mathbf{v}_{\text{base}}, \tilde{\mathbf{v}}_0) + \mathbf{N}(\tilde{\mathbf{v}}_0, \mathbf{v}_{\text{base}}) + \left[ \left( \mathcal{D}e^{i(\omega t + \varphi)} + \text{c.c.} \right) - \frac{\partial \mathcal{I}}{\partial \tau} \right] \tilde{\mathbf{v}}_0. \quad (25) $$

This $O(Ro)$ problem is linear, with a forcing term given by the RHS of (25). To avoid secular terms in the solution $\tilde{\mathbf{v}}_1$, the RHS must be orthogonal to the kernel of the LHS operator. This kernel being spanned by the free Kelvin modes, themselves given by the $O(1)$ problem solved above, a solvability condition is obtained by taking the dot product of (25) with $\mathbf{v}_l e^{i(\omega_l t + m_l \varphi)}$, $l = 1, 2$. Since the problem is self-adjoint, i.e. $\langle \mathbf{v}_l e^{i(\omega_l t + m_l \varphi)}, (\partial \mathcal{I}/\partial t + \mathcal{M}) \tilde{\mathbf{v}}_1 \rangle = 0$, we show in Appendix 3.6.2 that the slow time equations for $A_1$ and $A_2$ are

$$ \frac{dA_1}{d\tau} = c_1 A_2 \quad \text{with} \quad c_1 = \frac{\bar{d}_{12} + n_{1s} + \sum_j \bar{a}_j n_{1j}}{\langle \mathbf{v}_1, \mathcal{I} \mathbf{v}_1 \rangle}, \quad (26a) $$

$$ \frac{dA_2}{d\tau} = c_2 A_1 \quad \text{with} \quad c_2 = \frac{d_{21} + n_{2s} + \sum_j a_j n_{2j}}{\langle \mathbf{v}_2, \mathcal{I} \mathbf{v}_2 \rangle}. \quad (26b) $$

The terms $dA_l/d\tau$ come from the dot products $\langle \mathbf{v}_l e^{i(\omega_l t + m_l \varphi)}, \frac{\partial \mathcal{I}}{\partial \tau} \tilde{\mathbf{v}}_0 \rangle$. The terms $\bar{d}_{12}$ and $d_{21}$ represent the interaction between the resonant free Kelvin modes and the forcing due to precession. They come from the dot products

$$ \langle \mathbf{v}_l e^{i(\omega_l t + m_l \varphi)}, (\mathcal{D}e^{i(\omega t + \varphi)} + \text{c.c.}) \tilde{\mathbf{v}}_0 \rangle. \quad (27) $$

The terms $n_{1s}$ and $n_{2s}$ represent the nonlinear interactions between the resonant free Kelvin modes and the shear part of the base flow. They come from the dot products

$$ n_{ls} = \langle \mathbf{v}_l e^{i(\omega_l t + m_l \varphi)}, \mathbf{N}(\mathbf{v}_s e^{i(\omega t + \varphi)} + \text{c.c.}, \tilde{\mathbf{v}}_0) + \mathbf{N}(\tilde{\mathbf{v}}_0, \mathbf{v}_s e^{i(\omega t + \varphi)} + \text{c.c.}) \rangle. \quad (28) $$

where $\mathbf{v}_s$ may be the vertical shear given by (6) or the horizontal shear given by (7). Finally, the terms $n_{1j}$ and $n_{2j}$ represent the nonlinear interactions between the resonant free Kelvin modes and the $j$-th forced mode of the base flow. They come from the dot products

$$ n_{lj} = \langle \mathbf{v}_l e^{i(\omega_l t + m_l \varphi)}, \mathbf{N}(\mathbf{v}_j e^{i(\omega t + \varphi)} + \text{c.c.}, \tilde{\mathbf{v}}_0) + \mathbf{N}(\tilde{\mathbf{v}}_0, \mathbf{v}_j e^{i(\omega t + \varphi)} + \text{c.c.}) \rangle. \quad (29) $$

Simplified expressions for all these terms are given in Appendix 3.6.2.

Seeking solutions to the amplitude equations (26) as growing exponentials, $A_j \sim e^{\sigma t}$, yields an analytical prediction for the complex growth rate $\sigma$ of the instability

$$ \sigma = |Ro|\sqrt{c_1 c_2}. \quad (30) $$

The temporal growth rate $\sigma_r$ is then obtained by taking the real part of $\sigma$, i.e. $\sigma_r = |Ro|\text{Re}(\sqrt{c_1 c_2})$. A resonant combination of free Kelvin modes is unstable if $\sigma_r$ is positive. This is always the case for combinations issued from the intersection of dispersion relations with opposite slopes Fukumoto (2003).
### 3.3.4 Amplitude equations using the horizontal shear decomposition

Amplitude equations (26) apply to any decomposition, $\mathbf{v}_{\text{base}}^{\text{V}}$ or $\mathbf{v}_{\text{base}}^{\text{H}}$, of the base flow. If the vertical shear decomposition $\mathbf{v}_{\text{base}}^{\text{V}}$ is chosen for the base flow, the formula for the growth rate contains infinite summations (corresponding to $n_{1j}$ and $n_{2j}$). This is because, in this case, axial wavenumbers $k_j^H$ are not integer multiples of $\pi/h$ and the integral from $z = -\pi/h$ to $\pi/h$ in the dot products (29) cannot be easily simplified.

However, if the horizontal decomposition $\mathbf{v}_{\text{base}}^{\text{H}}$ is chosen, it is possible to achieve such a simplification. In this case, the forced modes have axial wavenumbers $k_j^H$ that are odd multiples of $\pi/h$. These forced modes will interact with the resonant free Kelvin modes if the dot product between $\mathbf{N}(\mathbf{v}_j^H e^{i(\omega t + \varphi)} + \text{c.c.}, \tilde{\mathbf{v}}_0)$ and $\mathbf{v}_l e^{i(\omega_l t + m_l \varphi)}$ lead to a non-zero integral over $z$. Since these terms have the following $z$-Fourier components

$$ \mathbf{v}_l : z\text{-Fourier components } \pm k_l, \quad (31a) $$

$$ \mathbf{N}(\mathbf{v}_j^H e^{i(\omega t + \varphi)} + \text{c.c.}, \tilde{\mathbf{v}}_0) : z\text{-Fourier components } k_j^H \pm k_l, -k_j^H \pm k_l, \quad (31b) $$

and since $k_1$ and $k_2$ have different parities and $k_j^H$ is an odd multiple of $\pi/h$, the axial wavenumber of the dot product only contains even multiple of $\pi/h$, i.e. multiple of $2\pi/h$. As a consequence, the integral over $z$ is non-zero only if the axial wavenumber of the dot product is equal to zero. It follows that only two forced modes of $\mathbf{v}_{\text{base}}^{\text{H}}$ will interact with the resonant free Kelvin modes and their axial wavenumbers are

$$ k_j^H = k_{j1} = |k_2 - k_1| \quad \text{and} \quad k_j^H = k_{j2} = |k_2 + k_1|. \quad (32) $$

It means that there are only two forced modes $j_1$ and $j_2$ which give non zero coefficients $n_{1j}$ and $n_{2j}$ such that the infinite summations in (26) becomes sum of only two terms. The coefficients $c_1$ and $c_2$ thus simplify to

$$ c_1 = \frac{\overline{d_{12}} + n_{1s} + \overline{a_{j_1}^H} n_{1j_1} + \overline{a_{j_2}^H} n_{1j_2}}{\langle \mathbf{v}_1, \mathcal{I} \mathbf{v}_1 \rangle}, \quad (33a) $$

$$ c_2 = \frac{d_{21} + n_{2s} + \overline{a_{j_1}^H} n_{2j_1} + \overline{a_{j_2}^H} n_{2j_2}}{\langle \mathbf{v}_2, \mathcal{I} \mathbf{v}_2 \rangle}. \quad (33b) $$

### 3.3.5 Case of a resonant cylinder

In a resonant cylinder, the precession forces a flow that is dominated by a single forced mode, which is a Kelvin mode (i.e. a forced mode which satisfies the no-slip boundary conditions on the cylinder walls). In that case, the flow can become unstable if this forced Kelvin mode interacts constructively with two free Kelvin modes. As we have shown in a previous work, amplitude equations can be derived for a resonant cylinder Lagrange et al. (2011). In this section, we want to explain how to recover these amplitude equations from the amplitude equations of a non-resonant precessing fluid cylinder derived here in (26a,b).
When a forced mode is resonant, its amplitude $a_j$ predicted by the linear theory diverges. As shown by the equation (45), it happens when the dispersion relation $D(1, \omega, k_j) = 0$ holds for an axial wavenumber $k_j$ which is an odd multiple of $\pi/h$. For large Ekman numbers, the viscous effects saturate the amplitude of the forced Kelvin mode to an order $Ek^{-1/2}$ larger than the amplitudes of the shear flow and the others forced modes, see refs. (Gans, 1970; Meunier et al., 2008). Thus, at main order the base flow $v_{\text{base}}$ is a single forced Kelvin mode with amplitude $|\varepsilon| = O(RoEk^{-1/2})$. It follows that the summations in (26) are truncated to the index of that mode and the amplitude $a_j$ must be replaced by $\varepsilon/Ro$. Since the terms $d_{12}$, $n_{1s}$, $d_{21}$ and $n_{2s}$ are of order $O(1)$, they are negligible compared to $\varepsilon/Ro$ and can be dropped in (26). The amplitude equations when the j-th forced mode is resonant thus become

$$ \frac{dA_1}{d\tau} = A_2 \frac{\bar{\varepsilon} n_{1j}}{\langle \mathbf{v}_1, \mathcal{I} \mathbf{v}_1 \rangle} = n_1 \bar{\varepsilon} A_2, \quad (34a) $$

$$ \frac{dA_2}{d\tau} = A_1 \frac{\varepsilon n_{2j}}{\langle \mathbf{v}_2, \mathcal{I} \mathbf{v}_2 \rangle} = n_2 \varepsilon A_1. \quad (34b) $$

As explained in §3.3.4, terms $n_1$ and $n_2$ are non-zeros only if the condition of resonance $|k_2 - k_1| = k_j = (2j - 1)\pi/h$ or $|k_2 + k_1| = k_j = (2j - 1)\pi/h$ is satisfied (together with similar conditions on $m$ and $\omega$: $m_2 - m_1 = 1$ and $\omega_2 - \omega_1 = \omega$).

Seeking solutions to the amplitude equations (34a,b) as growing exponentials $A_j \sim e^{\sigma t}$, yields an expression for the growth rate

$$ \sigma = |\varepsilon| \sqrt{n_1 n_2}, \quad (35) $$

similar to the one obtained in Lagrange et al. (2011). As expected (based on similarities with the elliptic instability), the growth rate scales as the amplitude of the forced Kelvin mode. For large Ekman number, the growth rate of a resonant precessing fluid cylinder is thus an order $Ek^{-1/2}$ larger than the growth rate of the non-resonant case.

For example, for $h = 1.62$ and $\omega = 1.18$, the first Kelvin mode is forced at its first resonance and we recover (see Table 3 in Lagrange et al. (2011)) that the resonant combination (6, 1, 1) has $n_1 = -1.672$, $n_2 = -2.456$, leading to a growth rate $\sigma = 2.026|\varepsilon|$.

### 3.3.6 Viscous effects

Amplitude equations (26a,b) have been derived under the assumption of an inviscid fluid. Accounting for viscosity, they modify to Lagrange et al. (2011)

$$ \frac{dA_1}{d\tau} = c_1 A_2 - \alpha_1 A_1, \quad (36a) $$

$$ \frac{dA_2}{d\tau} = c_2 A_1 - \alpha_2 A_2, \quad (36b) $$

with $\alpha_l = s_l Ek^{1/2} + v_l Ek$. The coefficients $s_l$ represent the surface viscous damping of the free Kelvin modes due to Ekman layers on the cylinder walls.
They are complex numbers with a positive real part Lagrange et al. (2011). The coefficients $v_l$ are real and represent the volume viscous damping of the free Kelvin modes. They come from the Laplace operator of the Navier-Stokes equations, and are proportional to $k_l^2 + \delta_l^2$, such that they strongly damp the free Kelvin modes with large axial and radial wavenumbers Lagrange et al. (2011).

The critical Rossby number for which the instability appears can be determined from the condition of a vanishing growth rate, leading to

$$|Ro_{\text{crit}}| = \left\{ \frac{\alpha_1^r \alpha_2^r}{c_1 c_2} \left[ 1 + \left( \frac{\alpha_1^i - \alpha_2^i}{\alpha_1^r + \alpha_2^r} \right)^2 \right] \right\}^{1/2}, \quad (37)$$

where $\alpha_l^r$ and $\alpha_l^i$ are respectively the real and imaginary parts of $\alpha_l$. For large $Ek$, the volume viscous effects (which scale as $Ek$) are larger than the surface viscous effects (which scale as $Ek^{1/2}$), such that $Ro_{\text{crit}}$ scales as $Ek^{3/2}$. On the contrary, for low $Ek$, surface viscous effects are dominant and $Ro_{\text{crit}}$ scales as $Ek$ Lagrange et al. (2011).

## 3.4 Numerical application far from a resonance

In this section we compute the inviscid growth rate of the precessing instability, given by (30), when the forced modes are not resonant. To do so, we first determine an aspect ratio $h$ and a frequency ratio $\omega$, such that the first 4 forced modes are away from a resonance, and such that combination of free Kelvin modes ($m_2 = 3, l_1 = 1, l_2 = 1$) forms a perfect triadic resonance with the first forced mode (see Fig. 3.6). These two conditions are illustrated in Fig. 3.7: the first 5 resonances of the first 4 forced modes are noted $\omega_{j,n}$ (n-th resonance of the j-th forced Kelvin mode) and are solutions to the dispersion relation

$$D \left( m_j, \omega_{j,n}, \frac{(4 - \omega_{j,n}^2)^{1/2}}{|\omega_{j,n}|} (2n - 1) \frac{\pi}{h} \right) = 0. \quad (38)$$

The thick blue curve in Fig. 3.7 gives the aspect ratio and the frequency ratio for which the combination (3, 1, 1) is resonant. This curve is the solution of the equation

$$\omega_2(k_2 = 2\pi/h, l_2 = 1, m_2 = 3) = \omega_1(k_1 = \pi/h, l_1 = 1, m_{1r} = 2) + \omega_{i,n}, \quad (39)$$

where $\omega_2(k_2 = 2\pi/h, l_2 = 1, m_2 = 3)$ means: the value of $\omega_2$ for $k_2 = 2\pi/h$ and $l_2 = 1$ (first branch of the dispersion relation $m_2 = 3$).

In the present numerical investigations, we choose $h = 2.3$ and $\omega = 1.34$. It can be checked that for these aspect and frequency ratios, (3, 1, 1) is a resonant point (Fig. 3.6). The theoretical axial vorticity fields of the corresponding free Kelvin modes, $m_1 = 2$ and $m_2 = 3$, are shown in Fig. 3.8. Since these modes correspond to the first branches of the dispersion relation ($l_1 = l_2 = 1$), their vorticity fields exhibit only one ring of 4 (for $m_1 = 2$) or 6 (for $m_2 = 3$) counter-rotating vortices.
<figure><img src="image_8.png" /><figcaption>Figure 3.7: Resonance location in the $\omega-h$ plane. The resonances of the 1st (black curves), 2nd (red curves), 3rd (pink curves), 4th (green curves) forced modes are shown. The notation $\omega_{j,n}$ refers to the $n$-th resonance of the $j$-th forced mode. The blue thick line indicates the values of $\omega$ and $h$ for which the unstable modes combination (3, 1, 1) forms a perfect triadic resonance with the first forced mode (see Fig. 3.6). The black circle ($\omega = 1.34, h = 2.3$) is the choice for the numerical application.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 3.8: Theoretical axial vorticity fields of the free Kelvin modes $m_1 = 2$ (a) and $m_2 = 3$ (b) which are expected to grow for $h = 2.3$ and $\omega = 1.34$.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 3.9: Stability diagram of the flow inside a precessing cylinder, for $h = 2.30$ and $\omega = 1.34$. The stable and unstable domains are separated by the solid line corresponding to the prediction (37).</figcaption></figure>

<table><caption>Table 3.1: Values of the parameters appearing in the amplitude equations. The aspect ratio and the frequency ratio are $h = 2.30$ and $\omega = 1.34$. For these values, the combination ($m_2 = 3, l_1 = 1, l_2 = 1$) is resonant and corresponds to a pair of free Kelvin modes ($m_1 = 2, \omega_1 = -0.466, k_1 = \pi/h$) and ($m_2 = 3, \omega_2 = 0.874, k_2 = 2\pi/h$). The forced modes are those with indices $j_1 = 1, j_2 = 2$ and have axial wavenumbers $k_{j1} = \pi/h$ and $k_{j2} = 3\pi/h$. Their amplitudes are $a_{j1}^H = 1.4983$ and $a_{j2}^H = 0.1130$.</caption><tbody><tr><td>⟨v<sub>1</sub>, Iv<sub>1</sub>⟩</td><td>d&#773;<sub>12</sub></td><td>n<sub>1s</sub></td><td>n<sub>1j1</sub></td><td>n<sub>1j2</sub></td></tr><tr><td>27.1230</td><td>6.8448</td><td>1.1623</td><td>-12.7524</td><td>3.6628</td></tr><tr><td></td><td></td><td>c<sub>1</sub></td><td>s<sub>1</sub></td><td>v<sub>1</sub></td></tr><tr><td></td><td></td><td>-0.3940</td><td>1.22 - 0.15i</td><td>34.37</td></tr><tr><td>⟨v<sub>2</sub>, Iv<sub>2</sub>⟩</td><td>d<sub>21</sub></td><td>n<sub>2s</sub></td><td>n<sub>2j1</sub></td><td>n<sub>2j2</sub></td></tr><tr><td>35.2671</td><td>-6.8448</td><td>21.8624</td><td>-23.9177</td><td>6.8697</td></tr><tr><td></td><td></td><td>c<sub>2</sub></td><td>s<sub>2</sub></td><td>v<sub>2</sub></td></tr><tr><td></td><td></td><td>-0.5683</td><td>1.50 + 0.027i</td><td>39.08</td></tr></tbody></table>
The values needed to compute the inviscid growth rate given by (30) are listed in Table 5.1. For $h = 2.3$ and $\omega = 1.34$, we obtain a growth rate $\sigma_r = |Ro|(\sqrt{c_1c_2}) = 0.4732|Ro|$. Since the combination (3, 1, 1) corresponds to free Kelvin modes with simple radial and axial structures, volume viscous damping (which scale as $k_l^2 + \delta_l^2$) is small, which make these modes the perfect candidates for an instability. The stability diagram of the resonant combination (3, 1, 1) is shown in Fig. 3.9 in the *Ek–Ro* plane: in this graph, the prediction from (37) is represented by a solid line that splits the *Ek–Ro* plane into a stable and an unstable domain.

## 3.5 Conclusion

In this paper, the instability of a fluid inside a precessing cylinder has been addressed theoretically. First, we have shown that the base flow can be written as a superposition of a vertical or horizontal shear flow and a sum of forced modes. The stability of this base flow has then been studied for a forcing at a non-resonant frequency, thus completing previous studies performed for resonant or near-resonant frequencies Lagrange *et al.* (2008, 2011). We have shown that the non-resonant base flow can trigger a triadic instability with two free Kelvin modes only if these modes have different axial parities. We also have obtained a prediction of the instability growth rate, showing that the inviscid growth rate is proportional to the Rossby number, *Ro*, and an order $Ek^{-1/2}$ smaller than the growth rate obtained for a resonant base flow. Finally, introducing viscous damping, we have given a prediction of *Ro*<sub>crit</sub>, the critical Rossby number for which the flow becomes unstable, as a function of *Ek*: for large (resp. low) *Ek* numbers, *Ro*<sub>crit</sub> scales as $Ek^{3/2}$ (resp. *Ek*).

We believe that the theoretical predictions provided in this paper will be valuable for future experimental and numerical studies of precessing flows performed at arbitrary frequencies. We shall note however that our analysis relies on the assumption of a small Rossby number, i.e. a small precession angle. For strong forcing, very different phenomena are expected (e.g. Kelvin-Helmholtz, centrifugal, or boundary layer instabilities), due to the generation of powerful zonal flows.

In closing, the precessional instability is typical of transition to turbulence in rotating flows. The presence of rotation ensures that energy is continuously provided to the flow. It also supports the existence of inertial waves that can lead to several instabilities (elliptic instability, libration instability, etc.). The structure of turbulence is also modified by the presence of the rotation because of the anisotropy induced. There is much more work to be done on this fascinating topic if we want to understand the mechanisms at play in turbulent rotating flows.

This study was carried out under CEA-CNRS contract No. 012171.
## 3.6 Appendix

### 3.6.1 Operators of the Euler equations, forced mode vector $\mathbf{v}_j$, amplitudes $a_j^V$ and $a_j^H$, free Kelvin modes vectors $\mathbf{v}_i^+$ and $\mathbf{v}_i^-$

Operators used for the matrix formulation (3) of the Euler equations are

$$ 
\begin{aligned}
\mathcal{I} &= \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}, \quad \mathcal{D} = \begin{pmatrix} 0 & 0 & -i & 0 \\ 0 & 0 & 1 & 0 \\ i & -1 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix}, \\
\mathcal{M} &= \begin{pmatrix} 0 & -2 & 0 & \frac{\partial}{\partial r} \\ 2 & 0 & 0 & \frac{1}{r} \frac{\partial}{\partial \varphi} \\ 0 & 0 & 0 & \frac{\partial}{\partial z} \\ \frac{\partial}{\partial r} + \frac{1}{r} & \frac{1}{r} \frac{\partial}{\partial \varphi} & \frac{\partial}{\partial z} & 0 \end{pmatrix},
\end{aligned} 
\quad (40) $$

and

$$ \mathbf{F}_0 = \begin{pmatrix} 0 \\ 0 \\ -r\omega \\ 0 \end{pmatrix}, \quad \mathbf{N}(\mathbf{v}_1, \mathbf{v}_2) = \begin{pmatrix} \mathbf{u}_1 \times (\nabla \times \mathbf{u}_2) \\ 0 \end{pmatrix}. \quad (41) $$

The forced modes quadri-vector $\mathbf{v}_j^V$ and $\mathbf{v}_j^H$ are given by the same expression

$$ \mathbf{v}_j^V = \mathbf{v}_j^H = 2 \begin{pmatrix} (\mathbf{w}_{1,\omega,k_j})_r i \sin(k_j z) \\ (\mathbf{w}_{1,\omega,k_j})_\varphi i \sin(k_j z) \\ (\mathbf{w}_{1,\omega,k_j})_z \cos(k_j z) \\ (\mathbf{w}_{1,\omega,k_j})_p i \sin(k_j z) \end{pmatrix}, \quad (42) $$

with

$$ \mathbf{w}_{m,\omega,k} = \begin{pmatrix} \frac{-1}{4-\omega^2} \left( \omega \delta J'_m(\delta r) + 2 \frac{m}{r} J_m(\delta r) \right) \\ \frac{-i}{4-\omega^2} \left( 2 \delta J'_m(\delta r) + \frac{\omega m}{r} J_m(\delta r) \right) \\ i \frac{k}{\omega} J_m(\delta r) \\ -i J_m(\delta r) \end{pmatrix}, \quad (43) $$

and where $J_m$ is the Bessel function of the first kind, $J'_m$ its derivative and $\delta = \sqrt{4 - \omega^2 k/|\omega|}$. However, they have been named differently since they have different boundary conditions and thus different axial wavenumbers $k_j^V$ or $k_j^H$. For the vertical shear, the radial wavenumber $\delta$ is imposed by the boundary conditions and $k$ is deduced from the dispersion relation whereas for the horizontal shear, the axial wavenumber is imposed by the boundary conditions and $\delta$ is given again by $\delta = \sqrt{4 - \omega^2 k/|\omega|}$.

The dispersion relation is

$$ D(m, \omega, k) = \omega \delta J'_m(\delta) + 2m J_m(\delta). \quad (44) $$

The amplitudes of the forced modes are

$$ a_j^V = \frac{\omega^2}{(\omega - 2)(k_j^2 + 1)k_j J_1(\delta_j) \cos(k_j h/2)}, \quad (45a) $$

$$ a_j^H = - \frac{2h(-1)^j \omega(2 + \omega)}{\pi^2(2j - 1)^2(\omega \delta_j J'_1(\delta_j) + 2J_1(\delta_j))}. \quad (45b) $$
Vectors $\mathbf{v}_l^+$ and $\mathbf{v}_l^-$ appearing in (20) are

$$ \mathbf{v}_l^+ = 2 \begin{pmatrix} \left( \mathbf{w}_{m_l, \omega_l, k_l^+} \right)_r \cos(k_l^+ z) \\ \left( \mathbf{w}_{m_l, \omega_l, k_l^+} \right)_\varphi \cos(k_l^+ z) \\ \left( \mathbf{w}_{m_l, \omega_l, k_l^+} \right)_z \sin(k_l^+ z) \\ \left( \mathbf{w}_{m_l, \omega_l, k_l^+} \right)_p \cos(k_l^+ z) \end{pmatrix}, \quad \mathbf{v}_l^- = 2 \begin{pmatrix} \left( \mathbf{w}_{m_l, \omega_l, k_l^-} \right)_r \sin(k_l^- z) \\ \left( \mathbf{w}_{m_l, \omega_l, k_l^-} \right)_\varphi \sin(k_l^- z) \\ \left( \mathbf{w}_{m_l, \omega_l, k_l^-} \right)_z \cos(k_l^- z) \\ \left( \mathbf{w}_{m_l, \omega_l, k_l^-} \right)_p \sin(k_l^- z) \end{pmatrix}. \quad (46) $$

### 3.6.2 Derivation of amplitude equations

In this Appendix, we derive the amplitude equations (26), starting from the order $R_0$ equation (25) that we report here

$$ \left( \frac{\partial \mathcal{I}}{\partial t} + \mathcal{M} \right) \tilde{\mathbf{v}}_1 = \mathbf{N}(\mathbf{v}_{\text{base}}, \tilde{\mathbf{v}}_0) + \mathbf{N}(\tilde{\mathbf{v}}_0, \mathbf{v}_{\text{base}}) + \left[ \left( \mathcal{D}e^{i(\omega t + \varphi)} + \text{c.c.} \right) - \frac{\partial \mathcal{I}}{\partial \tau} \right] \tilde{\mathbf{v}}_0. \quad (47) $$

As explained in the core of the manuscript, a solvability condition is obtained by taking the dot product of this equation with $\mathbf{v}_l e^{i(\omega_l t + m_l \varphi)}$, $l = 1, 2$. The problem being self adjoint, we have $\langle \mathbf{v}_l e^{i(\omega_l t + m_l \varphi)}, (\partial \mathcal{I}/\partial t + \mathcal{M}) \tilde{\mathbf{v}}_1 \rangle = 0$, so that, we are left with

$$ \begin{aligned} \left\langle \mathbf{v}_l e^{i(\omega_l t + m_l \varphi)}, \frac{\partial \mathcal{I}}{\partial \tau} \tilde{\mathbf{v}}_0 \right\rangle = & \left\langle \mathbf{v}_l e^{i(\omega_l t + m_l \varphi)}, \left( \mathcal{D} e^{i(\omega t + \varphi)} + \text{c.c.} \right) \tilde{\mathbf{v}}_0 \right\rangle \\ & + \left\langle \mathbf{v}_l e^{i(\omega_l t + m_l \varphi)}, \mathbf{N}(\mathbf{v}_{\text{base}}, \tilde{\mathbf{v}}_0) \right\rangle \\ & + \left\langle \mathbf{v}_l e^{i(\omega_l t + m_l \varphi)}, \mathbf{N}(\tilde{\mathbf{v}}_0, \mathbf{v}_{\text{base}}) \right\rangle. \end{aligned} \quad (48) $$

The computation of the LHS term is straightforward and follows from the linearity of the dot product and the orthogonality of the free Kelvin modes.
Introducing the expression of $\tilde{\mathbf{v}}_0$ given by (21) gives a LHS term equal to

$$
\begin{align}
& \left\langle \mathbf{v}_l e^{i(\omega_l t + m_l \varphi)}, \sum_{l=1}^{2} \frac{dA_j}{d\tau} \mathcal{I} \left( \mathbf{v}_j e^{i(\omega_j t + m_j \varphi)} \right) \right\rangle \nonumber \\
& + \left\langle \mathbf{v}_l e^{i(\omega_l t + m_l \varphi)}, \sum_{j=1}^{2} \frac{d\overline{A_j}}{d\tau} \mathcal{I} \left( \overline{\mathbf{v}_j} e^{-i(\omega_j t + m_j \varphi)} \right) \right\rangle \nonumber \\
& = \sum_{j=1}^{2} \frac{dA_j}{d\tau} \left\langle \mathbf{v}_l e^{i(\omega_l t + m_l \varphi)}, e^{i(\omega_j t + m_j \varphi)} \mathcal{I} \mathbf{v}_j \right\rangle \nonumber \\
& + \sum_{j=1}^{2} \frac{d\overline{A_j}}{d\tau} \left\langle \mathbf{v}_l e^{i(\omega_l t + m_l \varphi)}, e^{-i(\omega_j t + m_j \varphi)} \mathcal{I} \overline{\mathbf{v}_j} \right\rangle, \tag{49a}
\end{align}
$$

$$
\begin{align}
&= \sum_{j=1}^{2} \frac{dA_j}{d\tau} \underbrace{\left\langle \mathbf{v}_l, e^{i[(-\omega_l + \omega_j)t + (-m_l + m_j)\varphi]} \mathcal{I} \mathbf{v}_j \right\rangle}_{\neq 0 \text{ if } j=l} \nonumber \\
&\quad + \sum_{j=1}^{2} \frac{d\overline{A_j}}{d\tau} \underbrace{\left\langle \mathbf{v}_l, e^{i[(-\omega_l - \omega_j)t + (-m_l - m_j)\varphi]} \mathcal{I} \overline{\mathbf{v}_j} \right\rangle}_{0 \text{ because } \propto \int_0^{2\pi} e^{-i(m_l+m_j)\varphi} d\varphi} \tag{49b}
\end{align}
$$

$$
= \frac{dA_l}{d\tau} \langle \mathbf{v}_l, \mathcal{I} \mathbf{v}_l \rangle . \tag{49c}
$$

Eq. (49b) shows that it is not necessary to take into account the c.c. part of $A_j \mathbf{v}_j e^{i(\omega_j t+m_j\varphi)}$ in (48) since it leads to 0 integral terms. This observation still holds for computations with operators $\mathcal{D}$, $\overline{\mathcal{D}}$ and $\mathbf{N}$ since they do not change the wavenumbers in the exponential when applied to $\mathbf{v}_l e^{i(\omega_l t+m_l\varphi)}$. Therefore, the c.c. part of $A_j \mathbf{v}_j e^{i(\omega_j t+m_j\varphi)}$ will be omitted in the next computations.

Plugging (49c) into (48) yields the amplitude equations

$$
\begin{align}
\frac{dA_l}{d\tau} &= \frac{\langle \mathbf{v}_l e^{i(\omega_l t + m_l \varphi)}, (\mathcal{D} e^{i(\omega t + \varphi)} + \text{c.c.}) \tilde{\mathbf{v}}_0 \rangle}{\langle \mathbf{v}_l, \mathcal{I} \mathbf{v}_l \rangle} \nonumber \\
&\quad + \frac{\langle \mathbf{v}_l e^{i(\omega_l t + m_l \varphi)}, \mathbf{N}(\mathbf{v}_{\text{base}}, \tilde{\mathbf{v}}_0) + \mathbf{N}(\tilde{\mathbf{v}}_0, \mathbf{v}_{\text{base}}) \rangle}{\langle \mathbf{v}_l, \mathcal{I} \mathbf{v}_l \rangle}. \tag{50}
\end{align}
$$

We now proceed with the calculation of the RHS terms of (50). Computa-
tions are performed with $l = 1$ so that results for $l = 2$ will follow from the
permutation of indices $(1, 2) \rightarrow (2, 1)$. Also, from the operators properties
presented in §3.2.2, we know that the free Kelvin modes $\mathbf{v}_1 e^{i(\omega_1 t+m_1\varphi)}$ and
$\mathbf{v}_2 e^{i(\omega_2 t+m_2\varphi)}$ must have different z-parities to give nonzero coupling terms.
Thus, for $l = 1$ we can directly substitute in (50) vector $\tilde{\mathbf{v}}_0$ by $\mathbf{v}_2 e^{i(\omega_2 t+m_2\varphi)}$
and (50) writes

$$ \frac{dA_1}{d\tau} = \frac{\langle \mathbf{v}_1 e^{i(\omega_1 t+m_1 \varphi)}, (\mathcal{D} e^{i(\omega t+\varphi)} + \text{c.c.}) \mathbf{v}_2 e^{i(\omega_2 t+m_2 \varphi)} \rangle}{\langle \mathbf{v}_1, \mathcal{I} \mathbf{v}_1 \rangle} + \frac{\langle \mathbf{v}_1 e^{i(\omega_1 t+m_1 \varphi)}, \mathbf{N} (\mathbf{v}_{\text{base}}, \mathbf{v}_2 e^{i(\omega_2 t+m_2 \varphi)}) \rangle}{\langle \mathbf{v}_1, \mathcal{I} \mathbf{v}_1 \rangle} + \frac{\langle \mathbf{v}_1 e^{i(\omega_1 t+m_1 \varphi)}, \mathbf{N} (\mathbf{v}_2 e^{i(\omega_2 t+m_2 \varphi)}, \mathbf{v}_{\text{base}}) \rangle}{\langle \mathbf{v}_1, \mathcal{I} \mathbf{v}_1 \rangle}. \quad (51) $$

In what follows we compute each of the terms in the RHS of (51) and we assume that the resonance conditions (24) are fulfilled in order to drop the exponential terms.

$$ \text{3.6.3 Computation of} \quad \langle \mathbf{v}_1 e^{i(\omega_1 t+m_1 \varphi)}, (\mathcal{D} e^{i(\omega t+\varphi)} + \text{c.c.}) (\mathbf{v}_2 e^{i(\omega_2 t+m_2 \varphi)}) \rangle $$

Expanding the complex conjugate leads to 2 terms

$$ \langle \mathbf{v}_1 e^{i(\omega_1 t+m_1 \varphi)}, e^{i(\omega t+\varphi)} e^{i(\omega_2 t+m_2 \varphi)} \mathcal{D} \mathbf{v}_2 \rangle \langle \mathbf{v}_1 e^{i(\omega_1 t+m_1 \varphi)}, e^{-i(\omega t+\varphi)} e^{i(\omega_2 t+m_2 \varphi)} \bar{\mathcal{D}} \mathbf{v}_2 \rangle. \quad (52) $$

The first term vanishes because the azimuthal Fourier components are different on each side of the dot product such that the integral over $\varphi$ gives zero. In contrast, in the second term, the azimuthal Fourier components are equal and can thus be dropped. This term can thus be written as

$$ \bar{d}_{12} = \langle \mathbf{v}_1, \bar{\mathcal{D}} \mathbf{v}_2 \rangle. \quad (53) $$

$$ \text{3.6.4 Computation of } \langle \mathbf{v}_1 e^{i(\omega_1 t+m_1 \varphi)}, \mathbf{N} (\mathbf{v}_{\text{base}}, \mathbf{v}_2 e^{i(\omega_2 t+m_2 \varphi)}) \rangle $$

Here $\mathbf{v}_{\text{base}} = \mathbf{v}_b e^{i(\omega t+\varphi)} + \text{c.c.}$ is the base flow given either by (8) or (9) depending on which decomposition (vertical or horizontal shear) is used to express the base flow. We have

$$ \langle \mathbf{v}_1 e^{i(\omega_1 t+m_1 \varphi)}, \mathbf{N} (\mathbf{v}_{\text{base}}, \mathbf{v}_2 e^{i(\omega_2 t+m_2 \varphi)}) \rangle = \langle \mathbf{v}_1 e^{i(\omega_1 t+m_1 \varphi)}, \mathbf{N} (\mathbf{v}_b e^{i(\omega t+\varphi)} + \text{c.c.}, \mathbf{v}_2 e^{i(\omega_2 t+m_2 \varphi)}) \rangle, \quad (54a) $$

$$ = \langle \mathbf{v}_1 e^{i(\omega_1 t+m_1 \varphi)}, \mathbf{N} (\mathbf{v}_b e^{i(\omega t+\varphi)}, \mathbf{v}_2 e^{i(\omega_2 t+m_2 \varphi)}) \rangle + \langle \mathbf{v}_1 e^{i(\omega_1 t+m_1 \varphi)}, \mathbf{N} (\bar{\mathbf{v}}_b e^{-i(\omega t+\varphi)}, \mathbf{v}_2 e^{i(\omega_2 t+m_2 \varphi)}) \rangle, \quad (54b) $$

$$ = \langle \mathbf{v}_1 e^{i(\omega_1 t+m_1 \varphi)}, e^{i(\omega t+\varphi)} e^{i(\omega_2 t+m_2 \varphi)} \mathbf{N}_{im_2} (\mathbf{v}_b, \mathbf{v}_2) \rangle + \langle \mathbf{v}_1 e^{i(\omega_1 t+m_1 \varphi)}, e^{-i(\omega t+\varphi)} e^{i(\omega_2 t+m_2 \varphi)} \mathbf{N}_{im_2} (\bar{\mathbf{v}}_b, \mathbf{v}_2) \rangle, \quad (54c) $$

where $\mathbf{N}_{im_2}$ corresponds to operator $\mathbf{N}$ where $d/d\varphi$ has been replaced by $im_2$. As previously, the first term vanishes and the exponential can be dropped
from the second term. Introducing the expression of $\mathbf{v}_b$ makes this whole
term equal to

$$
\left\langle \mathbf{v}_1, \mathbf{N}_{im2} \left( \overline{\mathbf{v}_s} + \sum_{j=1}^{\infty} \overline{a_j \mathbf{v}_j}, \mathbf{v}_2 \right) \right\rangle = n_{1\bar{s}2} + \sum_{j=1}^{\infty} \overline{a_j} n_{1\bar{j}2}, \qquad (55a)
$$

$$
\begin{equation}
\begin{split}
& \text{with } n_{1\bar{s}2} = \langle \mathbf{v}_1, \mathbf{N}_{im2} (\overline{\mathbf{v}}_s, \mathbf{v}_2) \rangle \quad \text{and} \quad n_{1\bar{j}2} = \langle \mathbf{v}_1, \mathbf{N}_{im2} (\overline{\mathbf{v}}_j, \mathbf{v}_2) \rangle, \\
& \hfill (55\text{b})
\end{split}
\end{equation}
$$

3.6.5 Computation of $\langle \mathbf{v}_1 e^{i(\omega_1 t + m_1 \varphi)}, \mathbf{N} (\mathbf{v}_2 e^{i(\omega_2 t + m_2 \varphi)}, \mathbf{v}_{\text{base.}}) \rangle$

We have

$$
\begin{align}
& \left\langle \mathbf{v}_1 e^{i(\omega_1 t + m_1 \varphi)}, \mathbf{N} \left( \mathbf{v}_2 e^{i(\omega_2 t + m_2 \varphi)}, \mathbf{v}_{\text{base}} \right) \right\rangle \nonumber \\
&= \left\langle \mathbf{v}_1 e^{i(\omega_1 t + m_1 \varphi)}, \mathbf{N} \left( \mathbf{v}_2 e^{i(\omega_2 t + m_2 \varphi)}, \mathbf{v}_b e^{i(\omega t + \varphi)} + \text{c.c.} \right) \right\rangle, \tag{56a}
\end{align}
$$

$$
= \left\langle \mathbf{v}_1 e^{i(\omega_1 t + m_1 \varphi)}, \mathbf{N} \left( \mathbf{v}_2 e^{i(\omega_2 t + m_2 \varphi)}, \mathbf{v}_b e^{i(\omega t + \varphi)} \right) \right\rangle \\
+ \left\langle \mathbf{v}_1 e^{i(\omega_1 t + m_1 \varphi)}, \mathbf{N} \left( \mathbf{v}_2 e^{i(\omega_2 t + m_2 \varphi)}, \overline{\mathbf{v}_b} e^{-i(\omega t + \varphi)} \right) \right\rangle, \quad (56\text{b})
$$

$$
= \left\langle \mathbf{v}_1 e^{i(\omega_1 t + m_1 \varphi)}, e^{i(\omega_2 t + m_2 \varphi)} e^{i(\omega t + \varphi)} \mathbf{N}_i (\mathbf{v}_2, \mathbf{v}_b) \right\rangle \\
+ \left\langle \mathbf{v}_1 e^{i(\omega_1 t + m_1 \varphi)}, e^{i(\omega_2 t + m_2 \varphi)} e^{-i(\omega t + \varphi)} \mathbf{N}_{-i} (\mathbf{v}_2, \overline{\mathbf{v}}_b) \right\rangle, \quad (56c)
$$

where $\mathbf{N}_i$ and $\mathbf{N}_{-i}$ correspond to operator $\mathbf{N}$ where $d/d\varphi$ has been replaced by $i$ and $-i$, respectively. As previously, the first term vanishes and the ex-
ponentials can be dropped from the second term. Introducing the expression
of $\mathbf{v}_b$ makes this whole term equal to

$$
\left\langle \mathbf{v}_1, \mathbf{N}_{-i} \left( \mathbf{v}_2, \overline{\mathbf{v}}_s + \sum_{j=1}^{\infty} \overline{a_j \mathbf{v}_j} \right) \right\rangle = n_{12\bar{s}} + \sum_{j=1}^{\infty} \overline{a_j} n_{12\bar{j}}, \quad (57\text{a})
$$

$$
\begin{equation}
\begin{split}
& \text{with } n_{12\bar{s}} = \langle \mathbf{v}_1, \mathbf{N}_{-i} (\mathbf{v}_2, \overline{\mathbf{v}}_s) \rangle \quad \text{and} \quad n_{12\bar{j}} = \langle \mathbf{v}_1, \mathbf{N}_{-i} (\mathbf{v}_2, \overline{\mathbf{v}}_j) \rangle \\
& \hfill (57\text{b})
\end{split}
\end{equation}
$$

3.6.6 Conditions of resonance and amplitude equations

Collecting (53), (55b) and (57b) together, the amplitude equations (50)
rewrite

$$
\frac{dA_1}{d\tau} = A_2 \frac{\bar{d}_{12} + n_{1s} + \sum_{j=1}^{\infty} \bar{a}_j n_{1j}}{\langle \mathbf{v}_1, \mathcal{I} \mathbf{v}_1 \rangle}, \qquad (58a)
$$

$$
\frac{dA_2}{d\tau} = A_1 \frac{d_{21} + n_{2s} + \sum_{j=1}^{\infty} a_j n_{2j}}{\langle \mathbf{v}_2, \mathcal{I} \mathbf{v}_2 \rangle}, \quad (58\text{b})
$$

where coefficients in (58a) are

$$
\bar{d}_{12} = \langle \mathbf{v}_1, \bar{D} \mathbf{v}_2 \rangle, \qquad (59a)
$$

$$
n_{1s} = n_{1\bar{s}2} + n_{12\bar{s}} = \langle \mathbf{v}_1, \mathbf{N}_{im2} (\overline{\mathbf{v}}_s, \mathbf{v}_2) \rangle + \langle \mathbf{v}_1, \mathbf{N}_{-i} (\mathbf{v}_2, \overline{\mathbf{v}}_s) \rangle, \quad (59\text{b})
$$

$$
n_{1j} = n_{1\bar{j}2} + n_{12\bar{j}} = \langle \mathbf{v}_1, \mathbf{N}_{im2} (\overline{\mathbf{v}}_j, \mathbf{v}_2) \rangle + \langle \mathbf{v}_1, \mathbf{N}_{-i} (\mathbf{v}_2, \overline{\mathbf{v}}_j) \rangle. \quad (59c)
$$
Coefficients in (58b) are

$$d_{21} = \langle \mathbf{v}_2, \mathcal{D}\mathbf{v}_1 \rangle = -\bar{d}_{12}, \qquad (60a)$$

$$n_{2s} = n_{2s1} + n_{21s} = \langle \mathbf{v}_2, \mathbf{N}_{im_1}(\mathbf{v}_s, \mathbf{v}_1) \rangle + \langle \mathbf{v}_2, \mathbf{N}_i(\mathbf{v}_1, \mathbf{v}_s) \rangle, \qquad (60b)$$

$$n_{2j} = n_{2j1} + n_{21j} = \langle \mathbf{v}_2, \mathbf{N}_{im_1}(\mathbf{v}_j, \mathbf{v}_1) \rangle + \langle \mathbf{v}_2, \mathbf{N}_i(\mathbf{v}_1, \mathbf{v}_j) \rangle. \qquad (60c)$$
Bibliography

AGRAWAL, B. N. 1993 Dynamics characteristics of liquid motion in partially filled tanks of a spinning spacecraft. *J. Guid. Control Dynam.* **16** (4), 636–640.

BAO, G. W. & PASCAL, M. 1997 Stability of a spinning liquid filled spacecraft. *Appl. Mech.* **67**, 407–421.

CÉBRON, D., BARS, M. L. & MAUBERT, P. 2011 Magnetohydrodynamic simulations of the elliptical instability in triaxial ellipsoids. *Geophys. and Astrophys. Fluid Dyn.* **106**, 524–546.

FUKUMOTO, Y. 2003 The three dimensional instability of a strained vortex tube revisited. *J. Fluid Mech.* **493**, 287–318.

FULTZ, D. 1959 A note on overstability and elastoid-inertia oscillations of kelvin, solberg and bjerknes. *J. Meteorol.* **16**, 199–208.

GAILITIS, A., LIELAUSSIS, O., PLATACIS, E., STEFANI, F. & GERBETH, G. 2002 Laboratory experiments on hydromagnetic dynamos. *Rev. Mod. Phys.* **74**, 973.

GANS, R. F. 1970 On the precession of a resonant cylinder. *J. Fluid Mech.* **476**, 865–872.

GANS, R. F. 1984 Dynamics of a near-resonant fluid-filled gyroscope. *AIAA J.* **22**, 1465–1471.

GARG, S. C., FURUNOTO, N. & VANYO, J. P. 1986 Spacecraft nutational instability prediction by energy dissipation measurments. *J. Guid.* **9** (3), 357–361.

GREENSPAN, H. P. 1968 *The theory of rotating fluids*. Cambridge University Press.

GROTE, F. & BUSSE, H. 2001 Dynamics of convection and dynamos in rotating spherical fluid shells. *Fluid Dyn. Res.* **28**, 349–368.

GROTE, F., BUSSE, H. & TILGNER, A. 1999 Convection driven quadrupolar dynamos in rotating spherical shells. *Phys. Rev. E* **60**, R5025.

GROTE, F., BUSSE, H. & TILGNER, A. 2000 Regular and chaotic spherical dynamos. *Phys. Earth Planet. Inter.* **117**, 259.

HERBERT, T. 1986 Viscous fluid motion in a spinning and nutating cylinder. *J. Fluid Mech.* **167**, 181–198.

KARPOV, B. G. 1965 *Ballistic Research Labs. Maryland, U. S.*, Report BRL R 1302.

KELVIN, L. 1880 Vibrations of a columnar vortex. *Phil. Mag.* **10**, 155–168.
KERSWELL, R. R. 1996 Upper bounds on the energy dissipation in turbulent precession, *J. Fluid Mech.* **321**, 335–370.

KERSWELL, R. R. & BARENGHI, C. F. 1995 On the viscous decay rates of inertial waves in a rotating cylinder. *J. Fluid Mech.* **285**, 203–214.

KOBINE, J. J. 1995 Inertial wave dynamics in a rotating and precessing cylinder. *J. Fluid Mech.* **303**, 233–252.

KONG, D., CUI, Z., LIAO, X. & ZHANG, K. 2015 On the transition from the laminar to disordered flow in a precessing spherical-like cylinder. *Geophys. Astrophys. Fluid Dyn.* **109**, 62–83.

LAGRANGE, R., ELOY, C., NADAL, F. & MEUNIER, P. 2008 Instability of a fluid inside a precessing cylinder. *Physics of Fluids.* **20**(8), 081701.

LAGRANGE, R., MEUNIER, P., ELOY, C. & NADAL, F. 2009 Dynamics of a fluid inside a precessing cylinder. *Mechanics and Industry.* **10**, 187–194.

LAGRANGE, R., MEUNIER, P., NADAL, F. & ELOY, C. 2011 Precessional instability of a fluid cylinder. *Journal of Fluid Mechanics* **666**, 104–145.

LAMBELIN, J. P., NADAL, F., LAGRANGE, R. & SARTHOU, A. 2009 Non-resonant viscous theory for the stability of a fluid-filled gyroscope. *Journal of Fluid Mechanics* **639**, 167–194.

LOPER, D. E. 1975 Torque balance and energy budget for the precessionally driven dynamo. *Phys. Earth Planet. Inter.* **43**, 43.

MALKUS, W. V. R. 1989 An experimental study of global instabilities due to tidal (elliptical) distortion of a rotating elastic cylinder. *Geophys. Astrophys. Fluid Dynamics* **48**, 123–134.

MANASSEH, R. 1992 Breakdown regimes of inertia waves in a precessing cylinder. *J. Fluid Mech.* **243**, 261–296.

MANASSEH, R. 1994 Distorsions of inertia waves in a precessing cylinder forced near its fundamental mode resonance. *J. Fluid Mech.* **265**, 345–370.

MANASSEH, R. 1996 Nonlinear behaviour of contained inertia waves. *J. Fluid Mech.* **315**, 151–173.

MCEWAN, A. D. 1970 Inertial oscillations in a rotating fluid cylinder. *J. Fluid Mech.* **40** (3), 603–640.

MEUNIER, P., ELOY, C., LAGRANGE, R. & NADAL, F. 2008 A rotating fluid cylinder subject to weak precession. *J. Fluid Mech.* **599**, 405–440.

NORE, C., LÉORAT, J., GUERMOND, J. & LUDDENS, F. 2011 Nonlinear dynamo action in a precessing cylindrical container. *Phys. Rev. E* **84**, 016317.
POCHA, J. J. 1987 An experimental investigation of spacecraft sloshing. *Space Commun. Broadcasting* **5**, 323-332.

ROBERTS, P. H. & GUBBINS, D. 1987 Origin of the main field: Kinematics. *Geomagnetism*, Ed. *J. A. Jacobs*, Academic Press **2**, 185-249.

ROCHESTER, M. G., JACOBS, J. A., SMYLIE, D. E. & CHONG, K. F. 1975 Can precession power the geomagnetic dynamo? *Geophys. J. R. Astron. Soc.* **43**, 661.

STEWARTSON, K. 1958 On the stability of a spinning top containing liquid. *J. Fluid Mech.* **5**, 577-592.

THOMPSON, R. 1970 Diurnal tides and shear instabilities in a rotating cylinder. *J. Fluid Mech.* **40**, 737-751.

TILGNER, A. 2000 Oscillatory shear layer in source driven flows in an unbounded rotating fluid. *Phys. Fluids* **12**, 1101-1111.

TILGNER, A. 2005 Precession driven dynamos. *Phys. Fluids* **17**, 034104.

VANYO, J. P., WILDE, P. & CARDIN, P. 1995 Experiments on precessing flows in the earth's liquid core. *Geophys. J. Int.* **121**, 136-142.
## 4. INSTABILITÉS DE FLAMBAGE/PLISSEMENT

Une structure constituée d'un film mince (épaisseur $h_F$, module d'Young $E_F$, coefficient de Poisson $\nu_F$) attaché à un substrat plan semi-infini (épaisseur $h_S$, module d'Young $E_S$, coefficient de Poisson $\nu_S$), devient instable lorsqu'elle est soumise à une contrainte (axiale ou bi-axiale) supérieure à un chargement critique.

Lorsque le substrat est beaucoup plus épais que le film, des plis, voir Fig. 4.1, se forment à la surface de ce dernier : on parle alors d'une instabilité (locale) de plissement (wrinkling en Anglais). Le principal objectif des approches analytiques est de caractériser l'effet des paramètres, comme le rapport des épaisseurs $h_F/h_S$ ou le type de chargement (uniaxial, équi-bi-axial, bi-axial), sur la sélection du mode de plissement, sur sa longueur d'onde, sur son amplitude et sur sa stabilité. Généralement, les approches reposent sur la définition d'une énergie potentielle

$$ \mathcal{E} = \mathcal{E}_S + \mathcal{E}_B - W_{ext}, \quad (1) $$

incluant une énergie de déformation de membrane $\mathcal{E}_S$ (stretching energy en Anglais), de flexion $\mathcal{E}_B$ (bending energy) et un terme associé au travail $W_{ext}$ des efforts extérieurs, supposés conservatifs. La minimisation de cette énergie potentielle conduit alors, suivant la cinématique adoptée, à une ou plusieurs équations aux dérivées partielles non-linéaires pour le champ de déplacement et les états d'équilibre statiques associés. Les premiers modèles analytiques pour prédire les instabilités de plissement ont été principalement développés à partir des équations d'équilibre de Von-Karman (équations différentielles partielles non-linéaires, obtenues via la minimisation d'une énergie potentielle construite sous l'hypothèse de déformations planes, où toutes les non-linéarités liées aux composantes du déplacement dans le plan sont négligées). Une analyse linéaire de stabilité par Allen (1969) montre alors que la configuration plane devient instable par plissement pour un chargement critique

$$ \sigma_c^{(W)} = \frac{\overline{E_F}}{4} \left( 3 \frac{\overline{E_S}}{\overline{E_F}} \right)^{2/3}. \quad (2) $$

Dans cette équation, $\overline{E_S}^* = 4E_S/(3 - 4\nu_S)$ et $\overline{E_F} = E_F/(1 - \nu_F^2)$ représentent respectivement les rigidités effectives du substrat et du film. Le mode de l'instabilité est de la forme $w \sim \cos(2\pi/\lambda x)$ (mode cylindrique, voir Fig.
4.1) avec

$$ \lambda = 2\pi h_F \left( \frac{1}{3} \frac{\overline{E_F}}{\overline{E_S}^*} \right)^{1/3}, \qquad (3) $$

sa longueur d'onde. Ce mode cylindrique devient à son tour instable lorsque le chargement de compression dépasse un nouveau seuil critique. Le mode d'instabilité secondaire (mode ondulant ou carré, i.e. herringbone or square-checkerboard mode en Anglais) dépend alors de l'isotropie de la contrainte, voir le diagramme de stabilité de Audoly & Boudaoud (2008a), Fig. 4.1.

Lorsque l'épaisseur du substrat diminue, l'ensemble de la structure peut également fléchir selon un mode global de flambage, voir Fig. 4.2, apparaissant pour un chargement critique $\sigma_c^{(B)}$. Il existe alors une compétition (voire parfois une co-existence) entre le mode de plissement cylindrique et le mode de flambage. Dans le cas d'un film mince attaché à un substrat plan semi-infini, les limites de cette compétition peuvent être établies, voir Wang et al. (2008), en calculant la charge critique de flambage $\sigma_c^{(B)}$ du film et du substrat, modélisés sous la forme d'une poutre composite. La limite $\sigma_c^{(W)} = \sigma_c^{(B)}$ permet alors de définir un diagramme de phase entre l'instabilité de plissement et l'instabilité de flambage, voir Fig. 4.2.

L'extension de cette analyse au cas d'un substrat non-plan<fnref n="1" /> a été réalisée lors de mon séjour postdocoral au M.I.T., et est l'objet de la publication qui suit, parue dans la revue *Journal of the Mechanics and Physics of Solids*, sous le titre "From wrinkling to global buckling of a ring on a curved substrate".

<footnote n="1">¹ Vidéo M.I.T. de l'instabilité de plissement, cas d'un substrat semi-sphérique possédant une cavité intérieure dépressurisée : ici.</footnote>
<figure><img src="image_11.png" /><figcaption>FIGURE 4.1 : (a) Schéma d'un film mince attaché à un substrat élastique plan, sous compression bi-axiale. Modes de l'instabilité de plissement (par Cai *et al.* (2011)) : mode cylindrique (b), ondulant (c), carré (d). (e) Diagramme de phase par Audoly & Boudaoud (2008a) dans le cas général d'une compression bi-axiale non isotrope (compressions directionnelles adimensionnées $η_x$ et $η_y$). La configuration plane est notée (P).</figcaption></figure>
<figure><img src="image_12.png" /></figure>

FIGURE 4.2 : (a) Schéma d'un film mince attaché à un substrat élastique plan, sous compression axiale. (b) Diagramme de phase dans le plan $h_s/h_F$ (rapport des épaisseurs) - $L/h_F$ (rapport d'aspect du film). La courbe indique la transition entre instabilité de flambage et instabilité de plissement, établie théoriquement par Wang et al. (2008).
Abstract

We present a combined analytical approach and numerical study on the sta-
bility of a ring bound to an annular elastic substrate, which contains a circu-
lar cavity. The system is loaded by applying a pressure differential between
the inside of the cavity and the outside of the sample. The ring is modeled
as an Euler-Bernoulli beam and its equilibrium equations are derived from
the mechanical energy which takes into account both stretching and bend-
ing contributions. The curvature of the substrate is considered explicitly to
model the work done by its reaction force on the ring. We distinguish two
different instabilities: periodic wrinkling of the ring or global buckling of
the structure. Our model provides an expression for the critical pressure, as
well as a phase diagram that rationalizes the transition between instability
modes. Towards assessing the role of curvature, we compare our results for
the critical stress and the wrinkling wavelength to their planar counterparts.
We show that the critical stress is insensitive to the curvature of the sub-
strate, while the wavelength is only affected by the discrete values of the
azimuthal wavenumber imposed by the geometry of the problem. Through-
out, we contrast our analytical predictions against finite element simulations.

4.1 Introduction

Wrinkling is a stress-driven mechanical instability that occurs when a stiff
and slender surface layer bonded to a compliant substrate is subject to com-
pression. This universal instability phenomenon spans from numerous natu-
ral and technological/engineering examples. It can be observed over a wide
range of length scales, for example: carbon nanotubes (Lourie et al., 1998),
pre-stretched elastomers used in flexible electronics applications (Kim et al.,
2011), human skin (Chen & Yin, 2010), drying fruit (Yin et al., 2009), surface
morphology of the brain (Budday et al., 2014) and mountain topographies
generated due to tectonic stresses (Price & Cosgrove, 1990; Huddleston &
Lan, 1993).

Over the past decade, there has been an upsurge of interest in the study
of the mechanics of wrinkling, along with a change of paradigm in regarding
surface instabilities as an opportunity for functionality, instead of a first step
in the route to structural failure (Genzer & Groenewold, 2006; Li et al., 2012).
The first mechanical studies of wrinkling were motivated by the stability of
sandwich panels (Allen, 1969), used in lightweight structural applications,
in which the core acts as a soft substrate for the much stiffer skin. More
recently, Bowden et al. (1998) showed how the wrinkling of a thin film on
an elastomeric substrate can be used to produce complex self-organized pat-
terns. Their seminal work has instigated the realization of wrinkling through
several different actuation mechanisms, including thermal mismatch (Huck
et al., 2000), tissue growth/atrophy (Ben Amar & Goriely, 2005; Li et al.,
2011a; Budday et al., 2014), swelling of a liquid (Chan & Crosby, 2006b)
or vapor solvent (Breid & Crosby, 2009), and pneumatics (Terwagne et al.,
2014). The opportunities opened by such a wide range of external stimuli
have enabled the usage of wrinkling in various applications in photonics (Kim
et al., 2012), optics (Chan & Crosby, 2006a), self-assembly (Yoo et al., 2002), microfluidics (Yin et al., 2012) or morphogenesis (Efimenko et al., 2005).

In order to provide a theoretical background to these recent developments, several authors have built on the pioneering work of Allen (1969), who first provided close form solutions for the critical stress and wavelength obtained when an initially straight beam adhered to an infinite plane substrate is placed under a state of uniaxial compression. Chen & Hutchinson (2004) extended this work to consider the case of a plate adhered to a flat substrate under equi-biaxial compression and performed a nonlinear analysis of the Föppl-von Kármán equations (Landau & Lifshitz, 1959; Timoshenko & Gere, 1961). Huang et al. (2005) further refined these efforts by considering the effect of a finite substrate. Both studies showed the existence of multiple buckling modes associated with the same value of critical stress. The stability of these modes under different loadings conditions has been addressed by Audoly & Boudaoud (2008a,b,c), who produced a stability diagram covering the complete evolution from low to high values of overstress. However, experiments by Cai et al. (2011) found disagreement at low values of overstress, suggesting that a finite intrinsic curvature of the system, even if small, may play an important role in dictating the pattern selection.

Early studies of wrinkling on curved substrates, as in the flat configuration, were also motivated by a structural problem, in this case the stability of the outer shell of rockets (Kachman, 1959; Seide & Weingarten, 1961; Seide, 1962). More recent studies that consider instabilities as a possible source of functionality have led to applications of curved configurations in adhesion (Kundu et al., 2011), microfluidics (Mei et al., 2010), morphogenesis of microparticles (Yin et al., 2014), optics (Breid & Crosby, 2013) and aerodynamic drag reduction (Terwagne et al., 2014). Curvature also plays a relevant role in the growth of biological systems (Li et al., 2011b). Despite these important emerging applications, the mechanics of wrinkling on curved substrates remains weakly understood, when compared to the planar counterpart.

Systematic Finite Element simulations have been performed (Yin et al., 2009; Chen & Yin, 2010; Li et al., 2011c; Cao et al., 2012) that highlighted a complex pattern formation process and suggested the possibility of curvature to affect the selected patterns and modify the relevant characteristic length scales, which calls for a robust theoretical backing. Analytical predictions are challenged by the difficulty of modeling the stiffness of the substrate, even in two-dimensional configurations. Cheng (1996) and Cai et al. (2011) used the stiffness provided by Allen (1969) for the flat case, such that their model therefore neglects the contribution of curvature on the response of the substrate. Yin et al. (2009) used the prediction provided by Brush & Almroth (1975), which accounts for curvature but does not consider the influence of the wrinkling wavelength and their prediction does not converge to the classical planar case as the curvature tends to zero. As such, there is a need to quantify the effect of curvature on the stiffness of the substrate and its subsequent influence on wrinkling.

Here, to the best of our knowledge, we provide the first analytical work that accounts for both the curvature of a (2D) shell-substrate system, as
well as the finite size of the substrate. As an initial step, we focus our study on a curved film adhered to a cylindrical substrate, instead of dealing with non-zero Gaussian curvature geometries, which is left for a future study. We assume axial-symmetry to further simplify the system to the 2D problem of a ring on an annular substrate. Mechanical loading is applied by depressurizing a circular cavity inside the substrate, which places the system under a homogenous state of compression. This geometry is motivated by recent experiments (Terwagne et al., 2014) that demonstrated the usage of wrinkling on spherical samples for switchable and tunable aerodynamic drag reduction. In our simplified 2D system, we solve the elasticity problem for the substrate and derive a close form expression for its stiffness, which is then used in the stability analysis of the ring to quantify the buckling patterns.

The paper is organized as follows: in § 4.2 we describe our system along with its material and geometrical parameters. We also describe the possible instability modes, and present a simplified phase diagram, with the aim of providing physical insight on the problem. In § 4.3 we then introduce the kinematics of the ring attached to the substrate and determine the stiffness of the ring. We proceed by defining a strain energy that includes both bending and stretching of the ring, as well as the effect of the substrate. Energy minimization yields the equilibrium equations of the problem. An asymptotic expansion is then used to calculate the principal solution and the bifurcation at the onset of instability. In § 4.4, we describe the finite element simulations that we have performed for this same system.

The results of our investigation are presented in § 4.5. Throughout, we directly compare the analytical predictions to the numerical simulations. We start with the fundamental solution and the critical conditions that lead to instability. We then construct a phase diagram which rationalizes the dependence of the instability modes on the governing parameters. The results for our system are then quantitatively compared to those for wrinkling of a film on a planar substrate, highlighting the effect of curvature. Finally, § 4.6 summarizes our findings and provides perspectives for potential extensions of our work in future studies.

## 4.2 Definition of the problem

We study the stability of a thin elastic ring, bound to an equally curved 2D substrate that contains an inner cavity, a schematic diagram of which is presented in Fig. 4.3(a). The system is loaded by applying a pressure differential, $P$, between the inside of the cavity and the outside of the ring. The thickness of the ring is $H$, its Young's modulus $E_F$ and its Poisson's ratio $\nu_F$. We refer to $\overline{E}_F = E_F/(1-\nu_F^2)$ as the reduced Young's modulus of the film. The substrate is made of a linearly elastic material with Young's modulus $E_S$, Poisson's ratio $\nu_S$ and reduced Young's modulus $\overline{E}_S = E_S/(1-\nu_S^2)$. The thickness of the substrate is $R-R_0$, where $R_0$ is the radius of the inner cavity.

For convenience, we now introduce new rescaled quantities to reduce the number of parameters of the problem. As such, we use $R$ and $\overline{E}_F$ to
<figure><img src="image_13.png" /><figcaption>Figure 4.3: (a) Schematic diagram of our system: a ring is bound to a curved substrate which contains a circular cavity. The system is loaded by applying a pressure differential between the inside of the cavity and the outside of the ring. (b-d) Representative examples of the three possible instability modes of a ring on a curved substrate which contains a cavity that is depressurized. (e) is a zoom in of (d) that exhibits the deformation of the surface of the inner cavity in the Biot mode. The colorbar refers to the principal component of the strain tensor of the mode, which has been normalized.</figcaption></figure>

normalize lengths and pressures and define

$$h = \frac{H}{R}, \quad \beta = \frac{R_0}{R}, \quad \xi = \frac{\overline{E_F}}{\overline{E_S}}, \quad p = \frac{P}{\overline{E_F}}, \qquad (1)$$

as the dimensionless thickness, cavity size, stiffness ratio, and pressure, respectively.

In Fig. 4.3(b-e) we show representative results, obtained from Finite Element Modeling (FEM), of the three possible instability configurations of the ring-substrate system for different values of the dimensionless ring thickness h and ratio of stiffness ξ. All the other mechanical properties are kept constant: the Poisson's ratios are $ν_F = ν_S = 0.5$, and the cavity size is $β = 0.2$. In the simulations this is achieved by taking fixed values for R = 100 units of length and E_S = 1 units of pressure, while changing the values of H, R₀ and E_F accordingly. More details on our numerical simulations are provided in §4.4.

The first mode, a representative example of which is shown in Figure 4.3(b) for $h = 10^{-2}$ and $ξ = 10^3$, corresponds to periodic wrinkling of the film with a well defined wavelength. The strain in the bulk of the substrate is negligible, and no deformation of the inner cavity is observed. The second mode, for example at $h = 10^{-2}$ and $ξ = 10^6$ in Figure 4.3(c), corresponds to a global buckling of the structure, where both the ring and the cavity deform into an ellipse such that the wavelength is $λ = πR$.

In addition to these two instability modes (wrinkling and global buckling), we have also numerically observed an instability on the inner surface
<figure><img src="image_14.png" /><figcaption>Figure 4.4: Schematic phase diagram of the instability modes of the system, obtained for a critical value of the pressure differential: i) wrinkling of the ring or ii) global buckling of the structure (shaded regions). The primary parameters that govern this transition are the stiffness and thickness ratios: $\xi = \overline{E_F}/\overline{E_S}$ and $h = H/R$, respectively. The depicted examples from FEM simulations are for an incompressible film and substrate, $\nu_F = \nu_S = 0.5$. They were obtained for $R = 100$ units of length and $E_S = 1$ units of pressure, while varying $H$ and $E_F$. The colorbar refers to the principal component of the strain tensor of the mode, which has been normalized.</figcaption></figure>
of the cavity. However, this third mode does not affect the ring and is found
for low values of cavity size and stiffness ratio ($\xi = 10^2$ and $h = 10^{-2}$ in
Fig. 4.3c-d). This instability was first discussed by Biot (1965) and we there-
fore refer to it as the Biot mode; it is local in nature and only depends on
the compressive strain at the surface. This type of instability mode has been
recently studied in the case of elastomeric materials with voids by Michel
et al. (2007) and Cai et al. (2010). Understanding the specifics of this Biot
mode is however outside the scope of our work and we shall not take it into
account in our analytical model and systematic numerical investigation.

A schematic phase diagram of our system is provided in Fig. 4.4. For low
values of $h$ and $\xi$ the ring wrinkles with a short wavelength. As either $h$ or
$\xi$ are increased, the wavelength also increases. Once these parameters reach
a critical value, represented by the dashed line in Fig. 4.4, the instability
transitions from wrinkling to global buckling. In what follows, we focus on
rationalizing how the wavelength of the wrinkling mode, and the threshold
value for the transition to global buckling, evolve with the elastic and geo-
metrical parameters of the system.

## 4.3 Analytical approach

In this section we present our analytical approach. The ring is treated as an
Euler-Bernoulli beam. The effect of the substrate is modeled as a restoring
force that acts on the ring and is determined by solving the elasticity problem
of the substrate with adequate boundary conditions. Minimization of the
potential energy provides the equilibrium equations of the problem, which are
solved using an asymptotic expansion that yields the principal and bifurcated
solutions.

4.3.1 Kinematics, energy formulation and equations of equilibrium

We model the ring as an extensible Euler-Bernoulli beam made of an ho-
mogeneous and isotropic material. Polar coordinates are used to track the
position of the ring center-line, *C*. The initial configuration of the ring,
prior to depressurization, is assumed to be circular. The origin, *O*, is lo-
cated at the center of the cavity, and the initial and equilibrium configura-
tions of an arbitrary point of *C* are represented by *M*₀ and *M*, respectively.
Vectors are expressed in the physical base (*e*_r, *e*_θ), derived from the polar
coordinates (*r*, *θ*). The initial position of *C* is **OM**₀ = (*R*, 0), as shown
in the inset of Fig.4.3(a). When the system is loaded by depressurizing
the cavity, *C* deforms to a new configuration given by the position vector
**OM** = *R*(1 + *v*(*θ*), *u*(*θ*)), where *v* and *u* are the dimensionless radial and
orthoradial displacements, respectively.

The infinitesimal arclength of C in the initial and deformed configura-
tions are denoted by ds₀ = |dOM₀| and ds = |dOM|, respectively. More-
over, defining the tangent vector T = dOM/ds, allows us to express the
curvature of C in the deformed configuration as κ/R = |dT/ds|. Here, κ is
dimensionless and can be written in terms of $v$ and $u$ as

$$ \kappa = 1 + (-1 + 2u' + 2v') v'' + v^2 - v - u'^2 + (-u + v') u'' + \frac{1}{2} (v'^2 - u^2) + \text{h.o.t.,} \quad (2) $$

where the prime notation represents derivation with respect to $\theta$ and high order terms (h.o.t.) are neglected under the assumption of small displacements and moderate rotations. We now define the elongation of the ring as $e = ds/ds_0$ to express the stretching deformation, $\eta = (e^2 - 1)/2$, in terms of $v$ and $u$ as

$$ \eta = \frac{1}{2} [u'^2 + v^2 + (u - v')^2] + (1 + v) u' + v, \quad (3) $$

so that the hoop stress in the film is $\sigma_0 = \overline{E_F \eta}$.

Following Euler-Bernoulli beam theory (Timoshenko & Gere, 1961), the total energy of deformation $\mathcal{E}$ of the ring is the sum of a stretching energy $\mathcal{E}_S$ and a bending energy $\mathcal{E}_B$,

$$ \mathcal{E} = \mathcal{E}_S + \mathcal{E}_B = \int_0^{2\pi R} \overline{\mathcal{E}} ds_0, \quad (4) $$

with

$$ \mathcal{E}_S = \int_0^{2\pi R} \frac{\overline{E_F H}}{2} \eta^2 ds_0, \quad (5a) $$

$$ \mathcal{E}_B = \int_0^{2\pi R} \frac{\overline{E_F H^3}}{24R^2} (\kappa - 1)^2 ds_0, \quad (5b) $$

and $\overline{\mathcal{E}}$ the energy of deformation per unit length of the initial configuration of the ring. Assuming that the reaction force of the substrate derives from a potential $\int_0^{2\pi R} \overline{W} ds_0$, the equilibrium states of the ring are the solutions of

$$ \delta \mathcal{E}_S + \delta \mathcal{E}_B - \int_0^{2\pi R} \delta \overline{W} ds_0 = 0, \quad (6) $$

where $\delta A$ is the variation of quantity $A$, for an arbitrary displacement field $R(\delta v, \delta u)$, which is $2\pi$ periodic. The computation of the variations in Eq. (6) leads to the Euler-Lagrange equations for the equilibrium of the ring,

$$ \frac{\partial \overline{\mathcal{E}}}{\partial v} - \left( \frac{\partial \overline{\mathcal{E}}}{\partial v'} \right)' + \left( \frac{\partial \overline{\mathcal{E}}}{\partial v''} \right)'' - \left( \frac{\partial \overline{\mathcal{E}}}{\partial v'''} \right)''' + \left( \frac{\partial \overline{\mathcal{E}}}{\partial v'''} \right)'''' - \frac{\partial \overline{\delta W}}{\partial \delta v} = 0, \quad (7a) $$

$$ \frac{\partial \overline{\mathcal{E}}}{\partial u} - \left( \frac{\partial \overline{\mathcal{E}}}{\partial u'} \right)' + \left( \frac{\partial \overline{\mathcal{E}}}{\partial u''} \right)'' - \left( \frac{\partial \overline{\mathcal{E}}}{\partial u'''} \right)''' + \left( \frac{\partial \overline{\mathcal{E}}}{\partial u'''} \right)'''' - \frac{\partial \overline{\delta W}}{\partial \delta u} = 0, \quad (7b) $$

along with static boundary conditions that are naturally satisfied due to the $2\pi$ periodicity condition on the displacements $v$ and $u$.
### 4.3.2 Asymptotic expansion and reactive force of the substrate

We seek a solution of Eq. (7) as an expansion of the form

$$v = v_0 + \varepsilon A \sin(m\theta) + O(\varepsilon^2), \quad (8a)$$

$$u = \varepsilon B \cos(m\theta) + O(\varepsilon^2), \quad (8b)$$

where $(v_0, 0)$ corresponds to a radial pre-buckling deformation and the term $\varepsilon(A \sin(m\theta), B \cos(m\theta))$ represents an instability of azimuthal wavenumber $m$. From the requirement of $2\pi$ periodic functions $v$ and $u$, $m$ has to be an integer. Also, we consider $m > 1$ since $m = 1$ corresponds to a solid body translation of the ring. The instability displacement field has amplitudes $\varepsilon A$ and $\varepsilon B$, where $\varepsilon$ is a small parameter.

Before substituting Eq. (8) into Eq. (7) and solving at each order in $\varepsilon$, we first need to determine the work, $\delta W$, done by the reaction force, $\mathbf{F} = -(\sigma \mathbf{e}_{\perp} + \tau \mathbf{e}_{\parallel})$, that the substrate exerts on the ring. The normal stress $\sigma$ and the tangential stress $\tau$ at the interface between the ring and the substrate are computed by solving the corresponding two-dimensional elasticity problem using an Airy function, with pressure $P$ at $r = R_0$ and the displacement field given by Eq. (8), at $r = R$. It is worth to note that in our computation we thus enforce the continuity of the displacement field at the interface ring/substrate, contrarily to Seide (1962), who assumed zero shear stress and continuity of normal stress. Both computations can be carried out in a similar way using an Airy function, but they lead to significantly different estimates for the bending and stretching energy in the ring. Audoly & Boudaoud (2008a) have shown that it is preferable to enforce the continuity of the displacement, since the continuity of the stress yields a stretching energy in the ring much larger than the bending energy, in contradiction to what is expected in the wrinkling of a thin film. Our explicit solution for the boundary value problem of the substrate is reported in 4.7.1. In short, we find that the stresses $\sigma$ and $\tau$ at the interface are

$$\frac{\sigma}{E_F} = k_0 v_0 + k \varepsilon A \sin(m\theta) - \gamma p, \quad (9a)$$

$$\frac{\tau}{E_F} = \mu \varepsilon B \sin(m\theta), \quad (9b)$$

with the following governing parameters

$$k_0 = \frac{(1 - \nu_S) (1 - \beta^2) \frac{1}{\xi}}{1 - 2\nu_S + \beta^2}, \quad (10a)$$

$$k = \frac{(1 - \nu_S)}{2} \left( S_A + S_B \frac{B}{A} \right) \frac{1}{\xi}, \quad (10b)$$

$$\mu = \frac{(1 - \nu_S)}{2} \left( S_B \frac{A}{B} + T_B \right) \frac{1}{\xi}, \quad (10c)$$

$$\gamma = \frac{2\beta^2 (1 - \nu_S)}{1 - 2\nu_S + \beta^2}. \quad (10d)$$

Here, $k_0$, $k$, and $\mu$ are the dimensionless pre-wrinkling, wrinkling and shear stiffnesses of the substrate, respectively. The coefficient $\gamma$ quantifies the effect
of pressure, $p$, on the ring through the term $-\gamma p$ in Eq. (9a). The value of $\gamma$ is smaller than one, reflecting the fact that the transmitted pressure decreases within the substrate. There are however two exceptions for which $\gamma = 1$: the limit where there is no substrate, $\beta \to 1$, and the limit of a perfectly incompressible substrate, $\nu_S \to 0.5$, in which the volumetric pressure remains constant throughout the substrate. The quantities $S_A$, $S_B$ and $T_B$ used in Eq. (10) are functions of both the cavity size $\beta$ and the azimuthal wavenumber $m$, and their full expressions are reported in 4.7.1.

Depending on which configuration of the ring is used (initial or deformed), two cases need to be considered to define the normal and tangential vectors $\mathbf{e}_{\perp}$ and $\mathbf{e}_{\parallel}$, respectively. If we define these vectors from the initial configuration, then $(\mathbf{e}_{\perp}, \mathbf{e}_{\parallel}) = (\mathbf{e}_r, \mathbf{e}_{\theta})$, so that the reaction force $\mathbf{F}$ keeps a constant direction while the ring deforms. If we use the deformed configuration, then $\mathbf{F}$ is modeled as a follower force whose direction changes with the ring deformation. In this case, $(\mathbf{e}_{\perp}, \mathbf{e}_{\parallel}) = (\mathbf{N}, \mathbf{T})$, where $\mathbf{N}$ is the inward normal vector orthogonal to the ring center-line, in its deformed state. In order to take both of these options into account in the same model, we define a parameter, $\chi$, which can take the values $\chi = 1$ or $\chi = 0$, when either the undeformed or deformed configurations are used. Thus, we express the reaction force and its elementary work per unit length of the ring as

$$ \mathbf{F} = -[(\sigma\mathbf{e}_r + \tau\mathbf{e}_{\theta})\chi + (\sigma\mathbf{N} + \tau\mathbf{T})(\chi - 1)], \quad (11a) $$

$$ \bar{\delta W} = \mathbf{F} \cdot R (\delta v \mathbf{e}_r + \delta u \mathbf{e}_{\theta}). \quad (11b) $$

where $\cdot$ is the Euclidean dot product.

### The principal solution: order 0 in $\epsilon$

We now proceed to obtain the principal solution and the critical instability modes by substituting Eq. (8) into Eq. (7) and solving the resulting equations for each order of $\epsilon$. At order 0, Eq. (3) for the stretching deformation $\eta$ yields $\eta_0 = v_0 = \sigma_0/\bar{E}_F$ and the linear approximation of the Euler-Lagrange Eq. (7a)

$$ v_0 = \frac{\gamma p}{h + h^3/12 + k_0} = \frac{\sigma_0}{\bar{E}_F}, \quad (12) $$

relates the dimensionless pressure, $p$, and the dimensionless hoop stress, $\sigma_0/\bar{E}_F$, in the ring. Note that Eq. (7b) is automatically satisfied at order 0 in $\epsilon$. In the absence of a substrate, the ratio $\sigma_0/P$ reads

$$ \lim_{\beta \to 1} \left( \frac{\sigma_0}{P} \right) = \frac{1}{h} + O(h), \quad (13) $$

and in the limit of the ring that is much stiffer than the substrate we obtain

$$ \lim_{\xi \to \infty} \left( \frac{\sigma_0}{P} \right) = \frac{\gamma}{h} + O(h), \quad (14) $$

so that $\sigma_0/P < \lim_{\xi \to \infty} (\sigma_0/P) \le \lim_{\beta \to 1} (\sigma_0/P)$. For $\nu_S = 0.5$, we have $\gamma = 1$ through the definition in Eq. (10d) and both limits in Eqs. (13) and (14) are equal. Thus, the hoop stress in a very stiff ring lying on a soft incompressible substrate is the same as the hoop stress in a ring with no substrate, which serves as a verification of the rationale thus far.
The instability: order 1 in $\epsilon$

At order 1 in $\epsilon$, the Euler-Lagrange Eqs. (7) write

$$ (a_1 p + \tilde{a}_1) A + (b_1 p + \tilde{b}_1) B = 0, \quad (15a) $$

$$ (a_2 p + \tilde{a}_2) A + (b_2 p + \tilde{b}_2) B = 0, \quad (15b) $$

where $a_i$, $\tilde{a}_i$, $b_i$ and $\tilde{b}_i$ are functions of the azimuthal wavenumber $m$, the stiffness ratio $\xi$, the cavity size $\beta$ and the substrate Poisson's ratio $\nu_S$, reported in full in 4.7.2. The linear system in Eq. (15) has a nontrivial solution when its determinant vanishes, which for a given value of $m$ occurs at the dimensionless pressure

$$ p_m = \frac{P_m}{E_F} = \frac{-\tilde{a}_1 \tilde{b}_2 + \tilde{b}_1 \tilde{a}_2}{a_1 \tilde{b}_2 + \tilde{a}_1 b_2 - b_1 \tilde{a}_2 - \tilde{b}_1 a_2}. \quad (16) $$

The critical azimuthal wavenumber and the dimensionless critical pressure can now be obtained from Eq. (16) by minimizing over all possible values of $m$

$$ [m_c, p_c] = \min_{m=2,3,...} (p_m). \quad (17) $$

For the general case, this minimization must be performed numerically. However, in the absence of the substrate, we obtain that $m_c = 2$ and Eq. (16) simplifies to

$$ \lim_{\beta \to 1} (p_c) = \lim_{\beta \to 1} (p_2) = -\frac{1}{4-\chi}h^3 - \frac{4}{(4-\chi)^2}h^5 + O(h^7). \quad (18a) $$

The term of order $h^3$ corresponds to the classical dimensionless critical pressure for a ring with no substrate (Lévy, 1884; Carrier, 1947; Timoshenko & Gere, 1961; Combescure, 1981). The term in $h^5$ is a correction also reported by Atanackovic (1998). In the limit of a ring much stiffer than the substrate, we also have $m = 2$, with

$$ \begin{aligned} \lim_{\xi \to \infty} (p_c) &= \lim_{\xi \to \infty} (p_2) \\ &= -\frac{1}{2} \frac{(-1 + 2\nu_S - \beta^2)h^3}{\beta^2(-1 + \nu_S)(4 - \chi)} - \frac{2(-1 + 2\nu_S - \beta^2)h^5}{\beta^2(-1 + \nu_S)(4 - \chi)^2} + O(h^7), \end{aligned} \quad (19) $$

so that $\lim_{\xi \to \infty} (p_c) \le \lim_{\beta \to 1} (p_c) < 0$. Again, for $\nu_S = 0.5$ we observe that the limits equal, corroborating the physical intuition that as the stiffness ratio $\xi$ increases, the effect of the substrate becomes negligible.

### 4.3.3 Critical stress, wavelength and comparison to wrinkling on a planar substrate

To determine the effect of curvature on the instability, we compare the critical hoop stress, $\sigma_c$, in the ring, given by Eq. (12) with $p = p_c$, as well as the
wavelength $\lambda_c = 2\pi R/m_c$ of the wrinkling mode, against their counterparts for an initially planar film on an infinite plane substrate (Chen & Hutchinson, 2004)

$$ \sigma_{Plane} = \frac{\overline{E_F}}{4} \left( 3 \frac{\overline{E_S}}{\overline{E_F}} \right)^{2/3} = \frac{\overline{E_F}}{4} \left( \frac{4(1 - \nu_S)^2}{3 - 4\nu_S} \frac{3}{\xi} \right)^{2/3}, \quad (20a) $$

$$ \lambda_{Plane} = 2\pi H \left( \frac{1}{3} \frac{\overline{E_F}}{\overline{E_S}^*} \right)^{1/3} = 2\pi H \left( \frac{1}{3} \frac{3 - 4\nu_S}{4(1 - \nu_S)^2} \xi \right)^{1/3}, \quad (20b) $$

where $\overline{E_S}^* = 4\overline{E_S}(1 - \nu_S)^2/(3 - 4\nu_S)$ is the effective stiffness of the substrate. Together, Eqs. (20) and (12) yield

$$ \frac{\sigma_c}{\sigma_{Plane}} = \frac{4\gamma p_c}{h + h^3/12 + k_0} \left( \frac{12(1 - \nu_S)^2}{3 - 4\nu_S} \frac{1}{\xi} \right)^{-2/3}, \quad (21a) $$

$$ \frac{\lambda_c}{\lambda_{Plane}} = \frac{1}{hm_c} \left( \frac{3 - 4\nu_S}{12(1 - \nu_S)^2} \xi \right)^{-1/3}, \quad (21b) $$

and, in the limit of absence of the substrate, Eq. (21a) reduces to

$$ \lim_{\beta \to 1} \left( \frac{\sigma_c}{\sigma_{Plane}} \right) = \lim_{\beta \to 1} \left( \frac{\sigma_2}{\sigma_{Plane}} \right) = \frac{4}{4 - \chi} \left( \frac{12(1 - \nu_S)^2}{(3 - 4\nu_S)\xi} \right)^{-2/3} h^2 + O(h^4). \quad (22) $$

In § 4.5.4, we shall make use of Eqs. (21) and Eq. (22) to further quantify the sensitivity of the stress and wavelength to the curvature of the substrate.

## 4.4 Numerical simulations

In § 4.5, we will contrast the predictions from the above analysis with the results of a series of finite element simulations performed in the commercial package Abaqus, using the BUCKLE analysis, which provides the buckling load and the corresponding eigenmodes.

As in the analytical study, the cylindrical structure was modeled as an annulus of a soft substrate with a stiffer thin film adhered to its exterior, under plane strain conditions. The bonding between the substrate and the film is assumed to be perfect, such that both share nodes. Negative pressure is applied on the interior surface, to model the effect of the pressure differential between the inner surface of the annulus (the cavity) and the exterior of the system. The buckling analysis provides the value of the critical pressure, $p_c$, as well as the corresponding critical mode. Rigid body motions are removed by constraining the displacement of two points on the film.

Both substrate and film are modeled as incompressible linearly elastic materials, $\nu_S = \nu_F = 0.5$. These results were compared to additional simulations using a Neo-Hookean model but no difference was observed given the low values of strain involved. The substrate was modeled using quadrilateral plane strain elements. Due to the incompressibility, the corresponding hybrid element, CPE4H, was used. The film was modeled using B21 beam
elements. In order to account for the effect of plane strain, the stiffness of the beam is defined as $E_F / (1 - \nu_F^2)$.

All of the results presented were obtained using 1000 elements in the circumferential direction, and 150 ($R - R_0$) elements in the radial direction. The mesh size was validated with a convergence analysis. By way of example, differences less than 0.5% in critical pressure were obtained when comparing results for a mesh with twice the elements in each direction, even in the case of wrinkling with very short wavelengths. The deviations between the two meshes were, however, larger ($\sim 5\%$) for the case when the critical buckling mode is a Biot instability, due to the infinite number of wavelengths associated to the same buckling mode. However, as stated in § 4.2, this mode shall not be studied in detail as we focus on the wrinkling and global buckling loads.

## 4.5 Results

Having introduced our analytical and numerical methods, we now present the results of a systematic exploration of the mechanical response of our system for different geometric and material parameters. Throughout, we provide a direct comparison between analytical results and numerical simulations, finding good agreement, and instances of discrepancies will be discussed where appropriate.

For the geometric parameters, we have varied the cavity size, $\beta = R_0/R$, and the dimensionless thickness, $h = H/R$. Three representative values were chosen for $\beta$ (depicted in the insets of Fig. 4.5a-c): a small cavity, $\beta = 0.2$; a cavity with size half the external radius, $\beta = 0.5$; and a large cavity, $\beta = 0.8$. Moreover, $h$ was varied in the range $10^{-3}$ to $10^{-1}$. This parameter has two different physical interpretations. On one hand, for a substrate with given curvature, i.e. fixed $R$, increasing $h$ is equivalent to increasing the thickness of the ring. On the other hand, for a ring of given thickness $H$, the value of $h$ decreases with the curvature, $1/R$.

For the material properties, we have considered values for the stiffness ratio between the film and the substrate, $\xi = \bar{E}_F / \bar{E}_S$, spanning over five orders of magnitude, from $10^2$ to $10^7$. The substrate is taken to be incompressible, $\nu_S = 0.5$, since most of the relevant experiments that have motivated our study (Bowden et al., 1998; Huck et al., 2000; Yu & Jiang, 2010; Terwagne et al., 2014) use nearly incompressible elastomeric substrates. Moreover, it is important to note that, even if the analytical model has been derived assuming a general value of $\nu_S$, it is expected to be less accurate for increasing deviations from incompressibility, as discussed previously by Cai et al. (2011). The reason for these discrepancies is that, since the ring is modeled as a beam, it is not possible to enforce both displacement and traction continuity at the interface. In the present work, we impose displacement continuity, as discussed in § 4.3.2.

In the presentation of our results, we first consider the effect of $h$ and $\xi$ on the hoop-stress and on the critical pressure. Then, we rationalize the transition from wrinkling to global buckling, shown in Fig. 4.3. We finally compare the critical stress and wavelength of the wrinkling mode to their
planar substrate counterparts and discuss the effect of curvature.

### 4.5.1 Hoop stress prior to wrinkling

In Fig. 4.5, we plot the hoop stress in the ring, normalized by the pressure, $\sigma_0/P$, as a function of the dimensionless thickness, $h = H/R$. We find that the hoop stress decreases monotonically with the dimensionless thickness, and increases with both the stiffness ratio and the cavity size. When the ring is much stiffer than the substrate (i.e. $\xi \to \infty$), the hoop stress scales as $\sigma_0/P \sim (H/R)^{-1}$, with a prefactor given by Eq. (14). For a given thickness of the ring and a given cavity size, the hoop stress in the ring decreases as the curvature of the substrate increases. This observation is consistent with the classic result for the hoop stress, $\sigma_0 = PR/H$, for a depressurized thin-walled cylindrical pressure vessel. This result can be recovered from Eq. (12) by taking $\gamma = 1$, $k_0 = 0$ and performing a Taylor expansion in $h$, about 0.

### 4.5.2 Critical pressure

In Fig. 4.6, we plot the dimensionless critical pressure, $|p_c| = |P_c|/\overline{E_F}$, as a function of the dimensionless thickness, $h$, and observe two different regimes. For low values of $h$, the ring wrinkles with an azimuthal wavenumber $m_c \gg 2$, with a critical pressure $|p_c|$ that increases with $h$. When $h$ reaches a threshold value, $h^*$, the wavenumber decreases suddenly; the nature of the instability changes from a wrinkling to a global buckling mode with $m_c = 2$. In this regime, there is a clear asymptote as $\xi \to \infty$, given by Eq. (19), which corresponds to the buckling of a ring with no substrate.

Depending on the model used for the reaction force of the substrate (constant direction force, $\chi = 1$, or follower force, $\chi = 0$, see § 4.3.2), Eq. (11) yields two different analytical predictions, shown in Fig. 4.6 as a solid line for $\chi = 0$, and a dashed line for $\chi = 1$. We note that the agreement between FEM simulations and analytical predictions is superior for $\chi = 0$, in particular for the global mode. In other words, it is better to model the reaction force of the substrate as a pressure field which remains normal to the ring center-line $C$ as it deforms, than a force with constant direction. Thus, from now on, all analytical results will be only presented for $\chi = 0$.

Despite the overall good agreement between analytical and numerical predictions, there are noticeable discrepancies in the global modes for low values of the stiffness ratio, for example when $\xi = 10^2$. The reason is the impossibility of our model to enforce both displacement and traction continuity at the interface, as discussed at the beginning of § 4.5. For the wrinkling modes, the shear between ring and substrate is negligible, and consequently inaccuracies in modeling the interface do not affect the model (Audoly & Boudaoud, 2008a). For the global modes, however, the shear stress plays a crucial role, which is underestimated in our model. This leads to discrepancies in the critical pressure with the numerical results, where continuity of both displacement and stress is enforced. The disagreement is more pronounced when the ring and the substrate are of comparable stiffness (e.g. $\xi = 10^2$).
<figure><img src="image_15.png" /><figcaption>Figure 4.5: Pre-instability hoop stress, $\sigma_0$, in the ring, normalized by the pressure, $P$, as a function of the dimensionless ring thickness, $h = H/R$. Cavity sizes are: (a) $\beta = 0.2$, (b) $\beta = 0.5$, (c) $\beta = 0.8$. The Poisson's ratio of the substrate is $\nu_S = 0.5$. Analytical predictions are given by Eq. (12) as solid lines and FEM results are plotted as data points. The legend (bottom right) is common to all three plots.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 4.6: Dimensionless critical pressure $|p_c| = |P_c|/\overline{E_F}$ versus the dimensionless ring thickness $h = H/R$, for several stiffness ratios $\xi = \overline{E_F}/E_S$ and cavity sizes $\beta$: (a) $\beta = 0.2$, (b) $\beta = 0.5$, (c) $\beta = 0.8$. Theory is given by Eq. (17). The Poisson's ratio of the substrate is $\nu_S = 0.5$. The legend (bottom right) is common to all three plots.</figcaption></figure>
Another difference between the analytical predictions and the numerics is the appearance of Biot modes for $\beta = 0.2$ and $\xi = 10^2$. These modes, although possible, only occur in a small region of our parameter space and, as mentioned above, are beyond the scope of this work.

In short, our model exhibits limitations if the ring and the substrate have comparable stiffness and the system globally buckles, or when there is a Biot mode (small cavity). Apart from this extreme combination of parameters, rarely observed in experimental configurations, the model performs successfully.

### 4.5.3 Phase diagram

We proceed by focusing on the transition from wrinkling to global buckling, towards first constructing a phase diagram in the $(\xi, h)$ parameter space and then quantifying the dependence of the boundary, $h^*(\xi)$, between the two modes on the size of the cavity, $\beta$. Given that we do not have a closed form expression for $h^*$, we use a numerical method which tracks any jump from $m_c = 2$ to $m_c^* > 2$ in Eq. (17), when $h$ is decreased.

In Fig. 4.7(a), we plot $h^*$ as a function of $\xi$, for $\beta = 0.5$. The boundary between modes is consistent with a power-law, $h^* \sim \xi^{-3}$, which divides the phase diagram into a wrinkling domain ($h < h^*$) and a global buckling domain ($h > h^*$, shaded region). Again, there is good agreement between numerical and analytical solutions.

In Fig. 4.7(b), we extend the phase diagram to cavity sizes from $\beta = 0.2$ to $\beta = 0.8$. We predict that the power-law with exponent -3, mentioned above for $\beta = 0.5$, is still valid as $\beta$ increases, even if there are some deviations towards the higher values. To quantify the appropriateness of the -3 power-law, we fit a curve of the form $\xi = f_2h^{f_1}$ to the analytically calculated boundaries. In Fig. 4.7(c) we plot the fitting parameters $f_1$ and $f_2$ as a function of $\beta$. We find that the exponent is $f_1 \approx -3$ for $\beta < 0.7$, and then decreases rapidly to a value of $f_1 \approx -4.25$ for $\beta = 0.8$. The decrease of the prefactor $f_2$ is more pronounced and reflects the fact that $h^*$ decreases as $\beta$ increases.

The evolution of $m_c^*$ (i.e. azimuthal wavenumber of the wrinkling mode for $h = h^*$), as a function of $\beta$ and $\xi$ is presented in Fig. 4.7(d). Our prediction shows that $m_c^*$ increases with $\beta$, while it is nearly insensitive to $\xi$. For example, for $\beta = 0.5$, we find that $m_c^* = 9$ for $\xi > 10^2$ and $m_c^* = 10$ for $\xi = 10^2$. Predictions for $\xi > 10^2$ are in relatively good agreement with finite element simulations, which have shown that $m_c^* = 10$. For $\xi = 10^2$, FEM simulations yield $m_c^* = 15$, pointing out again at the limitation of our analytical approach in the limit when the stiffness of the ring and the substrate become comparable.

### 4.5.4 Comparison of wrinkling in our curved system with that on infinite planar substrate

Thus far, we have shown that the ring may wrinkle or buckle globally, depending on the curvature and stiffness of the substrate. We now focus on the wrinkling mode of the ring, with the aim of comparing the critical stress and
<figure><img src="image_17.png" /><figcaption>Figure 4.7: (a) Phase diagram in the $(\xi, h)$ parameter space, showing the transition boundary, $h^*$, from wrinkling to global buckling, for $\beta = 0.5$. (b) Phase diagram for $\beta = 0.2$ to $\beta = 0.8$. Solid lines are analytical predictions and dotted lines correspond to a numerical fit $\xi = f_2 h^{f_1}$. (c) Fitting coefficients, $f_1$ and $f_2$, as functions of $\beta$. (d) Azimuthal wavenumber of the wrinkling mode for $h = h^*$, as a function of $\beta$. The Poisson's ratio of the substrate is $\nu_S = 0.5$.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure 4.8: (a) Critical stress $\sigma_c$ normalized by its planar substrate counterpart, $\sigma_{Plane}$, as a function of the dimensionless ring thickness $h = H/R$. (b) Dimensionless critical normal force versus planar substrate counterpart. Analytical predictions (solid lines) are given by Eq. (21a) and FEM results are shown as data points. Insets are sketches of the wrinkling and global buckling modes. The cavity size is $\beta = 0.5$ and the Poisson's ratio of the substrate is $\nu_S = 0.5$. The legend is common to both plots.</figcaption></figure>

wavelength of our curved system to their counterparts for an infinite planar substrate. We shall center our discussion of this comparison for $\beta = 0.5$. In 4.7.3, we report the results for $\beta = 0.2$ and $\beta = 0.8$, which are qualitatively similar.

In Fig. 4.8(a), we plot $\sigma_c/\sigma_{Plane}$ given by Eq. (21a), as a function of $h$. For $h < h^*$ (wrinkling domain), we find that $\sigma_c/\sigma_{Plane} \approx 1$, in agreement with the FEM simulations. To further quantify how close to unity is this ratio, in Fig. 4.8(b) we plot the dimensionless critical normal force $\sigma_c h / \overline{E_F}$ as a function of the planar result, $\sigma_{Plane} h / \overline{E_F}$. We obtain a line with unit slope, indicating that the substrate curvature has no significative effect on the critical stress for wrinkling. From this observation, and using Eq. (12) with $\sigma_0 = \sigma_{Plane}$, we write the following approximation for the dimensionless critical pressure,

$$|p_c| = \frac{h + k_0 \sigma_{Plane}}{\gamma \overline{E_F}} + O(h^3), \qquad (23)$$

with $k_0$ and $\gamma$ given by Eq. (10). Finally, back to Fig. 4.8(a), for $h > h^*$ (global buckling domain), we find that $\sigma_c/\sigma_{Plane}$ first decreases with $h$, then reaches a local minimum and eventually increases as a power-law with slope 2, towards the asymptotic limit given by Eq. (22). The evolution of $\sigma_c/\sigma_{Plane}$ in the global buckling domain is also reproduced well by FEM simulations.

We now investigate the effect of curvature on the wavelength of the instability mode. In Fig. 4.9(a), we plot $\lambda_c/\lambda_{Plane}$, given by Eq. (21b), as a function of $h$. Focusing on the wrinkling domain, we find that $\lambda_c/\lambda_{Plane} \approx 1$, in agreement with FEM simulations. In Fig. 4.9(b), once again, we quantify how close to unity this ratio is by plotting the critical wavelength $\lambda_c$ as a function of the planar result, $\lambda_{Plane}$. To first approximation, noting the range of five orders of magnitude in both axes of the plot, we find a line with unit slope that passes through the origin, suggesting that the substrate cur-
<figure><img src="image_19.png" /><figcaption>Figure 4.9: (a) Critical wavelength, $λ_c$, normalized by its planar substrate counterpart, $λ_{Plane}$, as a function of the dimensionless ring thickness $h = H/R$. Analytical predictions (solid lines) are given by Eq. (21b). (b) Zoom of the wrinkling domain by plotting $λ_c$ versus $λ_{Plane}$. (c) Azimuthal wavenumber $m_c$ of the wrinkling mode as a function of $h$. Analytical predictions (solid lines) are given by Eq. (17). Dotted lines show $m_c$ for the planar substrate case. Insets on (a) and (b) are sketches of the wrinkling and global buckling modes. Inset on (c) is a zoom in showing the difference in $m_c$ for the curved and planar substrates. The cavity size is $β = 0.5$ and the Poisson's ratio of the substrate is $ν_S = 0.5$. FEM results shown as data points. The legend (bottom right) is common to all three plots.</figcaption></figure>
vature has no significative effect on the wavelength of the wrinkling mode. The only deviations arise from the discrete nature of the wavenumber, since the geometry of the ring enforces $2\pi$ periodic wrinkling modes. To highlight this phenomenon, in Fig. 4.9(c), we plot the critical azimuthal wavenumber $m_c$, given by Eq. (17), as a function of $h$ and we superimpose the planar substrate result (dashed line), $m_{Plane} = 2\pi R/\lambda_{Plane}$. We find that $m_c$ is a decreasing stair function of $h$. The deviation in wavelength between the curved and planar cases scales as $(|\lambda_c - \lambda_{Plane}|)/\lambda_{Plane} \sim 1/m$, which is maximum for $m_c^*$. As shown previously in Fig. 4.7(d), $m_c^*$ increases with $\beta$, hence the deviation $(|\lambda_c - \lambda_{Plane}|)/\lambda_{Plane}$ is maximum for a small cavity size. By way of example, for $\beta = 0.5$, Fig. 4.7(d) indicates that $m_c^* = 9$, leading to a maximum deviation $(|\lambda_c - \lambda_{Plane}|)/\lambda_{Plane} \approx 10\%$, that rapidly decreases as $m$ increases.

## 4.6 Conclusion

We have considered the two-dimensional problem of a ring bound to an elastic substrate which contains a cavity that is depressurized. We have used an energy formulation to derive the Euler-Lagrange equations that govern the equilibrium of the ring, and solved them via an asymptotic expansion. As an improvement to previous results in the literature, our analytical approach accounts for the effect of curvature in modeling the reaction force of the substrate. The results have been compared with numerical simulations.

We first studied the principal solution, obtaining an expression for the hoop stress in the ring as a function of the applied pressure. We then performed a stability analysis of the problem to determine the critical pressure, $P_c$, and the corresponding instability mode. Depending on the dimensionless thickness and stiffness ratio ($h$ and $\xi$) we have identified two different regimes: local wrinkling of the ring, and global buckling of the structure. The boundary between both regions of instability was described via a detailed phase diagram, which quantifies the value of $h$ and $\xi$ at which the transition between instabilities occurs and takes into account the cavity size, $\beta$. Our results can be used as a design guideline to target a desired mode. Finally, we have shown that the critical stress for wrinkling does not depend significantly on the curvature of the substrate. However, curvature imposes a discretization of the wrinkling wavelength due to the periodic closing conditions of the ring.

Our study focused on a 2D curved system which exhibits instability modes analogous to the cylindrical pattern found for uniaxial compression of a film on a flat infinite substrate. Considering more complex loading conditions (e.g. also introducing axial loading) or shells with non-zero Gaussian curvature, should lead to more complex patterns that deserve to be investigated further. Having validated the FEM analysis, as well as carefully considering the elastic response of the curved substrate, extending the study to these other scenarios and addressing the issue of pattern selection is an exciting avenue for future research.
# 4.7 Appendix

## 4.7.1 Response of the substrate

In this appendix we consider the boundary value problem of the substrate, subjected to the ring displacement $R(v(\theta), u(\theta))$, $v$ and $u$ given by Eq. (8), at the interface $r = R$ between the ring and the substrate and to the pressure $P$ at $r = R_0 = \beta R$. The substrate is assumed to be in a plane strain state. We note $\alpha = 3 - 4\nu_S$ and introduce the shear modulus $G = \overline{E_S}(1 - \nu_S)/2$. The components of the stress in the substrate are represented by $\sigma_{rr}$, $\sigma_{r\theta}$, and the displacement field is $(U_r, U_\theta)$. The 2D problem of elasticity is solved by finding an Airy function of the form (Michell, 1899),

$$ \phi_m (r, \theta) = B_1 r^2 + B_2 \ln (r) + \varepsilon (A_1 r^{m+2} + A_2 r^{-m+2} + A_3 r^m + A_4 r^{-m}) \sin (m\theta), \quad (24) $$

where $A_i$ and $B_i$ are unknown constants determined by the boundary conditions

$$ \sigma_{rr} (R_0, \theta) = -P, \quad (25a) $$

$$ \sigma_{r\theta} (R_0, \theta) = 0, \quad (25b) $$

$$ U_r (R, \theta) = Rv (\theta), \quad (25c) $$

$$ U_{\theta} (R, \theta) = Ru (\theta), \quad (25d) $$

which are assumed to apply at $r = R_0$ and $r = R$. The first two equations stand for the continuity of the stress at the boundary of the cavity, whereas the last two stand for the continuity of the displacement at the interface between the substrate and the ring.

The stress and displacement fields resulting from the Airy function Eq. (24) are (Barber, 2002)

$$ \sigma_{rr} (r, \theta) = 2B_1 + \frac{B_2}{r^2} + \begin{pmatrix} -A_1 (m+1) (m-2) r^m \\ -A_2 (m+2) (m-1) r^{-m} \\ -A_3 m (m-1) r^{m-2} \\ -A_4 m (m+1) r^{-m-2} \end{pmatrix} \sin (m\theta), \quad (26) $$

$$ \sigma_{r\theta} (r, \theta) = \begin{pmatrix} -A_1 m (m+1) r^m + A_2 m (m-1) r^{-m} \\ -A_3 m (m-1) r^{m-2} + A_4 m (m+1) r^{-m-2} \end{pmatrix} \cos (m\theta), \quad (27) $$

$$ U_r (r, \theta) = \frac{1}{2\mu} \left[ B_1 (\alpha - 1) r + \begin{pmatrix} A_1 (\alpha - m - 1) r^{m+1} \\ +A_2 (\alpha + m - 1) r^{-m+1} \\ -A_3 m r^{m-1} \\ +A_4 m r^{-m-1} \end{pmatrix} \sin (m\theta) \right], \quad (28) $$

$$ U_{\theta} (r, \theta) = \frac{1}{2\mu} \left[ B_1 (\alpha - 1) r + \begin{pmatrix} -A_1 (\alpha + m + 1) r^{m+1} \\ +A_2 (\alpha - m + 1) r^{-m+1} \\ -A_3 m r^{m-1} \\ -A_4 m r^{-m-1} \end{pmatrix} \cos (m\theta) \right]. \quad (29) $$
Applying the boundary conditions Eq. (25) to Eq. (26) yields a linear system for $A_i$ and $B_i$, with solution

$$B_1 = \frac{1}{2(\beta^2 + 1 - 2\nu_S)} [\overline{E_S v_0} (1 - \nu_S) - P\beta^2], \quad (30a)$$

$$B_2 = \frac{-R^2\beta^2}{\beta^2 + 1 - 2\nu_S} [\overline{E_S v_0} (1 - \nu_S) + P(1 - 2\nu_S)], \quad (30b)$$

$$A_1 = \frac{(m-1)\beta^{2m+2} + (\alpha - m + 1)\beta^{2m} + \beta^2}{2(-m^2 + 1)\beta^{2m+2} + m^2\beta^{2m+4} + \beta^{4m+2}\alpha + (m^2 - 1 + \alpha^2)\beta^{2m} + \alpha\beta^2} \times GR^{-m}A \\ + \frac{(m-1)\beta^{2m+2} + (-m - \alpha + 1)\beta^{2m} - \beta^2}{2(-m^2 + 1)\beta^{2m+2} + m^2\beta^{2m+4} + \beta^{4m+2}\alpha + (m^2 - 1 + \alpha^2)\beta^{2m} + \alpha\beta^2} \times GR^{-m}B, \quad (30c)$$

$$A_2 = \frac{-(m+1)\beta^{2m+2} + \beta^{4m+2} + \beta^{2m}(\alpha + m + 1)}{2(-m^2 + 1)\beta^{2m+2} + m^2\beta^{2m+4} + \beta^{4m+2}\alpha + (m^2 - 1 + \alpha^2)\beta^{2m} + \alpha\beta^2} \times GR^m A \\ + \frac{(m+1)\beta^{2m+2} + \beta^{4m+2} + \beta^{2m}(\alpha - m - 1)}{2(-m^2 + 1)\beta^{2m+2} + m^2\beta^{2m+4} + \beta^{4m+2}\alpha + (m^2 - 1 + \alpha^2)\beta^{2m} + \alpha\beta^2} \times GR^m B, \quad (30d)$$

$$A_3 = \frac{\beta^2 ((\alpha - m + 1)(m + 1)\beta^{2m} + m^2\beta^{2m+2} + \alpha + 1 + m)}{(2(-m^2 + 1)\beta^{2m+2} + m^2\beta^{2m+4} + \beta^{4m+2}\alpha + (m^2 - 1 + \alpha^2)\beta^{2m} + \alpha\beta^2)m} \times GR^{-m+2}A \\ - \frac{\beta^2 (m^2\beta^{2m+2} - 1 + (-m - \alpha + 1)(m + 1)\beta^{2m} - m + \alpha)}{(2(-m^2 + 1)\beta^{2m+2} + m^2\beta^{2m+4} + \beta^{4m+2}\alpha + (m^2 - 1 + \alpha^2)\beta^{2m} + \alpha\beta^2)m} \times GR^{-m+2}B, \quad (30e)$$

$$A_4 = \frac{(-m^2\beta^{2m+2} + (m-1)(\alpha + m + 1)\beta^{2m} - (\alpha - m + 1)\beta^{4m})\beta^2}{(2(-m^2 + 1)\beta^{2m+2} + m^2\beta^{2m+4} + \beta^{4m+2}\alpha + (m^2 - 1 + \alpha^2)\beta^{2m} + \alpha\beta^2)m} \times GR^{m+2}A \\ - \frac{(m^2\beta^{2m+2} + (m-1)(\alpha - m - 1)\beta^{2m} - (-m - \alpha + 1)\beta^{4m})\beta^2}{(2(-m^2 + 1)\beta^{2m+2} + m^2\beta^{2m+4} + \beta^{4m+2}\alpha + (m^2 - 1 + \alpha^2)\beta^{2m} + \alpha\beta^2)m} \times GR^{m+2}B. \quad (30f)$$
Substituting $A_i$ and $B_i$ into Eq. (26) yields the stress at the interface $r = R$

$$ \sigma_{rr}(R, \theta) = K_0 R v_0 + K R \epsilon A \sin(m\theta) - \gamma P, \quad (31a) $$

$$ \sigma_{r\theta}(R, \theta) = M R \epsilon B \cos(m\theta), \quad (31b) $$

where

$$ K_0 = \overline{E_S} \frac{1}{R} \frac{(1 - \nu_S)(1 - \beta^2)}{1 - 2\nu_S + \beta^2}, \quad (32a) $$

$$ K = \overline{E_S} \frac{1}{R} \frac{(1 - \nu_S)}{2} \left( S_A + S_B \frac{B}{A} \right), \quad (32b) $$

$$ M = \overline{E_S} \frac{1}{R} \frac{(1 - \nu_S)}{2} \left( S_B \frac{A}{B} + T_B \right), \quad (32c) $$

$$ \gamma = \frac{2\beta^2(1 - \nu_S)}{1 - 2\nu_S + \beta^2}, \quad (32d) $$

and

$$ \begin{aligned} \psi S_A = & 2(-m^2(\beta^4 + 2 + \alpha) + (3 + \alpha)(m^2 - 1)\beta^2 + 2(\alpha + 1))\beta^{2m} \\ & - (((\alpha + 1)m + \alpha - 1)\beta^{4m} - (\alpha + 1)m + \alpha - 1)\beta^2, \end{aligned} \quad (33a) $$

$$ \begin{aligned} \psi S_B = & -2m(-m^2(\beta^4 + 1) + 2(m^2 - 1)\beta^2 + 1 + \alpha)\beta^{2m} \\ & + (((\alpha - 1)m + \alpha + 1)\beta^{4m} + (\alpha - 1)m - \alpha - 1)\beta^2, \end{aligned} \quad (33b) $$

$$ \begin{aligned} \psi T_B = & 2(-m^2\beta^4 + (1 - \alpha)(m^2 - 1)\beta^2 + \alpha m^2)\beta^{2m} \\ & - (((\alpha + 1)m + \alpha - 1)\beta^{4m} - (\alpha + 1)m + \alpha - 1)\beta^2, \end{aligned} \quad (33c) $$

$$ \psi = (m^2\beta^4 - (m^2 - 1)(2\beta^2 - 1) + \alpha^2)\beta^{2m} + \alpha\beta^2(1 + \beta^{4m}). \quad (33d) $$

The dimensionless stress at the interface and the dimensionless stiffness parameters are obtained from Eq. (31), by dividing it with $\overline{E_F}$,

$$ \overline{\sigma} = k_0 v_0 + k \epsilon A \sin(m\theta) - \gamma p, \quad (34a) $$

$$ \overline{\tau} = \mu \epsilon B \sin(m\theta), \quad (34b) $$

$$ k_0 = \overline{\frac{K_0 R}{E_F}} = \frac{(1 - \nu_S)(1 - \beta^2)}{1 - 2\nu_S + \beta^2} \frac{1}{\xi}, \quad (34c) $$

$$ k = \overline{\frac{KR}{E_F}} = \frac{(1 - \nu_S)}{2} \left( S_A + S_B \frac{B}{A} \right) \frac{1}{\xi}, \quad (34d) $$

$$ \mu = \overline{\frac{MR}{E_F}} = \frac{(1 - \nu_S)}{2} \left( S_B \frac{A}{B} + T_B \right) \frac{1}{\xi}, \quad (34e) $$

as indicated in Eq. (10).

We note that in the case of an inextensible wrinkling mode $A/B = m$, the stiffness $K$ simplifies to $\tilde{K}$ given by

$$ \begin{aligned} \tilde{K} = & \overline{E_S} \frac{1}{R} \frac{(1 - \nu_S)}{2} \left( S_A + S_B \frac{1}{m} \right) \\ = & \overline{E_S} \frac{1}{mR} \frac{2(1 - \nu_S)^2 (\beta^{-2m} - \beta^{2m} + 2m(1 - \beta^{-2})) (m^2 - 1)}{(\beta^{-1} - \beta)^2 (m^2 - 1) + (\beta^{-m} - \beta^m)^2 (3 - 4\nu_S) + \left(\frac{3 - 4\nu_S}{\beta} + \beta\right)^2} \end{aligned} \quad (35) $$
which, for a substrate with no cavity, leads to the limiting case

$$ \lim_{\beta \to 0} (\tilde{K}) = \overline{E_S} \frac{1}{mR} \frac{2(1 - \nu_S)^2 (m^2 - 1)}{3 - 4\nu_S}, \qquad (36) $$

and, for an infinite plane substrate, yields

$$ \lim_{R \to \infty} \left( \lim_{\beta \to 0} (\tilde{K}) \right) = K_{Plane} = \overline{E_S} \frac{4(1 - \nu_S)^2 \pi}{3 - 4\nu_S}, \qquad (37) $$

in agreement with Audoly & Boudaoud (2008a).

### 4.7.2 Terms of the linear analysis of stability

The terms $a_i$, $\tilde{a}_i$, $b_i$ and $\tilde{b}_i$ that appear in Eq. (16) are

$$ a_1 = -4\xi h \left( \left( m^4 + \frac{3}{2} - \frac{11}{4}m^2 \right) h^2 - 9 - 3m^2 \right) \frac{\gamma}{h + h^3/12 + k_0}, \qquad (38a) $$

$$ \tilde{a}_1 = h \left( (m^2 - 1)^2 h^2 + 12 \right) \xi + 6 (1 - \nu_S) S_A, \qquad (38b) $$

$$ b_1 = -4\xi h \left( \frac{1}{4} h^2 m^3 + 12m \right) \frac{\gamma}{h + h^3/12 + k_0}, \qquad (38c) $$

$$ \tilde{b}_1 = -12\xi hm + 6 (1 - \nu_S) S_B, \qquad (38d) $$

$$ a_2 = - \left( h^3 m^2 + 48h + 12k_0 (1 - \chi) \right) m\xi \frac{\gamma}{h + h^3/12 + k_0} + 12\gamma (1 - \chi) m\xi, \qquad (38e) $$

$$ \tilde{a}_2 = -12hm\xi + 6 (1 - \nu_S) S_B, \qquad (38f) $$

$$ b_2 = \left( 36hm^2 + h^3 + 12h + 12k_0 (1 - \chi) \right) \xi \frac{\gamma}{h + h^3/12 + k_0} - 12\gamma (1 - \chi) \xi, \qquad (38g) $$

$$ \tilde{b}_2 = 12hm^2\xi + 6 (1 - \nu_S) T_B \qquad (38h) $$

where $\gamma$, $k_0$, $S_A$, $S_B$ and $T_B$ were given in 4.7.1.

### 4.7.3 Influence of the cavity size

In Figs. 4.10 and 4.11, we plot the $h$ dependence of $\sigma_c/\sigma_{Plane}$ and $\lambda_c/\lambda_{Plane}$, for cavity sizes $\beta = 0.2$ and $\beta = 0.8$, respectively. These plots are qualitatively similar to those obtained for $\beta = 0.5$ in Figs. 4.8 and 4.9, and discussed in 4.5.4 of the main text.
<figure><img src="image_20.png" /><figcaption>Figure 4.10: (a) Critical stress $\sigma_c$ normalized by its planar substrate counterpart $\sigma_{Plane}$, as a function of the dimensionless ring thickness $h = H/R$. Analytical prediction is given by Eq. (21a). (b) Critical wavelength $\lambda_c$ normalized by its planar substrate counterpart $\lambda_{Plane}$, as a function of $h$. Analytical prediction is given by Eq. (21b). Insets are sketches of the wrinkling and global buckling modes. The cavity size is $\beta = 0.2$ and the Poisson's ratio of the substrate is $\nu_S = 0.5$. The legend is common to both plots.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure 4.11: (a) Critical stress $\sigma_c$ normalized by its plane substrate counterpart $\sigma_{Plane}$, as a function of the dimensionless ring thickness $h = H/R$. Analytical prediction is given by Eq. (21a). (b) Critical wavelength $\lambda_c$ normalized by its plane substrate counterpart $\lambda_{Plane}$, as a function of $h$. Analytical prediction is given by Eq. (21b). Insets are sketches of the wrinkling and global buckling modes. The cavity size is $\beta = 0.8$ and the Poisson's ratio of the substrate is $\nu_S = 0.5$. The legend is common to both plots.</figcaption></figure>
Bibliography

ALLEN, H. G. 1969 *Analysis and Design of Structural Sandwich Panels*. Pergamon, New York.

ATANACKOVIC, T. M. 1998 Buckling of a compressible elastic ring. *Acta Mechanica* **127**, 121–134.

AUDOLY, B. & BOUDAOUD, A. 2008a Buckling of a thin film bound to a compliant substrate (part I). Formulation, linear stability of cylindrical patterns, secondary bifurcations. *Journal of the Mechanics and Physics of Solids* **56**, 2401–2421.

AUDOLY, B. & BOUDAOUD, A. 2008b Buckling of a thin film bound to a compliant substrate (part II). A global scenario for the formation of herringbone pattern. *Journal of the Mechanics and Physics of Solids* **56**, 2422–2443.

AUDOLY, B. & BOUDAOUD, A. 2008c Buckling of a thin film bound to a compliant substrate (part III). Herringbone solutions at large buckling parameter. *Journal of the Mechanics and Physics of Solids* **56**, 2444–2458.

BARBER, J. R. 2002 *Elasticity: 2nd Edition*. Kluwer Academic Publishers.

BEN AMAR, M. & GORIELY, A. 2005 Growth and instability in elastic tissues. *Journal of the Mechanics and Physics of Solids* **53**, 2284–2319.

BIOT, M. 1965 *Mechanics of Incremental Deformation*. Wiley, New York.

BOWDEN, N., BRITTAIN, S., EVANS, A. G., HUTCHINSON, J. W. & WHITESIDES, G. M. 1998 Spontaneous formation of ordered structures in thin films of metals supported on an elastomeric polymer. *Nature* **393**, 146–149.

BREID, D. & CROSBY, A. J. 2009 Surface wrinkling behavior of finite circular plates. *Soft Matter* **5**, 425–431.

BREID, D. & CROSBY, A. J. 2013 Curvature-controlled wrinkle morphologies. *Soft Matter* **9**, 3624–3630.

BRUSH, D. O. & ALMROTH, B. O. 1975 *Buckling of Bars, Plates and Shells*. McGraw–Hill, New York.

BUDDAY, S., STEINMANN, P. & KUHL, E. 2014 The role of mechanics during brain development. *Journal of the Mechanics and Physics of Solids* **72**, 75–92.

CAI, S., BERTOLDI, K., WANG, H. & SUO, Z. 2010 Osmotic collapse of a void in an elastomer: breathing, buckling and creasing. *Soft Matter* **6**, 5770–5777.

CAI, S., BREID, D., CROSBY, A. J., SUO, Z. & HUTCHINSON, J. W. 2011 Periodic patterns and energy states of buckled films on compliant substrates. *Journal of the Mechanics and Physics of Solids* **59**, 1094–1114.