CAO, Y.-P., LI, B. & FENG, X.-Q. 2012 Surface wrinkling and folding of core-shell soft cylinders. *Soft Matter* **8**, 556.

CARRIER, G. F. 1947 On the buckling of elastic rings. *Journal of Mathematics and Physics* **26**, 94–103.

CHAN, E. P. & CROSBY, A. J. 2006a Fabricating microlens arrays by surface wrinkling. *Advanced Materials* **18**(24), 3238–3242.

CHAN, E. P. & CROSBY, A. J. 2006b Spontaneous formation of stable aligned wrinkling patterns. *Soft Matter* **2**(4), 324–328.

CHEN, X. & HUTCHINSON, J. W. 2004 Herringbone buckling patterns of compressed thin films on compliant substrates. *Journal of Applied Mechanics* **71**(5), 597–603.

CHEN, X. & YIN, J. 2010 Buckling patterns of thin films on curved compliant substrates with applications to morphogenesis and three-dimensional micro-fabrication. *Soft Matter* **22**(6), 5667–5680.

CHENG, P. 1996 Weight optimization of cylindrical shells with cellular cores. Master of Science manuscript.

COMBESCURE, A. 1981 Calcul analytique de la déformation d'un anneau avec défauts. DEMT/SMTS/BAMS/81-142 F.E. : 5423-05-000-09.

EFIMENKO, K., RACKAITIS, M., MANIAS, E., VAZIRI, A., MAHADEVAN, L. & GENZER, J. 2005 Nested self-similar wrinkling patterns in skins. *Nature Materials* **4**, 293–297.

GENZER, J. & GROENEWOLD, J. 2006 Soft matter with hard skin: From skin wrinkles to templating and material characterization. *Soft Matter* **2**, 310–323.

HUANG, Z., HONG, W. & SUO, Z. 2005 Nonlinear analyses of wrinkles in a film bonded to a compliant substrate. *Journal of the Mechanics and Physics of Solids* **53**, 2101–2118.

HUCK, W. T. S., BOWDEN, N., ONCK, P., PARDOEN, T., HUTCHINSON, J. W. & WHITESIDES, G. M. 2000 Ordering of spontaneously formed buckles on planar surfaces. *Langmuir* **16**, 3497–3501.

HUDDLESTON, P. J. & LAN, L. 1993 Information from fold shapes. *Journal of Structural Geology* **15**, 253–264.

KACHMAN, D. 1959 Test report on buckling of propellant cylinders under compressive loads. Space Technology Labs, Inc., GM-59-7520.6-24.

KIM, D. H., LU, N., MA, R., KIM, Y. S., KIM, R. H., SHUODAO, W., WU, J., WON, S. M., KIM, T. I., CHOWDHURY, R., YING, M., XU, L., LI, M., CHUNG, H. J., KEUM, H., MCCORMICK, M., LIU, P., ZHANG, Y. W., OMENETTO, F. G., HUANG, Y., COLEMAN, T. & ROGERS, J. A. 2011 Epidermal electronics. *Science* **333**, 838.
KIM, J. B., KIM, P., PÉGARD, N. C., OH, S. J., KAGAN, C. R., FLEISCHER, J. W., STONE, H. A. & LOO, Y. L. 2012 Wrinkles and deep folds as photonic structures in photovoltaics. *Nature Photonics* **6**, 327–332.

KUNDU, S., DAVIS, C., LONG, T., SHARMA, R. & CROSBY, A. 2011 Adhesion of nonplanar wrinkled surfaces. *Journal of Polymer Science Part B: Polymer Physics* **49**, 179–185.

LANDAU, L. D. & LIFSHITZ, E. M. 1959 *Theory of Elasticity*. Pergamon, London.

LÉVY, M. M. 1884 Mémoire sur un nouveau cas intégrable du problème de l'élastique et de l'une de ses applications. *Journal de Mathématiques Pures et Appliquées* **3X**, 5–42.

LI, B., CAO, Y. P. & FENG, X. Q. 2011a Growth and surface folding of esophageal mucosa: a biomechanical model. *Journal of Biomechanics* **44**, 182–188.

LI, B., CAO, Y.-P., FENG, X.-Q. & GAO, H. 2011b Surface wrinkling of mucosa induced by volumetric growth : Theory, simulation and experiment. *Journal of the Mechanics and Physics of Solids* **59**, 758–774.

LI, B., CAO, Y.-P., FENG, X.-Q. & GAO, H. 2012 Mechanics of morphological instabilities and surface wrinkling in soft materials: a review. *Soft Matter* **8**, 5728.

LI, B., F., J., CAO, Y. P., FENG, X. Q. & GAO, H. 2011c Surface wrinkling patterns on a core–shell soft sphere. *Physical Review Letters* **106**, 234301.

LOURIE, O., COX, D. M. & WAGNER, H. D. 1998 Buckling and collapse of embedded carbon nanotubes. *Physical Review Letters* **81**, 1638–1641.

MEI, Y., KIRAVITTAYA, S., HARAZIM, S. & SCHMIDT, O. G. 2010 Principles and applications of micro and nanoscale wrinkles. *Materials Science and Engineering R* **70**, 209–224.

MICHEL, J., LOPEZ-PAMIES, O., PONTE CASTAÑEDA, P. & TRIANTAFYL-LIDIS, N. 2007 Microscopic and macroscopic instabilities in finitely strained porous elastomers. *Journal of the Mechanics and Physics of Solids* **55**, 900–938.

MICHELL, J. H. 1899 On the direct determination of stress in an elastic solid, with application to the theory of plates. *Proceedings of the London Mathematical Society* **s1-31(1)**, 100–124.

PRICE, N. J. & COSGROVE, J. W. 1990 *Analysis of Geological Structures*. Cambridge University Press, Cambridge.

SEIDE, P. 1962 The stability under axial compression and lateral pressure of circular-cylindrical shells with a soft elastic core. *Journal of the Aerospace Sciences* **29(7)**, 851–862.
SEIDE, P. & WEINGARTEN, V. I. 1961 Buckling of circular rings and long
cylinders enclosing an elastic material under uniform external pressure.
*American Rocket Society Journal* **32**(5), 680–687.

TERWAGNE, D., BROJAN, M. & REIS, P. 2014 Smart morphable surfaces for aerodynamic drag control. *Advanced Materials* p. In press.

TIMOSHENKO, S. P. & GERE, J. M. 1961 *Theory of Elastic Stability*.
McGraw-Hill, New York.

WANG, S., SONG, J., KIM, D., HUANG, Y. & ROGERS, J. 2008 Local ver-
sus global buckling of thin films on elastomeric substrates. *Applied Physics
Letters* **93**(2), 023126.

YIN, J., CHEN, X. & SHEINMAN, I. 2009 Anisotropic buckling patterns in
spheroidal film/substrate systems and their implications in some natural
and biological systems. *Journal of the Mechanics and Physics of Solids*
**57**, 1470–1484.

YIN, J., HAN, X., CAO, Y. & LU, C. 2014 Surface wrinkling on poly-
dimethylsiloxane microspheres via wet surface chemical oxidation. *Scientific Reports* **4**, 5710.

YIN, J., YAGÜE, J. L., EGGENSPIELER, D., GLEASON, K. K. & BOYCE, M. C. 2012 Deterministic order in surface micro-topologies through sequential wrinkling. *Advanced materials* **24**, 5441–5446.

YOO, P. J., SUH, K. Y., PARK, S. Y. & LEE, H. H. 2002 Physical self-
assembly of microstructures by anisotropic buckling. *Advanced materials*
**14**, 1383–1387.

YU, C. & JIANG, H. 2010 Forming wrinkled stiff films on polymeric sub-
strates at room temperature for strechable interconnects applications.
*Thin Solid Films* **519**, 818–822.
## 5. INSTABILITÉ FLUIDE-ÉLASTIQUE D'UN FAISCEAU DE TUBES

Les phénomènes d'interaction fluide-structure (I.F.S.) sont omniprésents dans les installations industrielles où des structures solides sont en contact avec un fluide qui exerce un chargement mécanique. Pour des structures élancées et souples, l'I.F.S. peut induire des phénomènes vibratoires et des instabilités mécaniques à l'origine de grandes amplitudes de déplacements. L'industrie nucléaire est confrontée à cette problématique, notamment au niveau des tuyauteries, des assemblages de combustible, des commandes de grappe ou encore des générateurs de vapeur (G.V.). Un réacteur à eau préssurisée comprend habituellement trois ou quatre générateurs de vapeur, permettant le transfert de chaleur du circuit primaire vers le circuit secondaire, les deux circuits restant isolés. L'échange thermique entre les deux circuits se fait au travers de tubes en forme de U renversé. L'eau chaude du circuit primaire, pressurisée (155 bar) et chauffée (320 °C) par le cœur du réacteur nucléaire s'écoule à l'intérieur des tubes alors que l'eau froide du circuit secondaire s'écoule à l'extérieur des tubes. Cette eau, portée à ébullition, se transforme en vapeur qui alimente les turbines génératrices d'électricité. Un G.V. est constitué de plusieurs milliers de tubes (faisceau de tubes) représentant une surface d'échange de l'ordre de 5000 m². Ces tubes élancés sont encastrés en partie basse sur une plaque tubulaire et guidés par des plaques entretoises composées d'orifices qui permettent à l'eau liquide du circuit secondaire de s'écouler et de se transformer progressivement en vapeur. En tant que structures élancées et flexibles, les tubes de G.V. sont sensibles aux mécanismes de l'I.F.S., en particulier au niveau de la partie supérieure cintrée (chignon) soumise à un écoulement transverse diphasique eau-vapeur. Les vibrations liées à cette interaction induisent des phénomènes (fatigue et usure par frottement et chocs) pouvant entrainer de graves conséquences pour la sûreté nucléaire. Pour faire face aux différents modes de dégradation des tubes de G.V. et éviter des fuites entre le circuit primaire radioactif et le circuit secondaire, des contrôles avec des sondes passant dans les différents tubes sont réalisés lors des arrêts des réacteurs. En fonction des résultats de ces contrôles, des bouchages de tubes via l'obturation de leur entrées et sorties sont réalisés de manière préventive, voir Fig. 5.1. Enfin, lorsque les dégradations des tubes deviennent trop importantes, des opérations lourdes de remplacement de l'ensemble des générateurs de vapeur associés à un réacteur sont réalisées. Dans le cadre de la prévention des risques vibratoires des tubes de G.V., et afin d'améliorer le dimensionnement et évaluer les marges en fonctionnement normal ou dégradé des G.V., des études expérimentales, numériques et théoriques sont ainsi conduites par les acteurs du nucléaire. Pour les tubes
de G.V., ces études portent essentiellement sur les mécanismes de l'I.F.S. en écoulement transverse monophasique (I.F.S. en partie basse des tubes) ou diphasique (I.F.S. au niveau du chignon). L'étude des mécanismes d'I.F.S. dans ces deux zones est alors réalisée à l'aide de faisceaux analytiques de tubes droits soumis à un écoulement transverse, monophasique ou diphasique.

Dans le but de mieux appréhender et quantifier les phénomènes vibra-
toires des tubes de G.V., des travaux expérimentaux sont menés au sein du
laboratoire d'Etudes de dynamique (DYN) du C.E.A., sur des maquettes à
échelle réduite. Afin de limiter le coût des installations et dans un souci de
sécurité, les études sont réalisées avec des fluides simulants, en l'occurrence
de l'eau et de l'air. Les deux principales maquettes exploitées au sein du
laboratoire DYN s'appellent DIVA et TITAN, voir Fig. 5.2. La maquette
DIVA est principalement dédiée à l'étude des vibrations dans les faisceaux à
pas carré. Elle reproduit ainsi à échelle réduite l'I.F.S. dans les générateurs
de vapeur qui équipent actuellement les 32 réacteurs nucléaires de 900 MWe
et les 20 réacteurs de 1300 MWe du parc Français. La maquette TITAN est,
quant à elle, dédiée à l'étude des vibrations dans les faisceaux à pas trian-
gulaire. Ces faisceaux permettent de placer un plus grand nombre de tubes
au sein des G.V. qui équipent actuellement les 4 réacteurs de 1450 MWe de
palier N4 du parc Français et équiperont les futurs EPR. Les deux maquettes
sont alimentées par une boucle équipée d'une pompe pouvant débiter de 5 à
300 l/s d'eau et d'un compresseur à air pouvant débiter jusqu'à 3000 m³/h
sous une pression absolue de 8 bar. La boucle peut ainsi fonctionner en écou-
lement d'eau, en écoulement d'air ou en écoulement diphasique eau-air. La
maquette DIVA est constituée de 15 tubes en acier inox placés sur 5 rangées
et de 10 demi-tubes pour limiter les effets de bord. Les tubes ont un diamètre
D = 30 mm, une longueur L = 300 mm et le faisceau un pas carré (distance
entre deux tubes) P = 45 mm. La maquette TITAN est, quant à elle, com-
posée de tubes similaires à ceux utilisés pour DIVA, mais disposés selon un
arrangement triangulaire parallèle et un pas P = 43.2 mm. Ces deux ma-
quettes permettent d'étudier plusieurs configurations flexibles, notamment
le cas d'un tube unique monté sur un bi-lame flexible et positionné au centre
d'un faisceau rigide, voir Fig. 5.2 (a). Deux caractéristiques fondamentales
de ces maquettes sont i) la possibilité de découpler le mouvement des tubes
selon les directions de portance (direction perpendiculaire à l'écoulement) et
de trainée (même direction que l'écoulement), ii) l'utilisation d'un système
de contrôle actif sans contact (actionneur électro-dynamique), voir Fig. 5.2
(b), permettant d'amortir les vibrations des tubes afin d'explorer leur dyna-
mique sous d'importantes vitesses d'écoulement. La mesure expérimentale
de la force fluide s'exerçant sur un tube flexible est réalisée par la méthode
directe ou indirecte. Dans la méthode directe, le déplacement du tube est
imposé par un pot vibrant (actionneur PRODERA), voir Fig. 5.2 (b), et la
force fluide est directement mesurée par un capteur piézoélectrique. Dans la
méthode indirecte, le mouvement du tube est libre, piloté par son interac-
tion naturelle avec l'écoulement incident. La force fluide est alors obtenue
indirectement à partir de la mesure, via un capteur laser, du déplacement et
de la vitesse du tube. Les variations de fréquence et d'amortissement de ces
<figure><img src="image_22.png" /><figcaption>FIGURE 5.1 : (a) Schéma d'un générateur de vapeur, source IRSN. (b) Photo d'un assemblage de tubes de G.V. avec zoom sur le chignon. Vidéo Areva sur les vibrations des tubes de G.V. : ici. (c) Evolution du nombre de tubes bouchés en fonction des années, source EDF.</figcaption></figure>
grandeurs par rapport aux références mesurées en fluide au repos permettent alors de remonter aux caractéristiques de la force fluide en écoulement.

En parallèle de son activité expérimentale, le laboratoire DYN conduit des études semi-analytiques pour déterminer le chargement mécanique exercé par un écoulement sur un tube flexible. La modélisation de ce chargement repose essentiellement sur l'hypothèse d'un découplage de la force fluide en deux phénomènes distincts

* les phénomènes aléatoires, modélisés par une force<fnref n="1" /> fluide indépendante du mouvement du tube. Cette force, aussi appelée excitation aléatoire, est liée à l'intensité turbulente générée par l'écoulement. L'augmentation de cette intensité avec la vitesse de l'écoulement conduit à une croissance de l'amplitude vibratoire de la structure, voir Fig. 5.3.

* les phénomènes couplés, modélisés par une force fluide-élastique $F_{FE}$ qui dépend du déplacement $X$, de la vitesse $\frac{dX}{dt}$ et de l'accélération $\frac{d^2X}{dt^2}$ du tube. Cette force est traditionnellement décomposée sous la forme

$$F_{FE} = -M_f \frac{d^2 X}{dt^2} - C_f \frac{dX}{dt} - K_f X, \quad (1)$$

où $(M_f, C_f, K_f)$ représentent des coefficients de masse, d'amortissement et de rigidité ajoutés.

Sous l'hypothèse de décomposition de (1), le système couplé fluide-structure est alors caractérisé par un triplé de coefficients $M_s + M_f, C_s + C_f, K_s + K_f$, où $M_s, C_s, K_s$ représentent la masse, l'amortissement et la rigidité du tube en air. Pour une vitesse d'écoulement critique $V_c$, l'amortissement total $C_T = C_s + C_f$ du système couplé devient nul, entrainant une instabilité fluide-élastique qui se manifeste sous la forme d'une croissance exponentielle de l'amplitude de vibration, voir Fig. 5.3. Ce type de phénomène peut conduire à une usure et un percement rapide des tubes de G.V. et a donc conditionné la plupart des études menées sur les faisceaux. Un adimensionnement est alors nécessaire pour transposer les résultats de ces études, réalisées sur des maquettes à échelle réduite, vers le cas du G.V. à échelle industrielle.

Les principaux nombres sans dimension permettant cette transposition d'échelle sont

* le rapport d'aspect, $l = L/D$, qui caractérise l'élancement des tubes sous écoulement,

* le pas réduit, $p = P/D$, qui caractérise la taille de l'espace libre utilisé par le fluide pour circuler entre les tubes,

* le nombre de Reynolds, $Re = DV/\nu$, avec $V$ la vitesse de l'écoulement au sein du faisceau de tubes et $\nu$ la viscosité cinématique du fluide de densité $\rho$. Ce nombre caractérise principalement la vitesse de l'écoulement,

<footnote n="1">¹ Une collaboration (thèse de C. Bazin, 2020-2023) entre le laboratoire DYN (R. Lagrange, encadrant C.E.A.), le Laboratoire de Modélisation et Simulation en mécanique des Fluides (M. G. Rodio, encadrant C.E.A.) et le LadHyX (C. Josserand, directeur) est en cours pour reproduire numériquement cette force observée dans la maquette expérimentale DIVA.</footnote>
<figure><img src="image_23.png" /><figcaption>FIGURE 5.2 : (a) Schéma de la maquette DIVA (faisceau à pas carré). (b) Photo de l'instrumentation. Vidéo C.E.A. de l'instabilité fluide-élastique en portance, configuration DIVA, écoulement monophasique : ici. (c) Photo de la maquette TITAN (faisceau à pas triangulaire).</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>FIGURE 5.3 : Evolution de l'amplitude vibratoire en fonction du débit de l'écoulement. L'excitation aléatoire est liée à l'intensité turbulente de l'écoulement. Le détachement tourbillonnaire, absent en écoulement diphasique, correspond à une organisation périodique de l'écoulement autour du tube. L'instabilité fluide élastique est caractérisée par une croissance exponentielle de l'amplitude vibratoire. Source : S. Caillaud.</figcaption></figure>

* le nombre de Stokes, $Sk = F_0 D^2 / \nu$, avec $F_0$ la fréquence du tube en fluide au repos. Ce nombre caractérise principalement la flexibilité du tube,

* la vitesse réduite, $Vr = V/F_0 D$, qui caractérise le rapport entre la vitesse de l'écoulement et une vitesse caractéristique de vibration du tube. Ce nombre sans dimension s'exprime en fonction de $Re$ et $Sk$ sous la forme $Vr = Re/Sk$,

* le nombre de Scruton, $Sc = C_0 / (2\rho D^2 L F_0)$ avec $C_0$ l'amortissement du tube en fluide au repos. Ce nombre caractérise principalement la masse et l'amortissement du tube et est parfois appelé paramètre masse-amortissement.

L'ensemble des nombres sans dimension caractérise alors le système couplé fluide-structure, en particulier l'évolution de son amortissement total $C_T = C_s + C_f$ et de son seuil de stabilité $C_T = 0$.

La détermination du seuil de stabilité $C_T = 0$ d'un tube flexible positionné au centre d'un faisceau rigide, à pas carré, soumis à un écoulement transverse monophasique et diphasique (eau-air), est l'objet de la publication qui suit, soumise dans la revue *Journal of Fluids and Structures*, sous le titre "A new criterion for the instability threshold of a square tube bundle subject to an air-water cross-flow".
Abstract

In this paper, we investigate the fluid-elastic instability of a square tube bundle subject to two-phase cross-flow. A dimensional analysis is carried out, leading to a new criterion of instability. This criterion establishes a direct link with the instability thresholds in single-phase flows. In parallel to the dimensional analysis, experimental work is carried out to i) determine the instability thresholds in single-phase flows (new relation between the Scruton, Stokes and Reynolds number), ii) to test the validity of the two-phase flow instability criterion, derived from the dimensional analysis. The experiments are carried out on a square tube bundle (pitch ratio of 1.5) consisting of 5 rows of 3 tubes (plus two end-rows of half tubes). The central tube is mounted on two flexible blades allowing a vibration in the transverse direction only, whereas all the other tubes are rigid. The instability threshold in single-phase (water) flow is obtained from a method of direct measurement of the fluid-elastic forces, in which the motion of the central tube is imposed. The instability threshold in two-phase (air-water) flow is obtained from a method of indirect measurement, with an active system stability control, of the fluid-elastic forces, in which the central tube vibrates freely. Three sets of blades with different stiffnesses are tested to investigate the stability of the central tube. The criterion of instability is in very good agreement with air-water experiments, predictive for all homogeneous void fractions, and so for all flow regimes (identified in our experiments with an optical probe). This new criterion is of theoretical interest for the understanding of complex two-phase flow excitations, as well as of practical significance for the predictive analysis of industrial components.

**5.1 Introduction**

The knowledge of the fluid force acting on a structure subject to a cross flow is a crucial information that must be accounted for when designing heat-exchanger tube bundles. The large vibrations resulting from a fluid-elastic instability may lead to some mechanical degradation of the concerned tube, which may affect the power plant operation and safety. This instability can be described as a self-excited feedback mechanism between the motion of the structure and the fluid forces.

The first work on fluid-elastic instability was published by Roberts (1962), who highlighted self-sustained oscillations within a row of cylinders subject to cross flow; the cylinders being alternately rigid and flexible in the direction of the flow. Connors (1970) studied the interaction between a row of flexible cylinders and a transverse flow and developed a quasi-static model following the work of Roberts (1962). He proposed an estimation of the critical instability velocity for a tube as a function of a mass damping parameter and a constant called the Connors constant. Several other theoretical models were subsequently proposed for the study of the fluid-elastic instability: the inviscid or velocity potential model by Dalton & Helfinstine (1971), the quasi-steady model by Price & Paidoussis (1982, 1983), the unsteady model by Tanaka & Takahara (1980, 1981); Chen (1983a,b, 1987), the
quasi-unsteady model by Granger & Paidoussis (1996) or the flow channel model by Lever & Weaver (1986).

Following the work of Roberts (1962) and Connors (1970), a wealth of ex-
perimental studies were also carried out, primarily to refine Connors' model
and determine the value of the Connors constant for various types of tube
bundles. Initially, most of the experiments were performed in single-phase
flows (liquid or gaseous) and it was not until the 1980s that studies on tube
bundles subject to two-phase flow became widespread, see for e.g. Axisa
et al. (1984, 1985); Nakamura et al. (1986a,b). Most experimental works
were done in air-water flows due to the very high cost of steam-water loops.
This however raised the question of the validity of experimental results for the
tubes of a steam generator that actually operate in a steam-water mixture.
To overcome this problem, some researchers as Pettigrew et al. (1995, 2002);
Pettigrew & Taylor (2009) used a mixture of liquid-vapor freon due to the
proximity of its density ratio to that of steam-water. The results have shown
very little difference in terms of critical velocity whether a mixture of freon,
an air-water or a steam-water flow is considered. It is therefore believed that
the results in air-water are representative of what happens in a steam gener-
ator. Still, it should be mentioned that, apart from experimental Axisa et al.
(1984, 1985); Nakamura et al. (1986a,b) and numerical Sadek et al. (2018,
2020) studies, most theoretical models have been developed for single-phase
flows, whereas steam generators operate mainly in two-phase flows. Re-
cently, Yasser et al. (2021) relied on the inviscid flow channel approach of
Hassan & Weaver (2015, 2016) to develop a model to predict fluid-elastic
instability in a square tube bundle subject to a two-phase bubbly cross-flow.
In his work, Yasser et al. (2021) treated the liquid phase as a continuous
medium and modeled gas bubbles as compressible dispersed particles. The
continuity and momentum equations are solved numerically to calculate the
fluid forces on the tube and the flow properties around it. The prediction
of the stability threshold obtained from this model showed promising results
compared to the experimental data. In our work, we present a very different
approach, in which the two-phase instability threshold is not obtained from
the implementation of a theoretical model in a numerical simulation code.
This work, an augmented version of Lagrange et al. (2022), is based solely
on dimensional analysis and on a hypothesis of similarity with single-phase
flow.

5.2 Position of the problem and dimensional analysis

We consider the lift vibration of a flexible tube (diameter *D*, length *L*), part of a square tube bundle (pitch *P*), subject to a two-phase cross-flow characterized by the mass densities (*ρ*₁, *ρ*₂), the kinematic viscosities (*ν*₁, *ν*₂), the superficial gap velocities (*J*₁, *J*₂) of its two constitutive phases, as well as a surface tension *σ*, see Fig. 5.4. The homogeneous void fraction is defined as *α* = *J*₂/(*J*₁ + *J*₂). We note *M*<sub>s</sub>, *C*<sub>s</sub> and *K*<sub>s</sub> the mass, damping and rigidity coefficients of the first mode of vibration of the flexible tube in vacuum. The modal displacement *X* in the lift direction is assumed to satisfy the equation
<figure><img src="image_25.png" /><figcaption>Figure 5.4: Sketch of a square tube bundle subject to two-phase cross-flow with superficial gap velocities $J_1$ and $J_2$. The two phases have mass densities $(\rho_1, \rho_2)$ and kinematic viscosities $(\nu_1, \nu_2)$. The central cylinder vibrates in the lift direction.</figcaption></figure>

of motion

$$ (M_s + M_f) \frac{d^2 X}{dT^2} + (C_s + C_f) \frac{dX}{dT} + (K_s + K_f) X = F_t, \quad (1) $$

with $T$ the time and where $M_f, C_f$ and $K_f$ are the fluid added mass, damping and rigidity terms, respectively. The right hand side term $F_t$ is an excitation fluid force, considered as independent on the tube motion. According to this forced oscillator model, the flexible tube turns unstable when the total damping term $C_T = C_s + C_f$ vanishes, i.e. $C_T = 0$. Here, we consider that $C_T$ is a function $H_{CT}^{(j)}$ of $(C_j, F_j, \rho_1, \rho_2, \nu_1, \nu_2, \sigma, J_1, J_2, P, D, L)$ so that we can write

$$ C_T = H_{CT}^{(j)} (C_j, F_j, \rho_1, \rho_2, \nu_1, \nu_2, \sigma, J_1, J_2, P, D, L), \quad (2) $$

where $C_j$ and $F_j$ are the damping coefficient and frequency of the flexible tube immersed in the single-phase fluid $j$ at rest. The relation (2) constitutes a minimal model, based on experimental observations, bibliography reporting and physical intuition. More advanced models would also consider the effect of some other parameters, for e.g. the roughness of the tubes. The dimensional analysis is based on the Vaschy-Buckingham theorem. The theorem states that an equation involving $n$ physical variables with $k$ fundamental units (usually $k = 3$ in classical mechanics) can be reduced to an equation involving $n-k$ dimensionless parameters. Thus, introducing a scale of length $\lambda$, mass $M$ and time $\tau$, the equation for the instability threshold $C_T = 0$ is physically meaningful if

$$ h_{CT}^{(j)} \left( \frac{C_j}{M\tau^{-1}}, \frac{F_j}{\tau^{-1}}, \frac{\rho_1}{M\lambda^{-3}}, \frac{\rho_2}{M\lambda^{-3}}, \frac{\nu_1}{\lambda^2\tau^{-1}}, \frac{\nu_2}{\lambda^2\tau^{-1}}, \frac{\sigma}{M\tau^{-2}}, \frac{J_1}{\lambda\tau^{-1}}, \frac{J_2}{\lambda\tau^{-1}}, \frac{P}{\lambda}, \frac{D}{\lambda}, \frac{L}{\lambda} \right) = 0. \quad (3) $$
As Eq. (3) involves twelve dimensional quantities with three fundamental dimensions (length, mass and time), it can be reduced to a relation between nine dimensionless quantities. These dimensionless quantities are not unique and derive from a specific choice for the characteristic length $\lambda$, mass $M$ and time $\tau$. Picking $\lambda = D$, $M = \rho_j D^2 L$ and $\tau = D/J_j$ with $j$ a free index, the dimensionless equation for the instability threshold can be rewritten as

$$h_{CT}^{(j)}(Sc_j, Sk_j, m, \varepsilon, We_j, Re_1, Re_2, p, l) = 0, \quad (4)$$

with $l = L/D$, $p = P/D$, $m = \rho_2/\rho_1$ and $\varepsilon = \nu_2/\nu_1$ the tube aspect ratio, the pitch ratio, the density and viscosity ratios, respectively. The other dimensionless numbers appearing in Eq. (4), i.e.

$$Re_j = \frac{DJ_j}{\nu_j}, \quad Sk_j = \frac{D^2 F_j}{\nu_j}, \quad Sc_j = \frac{C_j}{2\rho_j D^2 L F_j}, \quad We_j = \frac{\rho_j L \nu_j^2}{\sigma D^2}, \quad (5)$$

are the Reynolds, Stokes, Scruton (i.e. mass-damping parameter) and Weber numbers, constructed from the physical properties of phase $j$. Once again, the choice for the dimensionless numbers is purely arbitrary, and other combinations are possible. For example, instead of $Re_1$ or $Re_2$, it is also possible to introduce the homogeneous void fraction $\alpha$, via the relation $\alpha = \varepsilon Re_2 / (Re_1 + \varepsilon Re_2)$. Similarly, one can replace in Eq. (4) the Reynolds number $Re_j$ by a reduced velocity, $V_{rj} = Re_j/Sk_j$.

### 5.2.1 Single-phase flow instability threshold

A single-phase flow can be seen as a two-phase flow with two identical phases, i.e. $(\rho_1, \nu_1) = (\rho_2, \nu_2) = (\rho, \nu)$ and $\sigma = 0$. The density and viscosity ratios are then unitary, $m = 1$, $\varepsilon = 1$, while the Weber number diverges, $We \to \infty$. In this case, the stability threshold (4) simplifies to

$$h_{CT}^{(j)}(Sc, Sk, Re_1, Re_2, p, l) = 0, \quad (6)$$

with $Re_j = DJ_j/\nu$, $Sk = D^2 F/\nu$, $Sc = C / (2\rho D^2 L F)$, $C = C_1 = C_2$ and $F = F_1 = F_2$. As the two phases are identical, it is expected that Eq. (6) is in fact not a function of the two independent dimensionless numbers $Re_1$ and $Re_2$, but a function of the incident Reynolds number $Re = Re_1 + Re_2 = D(J_1 + J_2)/\nu$. Consequently, Eq. (6) reduces to an equation involving only five dimensionless numbers

$$h_{CT}(Sc, Sk, Re, p, l) = 0. \quad (7)$$

For a given geometry of the bundle and for some given material properties of the tube and the fluid, the dimensionless numbers $Sc$, $Sk$, $p$ and $l$ are fixed and the instability threshold is the set of Reynolds numbers $Re$ such that $Re = \text{constant}$. In the plane $(Re_1, Re_2)$, this threshold corresponds to a straight line of equation $Re = Re_1 + Re_2 = \text{constant} = Re_{\text{crit}}^{1\phi}$, i.e.

$$\frac{Re_1}{Re_{\text{crit}}^{1\phi}} + \frac{Re_2}{Re_{\text{crit}}^{1\phi}} = 1, \quad (8)$$
with $Re_{\text{crit}}^{1\phi} = H_{Re_{\text{crit}}^{1\phi}}(Sc, Sk, p, l)$ some function of $(Sc, Sk, p, l)$, according to Eq. (7). Note that the instability threshold proposed by Connors (1970) is actually a special case of Eq. (7), corresponding to $h_{C_T}(Sc, Sk, Re, p, l) = Re - K(p, l)^{-1/2}Sc^{1/2}Sk = 0$, which is also similar to $H_{Re_{\text{crit}}^{1\phi}}(Sc, Sk, p, l) = K(p, l)^{-1/2}Sc^{1/2}Sk$. In this expression $K$ is the Connors' constant, which is actually a function of $(p, l)$, and possibly some other dimensionless numbers characteristics of the geometry of the tube bundle.

### 5.2.2 Two-phase flow instability threshold

The two-phase flow instability threshold is given by the implicit equation (4). For a given geometry of the bundle and for some given material properties of the tube and the fluid, the dimensionless numbers $Sc_j, Sk_j, m, \epsilon, We_j, p$ and $l$ are fixed and the instability threshold is the set of Reynolds numbers $Re_1$ and $Re_2$ such that $h_{C_T}^{(j)} = 0$. In the plane $(Re_1, Re_2)$, this set corresponds to a curve whose shape is unknown. By similarity with single-phase flows, we suppose that this curve is also a straight line. The hypothesis of continuity of the threshold between single-phase and two-phase flows requires this line to pass through the points $(Re_1, Re_2) = (0, Re_{2, \text{crit}}^{1\phi})$ and $(Re_1, Re_2) = (Re_{1, \text{crit}}^{1\phi}, 0)$, where $Re_{j, \text{crit}}^{1\phi}$ is the instability threshold in single-phase flows for the phase labeled $j$. These two hypotheses imply that Eq. (4) is finally independent of $m, \epsilon$ and $We_j$, and is explicitly expressed as

$$ \frac{Re_1}{Re_{1, \text{crit}}^{1\phi}} + \frac{Re_2}{Re_{2, \text{crit}}^{1\phi}} = 1, \qquad (9) $$

with $Re_{j, \text{crit}}^{1\phi} = H_{Re_{\text{crit}}^{1\phi}}(Sc_j, Sk_j, p, l)$.

As expected (the stability threshold should not depend on how the label of the two phases), the relation (9) is invariant per index change $1 \leftrightarrow 2$ and establishes continuity with the single-phase thresholds, via the respective limits as $Re_2 \to 0$ and $Re_1 \to 0$

$$ Re_1 \to Re_{1, \text{crit}}^{1\phi} \quad \text{as} \quad Re_2 \to 0, \qquad (10a) $$

$$ Re_2 \to Re_{2, \text{crit}}^{1\phi} \quad \text{as} \quad Re_1 \to 0. \qquad (10b) $$

From the relation $\alpha = \epsilon Re_2 / (Re_1 + \epsilon Re_2)$, it is also possible to express the instability threshold in a form equivalent to Eq. (9), explicitly showing the homogeneous void fraction $\alpha$ and the viscosity ratio $\epsilon$

$$ Re_1 = \frac{Re_{1, \text{crit}}^{1\phi}}{1 + \frac{\alpha}{1-\alpha} \frac{1}{\epsilon} \frac{Re_{1, \text{crit}}^{1\phi}}{Re_{2, \text{crit}}^{1\phi}}}, \quad Re_2 = \frac{Re_{2, \text{crit}}^{1\phi}}{1 + \frac{1-\alpha}{\alpha} \epsilon \frac{Re_{2, \text{crit}}^{1\phi}}{Re_{1, \text{crit}}^{1\phi}}}. \qquad (11) $$

The above relationships remain invariant per index change $1 \leftrightarrow 2$ for $\alpha \to 1-\alpha$ and $\epsilon \to 1/\epsilon$. Similarly, the continuities with the single-phase thresholds
<figure><img src="image_26.png" /><figcaption>Figure 5.5: Left: 5 × 3 square tube bundle (pitch ratio $p = 1.5$) with a central tube mounted on two parallel blades, allowing vibrations in the lift direction only. Right: experimental setup and instrumentation. In the direct method of measurement of fluid-elastic forces, the PRODERA actuator is used to impose the harmonic vibration of the central tube.</figcaption></figure>

are maintained as $\alpha \to 0$ and $\alpha \to 1$

$$ Re_1 \to \begin{cases} Re_{1, \text{crit}}^{1\phi} & \text{as } \alpha \to 0 \\ 0 & \text{as } \alpha \to 1 \end{cases}, \quad Re_2 \to \begin{cases} 0 & \text{as } \alpha \to 0 \\ Re_{2, \text{crit}}^{1\phi} & \text{as } \alpha \to 1 \end{cases}. \qquad (12) $$

## 5.3 Single-phase flow experiments

The formulation (9) expresses the two-phase flow instability threshold as a function of the two single-phase flow thresholds, $Re_{1, \text{crit}}^{(1\phi)}$ and $Re_{2, \text{crit}}^{(1\phi)}$, images of $(Sc_1, Sk_1, p, l)$ and $(Sc_2, Sk_2, p, l)$ by the unknown function $H_{Re_{crit}^{1\phi}}$. To complete the formulation (9) it is therefore necessary to determine the function $H_{Re_{crit}^{1\phi}}$, in particular by establishing its dependence on the Scruton and Stokes numbers. To do so, we measure experimentally the fluid-elastic forces acting on a flexible tube that is imposed a harmonic displacement (direct method of Tanaka & Takahara (1980), see also Lagrange et al. (2019)) and subject to a single-phase (water) cross-flow. Measurements are performed for various flow rates $Q$, i.e. different Reynolds numbers, and different forcing frequencies, i.e. different Stokes numbers. From the measured fluid-forces, the total damping coefficient $C_T$ is extracted and the critical Scruton number, given by $C_T = 0$, is determined. This method makes it possible to map the instability threshold by independently varying the Reynolds number and the Stokes number.

### 5.3.1 Experimental setup and instrumentation

The experimental setup is sketched in Fig. 5.5 and described briefly in the following. Readers should refer to Caillaud et al. (1999, 2000, 2003); Piteau et al. (2012, 2019) for an extensive description of the experimental setup and to Lagrange et al. (2019) for a discussion on the methodology used to measure the fluid-elastic forces. The flexible tube is in the central
position of a 5 × 3 square bundle (pitch ratio $p = 1.5$) made of rigid tubes (plus two columns of 5 half-tubes at the boundaries). The tube bundle has immersed length $L = 300$ mm, diameter $D = 30$ mm, upstream cross section area $S = 6DL = 0.054$ m² and is subject to a water cross-flow (flow rate $Q$). The central tube is supported by two parallel flexible steel blades allowing vibrations in the lift direction. A PRODERA actuator is used to force the harmonic motion of the central tube (the forcing angular frequency $\omega$ and amplitude of the motion are imposed). A non-contact electrodynamic actuator is used for active speed control (system stability control of the first mode of vibration, via a feedback stabilizing force). The displacement and velocity of the flexible tube are provided by a KEYENCE laser sensor. The resultant force on the flexible tube is measured only in the lift direction by KISTLER sensors. This force, noted $F_m$, integrates the inertia, damping of the tube and the fluid-elastic forces, i.e.

$$F_m = (M_s + M_f) \frac{d^2 X}{dT^2} + (C_s + C_f) \frac{dX}{dT} + K_f X. \quad (13)$$

It follows that the transfer function between the measured force and the imposed displacement is

$$H(\omega) = -\omega^2 (M_s + M_f) + i\omega (C_s + C_f) + K_f, \quad (14)$$

so that $C_T = C_s + C_f = \Im\{H\}/\omega$, with $\Im$ the imaginary part operator. Let $C_0$ be the value of $C_T$ measured experimentally for a fluid at rest. Introducing $\bar{C}_f = C_T - C_0$, the critical Scruton number, defined by $C_T = 0$, is calculated as

$$Sc_{\text{crit}} = -\frac{\bar{C}_f}{2\rho F D^2 L}. \quad (15)$$

Varying the Reynolds number, $Re = DQ/(S\nu)$, and the Stokes number, $Sk = D^2F/\nu$, we eventually obtain a set of points ($Re, Sk, Sc_{\text{crit}}$) leading to a mapping of the instability threshold $h_{C_T}(Sc, Sk, Re, p, l) = 0$, and by extension of the unknown function $H_{Re}^{1\phi}$. In our experiments, measurements are performed for $Re \in \{27, 36, 45, 54, 72, 90, 117\} \times 10^3$. For each $Re$, 37 Stokes numbers, $Sk \in [3600, 36000]$, are tested (by continuous frequency sweep between 4 and 40 Hz).

### 5.3.2 Results and discussion

Figure 5.7 shows the evolution of the instability threshold, measured experimentally, represented in the planes ($Sk, Sc$) and ($Sc, V_r = Re/Sk$), for $Re \in \{27, 36, 45, 54, 72, 90, 117\} \times 10^3$. The experimental data show that (i) for a given Scruton number, a tube is all the more stable (in the sense that the Reynolds number must be increased to reach the instability threshold) as the Stokes number is large, (ii) for a given Stokes number, a tube is all the more stable as the Scruton number is large, (iii) for a given Reynolds number, if the Scruton number is increased, then the Stokes number has be decreased to turn the flexible tube unstable, (iv) the instability threshold function $H_{Re}^{1\phi}$ is discontinuous at $Sc \approx 0.28$. All these observations
<figure><img src="image_27.png" /><figcaption>Figure 5.6: Log-Log graph of the instability threshold for a water cross-flow. The log-log scale is used to identify different power laws for the instability threshold, depending on the Scruton number. The points correspond to the experimental measurements. The curves correspond to the regression (17).</figcaption></figure>

actually show that a tube is all the more stable as it is rigid and/or as its reduced mass-damping parameter (i.e. Scruton number) is large. Also, our experimental results show that the instability threshold is at least a function of three dimensionless parameters, for e.g. (Sc, Sk, Re), which can not be reduced to a function of only two parameters, for. e.g. Sc and $V_r = Re/Sk$, as suggested by Connors (1970) through the law $V_{rcrit} = K(p, l)Sc^{1/2}$. In other words, to correctly predict the instability threshold of a flexible tube subject to a single-phase flow, it is necessary to dissociate the fluid flow effects, represented by the Reynolds number, which will modify the fluid damping, and therefore affect the instability threshold, from structural frequency effects, represented by the Stokes number.

In order to obtain an analytical expression of the instability threshold, a non-linear regression of all the experimental points is carried out, in the form of two power laws, i.e.

$$ Sc(d, b_1, b_2, b_3, c_1, c_2, c_3) = \begin{cases} Sc_1 = b_1 Re^{b_2} Sk^{b_3} & \text{if } Sc < d \\ Sc_2 = c_1 Re^{c_2} Sk^{c_3} & \text{if } Sc > d \end{cases}, \quad (16) $$

with $b_i$, $c_i$ and $d$ the fitting parameters. The discontinuity at $Sc = d$ is a pure mathematical artifact introduced to fit as much data as possible. This discontinuity is only used to take into account the possibility of different instability mechanisms for small and large Scruton numbers, as suggested by Chen (1983b). Indeed, a log-log plot of the data, see Fig. 5.6, tends to indicate the existence of different power laws for the instability threshold, depending on the Scruton number. At first glance, it seems reasonable to identify two regimes, with two different slopes. Obviously, a continuous non-
<figure><img src="image_28.png" /><figcaption>Figure 5.7: Graphs of the instability threshold for a water cross-flow. The points correspond to the experimental measurements. The curves correspond to the regressions (17) and (19). Subfigures (c) and (d) are representative diagrams of subfigures (a) and (b), respectively.</figcaption></figure>
linear regression could be used to adjust the instability threshold, for example by imposing an additional constraint $Sc_1(d) = Sc_2(d)$ in the minimization process. This would make it possible to avoid introducing a non-physical discontinuity into the model, but this would also tend to unnecessarily complicate the expression of the instability threshold, in single-phase flow, and by extension in two-phase flow. Thus, in our study, guided by engineering applications, and following the example of Chen (1984), we stick to the model given by Eq. (16) (i.e. model based on power laws with a discontinuity), which has the merit of being simple. From the minimization process (least squares method), the discontinuity is obtained at $Sc = 0.28$ and the critical threshold writes

$$ Sc_{crit} = \begin{cases} 9.78 \times 10^5 Re^{4.94} Sk^{-7.29} & \text{if } Sc < 0.28 \\ 67.5 Re^{3/2} Sk^{-9/4} & \text{if } Sc > 0.28 \end{cases}, \quad (17) $$

or equivalently

$$ Re_{crit} = \begin{cases} 6.13 \times 10^{-2} Sc^{0.2} Sk^{1.48} & \text{if } Sc < 0.28 \\ 6.03 \times 10^{-2} Sc^{2/3} Sk^{3/2} & \text{if } Sc > 0.28 \end{cases}, \quad (18) $$

or equivalently

$$ V_{r_{crit}} = \begin{cases} 0.151 Sc^{0.137} Re^{0.322} & \text{if } Sc < 0.28 \\ 0.154 Sc^{4/9} Re^{1/3} & \text{if } Sc > 0.28 \end{cases}. \quad (19) $$

In Fig. 5.7, we show that the power laws (17) and (19) are in good agreement with experimental data. We conclude that Eq. (18) is a good prediction of the instability threshold in single-phase flow, that could be inserted in Eq. (9) to complete the criterion for the instability threshold in two-phase flow. Before going further, it should be noted that it is in fact possible to determine the total damping of the tube, $C_T$, for any Reynolds number, provided that $Re_{crit}$ is known. Indeed, let $k$ be an iso-value of the dimensionless coefficient $c_T = C_T/C_0$, so that

$$ c_T = \frac{C_0 + \overline{C_f}}{C_0} = k \Leftrightarrow \frac{(1-k) C_0 + \overline{C_f}}{C_0} = 0. \quad (20) $$

This equation indicates that the point of coordinates $(Re, c_T = k)$, associated to the tube with modal damping $C_0$, corresponds to the critical threshold of a tube with modal damping $(1-k)C_0$, noted $Re_{crit}$, i.e.

$$ Re(C_0, c_T = k) = Re_{crit}((1-k)C_0). \quad (21) $$

In Eq. (21), the notation $Re(C_0, c_T = k)$ is used to define the Reynolds number for which the dimensionless total damping coefficient $c_T$ of the tube with modal damping $C_0$ in still fluid equals $k$. According to Eq. (18), the instability threshold scales as $Re_{crit} \sim Sc^a$ with $Sc \propto C_0$ and $a = 0.2$ (resp. $2/3$) if $Sc < 0.28$ (resp. $Sc > 0.28$). Therefore, the right hand side term of Eq. (21) writes

$$ Re_{crit}((1-k)C_0) = (1-k)^a Re_{crit}(C_0), \quad (22) $$
and Eq. (21) simplifies to

$$Re(C_0, c_T = k) = (1 - k)^a Re_{\text{crit}}(C_0). \quad (23)$$

Replacing $k$ by $c_T$ and $Re(C_0, c_T = k)$ by $Re$ in Eq. (23), the coefficient $c_T$ can be written as

$$c_T = 1 - \left( \frac{Re}{Re_{\text{crit}}} \right)^{1/a} . \quad (24)$$

Finally, plugging Eq. (18) into Eq. (24) yields an expression of the total damping coefficient in single-phase flow

$$c_T = \begin{cases} 1 - \left( \frac{Re_1}{6.13 \times 10^{-2} Sc^{0.2} Sk^{1.48}} \right)^5 & \text{if } Sc < 0.28 \\ 1 - \left( \frac{Re_1}{6.03 \times 10^{-2} Sc^{2/3} Sk^{3/2}} \right)^{3/2} & \text{if } Sc > 0.28 \end{cases} . \quad (25)$$

Figure 5.8 shows the evolution of $c_T$ as a function of $Re$ for a set of blades with Scruton and Stokes numbers $Sc = 0.19$ and $Sk = 1.86 \times 10^4$ (blades with a frequency of 20.8 Hz in still fluid). The experimental points in Fig. 5.8 come from Lagrange et al. (2019) where $c_T$ is obtained from the modal response of the tube mounted on the blades (indirect method of measuring forces). In Lagrange et al. (2019), experiments are performed well beyond the instability boundary of the system, using active control, allowing exploration of a significant range of Reynolds numbers. For $Sc < 0.28$, the proposed law (25) predicts quite satisfactorily the variations of $c_T$ (as well as the prediction of the instability threshold $c_T = 0$) as a function of the Reynolds number, with the exception of the increasing part observed experimentally for the moderate Reynolds numbers, characteristic of tube stabilization. Since the approach is only valid for $k \le 1$, the law does not in fact make it possible to predict a coefficient $c_T > 1$. For the interest of the reader, we have also plotted in Fig. 5.8 the variations of Eq. (25), considering $Sc > 0.28$. We show that the variations of $c_T$ (as well as the prediction of the instability threshold $c_T = 0$) are not in good agreement with the experimental data, suggesting that it is better to fit the instability threshold with two power laws, possibly by introducing an artificial discontinuity, as done in our study.

# 5.4 Two-phase flow experiments

## 5.4.1 New criterion for the instability threshold

Introducing Eq. (18) into Eq. (9) yields a new criterion for the instability threshold in two-phase flow. This criterion writes

$$\frac{Re_1}{Re_{1, \text{crit}}^{(1\phi)}} + \frac{Re_2}{Re_{2, \text{crit}}^{(1\phi)}} = 1, \quad (26)$$

with

$$Re_{i, \text{crit}}^{(1\phi)} = \begin{cases} 6.13 \times 10^{-2} Sc_i^{0.2} Sk_i^{1.48} & \text{if } Sc_i < 0.28 \\ 6.03 \times 10^{-2} Sc_i^{2/3} Sk_i^{3/2} & \text{if } Sc_i > 0.28 \end{cases}. \quad (27)$$

To assess the validity of Eq. (26), along with Eq. (27), experimental measurements of the instability threshold in air-water cross-flow are carried out.
<figure><img src="image_29.png" /><figcaption>Figure 5.8: Evolution of the dimensionless total damping coefficient, $c_T$, as a function of the Reynolds number, for a set of blades characterized by $Sc = 0.19$ and $Sk = 1.86 \times 10^4$ (blades with a frequency of 20.8 Hz in still fluid). The points correspond to the experimental data of Lagrange et al. (2019), obtained with the indirect method and an active control system. The continuous (resp. dashed) line corresponds to the theoretical prediction (25) for $Sc < 0.28$ (resp. $Sc > 0.28$).</figcaption></figure>

### 5.4.2 Experimental setup and methodology

The experimental setup for measurements in two-phase flow is the same as the one used in single-phase flow, see Fig. 5.5, except that the central flexible tube is not imposed a harmonic motion. Here, the tube is free to vibrate and the fluid-elastic force is extracted from the changes in the modal frequency and damping of the tube. This indirect measurement method of the fluid-elastic force is discussed in details in Caillaud et al. (1999, 2000, 2003); Piteau et al. (2012, 2019). The modal frequency of the central tube is set by changing the thickness of the blades. Three sets of blades, denoted L1C, L4C and L2C, whose modal properties in still air and water are reported in Table 5.1, have been tested. The tube bundle is subject to an air-water cross-flow with flow rates in the ranges $Q_{water} = [3, 300]$ L/s and $Q_{air} = [20, 3000]$ m³/h. The measurement of the displacement and velocity of the flexible tube is provided by a laser transducer (Keyence LK-G500). A feedback control loop powered by an electro-dynamical shaker is used to ensure the stability of the tube, whatever the flow conditions. The principle of experimental determination of the instability threshold is as follows. For a given homogeneous void fraction, a modal identification is used to determine the reduced damping $\xi_1$ and the frequency $F_1$ of the coupled fluid-structure-shaker system. The total reduced damping $\xi_T$ of the fluid-structure system is obtained by subtracting the control damping $C_{ctrl}$ due to the electro-dynamical shaker, by defining a mass of reference.
<table><caption>Table 5.1: Characteristics of the three sets of blades L1C, L4C and L2C.</caption><thead><tr><th>Blades</th><th>L1C</th><th>L4C</th><th>L2C</th></tr></thead><tbody><tr><td>M<sub>air</sub> (kg)</td><td>1.29</td><td>1.25</td><td>1.22</td></tr><tr><td>M<sub>water</sub> (kg)</td><td>1.67</td><td>1.64</td><td>1.57</td></tr><tr><td>ξ<sub>air</sub> (%)</td><td>0.56</td><td>0.46</td><td>0.38</td></tr><tr><td>ξ<sub>water</sub> (%)</td><td>1.14</td><td>0.89</td><td>0.77</td></tr><tr><td>F<sub>air</sub> (Hz)</td><td>14.11</td><td>22.06</td><td>29.18</td></tr><tr><td>F<sub>water</sub> (Hz)</td><td>11.19</td><td>18.53</td><td>24.9</td></tr><tr><td>Sk<sub>air</sub> (×10<sup>4</sup>)</td><td>0.08140</td><td>0.12727</td><td>0.1683</td></tr><tr><td>Sk<sub>water</sub> (×10<sup>4</sup>)</td><td>1.000</td><td>1.656</td><td>2.225</td></tr><tr><td>Sc<sub>air</sub></td><td>140.1</td><td>111.5</td><td>89.90</td></tr><tr><td>Sc<sub>water</sub></td><td>0.4430</td><td>0.3397</td><td>0.2813</td></tr></tbody></table>

In our case, two choices for this mass have been tested. The first choice is based on the homogeneous mass $M_h = M_{air} + \rho_h (M_{water} - M_{air}) / \rho_{water}$ with $\rho_h = \alpha \rho_{air} + (1 - \alpha) \rho_{water}$. The second choice is based on a mass measurement $M_{2\phi}$ obtained from the frequency $F_{2\phi}$ of the tube motion for very low flow rates of water $Q_{water}$ and air $Q_{air}$. Figure 5.9 shows that at first order the mass ratio $M_{ratio} = (M_{2\phi} - M_{air})/(M_{water} - M_{air})$ decreases linearly with the homogeneous void fraction, $\alpha$. Similar trends were observed by Pettigrew et al. (1989); Carlucci (1980). Just like Pettigrew et al. (1989), however, we observe an oscillation of the mass ratio, which remains lower than $\rho_h / \rho_{eau} = (M_h - M_{air})/(M_{water} - M_{air})$ for $\alpha < 0.5$, then becomes higher for higher void fractions. The nature of this oscillation remains unknown but we suspect that it may be related to a non-uniform distribution of the local void fraction close to the tube and to a change in the flow regime for $\alpha \approx 0.5$. The two choices for a mass reference, $M_h$ or $M_{2\phi}$, give two estimates of reduced damping, $\xi_T = \xi_{Th} = \xi_1 - C_{ctrl}/(2M_h(2\pi F_1))$ or $\xi_T = \xi_{TMm} = \xi_1 - C_{ctrl}/(2M_{2\phi}(2\pi F_1))$, which ultimately differ little from each other, i.e. $\xi_T \approx \xi_{Th} \approx \xi_{TMm}$. The instability threshold is then defined as the value of $Re = (Re_{water}^2 + Re_{air}^2)^{1/2}$ for which $\xi_{TMm} = 0$. In order to determine this critical value $Re_{crit}$, a linear interpolation between the last stable measuring point ($\xi_{TMm} > 0$) and the first unstable measuring point ($\xi_{TMm} < 0$) is performed. The repetition of this process for different void fractions $\alpha$ yields a set of points $(\alpha, Re_{crit})$. Once $Re_{crit}$ has been determined, the results are reported in the planes $(Re_{air}, Re_{water})$, $(\alpha, Re_{air})$ and $(\alpha, Re_{water})$ with $\alpha = \varepsilon Re_{air} / (Re_{water} + \varepsilon Re_{air})$.

### 5.4.3 Results and discussion

Figure 5.10 shows graphs of the instability threshold in air-water cross-flow, represented in planes $(Re_{air}, Re_{water})$, $(Re_{air}/Re_{air,crit}^{1\phi}, Re_{water}/Re_{water,crit}^{1\phi})$, $(\alpha, Re_{air})$ and $(\alpha, Re_{water})$, for the three sets of blades, L1C, L4C and L2C. These blades are characterized by Scruton numbers $Sc_{air}$ and $Sc_{water}$ greater
<figure><img src="image_30.png" /><figcaption>Figure 5.9: Evolution of the mass ratio $M_{ratio} = (M_{2\phi} - M_{air})/(M_{water} - M_{air})$ as a function of the homogeneous void fraction, $\alpha$, for the three sets of blades tested, L1C, L4C and L2C. The straight line corresponds to the homogeneous theory, $M_{ratio} = \rho_h/\rho_{water}$, with $\rho_h = \alpha\rho_{air}+(1-\alpha)\rho_{water}$ the homogeneous mass density. Square symbols are experimental data from Pettigrew et al. (1989). Triangle symbols are experimental data from Carlucci (1980).</figcaption></figure>

than 0.28, see Table 5.1, so that the criterion (26) reduces to

$$ \frac{Re_{water}}{Sc_{water}^{2/3} Sk_{water}^{3/2}} + \frac{Re_{air}}{Sc_{air}^{2/3} Sk_{air}^{3/2}} = 6.03 \times 10^{-2}, \quad (28) $$

or equivalently

$$ Re_{air,crit} = 6.03 \times 10^{-2} \frac{Sc_{air}^{2/3} Sk_{air}^{3/2}}{1 + \frac{1-\alpha}{\alpha} \frac{\nu_{air}}{\nu_{water}} \left( \frac{Sc_{air}}{Sc_{water}} \right)^{2/3} \left( \frac{Sk_{air}}{Sk_{water}} \right)^{3/2}}, \quad (29) $$

or equivalently

$$ Re_{water,crit} = 6.03 \times 10^{-2} \frac{Sc_{water}^{2/3} Sk_{water}^{3/2}}{1 + \frac{\alpha}{1-\alpha} \frac{\nu_{water}}{\nu_{air}} \left( \frac{Sc_{water}}{Sc_{air}} \right)^{2/3} \left( \frac{Sk_{water}}{Sk_{air}} \right)^{3/2}}. \quad (30) $$

In the plane $(Re_{air}, Re_{water})$, experimental results nearly collapse on a straight line, as suggested by Eq. (28). Lines for the three sets of blades are staggered according to the frequency of the blades, reflecting the effect of the Stokes number on the single-phase flow critical Reynolds number, see Eq. (27). In the plane $(Re_{air}/Re_{air,crit}^{1\phi}, Re_{water}/Re_{water,crit}^{1\phi})$, all experimental results align into a single line, as predicted by Eq. (26). Finally, in the plane $(\alpha, Re_{air})$ (resp. $(\alpha, Re_{water})$), the graph of the instability threshold is an increasing (resp. decreasing) function of $\alpha$, with $Re_{air,crit} \to 0$ as $\alpha \to 0$ and
<figure><img src="image_31.png" /><figcaption>Figure 5.10: Graphs of the instability threshold for an air-water cross-flow. The points correspond to the experimental results. The curves correspond to the predictions from the criterion (26), which reduces to Eq. (28) for the three sets of blades L1C, L4C and L2C.</figcaption></figure>
$Re_{air,crit} \rightarrow Re_{air,crit}^{1\phi}$ as $\alpha \rightarrow 1$ (resp. $Re_{water,crit} \rightarrow Re_{water,crit}^{1\phi}$ as $\alpha \rightarrow 0$ and $Re_{water,crit} \rightarrow 0$ as $\alpha \rightarrow 1$), as predicted by Eq. (29) (resp. Eq. (30)).

At first sight, the criterion for the instability threshold seems to perform worse for $\alpha < 0.80$ than for higher $\alpha$. But this first impression actually depends on the plane of observation, the definition of the relative deviation, and the set of blades under consideration. To be more specific, in our experiments, the instability threshold is defined as the value of $Re = (Re_{water}^{2} + Re_{air}^{2})^{1/2}$ for which the total damping of the system vanishes. Therefore, to correctly assess the performance of the criterion, it is necessary to compare the experimental value of $Re_{crit}$, noted $Re_{crit}^{\text{Exp}}$, with its prediction, $Re_{crit}^{\text{Th}}$. In Fig. 5.11, we then show the evolution of the relative deviation, $\iota (\%) = 100|Re_{crit}^{\text{Exp}} - Re_{crit}^{\text{Th}}|/Re_{crit}^{\text{Exp}}$, as a function of the homogeneous void fraction. For the stiffest blades L2C and L4C, we observe that the maximum relative deviation (approximately 20%) is reached for intermediate homogeneous void fractions $\alpha \in [0.3, 0.6]$. For such values of $\alpha$, the regime map shown in Fig. 5.12 indicates that the instability threshold of these two sets of blades is located near a transition between a bubbly regime (identified by the probability density function (PDF) of a local void fraction, with a PDF characterized by a marked peak close to $\alpha = 0$) and a slug regime (PDF characterized by a peak close to $\alpha = 0$, which extends to higher values of $\alpha$). On the contrary, for the most flexible set of blades L1C, we observe that the relative deviation decreases and hits a minimum in the range $\alpha \in [0.3, 0.6]$. For this set of blades, the regime map indicates that the instability threshold is located near a transition between a slug regime (PDF characterized by a peak close to $\alpha = 0$, which extends to higher values of $\alpha$) and a churn regime (PDF characterized by a peak close to $\alpha = 1$, which extends to lower values of $\alpha$). Interestingly, we also observe that the relative deviation increases in the range $\alpha \in [0.6, 1]$, which, according to the regime map, could be related to a transition from a churn to an annular regime (PDF characterized by a marked peak close to $\alpha = 1$). Let us note that the most flexible sets of blades, L1C and L4C, present a relative deviation greater than 15% for $\alpha = 0.99\%$. In fact, the existence of this peak is explained by the experimental difficulties associated with controlling a very low water flow (a few liters per second), while ensuring a high air flow. Also, additional measurements seem to indicate the possible existence of a cascade of bifurcations and hysteresis cycles, making the precise determination of the instability threshold very difficult for such high values of $\alpha$ and such flexible blades. Finally, despite some peaks in the relative deviation, the criterion (26), which does not take into account any tension surface effect, by including a Webber number for example, yields a rather accurate prediction of the instability threshold for all homogeneous voids fractions $\alpha \in [0.1]$.

## 5.5 Conclusion

In this work, we have studied the fluid-elastic instability of a tube bundle subject to a two-phase cross-flow. From a dimensional analysis, we have proposed a new criterion for the instability threshold, involving the
<figure><img src="image_32.png" /><figcaption>Figure 5.11: Left: geometrical interpretation of the relative deviation, $\iota (\%) = 100|Re_{crit}^{Exp} - Re_{crit}^{Th}|/Re_{crit}^{Exp}$. Right: evolution of the relative deviation, $\iota$, as a function of the homogeneous void fraction, $\alpha$.</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>Figure 5.12: Left: probability density functions (PDF) of the void fraction, $\alpha$, measured locally with an optical probe, in a rigid square tube bundle (5 x 3 full rigid square tube bundle, pitch ratio $p = 1.5$) subject to an air-water cross-flow. The optical probe is positioned halfway between the central cylinder (3, 2) and the downstream cylinder (2, 3). The superficial gap velocities are ($J_{water}, J_{air}$) = (5.4, 1.79) m/s (bubbly, $\alpha = 0.25$), ($J_{water}, J_{air}$) = (2.5, 5.86) m/s (slug, $\alpha = 0.70$), ($J_{water}, J_{air}$) = (1.7, 4) m/s (churn, $\alpha = 0.70$), ($J_{water}, J_{air}$) = (0.8, 22.1) m/s (annular, $\alpha = 0.965$). Right: regime map. The color of points is related to the experimental identifications of flow regimes: bubbly (purple), slug (light blue), churn (dark blue), annular (yellow). The curves correspond to the instability threshold given by Eq. (26), for the three sets of blades L1C, L2C and L4C.</figcaption></figure>
Reynolds numbers of the two phases and the critical Reynolds number in single-phase flow. This new criterion has been established by similarity with single-phase flows, assuming that its representative curve is a straight line in the plane generated by the two Reynolds numbers. Experiments in water cross-flow (tube bundle with central tube mounted on two flexible blades, direct method of measurement of fluid-elastic forces) have shown that the critical Reynolds number in single-phase flow can be expressed in the form of a power law involving the Scruton and Stokes numbers. The instability criterion in two-phase flow has been corroborated with air-water experiments (indirect method of measurement of fluid-elastic forces). This criterion is shown to be predictive for homogeneous void fractions from 0 to 1, and so for all flow regimes (identified in experiments with an optical probe). The validity of the criterion has now to be tested for different tube bundle geometries (e.g. triangular pitch) and mixtures (e.g. water-freon and ultimately steam-water).

## Acknowledgements

This work was performed in the framework of a joint research program co-funded by FRAMATOME, EDF and CEA (France). A valuable contribution by M. Bellanger and T. Valin (CEA-Saclay) to the experimental work is gladly acknowledged.
Bibliography

AXISA, F., BOHEAS, M. A. & VILLARD, B. 1985 Vibration of tube bundles subjected to steam-water cross-flow: A comparative study of square and triangular arrays. In *8th International Conference on Structural Mechanics in Reactor Technology, Paper No. B1/2*.

AXISA, F., VILLARD, B., GIBERT, R. J., HETSRONI, G. & SUNDHEIMER, P. 1984 Vibration of tube bundles subjected to air-water and steam-water cross flow: Preliminary results on fluidelastic instability. In *Symposium on Flow-Induced Vibrations, Volume 2: Vibration of Arrays of Cylinders in Cross Flow. Presented at the 1984 ASME Winter Annual Meeting. ASME New York NY USA*.

CAILLAUD, S., DE LANGRE, E. & BAJ, F. 2003 Active vibration control for the measurement of fluidelastic effects. *ASME Journal of Pressure Vessel Technology* **125**, 165-170.

CAILLAUD, S., DE LANGRE, E. & PITEAU, P. 1999 The measurement of fluidelastic effects at low reduced velocities using piezoelectric actuators. *ASME Journal of Pressure Vessel Technology* **121**, 232-238.

CAILLAUD, S., DE LANGRE, E. & PITEAU, P. 2000 Measurement of critical velocities for fluidelastic instability using vibration control. *ASME Journal of Pressure Vessel Technology* **122**, 341-345.

CARLUCCI, N. L. 1980 Damping and hydrodynamic mass of a cylinder in simulated two-phase flow. *Journal of Mechanical Design, Transactions of the ASME* **102**, 597-602.

CHEN, S. S. 1983a Instability mechanisms and stability criteria of a group of circular cylinders subjected to cross-flow. part 1: Theory. *Journal of Vibration, Acoustics, Stress, and Reliability in Design* **105**, 51-58.

CHEN, S. S. 1983b Instability mechanisms and stability criteria of a group of circular cylinders subjected to cross-flow. part ii: Numerical results and discussion. *Journal of Vibration, Acoustics, Stress, and Reliability in Design* **105**, 253-260.

CHEN, S. S. 1984 Guidelines for the instability flow velocity of tube arrays in crossflow. *Journal of sound and Vibration* **93**, 439-455.

CHEN, S. S. 1987 A general theory for dynamic instability of tube arrays in cross-flow. *Journal of Fluids Engineering, Transaction of the ASME* **1**, 35-53.

CONNORS, H. 1970 Fluidelastic vibration of tube arrays excited by crossflows. *Flow-Induced Vibration in Heat Exchangers* pp. 42-56.

DALTON, C. & HELFINSTINE, R. A. 1971 Potential flow past a group of circular cylinders. *ASME Journal of Basic Engineering* **93**, 636-642.
GRANGER, S. & PAIDOUSSIS, M. P. 1996 An improvement to the quasi-
steady model with application to cross-flow-induced vibration of tube ar-
rays. *Journal of Fluid Mechanics* **320**, 163–184.

HASSAN, M. & WEAVER, D. S. 2015 The effect of flat bar supports on streamwise fluidelastic instability in heat exchanger tube arrays. *ASME. Journal of Pressure Vessel Technology* **137**, 061302.

HASSAN, M. & WEAVER, D. S. 2016 Modeling of streamwise and transverse fluidelastic instability in tube arrays. *ASME. Journal of Pressure Vessel Technology* **138**, 051304.

LAGRANGE, R., PANUNZIO, D., PITEAU, P., DELAUNE, X. & ANTUNES, J. 2022 A new criterion for the instability threshold of a square tube bundle subject to an air-water cross-flow. In *Proceedings of the 12th International Conference on Flow-Induced Vibration*, L. Zimmer, Laboratoire EM2C, CNRS, CentraleSupélec, Université Paris-Saclay edn. (ed. C. Habchi & P. Moussou), pp. 391–399.

LAGRANGE, R., PITEAU, P., DELAUNE, X. & ANTUNES, J. 2019 Fluid-elastic coefficients in single phase cross flow: Dimensional analysis, direct and indirect experimental methods. In *Pressure Vessels and Piping Conference*, vol. 4: Fluid-Structure Interaction, p. 93984.

LEVER, J. & WEAVER, D. S. 1986 On the stability of heat exchanger tube bundles, part 1: Modified theoretical model. *Journal of Sound and Vibration* **107**, 375–392.

NAKAMURA, T., FUJITA, K., KAWANISHI, K. & SAITO, I. 1986a Study on the flow induced vibration of a tube array by a two-phase flow (2nd report, large amplitude vibration by steam - water flow). *Japan Society of Mechanical Engineers Part C.* **52**, 2790–2795.

NAKAMURA, T., FUJITA, K., YAMAGUCHI, N., SAKATA, K., TSUGE, A. & SAITO, I. 1986b Study on flow induced vibration of a tube array by a two-phase flow. (1st report: large amplitude vibration by air-water flow). *Japan Society of Mechanical Engineers Part C.* **52**, 252–257.

PETTIGREW, M. J. & TAYLOR, C. E. 2009 Vibration of a normal triangular tube bundle subjected to two-phase freon cross flow. *Journal of Pressure Vessel Technology, Transactions of the ASME* **131**, 051302.

PETTIGREW, M. J., TAYLOR, C. E., JANZEN, V. P. & WHAN, T. 2002 Vibration behavior of rotated triangular tube bundles in two-phase cross flows. *Journal of Pressure Vessel Technology, Transactions of the ASME* **124**, 144–153.

PETTIGREW, M. J., TAYLOR, C. E., JONG, J. H. & CURRIE, I. G. 1995 Vibration of a tube bundle in two-phase freon cross-flow. *Journal of Pressure Vessel Technology, Transactions of the ASME* **117**, 321–329.
PETTIGREW, M. J., TAYLOR, C. E. & KIM, B. S. 1989 Vibration of tubes bundles in two-phase cross-flow: Part 1-hydrodynamic mass and damping. *Journal of Pressure Vessel Technology, Transactions of the ASME* **111**, 466-477.

PITEAU, P., DELAUNE, X., ANTUNES, J. & BORSOI, L. 2012 Experiments and computations of a loosely supported tube in a rigid bundle subjected to single-phase flow. *Journal of Fluids and Structures* **28**, 56-71.

PITEAU, P., DELAUNE, X., BORSOI, L. & ANTUNES, J. 2019 Experimental identification of fluid-elastic coupling forces on a flexible tube within a rigid square bundle subjected to single-phase cross-flow. *Journal of Fluids and Structures* **86**, 156-169.

PRICE, S. J. & PAIDOUSSIS, M. P. 1982 A theoretical investigation of the parameters affecting the fluidelastic instability of a double row of cylinders subject to a cross-flow. In *Vibration in Nuclear Plant. Proceedings of the 3rd International Conference. British Nucl. Energy Soc, London UK*.

PRICE, S. J. & PAIDOUSSIS, M. P. 1983 Fluidelastic instability of an infinite double row of circular cylinders subject to a uniform cross-flow. *Journal of Vibration Acoustics Stress Reliability in Design* **105**, 59-66.

ROBERTS, B. W. 1962 Low frequency, self-excited vibration in a row of circular cylinders mounted in an airstream. *Ph.D Thesis, University of Cambridge, Cambridge, UK*.

SADEK, O., MOHANY, A. & HASSAN, M. 2018 Numerical investigation of the cross flow fluidelastic forces of two-phase flow in tube bundle. *Journal of Fluids and Structures* **79**, 171-186.

SADEK, O., MOHANY, A. & HASSAN, M. 2020 Numerical simulation of streamwise fluidelastic instability of tube bundles subjected to two-phase cross flow. *Journal of Fluids and Structures* **92**, 102816.

TANAKA, H. & TAKAHARA, S. 1980 Unsteady fluid dynamic force on tube bundle and its dynamic effect on vibration. *American Society of Mechanical Engineers, Pressure Vessels and Piping Division (Publication) PVP* pp. 77-92.

TANAKA, H. & TAKAHARA, S. 1981 Fluid elastic vibration of tube array in cross-flow. *Journal of Sound and Vibration* **77**, 19-37.

YASSER, S., HASSAN, M., MOHANY, A. & AHMED, W. 2021 Modelling of fluidelastic instability in tube bundles under two-phase bubbly flow conditions. *Journal of Fluids and Structures* **103**, 103256.
# 6. PERSPECTIVES DE RECHERCHE

Dans ce mémoire ont été présentés quelques un de mes travaux de recherche post-doctoraux sur les instabilités en mécanique des fluides et des solides. Afin de conclure ce mémoire, je mets l’accent sur quelques projets qui seront privilégiés dans les années à venir. En tant qu’ingénieur-chercheur au Commissariat à l’Energie Atomique (C.E.A.), ces projets s’inscrivent dans le cadre d’une problématique industrielle, en partenariat avec les grands acteurs Français du nucléaire (E.D.F., Framatome), avec pour ligne directrice la modélisation fine des phénomènes d’interaction fluide-structure (I.F.S.). D’un point de vue recherche amont, ces travaux visent à améliorer la compréhension des mécanismes physiques de l’I.F.S.. D’un point de vue applicatif, ces travaux doivent alimenter les codes de calcul (par exemple TrioCFD<fnref n="1" />, Neptune-CFD<fnref n="2" />) développés par le C.E.A. et ses partenaires, afin de disposer d’outils numériques fiables, qualifiés et prédictifs. In fine, ces outils doivent permettre de décrire efficacement le comportement mécanique de composants structurels utilisés dans le domaine nucléaire, pour améliorer leur dimensionnement, évaluer leurs marges et leurs incertitudes de fonctionnement. Dans ce contexte, j’ambitionne donc de poursuivre mon activité de recherche sur l’I.F.S. des écoulements transverses (monophasiques ou diphasiques), dans la continuité du travail présenté dans le chapitre 5, sur l’instabilité fluide-élastique d’un faisceau de tubes. Je souhaite également développer une activité sur l’I.F.S. des écoulements monophasiques axiaux, présents en partie basse des tubes de G.V., dans les assemblages combustibles, dans les in-

<footnote n="1">¹ TrioCFD est un code de calcul numérique développé au C.E.A. par le Laboratoire de Modélisation et Simulation en mécanique des Fluides (L.M.S.F.). Ce code, massivement parallèle, est principalement dédié à la thermohydraulique monophasique des réacteurs nucléaires, via la résolution des équations de Navier-Stokes incompressibles ou faiblement compressibles, éventuellement couplées à un modèle de turbulence. Une collaboration entre le laboratoire DYN (R. Lagrange) et le L.M.S.F. (M. A. Puscas) est en cours pour développer et éprouver les capacités de TrioCFD et de son module A.L.E. (Arbitrary Lagrangian-Eulerian) pour la simulation de problèmes de mécanique impliquant un couplage entre des structures vibrantes et un fluide monophasique (au repos ou en écoulement).</footnote>

<footnote n="2">² Neptune-CFD est un code de calcul numérique développé par E.D.F., le C.E.A., Framatome et l’Institut de Radioprotection et de Sûreté Nucléaire, pour la simulation des écoulements multiphasiques. Ce code repose sur la résolution des équations de Navier-Stokes pour chaque phase constitutive du fluide. Ces équations sont complétées par des modèles théoriques et/ou phénoménologiques permettant de décrire l’interaction entre les phases. Une collaboration (thèse de C. Bazin, 2020-2023) entre le laboratoire DYN (R. Lagrange, encadrant C.E.A.), le L.M.S.F. (M. G. Rodio, encadrant C.E.A.) et le LadHyX (C. Josserand, directeur) est en cours pour assurer la validation de ces modèles, via la simulation des régimes d’écoulement et des phénomènes aléatoires observés dans la maquette expérimentale DIVA, voir chapitre 5.</footnote>
ternes de cuve<fnref n="3" />, mais également responsables des vibrations des tubes guides du réacteur Jules Horowitz<fnref n="4" />. Pour cela, je compte m’appuyer sur la résolution de problèmes fondamentaux (couplage hydrodynamique entre deux cylindres concentriques ou parallèles, vibrant dans un fluide au repos ou en écoulement), exigeants d’un point de vue théorique mais aussi numérique. Ces problèmes modèles doivent permettre de tester les capacités des codes numériques, définir des pistes de développement, mais aussi de réaliser des études paramétriques pour établir de nouvelles formulations analytiques, par exemple pour les coefficients de masse, d’amortissement, de rigidité ajoutés, qui déterminent la stabilité des systèmes couplés. Une présentation synthétique d’un programme de recherche pluriannuel est proposée dans les tables 6.1 (I.F.S. en écoulement axial monophasique) et 6.2 (I.F.S. en écoulement transverse monophasique ou diphasique). Les premiers résultats de ce programme de recherche sont illustrés sur les figures 6.1, 6.2 et des publications représentatives sont reportées en annexes C et D.

<table><caption>TABLE 6.1 : Programme de recherche sur l’I.F.S. en écoulement axial monophasique.</caption><thead><tr><th colspan="2">Hydrodynamic interaction between two coaxial cylinders</th></tr><tr><th colspan="2">Free vibrations</th></tr></thead><tbody><tr><th scope="row">Goal</th><td>- Simulate and model the free oscillations of a slender and confined structure vibrating in a still fluid</td></tr><tr><th scope="row">References</th><td>- R. Lagrange et al., J. Appl. Mech., 89(8) : 081006, 2022.<br>- R. Lagrange et al., J. Fluids Struct., 114, 103754, 2022.<br>- R. Lagrange et al., J. Fluids Struct., 120, 103899, 2023.<br>- R. Lagrange et al., J. Appl. Mech., 90(6) : 061009, 2023.</td></tr><tr><th scope="row">Actions</th><td>- Perform Helmholtz expansion and modal projection of the dynamical equations<br>- Implement coupling models for beams in TrioCFD<br>- Carry out parametric studies to determine new scaling laws for the frequency and damping characteristics of the free oscillations<br>- Compare TrioCFD simulations to theoretical predictions</td></tr></tbody></table>

<footnote n="3">³ Une collaboration (thèse de M. Kocher, 2021-2024) entre le laboratoire DYN (R. Lagrange et D. Panunzio, encadrants C.E.A.) et E.D.F. (Aurélien Joly, encadrant, et Pierre Moussou, directeur) est en cours afin de répondre aux objectifs conjoints du C.E.A. et d’E.D.F. sur le sujet des vibrations de circuit primaire principal. Le travail de thèse consiste à construire, avec les outils habituels (simulations numériques, essais), un modèle acoustique-mécanique représentatif de la cuve et de ses internes (modèle de transmission des ondes acoustiques entre les entrées et sorties de la cuve, via la prise en compte de la turbulence dans l’espace annulaire et des oscillations de l’enveloppe de cœur).</footnote>

<footnote n="4">⁴ Le réacteur Jules Horowitz est un réacteur d’essai principalement dédié aux études sur les comportements des matériaux et des combustibles utilisés dans le domaine nucléaire. Les tubes guides sont des structures cylindriques, en aluminium, qui assurent un positionnement précis des éléments à tester, permettant de contrôler et d’étudier leurs interactions avec les éléments combustibles. Les tubes guides sont insérés dans l’espace libre cylindrique entre les assemblages combustibles. Comme les éléments combustibles, ces tubes sont soumis à un écoulement axial monophasique provenant du circuit primaire, avec une vitesse caractéristique de 15 m/s dans la partie centrale.</footnote>
<table><tr><td>Expected</td><td>- 2024</td></tr><tr><td colspan="2">Abrupt change of section</td></tr><tr><td>Goal</td><td>- Simulate and model the fluid forces acting on a rigid slender structure confined in an environment with an abrupt change of section</td></tr><tr><td>References</td><td>- R. Lagrange et al., J. Appl. Mech., 90(6): 061009, 2023.<br>- Hydro note by B. Molin (inviscid)</td></tr><tr><td>Actions</td><td>- Perform Helmholtz expansion with matching conditions<br>- Perform simulations with the A.L.E. module and with the beam model (internal coupling) in TrioCFD<br>- Carry out parametric studies to determine the effect of a change of section on the added-coefficients<br>- Compare TrioCFD simulations to theoretical predictions</td></tr><tr><td>Expected</td><td>- 2024</td></tr><tr><td colspan="2">Stability under axial flow</td></tr><tr><td>Goal</td><td>- Simulate and model the coupled dynamics of slender structures subject to axial flow</td></tr><tr><td>References</td><td>- Internship E. Diaferia (2021)<br>- R. Lagrange et al., J. Appl. Mech., 89(8): 081006, 2022.<br>- R. Lagrange et al., J. Fluids Struct., 114, 103754, 2022.<br>- R. Lagrange et al., J. Fluids Struct., 120, 103899, 2023.<br>- R. Lagrange et al., J. Appl. Mech., 90(6): 061009, 2023.</td></tr><tr><td>Actions</td><td>- Couple the fluid equations (potential theory) with the Euler-Bernoulli beam equation<br>- Perform simulations with the A.L.E. module and with the beam model (internal coupling) in TrioCFD<br>- Carry out parametric studies to establish new Argand diagrams<br>- Compare TrioCFD simulations to theoretical predictions</td></tr><tr><td>Expected</td><td>- 2025</td></tr></table>
<figure><img src="image_34.png" /><figcaption>FIGURE 6.1 : Couplage hydrodynamique entre deux cylindres concentriques flexibles séparés par un fluide Newtonien au repos. (a) Schéma de principe de deux cylindres en vibrations imposées. (b) Exemple de champ adimensionné de pression généré par le déplacement rigide du cylindre intérieur, simulation numérique TrioCFD (figure de gauche), théorie potentielle 2D (figure de droite). Vidéo de l'I.F.S. entre deux cylindres concentriques (simulation TrioCFD) : ici. (c) Evolution des coefficients modaux adimensionnés ajoutés (masse et amortissement en fluide au repos) en fonction de l'élancement des cylindres. Le cylindre intérieur vibre selon le premier mode de flexion d'une poutre encastrée-libre, théorie visqueuse 3D.</figcaption></figure>
<table><caption>TABLE 6.2 : Programme de recherche sur l'I.F.S. en écoulement transverse monophasique ou diphasique.</caption><thead><tr><th colspan="2">Hydrodynamic interaction between two parallel cylinders (single phase)</th></tr></thead><tbody><tr><th colspan="2">Frequency effect</th></tr><tr><th>Goal</th><td>- Simulate and model the Stokes effect, i.e. frequency and viscosity, on the forced oscillations of two slender structures immersed in still fluid</td></tr><tr><th>References</th><td>- Internships M. Monteverde and D. Zegarra (2020)<br>- R. Lagrange et al., J. Fluids Struct., 77, 102-114, 2018.<br>- R. Lagrange et al., J. Fluids Struct., 92, 102818, 2020.</td></tr><tr><th>Actions</th><td>- Perform Helmholtz expansion with 2D bipolar coord.<br>- Perform simulations with the A.L.E. module of TrioCFD<br>- Carry out parametric studies to determine new scaling laws for the added-coefficients<br>- Compare TrioCFD simulations to theoretical predictions</td></tr><tr><th>Expected</th><td>- 2025</td></tr><tr><th colspan="2">Three-dimensional and stiffness effects</th></tr><tr><th>Goal</th><td>- Simulate and model the 3D and stiffness effects on the forced oscillations of two slender structures immersed in still fluid</td></tr><tr><th>References</th><td>- Internships M. Monteverde and D. Zegarra (2020)<br>- R. Lagrange et al., J. Fluids Struct., 77, 102-114, 2018.<br>- R. Lagrange et al., J. Fluids Struct., 92, 102818, 2020.</td></tr><tr><th>Actions</th><td>- Perform Helmholtz expansion with 3D bipolar coord.<br>- Perform simulations with the A.L.E. module of TrioCFD<br>- Carry out parametric studies to determine new scaling laws for the added-coefficients<br>- Compare TrioCFD simulations to theoretical predictions</td></tr><tr><th>Expected</th><td>- 2025</td></tr><tr><th colspan="2">Free vibrations</th></tr><tr><th>Goal</th><td>- Simulate and model the free oscillations of a slender structure vibrating in a still fluid</td></tr><tr><th>References</th><td>- Internships M. Monteverde and D. Zegarra (2020)<br>- R. Lagrange et al., J. Fluids Struct., 77, 102-114, 2018.<br>- R. Lagrange et al., J. Fluids Struct., 92, 102818, 2020.</td></tr><tr><th>Actions</th><td>- Couple the fluid equations (potential theory) with the Euler-Bernoulli beam equation<br>- Perform simulations with the beam model of TrioCFD<br>- Carry out parametric studies to determine new scaling laws for the frequency and damping characteristics of free oscillations<br>- Compare TrioCFD simulations to theoretical predictions</td></tr><tr><th>Expected</th><td>- 2026</td></tr><tr><th colspan="2">Stability under cross-flow</th></tr><tr><th>Goal</th><td>- Simulate and model the coupled dynamics of slender structures subject to cross-flow</td></tr><tr><th>References</th><td>- Internships M. Monteverde and D. Zegarra (2020)</td></tr></tbody></table>
<table><tr><td></td><td>- Internship E. Beaugendre (2021)<br/>- R. Lagrange et al., J. Fluids Struct., 77, 102-114, 2018.<br/>- R. Lagrange et al., J. Fluids Struct., 92, 102818, 2020.</td></tr><tr><td>Actions</td><td>- Carry out potential flow theory with incident flow<br/>- Implement new turbulence models in TrioCFD<br/>- Carry out parametric studies to establish parametric Argand diagram<br/>- Compare TrioCFD simulations to theoretical predictions</td></tr><tr><td>Expected</td><td>- 2026</td></tr><tr><td colspan="2">Nonlinear vibrations (large amplitudes)</td></tr><tr><td>Goal</td><td>- Study the effect of large Keulegan-Carpenter numbers on the fluid forces and derive new power laws</td></tr><tr><td>References</td><td>- M. A. Puscas et al.,<br/>Eur. J. Mech. B Fluids, 101, 106-117, 2023.</td></tr><tr><td>Actions</td><td>- Perform asymptotic development of the Navier-Stokes equations with the Keulegan-Carpenter number<br/>- Extend the amplitudes of the vibrations possible with the A.L.E. module of TrioCFD<br/>- Carry out parametric studies to determine new scaling laws for the fluid-elastic forces<br/>- Compare TrioCFD simulations to theoretical predictions</td></tr><tr><td>Expected</td><td>- 2027</td></tr><tr><td colspan="2">Hydrodynamic interactions in tube bundles (single and two-phase)</td></tr><tr><td>Goal</td><td>- Derive new scaling laws for the turbulent fluid forces and the fluid-elastic forces</td></tr><tr><td>References</td><td>- Technical notes on DIVA and TITAN<br/>- C. Bazin et al., submitted to Int. J. Multiph. Flow, 2023.<br/>- R. Lagrange et al., submitted to J. Fluids Struct., 2023.<br/>- P. Piteau et al., submitted to J. Fluids Struct., 2023.<br/>- Internships D. Panunzio (2020), E. Beaugendre (2021)<br/>- PhD, C. Bazin (2020-2023)</td></tr><tr><td>Actions</td><td>- Extend the theoretical studies performed with two parallel cylinders to the case of multiple cylinders, using complex analysis, symmetries and conformal mappings<br/>- Develop new law of walls in TrioCFD<br/>- Develop the all regimes model in Neptune-CFD<br/>- Carry out parametric studies to determine new scaling laws for the fluid forces (Stokes, Reynolds, Scruton effects)<br/>- Compare theoretical, numerical and experimental results<br/>- Extend the results to steam-water flows and to the industrial scale</td></tr><tr><td>Expected</td><td>- 2028</td></tr></table>
<figure><img src="image_35.png" /><figcaption>FIGURE 6.2 : Couplage hydrodynamique entre deux cylindres parallèles rigides immersés dans un fluide au repos. (a) Transformation conforme des coordonnées cartésiennes (x, y) en coordonnées bipolaires (σ, τ). (b) Exemple de champ de vorticité axiale d'un mode d'Helmholtz, théorie visqueuse 2D. (c) Evolution des coefficients adimensionnés de masse ajoutée (traits continus : $m_{11}$, traits discontinus : $m_{12}$) en fonction de la distance adimensionnée $\varepsilon$ entre les deux cylindres et de leur rapport d'aspect $r$, théorie potentielle 2D. Vidéo de l'I.F.S. entre deux cylindres parallèles, théorie potentielle 2D : ici.</figcaption></figure>
ANNEXES
A. FLAMBAGE D'UNE POUTRE
IMPARFAITE SUR FONDATION
NON-LINÉAIRE

Dans cette annexe, la publication parue dans la revue *International Journal of Mechanical Sciences*, sous le titre "Solution methods for the growth of a repeating imperfection in the line of a strut on a nonlinear foundation" est reproduite.
<figure><img src="image_36.png" /></figure>

ELSEVIER

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

# Solution methods for the growth of a repeating imperfection in the line of a strut on a nonlinear foundation

R. Lagrange*, D. Averbuch**

IFP Énergies nouvelles, Rond-point de l'échangeur de Solaize, BP3, 69360 Solaize, France

## ARTICLE INFO

**Article history:**

Received 1 June 2011

Received in revised form

9 June 2012

Accepted 23 June 2012

Available online 29 June 2012

**Keywords:**

Buckling

Nonlinear elastic foundation

Imperfection

Stability and bifurcation

## ABSTRACT

This paper is a theoretical and numerical study of the uniform growth of a repeating sinusoidal imperfection in the line of a strut on a nonlinear elastic Winkler-type foundation. The imperfection is introduced by considering an initially deformed shape which is a sine function with an half wavelength. The restoring force is either a bi-linear or an exponential profile. Periodic solutions of the equilibrium problem are found using three different approaches: a semi-analytical method, an explicit solution of a Galerkin method and a direct numerical resolution. These methods are found in very good agreement and show the existence of a maximum imperfection size which leads to a limit point in the equilibrium curve of the system. The existence of this limit point is very important since it governs the appearance of localization phenomena.

Using the Galerkin method, we then establish an exact formula for this maximum imperfection size and we show that it does not depend on the choice of the restoring force. We also show that this method provides a better estimate with respect to previous publications. The decrease in the maximum compressive force supported by the beam as a function of the imperfection magnitude is also determined. We show that the leading term of the development has a different exponent than in subcritical buckling of elastic systems, and that the exponent values depend on the choice of the restoring force.

© 2012 Elsevier Ltd. All rights reserved.

## 1. Introduction

The subject of beam buckling can be found in several situations in industrial applications. Among the most studied are the thermal track buckling and the buckling of subsea pipelines under the effect of temperature and/or pressure. In the latter case, buckling can appear in both the vertical and horizontal planes, according to the existing restraints imposed by the environment and the backfill. Numerous authors have studied these two applications due to their high practical importance, and have proposed solution methods to determine buckling loads and post-buckling situations. Along time, the techniques used have progressed, based firstly on analytical analyses and latter on numerical methods mostly derived from finite elements models. Thus, based on some early work by Kerr [9,10] studying the stability of railway tracks subjected to thermal buckling, several authors such as Bournazel [1], Hobbs [4,5] have proposed solution methods where the equilibrium equations were solved in

post-buckling configurations to establish relevant buckling loads. In these works, the soil was supposed to be rigid, while the external forces acting on the beam were assumed as constant as a dead weight or constant friction force. One of the key features of these theories is the fact that the loss of contact (or movement) induces a loss of global stiffness of the structure which leads to subcritical buckling and infinite buckling loads if no imperfection is assumed. Using slightly different arguments, other models were proposed in Croll [2] and Maltby and Calladine [14]. In the latter work, equilibrium equations were obtained by assuming sine deflections in the post-buckling situations and using an approximate Galerkin solution method. This method was compared with numerical solutions in Maltby and Calladine [13] and against experimental results in Maltby and Calladine [14]. Though using an approximate solution method, the approach showed good results with respect to the numerical simulation. In order to improve the earlier methods, numerical models were developed in Ju and Kyriakides [8], Klever et al. [11], Leroy and Putot [12], Yun and Kyriakides [25] to incorporate for instance additional nonlinear effects in the models, such as non linear geometric and material models. One of the key aspects of work related to the study of the upheaval buckling is the study of the localization phenomenon, which was suspected early in the pioneering work of Tveergard and Needleman [20,21]. This aspect was analyzed

<footnote>* Principal corresponding author.</footnote>

<footnote>** Corresponding author. Tel.: +33 437 702 000.</footnote>

<footnote>E-mail addresses: romain.g.lagrange@gmail.com (R. Lagrange),</footnote>

<footnote>daniel.averbuch@ifpen.fr (D. Averbuch).</footnote>

<footnote>URL: http://rlagrange.perso.centrale-marseille.fr/visible/Site/ (R. Lagrange).</footnote>
through numericalsimulations in these first papers and in Hunt and Wadee [7], Hunt and Blackmore [6], and through analytical approaches based on a double-scale expansion of the equilibrium equations in Potier-Ferry [16], for a beam resting on an elastic non-linear foundation. Based on these results, the study of the localization phenomenon was continued by using Galerkin techniques (see [22,24]) using the displacement envelopes obtained through the double-scale expansion of Potier-Ferry [16]. A lot of attention has also been paid to the estimation of the mechanical restraint induced by the soil friction and to the effect of backfill on the pipeline behavior, since the corresponding forces were found to highly influence the mechanical behavior of the pipe. In order to feed the corresponding models, experiments were performed (see [15,17–19]) either through centrifuge testing of small-scale pipeline models or through direct testing of buried full-scale pipe sections.

The present paper is an attempt to provide additional solutions for the study of the growth of a repeating imperfection in the line of a strut on a nonlinear foundation. In this work, the foundation is supposed to act through an either bi-linear or exponential regularized friction model relating the interaction line force to the transverse displacement (see Section 2). These two kinds of models are indeed found in the above mentioned papers describing the soil-pipe interaction models. In the former case, a solution method (piecewise solution) in Section 3 is proposed by explicitly solving the equilibrium equation in the regions where the foundation acts linearly and where the friction force is constant, and by connecting the two solutions by adequate boundary solutions. Alternatively, a Galerkin approach of the same problem is developed in Section 3 and leads to an explicit solution of the problem, which is developed for the two regularization models. The piecewise solution and the Galerkin approach are consequently compared together and with numerical solutions of the problem in Section 4. The post-buckling problem is then studied through the Galerkin approach which provides precise analytical solutions, focusing on the characteristics of a limit point (see Section 5) in the equilibrium curve depending on the magnitude of the initial imperfection.

## 2. Formulation of the differential equation

This section formulates the differential equation for the growth of an imperfection in the line of a strut on a nonlinear Winkler-type foundation, see Fig. 1. The imperfection is introduced by supposing an initially deformed shape $W_0$ whose form is

$$W_0 = A_0 \sin\left(\frac{\pi}{L}x\right), \qquad (1)$$

with $A_0$ the amplitude of the imperfection, $L$ its length and $x$ the longitudinal coordinate. The compressive load is $P$ and the restoring force per unit length is $\bar{P}$. These two forces are assumed to be conservative. The differential equation governing the

<figure><img src="image_38.png" /><figcaption>Fig. 1. Half-wavelength imperfection in the line of a strut resting on an elastic foundation. The imperfection is $W_0$, the compressive load is $P$ and the buckling displacement is $W$. Springs have a nonlinear elastic force-displacement relationship $\bar{P}(W)$.</figcaption></figure>

deflection $W$ may be derived either: directly by equilibration of forces; or from the Principal of Virtual Work; or using an energy formulation. The latter approach is adopted here; the total potential energy at first order being

$$V = \int_{0}^{L} \left[ \frac{1}{2} EI W'^2 - P \left( \frac{1}{2} W'^2 + W_0' W \right) - \int_{0}^{W} \bar{P}(t) dt \right] dX, \qquad (2)$$

where a prime indicates differentiation with respect to $X$. The first term is the strain energy of bending ($EI$ is the bending stiffness of the strut), the second is the work done by the load $P$ and the remainder is the energy stored in the elastic foundation. Equilibrium is given by stationary values of $V$. In what follows the strut is assumed to be simply supported, such that the conditions at the boundaries are $W(0) = W(L) = 0$. The calculus of variations on (2) gives, for a virtual displacement $\delta W$ such that $\delta W(0) = \delta W(L) = 0$

$$\delta V = EI[\delta W'W'']_0^L + \int_0^L [EIW'''' + P(W'' + W_0'') - \bar{P}(W)]\delta W \, dX, \qquad (3)$$

so that the Euler-Lagrange equation and the conditions at the boundaries for a simply supported strut are

$$EIW'''' + PW'' - \bar{P}(W) = -PW_0'', \qquad (4a)$$

$$W(0) = 0, \qquad (4b)$$

$$W''(0) = 0, \qquad (4c)$$

$$W(L) = 0, \qquad (4d)$$

$$W''(L) = 0. \qquad (4e)$$

In Eq. (4a), $P$ is the compressive load before buckling. The compressive load after buckling, considering a strut of section $S$, should be written as $N = P - (ES/2L) \int_0^L (W_x)^2 dX$, last term of this expression being a geometric shortening which allows for the additional length introduced by the lateral movement. Therefore, $N$ should be used in the equation for equilibrium. However, Tveergard and Needleman [21] have shown that the buckle will only become unstable if $N(y)$ has a maximum is correct for an isolated half-wave but is not correct for a long strut which contains a sequence of half-waves end-to-end. In such a case the key point is that a localization of the buckling, in which one particular half-wave grows at the expenses of its neighbors, can occur whenever the curve $P(y)$ has a maximum. Under this consideration, and as Maltby and Calladine [14] did, we use $P$ instead of $N$ as the load parameter.

### 2.1. The restoring force

The restoring force per unit length is assumed to be nonlinear and two particular $\bar{P}$ functions are considered. The first one is referred as the bi-linear function and is defined by

$$\bar{P}(W) = \begin{cases} -KW & \text{if } |W| < \Delta, \\ -K\Delta & \text{if } W > \Delta, \\ K\Delta & \text{if } W < -\Delta, \end{cases} \qquad (5)$$

where $K$ is the linear stiffness and $\Delta$ the mobilization. These two constants are positive. The second $\bar{P}$ function considered in this paper is referred as the exponential profile and is defined by

$$\bar{P}(W) = \begin{cases} -K\Delta(1 - e^{-W/\Delta}) & \text{if } W > 0, \\ K\Delta(1 - e^{W/\Delta}) & \text{if } W < 0. \end{cases} \qquad (6)$$

The two $\bar{P}$ functions share the same initial slope and limiting force (see Fig. 2(a)).
<figure><img src="image_39.png" /><figcaption>**Fig. 2.** (a) Restoring force. (b) Dimensionless restoring force. Bi-linear profile (solid line), exponential profile (dashed line). Dotted lines: limiting plateau (horizontal line), mobilization (vertical line).</figcaption></figure>

## 2.2. Nondimensionalization

Let us introduce a characteristic length $L_{char} = (EI/K)^{1/4}$ and nondimensional quantities

$$l = \frac{L}{L_{char}}, \quad x = \frac{X}{L_{char}}, \quad w = \frac{W}{\Delta}, \quad w_0 = \frac{W_0}{\Delta}, \quad \lambda = \frac{P}{KL_{char}^2}, \quad p = \frac{\bar{P}}{K\Delta}. \qquad (7)$$

Hence, from (4) the deflection $w$ is the solution of the differential problem

$$w''' + \lambda w'' - p(w) = -\lambda w_0'' \qquad (8a)$$

$$w(0) = 0, \qquad (8b)$$

$$w''(0) = 0, \qquad (8c)$$

$$w(l) = 0, \qquad (8d)$$

$$w''(l) = 0, \qquad (8e)$$

the dimensionless imperfection being

$$w_0 = a_0 \sin\left(\frac{\pi}{l}x\right), \qquad (9)$$

with $a_0 = A_0/\Delta$.

After nondimensionalization, the bi-linear function rewrites

$$p(w) = \begin{cases} -w & \text{if } |w| < 1, \\ -1 & \text{if } w > 1, \\ 1 & \text{if } w < -1, \end{cases} \qquad (10)$$

such that the dimensionless mobilization and limiting force equal 1. The nondimensional exponential profile is given by

$$p(w) = \begin{cases} -(1-e^{-w}) & \text{if } w > 0, \\ 1-e^{w} & \text{if } w < 0. \end{cases} \qquad (11)$$

Fig. 2(b) shows the evolution of these two dimensionless restoring forces for $w > 0$.

## 3. Theoretical resolution

In this paper, imperfections with dimensionless lengths $l < \sqrt{2}\pi$ are considered. In order to give a practical meaning to this inequation, let us consider the case of an imperfect railway. The most common rails in France (rail type 50E6) have a weight per meter of about $m=500 \text{ N m}^{-1}$, a length from 10 to 400 m

(length between two joints connection) and a bending stiffness of about $EI = 4 \times 10^6 \text{ N m}^2$. Considering a coefficient of friction $\phi$ between steel and concrete of 0.4 and a mobilization $\Delta$ from 0.1 to 1 mm it comes a linear stiffness $K = mg\phi/\Delta$ (g being the gravity) from $2 \times 10^5$ to $2 \times 10^6 \text{ N m}^{-2}$ and a characteristic length $L_{char}$ from 1 to 2 m. Therefore $l < \sqrt{2}\pi$ would correspond to a repeating sinusoidal imperfection in the railway whose half-wavelength is no more than 10 m. In a more general way, $l < \sqrt{2}\pi$ deals with imperfections whose length does not exceed some meters. Therefore the specific calculations described in this paper have been made in the context of a small-scale experimental setup.

For $l < \sqrt{2}\pi$, the first buckling mode predicted by the linear analysis is excited when

$$\lambda = \lambda_c = \left(\frac{\pi}{l}\right)^2 + \left(\frac{\pi}{l}\right)^{-2}, \qquad (12)$$

and it has the same shape as the imperfection. In what follows we will also introduce the Euler load $\lambda_e$

$$\lambda_e = \left(\frac{\pi}{l}\right)^2, \qquad (13)$$

which is the buckling load of the first mode when the restoring force equals 0.

Two theories are developed to solve the equilibrium problem. The first one, named piecewise solution theory is an exact resolution of the equilibrium problem when the bi-linear restoring force is considered. The second theory is based on a Galerkin method: it leads to an approximate resolution of the equilibrium problem by considering equally the bi-linear or the exponential restoring force. To initiate this method, the deflection shape is assumed to be a sinusoid, as the imperfection. Explicit solutions of the Galerkin equation are obtained without any assumptions.

### 3.1. Piecewise solution theory

The principle of the piecewise solution theory is to solve (8a) on each piece of the bi-linear function. Then, the solutions are connected thanks to the boundary conditions and assuming the continuity of $w, w', w''$ and $w'''$ at two connecting points $x_1$ and $x_2$.

Substituting $p(w)$ by $-w$ and by $-1$ in (8a) yields respectively

$$w_1''' + \lambda w_1'' + w_1 = -\lambda w_0'', \qquad (14a)$$

$$w_2''' + \lambda w_2'' = -1 - \lambda w_0'', \qquad (14b)$$
The solutions $w_1$ and $w_2$ belong to two affine spaces of dimension 4 given by

$$w_1 = w_{1h} + w_{1part}, \qquad (15a)$$

$$w_2 = w_{2h} + w_{2part}, \qquad (15b)$$

where $w_{1h}$ and $w_{2h}$ satisfy the homogeneous equations

$$w_1''' + \lambda w_1'' + w_{1h} = 0, \qquad (16a)$$

$$w_2''' + \lambda w_2'' + w_{2h} = 0. \qquad (16b)$$

Inserting $w_{1h} = Ae^{r_1 x}$ and $w_{2h} = Be^{r_2 x}$ in (16) yields two algebraic equations

$$r_1^4 + \lambda r_1^2 + 1 = 0, \qquad (17a)$$

$$r_2^4 + \lambda r_2^2 = 0, \qquad (17b)$$

whose solutions are

$$r_1^{(1)} = \frac{1}{2}[-2\lambda + 2(\lambda^2 - 4)^{1/2}]^{1/2}, \qquad (18a)$$

$$r_1^{(2)} = -\frac{1}{2}[-2\lambda + 2(\lambda^2 - 4)^{1/2}]^{1/2}, \qquad (18b)$$

$$r_1^{(3)} = \frac{1}{2}[-2\lambda - 2(\lambda^2 - 4)^{1/2}]^{1/2}, \qquad (18c)$$

$$r_1^{(4)} = -\frac{1}{2}[-2\lambda - 2(\lambda^2 - 4)^{1/2}]^{1/2}, \qquad (18d)$$

and

$$r_2 = 0, \quad \pm \sqrt{\lambda} i. \qquad (19)$$

Let us introduce $\alpha_i$ and $\omega_i$ the real and imaginary parts of $r_1^{(i)}$. For $\lambda < 2$, the roots $r_1$ are complex numbers and the function $w_{1h}$ is

$$w_{1h} = A_1 e^{\alpha_1 x} \cos(\omega_1 x) + A_2 e^{\alpha_2 x} \cos(\omega_2 x) + A_3 e^{\alpha_3 x} \sin(\omega_1 x) + A_4 e^{\alpha_4 x} \sin(\omega_2 x), \qquad (20)$$

with $A_1, A_2, A_3$ and $A_4$ four real constants. For $\lambda = 2$, the roots $r_1$ are double imaginary numbers ($r_1^{(1)} = r_1^{(3)}$ and $r_1^{(2)} = r_1^{(4)}$). The function $w_{1h}$ is

$$w_{1h} = (A_1 x + A_2) \cos(x) + (A_3 x + A_4) \sin(x). \qquad (21)$$

For $\lambda > 2$, the roots $r_1$ are imaginary numbers. The solution $w_{1h}$ is

$$w_{1h} = A_1 \cos(\omega_1 x) + A_2 \cos(\omega_2 x) + A_3 \sin(\omega_1 x) + A_4 \sin(\omega_2 x). \qquad (22)$$

For $\lambda \neq 0$ the function $w_{2h}$ is given by

$$w_{2h} = B_1 + B_2 x + B_3 \cos(\sqrt{\lambda} x) + B_4 \sin(\sqrt{\lambda} x), \qquad (23)$$

with $B_1, B_2, B_3$ and $B_4$ four real constants.

Since $w_{1h}$ and $w_{2h}$ depend on four undetermined constants, these functions will be noted $w_{1h}(x, A_1, A_2, A_3, A_4)$ and $w_{2h}(x, B_1, B_2, B_3, B_4)$.

The functions $w_{1part}$ and $w_{2part}$ appearing in (15a) and (15b) are particular solutions of (14a) and (14b), with $w_0$ given by (9). Searching $w_{1part}$ and $w_{2part}$ as $w_{1part} = Aw_0$ and $w_{2part} = Cx^2 + Aw_0$ yields

$$w_{1part} = \frac{\lambda}{\lambda_c - \lambda} w_0, \qquad (24a)$$

$$w_{2part} = -\frac{1}{\lambda} \frac{x^2}{2} + \frac{\lambda}{\lambda_c - \lambda} w_0. \qquad (24b)$$

Finally, the functions $w_1$ and $w_2$, defined by (15a) and (15b), are

$$w_1 = w_{1h}(x, A_1, A_2, A_3, A_4) + \frac{\lambda}{\lambda_c - \lambda} w_0, \qquad (25a)$$

$$w_2 = w_{2h}(x, B_1, B_2, B_3, B_4) - \frac{1}{\lambda} \frac{x^2}{2} + \frac{\lambda}{\lambda_c - \lambda} w_0. \qquad (25b)$$

In what follows, we will introduce the function $w_3$ defined as

$$w_3 = w_{1h}(x, C_1, C_2, C_3, C_4) + \frac{\lambda}{\lambda_c - \lambda} w_0. \qquad (26)$$

This function is a solution of (14a), as $w_1$, but with four different undetermined constants $C_1, C_2, C_3$ and $C_4$.

The aim of the piecewise solution theory is to search for a deflection $w$ solution of (8) by connecting the functions $w_1, w_2$ and $w_3$ at two unknown points $x_1$ and $x_2$ such that

$$w = \begin{cases} w_1 & \text{if } x \in [0, x_1[, \\ 1 & \text{if } x = x_1, \\ w_2 & \text{if } x \in ]x_1, x_2[, \\ 1 & \text{if } x = x_2, \\ w_3 & \text{if } x \in ]x_2, l[. \end{cases} \qquad (27)$$

With these notations, the boundary conditions (8b), (8c), (8d) and (8e) rewrite as respectively

$$w_1(0) = 0, \qquad (28a)$$

$$w_1''(0) = 0, \qquad (28b)$$

$$w_3(l) = 0, \qquad (28c)$$

$$w_3''(l) = 0. \qquad (28d)$$

The continuity of the displacement, the tangent, the curvature and the shear at $x_1$ yields

$$w_1(x_1) - w_2(x_1) = 0, \qquad (29a)$$

$$w_1'(x_1) - w_2'(x_1) = 0, \qquad (29b)$$

$$w_1''(x_1) - w_2''(x_1) = 0, \qquad (29c)$$

$$w_1'''(x_1) - w_2'''(x_1) = 0, \qquad (29d)$$

and at $x_2$

$$w_2(x_2) - w_3(x_2) = 0, \qquad (30a)$$

$$w_2'(x_2) - w_3'(x_2) = 0, \qquad (30b)$$

$$w_2''(x_2) - w_3''(x_2) = 0, \qquad (30c)$$

$$w_2'''(x_2) - w_3'''(x_2) = 0. \qquad (30d)$$

Eqs. (28), (29) and (30) lead to a linear system with 12 equations and 12 unknowns (i.e. the amplitudes $A_i, B_i$ and $C_i$). A matrix representation of this system is

$$G(x_1, x_2) = b(w_0, x_1, x_2), \qquad (31)$$

with $G$ a 12 by 12 real matrix and $b$ the vector of the unknown amplitudes. The vector $b$ contains the particular solutions $w_{1part}$ and $w_{2part}$ which depend on $w_0, x_1$ and $x_2$.

If $\det(G) \neq 0$ then

$$a(w_0, x_1, x_2) = G^{-1}(x_1, x_2)b(w_0, x_1, x_2). \qquad (32)$$

Eq. (32) expresses the amplitudes as functions of $x_1$ and $x_2$. These two connecting points are obtained by solving the nonlinear system

$$f_1(w_0, x_1, x_2) = w_1(x_1) - 1 = 0, \qquad (33a)$$

$$f_2(w_0, x_1, x_2) = w_2(x_2) - 1 = 0. \qquad (33b)$$

The numerical resolution of this system has been carried out with Matlab, using the `fzero` function. This function tries to find a zero of (33) near $X_0, X_0$ being a vector of length two. Depending on $\lambda$ (33) has 0 or several solutions. For small $\lambda$, the restoring force $p$ remains linear, such that $w = w_1$ for any $x \in [0, l]$. Thus, (33) has no solution. For high $\lambda$, the restoring force $p$ is nonlinear: the existence and the uniqueness of a solution for (33) is not trivial.
<figure><img src="image_40.png" /><figcaption>**Fig. 3.** Equilibrium paths predicted by (a) the piecewise solution theory, (b) the Galerkin method, (c) the numerical resolution of (8), case of a bi-linear restoring force. On each graph, the equilibrium paths are plotted (from top to bottom) for $a_0 = 0$, $a_0 = 0.0238$, $a_0 = 0.595$, $a_0 = 1.19$ and $l=3$. Dotted lines: critical load (upper line), Euler load (lower line), mobilization (vertical line).</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>**Fig. 4.** Deflection predicted by the piecewise solution theory. (a) is a zoom of (a). Function $w_1$ (solid line), function $w_2$ (dashed line), function $w_3$ (dash-dotted line). The functions $w_1$ and $w_2$ are connected at $x = x_1 \approx 0.618$ (first vertical dotted line). The functions $w_2$ and $w_3$ are connected at $x = x_2 \approx 2.382$ (second vertical dotted line). The horizontal dotted line represents the mobilization. $l=3$, $a_0 = 0$, $\lambda = 1.75$.</figcaption></figure>

In Fig. 5 the equilibrium paths predicted for a bi-linear restoring force have been reported in order to discuss the influence of the restoring force. It appears that the restoring force (bi-linear or exponential) has no influence on the shape of the

equilibrium paths (the variations are preserved, it always exists a limit point for small imperfection sizes, the same asymptotes are recovered). Nevertheless, the equilibrium paths for an exponential profile are below the equilibrium paths for a bi-linear profile.
<figure><img src="image_42.png" /><figcaption>**Fig. 5.** Equilibrium paths predicted by the Galerkin method. Bi-linear restoring force (solid line), exponential restoring force (dashed lines). Galerkin method (upper dashed line), Galerkin method with the [14] assumption (lower dashed line). Dotted lines: critical load (upper line), Euler load (lower line), mobilization (vertical line). $l=3$, $a_0 = 0$ (a), $a_0 = 0.0238$ (b), $a_0 = 0.595$ (c), $a_0 = 1.19$ (d).</figcaption></figure>

Therefore, the choice of the restoring force has a non-negligible influence on the maximal load acceptable by the system. This point is the purpose of the next section.

## 5. Limit point

A limit point corresponds to a maximum of $\lambda$. Therefore, if it exists, this point $(y_m, \lambda_m)$ satisfies $(d\lambda/dy)(y_m) = 0$, that is to say using (37)

$$ (a_0 + y_m) \frac{dQ}{dy}(y_m) - Q(y_m) = -a_0 \lambda_c \lambda_e. \quad (41) $$

This equation yields $Q(y_m)/\lambda_e = a_0\lambda_c - (a_0 + y_m)f(a_0)$ with $f$ a function of $a_0$. Inserting this last expression in (37) yields

$$ \lambda_m = \lambda(y_m) = \lambda_c - f(a_0). \quad (42) $$

In practice, the coordinates of a limit point are obtained by varying $y_m$ and looking for an amplitude $a_0$ which satisfies (41). Fig. 7(a) represents the evolution of $y_m$ as a function of $a_0$. For the two restoring forces it appears that $y_m$ is an increasing function of $a_0$ and diverges for the same critical amplitude. This amplitude is determined by enforcing $y_m \to \infty$ in (41). It yields

$$ a_{0c} = \frac{4}{\pi} \lambda_e^{-2}, \quad (43) $$

whose dimensional equivalent form is

$$ A_{0c} = \frac{4 L^4}{\pi \pi^4 EI} K \Delta. \quad (44) $$

This result is of great importance since it governs the appearance of localized solutions, as demonstrated in Tveergard and Needleman [21]. In that paper, they indeed showed the close link between the existence of a limit point in the $\lambda(y)$ curve and localization, the bifurcation into a localized form occurring a little way beyond this limit point.

Considering the bi-linear restoring force, when $a_0 \to 0$, Eq. (41) has an infinity of solutions $y_m$ in $[0, 1]$. This observation is coherent with the equilibrium paths depicted in Fig. 3 since the points lying on the plateau $\lambda = \lambda_c$ are a maximum.

Considering the exponential restoring force, $y_m \to 0$ when $a_0 \to 0$. This result is coherent with the equilibrium paths depicted in Figs. 5(a) and 6(a) since $(0, \lambda_c)$ is a maximum.

Substituting $Q(y_m)$ by $N(y_m) = -y_m + 1 - e^{-y_m}$ in (41) (assumption made by [14]) and enforcing $y_m \to \infty$ yields

$$ a_{0c} = \lambda_e^{-2}, \quad (45) $$

whose dimensional equivalent form is

$$ A_{0c} = \frac{L^4}{\pi^4 EI} K \Delta, \quad (46) $$

which is $4/\pi$ times smaller than the exact prediction (43). This result means that a stable equilibrium (i.e. an equilibrium path without a limit point) predicted with the assumption introduced by Maltby and Calladine [14], could in fact be unstable (see Fig. 9).

Fig. 7(b) represents the evolution of the limit point $\lambda_m$ as a function of $a_0$. Considering the two restoring forces, the same asymptotes are observed: $\lambda_m \to \lambda_e$ when $a_0 \to a_{0c}$ and $\lambda_m \to \lambda_c$ when $a_0 \to 0$. We also observe that the assumption introduced by [14]
<figure><img src="image_43.png" /><figcaption>**Fig. 6.** Equilibrium paths obtained via a numerical resolution of (8), case of an exponential restoring force. Dotted lines: critical load (upper line), Euler load (lower line), mobilization (vertical line). $l=3$, $a_0=0$ (a), $a_0=0.0238$ (b), $a_0=0.595$ (c), $a_0=1.19$ (d).</figcaption></figure>

<figure><img src="image_44.png" /><figcaption>**Fig. 7.** Limit point as a function of the amplitude of the imperfection. Bi-linear restoring force (solid line), exponential restoring force (dashed line), exponential restoring force with the assumption introduced by [14] (dash-dotted line). Dotted lines: (a) critical amplitudes, (b) critical load (upper line), Euler load (lower line), critical amplitudes (vertical lines). $l=3$.</figcaption></figure>

(see Section 3.2) leads to an underestimation of the limit load. These observations are coherent with the equilibrium paths depicted in Figs. 3, 5 and 6.

Fig. 8 shows the evolution of $\lambda_c - \lambda_m$ as a function of $a_0$, in a logarithmic scale. This picture allows to determine a

scaling for the function $f$ appearing in (42). Considering the bi-linear restoring force, for small amplitudes $a_0$, $\lambda_m - \lambda_c$ scales as

$$ \lambda_m - \lambda_c \sim -a_0. \quad (47) $$
<figure><img src="image_45.png" /><figcaption>Fig. 8. Scaling of the limit point. Bi-linear restoring force (solid line), exponential restoring force (dashed line), exponential restoring force with the assumption introduced by Maltby and Calladine [14] (dash-dotted line), critical amplitudes (dotted lines). l=3.</figcaption></figure>

<figure><img src="image_46.png" /><figcaption>Fig. 9. Equilibrium paths predicted by the exact Galerkin method (solid curve) and the Galerkin method with the assumption introduced by Maltby and Calladine [14] (dashed line). The dotted line corresponds to the Euler load. a₀=0.9, l=3.</figcaption></figure>

When the exponential restoring force is considered this scaling becomes

$$ \lambda_m - \lambda_c \sim -a_0^{1/2}. \qquad (48) $$

We observe in Fig. 8 that the assumption introduced by Maltby and Calladine [14] has no influence on this last scaling.

## 6. Conclusion

In this paper, the growth of a repeating sinusoidal imperfection in the line of a strut on a nonlinear elastic Winkler-type foundation is considered. The imperfection is introduced by considering an initially sinusoidal deformed shape with an half wavelength. The imperfection length is chosen such that the buckle mode predicted by the linear theory has the same shape as the imperfection (first buckle mode). The nonlinearities are only due to the restoring force provided by the foundation. This restoring force is expressed as a force-displacement relationship which is either a bi-linear or an exponential function. The equilibrium problem is solved using three different methods. The first one, named piecewise solution theory, is dedicated to the bi-linear profile and leads to an exact resolution of the equilibrium problem. The second one is available whatever the restoring force and is based on a Galerkin procedure. This procedure is initiated with a test function which has the same shape as the imperfection. It yields an explicit relation between the compressive load and the amplitude of the test function. This expression is an exact solution of the Galerkin equation and gives an approximate solution of the equilibrium problem. The last method is a numerical resolution of the equilibrium problem, using the ODE45 solver from Matlab. These three solving methods yield the same results: whatever the restoring force (bi-linear or exponential), the bifurcation is subcritical, the system is imperfection sensitive and the deformed shape is an amplification of the default. Moreover, it exists a critical imperfection size $a_{0c} = (4/\pi)\lambda_e^{-2}$ ($\lambda_e$ being the Euler load) which does not depend on the restoring force and such that

* if $a_0 > a_{0c}$, then the equilibrium path shows the load increasing monotonically and remains asymptotic to the Euler load.

* if $a_0 < a_{0c}$, then the equilibrium path shows the load increasing at first but then hits a limit point and the rest of the path is asymptotic to the Euler load.

This paper provides a better estimate of $a_{0c}$ with respect to previous publications.

For each restoring force, an approximate mathematical rule is derived relating the imperfection size $a_0$ to the corresponding limit load $\lambda_m$. Considering the bi-linear profile (resp. the exponential profile) the limit point scales as $\lambda_m - \lambda_c \sim -a_0$ (resp. $\lambda_m - \lambda_c \sim -a_0^{1/2}$), where $\lambda_c$ is the critical load issued from the classical linear analysis. Therefore, the scaling of the limit point depends on the regularization method.

In this paper, the restoring force and the compressive load are independent. Nevertheless, in some industrial applications (such as in drilling problems) the restoring force slightly depends on the axial compressive load. Therefore, we are currently carrying out a study with a bi-linear restoring force proportional to the axial load. First results issued from the Galerkin approach indicate that it is necessary to redefine the dimensionless parameters, leading to new scalings for the critical imperfection size and the limit load.

## Appendix A. Function Q, bi-linear restoring force

The aim of this appendix is to calculate the function Q appearing in (39) when the bi-linear restoring force is considered.

Introduce H the Heaviside function defined as $H(x) = 0$ if $x < 0$ and $H(x) = 1$ if $x \ge 0$. The restoring force p rewrites

$$ p = -w(Sgn(w) - w)H(|w| - 1), \qquad (A.1) $$
where $Sgn$ is the sign function. Since $p = -w-N(w)$ it comes

$$N(w) = (Sgn(w) - w)H(|w| - 1). \qquad (A.2)$$

For an imperfection with an half-wavelength, it can be assumed that $w > 0$. Then, the $Q$ function is (see (38))

$$Q(y) = -\frac{2}{l} \int_{0}^{l} \left[ y \sin\left(\frac{\pi}{l}x\right) - 1 \right] \sin\left(\frac{\pi}{l}x\right) H\left(y \sin\left(\frac{\pi}{l}x\right) - 1\right) dx. \quad (A.3)$$

If $y < 1$ then $y \sin((\pi/l)x) - 1 < 0$ so $H(y \sin((\pi/l)x) - 1) = Q = 0$. Therefore, the function $Q$ can be written as

$$Q(y) = -\frac{2H(y-1)}{l} \int_{0}^{l} \left[ y \sin\left(\frac{\pi}{l}x\right) - 1 \right] \sin\left(\frac{\pi}{l}x\right) H\left(y \sin\left(\frac{\pi}{l}x\right) - 1\right) dx. \quad (A.4)$$

The change of variable $t = (\pi/l)x$ gives

$$Q(y) = -\frac{2}{\pi} H(y-1) \int_{0}^{\pi} [y \sin(t) - 1] \sin(t) H(y \sin(t) - 1) dt. \quad (A.5)$$

The argument of the Heaviside function under the integral sign equals 0 when $\sin(t) = 1/y$, that is to say for $t = t_1 = \arcsin(1/y)$ and $t = \pi - t_1$. It comes that the function $Q$ is non-zero between $t_1$ and $\pi - t_1$

$$Q(y) = -\frac{2}{\pi} H(y-1) \int_{t_1}^{\pi-t_1} [y \sin(t) - 1] \sin(t) dt. \quad (A.6)$$

Finally, this integral yields

$$Q(y) = \frac{2}{\pi} H(y-1) \left\{ y \left[ \arcsin\left(\frac{1}{y}\right) - \frac{\pi}{2} \right] + \left(\frac{y^2-1}{y^2}\right)^{1/2} \right\}, \quad (A.7)$$

and the result from (39) is recovered.

## Appendix B. Function $Q$, exponential restoring force

The aim of this Appendix is to calculate the function $Q$ appearing in (40) when the exponential restoring force is considered. For this calculus, we recall that the modified Struve and Bessel functions of parameter 1 can be expended as power series

$$L_1(y) = \sum_{p=1}^{+\infty} \frac{2 (p!)^2}{\pi (2p+1)[(2p)!]^2} (2y)^{2p}, \quad (B.1a)$$

$$I_1(y) = \sum_{p=0}^{+\infty} \frac{1}{p!(p+1)!} \left(\frac{y}{2}\right)^{2p+1}. \quad (B.1b)$$

Substituting the exponential function in (11) by its power series yields

$$p(w) = -w + \sum_{n=2}^{+\infty} \frac{(-1)^n w^n}{n!}, \quad (B.2)$$

so that

$$N(w) = \sum_{n=2}^{+\infty} \frac{(-1)^{n+1} w^n}{n!}. \quad (B.3)$$

Therefore, (38) gives

$$Q(y) = \frac{2}{l} \int_{0}^{l} \sum_{n=2}^{+\infty} \frac{(-1)^{n+1} \sin\left(\frac{\pi}{l}x\right)^{n+1}}{n!} y^n dx. \quad (B.4)$$

Inverting the sum and integral signs and introducing the change of variable $t = (\pi/l)x$ yields

$$Q(y) = \frac{4}{\pi} \sum_{n=2}^{+\infty} \frac{(-1)^{n+1} W_{n+1} y^n}{n!}, \quad (B.5)$$

with $W_n$ the Wallis integral

$$W_n = \int_{0}^{\pi/2} \sin^n(x) dx. \quad (B.6)$$

The terms $W_n$ are classical to calculate. For $n = 2p$ and $n = 2p+1$ it yields

$$W_{2p} = \frac{(2p)!}{2^{2p}(p!)^2} \frac{\pi}{2}, \quad (B.7a)$$

$$W_{2p+1} = \frac{2^{2p}(p!)^2}{(2p+1)!}, \quad (B.7b)$$

Splitting the series appearing in (B.5) into odd and even indices gives $Q(y) = \Sigma_1(y) + \Sigma_2(y)$ with

$$\Sigma_1(y) = -\frac{4}{\pi} \sum_{p=1}^{+\infty} \frac{W_{2p+1}}{(2p)!} y^{2p}, \quad (B.8a)$$

$$\Sigma_2(y) = \frac{4}{\pi} \sum_{p=1}^{+\infty} \frac{W_{2(p+1)}}{(2p+1)!} y^{2p+1}. \quad (B.8b)$$

Inserting (B.7a) and (B.8a) yields

$$\Sigma_1(y) = -2 \sum_{p=1}^{+\infty} \frac{2 (p!)^2}{\pi (2p+1)[(2p)!]^2} (2y)^{2p}, \quad (B.9a)$$

$$\Sigma_2(y) = 2 \sum_{p=0}^{+\infty} \frac{1}{p!(p+1)!} \left(\frac{y}{2}\right)^{2p+1} - y. \quad (B.9b)$$

Eq. (B.1) gives

$$\Sigma_1(y) = -2L_1(y), \quad (B.10a)$$

$$\Sigma_2(y) = 2I_1(y) - y, \quad (B.10b)$$

with respectively $L_1$ and $I_1$ the modified Struve and Bessel functions of parameter 1. Finally

$$Q(y) = \Sigma_1(y) + \Sigma_2(y) = 2[L_1(y) - L_1(y)] - y, \quad (B.11)$$

and the result from (40) is recovered.

## References

[1] Bournazel C. Vertical buckling of buried pipes. Revue de l'Institut Francais du Pétrole 1982;37(1):113-22.

[2] Croll JGA. A simplified model of upheaval thermal buckling of subsea pipelines. Thin Wall Struct 1997;29:59-78.

[3] Fox C. An introduction to the calculus of variations. New York: Dover; 1987.

[4] Hobbs RE. Pipeline buckling caused by axial loads. J Constr Steel Res 1981;1(2):2-10.

[5] Hobbs RE. In-service buckling of heated pipelines. J Transp Eng 1984;110:175-89.

[6] Hunt GW, Blackmore A. Principles of localized buckling for a strut on an elastoplastic foundation. J Appl Mech 1996;63(1):234-9.

[7] Hunt GW, Wadee MK. Comparative lagrangian formulations for localized buckling. Proc R Soc Lond A Math Phys Eng Sci 1991;434:485-502.

[8] Ju GT, Kyriakides S. Thermal buckling of offshore pipelines. J Offshore Mech Arctic Eng 1988;110:355-64.

[9] Kerr AD. On the stability of the railroad track in the vertical plane. Rail Int 1974;5:131-42.

[10] Kerr AD. Analysis of thermal track buckling in the lateral plane. Acta Mech 1978;30:17-50.

[11] Klever FJ, Van Helvoirt LC, Sluyterman AC. A dedicated finite-element model for analyzing buckling response of submarine pipelines. In: Proceedings from the offshore technology conference, Houston; 1990. 6333-MS.

[12] Leroy JM, Putot CJM. Behavior of buried flexible pipelines. In: Proceedings from the offshore mechanics and arctic engineering conference, Calgary; 1992.

[13] Maltby TC, Calladine CR. An investigation into upheaval buckling of buried pipelines: experimental apparatus and some observations. Int J Mech Sci 1995;37:943-63.

[14] Maltby TC, Calladine CR. An investigation into upheaval buckling of buried pipelines: theory and analysis of experimental observations. Int J Mech Sci 1995;37:965-83.

[15] Palmer AC, White DJ, Baumgard AJ, Bolton MD, Barefoot AJ, Finch M, Powell T, Faranski AS, Baldry JAS. Uplift resistance of buried submarine pipelines:
comparison between centrifuge modeling and full-scale tests. Geotechnique 2003;53(10):877-83.

[16] Potier-Ferry M. Foundations of elastic postbuckling theory. In: Lecture notes in physics, Buckling and post-buckling, vol. 288; 1987. p. 1–82.

[17] Schaminee PEL, Zorn NF, Schotman GJM. Soil response for pipeline upheaval buckling analyses: full-scale laboratory tests and modelling. In: Proceedings from the offshore technology conference, Houston; 1990. 6486-MS.

[18] Trautmann CH, O'Rourke TD, Kulhawy FH. Lateral force-displacement response of buried pipe. J Geotech Eng 1985;111(9):1077-92.

[19] Trautmann CH, O'Rourke TD, Kulhawy FH. Uplift force-displacement response of buried pipe. J Geotech Eng 1985;111(9):1061-76.

[20] Tveergard V, Needleman A. On the localization of buckling patterns. J Appl Mech 1981;47(3):613-9.

[21] Tveergard V, Needleman A. On localized thermal track buckling. Int J Mech Sci 1981;23:577-87.

[22] Wadee MA. Effects of periodic and localized imperfections on struts on nonlinear foundations and compression sandwich panels. Int J Solids Struct 2000;37:1191-209.

[23] Wadee MK, Hunt GW, Whiting AIM. Asymptotic and Rayleigh-Ritz routes to localized buckling solutions in an elastic instability problem. Proc R Soc Lond A Math Phys Eng Sci 1997;453:2085-107.

[24] Whiting AIM. A Galerkin procedure for localized buckling of a strut on a nonlinear elastic foundation. Int J Solids Struct 1997;34:727-39.

[25] Yun H, Kyriakides S. Model for beam-mode buckling of buried pipelines. J Eng Mech 1985;111(2):235-53.
B. PLISSEMENT D'UN FILM MINCE
ATTACHÉ À UN SUBSTRAT NON-PLAN

Dans cette annexe, la publication parue dans la revue *Nature Materials*, sous le titre "Curvature-induced symmetry breaking selects elastic surface patterns" est reproduite.
# Curvature-induced symmetry breaking determines elastic surface patterns

Norbert Stoop<fnref n="1" />, Romain Lagrange<fnref n="1" />, Denis Terwagne<fnref n="2" />†, Pedro M. Reis<fnref n="2" />,³ and Jörn Dunkel<fnref n="1" />*

Symmetry-breaking transitions associated with the buckling and folding of curved multilayered surfaces—which are common to a wide range of systems and processes such as embryogenesis, tissue differentiation and structure formation in heterogeneous thin films or on planetary surfaces—have been characterized experimentally. Yet owing to the nonlinearity of the underlying stretching and bending forces, the transitions cannot be reliably predicted by current theoretical models. Here, we report a generalized Swift-Hohenberg theory that describes wrinkling morphology and pattern selection in curved elastic bilayer materials. By testing the theory against experiments on spherically shaped surfaces, we find quantitative agreement with analytical predictions for the critical curves separating labyrinth, hybrid and hexagonal phases. Furthermore, a comparison to earlier experiments suggests that the theory is universally applicable to macroscopic and microscopic systems. Our approach builds on general differential-geometry principles and can thus be extended to arbitrarily shaped surfaces.

Symmetry breaking and structure formation are intrinsically linked. Symmetry-breaking transitions encompass a diverse range of phenomena, from the emergence of large-scale cosmological structures<fnref n="1" /> or the formation of sand dunes<fnref n="2" /> to crystallization of solids³ and the development of form and function in living organisms⁴,⁵. Theoretical analysis of symmetry breaking typically builds on effective nonlinear field equations that describe complex many-particle systems⁶ by a few macroscopic field variables. This approach has proved fruitful in identifying generic aspects of structure formation, as exemplified by the Ginzburg-Landau theory of phase transitions⁷ and Turing's description of reaction-diffusion patterns⁸,⁹. Generally, however, it is challenging to derive nonlinear field theories systematically¹⁰⁻¹³ from the underlying microscopic dynamics. Instead, effective field equations are often inferred from abstract symmetry considerations and bifurcation theory¹⁴, resulting in a large number of undetermined parameters that limit the predictive power and complicate comparison with experimental data. Here, we systematically derive and experimentally test an effective field theory that predicts quantitatively the surface-pattern selection in curved bilayer systems consisting of a stiff film on a soft substrate (Figs 1 and 2).

Buckling of thin films plays a prominent role in the morpho-
genesis of multilayered soft tissues, governing the wrinkling of
skin¹⁵, fingerprint formation¹⁶ and the development of brain
convolutions¹⁷. In addition to their biological relevance, wrinkling
processes under curvature constraints are attracting considerable
interest as promising techniques for nanoscale surface patterning¹⁸,
microlens array fabrication¹⁹ and adaptive aerodynamic drag
control²⁰. Recent experiments and simulations suggest that
wrinkling patterns may vary strongly with applied stress²¹⁻²⁵
and substrate curvature²⁰,²⁶⁻³⁰. However, so far, the complexity of
the numerically implemented tensor equations has prevented a
detailed analytical understanding. Despite substantial progress

in the theoretical description of planar bilayer membranes³¹,³²,
it is unclear how curvature controls pattern selection in non-
planar geometries.

The scalar field theory presented below solves this longstand-
ing problem by providing detailed quantitative predictions for
curvature- and stress-induced pattern-formation transitions. Start-
ing from Koiter's shell theory³³, we derive a generalized fourth-order
Swift-Hohenberg (GSH) equation for the normal displacement field
of a film bound to an arbitrarily curved surface (Supplementary
Information). In the case of a spherical geometry, our GSH theory
reveals that curvature triggers a transition from labyrinth-like to
hexagonal wrinkling patterns through a curvature-induced symme-
try breaking in the field equation. The theory further predicts a
coexistence region separating two 'pure' phases (Fig. 1a-c). Both
the theoretically predicted surface patterns and the analytically
predicted phase diagram agree quantitatively with data from our
macroscale experiments (Figs 1d-f and 3). The GSH model implies
that analogous transitions occur when the compressive stress in
the film is increased. This prediction is in agreement with recent
microscale experiments²⁸,³⁴ (Fig. 1g-i), suggesting that the theory
is universally applicable to both microscopic and macroscopic
systems. As our derivation of the GSH model builds on general
differential-geometric principles, it can be extended to arbitrarily
shaped surfaces, thus providing a generic framework for future
studies of curvature-controlled wrinkling in physical, biological and
chemical systems.

Theory of thin-film deformation on soft substrates

Our derivation starts from the covariant Koiter shell equations³³,
obtained from three-dimensional elasticity theory through an
expansion in the film thickness *h* → 0 (Fig. 2). Koiter's model
expresses the elastic energy of a freestanding curved shell in terms
of deformations of its central surface (Supplementary Information).
Although the Koiter equations have been successfully used in

<footnote n="1">¹Department of Mathematics, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, Massachusetts 02139-4307, USA.</footnote>

<footnote n="2">²Department of Civil & Environmental Engineering, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, Massachusetts 02139-4307, USA. ³Department of Mechanical Engineering, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, Massachusetts 02139-1713, USA. †Present address: Faculté des Sciences, Université Libre de Bruxelles (ULB), Bruxelles 1050, Belgium.</footnote>

<footnote>*e-mail: dunkel@mit.edu</footnote>
<figure><img src="image_47.png" /><figcaption>**Figure 1** | Macroscopic and microscopic wrinkling morphologies of stiff thin films on spherically curved soft substrates. a-c, Theoretical predictions based on numerical steady-state solutions of equation (1). Colour red (blue) signals inward (outward) wrinkles. Simulation parameters: (a) $γ_0 = -0.029$, $a=0.00162$, $c=0.0025$; (b) $γ_0 = -0.04$, $a=-1.26 \times 10^{-6}$, $c=0.002$; (c) $γ_0 = -0.02$, $a=1.49 \times 10^{-4}$, $c=0.0025$ (see Table 1). d-f, Experimentally observed patterns confirm the transition from hexagonal (d) to labyrinth-like wrinkles (f) via a bistable region (e) when the radius-to-thickness ratio $R/h$ (see Fig. 2) is increased. Scale bars, 10 mm. Parameters: $E_f = 2,100$ kPa, $R = 20$ mm, $v = 0.5$ and (d) $E_s = 230$ kPa, $h = 0.630$ mm; (e) $E_s = 29$ kPa, $h = 0.14$ mm; (f) $E_s = 63$ kPa, $h = 0.10$ mm. g-i, Oxide layers on microscopic PDMS hemispheres exhibit a similar transition from hexagonal to labyrinth patterns when the excess film stress is increased through changes in the ambient ethanol concentration (indicated in per cent). Scale bars, 250 µm. Micrographs courtesy of D. Breid and A. Crosby$^{28}$.</figcaption></figure>

<figure><img src="image_48.png" /><figcaption>**Figure 2** | Notation and experimental system. a, Schematic of a curved thin film adhering to a soft spherical substrate of outer radius $R$. b, The film (thickness $h$) is driven towards a wrinkling instability by the compressive film stress $\sigma$, leading to a wrinkling pattern with wavelength $\lambda$ and radial displacement $u$. c, The experimental system consists of two merged hemispherical caps. An air channel allows one to tune the film stress $\sigma$ through the pressure difference $\Delta p = p_e - p_i$.</figcaption></figure>

computational wrinkling studies$^{26,27}$, their nonlinear tensorial structure offers limited insight beyond linear stability analysis. We found, however, that substantial analytical simplifications are possible when a stiff film (Young modulus $E_f$) is adhered to a soft substrate with Young modulus $E_s \ll E_f$.

As relevant to our experiments, which are described in detail below, we consider a spherical geometry with radius $R/h \gg 1$ and

assume that film and substrate have the same Poisson ratio $v$. The generalization to non-spherical surfaces is obtained by replacing the metric tensor appropriately (Supplementary Information). Continuity across the film-substrate interface favours deformations that are dominated by the radial displacement $u$ (Fig. 2; from here onwards all lengths are normalized by $h$). Neglecting secondary lateral displacements, one can systematically expand the strain energy, which contains the original Koiter shell energy density as well as additional substrate coupling and overstress contributions, in terms of the covariant surface derivative $\nabla u$ and powers of $u$ (Supplementary Information). Functional variation of the elastic energy with respect to $u$ then yields a nonlinear partial differential equation for the wrinkled equilibrium state of the film. Assuming overdamped relaxation dynamics, one thus obtains the following GSH equation (Supplementary Information)

$$ \partial_t u = \gamma_0 \Delta u - \gamma_2 \Delta^2 u - au - bu^2 - cu^3 + \Gamma_1 [(\nabla u)^2 + 2u\Delta u] + \Gamma_2 [u(\nabla u)^2 + u^2 \Delta u] \quad (1) $$

Here, $\Delta$ denotes the Laplace-Beltrami operator, involving the surface metric tensor of the sphere and Christoffel symbols of the second kind, and $\Delta^2$ is the surface biharmonic operator$^{35}$. The ($\gamma_0, \gamma_2$) terms describe stress and bending, the (a, b, c) terms comprise local film-substrate interactions and stretching contributions, and the ($\Gamma_1, \Gamma_2$) terms account for higher-order stretching forces. For $\Gamma_1 = \Gamma_2 = 0$, equation (1) reduces to the standard Swift-Hohenberg equation, as originally derived in the context of Rayleigh-Bénard convection$^{10,36}$. The additional ($\Gamma_1, \Gamma_2$) terms will prove crucial below when matching theory and experiments. The generalization of equation (1) for arbitrary surfaces is given in Supplementary Equation (34).

The detailed derivation (Supplementary Information), combined with systematic asymptotic analysis of the planar limit $R/h \to \infty$, allows us to express the coefficients in equation (1) in terms of the standard material parameters: Poisson ratio of the film $v$,
<figure><img src="image_49.png" /><figcaption>**Figure 3** | Phase diagram of wrinkling morphologies. Experimental data points for hexagonal (blue), bistable (yellow) and labyrinth (red) patterns are shown for different values of curvature radius $R/h = \kappa^{-1}$ and excess film stress $\Sigma_e$. Symbols indicate the elastic moduli ratio $\eta = 3E_s/E_f$ (square, $\eta = 0.019$; circle, $\eta = 0.036$; clubsuit, $\eta = 0.041$; downtriangle, $\eta = 0.055$; diamond, $\eta = 0.09$; triangle, $\eta = 0.328$). The data suggest that phase boundaries are independent of $\eta$ in the experimentally tested range. Only the largest vertical error bars are shown (standard deviation of 12 amplitude measurements; see Methods). Horizontal error bars are smaller than the symbol size. Solid lines are theoretically predicted phase boundaries, obtained from equation (3) with parameter $c_1 = 0.0188$ (Table 1).</figcaption></figure>

effective curvature $\kappa = h/R$, Young ratio $\eta = 3E_s/E_f$, and excess stress $\Sigma_e = (\sigma/\sigma_c) - 1$ (Table 1). The theory contains only a single fitting parameter, $c_1$, related to the cubic stretching force term $cu^3$. Equation (1) predicts that the unbuckled solution $u=0$ is stable for negative excess stresses $\Sigma_c < 0$, whereas wrinkling occurs for $\Sigma_c \ge 0$. Linear stability analysis at $\Sigma_c = 0$ and $\kappa = 0$ reproduces the classical$^{37}$ pattern-wavelength relation for planar wrinkling, $\lambda/h = 2\pi\eta^{-1/3}$ (Supplementary Information).

Numerical simulation of equation (1) is non-trivial owing to the metric dependence of the biharmonic operator $\Delta^2$ (ref. 35). To compute the stationary wrinkling patterns (Fig. 1a-c) predicted by equation (1), we implemented a $C^1$-continuous finite-element algorithm specifically designed for covariant fourth-order problems (Methods). A main benefit of equation (1), however, is that it enables analytical prediction of the various pattern-formation regimes.

## Pattern selection

Pattern selection in the wrinkling regime $\Sigma_c \ge 0$ is a nonlinear process and, therefore, cannot be inferred from linear stability analysis. Numerical parameter scans of equation (1) yield a variety of qualitatively different stationary states that can be classified as representatives of a hexagonal phase (Fig. 1a), labyrinth phase (Fig. 1c) or intermediate coexistence phase (Fig. 1b). Qualitatively, the transition from hexagons to labyrinths can be understood through a symmetry argument: the ($b, \Gamma_1$) terms in equation (1) break the radial reflection invariance of its solutions under the transformation $u \to -u$, as also evident from the corresponding energy functional that is given in Supplementary Equation (37). As $b$ and $\Gamma_1$ are controlled by $\kappa = h/R$ (Table 1), we expect a curvature-induced symmetry-breaking transition at some critical value of $\kappa$. Furthermore, recalling that the inclusion of similar symmetry-breaking terms causes a transition from labyrinths to hexagonal

<caption>**Table 1** | List of parameters for equation (1) in units $h=1$, with $\eta = 3E_s/E_f$, $\gamma_2 = 1/12$, $\Sigma_e = (\sigma/\sigma_c) - 1$ and $\kappa = h/R$.</caption>

$$ 
\begin{aligned}
\gamma_0 &= -\frac{\eta^{2/3}}{6} - \left[ \frac{2(1+\nu)}{\eta^{2/3}} - \frac{1}{3} \right] \kappa^2 \\
a &= \frac{\eta^{4/3}}{12} + \frac{6(1+\nu) - \eta^{2/3}}{3} \kappa^2 + \bar{a}_2 \Sigma_e \\
b &= 3(1+\nu)\kappa^3 \\
c &= \frac{2(1+\nu)\eta^{2/3}}{3} c_1 \\
\Gamma_1 &= \frac{1+\nu}{2} \kappa \\
\Gamma_2 &= \frac{1+\nu}{2} \kappa^2 \\
\bar{a}_2 &= -\frac{\eta^{4/3}(c+3|\gamma_0|\Gamma_2)}{48\gamma_0^2}
\end{aligned}
$$

The only remaining fitting parameter of the model is $c_1$.

patterns in the classical Swift-Hohenberg model$^{36}$, it is plausible to expect a hexagonal phase at large curvatures $\kappa$ and labyrinths at smaller values of $\kappa$ in our system.

To obtain a quantitative prediction for the phase boundaries, we approximate equation (1) through a standard Swift-Hohenberg equation and make use of established results from nonlinear stability analysis$^{38}$. Assuming plane-wave solutions with amplitude $A$ and wavevector $k$, the $\Gamma_1$ term exerts an average force $\Gamma_1 \langle (\nabla u)^2 + 2u\Delta u \rangle_\lambda = -\Gamma_1 A^2 k^2/2$ per wavelength $\lambda$. One may therefore approximate the $\Gamma_1$ term by an effective quadratic force $-\Gamma_1 k^2 u^2$, and similarly the $\Gamma_2$ term by an effective cubic force $\Gamma_2 k^2 u^3/2$ (Supplementary Information). Inserting for $k$ the most unstable mode, $k_* = \sqrt{|\gamma_0|/(2\gamma_2)}$, equation (1) can be approximated by the standard Swift-Hohenberg equation

$$ \partial_t \phi = -2\Delta\phi - \Delta^2\phi - A\phi - B\phi^2 - \phi^3 \quad (2) $$

where $\phi = u/u_*$, $u_* = |\gamma_0|/\sqrt{(c/3) + \Gamma_2|\gamma_0|}$, $A = 3a/\gamma_0^2$, and $B = u_* [(b/3) + 2|\gamma_0|\Gamma_1]/\gamma_0^2$. Nonlinear stability analysis of equation (2) yields the critical phase transition curves as functions of $A$ and $B$ (ref. 38). Note that the coefficients in equation (2) can be directly traced back to geometric and material parameters, whereas in many other pattern formation processes Swift-Hohenberg-type equations have been applied only in a purely phenomenological manner$^6$. In terms of the original system parameters, one finds the stability criteria (Supplementary Information)

$$ 
\begin{aligned}
\text{Hexagonal phase:} \quad & -\kappa^2/(20c_1^2) < \Sigma_e < \kappa^2/c_1^2 \\
\text{Bistable phase:} \quad & \kappa^2/c_1^2 < \Sigma_e < 4\kappa^2/c_1^2 \quad (3) \\
\text{Labyrinth phase:} \quad & 4\kappa^2/c_1^2 < \Sigma_e
\end{aligned}
$$

where the parameter $c_1$ sets the strength of the cubic stretching force (Table 1). In the bistable coexistence phase, both hexagon and labyrinth solutions are stable, suggesting a strong dependence on initial conditions in this regime (Fig. 4).

Equation (3) confirms our qualitative symmetry argument and implies, moreover, that the pattern-formation transitions can be controlled not only by curvature, but also through the excess film stress $\Sigma_e$, in agreement with recent experimental results$^{28}$ (Fig. 1g-i).
<figure><img src="image_50.png" /><figcaption>**Figure 4 | Bifurcation diagram of wrinkling patterns.** Stability analysis of equation (2) predicts two hysteresis cycles, shown here for R/h = 40. Solid (dashed) lines correspond to stable (unstable) amplitude solutions; see equation (4) and Supplementary Information. **a.** The hysteresis path across the bistable phase is realized by first decreasing (red) and subsequently increasing (blue) the excess stress $\Sigma_e$. **b.** Enlarged view of the second curvature-dependent hysteresis cycle near $\Sigma_e = 0$, corresponding to the highlighted region in **a**. Starting from a stable unwrinkled solution at $\Sigma_e < 0$, the system switches to a hexagonal state at $\Sigma_e = 0$ (blue path). When decreasing the excess stress again to negative values (red), the hexagons remain stable in the subcritical region until a critical value $\Sigma_e = -\delta\Sigma_e$ is reached. The width $\delta\Sigma_e$ and height $\delta U_H$ of this hysteresis loop depend on curvature $\kappa = h/R$; see equation (4).</figcaption></figure>

## Comparison with experiments

We test the theoretical predictions, obtained from equations (1)–(3), by studying the wrinkling of centimetre-sized polydimethylsiloxane (PDMS)-coated elastomer hemispheres (Methods). In our experiments, wrinkling is controlled by the swelling of the film during fabrication and by manual depressurization after fabrication (Fig. 2a,c). The displacement field $u$, from which the excess film stress $\Sigma_e$ can be estimated through amplitude measurements³², is obtained from three-dimensional (3D) surface scans (Methods).

The experimental data confirm quantitatively the theoretically predicted curvature-induced phase transitions from hexagons to labyrinths (Figs 1 and 3). At high values of curvature $\kappa = h/R$, we find the hexagonal phase, characterized by localized spherical depressions that are typically surrounded by 6 neighbours (Fig. 1a,d), although occasional topological defects with 5 or 7 neighbours exist as required by Euler's polyhedral theorem³⁹. As predicted by equation (1), experimentally observed hexagons always buckle inwards. For intermediate values of $\kappa$, the experiments further confirm coexisting domains of hexagonal and labyrinth-like patterns (Fig. 1b,e). In our simulations of equation (1), we find that the energy of such hybrid patterns remains constant asymptotically, suggesting that they are not transient but correspond to local energy minima. When the curvature is decreased, $\kappa \to 0$, at constant stress $\Sigma_e$, the experimental system transitions into the labyrinth phase (Fig. 1c,f), characterized by a network of connected ridges and extended but disconnected valleys (Fig. 1f). Equation (1) shows that this ridge-valley asymmetry is due to the small but non-vanishing symmetry-breaking effect of curvature.

Moreover, in agreement with previous microscale experiments²⁸ (Fig. 1g-i), equations (1) and (3) imply that the phase transition from hexagons to labyrinths can also be triggered by increasing the excess film stress $\Sigma_e = (\sigma/\sigma_c) - 1$ at constant surface curvature. The morphological phase diagram constructed from our macroscale data confirms this prediction (Fig. 3). In particular, by fixing just a single fitting parameter $c_1 = 0.0188 \pm 0.0002$, the analytical results for the two critical curves in equation (3) are in good quantitative agreement with the experimental data for a wide range of Young modulus ratios $\eta = 3E_s/E_f$ (Fig. 3). Strikingly, we find that the phase boundaries are independent of $\eta$ over the range $0.019 < \eta < 0.328$ realized in our experiments, suggesting that the parameter $c_1$ may be a universal numerical constant independent of material properties.

## Predictions for future experiments

The good agreement between theory and available experimental data encourages additional predictions that ought to be tested in

future experiments. The nonlinear stability analysis of equation (2) suggests that, for sufficiently small overstress $\Sigma_e$, the hexagonal phase continues to exist even for weakly curved substrates³² with $\kappa \ll 1$ (Fig. 3). Simulations of equation (1) for time-varying overstress $\Sigma_e(t)$ confirm that, owing to the presence of symmetry-breaking terms for $\kappa \neq 0$, hexagonal patterns always appear first after crossing the wrinkling threshold $\Sigma_e = 0$ from below. Once the hexagons have been formed, they remain stable throughout the bistable phase when the film stress is slowly increased. A similar reverse effect is observed when the film stress is slowly decreased in simulations that start from the labyrinth phase. In this case, the labyrinths persist throughout the bistable region. Equation (1) makes it possible to understand such memory effects analytically (Fig. 4).

Specifically, the above bifurcation analysis of equations (1)–(3) predicts two hysteresis cycles. The first cycle relates to the onset of wrinkling at $\Sigma_e = 0$ (Fig. 4b), whereas the second encompasses the bistable phase (Fig. 4a). The amplitude $U_H = \max u_H - \min u_H$ of the hexagonal solutions $u_H$ grows according to a square-root law, shifted by the coefficient of the symmetry-breaking term in equation (2) (Supplementary Information),

$$U_H = \frac{3}{5} \left[ Bu_* + \sqrt{(Bu_*)^2 + \frac{45\Sigma_e}{4}} \right] \quad (4)$$

where $Bu_* \simeq 3\kappa/(4c_1)$ to leading order in $\kappa$, with $B$ and $u_*$ as defined in equation (2). Equation (4) implies that, for $\kappa > 0$, the hexagonal phase is stable subcritically: on reducing the excess film stress from the hexagonal phase, hexagons remain stable even when the film stress is below the critical wrinkling stress $\sigma_c$ (Fig. 4b). The width of the subcritical region, $\delta\Sigma_e = \kappa^2/(20c_1^2)$, and the amplitude at onset, $\delta U_H = U_H(\Sigma_e = 0) = 9\kappa/(10c_1)$, scale with $\kappa$. The bifurcation at $\Sigma_e = 0$ is transcritical, corresponding to a Lifshitz point³⁸. Such bifurcations are typical of Swift-Hohenberg-type models, and have been predicted and observed in optics⁴⁰ and nonlinear biological and chemical systems⁴¹. For values of $\Sigma_e$ in the subcritical hysteresis region, the hexagonal and the flat state are simultaneously stable in a narrow parameter range, potentially allowing for localized hexagonal patterns as found for the standard Swift-Hohenberg equation⁴²,⁴³.

The detailed analysis of the second hysteresis cycle (Fig. 4a) shows that the amplitude $U_L$ of the labyrinth solutions follows a square-root law (Supplementary Information). Starting from the labyrinth phase, the system remains in a labyrinth state when the film stress is lowered across the bistable region until one reaches the
instability threshold, located at $\Sigma_c \approx 1.75$ in the depicted example with $R/h = 40$ (red path in Fig. 4a). At that point, the system transitions into a hexagonal state. As $\Sigma_c$ is increased again, the film maintains the hexagonal configuration until the stress exceeds the upper instability threshold $\Sigma_c \approx 7.5$ (blue path in Fig. 4a).

The direct verification of the two predicted hysteresis cycles poses a substantial experimental challenge, requiring high accuracy in the amplitude measurements and precise reversible tuning mechanisms for the excess film stress $\Sigma_c$. For instance, the large stress variations needed to trace out the hysteresis loops with a single sample cannot be realized with the present depressurization set-up²⁰. Some preliminary experimental support for the hysteresis predictions comes from a recent study³² of low-stress films, which pointed out the frequent appearance of hexagonal patterns when the excess stress is slowly varied from negative to positive values (see also Fig. 1g-i). These findings are consistent with the results of the above bifurcation analysis (Fig. 4a). We hope that our detailed theoretical predictions will stimulate further experimental work.

In closing, we showed that a systematically derived effective field theory provides a comprehensive quantitative description of surface-pattern formation in non-planar elastic media (see Supplementary Fig. 2 for additional examples with spatially varying curvature). The observation of similar pattern transitions in systems ranging from a few micrometres²⁶,²⁸,³⁴ to several centimetres²⁰,²⁷, combined with the fact that curvature-induced pattern selection can now be understood in terms of a symmetry breaking in the effective field equations, suggests that such processes form common universality classes. The generic differential-geometric framework developed here enables a systematic classification of wrinkling phenomena in complex geometries, by examining the symmetry properties of effective higher-order differential operators built from the surface metric and film-substrate coupling forces. Moreover, equation (1) and its generalization to arbitrary shapes (Supplementary Information) provide a basis for studying weakly time-dependent phenomena such as the nucleation of wrinkling patterns under adiabatic (slow) increase of stress. Thus, in practice, the above analytical approach can help us to predict and control wrinkling processes under natural conditions, promising improved microfabrication techniques and, perhaps, even a better understanding of tissue mechanics and developmental morphogenesis.

## Methods

**Algorithm.** The fourth-order covariant derivatives in the metric-dependent biharmonic operator $\Delta^2$ make it challenging to solve equation (1) numerically. We simulate equation (1) by employing an extension of the finite-element scheme, based on subdivision surface basis functions. Previous studies⁴⁴,⁴⁵ show that this method yields high accuracy and excellent performance for related problems in nonlinear elasticity. The underlying algorithm ensures the C¹-continuity of the basis functions, as required for the numerical integration of fourth-order equations. The method also allows for a direct computation of the various covariant derivatives. As in standard finite-element algorithms, we discretize the spherical surface with a mesh consisting of up to 50,000 triangular elements. A solution coefficient $u_i$ is assigned to each of the $i=1, \dots, N$ triangle vertices, such that the system state is interpolated by $u(r) = \sum_{i=1}^{N} u_i N_i(r)$, where the $N_i$ terms are the finite-element basis functions. For each time step $[t, t+\Delta t]$, we solve the weak form of equation (1) as a sum over individual element-wise contributions. Starting from random initial conditions with $\|u_i\|^2 < 1-a$, we integrate the dynamics of the system in time using a standard explicit Euler scheme, to obtain the system state $\{u_i\}$ at time $t+\Delta t$. We determine steady-state solutions by continuously monitoring the associated free energy. We consider the system to be in a steady state if its relative change remains below $10^{-5}$ for more than 10,000 successive time steps.

**Experiments.** Hemispherical samples were fabricated using rapid digital prototyping techniques that allow flexibility in the choice of geometrical and material parameters²⁰. Samples were casted and coated using silicone-based elastomers, such as PDMS (Sylgard 184, Dow Corning), Ecoflex (Smooth-on) and vinylpolysiloxane (VPS, Zhermack), allowing us to examine a wide range of elastic moduli for film ($E_f$) and substrate ($E_s$), spanning $9 \le E_f/E_s \le 162$. Typical sample parameters are: outer radius $R=20$ mm, radius of the inner cavity 9 mm,

and film thickness $20 \mu m \le h \le 1,000 \mu m$. The inner cavity of the samples was depressurized to create a state of homogeneous compression and to trigger wrinkling of the stiff surface film. Surface profiles were measured using a NextEngine 3D Laser scanner. The excess film stress $\Sigma_c$ was estimated from the pattern amplitude $A$ determined in the 3D surface scans, using the established amplitude versus stress relations $A = k\sqrt{\Sigma_c}$, from classical wrinkling theory³², with $k=1$ for labyrinths and $k=2/\sqrt{11+6v-5v^2}$ for hexagonal patterns.

Received 17 August 2014; accepted 19 December 2014; published online 2 February 2015

## References

1. Bond, J. R., Kofman, L. & Pogosyan, D. How filaments of galaxies are woven into the cosmic web. *Nature* **380**, 603–606 (1996).

2. Parteli, E., Durán, O., Tsoar, H., Schwämmle, V. & Herrmann, H. J. Dune formation under bimodal winds. *Proc. Natl Acad. Sci. USA* **106**, 22085–22089 (2009).

3. Onuki, A. *Phase Transition Dynamics* (Cambridge Univ. Press, 2002).

4. Palmer, A. R. Symmetry breaking and the evolution of development. *Science* **306**, 828–833 (2004).

5. Chirat, R., Moulton, D. E. & Goriely, A. Mechanical basis of morphogenesis and convergent evolution of spiny seashells. *Proc. Natl Acad. Sci. USA* **110**, 6015–6020 (2013).

6. Aranson, I. S. & Tsimring, L. S. Patterns and collective behavior in granular media: Theoretical concepts. *Rev. Mod. Phys.* **78**, 641–692 (2006).

7. Ginzburg, V. L. & Landau, L. D. On the theory of superconductivity. *Zh. Eksp. Teor. Fiz.* **20**, 1064–1082 (1950).

8. Turing, A. M. The chemical basis of morphogenesis. *Phil. Trans. R. Soc. Lond. B* **237**, 37–72 (1952).

9. Ouyang, Q. & Swinney, H. L. Transition from a uniform state to hexagonal and striped Turing patterns. *Nature* **352**, 610–612 (1991).

10. Swift, J. & Hohenberg, P. C. Hydrodynamic fluctuations at the convective instability. *Phys. Rev. A* **15**, 319–328 (1977).

11. Paczuski, M., Kardar, M. & Nelson, D. R. Landau theory of the crumpling transition. *Phys. Rev. Lett.* **60**, 2638–2640 (1988).

12. Delprato, A. M., Samadani, A., Kudrolli, A. & Tsimring, L. S. Swarming ring patterns in bacterial colonies exposed to ultraviolet radiation. *Phys. Rev. Lett.* **87**, 158102 (2001).

13. Brenner, M. P. Chemotactic patterns without chemotaxis. *Proc. Natl Acad. Sci. USA* **107**, 11653–11654 (2010).

14. Wiggins, S. & Golubitsky, M. *Introduction to Applied Nonlinear Dynamical Systems and Chaos* Vol. 2 (Springer-Verlag, 1990).

15. Cerda, E. & Mahadevan, L. Geometry and physics of wrinkling. *Phys. Rev. Lett.* **90**, 074302 (2003).

16. Efimenko, K. et al. Nested self-similar wrinkling patterns in skins. *Nature Mater.* **4**, 293–297 (2005).

17. Richman, D. P., Stewart, R. M., Hutchinson, J. W. & Caviness, V. S. Mechanical model of brain convolutional development. *Science* **189**, 18–21 (1975).

18. Bowden, N., Brittain, S., Evans, A. G., Hutchinson, J. W. & Whitesides, G. M. Spontaneous formation of ordered structures in thin films of metals supported on an elastomeric polymer. *Nature* **393**, 146–149 (1998).

19. Chan, E. P. & Crosby, A. J. Fabricating microlens arrays by surface wrinkling. *Adv. Mater.* **18**, 3238–3242 (2006).

20. Terwagne, D., Brojan, M. & Reis, P. M. Smart morphable surfaces for aerodynamic drag control. *Adv. Mater.* **26**, 6608–6611 (2014).

21. Jiang, H. et al. Finite deformation mechanics in buckled thin films on compliant supports. *Proc. Natl Acad. Sci. USA* **104**, 15607–15612 (2007).

22. Pocivavsek, L. et al. Stress and fold localization in thin elastic membranes. *Science* **320**, 912–916 (2008).

23. Brau, F. et al. Multiple-length-scale elastic instability mimics parametric resonance of nonlinear oscillators. *Nature Phys.* **7**, 56–60 (2011).

24. Xu, F., Potier-Ferry, M., Belouettar, S. & Cong, Y. 3d finite element modeling for instabilities in thin films on soft substrates. *Int. J. Solids Struct.* **51**, 3619–3632 (2014).

25. Brau, F., Damman, P., Diamant, H. & Witten, T. A. Wrinkle to fold transition: Influence of the substrate response. *Soft Matter* **9**, 8177–8186 (2013).

26. Cao, G., Chen, X., Li, C., Ji, A. & Cao, Z. Self-assembled triangular and labyrinth buckling patterns of thin films on spherical substrates. *Phys. Rev. Lett.* **100**, 036102 (2008).

27. Li, B., Jia, F., Cao, Y-P., Feng, X-Q. & Gao, H. Surface wrinkling patterns on a core-shell soft sphere. *Phys. Rev. Lett.* **106**, 234301 (2011).

28. Breid, D. & Crosby, A. J. Curvature-controlled wrinkle morphologies. *Soft Matter* **9**, 3624–3630 (2013).

29. Schroll, R. D. et al. Capillary deformations of bendable films. *Phys. Rev. Lett.* **111**, 014301 (2013).
30. King, H., Schroll, R. D., Davidovitch, B. & Menon, N. Elastic sheet on a liquid drop reveals wrinkling and crumpling as distinct symmetry-breaking instabilities. *Proc. Natl Acad. Sci. USA* **109**, 9716–9720 (2012).

31. Audoly, B. & Boudaoud, A. Buckling of a thin film bound to a compliant substrate, part i: Formulation, linear stability of cylindrical patterns, secondary bifurcations. *J. Mech. Phys. Solids* **56**, 2401–2421 (2008).

32. Cai, S., Breid, D., Crosby, A. J., Suo, Z. & Hutchinson, J. W. Periodic patterns and energy states of buckled films on compliant substrates. *J. Mech. Phys. Solids* **59**, 1094–1114 (2011).

33. Ciarlet, P. G. *Mathematical Elasticity* Vol. 3 (North Holland, 2000).

34. Yin, J., Han, X., Cao, Y. & Lu, C. Surface wrinkling on polydimethylsiloxane microspheres via wet surface chemical oxidation. *Sci. Rep.* **4**, 5710 (2014).

35. Jost, J. *Riemannian Geometry and Geometric Analysis* (Springer, 2008).

36. Cross, M. C. & Hohenberg, P. C. Pattern formation outside of equilibrium. *Rev. Mod. Phys.* **65**, 851–1112 (1993).

37. Allen, H. G. *Analysis and Design of Structural Sandwich Panels* (Pergamon, 1969).

38. Golovin, A. A. & Nepomnyashchy, A. A. *Self-Assembly, Pattern Formation and Growth Phenomena in Nano-Systems* (Springer, 2006).

39. Coxeter, H. S. M. *Regular Polytopes* (Courier Dover Publications, 1973).

40. Clerc, M. G., Petrossian, A. & Residori, S. Bouncing localized structures in a liquid-crystal light-valve experiment. *Phys. Rev. E* **71**, 015205(R) (2005).

41. Kozyreff, G. & Tlidi, M. Nonvariational real Swift–Hohenberg equation for biological, chemical, and optical systems. *Chaos* **17**, 037103 (2007).

42. Lloyd, D. J. B., Sandstede, B., Avitabile, D. & Champneys, A. R. Localized hexagon patterns of the planar Swift–Hohenberg equation. *SIAM J. Appl. Dyn. Syst.* **7**, 1049–1100 (2008).

43. Burke, J. & Knobloch, E. Localized states in the generalized Swift–Hohenberg equation. *Phys. Rev. E* **73**, 056211 (2006).

44. Cirak, F., Ortiz, M. & Schröder, P. Subdivision surfaces: A new paradigm for thin-shell finite-element analysis. *Int. J. Numer. Methods Eng.* **47**, 2039–2072 (2000).

45. Stoop, N., Wittel, F. K., Amar, M. B., Müller, M. M. & Herrmann, H. J. Self-contact and instabilities in the anisotropic growth of elastic membranes. *Phys. Rev. Lett.* **105**, 068101 (2010).

## Acknowledgements

This work was supported by the Swiss National Science Foundation grant No. 148743 (N.S.), by the National Science Foundation, CAREER CMMI-1351449 (P.M.R.) and by an MIT Solomon Buchsbaum Award (J.D.).

## Author contributions

N.S., R.L. and J.D. developed the theory. N.S. and R.L. performed analytical calculations. N.S. implemented and performed the numerical simulations. D.T. and P.M.R. developed the experiments. N.S., R.L. and D.T. analysed data. All authors discussed the results and contributed to writing the paper.

## Additional information

Supplementary information is available in the online version of the paper. Reprints and permissions information is available online at www.nature.com/reprints. Correspondence and requests for materials should be addressed to J.D.

## Competing financial interests

The authors declare no competing financial interests.
Curvature-induced symmetry breaking determines elastic surface patterns

Norbert Stoop¹, Romain Lagrange¹, Denis Terwagne²*, Pedro M. Reis²,³, and Jörn Dunkel¹

¹Department of Mathematics, Massachusetts Institute of Technology,
77 Massachusetts Avenue, Cambridge, MA 02139-4307, USA

²Department of Civil & Environmental Engineering, Massachusetts Institute of Technology,
77 Massachusetts Avenue, Cambridge, MA 02139-4307, USA

³Department of Mechanical Engineering, Massachusetts Institute of Technology,
77 Massachusetts Avenue, Cambridge, MA 02139-1713, USA

(Dated: December 14, 2014)

We derive an effective theory for the wrinkling of thin hard films bound to arbitrarily curved soft substrates. Starting from the nonlinear Koiter shell equations, we show that the elastic equations can be reduced to a generalized Swift-Hohenberg theory, yielding Eq. (1) in the Main Text for the special case of a spherical surface geometry. Using nonlinear analysis of this effective fourth-order equation, we derive predictions for hexagonal and labyrinth-like wrinkling patterns in dependence on the film stress and the substrate curvature. To illustrate the effects of spatially varying substrate curvature on wrinkling, we present additional numerical results for a toroidal geometry.

DEFINITIONS

Let $\mathbf{S} = \mathbf{\Theta}(\theta_1, \theta_2)$ be a two-dimensional surface in $\mathbb{R}^3$, parameterized by $y = (\theta_1, \theta_2) \in \omega \subset \mathbb{R}^2$. Throughout, Greek indices $\alpha, \beta, \dots$ take values in $\{1, 2\}$, whereas Latin indices $i, j, \dots$ run from 1 to 3. The induced metric $a_{\alpha\beta}$ (first fundamental form) on the surface $\mathbf{S} = \mathbf{\Theta}(\theta_1, \theta_2)$ is given by

$$a_{\alpha\beta} = \mathbf{a}_{\alpha} \cdot \mathbf{a}_{\beta} = a_{\beta\alpha}, \qquad (1)$$

where

$$\mathbf{a}_{\alpha} = \mathbf{\Theta}_{,\alpha} \equiv \partial_{\alpha} \mathbf{\Theta} \equiv \frac{\partial \mathbf{\Theta}}{\partial \theta_{\alpha}} \qquad (2)$$

are the tangent vectors, and $\cdot$ denotes the Euclidean in- ner product on $\mathbb{R}^3$. The unit-length normal vector $\mathbf{n}$ is defined by

$$\mathbf{n} \equiv \mathbf{a}^3 = \frac{\mathbf{a}_1 \times \mathbf{a}_2}{|\mathbf{a}_1 \times \mathbf{a}_2|} \qquad (3)$$

and characterized by the properties

$$\begin{align*}
\mathbf{n} \cdot \mathbf{n} &= 1, & \mathbf{n} \cdot \mathbf{a}_{\alpha} &= 0, & \mathbf{n}_{,\alpha} \cdot \mathbf{n} &= 0 \\
\mathbf{n}_{,\alpha} \cdot \mathbf{a}_{\beta} &= -\mathbf{a}_{\alpha,\beta} \cdot \mathbf{n}, & \mathbf{n} \cdot \mathbf{n}_{,\alpha\beta} &= -\mathbf{n}_{,\alpha} \cdot \mathbf{n}_{,\beta}
\end{align*}$$

The surface element is

$$d\omega = \sqrt{|\det(a_{\alpha\beta})|} dy \qquad (4)$$

We also introduce the second and third fundamental forms $b_{\alpha\beta}$, $c_{\alpha\beta}$ with components given by

$$b_{\alpha\beta} = \mathbf{n} \cdot \mathbf{a}_{\alpha,\beta} \qquad (5a)$$

$$c_{\alpha\beta} = \mathbf{n}_{,\alpha} \cdot \mathbf{n}_{,\beta} \qquad (5b)$$

The second fundamental form, $b_{\alpha\beta}$, is often also referred to as the curvature tensor. The Christoffel symbols are

$$\Gamma^\sigma_{\alpha\beta} = \frac{1}{2} a^{\sigma\gamma} (a_{\gamma\alpha,\beta} + a_{\gamma\beta,\alpha} - a_{\alpha\beta,\delta}) \qquad (6)$$

where $a^{\alpha\beta}$ are the components of the contravariant metric tensor, defined by $a^{\alpha\gamma}a_{\beta\gamma} = \delta^\alpha_\beta$. Introduc- ing $\mathbf{a}^\alpha = a^{\alpha\beta}\mathbf{a}_\beta$, the following identities will be useful later:

$$\mathbf{n}_{\alpha} = -b_{\alpha\beta}\mathbf{a}^{\beta} = -b_{\alpha}^{\sigma}\mathbf{a}_{\sigma} \qquad (7a)$$

$$b_{\alpha}^{\gamma} b_{\beta\gamma} = \mathbf{n}_{,\alpha} \cdot \mathbf{n}_{,\beta} = c_{\alpha\beta} \qquad (7b)$$

The covariant derivative of a scalar function $\psi$ is

$$\nabla_{\alpha} \psi = \psi_{,\alpha} \qquad (8)$$

The gradient of $\psi$ on the surface $\mathbf{S}$ has components

$$\nabla^{\alpha} \psi = a^{\alpha\beta} \nabla_{\alpha} \psi \qquad (9)$$

The action of the Laplace-Beltrami operator $\Delta$ on scalar functions is defined as

$$\Delta \psi = \nabla_{\alpha} \nabla^{\alpha} \psi = a^{\gamma\delta} \psi_{,\gamma\delta} - a^{\gamma\delta} \Gamma^{\lambda}_{\gamma\delta} \psi_{,\lambda} \qquad (10)$$

For a vector field $V^\alpha$ or for a $(0, 1)$-tensor field $V_\alpha$, the covariant derivative involves the Christoffel symbols,

$$\nabla_{\alpha} V^{\beta} = V_{,\alpha}^{\beta} + \Gamma^{\alpha}_{\beta\gamma} V^{\gamma} \qquad (11a)$$

$$\nabla_{\alpha} V_{\beta} = V_{,\beta,\alpha} - \Gamma^{\gamma}_{\alpha\beta} V_{\gamma} \qquad (11b)$$

All gradients $\nabla$ and Laplacians $\Delta$ below refer to these surface-specific differential operators.

<footnote>*Current address: Faculté des Sciences, Université Libre de Brux-elles (ULB), Bruxelles 1050, Belgium</footnote>
ENERGY FUNCTIONAL

We extend the classical Koiter shell (KS) energy func-
tional to account for film-substrate coupling and excess
stresses.

Koiter shell equations

The KS equations describe the equilibrium of a thin
shell (precurved plate) when the thickness *h* of the shell
is small compared to its curvature in undeformed and
deformed configurations. The KS equations follow rigor-
ously by means of Γ-convergence from the full 3D elas-
ticity problem in the limit *h* → 0 [1].

In the absence of forces and boundary conditions, we
assume that the shell adopts a stress-free equilibrium
configuration which we call the reference configuration,
parametrized by the surface map Θ, with fundamental
forms a<sub>αβ</sub>, b<sub>αβ</sub> etc. as specified above. Under the influ-
ence of forces and boundary conditions, the shell adopts a
new, deformed configuration characterized by a displace-
ment field Ψ defined with respect to the curved reference
state Θ,

$$
\Psi = \Psi_1 a^1 + \Psi_2 a^2 + \Psi_3 a^3 \equiv \Psi_i a^i \quad (12)
$$

For any point $y \in \omega$, its displaced position is given by
$\Theta(y) + \Psi(y)$, and the respective surface geometry will be
denoted as $a_{\alpha\beta}(\Psi)$, $b_{\alpha\beta}(\Psi)$ etc.

The KS energy of the shell is given by [1]

$$
\mathcal{E}_{\text{KS}}(\mathbf{\Psi}) = \mathcal{E}_b(\mathbf{\Psi}) + \mathcal{E}_s(\mathbf{\Psi}) + \mathcal{E}_f(\mathbf{\Psi}) \quad (13\text{a})
$$

with bending energy

$$
\mathcal{E}_b = \frac{E_f}{2(1 - \nu^2)} \int_{\omega} d\omega \frac{h^3}{24} C^{\alpha\beta\gamma\delta} R_{\gamma\delta}(\mathbf{\Psi}) R_{\alpha\beta}(\mathbf{\Psi}) \quad (13\text{b})
$$

stretching energy

$$
\mathcal{E}_s = \frac{E_f}{2(1 - \nu^2)} \int_{\omega} d\omega \frac{h}{2} C^{\alpha\beta\gamma\delta} G_{\gamma\delta}(\Psi) G_{\alpha\beta}(\Psi) \quad (13c)
$$

and energy contributions

$$
\mathcal{E}_f = - \int_{\omega} d\omega \, p^i \Psi_i \qquad (13\text{d})
$$

due to external forces $p^i$ (pressure, body loads, etc.). In
Eqs. (13), $E_f$ denotes the Young modulus of the film,
$\nu$ its Poisson ratio and $C^{\alpha\beta\gamma\delta}$ the constitutive tensor.
We focus on the case of a Kirchhoff-St. Venant material,
corresponding to an extension of Hook's law to large de-
formations, described by

$$
C^{\alpha\beta\gamma\delta} = (1 - \nu)(a^{\alpha\delta}a^{\beta\gamma} + a^{\alpha\gamma}a^{\beta\delta}) + 2\nu a^{\alpha\beta}a^{\gamma\delta} \quad (14)
$$

The nonlinear membrane bending and stretching strains
$R_{\alpha\beta}$ and $G_{\alpha\beta}$ are given by [1]

$$
R_{\alpha\beta} = b_{\alpha\beta}(\mathbf{\Psi}) - b_{\alpha\beta} \quad (15a)
$$

$$
G_{\alpha\beta} = \frac{1}{2} [a_{\alpha\beta}(\mathbf{\Psi}) - a_{\alpha\beta}] \quad (15\text{b})
$$

The bending energy scales with $h^3$ and will be small
compared to the stretching contributions. We will there-
fore linearize $R_{\alpha\beta}(\Psi)$ in the bending energy, but keep
higher-order terms in the stretching strains $G_{\alpha\beta}$.

*Normal component of bending strains.* The linearized bending strains read [1]

$$
R_{\alpha\beta} \approx (\Psi_{,\alpha\beta} - \Gamma_{\alpha\beta}^{\sigma} \Psi_{,\sigma}) \cdot \mathbf{n} \quad (16)
$$

With this approximation, the normal displacement com-
ponent $\Psi_3$ decouples from the in-plane components.
Since the dominant bending contribution comes from the
normal displacement $\Psi_3$, we may neglect the in-plane
components

$$
\begin{align}
R_{\alpha\beta} \approx \rho_{\alpha\beta}(\Psi_3) &\equiv \Psi_{3,\alpha\beta} - \Gamma_{\alpha\beta}^{\sigma} \Psi_{3,\sigma} - \Psi_3 c_{\alpha\beta} \nonumber \\
&= \nabla_{\alpha} \nabla_{\beta} \Psi_3 - \Psi_3 c_{\alpha\beta} \tag{17}
\end{align}
$$

Normal component of stretching strains. For the
stretching strains, one has [1]

$$
\begin{align}
G_{\alpha\beta} &= \frac{1}{2} (\nabla_{\beta}\Psi_{\alpha} + \nabla_{\alpha}\Psi_{\beta} + \Psi_{,\alpha} \cdot \Psi_{,\beta}) - b_{\alpha\beta}\Psi_3 \nonumber \\
&= \frac{1}{2} (\Psi_{\alpha,\beta} + \Psi_{\beta,\alpha} + \Psi_{,\alpha} \cdot \Psi_{,\beta}) - \Gamma_{\alpha\beta}^{\sigma} \Psi_{\sigma} - b_{\alpha\beta}\Psi_3 \tag{18}
\end{align}
$$

The displacement derivative can be split into an in-plane
and normal part,

$$
\mathbf{\Psi}_{,\alpha} = (\Psi_{\gamma,\alpha} - \Gamma_{\alpha\gamma}^{\sigma} \Psi_{\sigma} - b_{\alpha\gamma} \Psi_3) \mathbf{a}^{\gamma} + (\Psi_{3,\alpha} + b_{\alpha}^{\gamma} \Psi_{\gamma}) \mathbf{a}^3
$$

Due to the orthogonality of **a**<sup>α</sup> and **n** ≡ **a**<sup>3</sup>, the nonlinear
term in the stretching strains becomes

$$
\mathbf{\Psi}_{,\alpha} \cdot \mathbf{\Psi}_{,\beta} = t_{\alpha}^{\delta} t_{\beta\delta} + s_{\alpha} s_{\beta} \quad (19\text{a})
$$

where

$$
t_{\alpha\beta} = \Psi_{\beta,\alpha} - \Gamma_{\alpha\beta}^{\sigma} \Psi_{\sigma} - b_{\alpha\beta} \Psi_3 \quad (19\text{b})
$$

$$
s_{\alpha} = \Psi_{3,\alpha} + b_{\alpha}^{\sigma} \Psi_{\sigma} \quad (19c)
$$

Expanding Eqs. (19) for small in-plane displacements,
$|\Psi_\alpha| \ll |\Psi_3|$, one finds to leading order

$$
\mathbf{\Psi}_{,\alpha} \cdot \mathbf{\Psi}_{,\beta} \approx \Psi_{3,\alpha} \Psi_{3,\beta} + b_{\alpha}^{\delta} b_{\beta\delta} (\Psi_3)^2
$$

Using Eq. (7) and (18), we obtain

$$
G_{\alpha\beta} \approx \gamma_{\alpha\beta}(\Psi_3) = \frac{1}{2} [\Psi_{3,\alpha}\Psi_{3,\beta} + c_{\alpha\beta}(\Psi_3)^2] - b_{\alpha\beta}\Psi_3 \quad (20)
$$

We note that the obtained bending and stretching
strains are symmetric,

$$
\rho_{\alpha\beta} = \rho_{\beta\alpha}, \quad \gamma_{\alpha\beta} = \gamma_{\beta\alpha} \tag{21}
$$
Additional remarks. In our and previous [2] experi- ments, the film stress is imposed in two ways: (i) depres- surization of the shell-substrate system and (ii) swelling of the film. Before the onset of the buckling transition, either technique creates a pre-stress in the film. The asso- ciated prestrain can be decomposed into an in-plane part ¯ γ$_{αβ}$ and a bending part ¯ ρ$_{αβ}$ , with ¯ γ$_{αβ}$ ≫ ¯ ρ$_{αβ}$ . Adopt- ing the same approximations as above (Eqs. 17, 20), the prestrains remain symmetric. Restricting ourselves fur- thermore to equi-biaxial pre-stress without shear, ¯ γ$_{αβ}$ is proportional to the metric tensor a$_{αβ}$ , and we write this as

$$
\bar{\gamma}_{\alpha\beta} = \bar{\gamma} a_{\alpha\beta} \qquad (22)
$$

where from now on ¯ γ denotes the constant of proportion- ality.

If the film stress due to depressurization or swelling becomes too large, the film will buckle, with an inho- mogeneous deformation u around the pre-stressed state. The strains of a buckled configuration can then be ex- pressed as

$$
\gamma_{\alpha\beta}(u_0 + u) = \bar{\gamma}_{\alpha\beta} + \frac{1}{2} (u_{,\alpha}u_{,\beta} - 2b_{\alpha\beta}u + c_{\alpha\beta}u^2) \quad (23\text{a})
$$

$$
\rho_{\alpha\beta}(u_0 + u) = \bar{\rho}_{\alpha\beta} + \nabla_{\alpha} \nabla_{\beta} u - c_{\alpha\beta} u \quad (23\text{b})
$$

To contract the strains with the constitutive tensor, we use the fact that the bending and stretching strains, Eqs. (23), are symmetric. For any symmetric (0,2)- tensor τ$_{αβ}$ , the contraction C ( τ ) of C αβγδ with τ$_{αβ}$ can be written as

$$
\begin{align}
C(\tau) &= C^{\alpha\beta\gamma\delta} \tau_{\alpha\beta} \tau_{\gamma\delta} \nonumber \\
&= 2 \left[ (1-\nu)\tau^{\alpha\beta} \tau_{\alpha\beta} + \nu(\tau_{\gamma}^{\gamma})^2 \right] \nonumber \\
&= 2(\tau_M + \tau_G) \tag{24}
\end{align}
$$

with

$$
\begin{align}
\tau_M &= (\tau_\gamma^\gamma)^2 \nonumber \\
\tau_G &= (1 - \nu) [\tau^{\alpha\beta} \tau_{\alpha\beta} - (\tau_\gamma^\gamma)^2] \tag{25}
\end{align}
$$

For later use, it is convenient to introduce the identities

$$
\mathcal{H} = \frac{1}{2} b_{\gamma}^{\gamma} \qquad (26a)
$$

$$
\mathcal{K} = \det(\{a_{\alpha\beta}\}) / \det(\{b_{\alpha\beta}\}) \qquad (26b)
$$

$$
\mathcal{R} \equiv b^{\alpha\beta} b_{\alpha\beta} = c_{\gamma}^{\gamma} = 4\mathcal{H}^2 - 2\mathcal{K} \qquad (26c)
$$

$$
S \equiv b^{\alpha\beta} c_{\alpha\beta} = 2\mathcal{H}(4\mathcal{H}^2 - 3\mathcal{K}) \quad (26d)
$$

$$
\mathcal{T} \equiv c^{\alpha\beta} c_{\alpha\beta} = 16\mathcal{H}^2(\mathcal{H}^2 - \mathcal{K}) + 2\mathcal{K}^2 \quad (26e)
$$

H is the mean curvature and K the Gaussian curvature. Most of these expressions follow directly from the con- tracted Gauss-Codazzi-Mainardi equations [3].

Relevant energy contributions

Bending energy density. Using the definition of the Laplace-Beltrami operator, Eq. (10), and Eq. (24) with τ$_{αβ}$ = ρ$_{αβ}$ , the mean-curvature contribution ρ$_{M}$ can be written as

$$
\begin{equation}
\begin{split}
\rho_M = (\rho_\gamma^\gamma)^2 = (\Delta u)^2 & - 2\mathcal{R}u\Delta u + \mathcal{R}^2 u^2 + \\
& 2(\bar{\rho}_\gamma^\gamma \Delta u - \mathcal{R}\bar{\rho}_\gamma^\gamma u) + (\bar{\rho}_\gamma^\gamma)^2
\end{split}
\tag{27}
\end{equation}
$$

The ρ$_{G}$ -term in Eq. (24) accounts for the energy cost due to a change of Gaussian curvature. This term is negligible if the typical wrinkling wavelengths are small compared to the local radii of curvature of the underlying surface. More precisely, in this case, one finds that Eqs. (15a,17) reduce to ρ β α ≃ b β α ( u ) in leading order, where b β α ( u ) is the curvature tensor of the deformed configuration. Accord- ingly, Eq. (25) then yields ρ$_{G}$ = (1 − ν ) K ( u ). Since the Gaussian curvature K ( u ) of the deformed configuration integrates to a topological invariant for a closed surface (or when the geodesic curvature of the boundary curve is fixed), the ρ$_{G}$ -term can be neglected in the variational formulation, just as in the classical Helfrich model [4]. Under these assumptions, the resulting bending energy takes the form

$$
\begin{equation}
\begin{split}
\mathcal{E}_b = \frac{E_f}{2(1 - \nu^2)} \int_{\omega} d\omega \frac{h^3}{12} \biggl[ & (\Delta u)^2 + 2\mathcal{R}(\nabla u)^2 + \\
& (\mathcal{R}^2 - \Delta \mathcal{R})u^2 + \\
& 2(\Delta \bar{\rho}_{\gamma}^{\gamma} - \mathcal{R}\bar{\rho}_{\gamma}^{\gamma})u + (\bar{\rho}_{\gamma}^{\gamma})^2 \biggr]
\end{split}
\tag{28}
\end{equation}
$$

where from now on the product symbol · as in ( ∇ u ) 2 = ( ∇ u ) · ( ∇ u ) denotes the scalar product with respect to the surface metric. To obtain Eq. (28), we used the gener- alized Stokes theorem 1 to rewrite the second and fourth term on the rhs. of Eq. (27).

Stretching energy density. For the stretching energy

$$
\mathcal{E}_s = \frac{E_f}{2(1 - \nu^2)} \int_{\omega} d\omega \frac{h}{2} \mathcal{C}(\gamma) \quad (29)
$$

we find by using Eq. (24) with τ$_{αβ}$ = γ$_{αβ}$ , and combining with Eqs. (22), (23) and (26),

$$
\begin{align}
\mathcal{C}(\gamma) ={}& 2(1 + \nu)\bar{\gamma}(\nabla u)^2 + 2(1 + \nu)(\nabla u)^4 - \nonumber \\
& 2 \left[ (1 - \nu)b^{\alpha\beta} \nabla_{\alpha} u \nabla_{\beta} u + 2\nu\mathcal{H}(\nabla u)^2 \right] u + \nonumber \\
& \left[ (1 - \nu)c^{\alpha\beta} \nabla_{\alpha} u \nabla_{\beta} u + \nu\mathcal{R}(\nabla u)^2 \right] u^2 + \nonumber \\
& 2 \left[ (1 - \nu + \bar{\gamma} + \nu\bar{\gamma})\mathcal{R} + 4\nu\mathcal{H}^2 \right] u^2 - \nonumber \\
& 2 \left[ (1 - \nu)\mathcal{S} + 2\nu\mathcal{H}\mathcal{R} \right] u^3 + \nonumber \\
& \frac{1}{2} \left[ (1 - \nu)\mathcal{T} + \nu\mathcal{R}^2 \right] u^4 + \mathcal{O}(1) + \mathcal{O}(u) \tag{30}
\end{align}
$$

<footnote n="1">1 Note that − ∫ dω R u Δ u = ∫ dω u ( ∇ R ) · ( ∇ u ) + ∫ dω R ( ∇ u ) $^{2}$. Furthermore, repeated application of the Stokes theorem gives ∫ dω u ( ∇ R ) · ( ∇ u ) = 1 2 ∫ dω Δ R u $^{2}$.</footnote>
where we did not explicitly write down the terms linear and constant in $u$, as they will not be relevant for the later analysis (see detailed remarks in Total energy density below).

**Substrate coupling energy.** In our experiments, the thin film is coupled to a curved soft substrate. To simplify further analysis, we assume from now on that the substrate has the same Poisson ratio $\nu$ as the film, as is the case in our experiments. We model the substrate coupling as a nonlinear spring by adding a substrate energy $\mathcal{E}_{sub}$ to the KS energy from Eq. (13a), where

$$ \mathcal{E}_{sub} = \frac{E_s}{2} \int_{\omega} d\omega \left( Au + \frac{\tilde{a}}{h} u^2 + \frac{\tilde{c}}{h^3} u^4 \right) \quad (31) $$

with $E_s$ denoting the Young modulus of the substrate. The constant film thickness $h$ could have been absorbed into the coefficients $\tilde{a}$ and $\tilde{c}$, but simplifies subsequent formulas. Note that $\mathcal{E}_{sub}$ contains a term linear in $u$ because we are considering the state of the film-substrate system around a flat but displaced equilibrium solution $u_0$. This linear term gives rise to a corresponding constant normal force that is needed to balance the internal normal forces of the film.

**Energy due to excess film stress.** Finally, we still have to account for the excess film stress

$$ \Sigma_e \equiv \frac{\sigma}{\sigma_c} - 1 \quad (32a) $$

where $\sigma$ is the film stress and $\sigma_c$ the critical stress needed for wrinkling. In our model, the energy due to excess film stress is included by adding a term

$$ \mathcal{E}_{\sigma} = \frac{E_f}{2(1 - \nu^2)} \int_{\omega} d\omega \frac{\tilde{a}_2}{h} \Sigma_e u^2 \quad (32b) $$

to the KS energy from Eq. (13a). The energy contribution $\mathcal{E}_{\sigma}$ is crucial for capturing the system behavior beyond the wrinkling instability. We discuss below how the dimensionless parameter $\tilde{a}_2$ is related to the elastic properties of the substrate. The $u^2$-dependence of $\mathcal{E}_{\sigma}$ is a classical result from elastic wrinkling theory [5], ensuring that the amplitude-stress relationship in the effective model agrees with classical wrinkling theory, as is shown in detail further below.

**Total energy density.** Adding the contributions due to substrate coupling, external forces and excess stress to the KS energy (13a), we obtain the total elastic energy

$$ \mathcal{E} = \frac{E_f}{1 - \nu^2} \bar{\mathcal{E}} \quad (33a) $$

where to leading order

$$ \bar{\mathcal{E}} = \int_{\omega} d\omega \left[ \frac{\gamma_0}{2} (\nabla u)^2 + \frac{\gamma_2}{2} (\Delta u)^2 + \frac{a}{2} u^2 + \frac{b}{3} u^3 + \frac{c}{4} u^4 - \frac{h}{2} \left[ (1 - \nu) b^{\alpha\beta} \nabla_{\alpha} u \nabla_{\beta} u + 2\nu \mathcal{H}(\nabla u)^2 \right] u + \frac{h}{4} \left[ (1 - \nu) c^{\alpha\beta} \nabla_{\alpha} u \nabla_{\beta} u + \nu \mathcal{R}(\nabla u)^2 \right] u^2 \right] \quad (33b) $$

with coefficients

$$ \gamma_0 = h\bar{\gamma}(1 + \nu) + \frac{h^3 \mathcal{R}}{6} $$

$$ \gamma_2 = \frac{h^3}{12} $$

$$ a = h \left[ (1 - \nu + \bar{\gamma} + \nu \bar{\gamma}) \mathcal{R} + 4\nu \mathcal{H}^2 \right] + \frac{h^3 (\mathcal{R}^2 - \Delta \mathcal{R})}{12} + \frac{E_s (1 - \nu^2)}{E_f} \frac{\tilde{a}}{h} + \frac{\tilde{a}_2}{h} \Sigma_e $$

$$ b = -\frac{3h}{2} \left[ (1 - \nu) \mathcal{S} + 2\nu \mathcal{H} \mathcal{R} \right] $$

$$ c = \frac{h}{2} \left[ (1 - \nu) \mathcal{T} + \nu \mathcal{R}^2 \right] + 2 \frac{E_s (1 - \nu^2)}{E_f} \tilde{c} \quad (33c) $$

Note that for compressive stresses $\bar{\gamma} < 0$. Thus, for sufficiently large film pre-stress, $\gamma_0 < 0$. To obtain the effective energy functional (33), the following additional simplifications and assumptions were adopted:

* **Constant terms.** We neglected all constant terms in the energy, as they will not contribute to the equations of motion, obtained by variation of the energy with respect to $u$.

* **Terms linear in $u$.** We note that the term linear in $u$ gives rise to a inhomogeneous, constant term in the equation of motion. However, $u = 0$ always is an equilibrium solution by construction. More precisely, $u = 0$ means that the film is radially displaced by $u_0$, which is a fundamental solution of the problem. Therefore, the inhomogeneous term in the equation of motion has to vanish, implying that the coefficient of the energy term linear in $u$ must be zero. This condition can be interpreted as follows: For $u = 0$ to be an equilibrium solution, the sum of all normal forces acting on the film must vanish.

* **Quartic terms.** The quartic terms in $u$ and $\nabla u$ ensure that the effective theory remains stable above the wrinkling threshold, as these terms limit the growth of the most unstable modes. To keep the theory as simple as possible, we only include the dominant $u^4$-contribution and neglect terms $\propto (\nabla u)^4$.

Taking the variation of $\bar{\mathcal{E}}$ with respect to $u$, we obtain

$$ -\frac{\delta \bar{\mathcal{E}}}{\delta u} = \gamma_0 \Delta u - \gamma_2 \Delta^2 u - au - bu^2 - cu^3 + \frac{h}{2} \left\{ (\nu - 1) \left[ b^{\alpha\beta} \nabla_{\alpha} u \nabla_{\beta} u + 2u \nabla_{\beta} (b^{\alpha\beta} \nabla_{\alpha} u) \right] + 2\nu \left[ \mathcal{H}(\nabla u)^2 - 2\nabla \cdot (\mathcal{H} u \nabla u) \right] \right\} + \frac{h}{2} \left[ (1 - \nu) u \nabla_{\beta} (uc^{\alpha\beta} \nabla_{\alpha} u) - \nu \mathcal{R} u (\nabla u)^2 + \nu \nabla \cdot (\mathcal{R} u^2 \nabla u) \right] \quad (34) $$
where $\nabla \cdot$ denotes the surface divergence. We note that, for a flat Euclidean metric, the first line of Eq. (34) coincides with the Swift-Hohenberg (SH) equation, as originally derived in the context of Rayleigh-Bénard convection [6, 7]. The b-term and the first term $\propto h$ (second and third line) break the symmetry $u \to -u$, which is known to lead to a transition from labyrinth-like patterns to hexagons for SH-like equations. As both terms depend on the curvature tensor $b_{\alpha\beta}$, we can infer that regions of high curvature will show different wrinkling patterns than regions of low curvature.

We next apply Eq. (34) to derive quantitative predictions for wrinkling patterns transitions on spherical geometries, which are then compared with our experiments (Main Text). Subsequently, we still present numerical solutions for a toroidal geometry, as an example of a surface with locally varying curvature.

APPLICATION TO SPHERICAL GEOMETRIES

Using spherical coordinates $(\theta_1, \theta_2) \in [0, 2\pi) \times [0, \pi]$, a spherical surface of radius $R$ is described by the metric tensor

$$ (a_{\alpha\beta}) = \begin{pmatrix} (R \sin \theta_2)^2 & 0 \\ 0 & R^2 \end{pmatrix} \qquad (35) $$

In this case, we have $\mathcal{H} = -1/R$, $\mathcal{K} = 1/R^2$. Thus

$$ \mathcal{R} = 2/R^2, \quad \mathcal{S} = -2/R^3, \quad \mathcal{T} = 2/R^4 \qquad (36) $$

*Total energy density.* Noting that the sphere has constant mean and Gaussian curvature $\mathcal{H}$ and $\mathcal{K}$, Eq. (36) allows to simplify the total energy density (33) considerably,

$$ \bar{\mathcal{E}} = \int_{\omega} d\omega \left[ \frac{\gamma_0}{2} (\nabla u)^2 + \frac{\gamma_2}{2} (\Delta u)^2 + \frac{a}{2} u^2 + \frac{b}{3} u^3 + \frac{c}{4} u^4 + \Gamma_1 (\nabla u)^2 u + \frac{\Gamma_2}{2} (\nabla u)^2 u^2 \right] \qquad (37) $$

with coefficients

$$ \begin{aligned} \gamma_0 &= h\bar{\gamma}(1 + \nu) + \frac{h^3}{3R^2} \\ \gamma_2 &= \frac{h^3}{12} \\ a &= \frac{2h(1 + \bar{\gamma})(1 + \nu)}{R^2} + \frac{h^3}{3R^4} + \frac{E_s(1 - \nu^2)}{E_f}\frac{\tilde{a}}{h} + \frac{\tilde{a}_2}{h}\Sigma_e \\ b &= \frac{3h(1 + \nu)}{R^3} \\ c &= \frac{h(1 + \nu)}{R^4} + 2\frac{E_s(1 - \nu^2)}{E_f}\tilde{c} \\ \Gamma_1 &= \frac{h(1 + \nu)}{2R} \\ \Gamma_2 &= \frac{h(1 + \nu)}{2R^2} \end{aligned} \qquad (38) $$

Equations of motions

To identify the equilibrium configurations, we assume that the film exhibits an overdamped relaxation dynamics. Then, the equations of motion follow by functional variation of the elastic energy (37) with respect to the displacement field $u$,

$$ \frac{\rho}{\tau_0} \partial_t u = - \frac{\delta \mathcal{E}}{\delta u} \qquad (39) $$

where $\rho$ is the constant surface mass density of the film and $\tau_0$ the damping-time scale. The relaxation dynamics (39) can be written in the equivalent form

$$ \mu \partial_t u = - \frac{\delta \bar{\mathcal{E}}}{\delta u} \qquad (40a) $$

where the coefficient

$$ \mu = \frac{\rho(1 - \nu^2)}{\tau_0 E_f} \qquad (40b) $$

is the inverse relaxation speed. Calculating the functional derivative $\delta \bar{\mathcal{E}}/\delta u$ gives

$$ \mu \partial_t u = \gamma_0 \Delta u - \gamma_2 \Delta^2 u - au - bu^2 - cu^3 + \Gamma_1 [(\nabla u)^2 + 2u\Delta u] + \Gamma_2 [u(\nabla u)^2 + u^2\Delta u] \qquad (41) $$

Since we are only interested in the steady-state solutions, the exact value of $\mu$ is not relevant for our analysis. It is convenient to rewrite Eq. (41) in dimensionless form by measuring length in units of the film thickness $h$ and time in units of $\tau_h = \mu h$. Introducing the dimensionless curvature parameter

$$ \kappa = h/R \qquad (42) $$

Eq. (41) reduces to

$$ \partial_t u = \gamma_0 \Delta u - \gamma_2 \Delta^2 u - au - bu^2 - cu^3 + \Gamma_1 [(\nabla u)^2 + 2u\Delta u] + \Gamma_2 [u(\nabla u)^2 + u^2\Delta u] \qquad (43) $$

with rescaled dimensionless parameters

$$ \begin{aligned} \gamma_0 &= \bar{\gamma}(1 + \nu) + \frac{\kappa^2}{3} < 0 \\ \gamma_2 &= \frac{1}{12} \\ a &= 2(1 + \bar{\gamma})(1 + \nu)\kappa^2 + \frac{\kappa^4}{3} + \frac{E_s(1 - \nu^2)}{E_f}\tilde{a} + \tilde{a}_2\Sigma_e \\ b &= 3(1 + \nu)\kappa^3 \\ c &= (1 + \nu)\kappa^4 + 2\frac{E_s(1 - \nu^2)}{E_f}\tilde{c} \\ \Gamma_1 &= \frac{(1 + \nu)\kappa}{2} \\ \Gamma_2 &= \frac{(1 + \nu)\kappa^2}{2} \end{aligned} \qquad (44) $$
Note that the covariant derivatives $\nabla$ and $\Delta$ in Eq. (43) are now also defined with respect to the rescaled dimensionless sphere of radius $\kappa^{-1} = R/h$. In Eq. (44), we kept the term $\sim \kappa^4$ in the coefficient $a$. Although this higher order term is negligible, its inclusion will lead to simpler expressions when matching our model with experiments (see section *Curvature-dependence of the critical strain* $\bar{\gamma}$ below). Given the dimensionless parameters in Eq. (44), the corresponding values in physical units are recovered through the transformations

$$ 
\begin{aligned}
u &\rightarrow hu, & R &\rightarrow h/\kappa, & t &\rightarrow \mu h t \\
\gamma_0 &\rightarrow h\gamma_0, & \gamma_2 &\rightarrow h^3\gamma_2, & & \\
a &\rightarrow a/h, & b &\rightarrow b/h^2, & c &\rightarrow c/h^3 \\
\Gamma_1 &\rightarrow \Gamma_1, & \Gamma_2 &\rightarrow \Gamma_2/h & &
\end{aligned} 
\quad (45) 
$$

As evident from Eq. (44), the model is specified through dimensionless parameters

$$ (\kappa, \bar{\gamma}, E_s/E_f, \nu, \tilde{a}, \tilde{a}_2, \tilde{c}) \quad (46) $$

### Parameter determination

The parameters ($h, R, \nu, E_s, E_f$) can be directly measured for our experimental system. To determine the remaining parameters ($\bar{\gamma}, \tilde{a}, \tilde{a}_2, \tilde{c}$), we proceed as follows:

1. Linear stability analysis will enable us to relate $\gamma_0$ at the onset of wrinkling with the wavelength $\lambda$, which yields a relation between the critical buckling strain $\bar{\gamma}$ and the ratio $E_s/E_f$.

2. The value of the substrate parameter $\tilde{a}$ can be estimated from known results for the critical bulking stress in planar elasticity theory [8]. Below, we will extend the classical derivation to the weakly curved case to confirm that our model predictions agree with recent results by Cai et al. [9]

3. By means of nonlinear stability analysis and comparison with analytical results for the standard Swift-Hohenberg equation, we will express the parameter $\tilde{a}_2$ in terms of $c$, leaving $\tilde{c}$ as the only remaining fit parameter. We estimate $\tilde{c}$ by comparing our numerical simulations with the experimentally measured surface morphologies.

### Critical stress $\bar{\gamma}$ and $\tilde{a}$

We estimate $\bar{\gamma}$ by comparing our effective theory with known results for the full elastic equations in the planar limit case $R \rightarrow \infty$. Letting $\kappa \rightarrow 0$ and linearizing Eq. (43) for a small perturbation $\epsilon e^{ikx}$ of the unbuckled

homogeneous solution, one finds the dominant unstable wave-mode

$$ |k| = \sqrt{\frac{|\gamma_0|}{2\gamma_2}} = \sqrt{6|\gamma_0|} \quad (47a) $$

which selects the characteristic wavelength

$$ \lambda_c = \frac{2\pi}{|k|} = \frac{2\pi}{\sqrt{6|\gamma_0|}} \quad (47b) $$

Equating $\lambda_c$ with the known wrinkling wavelength $\lambda_{el}$ of a planar elastic film-substrate system, which in units $h = 1$ is given by [8]

$$ \lambda_{el} = 2\pi \left( \frac{E_f}{3E_s} \right)^{1/3} \quad (48) $$

we obtain

$$ \gamma_0 = -\frac{1}{6} \left( \frac{3E_s}{E_f} \right)^{2/3} \quad (49) $$

The sign indicates a compressive strain, which in our terminology is negative. From Eq. (44) with $\kappa \rightarrow 0$, we find the planar estimate

$$ \bar{\gamma} \approx -\frac{1}{6(1+\nu)} \left( \frac{3E_s}{E_f} \right)^{2/3} \equiv \bar{\gamma}_p \quad (50) $$

A correction due to curvature will be discussed below.

**Critical stress.** The strain $\bar{\gamma}_p$ can be associated with the critical stress $\sigma_c$ at the wrinkling threshold. As expressed by Eq. (22), our system is in a state of equibiaxial strains $\bar{\gamma}$ implying that, in a locally orthogonal coordinate frame, the in-plane elasticity tensor reduces to $\epsilon_{11} = \epsilon_{22} = \epsilon$, $\epsilon_{12} = \epsilon_{21} = 0$. The usual stress-strain relationship of a Hookean material then reads [10]

$$ 
\begin{aligned}
\sigma_{11} &= \sigma_{22} = \frac{E}{1 - \nu^2}(\epsilon_{11} + \nu\epsilon_{22}) = \frac{E}{1 - \nu}\epsilon \equiv \sigma \\
\sigma_{12} &= \sigma_{21} = 0
\end{aligned} 
$$

Assuming a standard linear relation between stress and strain, we expect

$$ \sigma_c = k_\sigma \frac{E_f}{1 - \nu} \bar{\gamma}_p \quad (51) $$

with some constant prefactor $k_\sigma$. Inserting Eq. (50), we obtain

$$ \sigma_c = -\frac{k_\sigma E_f}{6(1-\nu^2)} \left( \frac{3E_s}{E_f} \right)^{2/3} \quad (52) $$

which for $k_\sigma = 3/2$ agrees with the known critical stress of elastic wrinkling analysis [8]

$$ \sigma_{c_{el}} = -\frac{E_f}{4(1-\nu^2)} \left( \frac{3E_s}{E_f} \right)^{2/3} \quad (53) $$
Estimation of $\tilde{a}$ near the critical value $\sigma_c$. In the planar limit $\kappa \to 0$, Eq. (43) exhibits a bifurcation from a uniform state to nontrivial pattern formation only if $a < a_c$, where

$$a_c = \frac{\gamma_0^2}{4\gamma_2} = \frac{1}{12} \left( \frac{3E_s}{E_f} \right)^{4/3} \quad (54)$$

At the wrinkling onset, corresponding to $a = a_c$, the film stress $\sigma$ equals the critical stress $\sigma_c$ so that $\Sigma_e = 0$. For the planar case, Eq. (54) thus determines the substrate parameter $\tilde{a}$ as

$$\tilde{a} = \frac{1}{4(1 - \nu^2)} \left( \frac{3E_s}{E_f} \right)^{1/3} \quad (55)$$

### Curvature-dependence of the critical strain $\tilde{\gamma}$

Recent simulations of the full coupled elasticity equations [11] and experiments with polymer colloids [12] report wave-length reduction of pattern on curved substrates compared with the planar case. It is therefore interesting to study how the critical strain $\tilde{\gamma} < 0$ depends on the curvature parameter $\kappa = h/R$ in our model.

Similar to the planar case, cf. Eq. (54), the wrinkling bifurcation occurs when

$$\frac{\gamma_0^2(\kappa)}{4\gamma_2 a_c(\kappa)} = 1. \quad (56)$$

Recalling that $\Sigma_e = 0$ at the transition point and using the above result for $\tilde{a}$, we can solve Eq. (56) for the critical strain $\tilde{\gamma}$. Using Eq. (44) we then obtain for $\gamma_0$

$$\gamma_0 = \frac{\kappa^2}{3} - \frac{1}{6} \sqrt{\left(\frac{3E_s}{E_f}\right)^{4/3} + 24(1 + \nu)\kappa^2} \quad (57)$$

which reduces to Eq. (49) in the planar case ($\kappa = 0$). For $\tilde{\gamma}$, we find to $O(\kappa^4)$

$$\frac{\tilde{\gamma}}{\tilde{\gamma}_p} = 1 + 12\kappa^2(1 + \nu) \left(\frac{E_f}{3E_s}\right)^{4/3} \quad (58)$$

where $\tilde{\gamma}_p$ is the critical buckling strain for the planar case, given in Eq. (50). This asymptotic scaling behavior in $\kappa$ is similar to the results of Cai et al. [9], although the numerical prefactors and the dependence on the Poisson ratio $\nu$ differ.

According to Eq. (57), the absolute value $|\gamma_0|$ increases with curvature with a leading order correction $\sim \kappa^2$. Equation (47b) then implies that the wavelength $\lambda_c$ decreases with increasing curvature, in qualitative agreement with experimental and numerical findings [11, 12]. It is interesting to note that the wavelength reduction in our model is due to the nonlinear, curvature-dependent stretching terms of the underlying KS energy, unlike the

<table><caption>TABLE I: List of parameters for Eqs. (43) and (64) as obtained by systematic asymptotic matching to classical elastic wrinkling theory, with $\eta = 3E_s/E_f$, $\gamma_2 = 1/12$, $\Sigma_e = (\sigma/\sigma_c) - 1$, and $\kappa = h/R$ where $h$ is the film thickness and $R$ the radius (see Fig. 2 of Main Text). We substituted $\tilde{a}$ and $\tilde{c}$ by Eqs. (55) and (77). Focusing on the leading order contribution, we only kept terms up to $O(\kappa^3)$, cf. Eq. (57). The only remaining fitting parameter of the model is $c_1$.</caption><tbody><tr><td>γ<sub>0</sub></td><td>=</td><td>-η<sup>2/3</sup>/6 - [2(1 + ν) - 1/3]κ<sup>2</sup></td></tr><tr><td>a</td><td>=</td><td>η<sup>4/3</sup>/12 + (6(1 + ν) - η<sup>2/3</sup>)κ<sup>2</sup> + ˜a<sub>2</sub>Σ<sub>e</sub></td></tr><tr><td>b</td><td>=</td><td>3(1 + ν)κ<sup>3</sup></td></tr><tr><td>c</td><td>=</td><td>2(1 + ν)η<sup>2/3</sup>/3 c<sub>1</sub></td></tr><tr><td>Γ<sub>1</sub></td><td>=</td><td>(1 + ν)/2 κ</td></tr><tr><td>Γ<sub>2</sub></td><td>=</td><td>(1 + ν)/2 κ<sup>2</sup></td></tr><tr><td>˜a<sub>2</sub></td><td>=</td><td>-η<sup>4/3</sup>(c + 3|γ<sub>0</sub>|Γ<sub>2</sub>)/48γ<sub>0</sub><sup>2</sup></td></tr></tbody></table>

model of Yin et al. [13], where a wavelength reduction was obtained by assuming a curvature-dependent substrate model. Unfortunately, for the range of parameters realized in our experiments the curvature-induced wavelength reduction is below the detection threshold. We therefore did not include any curvature-dependent substrate response in our model.

### Nonlinear behavior above onset

Having determined estimates for the parameters ($\tilde{\gamma}, \tilde{a}$) by analyzing the onset of wrinkling, the two remaining unknown parameters are ($\tilde{a}_2, \tilde{c}$). Aiming to further reduce the number of free parameters, we now turn to the regime beyond the wrinkling threshold, where patterns are selected by nonlinear effects. To this end, we first reduce the generalized Swift-Hohenberg equation (43) to a standard Swift-Hohenberg (SH) equation by approximating mixed $\Gamma_{1,2}$-terms in Eq. (43) that contain both $u$ and $\nabla u$ through effective expressions that only contain $u$. Assuming a typical relation between pattern amplitude and excess film stress $\Sigma_e = (\sigma/\sigma_c) - 1$, we can then exploit existing results for the stability of patterns in the SH equation to predict the morphological phase diagram of the wrinkling patterns in our experimental system.

**Swift-Hohenberg approximation.** To approximate Eq. (43) by a standard SH equation, we recall that $\gamma_0$ and $\gamma_2$ select the dominant (most unstable) wave number vector $k_c = \pm\sqrt{6|\gamma_0|}$, see Eq. (47a). Considering the