# Mathematical Aspects of ANM/FEM Numerical Model, Applied to Nonlinear Elastic, and Thermo Elastic Analysis of Wrinkles in Film/Substrate Systems, and a New Implementation in the FreeFEM++ Language

Pascal Ventura, Frédéric Hecht, Michel Potier-Ferry, Hamid Zahrouni, Fan Xu, Hamza Azzayani, Michael Brun, Anh-Khoa Chau

► To cite this version:

Pascal Ventura, Frédéric Hecht, Michel Potier-Ferry, Hamid Zahrouni, Fan Xu, et al., Mathematical Aspects of ANM/FEM Numerical Model, Applied to Nonlinear Elastic, and Thermo Elastic Analysis of Wrinkles in Film/Substrate Systems, and a New Implementation in the FreeFEM++ Language. Mathematics, 2025, 13 (19), pp.3063. 10.3390/math13193063. hal-05284277

HAL Id: hal-05284277

https://hal.univ-lorraine.fr/hal-05284277v1

Submitted on 29 Sep 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Mathematical Aspects of ANM/FEM Numerical Model, Applied to Nonlinear Elastic, and Thermo Elastic Analysis of Wrinkles in Film/Substrate Systems, and a New Implementation in the FreeFEM++ Language

Pascal Ventura <fnref n="1" />,*,† ID, Frédéric Hecht <fnref n="2" />, Michel Potier-Ferry <fnref n="1" />,† ID, Hamid Zahrouni <fnref n="1" />, Fan Xu <fnref n="3" /> ID, Hamza Azzayani <fnref n="1" /> ID, Michael Brun <fnref n="1" /> and Anh-Khoa Chau <fnref n="1" />

## Abstract

The main purposes of the present paper are to present the mathematical and algorithmic aspects of the ANM/FEM numerical model and to show how it is applied to analyze elastic and thermo-elastic nonlinear solid mechanical problems. ANM is a robust continuation method based on a perturbation technique for solving nonlinear problems dependent on a loading parameter. Historically, this technique has been successfully applied to problems in various fields of solid and fluid mechanics. This paper shows how ANM is used to solve nonlinear elastic and nonlinear thermo-elastic problems involving elastic behavior and geometrical nonlinearities. The implementation of ANM for FEM in the FreeFEM++ language is then presented. The FEM software development platform, called FreeFEM++, is structured to work with variational formulations and, therefore, is well adapted to implement ANM for instability problems in solid mechanics. In order to illustrate the great efficiency of FreeFEM++, scripts will be presented for computing the different steps of ANM continuation for solid elastic structures, considering simple geometries subjected to conservative loading. For the purpose of validation, the problem of a cantilever subjected to an applied force is presented. Next, the new numerical model is applied to study wrinkles appearing in a planar film/substrate system that is subjected to compressive surface forces at the lateral faces of the film. Finally, the model is applied to a spherical film/substrate system subjected to thermo-elastic shrinkage. In both cases, the ANM/FEM prediction method, together with a Newton-Riks correction (if needed), identifies the equilibrium paths efficiently, especially after the post-buckling regime.

**Keywords:** implicit function theorem; partial differential equations; ANM; FEM; Domain-Specific Language; FreeFEM++; film/substrate system; wrinkling

**MSC:** 35B32; 35B36; 35B38; 47J25

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

check for
updates

Academic Editor: Andrey Amosov

Received: 30 June 2025

Revised: 5 September 2025

Accepted: 15 September 2025

Published: 23 September 2025

**Citation:** Ventura, P.; Hecht, F.; Potier-Ferry, M.; Zahrouni, H.; Xu, F.; Azzayani, H.; Brun, M.; Chau, A.-K. Mathematical Aspects of ANM/FEM Numerical Model, Applied to Nonlinear Elastic, and Thermo Elastic Analysis of Wrinkles in Film/Substrate Systems, and a New Implementation in the FreeFEM++ Language. Mathematics 2025, 13, 3063. https://doi.org/10.3390/math13193063

**Copyright:** © 2025 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

<footnote n="1">¹ Arts et Métiers Institute of Technology, Université de Lorraine, CNRS, LEM3, F-57000 Metz, France; michel.potier-ferry@univ-lorraine.fr (M.P.-F.); hamid.zahrouni@univ-lorraine.fr (H.Z.); hamza.azzayani@univ-lorraine.fr (H.A.); michael.brun@univ-lorraine.fr (M.B.); anh-khoa.chau@univ-lorraine.fr (A.-K.C.)</footnote>

<footnote n="2">² Université Paris Cité, Sorbonne Université, LJLL, CNRS, Inria, Alpines, F-75005 Paris, France; frederic.hecht@sorbonne-universite.fr</footnote>

<footnote n="3">³ Department of Aeronautics and Astronautics, Institute of Mechanics and Computational Engineering, College of Intelligent Robotics and Advanced Manufacturing, Fudan University, Shanghai 200433, China; fanxu@fudan.edu.cn</footnote>

<footnote>* Correspondence: pascal.ventura@univ-lorraine.fr; Tel.: +33-682349360</footnote>

<footnote>† These authors contributed equally to this work.</footnote>
# 1. Introduction

Film/substrate (bi-layer) systems composed of a thin film deposited on a thicker substrate are present in numerous applications [1]. For example, under compression loading, these systems can exhibit instabilities such as wrinkling [2], folding [3], and creasing [4]. Even though these instabilities are often avoided, they can be useful in determining mechanical properties. The bi-layer systems usually used in engineering applications are characterized by a stiff film deposited on a more compliant substrate, as for sandwich panels [5]. This paper will focus on wrinkling instability. Wrinkling is a mechanical instability due to stresses appearing when a stiff thin layer deposited on a compliant substrate is compressed. When subjected to such a loading, the initially smooth surface bifurcates towards a periodic pattern with a short wavelength. This phenomenon occurs both in nature and in technology. At the nanometric scale, the wrinkling of carbon nanotubes appears under loading conditions [6]. In the case of materials used in flexible electronics, pre-stressed stiff layers of elastometer will wrinkle under compressive loads [7]. In the field of biology, wrinkles can happen on the skin [8] during the aging process, or as a consequence of multiple movements, and also in the brain tissues [9]. The interest in wrinkling mechanisms has greatly increased. It is now considered more as a functional feature than a failure [10].

The research has focused on the exploration of complex auto-organized patterns resulting from the wrinkling of thin films on an elastomer substrate. For example, it has been shown that wrinkling patterns can be very useful as photonic structures for solar panels [11].

More precisely, the surface morphological instabilities of a soft material with a stiff thin surface layer have been extensively studied over the past years. Depositing a stiff thin film on a polymeric soft substrate may result in residual compressive stresses in the film during the cooling process. This can occur when there is a large mismatch in the thermal expansion coefficients of the film and the substrate. These stresses are relieved by buckling with wrinkling patterns, pioneered by Bowden et al. [12] in 1998. Note that wrinkles of a thin stiff film deposited on a soft substrate can be widely used in industry for many applications, as in the micro/nano-fabrication of flexible electronic devices with controlled morphological patterns [13,14], or the mechanical property measurement of material characteristics [15].

In most research, the resolution of the problem in 2D and 3D is achieved either using spectral methods or with the Fast Fourier Transform algorithm. Because of the need for periodic boundary conditions, these methods face some limitations. Later, Chen et al. [16] showed that film/substrate systems can be analyzed using the Finite Element Method, allowing for complex geometries.

Most of the methods used to follow the equilibrium path of a nonlinear system of equations are iterative from an initial solution guess. The best-known are the Newton-Raphson, or pseudo arc-length methods that use the Jacobian matrix interpolation. When studying wrinkling phenomena, the numerical modeling of the instabilities may be complicated due to the existence of many equilibrium solutions and the low stiffness of unloaded membranes. In such situations, it is not easy to apply a classical continuation algorithm that may fail [17,18] or undergo back-and-forth loops. As explained below, because of the use of high-order Taylor expansion, the ANM continuation algorithm may be more efficient for following equilibrium paths when other techniques fail to converge. Such a difficulty may also be avoided by using a pseudo-dynamic algorithm, e.g., [18-20]. ANM is a continuation method based on the Taylor series. A key property of ANM is the a posteriori determination of the domain of validity of the step length for each ANM branch. It is often associated with classical discretization techniques, generally incorporating the Finite Element Method (FEM). Most of the explanation of the properties of ANM is given in [21].
Recently, Potier-Ferry [22] wrote a review article that summarizes the key points of ANM and its applications to hyperelasticity and plasticity. In his article, the main mathematical aspects of ANM are presented. In this paper, these are developed (for example, for the needed mathematical spaces) and applied to the resolution (with FEM) of nonlinear elasticity or nonlinear thermo-elasticity problems using the formulation of internal virtual work based on the second Piola–Kirchhoff stress tensor and the Green–Lagrange strain tensor.

The coupling between a perturbation technique and FEM was first proposed by Thompson and Walker [23]. In [24], a continuation method was introduced, together with an efficient procedure to build the Taylor series and the convergence acceleration using Padé approximants. This has been successfully applied to various engineering and scientific problems in solid and fluid mechanics (for example, [25]). It is different from classical incremental-iterative methods such as Riks strategy [26] in that full solution paths are computed, rather than pointwise solutions, and because the a posteriori estimation of the radius of convergence ensures an optimal loading step size.

ANM has been implemented using various computer languages: FORTRAN 77, FORTRAN 90, MATLAB 2024, and C++98. Many home-made implementations were established. For example, the open-source software MANLAB [27] is widely used in the field of nonlinear vibrations. Lejeune et al. [28] described a full object-oriented implementation within C++. Nonetheless, the home-made software is difficult to maintain and often provides a framework that is not sufficiently user-friendly. Application in fluid mechanics within the open-source finite-element platform ELMER was made with results presented in [29].

FreeFEM++ [30] is an efficient environment for implementing variational formulations. As we will show in this paper, FreeFEM++ uses a high-level language that makes it easy to derive a variational formulation and implement the ANM algorithm. The goal is to develop a general, and easy-to-use, environment to implement the ANM algorithm in order to solve solid mechanical problems. FreeFEM++ also offers MPI parallel solvers and multigrid or domain decomposition tools that reduce the CPU time and allow a large number of degrees of freedom to be taken into account. FreeFEM++ has been used to develop a numerical tool for the analysis of infinite periodic surface acoustic wave transducers and more recently to solve numerically fluid-structure interaction problems. In this paper, the implementation of ANM in the FreeFEM++ code is discussed.

Section 2 will present the main mathematical aspects of ANM, the application of ANM in the cases of nonlinear solid elastic problems and thermo-elastic nonlinear problems. In this section, it is explained what makes FreeFEM++ an efficient tool to implement ANM for various nonlinear mechanical problems, and the implementation of the ANM algorithm in FreeFEM++ is described. Sections 3 and 4 present the simulation results, using the FreeFEM++ implementation of ANM for the nonlinear mechanical behavior of a clamped cantilever and subjected to a vertical surface force at its right surface, for a planar film/substrate system under uni-axial loading, and for a spherical film/substrate system subjected to a thermo-mechanical shrinkage of its core [18].

## 2. Methods

### 2.1. Theoretical Aspects of ANM

Beginning with the resolution of Equation (1) in the general case of a Banach space $B$:

$$ \mathbf{R}(\mathbf{U}, \lambda) = 0 \qquad (1) $$
with, $\mathbf{R} : B \times \mathbb{R} \to B$, $\mathbf{U} \in B$, and $\lambda \in \mathbb{R}$

Since the problem (1) is formulated for $\mathbf{U}$ belonging to a Banach space, the following results are valid for Sobolev spaces used in the weak formulations of the Finite Element Method.
For the purpose of simplification, problems from engineering and physics will be considered and formulated in a finite dimensional space using a numerical method like the Finite Element Method. In the following, we define $\mathbf{R}: \mathbb{R}^n \times \mathbb{R} \to \mathbb{R}^n$, $\mathbf{U} \in \mathbb{R}^n$, and $\lambda \in \mathbb{R}$.

Assume that the function $\mathbf{R}(\mathbf{U}, \lambda)$ is a class-$C^p$ function (with $p \ge 1$) from $H \times I$, where $H$ is an open set of $\mathbb{R}^n$, and $I$ is an open segment of $\mathbb{R}$, and assume that, for $(\mathbf{U}_0, \lambda_0) \in H \times I$, a solution of (1), the partial differential $\frac{\partial \mathbf{R}}{\partial \mathbf{U}}|_{(\mathbf{U}_0, \lambda_0)}$ is a bi-continuous isomorphism. The implicit function theorem [31] states that Equation (1) has a unique solution branch, $\phi: \lambda \to \mathbf{U} = \phi(\lambda)$, from a neighborhood of $\lambda_0$ in $I$ to $H$, which is also a class-$C^p$ function.

Furthermore, if we assume that $\mathbf{R}(\mathbf{U}, \lambda)$ is an analytic function, the function $\phi(\lambda)$ is also analytic.

Note that, because of the nonlinearity of the problem, the existence and uniqueness of the solution $\mathbf{U} = \phi(\lambda)$ is established only locally by the implicit function theorem, but not globally. In the presence of bifurcations, the Jacobian matrix $\frac{\partial \mathbf{R}}{\partial \mathbf{U}}|_{(\mathbf{U}_0, \lambda_0)}$ will be not invertible, at the critical point, and several forked branches may exist.

By the definition of an analytic function, it is possible to expand $\mathbf{U}$ in series, according to $\lambda$. A path parameter, $a$, is now defined, and both $\mathbf{U}$ and $\lambda$ are expanded into Taylor series at the order $N$, according to $a$ within a given radius of convergence [22,32,33]:

$$ \left\{ \begin{array}{l} \mathbf{U}(a) = \mathbf{U}_0 + \sum_{i=1}^{N} a^i \mathbf{U}_i \\ \\ \lambda(a) = \lambda_0 + \sum_{i=1}^{N} a^i \lambda_i \end{array} \right. \qquad (2) $$

Let us explain the algebraic properties of the ANM. If we consider a given function, $\hat{\mathbf{Y}}: \mathbf{X} \in \mathbb{R}^m \to \mathbf{Y} = \hat{\mathbf{Y}}(\mathbf{X}) \in \mathbb{R}^n$, as detailed in [22], if $\mathbf{X}$ is a function of the path parameter $a$, then because of the chained rule, the i-th derivative $\mathbf{Y}^{(i)}$ of $\mathbf{Y}(a)$ with respect to $a$ is a function of all the k-th derivatives $\mathbf{X}^{(k)}$ of $\mathbf{X}(a)$ with respect to $a$ for $k \in \{0, 1, \dots, i\}$. Furthermore, we have the following general expression:

$$ \mathbf{Y}^{(i)} = \mathbf{J}(\mathbf{X}_0) \cdot \mathbf{X}^{(i)} + \mathbf{Y}_i^{nl} (\mathbf{X}^{(0)}, \mathbf{X}^{(1)}, \dots, \mathbf{X}^{(i-1)}) \qquad (3) $$

Note that there is a linear relationship between $\mathbf{Y}^{(i)}$ and $\mathbf{X}^{(i)}$ with the same Jacobian matrix $\mathbf{J}(\mathbf{X}_0)$ for each order $i$ and that $\mathbf{Y}_i^{nl}$ is a multi-linear function of previous orders, $k \in \{0, 1, \dots, i-1\}$. As a consequence, the matrix $\mathbf{J}(\mathbf{X}_0) = D_{\mathbf{X}}(\mathbf{Y})$ can be obtained from (3) with $i=1$ since $\mathbf{Y}_1^{nl} = 0$, and, at order $i$ ($1 < i \le N$), with the assumption $\mathbf{X}^{(i)} = 0$, the expression $\mathbf{Y}_i^{nl}(\mathbf{X}^{(0)}, \mathbf{X}^{(1)}, \dots, \mathbf{X}^{(i-1)})$ is given by $\mathbf{Y}^{(i)}$.

## 2.2. ANM for Elastic Nonlinear Problems

### 2.2.1. Theory

Let us explain how the ANM can be applied in the case of a quasi-static geometrically nonlinear elastic problem [24,34].

In this document, the Einstein summation convention is used, and the double contracted product, operating on second-rank tensors is defined, in $\mathbb{R}^3$, by $A: B = \sum_{i=1}^3 \sum_{j=1}^3 A_{i,j}B_{i,j}$. A Lagrangian formulation is used. $\Omega_0$ is the reference configuration (at time $t=0$) for a material domain, and $\Omega$ is its deformed configuration (at current parameter $\lambda$). If $\mathbf{X}$ is the position of a particle in the reference configuration, and $\mathbf{x} = \chi(\mathbf{X}, \lambda)$, the position of the same particle in the deformed configuration, with $\chi(\mathbf{X}, \lambda)$, is a sufficiently regular function from $\Omega_0 \times \mathbb{R} \to \Omega$. For a given parameter, $\lambda$, the displacement field, defined on the open set $\Omega_0$, is $\mathbf{u} = \mathbf{x}(\mathbf{X}, \lambda) - \mathbf{X}$, and it belongs to the Sobolev space $H_1(\Omega_0)$, satisfying the kinematic boundary conditions.
The Green–Lagrange strain tensor $\gamma(\mathbf{u})$ and the second Piola–Kirchhoff stress tensor $\mathbf{S}$, both second-rank tensors, are very well defined in classical continuum mechanics books [21,35]. Here, the material is supposed to obey the Saint Venant–Kirchhoff law, which is characterized by a linear relation between $\mathbf{S}$ and $\gamma(\mathbf{u})$. The expression of $\gamma(\mathbf{u})$ is given as follows:

$$ \gamma(\mathbf{u}) = \underbrace{\frac{1}{2}(\nabla\mathbf{u} + \nabla^T\mathbf{u})}_{\gamma_l(\mathbf{u})} + \underbrace{\frac{1}{2}(\nabla^T\mathbf{u} \cdot \nabla\mathbf{u})}_{\gamma_{nl}(\mathbf{u},\mathbf{u})} \quad (4) $$

and its variation is given as follows:

$$ \delta\gamma(\mathbf{u}) = \gamma_l(\delta\mathbf{u}) + 2\gamma_{nl}(\mathbf{u}, \delta\mathbf{u}) \quad (5) $$

The fourth-rank elastic stiffness tensor $\mathbf{D}$ is used to define a linear relationship between stress and strain:

$$ \mathbf{S} = \mathbf{D} : \gamma(\mathbf{u}) \quad (6) $$

The formulation of the full boundary value problem is quite classical and combines the elastic constitutive law with the virtual work equation

$$ \left\{ \begin{array}{l} \mathbf{S} = \mathbf{D} : \gamma(\mathbf{u}) \\ \int_{\Omega} \mathbf{S} : \delta\gamma(\mathbf{u}) d\Omega - \lambda \int_{\Gamma_N} \mathbf{t} \cdot \delta\mathbf{u} d\Gamma = 0 \\ + \text{Dirichlet boundary conditions} \end{array} \right. \quad (7) $$

where $\lambda$ is a scalar load parameter, and $\mathbf{t}$ represents the given force on the part $\Gamma_N$ of the boundary $\partial\Omega$.

Each step of ANM relies on a Taylor series, at the order $N$, according to a well-chosen path parameter, $a$:

$$ \left\{ \begin{array}{l} \mathbf{u}(a) = \mathbf{u}_0 + a\mathbf{u}_1 + a^2\mathbf{u}_2 + \dots + a^N\mathbf{u}_N \\ \mathbf{S}(a) = \mathbf{S}_0 + a\mathbf{S}_1 + a^2\mathbf{S}_2 + \dots + a^N\mathbf{S}_N \\ \lambda(a) = \lambda_0 + a\lambda_1 + a^2\lambda_2 + \dots + a^N\lambda_N \end{array} \right. \quad (8) $$

Several efficient choices of the path parameter are available [36]. The most popular one is a linearized arc-length parameter that permits us to easily bypass all the extrema of the response curves. It is given as follows [21]:

$$ a = \frac{1}{\bar{\mathbf{u}}^2} \langle \mathbf{u} - \mathbf{u}_0, \mathbf{u}_1 \rangle + \alpha(\lambda - \lambda_0)\lambda_1 \quad (9) $$

where $\bar{\mathbf{u}}$ (resp. $\alpha$) are normalized displacement vectors (resp. loading parameter). Here, by default, we will use $\bar{\mathbf{u}}^2 = 1$ (resp. $\alpha = 1$).

When the Taylor series (8) are substituted into Equations (7) and (9), the identifying mathematical expressions related to a given order, $p$, in $a$ results in a linear system. Here, the boundary value problem in Equation (7) has a simple structure, with the equations involving only linear and quadratic terms with respect to the unknown fields $(\mathbf{u}, \mathbf{S})$, and this system can be deduced from the usual Leibniz rule to compute high-order derivatives or, equivalently, the Taylor coefficients of a product: $(ab)_p = \sum_{r=0}^{p} a_r b_{p-r}$. First of all, let us identify the coefficients at order 1 in $a$, which yield the classical tangent variational problem, as in Newton–Raphson’s method:

$$ k(\mathbf{u}_1, \delta\mathbf{u}) = \lambda_1 P_e(\delta\mathbf{u}) \quad (10) $$

where the bilinear form $k(\cdot, \cdot)$ is the virtual work of internal elastic forces, given as follows:

$$ k(\mathbf{u}, \delta\mathbf{u}) = \int_{\Omega} (\delta\gamma(\mathbf{u}_0) : \mathbf{D} : (\gamma_l(\mathbf{u}) + 2\gamma_{nl}(\mathbf{u}_0, \mathbf{u})) + \mathbf{S}_0 : 2\gamma_{nl}(\mathbf{u}, \delta\mathbf{u})) d\Omega \quad (11) $$

And the linear form is the virtual work of external forces, given as follows:
$$P_e(\delta \mathbf{u}) = \int_{\Gamma_D} \mathbf{t} \cdot \delta \mathbf{u} d\Gamma \qquad (12)$$

The definition in Equation (9) of the path parameter gives an additional equation:

$$\frac{1}{\mathbf{u}^2} \langle \mathbf{u}_1, \mathbf{u}_1 \rangle + \alpha \lambda_1 \lambda_1 = 1 \qquad (13)$$

At the order $p \ge 2$ in $a$, we also get a variational formulation involving the same bilinear form:

$$k(\mathbf{u}_p, \delta \mathbf{u}) = \lambda_p P_e(\delta \mathbf{u}) + \langle \mathbf{F}_p^{nl}, \delta \mathbf{u} \rangle \qquad (14)$$

The last term is a linear form that depends on the variables at orders $k \le (p-1)$. It is given as follows:

$$\langle \mathbf{F}_p^{nl}, \delta \mathbf{u} \rangle = -2 \sum_{r=1}^{p-1} \int_{\Omega} \mathbf{S}_r : \gamma_{nl}(\mathbf{u}_{p-r}, \delta \mathbf{u}) d\Omega \\ - \sum_{r=1}^{p-1} \int_{\Omega} \gamma_{nl}(\mathbf{u}_r, \mathbf{u}_{p-r}) : \mathbf{D} : \delta \gamma(\mathbf{u}_0) d\Omega \qquad (15)$$

while the stress on order $p$ is as follows:

$$\mathbf{S}_p = \mathbf{D} : \left[ \gamma_l(\mathbf{u}_p) + 2\gamma_{nl}(\mathbf{u}_0, \mathbf{u}_p) + \sum_{r=1}^{p-1} \gamma_{nl}(\mathbf{u}_r, \mathbf{u}_{p-r}) \right] \qquad (16)$$

The last equation is the arc-length condition in Equation (9), which yields the load parameter $\lambda_p$ as a function of the displacement $\mathbf{u}_p$:

$$\frac{1}{\mathbf{u}^2} \langle \mathbf{u}_p, \mathbf{u}_1 \rangle + \alpha \lambda_p \lambda_1 = 0 \qquad (17)$$

Usually, the Taylor expansion order is chosen between 15 and 50 to take advantage of the exponential convergence of the power series. Of course, the path parameter $a$ is required to be inferior to $a_{rc}$, the convergence radius of the Taylor series (8). The truncated Taylor series is tested to verify a given accuracy: the quotient of the norm of the difference between the two last orders in the Taylor series by the norm of the whole Taylor series is tested to be inferior to a given parameter, $\delta$, which results for the following inequality to be verified: $||a^N \mathbf{u}_N|| \le \delta ||a \mathbf{u}_1||$. This will lead to the expression of the path parameter $a_{max}$, called the validity parameter:

$$a_{max} = \left( \delta \frac{\| \mathbf{u}_1 \|}{\| \mathbf{u}_N \|} \right)^{\frac{1}{N-1}} \qquad (18)$$

In this way, the Taylor expansion (8) gives a part of the solution path in the interval $[0, a_{max}]$. The continuation procedure may be very simple, by chaining several series, with the end point of this path being the starting point of the next step. It is also possible to improve the accuracy of the end point by using a corrective step via the Newton method or a high-order method or by applying a convergence acceleration technique [37]. The Newton-Riks correction method can be easily implemented in the FreeFEM++ environment. In conclusion, the Taylor-series computation method was detailed here in the simple case of geometrically nonlinear elasticity, but it is easily extended to more complex models, such as finite-strain elastoplasticity.

### 2.2.2. Numerical Algorithm

Let us now describe how to derive a numerical model from the variational formulation established in the previous chapter using the Finite Element Method (FEM). The finite element implementation of the ANM is described in Section 7 of [21].
The domain $\Omega$ is meshed into tetrahedral finite elements, and a Lagrangian $P_2$ interpolation is used. The current finite element is called $\Omega_e$. Let us write the expression of the bilinear form (11) on the discretized domain:

$$k_h(\mathbf{u}, \delta\mathbf{u}) = \sum_e \int_{\Omega_e} (\delta\gamma(\mathbf{u}_0) : \mathbf{D} : (\gamma_l(\mathbf{u}) + 2\gamma_{nl}(\mathbf{u}_0, \mathbf{u})) + \mathbf{S}_0 : 2\gamma_{nl}(\mathbf{u}, \delta\mathbf{u})) d\Omega \quad (19)$$

Using the well-known $P_2$ basis functions, and following all the steps detailed in Section 7 of [21], it is straightforward to obtain the expression for the elementary stiffness matrix, and the FEM stiffness matrix $[K_t(\mathbf{u}_0, \mathbf{S}_0)]$ is the assembly of the elementary stiffness matrices. The FEM right-hand side vector $\{\mathbf{F}\}$ is the implementation of the linear form $P_e(\delta\mathbf{u})$ (12), and the FEM right-hand side vector $\{\mathbf{F}_p^{nl}\}$ is the implementation of the linear form $\langle \mathbf{F}_p^{nl}, \delta\mathbf{u} \rangle$ (15).

The details of the implementation in FreeFEM++ are given in Appendix A.

The main steps for the implementation of the resolution of a mechanical problem with geometrical nonlinearities using the ANM algorithm in the FreeFEM++ language are described in Section 2.4. The ANM algorithm is detailed below (Algorithm 1):

<caption>**Algorithm 1** Algorithm ANM</caption>

1: **Initialize** $\mathbf{u}_0$ AND $\lambda_0$
2: **for** $i = 1$ until $N$ **do**
3: &nbsp;&nbsp;&nbsp;&nbsp;**Assemble** $[K_t(\mathbf{u}_0, \mathbf{S}_0)]$ AND $\{\mathbf{F}\}$
4: &nbsp;&nbsp;&nbsp;&nbsp;$\{\hat{\mathbf{u}}\} \leftarrow [K_t(\mathbf{u}_0, \mathbf{S}_0)]^{-1}\{\mathbf{F}\}$
5: &nbsp;&nbsp;&nbsp;&nbsp;$$\lambda_1 \leftarrow \frac{1}{\sqrt{1+\{\hat{\mathbf{u}}\}^T\{\hat{\mathbf{u}}\}}}$$
6: &nbsp;&nbsp;&nbsp;&nbsp;$$\{\mathbf{u}_1\} \leftarrow \lambda_1\{\hat{\mathbf{u}}\}$$
7: &nbsp;&nbsp;&nbsp;&nbsp;**FOR EACH GAUSS POINT:** $\{\mathbf{S}_1\} \leftarrow [\mathbf{D}]\{\{\gamma_l(\mathbf{u}_1)\} + \{\gamma_{nl}(\mathbf{u}_0, \mathbf{u}_1)\}\}$
8: &nbsp;&nbsp;&nbsp;&nbsp;**FOR EACH GAUSS POINT:** $\{\mathbf{S}_2^{nl}\} \leftarrow [\mathbf{D}]\{\gamma_{nl}(\mathbf{u}_1, \mathbf{u}_1)\}$
9: &nbsp;&nbsp;&nbsp;&nbsp;**FOR EACH GAUSS POINT:** $\{\mathbf{S}_2^s\} \leftarrow [\mathbf{A}(\mathbf{u}_1)]^T\{\mathbf{S}_1\}$
10: &nbsp;&nbsp;&nbsp;&nbsp;**Assemble** $\{\mathbf{F}_2^{nl}\}$
11: &nbsp;&nbsp;&nbsp;&nbsp;**for** $p = 2$ until $N$ **do**
12: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$\{\mathbf{u}_p^{nl}\} \leftarrow [K_t(\mathbf{u}_0, \mathbf{S}_0)]^{-1}\{\mathbf{F}_p^{nl}\}$$
13: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$\lambda_p \leftarrow -\lambda_1\{\mathbf{u}_p^{nl}\}^T\{\mathbf{u}_1\}$$
14: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$$\{\mathbf{u}_p\} \leftarrow \lambda_p\{\hat{\mathbf{u}}\} + \{\mathbf{u}_p^{nl}\}$$
15: &nbsp;&nbsp;&nbsp;&nbsp;**FOR EACH GAUSS POINT:** $\{\mathbf{S}_p\} \leftarrow [\mathbf{D}]\{\{\gamma_l(\mathbf{u}_p)\} + \{\gamma_{nl}(\mathbf{u}_0, \mathbf{u}_p)\}\} + \{\mathbf{S}_p^{nl}\}$
16: &nbsp;&nbsp;&nbsp;&nbsp;**FOR EACH GAUSS POINT:** $\{\mathbf{S}_{p+1}^{nl}\} \leftarrow [\mathbf{D}]\{\sum_{r=1}^p\{\gamma_{nl}(\mathbf{u}_r, \mathbf{u}_{p+1-r})\}\}$
17: &nbsp;&nbsp;&nbsp;&nbsp;**FOR EACH GAUSS POINT:** $\{\mathbf{S}_{p+1}^s\} \leftarrow \sum_{r=1}^p ([\mathbf{A}(\mathbf{u}_{p+1-r})]^T\{\mathbf{S}_r\})$
18: &nbsp;&nbsp;&nbsp;&nbsp;**Assemble** $\{\mathbf{F}_{p+1}^{nl}\}$
19: &nbsp;&nbsp;&nbsp;&nbsp;**end for**
20: **Compute** $a_{max}$
21: **Compute** $\lambda(a)$ AND $\mathbf{u}(a)$ WITH $a \in [0, a_{max}]$
22: **Compute the normalized residual error**
23: **Actualize** $\mathbf{u}_0$ AND $\lambda_0$
24: **end for**

Algorithm A2 of the Newton-Riks correction is given in Appendix B, as well as Algorithm A1 of Newton-Raphson predictions, which will be used later for comparison purposes.

## 2.3. ANM for Thermo-Elastic Nonlinear Problems

### 2.3.1. Theory

In this section, we show that it is possible to take into account an isotropic thermal expansion (or shrinkage) of an elastic domain with the ANM. The thermal expansion (or shrinkage) is driven by the scalar parameter $\lambda$. As in the previous section, we assume that the domain $\Omega$ is a Saint Venant-Kirchhoff material, and we will consider nonlinear geometric elasticity, as in the previous section.
We will assume that the thermomechanical effects result in an additional contribution to the Green–Lagrange strain tensor (4), and its expression is $\lambda I$ (where I is the identity matrix). Notice that the expression of $\delta\gamma(\mathbf{u})$ (5) is not changed.

Incorporating thermal effects, the virtual work of Equation (7) results in the following:

$$ 
\left\{ 
\begin{array}{l}
\mathbf{S} = \mathbf{D} : (\gamma_l(\mathbf{u}) + \gamma_{nl}(\mathbf{u}, \mathbf{u}) - \lambda I) \\
\int_{\Omega} \mathbf{S} : \delta\gamma(\mathbf{u}) d\Omega = 0 \\
+ \text{Dirichlet boundary conditions}
\end{array}
\right. 
\qquad (20) 
$$

$\lambda > 0$ means a thermo-elastic shrinkage.

Let us present how to incorporate the thermo-elastic formulation within the ANM algorithm. We can easily derive the expression of $S_p$ at the order $p \ge 1$:

$$ S_p = D : \left( \gamma_l(\mathbf{u}_p) + 2\gamma_{nl}(\mathbf{u}_0, \mathbf{u}_p) + \sum_{k=1}^{p-1} \gamma_{nl}(\mathbf{u}_k, \mathbf{u}_{p-k}) - \lambda_p I \right) \qquad (21) $$

Inserting (21) into (20), we obtain at the order $p = 1$:

$$ k(\mathbf{u}_1, \delta\mathbf{u}) = \lambda_1 P_{\text{th}}(\delta\mathbf{u}) \qquad (22) $$

With $k(.,.)$, the same bilinear form as in Section 2.2. The additional equation for the path parameter is the same as Equation (13).

The bilinear form $k(.,.)$ is the classical tangent stiffness operator already given in (11), and the linear form $P_{\text{th}}(.)$ that accounts for the thermomechanical shrinkage is given as follows:

$$ P_{\text{th}}(\delta\mathbf{u}) = \int_{\Omega} I : \mathbf{D} : \delta\gamma(\mathbf{u}_0) d\Omega \qquad (23) $$

Notice that the integrand on the right member of (23) is the trace of the second-rank tensor $(\mathbf{D} : (\gamma_l(\delta\mathbf{u}) + \gamma_{nl}(\mathbf{u}_0, \delta\mathbf{u})))$.

At order $p \ge 2$, we obtain a variational form involving the same bilinear form, and the linear form $\langle \mathbf{F}_p^{nl}, \delta\mathbf{u} \rangle$, accounting for the nonlinear second member, has the same expression as in Section 2.2 (see Equation (15)).

$$ k(\mathbf{u}_p, \delta\mathbf{u}) = \lambda_p P_{\text{th}}(\delta\mathbf{u}) + \langle \mathbf{F}_p^{nl}, \delta\mathbf{u} \rangle \qquad (24) $$

The additional equation for the path parameter is the same as that of Equation (17).

### 2.3.2. Numerical Algorithm

The ANM numerical algorithm for the thermo-elastic nonlinear problem is very close to Algorithm 1 presented in Section 2.3.2. The only difference lies in the FEM right-hand side vector {$\mathbf{F}$}, which is now the implementation of the linear form (23).

### 2.4. Some Aspects of the FreeFEM++ Implementation

The main details of the implementation in FreeFEM++ are given in Appendix A.

By comparison with Matlab, which uses matrices, and C++, which uses arrays, FreeFEM++ manipulates variational formulations as abstract objects. FreeFEM++ achieves a high level of abstraction, optimized for finite elements, linked with up-to-date mathematical libraries (MUMPS, PETSc, ...) that allow it to improve the CPU time.

Note that the FreeFEM++ scripts for the FEM simulation of nonlinear solid mechanical problems using ANM will be available soon on the FreeFEM++ website (https://freefem.org/).

### 2.5. Methods for the Numerical Experiments

This section aims to describe two important numerical experiments in the context of the study of mechanical or thermo-mechanical instabilities of film/substrate systems. The numerical results will be presented in the next section.
### 2.5.1. Film/Substrate Systems Under Uniaxial Compression

We are now going to consider a more complex nonlinear mechanical problem: the film/substrate system. The film/substrate system being studied consists of a stiff film deposited on a soft rectangular parallelepiped substrate (Figure 1). Lateral uniaxial forces are applied on the lateral $Σ_l$ surfaces of the film (Figure 2).

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Schematic of the film/substrate system. $Σ_l$ is the lateral surface (right and left) of the film, and $Σ_d$ is the bottom interface of the substrate.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Schematic of the film/substrate system with an upper view showing the applied uniaxial surface force $λt$ (resp. $-λt$) on the right (resp. left) part of $Σ_l$ surfaces of the film, with $λ$ the loading parameter.</figcaption></figure>

Film/substrate systems have already been studied with the help of ANM using a coupling between shell finite element (Büchter-Ramm [38]) and hexahedron finite elements by Fan Xu. In this section, we will present a numerical model with ANM using only a tetrahedron finite element (due to the FreeFEM++ implementation). As in Fan Xu's research, we will assume geometrical nonlinear elasticity for the stiff film and linear elasticity for the soft substrate. The planes $(Oxz)$, $(Oyz)$ are symmetry planes of the modeled film/substrate system (due to symmetry, only one-fourth of the structure is modeled). Moreover, on the interfaces $Σ_d$ and $Σ_l$, the vertical displacements are assumed to be zero. The length, $L_x$, width, $L_y$, and height, $h_t$ of the film/substrate are, respectively, 1.5 mm, 0.75 mm, and 0.1 mm. The thickness of the film, $h_f$, is 0.001 mm. The substrate and the film are assumed to be homogeneous isotropic and elastic materials. Young's modulus of the substrate and the film are, respectively, $E_s = 1.8$ MPa and $E_f = 1.3 \times 10^5$ MPa). The Poisson ratios of the substrate and the film are, respectively, $v_s = 0.48$ and $v_f = 0.3$. Since the Young's
modulus of the film is much larger than that of the substrate ($E_f/E_s = 0.72 \times 10^5$), the
substrate domain is assumed to have linear elasticity, while the film domain has nonlinear
geometrical elasticity [16].

The number of elements for both the length and width of the film/substrate system is 100. We consider one element in the film thickness and 5 elements in the substrate thickness. P2 Lagrange finite element interpolation is assumed. The mesh consists of 360,000 tetrahedra and 1,575,639 degrees of freedom.

2.5.2. Spherical Film/Substrate Under Thermal Loading

In this section, we study film/substrate systems with a spherical curvature for the case of a thermo-mechanical shrinkage of their compliant core. The details of the continuation algorithm ANM used to study thermo-mechanical problems were presented in Section 2.3. In [18], Xu et al. presented an analytical model, experimental results, and numerical results. The numerical results were obtained with the commercial software ABAQUS using both the Riks method and the pseudo-dynamic method. Mechanical experiments have been made for spherical film/substrate systems in [18]. In order to characterize surface morphological pattern formation of core-shell spheres upon the shrinkage of the core, Xu et al. introduced the relevant parameter $C_s = (E_s/E_f)(R/h_f)^{3/2}$, where $E_s$ and $E_f$ are, respectively, the Young's modulus for the substrate and the film. R the radius of the sphere ($R = 0.5$ mm, including the film thickness), and $h_f$ the thickness of the film ($h_f = 0.01$ mm).

Figure 3 is a plot of the mesh of one-eighth of the spherical film/substrate system (three symmetry planes are used to obtain the entire spherical structure). The three symmetry planes allow the computational time to be reduced, resulting in a restriction to only symmetric bulking modes. The tetrahedron geometry of the mesh and the presence of smaller elements within the film part of the mesh are apparent in the figure. The spherical film/substrate geometry is meshed with two elements in the film thickness, having an average element $h = 0.024$ mm, and it consists of 164,528 tetrahedra (689,358 degrees of freedom). We choose four different values of the parameter $C_s$ (0.2, 2.6, 4.9, 21.2); the mesh is the same for the four parameters.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Three-dimensional mesh of the film/substrate system with a spherical curvature. One-eighth of the whole structure is meshed, using symmetry planes. The mesh has been densified in the film thickness.</figcaption></figure>
## 3. Results

First, the FreeFEM++ numerical model is validated with a classical simulation example: a cantilever subjected to a vertical force. This numerical experiment gives satisfactory results and is presented in Appendix C. Second, in Section 3.1, the numerical results for planar film/substrate systems under uniaxial compression are presented. We took the opportunity to present impressive improvements in the convergence acceleration of ANM. The bifurcation curves and deformation profile of the film with respect to the loading forces are presented. Finally, in Section 3.2, we present bifurcation curves and deformation profile of the film in the case of a spherical film substrate subjected to a thermo-mechanical shrinkage of the substrate.

### 3.1. Planar Film/Substrate Systems Under Uniaxial Compression

The details of the geometry, materials, and the FEM interpolation are given in Section 2.5.1. According to the ANM algorithm, an ANM strategy with quite small steps has been chosen ($\delta = 10^{-6}$), and the ANM order is 15. In the case of the mesh described in Section 2.5.1, one ANM step requires almost 40 minutes of CPU time when the FreeFEM++ program is parallelized on 6 MPI processors while using the multifrontal MUMPS solver.

#### 3.1.1. Convergence Acceleration Improvements

Let us clarify that, according to recent works [37], acceleration convergence improvements of the ANM algorithm, including MMPE extrapolation, adaptive steps, and Newton-Riks correction, have been used. Let us take the opportunity in this subsection to present some important achievements regarding the convergence acceleration of the ANM algorithm: First, the total number of Newton's corrections according to the parameter $\delta$ is presented in Table 1. We can see that both acceleration convergence and step adaptation are very efficient in reducing Newton-Riks's corrections.

<table><caption>**Table 1.** Number of Newton corrections during 100 steps, according to the accuracy parameter $\delta$. The full algorithm includes convergence acceleration, step length adaptation, and Newton correction(s). Required accuracy $\epsilon_1 = 10^{-5}$. The first line reports the results of ANM is only followed by Newton corrections when necessary.</caption><thead><tr><th>δ</th><th>10<sup>-5</sup></th><th>10<sup>-6</sup></th><th>10<sup>-8</sup></th></tr></thead><tbody><tr><td>Pure Newton</td><td>152</td><td>87</td><td>53</td></tr><tr><td>Full algorithm</td><td>5</td><td>2</td><td>0</td></tr></tbody></table>

The residual before and after the application of the acceleration MMPE are presented in Table 2. Seven starting points $\mathbf{u}(ra_{max})$ with $r \in [0.7, 1.3]$ are considered. Step 34 has been chosen, with the full algorithm and $\delta = 10^{-6}$. We can notice again the great efficiency of the MMPE convergence acceleration in the normalized residual error.

<table><caption>**Table 2.** Convergence acceleration by MMPE: a typical example. The table presents the residual before and after the application of MMPE. Seven starting points $\mathbf{u}(ra_{max})$, $r \in [0.7, 1.3]$ are considered. Step 34 with the full algorithm and $\delta = 10^{-6}$.</caption><thead><tr><th>r</th><th>0.7</th><th>0.8</th><th>0.9</th><th>1</th><th>1.1</th><th>1.2</th><th>1.3</th></tr></thead><tbody><tr><td>before (×10<sup>-5</sup>)</td><td>0.9012</td><td>7.404</td><td>48.96</td><td>265.6</td><td>1224</td><td>4925</td><td>17,660</td></tr><tr><td>after (×10<sup>-6</sup>)</td><td>0.9346</td><td>0.9843</td><td>2.359</td><td>13.82</td><td>75.7</td><td>364.7</td><td>1584</td></tr></tbody></table>
3.1.2. Bifurcation Curves and Deformation Profile of the Film

The bifurcation diagram is shown in Figure 4, which represents the vertical displacement of the center of the film/substrate system (upper face) as a function of the load parameter $\lambda$. The accumulation of small steps in the diagram shows the presence of a bifurcation. Figure 5 shows the comparison between the bifurcation diagram obtained with ANM (full algorithm) and a classical Newton-Raphson continuation algorithm with Newton-Riks correction.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Bifurcation diagram: vertical displacement of the center of the film/substrate system (upper face) with $E_f/E_s = 0.72 \times 10^5$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Bifurcation diagram: vertical displacement of the center of the film/substrate system (upper face) with $E_f/E_s = 0.72 \times 10^5$: comparison ANM full algorithm and Newton-Raphson continuation algorithm with Newton-Riks-like corrections.</figcaption></figure>

Many bifurcations can be seen on the diagram of Figure 4. At the first bifurcation for $\lambda = 0.048$, uniformly distributed one-dimensional wrinkles appears. The left part of Figure 6 shows such periodic wrinkles appearing at the ANM step 20. Next, the curve
goes back for $\lambda = 0.09$, and five limit points are obtained that correspond to the increase in a single wrinkle near the boundary and to a change of the number of wrinkles [37]. The right part of Figure 6 shows the appearance of a localization phenomenon of the periodic wrinkles at ANM step 34. Small steps accumulate close to the limit points, which indicates that these singularities are, in fact, perturbed bifurcation points [39].

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Deformation of the film profile at the upper face (only the right half because of the symmetric plane). The color code is the displacement on the z-axis in mm. The left figure (ANM step 20, just before the first bifurcation), exhibits a periodic pattern with a few wavelengths, and the right figure (ANM step 34, just before the second bifurcation), exhibits a localization at the right edge superposed to the periodic pattern.</figcaption></figure>

## 3.2. Spherical Film/Substrate Under Thermal Loading

The details of the geometry, materials, and the FEM interpolation are given in Section 2.5.2. Let us first make the numerical simulation for $C_s = 0.2$. The top component of Figure 7 is the plot of the maximum normalized deflection versus the thermal loading parameter $\lambda$. Two limit or turning points can be observed, the bifurcation curve being subcritical: the first one is close to step 10, and the second is close to step 40. Then, close to step 100, another bifurcation happens, and the curve goes backward. The ANM continuation does not allow a jump to another probable existing branch: future research will be done to address this issue. The plots of the deformation of the film profile inserted in the figure at various ANM steps show isolated dimples, with a greater localization when the deflection amplitude increased.

Second, let us now consider the numerical simulation for $C_s = 2.6$. The bottom component of Figure 7 presents a plot of the maximum normalized deflection versus the thermal loading parameter $\lambda$. First, a limit or turning point can be observed, with the bifurcation being weakly subcritical. Small back-and-forth motions can be observed, and the ANM continuation still persists. The plots of the deformation of the film profile are inserted into the figure at various ANM steps, and they show buckyball deformation patterns.

Third, let us consider the case of $C_s = 4.9$. The top figure of Figure 8 is the plot of the maximum normalized deflection versus the thermal loading parameter $\lambda$. The first bifurcation appears close to step 10 and turns out to be supercritical. After, the curve shows a zigzag, and the progression slows down. The plots of the deformation of the film profile, which are inserted into the figure at various ANM steps, show a periodicity in the pattern, a transition between the buckyball, and a symmetry-breaking of the buckyball pattern around step 120 (the deflection amplitude of some dimples decreases).
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Bifurcation diagram (maximum normalized deflection at the upper face of the film versus the thermal loading parameter $\lambda$) for a spherical film/substrate system with a parameter $C_s = 0.2$ (defined in Section 2.5.2) in the **top** figure and the parameter $C_s = 2.6$ in the **bottom** figure. Plots of the deformation of the film profile at different steps are inserted into the figures.</figcaption></figure>

Finally, let us consider $C_s = 21.2$. The bottom component of Figure 8 is the plot of the maximum normalized deflection versus the thermal loading parameter $\lambda$. The first bifurcation is also supercritical as previously, and it appears close to step 10. Then, the bifurcation curve behavior is almost linear and reaches greater thermal loading values than in the previous case. The plots of the deformation of the film profile, which are inserted into the figure at various ANM steps, show a smaller periodicity in the pattern and also a transition between the buckyball and a symmetry breaking of the buckyball pattern.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Bifurcation diagram (maximum normalized deflection at the upper face of the film versus the thermal loading parameter $\lambda$) for a spherical film/substrate system with a parameter $C_s = 4.9$ (defined in Section 2.5.2) in the top figure and the parameter $C_s = 21.2$ in the bottom figure. Plots of the deformation of the film profile at different steps are inserted into the figures.</figcaption></figure>

**4. Discussion**

This paper has shed new light on the mathematical and algorithmic aspects of ANM based on M. Potier-Ferry's review article [22]. It has been shown that ANM traces its roots back to the implicit function theorem on a Banach space. Since ANM employs the Taylor-series expansion of both the unknown vectorial field $U$ and the scalar load $\lambda$, the chain rule for the derivation of function composition is a key point for deriving an efficient numerical algorithm for the ANM prediction part.

A new implementation of the ANM has been presented within the user-friendly finite element environment FreeFEM++. The ANM continuation method, based on a perturbation technique, has proven historically to be a very efficient technique for the numerical simulation of nonlinear instability problems. The FreeFEM++ implementation of ANM has shown many advantages compared to a previously implemented version
of ANM in FORTRAN, or in Matlab [27], and C++ [28], to solve various problems of
fluid or solid mechanics. Its high-level Design Specific Language dedicated to variational
formulations and finite elements makes it easier to write scripts for solving nonlinear
solid or fluid mechanics problems using ANM, and it is a key feature of this work. This
new numerical tool written in FreeFEM++ has been very helpful for new developments
like those concerning convergence acceleration [37]. Up to now, there are only a limited
number of applications of ANM for large-scale problems, generally for fluid flows ([29,40]).
The new implementation in FreeFEM++, which is linked to the library PETSc can be very
helpful for this purpose.

The ANM prediction, together with classical Newton-Riks corrections, has been shown
to be a very efficient tool to solve various nonlinear problems of fluid or solid mechanics. A
first computation, for validation purposes, using the newly developed FreeFEM++ code,
consisting in a cantilever (Saint Venant-Kirchhoff material) submitted to a vertical surface
force, in which geometrical nonlinearities are taken into account, has shown a very good
comparison with the commercial software Abaqus.

ANM was applied to two film-substrate systems. There were a few papers in the literature to analyze the post-bifurcation response of such systems. These authors have underlined the difficulty of these computations via a continuation procedure [18], and often, they switched to dynamical approaches. The first case is the study of a planar film/substrate system subjected to axial compression, and the second case is a spherical film/substrate system subjected to the thermo-mechanical shrinkage of its core. Both film/substrate numerical simulations have demonstrated the efficiency of the implementation of ANM in FreeFEM++ for the following of equilibrium paths when solving solid mechanical problems. In the case of the spherical film/substrate system, it has been possible to improve the following aspect of the equilibrium path for large thermal loading with the ANM compared to the previous work [18].

As emphasized in [41], the main strengths of the ANM are the high-order Taylor
expansion, which offers better interpolation than only a first differential order like the
classical Newton-Raphson method, and also the a posteriori estimation of the domain of
validity for the path parameter, which results in a consistent chaining of the ANM branches.

**5. Conclusions**

This paper has presented the mathematical and algorithmic basis of ANM supported
by two numerical experiments about film/substrate systems, which have illustrated the
efficiency of ANM to solve nonlinear solid mechanical problems. Of course, nonlinear
fluid mechanical problems can also be addressed through ANM. Since ANM continuation,
together with MMPE convergence acceleration, minimizes the number of required Newton-
Riks corrections, there is a lot of computational time saved compared to other continuation
techniques, while preserving very good accuracy.

In the field of perturbation theory, one has to remember Milton Van Dyke’s famous
sentence: “A reasonable number of terms of a perturbation series will reveal part of the
analytic structure of the solution”. It will be interesting, in the future, to better study the
properties of the Taylor series coefficients, for the purpose of having more information on
the bifurcation points, bifurcated branches, stability...

Also, recently, one of the co-authors of this paper (M. Potier-Ferry) claimed that it
would be worth looking for a separated equilibrium branch difficult to obtain with a
continuation algorithm, like ANM.

The numerical experiments presented in this paper have used a Saint Venant-Kirchhoff
law; however, ANM can also be easily implemented in the case of hyperelasticity or
plasticity [22].
Finally, following previous research works [42], the new implementation of ANM in the FreeFEM++ numerical development tool can handle large-scale problems, thanks to its interface with the PETSc library (multigrid, FETI, or domain decomposition techniques).

**Author Contributions:** Conceptualization, P.V., M.P.-F. and H.Z.; Methodology, P.V., M.P.-F., H.Z., F.X. and M.B.; Software, F.H. and A.-K.C.; Validation, P.V. and H.A.; Formal analysis, P.V., F.H., M.P.-F., M.B. and A.-K.C.; Investigation, P.V.; Data curation, H.A.; Writing—original draft, P.V. and M.P.-F.; Supervision, P.V.; Funding acquisition, M.P.-F. and F.X. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the French State through the program "Investment in the future" operated by the National Research Agency (ANR) and referenced by ANR-11-LABX-0008-01 (LabEx DAMAS). Fan Xu acknowledges financial support from the National Natural Science Foundation of China (Grants Nos. 12425204, 12122204, and 12372096), the Shanghai Pilot Program for Basic Research-Fudan University (Grant No. 21TQ1400100-21TQ010), and the Shanghai Municipal Education Commission (Grant No. 24KXZNA14).

**Data Availability Statement:** Data are contained within the article. The original contributions presented in this study are included in the article. Further inquiries can be directed to the corresponding author.

**Acknowledgments:** The authors thank the EXPLOR HPC center for computational resources.

**Conflicts of Interest:** The authors declare no conflicts of interest.

## Appendix A. The FreeFEM++ Implementation

FreeFEM++ is based on an efficient Domain-Specific Language (DSL), which allows partial differential equations and, more generally, nonlinear multiphysics systems in 2D, 3D, and the surface to be solved. The FreeFEM++ language is typed, polymorphic, and reentrant with macro-generation. The kernel of FreeFEM++ is written in C++.

The following are the main characteristics of FreeFEM++:

1. FreeFEM++ involves a high-level, user-friendly typed input language with an algebra of analytic and finite element functions.

2. The finite element problems are described according to their variational formulations, and they are easily implemented. It is possible to gain access to the internal vectors and matrices if needed.

3. Many kinds of problems can be solved: multi-variables, multi-equations, and one-, two-, or three-dimensional static or time-dependent, linear, or nonlinear coupled systems.

4. An automatic mesh generator is available, based on the Delaunay-Voronoï algorithm.

5. A large variety of triangular finite elements: linear and quadratic Lagrangian finite elements and more, discontinuous P1, and Raviart-Thomas elements.

6. A large variety of linear direct and iterative solvers (LU, Cholesky, Crout, CG, GMRES, UMFPACK, MUMPS, SuperLU, and ...), eigenvalue, and eigenvector solvers (ARPACK) are available.

Most of these characteristics are of great interest for simulating nonlinear mechanical problems using the ANM continuation algorithm. The following demonstrates some of the advantages of code implementation with FreeFEM++.

1. The mesh generator is very easy to use for simple geometries.

2. The macro-generation enables easy implementation of the differential operators useful in nonlinear solid mechanics.

3. The variational formulation oriented commands in FreeFEM++ allow us to implement the weak form of nonlinear solid mechanical differential equations very efficiently, discretized using the Finite Element Method.
The main steps of the Finite Element Method discretization in the FreeFEM++ implementation are presented next. They begin with the Laplacian problem, a basic problem for applied mathematics.

We would like to find $u$, a real function defined on a domain, $\Omega$, of $\mathbb{R}^3$, such that

$$ egin{cases} -\Delta u = f \text{ in } \Omega \\ u = 0 \text{ on } \Gamma_D \\ \frac{\partial u}{\partial n} = g \text{ on } \Gamma_N \end{cases} \tag{A1} $$

with $\Gamma_D \cup \Gamma_N = \partial\Omega$. The boundary condition on $\Gamma_D$ (resp. on $\Gamma_N$) is called Dirichlet (resp. Neumann) boundary conditions.

Let us define the functional space, $V_0$, as follows:

$$ V_0(\Omega) = \{v \in H^1(\Omega), v|_{\Gamma_D} = 0\} \tag{A2} $$

Using Green's formula that generalizes the integration by part, an equivalent weak form (or variational form) of the previous problem (A1) is obtained: Find $u \in V_0(\Omega)$, such that

$$ \int_{\Omega} \nabla u \nabla v d\Omega = \int_{\Omega} f v d\Omega + \int_{\Gamma_D} v \frac{\partial u}{\partial n} ds, \quad \forall v \in V_0(\Omega) \tag{A3} $$

The Finite Element Method starts with the partitioning of the initial domain, $\Omega$, into subdomains, $\Omega_e$, called elements. It results in the spatial discretization of the domain (mesh). In FreeFEM++, the 2D (resp. 3D) mesh, is partitioned into triangles (resp. tetrahedra). FEM approximates all functions, $w$, as follows:

$$ w(x, y, z) = \sum_{i=1}^{N} w_i \phi_i(x, y, z) \tag{A4} $$

where $N$ is the number of nodes of the mesh, depending on the FEM interpolation, detailed below. The restriction of the FEM basis functions, $\phi_k(x, y, z)$, to each tetrahedron is called the shape function.

The finite element space, $V_h$, is defined as follows:

$$ V_h = \left\{ w = \sum_{i=1}^{N} w_i \phi_i, w_i \in \mathbb{R} \right\} \tag{A5} $$

The list of finite elements available in the FreeFEM++ environment is wide: $P_0$ (piecewise constant discontinuous finite element), $P_1$ (piecewise linear continuous finite element), $P_{1dc}$ (piecewise linear discontinuous finite element), $P_2$ (piecewise $P_2$ continuous finite element), ...

First, let us focus on the elementary Laplacian problem (A1). We will assume that the 3D mesh has already been generated and saved in the file "Th3D.msh". Let us load the mesh and build the $P_1$ finite element interpolation space on this mesh. This is accomplished in FreeFEM++ with the following two command lines:

```
msh3 Th("Th3D.msh ");
fespace Vh(Th,P1);
```

Then, the differential operator gradient will be defined with the following command macro:

```
macro Grad(u) [dx(u), dy(u), dz(u)] //
```
Let us define $u$, the solution function, and $v$, the test function, belonging to the finite element space $V_h$, and implement the equivalent variational formulation of the Laplacian problem:

$$
\begin{array}{l}
V_h u, v; \\
\\
\text{problem Poisson}(u, v) = \\
\quad \text{int3d}(\text{Th}) (\text{Grad}(u)' * \text{Grad}(v)) - \text{int3d}(\text{Th}) (f * v) - \text{int2d}(\text{Th}, \text{LN}) (g * v) \\
\quad + \text{on}(\text{LD}, u=0);
\end{array}
$$

The first (resp. second) operator int3d is used for the bilinear part (resp. linear part) of
the variational formulation (A3). The operator int2d allows the Neumann boundary condi-
tions (LN is the label of the mesh concerned with the Neumann boundary conditions) to be
taken into account, and the operator on(.) is used to describe the Dirichlet boundary condi-
tions (here, LD is the label of the mesh concerned with the Dirichlet boundary conditions).

The numerical solution of the problem (A3) is obtained using the following command:

Poisson;

It is also possible to plot the solution field *u* using the following command:

plot(u, wait = 1);

It is easy to generalize the FreeFEM++ implementation of the Laplacian problem
to the resolution of a linear or nonlinear elastic problem. Well-documented FreeFEM++
implementations can be found on the FreeFEM++ website: https://modules.freefem.
org/modules/elasticity/ and https://modules.freefem.org/modules/nonlinear-elasticity/
(both accessed on 1 September 2025).

FreeFEM++ enables both the sparse Finite Element matrices and the right-hand sides to be generated directly. Both are useful for the implementation of the ANM algorithm.

Let us now describe the main steps for the implementation of the resolution of a
mechanical problem with geometrical nonlinearities using the ANM algorithm in the
FreeFEM++ language. We will assume that the material is elastic with a Saint Venant-
Kirchhoff constitutive law.

The FEM stiffness matrix [K<sub>f</sub>(u<sub>0</sub>, S<sub>0</sub>)] is the implementation of the bilinear form (11),
the FEM right-hand side vector {F} is the implementation of the linear form P<sub>e</sub>(δu) (12),
and the FEM right-hand side vector {F<sub>p</sub><sup>nl</sup>} is the implementation of the linear form
⟨F<sub>p</sub><sup>nl</sup>, δu⟩ (15).

Differential operators, needed to define Green–Lagrange tensor and its variation, are
easily implemented in the FreeFEM++ language using a command macro. It represents
one of the strengths of FreeFEM++ for variational formulation implementation. For this
purpose, let us present the main macro commands:

macro GammaL(u,v,w) [dx(u),dy(v),dz(w),(dy(u)+dx(v)),
(dz(u)+dx(w)),(dz(v)+dy(w))] //

macro GammaNL(u1,v1,w1,u2,v2,w2)
[ (dx(u1)*dx(u2)+dx(v1)*dx(v2)+dx(w1)*dx(w2))*0.5 ,
( dy(u1)*dy(u2)+dy(v1)*dy(v2)+dy(w1)*dy(w2))*0.5 ,
( dz(u1)*dz(u2)+dz(v1)*dz(v2)+dz(w1)*dz(w2))*0.5 ,
( dy(u1)*dx(u2)+dx(u1)*dy(u2)+dy(v1)*dx(v2)
+dx(v1)*dy(v2)+dy(w1)*dx(w2)+dx(w1)*dy(w2))*0.5 ,
( dz(u1)*dx(u2)+dx(u1)*dz(u2)+dz(v1)*dx(v2)
+dx(v1)*dz(v2)+dz(w1)*dx(w2)+dx(w1)*dz(w2))*0.5 ,
( dz(u1)*dy(u2)+dy(u1)*dz(u2)+dz(v1)*dy(v2)
+dy(v1)*dz(v2)+dz(w1)*dy(w2)+dy(w1)*dz(w2))*0.5 ] //
```c
macro Gamma(u,v,w) (GammaL(u,v,w)+GammaNL(u,v,w,u,v,w)) //
macro dGammaNL(u,v,w,uu,vv,ww) (2.0*GammaNL(u,v,w,uu,vv,ww)) //
macro dGamma(u,v,w,uu,vv,ww) (GammaL(uu,vv,ww)+
                                dGammaNL(u,v,w,uu,vv,ww)) //
```

Here, the classical Voigt notation has been used for the Green-Lagrange tensor.

GammaL(u,v,w) represents the linear part of the Green-Lagrange tensor, and
GammaNL(u1,v1,w1,u2,v2,w2)) represents its nonlinear part (4).

The tensor dGammaNL(u,v,w,uu,vv,ww) is the variation in the non-linear part of the Green-Lagrange tensor (5).

Then, the elasticity matrix is easily implemented:

```c
real E = 1.e5; nu = 0.;
real lambda = E*nu/(1.0+nu)/(1.0-2.0*nu); mu = E/2.0/(1.0+nu);
func D = [
    (lambda+2.0*mu), lambda, lambda, 0.0, 0.0, 0.0 ],
    [ lambda, (lambda+2.0*mu), lambda, 0.0, 0.0, 0.0 ],
    [ lambda, lambda, (lambda+2.0*mu), 0.0, 0.0, 0.0 ],
    [ 0.0, 0.0, 0.0, mu, 0.0, 0.0 ],
    [ 0.0, 0.0, 0.0, 0.0, mu, 0.0 ],
    [ 0.0, 0.0, 0.0, 0.0, 0.0, mu ]
];
```

E, nu are Young's modulus and the Poisson ratio, lambda and mu are Lamé's coefficients, and D is the elasticity matrix.

Let us now describe how to implement in the language FreeFEM++ the variational formulation that will lead to the computation of both the stiffness matrix and the right-hand-side vector.

P2 Lagrange FEM interpolation is assumed and defined on the mesh Th3D via the fespace command:

```c
fespace Vh(Th3D, [P2,P2,P2]);
```

Let us declare the displacement fields for each ANM order:

```c
Vh[int] [u,v,w](Norder+1);
```

Here, we will not use the previous command problem but the key command varf, which facilitates the computation of the tangent matrix and the linear right-hand side while taking into account the Dirichlet boundary conditions (using the keyword on).

```c
varf PbTg ([u1,v1,w1],[uu,vv,ww]) =
    int3d (Th3D) ( (dGamma(u[0],v[0],w[0],uu,vv,ww))'*
    (D*(GammaL(u1,v1,w1)+2*GammaNL(u[0],v[0],w[0],u1,v1,w1)))
    +(dGammaNL(u1,v1,w1,uu,vv,ww))'* (D*(GammaL(u[0],v[0],w[0])
    +GammaNL(u[0],v[0],w[0],u[0],v[0],w[0])) )
    +on(lencasmid,u1=0.,v1=0.,w1=0.);
varf PbF ([u1,v1,w1],[uu,vv,ww]) = int2d (Th3D, lrigthmid) (Pa*ww)
    +on(lencasmid,u1=0.,v1=0.,w1=0.);
```

lrigthmid (resp. lencasmid) is the label of the mesh concerned with the vertical pressure (resp. Dirichlet boundary conditions).

The first varf PbTg is used to assemble the tangent matrix:

```c
matrix Kt = PbTg(Vh, Vh, solver(sparsesolver));
```

The second varf PbF is used to assemble the nodal vector of the linear right-hand side:
Vh [Fu, Fv, Fw];
Fu[] = PbF(0, Vh);

Let us now mention the important points appearing in the loop on the ANM orders for one ANM step.

First, the finite element displacement field {$\hat{u}$} is computed using a linear system resolution (we used the MUMPS direct solver, which has proven its efficiency).

Vh [u1hat, v1hat, w1hat];
u1hat[] = Kt^-1*Fu[];

It is now necessary to assemble the nonlinear right-hand-side vector {$F_2^{nl}$}. As is done in the ANM Algorithm 1, we define intermediate tensors (written in a vectorial form) {$S_1$} and {$S_2^{nl}$} defined at the Gauss points:

fespace QFh6(ThL3D,
[FEQF53d, FEQF53d, FEQF53d, FEQF53d, FEQF53d, FEQF53d]);
QFh6[int] [Sxx, Syy, Szz, deuxSxy, deuxSxz, deuxSyz](Norder+1);
QFh6[Sn1xx, Sn1yy, Sn1zz, deuxSnlxy, deuxSnlxz, deuxSnlyz](Norder+1);

The first step is to compute {$S_1$}:

[Sxx[1], Syy[1], Szz[1], deuxSxy[1], deuxSxz[1], deuxSyz[1]] =
DL*(GammaL(u[1], v[1], w[1]) +
2*GammaNL(u[0], v[0], w[0], u[1], v[1], w[1]));

Then, we obtain {$S_2^{nl}$}:

[Sn1xx, Sn1yy, Sn1zz, deuxSnlxy, deuxSnlxz, deuxSnlyz] =
DL*(GammaNL(u[1], v[1], w[1], u[1], v[1], w[1]));

It is now possible to assemble the nonlinear second member, {$F_2^{nl}$}:

varf PbFn12 ([u1,v1,w1],[uuu,vvv,www]) = - int3d(Th3D) (
(dGammaNL(u[1],v[1],w[1],uuu,vvv,www))'*
([Sxx[1],Syy[1],Szz[1],deuxSxy[1],deuxSxz[1],deuxSyz[1]])
+ (dGamma(u[0],v[0],w[0],uuu,vvv,www))'*
([Sn1xx,Sn1yy,Sn1zz,deuxSnlxy,deuxSnlxz,deuxSnlyz]))
+ on(lencasmid,u1=0.,v1=0.,w1=0.);
Fnlu[] = PbFn12(0,Vh);

Let us now enter the loop with respect to the order $p \ge 2$.

Let us compute {$\mathbf{u}_p$}:

unl[] = Kt^-1*Fnlu[];
lambda[p] = -lambda[1] * (unl[]'*u[1][]) ;
u[p][] = lambda[p]/lambda[1]*u[1][] + unl[];

Then, we compute {$S_p$}:

[Sxx[p], Syy[p], Szz[p], deuxSxy[p], deuxSxz[p], deuxSyz[p]] =
DL*(GammaL(u[p], v[p], w[p]) +
2*GammaNL(u[0], v[0], w[0], u[p], v[p], w[p]))
+[Sn1xx, Sn1yy, Sn1zz, deuxSnlxy, deuxSnlxz, deuxSnlyz];

a loop is necessary to compute {$S_{p+1}^{nl}$} for the next ANM order;

[Sn1xx, Sn1yy, Sn1zz, deuxSnlxy, deuxSnlxz, deuxSnlyz] =
[0., 0., 0., 0., 0., 0.];
```c
for (int ir = 1; ir < (p + 1); ir++)
{
    [Snlxx, Snlyy, Snlzz, deuxSnlxy, deuxSnlxz, deuxSnlyz] =
        [Snlxx, Snlyy, Snlzz, deuxSnlxy, deuxSnlxz, deuxSnlyz]
        + GammaNL(u[p+1-ir], v[p+1-ir], w[p+1-ir], u[ir], v[ir], w[ir]);
}

[Snlxx, Snlyy, Snlzz, deuxSnlxy, deuxSnlxz, deuxSnlyz] = D*
[Snlxx, Snlyy, Snlzz, deuxSnlxy, deuxSnlxz, deuxSnlyz];

Let us notice that {S<sup>nl</sup><sub>p</sub>} is overwritten at each order.
It is now possible to assemble the nonlinear second member, {F<sup>nl</sup><sub>p+1</sub>}:

Vh [Fnlutmp, Fnlvtmp, Fnlwtmp];
    [Fnlu, Fnlv, Fnlw] = [0., 0., 0.];
    [Fnlutmp, Fnlvtmp, Fnlwtmp] = [0., 0., 0.];

for (int ir = 1; ir < (p + 1); ir++)
{
    varf PbFnla ([utmp, vtmp, wtmp], [uuu, vvv, www]) = - int3d (Th3D)
        (dGammaNL(u[p+1-ir], v[p+1-ir], w[p+1-ir], uuu, vvv, www))';
    ([Sxx[ir], Syy[ir], Szz[ir], deuxSxy[ir], deuxSxz[ir], deuxSyz[ir])
        )
        + on (lencasmid, utmp = 0., vtmp = 0., wtmp = 0.);
    Fnlutmp[] = PbFnla(0, Vh);
    Fnlu[] = Fnlu[] + Fnlutmp[];
}

Vh [Fnlutmp, Fnlvtmp, Fnlwtmp];
Fnlutmp[] = 0.0;
varf PbFnlc ([utmp, vtmp, wtmp], [uuu, vvv, www]) =
    - int3d (Th3D) (dGamma(u[0], v[0], w[0], uuu, vvv, www))';
    ([Snlxx, Snlyy, Snlzz, deuxSnlxy, deuxSnlxz, deuxSnlyz])
        )
        + on (lencasmid, utmp = 0., vtmp = 0., wtmp = 0.);
    [Fnlutmp, Fnlvtmp, Fnlwtmp] = [0., 0., 0.];
    Fnlutmp[] = PbFnlc(0, Vh);
    Fnlu[] = Fnlu[] + Fnlutmp[];
```

Let us now compute the norm

$L_2$ of $u_1$ and $u_N$ to determine $a_{max}$ (18):

```c
real Norm1 = sqrt(int3d(Th3D) (u[1]'*u[1]+v[1]'*v[1]+w[1]'*w[1])));
real NormNo = sqrt(int3d(Th3D) (u[Norder]'*u[Norder]+
                                v[Norder]'*v[Norder]+
                                w[Norder]'*w[Norder])));
amax = (delta*Norm1/NormNo)^(1.0/(Norder-1.0));
```

The final task consists of evaluating the Taylor series of the displacement nodal vector, computing the residual error, and actualizing both $u_0$ and $\lambda_0$.

The FreeFEM++ scripts for the FEM simulation of non-linear solid mechanics using ANM will be available soon on the FreeFEM++ website (https://freefem.org/).

## Appendix B. Newton-Raphson Continuation and Newton-Riks Correction Algorithms

Algorithm A1 (resp. Algorithm A2) presents the several steps of the classical Newton-Raphson (resp. Newton-Riks correction) algorithm.
**Algorithm A1 Algorithm: Newton-Raphson prediction iteration**

1: **Initialize** $\Delta s$ (ARC LENGTH)

2: **Initialize** $\{U_0\}$ AND, $\{S_0\}$ BEFORE NEWTON-RAHPSON PREDICTION

3: **Assemble** $[K_t(U_0, S_0)]$ AND $\{F\}$

4: $\{\Delta\hat{U}\} \leftarrow [K_t(U_0, S_0)]^{-1}\{F\}$

5: $\Delta\lambda_0 \leftarrow \frac{\pm\Delta s}{\sqrt{1+\{\Delta\hat{U}\}^2}}$

6: $\{\Delta U_0\} \leftarrow \Delta\lambda_0\{\Delta\hat{U}\}$

**Algorithm A2 Algorithm: Newton-Riks correction iteration**

1: **Initialize** $U, S$ AND $\lambda$ BEFORE NEWTON-RIKS CORRECTIONS

2: **Initialize** $\Delta U_0$ AND $\Delta\lambda_0$ THE PREDICTION INCREMENT

3: **Assemble** $[K_t(U, S)]$ AND $\{F\}$

4: **Compute** THE RESIDUAL VECTOR $\{R\}$

5: $\{\Delta U_R\} \leftarrow -[K_t(U, S)]^{-1}\{R\}$

6: $\{\Delta U_F\} \leftarrow [K_t(U, S)]^{-1}\{F\}$

7: $\Delta\lambda_c \leftarrow -\frac{\{\{\Delta U_0\}, \{\Delta U_R\}\}}{\{\{\Delta U_0\}, \{\Delta U_F\}\} + \Delta\lambda_0}$

8: $\{\Delta U_c\} \leftarrow \{\Delta U_R\} + \Delta\lambda_c\{\Delta U_F\}$

**Appendix C. Clamped Beam Subjected to Conservative Vertical Surface Traction**

Let us consider a clamped beam subjected to vertical conservative surface traction. This problem has already been studied by Zahrouni [25] in the case of shell finite elements and a concentrated force. For validation purposes, we will add some simulation results obtained with the commercial software Abaqus.

Consider a clamped beam subjected to a conservative load at its free end. The length of the beam is 10 mm, its width is 1 mm, and the thickness is 1 mm. The left face of the beam is clamped, and conservative vertical traction is applied to the left face (see Figure A1).

<figure><img src="image_11.png" /><figcaption>**Figure A1.** Clamped beam submitted to a vertical conservative surface traction.</figcaption></figure>

The beam has been meshed using tetrahedral P2 Lagrangian finite elements (see Figure A2).

The beam is supposed to be a homogeneous isotropic elastic material with a Young's modulus of $10^5$ MPa and a 0 Poisson's ratio. A FreeFEM++ script has been developed in order to simulate the non-linear elastic behavior of the clamped beam. To build the finite element mesh of the beam, we chose 40 elements for their length and 4 elements for their width and height. The ANM order is 15, and the ANM branch length is determined by $\delta = 10^{-8}$. The amplitude of the surface traction is $t_{0z} = 1$ MPa. Let us recall that the loading parameter $\lambda$ appears in the weak form related to the applied force (the right term
of the second equation of (7)). Figure A3 shows the plot of the load parameter $\lambda$ versus the vertical displacement of the center of the right face. We can see that FreeFEM++ and Abaqus computations are in very good agreement. Large deflection up to $u_z = 8.5\text{mm}$ were obtained with only 7 ANM steps. At the end of each step, the normalized residual error is close to $10^{-6}$. Figure A3 shows the plot of the FreeFEM++ computations for the ANM orders 15 and 30. We can notice that, for the ANM order 30, the step length is greater than for the ANM order 15. Only 4 steps are needed for the whole curve when the ANM order is equal to 30.

<figure><img src="image_12.png" /><figcaption>Figure A2. Mesh of the beam using tetrahedral P2 Lagrangian finite elements.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure A3. Plot of $\lambda$ versus the vertical displacement of the center of the right face. Dots mark the end of ANM steps.</figcaption></figure>

When we choose $\delta = 10^{-5}$, the step length increases, but the residual error increases roughly to $10^{-3}$.

## References

1. Li, B.; Cao, Y.P.; Feng, X.Q.; Gao, H. Mechanics of morphological instabilities and surface wrinkling in soft materials: A review. *Soft Matter* **2012**, *8*, 5728–5745. [CrossRef]

2. Cao, Y.; Hutchinson, J. Wrinkling phenomena in neo-Hookean film/substrate bilayers. *Appl. Mech. 2012*, *79*, 031019-1-031019-9. [CrossRef]

3. Sun, J.Y.; Xia, S.; Moon, M.-W.; Oh, K.H.; Kim, K.S. Folding wrinkles of a thin stiff layer on a soft substrate. *Proc. R. Soc. A Math. Phys. Eng. Sci.* **2012**, *468*, 932–954. [CrossRef]
4. Cao, Y.; Hutchinson, J. From wrinkles to creases in elastometers: The instability and imperfection-sensitivity wrinkling. *Proc. R. Soc. A Math. Phys. Eng. Sci.* **2012**, *468*, 94–115. [CrossRef]

5. Allen, H. *Analysis and Design of Structural Sandwich Panels*; Elsevier: Amsterdam, The Netherlands, 1969.

6. Lourie, O.; Cox, D.; Wagner, H. Buckling and collapse of embedded carbon nanotubes. *Phys. Rev. Lett.* **1998**, *81*, 1638–1641. [CrossRef]

7. Rogers, J.; Someya, T.; Huang, Y. Materials and mechanics for stretchable electronics. *Science* **2010**, *37*, 207–213. [CrossRef] [PubMed]

8. Efimenko, K.; Rackaitis, M.; Manias, E.; Vaziri, A.; Mahadevan, L.; Genzer, J. Nested self-similar wrinkling patterns in skin. *Nat. Mater.* **2005**, *4*, 293–297. [CrossRef] [PubMed]

9. Budday, S.; Steinmann, P.; Kuhl, E. The role of mechanics during brain development. *J. Mech. Phys. Solids* **2014**, *72*, 75–92. [CrossRef]

10. Schweikart, A.; Fery, A. Controlled wrinkling as a novel method for the fabrication of patterned surfaces. *Microchim. Acta* **2009**, *165*, 249–263. [CrossRef]

11. Lin, J.; Zhong, J.; Zhong, S.; Li, H.; Zhang, H.; Chen, W. Modulating electronic transport properties of MoS2 field effect transistor by surface overlayers. *Appl. Phys. Lett.* **2013**, *103*, 063109. [CrossRef]

12. Bowden, N.; Brittain, S.; Evans, E.; Hutchinson, J.; Whitesides, G. Spontaneous formation of ordered structures in thin films of metals supported on an elastomeric polymer. *Nature* **1998**, *393*, 146–149. [CrossRef]

13. Bowden, N.; Huck, W.; Paul, K.; Whitesides, G.M. The control formation of ordered, sinusoidal structures by plasma oxidation of an elastometric polymer. *Appl. Phys. Lett.* **1999**, *75*, 2557–2559. [CrossRef]

14. Wang, S.; Song, J.; Kim, D.H.; Huang, Y.; Rogers, J.A. Local versus global buckling of thin films on elastometric substrates. *Appl. Phys. Lett.* **2008**, *93*, 023126. [CrossRef]

15. Howarter, J.; Stafford, C. Instabilities as a measurement tool for soft materials. *Soft Mater* **2010**, *6*, 5661–5666. [CrossRef]

16. Chen, X.; Hutchinson, J. Herringbone buckling patterns of compressed thin films on compliant substrates. *J. Appl. Mech.* **2004**, *71*, 597–603. [CrossRef]

17. Wong, W.; Pellegrino, S. Wrinkled membranes III: Numerical simulations. *J. Mech. Mater. Struct.* **2006**, *1*, 63–95. [CrossRef]

18. Xu, F.; Zhao, S.; Conghua, L.; Potier-Ferry, M. Pattern selection in core-shell spheres. *J. Mech. Phys. Solids* **2020**, *137*, 103892. [CrossRef]

19. Rodriguez, J.; Rio, G.; Cadou, J.M.; Troufflard, J. Numerical study of dynamic relaxation with kinetic damping applied to inflatable fabric structures with extensions for 3D solid element and non-linear behavior. *Thin-Walled Struct.* **2011**, *49*, 1468–1474. [CrossRef]

20. Taylor, M.; Bertoldi, K.; Steigmann, D.J. Spatial resolution of wrinkle patterns in thin elastic sheets at finite strain. *J. Mech. Phys. Solids* **2014**, *62*, 163–180. [CrossRef]

21. Cochelin, B.; Damil, N.; Potier-Ferry, M. *Méthode Asymptotique Numérique*; Hermès-Lavoisier: Paris, France, 2007.

22. Potier-Ferry, M. Asymptotic numerical method for hyperelasticity and elastoplasticity: A review. *Proc. R. Soc.* **2024**, *480*, 20230714. [CrossRef]

23. Thompson, J.M.T.; Walker, A.C. The non-linear perturbation analysis of discrete structural systems. *Int. J. Solids Struct.* **1968**, *4*, 757–768. [CrossRef]

24. Cochelin, B. A path following technique via an asymptotic-numerical method. *Comput. Struct.* **1994**, *53*, 1181–1192. [CrossRef]

25. Zahrouni, H.; Cochelin, B.; Potier-Ferry, M. Computing finite rotations of shells by an asymptotic-numerical method. *Comput. Methods Appl. Mech. Eng.* **1999**, *175*, 71–85. [CrossRef]

26. Riks, E. An incremental approach to the solution of snapping and buckling problems. *Int. J. Solids Struct.* **1979**, *15*, 529–551. [CrossRef]

27. R. Arquier, B. Cochelin, C. Manlab. Available online: http://manlab.lma.cnrs-mrs.fr/ (accessed on 1 September 2025).

28. Lejeune, A.; Béchet, F.; Boudaoud, H.; Mathieu, N.; Potier-Ferry, M. Object-oriented design to automate a high order non-linear solver based on asymptotic numerical method. *Adv. Eng. Softw.* **2012**, *48*, 70–88. [CrossRef]

29. Guevel, Y.; Allain, T.; Girault, G.; Cadou, J.M. Numerical bifurcation analysis for 3-dimensional sudden expansion fluid dynamic problem. *Int. J. Numer. Methods Fluids* **2018**, *87*, 1–26. [CrossRef]

30. Hecht, F. New development in FreeFem++. *J. Numer. Math.* **2012**, *20*, 251–265. [CrossRef]

31. Gourdon, X. *Analyse*; Ellipse: Boulogne-Billancourt, France, 2020.

32. Berger, M.S. *Nonlinearity and Functional Analysis: Lectures on Nonlinear Problems in Mathematical Analysis*; Academic Press: New York, NY, USA, 1977; Volume 74.

33. Berger, M.S. Method of Lyapunov and Schmidt in the study of non-linear equation and their further development. *Russ. Math. Surv.* **1962**, *17*, 1–60.

34. Cochelin, B.; Damil, N.; Potier-Ferry, M. The asymptotic numerical method: An efficient technique for non linear structural mechanics. *Rev. Eur. Eléments Finis* **1994**, *3*, 281–297. [CrossRef]
35. Holzapfel, G. *Nonlinear Solid Mechanics: A Continuum Approach for Engineering*; Wiley: Hoboken, NJ, USA, 1969.

36. Mottaqui, H.; Braikat, B.; Damil, N. Discussion about parameterization in the asymptotic numerical method: Application to nonlinear elastic shells. *Comput. Methods Appl. Mech. Eng.* **2010**, *199*, 1701–1709. [CrossRef]

37. Ventura, P.; Potier-Ferry, M.; Zahrouni, H. A secure version of asymptotic numerical method via convergence acceleration. C. R. **Mécaniques** **2020**, *348*, 361–374. [CrossRef]

38. Büchter, N.; Ramm, E.; Roehl, D. Three dimensional extension of non-linear shell formulation based on the enhanced assumed strain concept. *Int. J. Numer. Methods Eng.* **1994**, *37*, 2551–2568. [CrossRef]

39. Baguet, S.; Cochelin, B. On the behaviour of the ANM continuation in the presence of bifurcation. *Commun. Numer. Methods Eng.* **2003**, *19*, 459–471. [CrossRef]

40. Medale, M.; Cochelin, B. A parallel computer implementation of the asymptotic numerical method to study thermal convection instabilities. *J. Comput. Phys.* **2009**, *228*, 8249–8262. [CrossRef]

41. Cochelin, B.; Medale, M. Power series analysis as a major breakthrough to improve the efficiency of Asymptotic Numerical Method in the vicinity of bifurcations. *J. Comput. Phys.* **2013**, *236*, 594–607. [CrossRef]

42. Galliet, I.; Cochelin, B. Une version parallèle des MAN par décomposition de domaine. *Rev. Eur. Éléments Finis* **2012**, *13*, 177–195. [CrossRef]

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.