Nonlinear Analysis of Wrinkles in Film-Substrate Systems by Finite Element Method and Asymptotic Numerical Method

1. Introduction

The wrinkle phenomena of thin membranes [1] or film-substrate systems [2] are very common in the recent literature. The numerical modeling of these instabilities may be highly complicated due to the existence of many equilibrium solutions and the poor stiffness of unloaded membranes. In such a situation, it is not easy to apply a continuation algorithm that may fails [1,3] or undergoes back-and-forth loops [4]. Such a difficulty may be avoided by using a pseudo-dynamic algorithm, e.g. [3,5,6]. In order to assess the limit of validity of continuation techniques in the quasi-static regime, we consider in this paper ANM with applications to film-substrate systems.

ANM is a continuation method based on Taylor series. It is associated with classical discretization techniques, generally the Finite Element Method (FEM). The coupling between a perturbation technique and FEM was first proposed by Thompson and Walker [7]. In [8,9,10], a continuation method was introduced together with an efficient procedure to build the Taylor series and the convergence acceleration using Padé approximants. Then ANM has been successfully applied to various engineering and scientific problems in solid and fluid mechanics [11,12,13,14,15,16,17]. It is rather different from classical incremental-iterative methods such as Riks strategy [18], because full solution paths are computed and not only 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, C++. Many home-made implementations were established. For example, the open source software MANLAB [19] is widely used in the field of nonlinear vibrations. Lejeune et al. [20] 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 enough user-friendly. Application in fluid mechanics within the open source finite element platform ELMER has been done and results were presented in [21].

FreeFEM++ [22] is a very efficient environment to implement variational formulations. As we will show in this paper, FreeFEM++ uses a high level language that allows to derive very easily variational formulation and to implement the ANM algorithm. The main idea is to develop a general, and easy to use, environment to implement the ANM algorithm for solving solid mechanical problems. FreeFEM++ also offers MPI parallel solvers and multigrid or domain decomposition tools that enables to reduce the CPU time and gives the opportunity to take into account large number of degrees of freedom. In this paper, the implementation of ANM in the FreeFEM++ code is discussed. Let us mention a publication regarding the analysis of infinite periodic surface acoustic waves transducers using FreeFEM++ [23] and a recent publication about fluid-structure interaction problems solved using FreeFEM++ [24].

The second section will present the main mathematical aspects of ANM. The third section the application of ANM in the case of nonlinear solid elastic problems. And the fourth section is dedicated to thermo-elastic nonlinear problems. The fifth section we will show what makes FreeFEM++ an efficient tool to implement ANM for various nonlinear mechanical problems, and will describe the implementation of the ANM algorithm in FreeFEM++. The sixth section presents the simulation results, using the FreeFEM++ implementation of ANM, of the nonlinear mechanical behavior of a clamped cantilever and subjected to a vertical surface force at its right surface, of a planar film/substrate system under uni-axial loading [25,26], and of a spherical film/substrate system subjected to a thermo-mechanical shrinkage of its core [3].

2. Theoretical Aspects of ANM

Let us recall the mathematical aspects of the Asymptotic Numerical Method (ANM) [27][28]. Let us begin with the resolution of the equation (1) in the case of a finite n-dimensional vectorial space:

$$\mathbf{R}\left(\mathbf{U},\lambda \right)=0$$

(1)

with, $\mathbf{R}:\left({\mathbb{R}}^n\times \mathbb{R}\right)⟶{\mathbb{R}}^n$, $\mathbf{U}\in {\mathbb{R}}^n$, and $\lambda \in \mathbb{R}$

Generally, the problem (1) is formulated for $\mathbf{U}$ belonging to a Banach space, but here we will considered problems coming from engineering, and physics, formulated in a discretized space using a numerical method like the Finite Element Method or using only a finite number of degrees of freedom.

Let us now assume that the function $\mathbf{R}\left(\mathbf{U},\lambda \right)$ is a class $C^p$ function (with $p>=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 that for $\left({\mathbf{U}}_0,{\lambda }_0\right)\in H\times I$, the partial differential ${\displaystyle \frac{\partial \mathbf{R}}{\partial \mathbf{U}}}{|}_{\left({\mathbf{U}}_0,{\lambda }_0\right)}$ is a bi-continuous isomorphism. The implicit function theorem [29] states that the equation (1) has a unique solution branch $\varphi :\lambda ⟶\mathbf{U}=\varphi \left(\lambda \right)$ 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}\left(\mathbf{U},\lambda \right)$ is an analytic function, the function $\varphi \left(\lambda \right)$ is also analytic and it is possible to expand both $\left(\mathbf{U}\right)$ and $\lambda$ into serial using the parameter a called path parameter within a given radius of convergence [28][30][31]:

$$\left\{\begin{array}{cc}\hfill \mathbf{U}\left(a\right)& ={\mathbf{U}}_0+\sum _{i=1}^n{a}^i{\mathbf{U}}_i\hfill \\ \hfill \lambda \left(a\right)& ={\lambda }_0+\sum _{i=1}^n{a}^i{\lambda }_i\hfill \end{array}\right.$$

(2)

Let us explain the algebraic properties of the ANM. If we consider a given function $\widehat{\mathbf{Y}}:\mathbf{X}\in {\mathbb{R}}^m⟶\mathbf{Y}=\widehat{\mathbf{Y}}\left(\mathbf{X}\right)\in {\mathbb{R}}^n$, as detailed in [28], the i-th derivative of $\mathbf{Y}\left(a\right)$ is a function of all the derivative of $\mathbf{X}\left(a\right)$ up to order i. Furthermore, we have the general expression:

$${\mathbf{Y}}^{\left(i\right)}=\mathbf{J}\left({\mathbf{X}}_0\right)·{\mathbf{X}}^{\left(i\right)}+{\mathbf{Y}}_i^{nl}\left({\mathbf{X}}^{\left(0\right)},{\mathbf{X}}^{\left(1\right)},\dots ,{\mathbf{X}}^{\left(i-1\right)}\right)$$

(3)

Let us notice that there is a linear relationship between ${\mathbf{Y}}^{\left(i\right)}$ and ${\mathbf{X}}^{\left(i\right)}$ with the same Jacobian matrix $\mathbf{J}\left({\mathbf{X}}_0\right)$ for each order, and, that ${\mathbf{Y}}_i^{nl}$ is a multilinear function of previous orders. As a consequence, the expression (3) for $i=1$ allows to obtain the matrix $\mathbf{J}\left({\mathbf{X}}_0\right)$ since ${\mathbf{Y}}_1^{nl}=0$, and, at order $i>1$, by assuming ${\mathbf{X}}^{\left(i\right)}=0$, the expression ${\mathbf{Y}}_i^{nl}\left({\mathbf{X}}^{\left(0\right)},{\mathbf{X}}^{\left(1\right)},\dots ,{\mathbf{X}}^{\left(i-1\right)}\right)$ is given by ${\mathbf{Y}}^{\left(i\right)}$.

3. ANM for Elastic Nonlinear Problems

3.1. Theory

Let us explain how the Asymptotic Numerical Method (ANM) can be applied in the case of quasi static geometrically nonlinear elastic problems [8,9,10].

The considered elasticity model, known as the Saint Venant-Kirchhoff law, is a linear relation between the second Piola-Kirchhoff stress tensor $\mathbf{S}$ and the Green-Lagrange strain tensor $\mathbf{\gamma }\left(\mathbf{u}\right)$ defined by the following:

$$\gamma \left(\mathbf{u}\right)=\underset{{\gamma }_l\left(\mathbf{u}\right)}{\underbrace{{\displaystyle \frac{1}{2}}\left(\nabla \mathbf{u}+{\nabla }^T\mathbf{u}\right)}}+\underset{{\gamma }_{nl}\left(\mathbf{u},\mathbf{u}\right)}{\underbrace{{\displaystyle \frac{1}{2}}\left({\nabla }^T\mathbf{u}·\nabla \mathbf{u}\right)}}$$

(4)

and its variation is given by:

$$\delta \mathbf{\gamma }\left(\mathbf{u}\right)={\mathbf{\gamma }}_l\left(\mathbf{\delta }\mathbf{u}\right)+2{\mathbf{\gamma }}_{nl}\left(\mathbf{u},\mathbf{\delta }\mathbf{u}\right)$$

(5)

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}{c}\mathbf{S}=\mathbf{D}:\mathbf{\gamma }\left(\mathbf{u}\right)\hfill \\ \\ {\int }_{\mathrm{\Omega }}\mathbf{S}:\delta \mathbf{\gamma }\left(\mathbf{u}\right)d\mathrm{\Omega }-\mathbf{\lambda }{\int }_{{\mathrm{\Gamma }}_N}\mathbf{t}·\mathbf{\delta }\mathbf{u}d\mathrm{\Gamma }=0\hfill \end{array}\right.$$

(6)

where $\mathbf{D}$ is the elastic stiffness tensor, $\lambda$ is a scalar load parameter and $\mathbf{t}$ represents the given force on the part ${\mathrm{\Gamma }}_N$ of the boundary $\partial \mathrm{\Omega }$.

Each step of ANM relies on a Taylor series according to a well chosen path parameter a:

$$\left\{\begin{array}{c}\mathbf{u}\left(a\right)={\mathbf{u}}_0+a{\mathbf{u}}_1+a^2{\mathbf{u}}_2+\dots \hfill \\ \mathbf{S}\left(a\right)={\mathbf{S}}_0+a{\mathbf{S}}_1+a^2{\mathbf{S}}_2+\dots \hfill \\ \lambda \left(a\right)={\lambda }_0+a{\lambda }_1+a^2{\lambda }_2+\dots \hfill \end{array}\right.$$

(7)

Several efficient choices of the path parameter are available [32]. 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 by [27]:

$$a={\displaystyle \frac{1}{{\overline{\mathbf{u}}}^2}}〈\mathbf{u}-{\mathbf{u}}_0,{\mathbf{u}}_1〉+\alpha \left(\lambda -{\lambda }_0\right){\lambda }_1$$

(8)

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

When the Taylor series (7) are substituted into Eqs.(6-8). Identifying mathematical expressions related to a given order p in a results in a linear system. Here, the boundary value problem in Eq.(6) has a simple structure, 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: ${\displaystyle {\left(ab\right)}_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\left(\delta \mathbf{u}\right)$$

(9)

where the bilinear form $k(.,.)$ is the classical tangent stiffness operator given by:

$$\begin{array}{c}\hfill k(\mathbf{u},\delta \mathbf{u})={\int }_{\mathrm{\Omega }}\left(\delta \gamma \left({\mathbf{u}}_0\right):\mathbf{D}:\left({\gamma }_l\left(\mathbf{u}\right)+2{\gamma }_{nl}\left({\mathbf{u}}_0,\mathbf{u}\right)\right)+{\mathbf{S}}_0:2{\gamma }_{nl}\left(\mathbf{u},\delta \mathbf{u}\right)\right)d\mathrm{\Omega }\end{array}$$

(10)

And the linear form related to the virtual work of external forces is given by:

$$P_e\left(\delta \mathbf{u}\right)={\int }_{{\mathrm{\Gamma }}_D}\mathbf{t}·\delta \mathbf{u}d\mathrm{\Gamma }$$

(11)

The definition in Eq. (8) of the path parameter gives an additional equation:

$${\displaystyle \frac{1}{{\overline{\mathbf{u}}}^2}}〈{\mathbf{u}}_1,{\mathbf{u}}_1〉+\alpha {\lambda }_1{\lambda }_1=1$$

(12)

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

$$\begin{array}{c}\hfill k({\mathbf{u}}_p,\delta \mathbf{u})={\lambda }_p{P}_e\left(\delta \mathbf{u}\right)+〈{\mathbf{F}}_p^{nl},\delta \mathbf{u}〉\end{array}$$

(13)

The last term is a linear depends on the variables at orders $k<=\left(p-1\right)$. It is given by:

$$\begin{array}{ccc}\hfill 〈{\mathbf{F}}_p^{nl},\delta \mathbf{u}〉& =& {\displaystyle -2\sum _{r=1}^{p-1}{\int }_{\mathrm{\Omega }}{\mathbf{S}}_r:{\gamma }_{nl}\left({\mathbf{u}}_{p-r},\delta \mathbf{u}\right)d\mathrm{\Omega }}\hfill \\ & -& {\displaystyle \sum _{r=1}^{p-1}{\int }_{\mathrm{\Omega }}{\gamma }_{nl}\left({\mathbf{u}}_r,{\mathbf{u}}_{p-r}\right):\mathbf{D}:\delta \gamma \left({\mathbf{u}}_0\right)d\mathrm{\Omega }}\hfill \end{array}$$

(14)

while the stress on order p is:

$${\mathbf{S}}_p=\mathbf{D}:\left[{\gamma }_l\left({\mathbf{u}}_p\right)+2{\gamma }_{nl}\left({\mathbf{u}}_0,{\mathbf{u}}_p\right)+\sum _{r=1}^{p-1}{\gamma }_{nl}\left({\mathbf{u}}_r,{\mathbf{u}}_{p-r}\right)\right]$$

(15)

The last is the arc-length condition in Eq. (8) yields the load parameter ${\lambda }_p$ as a function of the displacement ${\mathbf{u}}_p$:

$${\displaystyle \frac{1}{{\overline{\mathbf{u}}}^2}}〈{\mathbf{u}}_p,{\mathbf{u}}_1〉+\alpha {\lambda }_p{\lambda }_1=0$$

(16)

Usually, the Taylor expansion order is chosen between 15 and 50 to take advantage of the exponential convergence of power series. Of course, the path parameter a is required to be inferior to $a_{rc}$ the convergence radius of the series 7. In fact, We will ask the truncated Taylor series to verify a given accuracy: the normalized norm of the difference between the two last order in the Taylor series is asked to be inferior to a given parameter $\delta$ which results for the following inequality to be verified: $\left|\right|a^N{\mathbf{u}}_N\left|\right|\le \delta \left|\right|a{\mathbf{u}}_1\left|\right|$. This will leads to the expression of the path parameter $a_{max}$ called validity parameter:

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

(17)

In this way, the Taylor expansion (7) gives a part of the solution path in the interval $\left[0,a_{max}\right]$. The continuation procedure may be very simple by chaining several series, 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 Newton method or the high-order method [33] or by applying a convergence acceleration technique [34]. 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 [35,36].

3.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 the section 7 of [27].

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. The ANM algorithm is detailed below (Algorithm 1):

Algorithm 1 Algorithm ANM

1: Initialize ${\mathbf{u}}_0$ and ${\lambda }_0$ 2: for $step=1$ until $Nstep$ do 3:    Assemble $\left[{\mathbf{K}}_t\left({\mathbf{u}}_0,{\mathbf{S}}_0\right)\right]$ and $\left\{\mathbf{F}\right\}$ 4:    Solve $\left\{\widehat{\mathbf{u}}\right\}={\left[{\mathbf{K}}_t\left({\mathbf{u}}_0,{\mathbf{S}}_0\right)\right]}^{-1}\left\{\mathbf{F}\right\}$ 5:    Compute ${\lambda }_1={\displaystyle \frac{1}{\sqrt{1+{\left\{\widehat{\mathbf{u}}\right\}}^T\left\{\widehat{\mathbf{u}}\right\}}}}$ 6:    Compute $\left\{{\mathbf{u}}_1\right\}={\lambda }_1\left\{\widehat{\mathbf{u}}\right\}$ 7:    Compute (for each Gauss point) $\left\{{\mathbf{S}}_1\right\}=\left[\mathbf{D}\right]\left(\left\{{\gamma }_l\left({\mathbf{u}}_1\right)\right\}+\left\{{\gamma }_{nl}\left({\mathbf{u}}_0,{\mathbf{u}}_1\right)\right\}\right)$ 8:    Compute (for each Gauss point) $\left\{{\mathbf{S}}_2^{nl}\right\}=\left[\mathbf{D}\right]\left\{{\gamma }_{nl}\left({\mathbf{u}}_1,{\mathbf{u}}_1\right)\right\}$ 9:    Compute (for each Gauss point) $\left\{{\mathbf{S}}_2^{*}\right\}={\left[\mathbf{A}\left({\mathbf{u}}_1\right)\right]}^T\left\{{\mathbf{S}}_1\right\}$ 10:    Compute $\left\{{\mathbf{F}}_2^{nl}\right\}$ 11:    for $p=2$ until N do 12:      Compute $\left\{{\mathbf{u}}_p^{nl}\right\}={\left[{\mathbf{K}}_t\left({\mathbf{u}}_0,{\mathbf{S}}_0\right)\right]}^{-1}\left\{{\mathbf{F}}_p^{nl}\right\}$ 13:      Compute ${\lambda }_p=-{\lambda }_1{\left\{{\mathbf{u}}_p^{nl}\right\}}^T\left\{{\mathbf{u}}_1\right\}$ 14:      Compute $\left\{{\mathbf{u}}_p\right\}={\lambda }_p\left\{\widehat{\mathbf{u}}\right\}+\left\{{\mathbf{u}}_p^{nl}\right\}$ 15:      Compute (for each Gauss point) $\left\{{\mathbf{S}}_p\right\}=\left[\mathbf{D}\right]\left(\left\{{\gamma }_l\left({\mathbf{u}}_p\right)\right\}+\left\{{\gamma }_{nl}\left({\mathbf{u}}_0,{\mathbf{u}}_p\right)\right\}\right)+\left\{{\mathbf{S}}_p^{nl}\right\}$ 16:      Compute (for each Gauss point) $\left\{{\mathbf{S}}_{p+1}^{nl}\right\}=\left[\mathbf{D}\right]\left({\sum }_{r=1}^p\left\{{\gamma }_{nl}\left({\mathbf{u}}_r,{\mathbf{u}}_{p+1-r}\right)\right\}\right)$ 17:      Compute (for each Gauss point) $\left\{{\mathbf{S}}_{p+1}^{*}\right\}={\sum }_{r=1}^p\left({\left[\mathbf{A}\left({\mathbf{u}}_{p+1-r}\right)\right]}^T\left\{{\mathbf{S}}_r\right\}\right)$ 18:      Compute $\left\{{\mathbf{F}}_{p+1}^{nl}\right\}$ 19:    end for 20:    Compute $a_{max}$ 21:    Compute $\lambda \left(a\right)$ and $\mathbf{u}\left(a\right)$ with $a\in \left[0,a_{max}\right]$ 22:    Compute the normalized residual error 23:    Actualize ${\mathbf{u}}_0$ and ${\lambda }_0$ 24: end for

The FEM stiffness matrix $\left[{\mathbf{K}}_t\left({\mathbf{u}}_0,{\mathbf{S}}_0\right)\right]$ is the implementation of the bilinear form (10), the FEM right-hand side vector $\left\{\mathbf{F}\right\}$ is the implementation of the linear form $P_e\left(\delta \mathbf{u}\right)$ (11), and the FEM right-hand side vector $\left\{{\mathbf{F}}_p^{nl}\right\}$ is the implementation of the linear form $〈{\mathbf{F}}_p^{nl},\delta \mathbf{u}〉$ (14). The main details of the implementation in FreeFEM++ will be given in the Section 5. The algorithm A3 of the Newton-Riks correction is given in the Appendix (B) as well as the algorithm A2 of Newton-Raphson predictions which will be used later for comparison purpose.

4. ANM for Thermo-Elastic Nonlinear Problems

4.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. As in the previous section, we assume that the domain $\mathrm{\Omega }$ is a Saint Venant-Kirchhoff material and we will consider nonlinear geometric elasticity, as in previous section. This domain will be subjected to isotropic thermal expansion (or shrinkage), driven by the scalar parameter $\lambda$.

We will assume that the thermomechanical effects result in an additional contribution to the Green-Lagrange strain tensor (4)), and its expression is $\lambda \mathbf{I}$ (where $\mathbf{I}$ is the identity matrix). Notice that the expression of $\delta \mathbf{\gamma }\left(\mathbf{u}\right)$ (5) is not changed.

Incorporating thermal effects, the virtual work equation (6) results in the following:

$$\left\{\begin{array}{c}\mathbf{S}=\mathbf{D}:\left({\mathbf{\gamma }}_l\left(\mathbf{u}\right)+{\mathbf{\gamma }}_{nl}\left(\mathbf{u},\mathbf{u}\right)-\mathbf{\lambda }\mathbf{I}\right)\hfill \\ \\ {\int }_{\mathrm{\Omega }}\mathbf{S}:\delta \mathbf{\gamma }\left(\mathbf{u}\right)d\mathrm{\Omega }=0\hfill \end{array}\right.$$

(18)

$\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 ${\mathbf{S}}_p$ at the order $p\ge 1$:

$${\mathbf{S}}_p=\mathbf{D}:\left({\mathbf{\gamma }}_l\left({\mathbf{u}}_p\right)+2{\mathbf{\gamma }}_{nl}\left({\mathbf{u}}_0,{\mathbf{u}}_p\right)+\sum _{k=1}^{p-1}{\mathbf{\gamma }}_{nl}\left({\mathbf{u}}_k,{\mathbf{u}}_{p-k}\right)-{\mathbf{\lambda }}_p\mathbf{I}\right)$$

(19)

Inserting (19) into (18), we obtain at the order $p=1$:

$$k({\mathbf{u}}_1,\delta \mathbf{u})={\lambda }_1{P}_{\mathrm{th}}\left(\delta \mathbf{u}\right)$$

(20)

The additional equation for the path parameter is the same than Eq. 12.

The bilinear form $k(.,.)$ is the classical tangent stiffness operator already given in (10), and the linear form $P_{\mathrm{th}}\left(.\right)$ that accounts for the thermomechanical shrinkage is given by:

$$P_{\mathrm{th}}\left(\delta \mathbf{u}\right)={\int }_{\mathrm{\Omega }}I:\mathbf{D}:\delta \mathbf{\gamma }\left({\mathbf{u}}_{\mathbf{0}}\right)d\mathrm{\Omega }$$

(21)

Notice that the integrand on the right member of (21) is the trace of the second-rank tensor $\left(\mathbf{D}:\left({\mathbf{\gamma }}_l\left(\mathbf{\delta }\mathbf{u}\right)+{\mathbf{\gamma }}_{nl}\left({\mathbf{u}}_0,\mathbf{\delta }\mathbf{u}\right)\right)\right)$.

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

$$\begin{array}{c}\hfill k({\mathbf{u}}_p,\delta \mathbf{u})={\lambda }_p{P}_{\mathrm{th}}\left(\delta \mathbf{u}\right)+〈{\mathbf{F}}_p^{nl},\delta \mathbf{u}〉\end{array}$$

(22)

The additional equation for the path parameter is the same than Eq. 16.

4.2. Numerical Algorithm

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

5. Some Aspects of 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 on surface to be solved. The FreeFEM++ language is typed, polymorphic and reentrant with macro generation. The kernel of FreeFEM++ is written in C++.

Let us now resume the main characteristics of FreeFEM++:

• FreeFEM++ has a high level user friendly typed input language with an algebra of analytic and finite element functions.
• The finite element problems are described by their variational formulations, and are easily implemented. It is possible to have access to the internal vectors and matrices if needed.
• Many kinds of problems can be solved: multi-variables, multi-equations, one, two or three dimensional static or time dependent, linear or nonlinear coupled systems.
• An automatic mesh generator is available, based on the Delaunay-Voronoï algorithm.
• A large variety of triangular finite elements: linear and quadratic Lagrangian finite elements and more, discontinuous P1, and Raviart-Thomas elements.
• A large variety of linear direct and iterative solvers (LU, Cholesky, Crout, CG, GMRES, UMFPACK, MUMPS, SuperLU, ...), 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. As we will see in more detail in the following, the mesh generator is very easy to use for simple geometries, the macro generation enables us an easy implementation of the differential operators useful in nonlinear solid mechanics, and the variational formulation oriented commands in FreeFEM++ allows us to implement very efficiently the weak form of nonlinear solid mechanical differential equations, discretized using the Finite Element Method.

Let us present the main steps of the Finite Element Method discretization by getting closer to the FreeFEM++ implementation. In order to set these ideas down, let us begin with the Laplacian problem, a basic problem for applied mathematics.

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

$$-\Delta u=f\mathrm{in}\mathrm{\Omega }$$

(23a)

$$u=0\mathrm{on}{\mathrm{\Gamma }}_D$$

(23b)

$${\displaystyle \frac{\partial u}{\partial n}}=g\mathrm{on}{\mathrm{\Gamma }}_N$$

(23c)

with ${\mathrm{\Gamma }}_D\cup {\mathrm{\Gamma }}_N=\partial \mathrm{\Omega }$. The boundary condition on ${\mathrm{\Gamma }}_D$ (resp. on ${\mathrm{\Gamma }}_N$) is called Dirichlet (resp. Neumann) boundary conditions.

Let us define the functional space $V_0$ by:

$$V_0\left(\mathrm{\Omega }\right)=\left\{v\in H^1\left(\mathrm{\Omega }\right),v_{|{\mathrm{\Gamma }}_D}=0\right\}$$

(24)

Using the Green’s formula that generalizes the integration by part, an equivalent weak form (or variational form) of the previous problem (23) is obtained: Find $u\in V_0\left(\mathrm{\Omega }\right)$ such that:

$${\int }_{\mathrm{\Omega }}\nabla u\nabla vd\mathrm{\Omega }={\int }_{\mathrm{\Omega }}fvd\mathrm{\Omega }+{\int }_{{\mathrm{\Gamma }}_D}v{\displaystyle \frac{\partial u}{\partial n}}ds,\forall v\in V_0\left(\mathrm{\Omega }\right)$$

(25)

The Finite Element Method starts with the partitioning of the initial domain $\mathrm{\Omega }$ into subdomains ${\mathrm{\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:

$$w\left(x,y,z\right)=\sum _{i=1}^N{w}_i{\varphi }_i\left(x,y,z\right)$$

(26)

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

The finite element space $V_h$ is defined by:

$$V_h=\left\{w=\sum _{i=1}^N{w}_i{\varphi }_i,w_i\in \mathbb{R}\right\}$$

(27)

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 (23). 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 done in FreeFEM++ with the two command lines:

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

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

• macro Grad(u) [dx(u), dy(u), dz(u)] //

Let us define u the solution function, v the test function, belonging to the finite element space $V_h$, and implement the equivalent variational formulation of the Laplacian problem:

• Vh u,v;
• problem Poisson(u,v) =
• int3d(Th)(Grad(u)’*Grad(v))-int3d(Th)(f*v)
• -int2d(Th,LN)(g*v)
• +on(LD,u=0);

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

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

• Poisson;

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

• plot(u,wait=1);

It is easy to generalize the FreeFEM++ implementation of the Laplacian problem to the resolution of nonlinear mechanical problem using ANM. FreeFEM++ enables to generate directly the sparse Finite Element matrices, and the right-hand sides, and this is usefull regarding the implementation of 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 $\left[{\mathbf{K}}_t\left({\mathbf{u}}_0,{\mathbf{S}}_0\right)\right]$ is the implementation of the bilinear form (10), the FEM right-hand side vector $\left\{\mathbf{F}\right\}$ is the implementation of the linear form $P_e\left(\delta \mathbf{u}\right)$ (11), and the FEM right-hand side vector $\left\{{\mathbf{F}}_p^{nl}\right\}$ is the implementation of the linear form $〈{\mathbf{F}}_p^{nl},\delta \mathbf{u}〉$ (14).

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 strength 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] //
• 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)) ) its nonlinear part (4). The tensor dGammaNL(u,v,w,uu,vv,ww) is the variation of the non linear part of the Green Lagrange tensor (5).

Then, the elasticity matrix is easily implemented:

• 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 Poisson ratio, lambda, 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 which will lead to the computation of both stiffness matrix, and right-hand side vector.

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

• fespace Vh(Th3D,[P2,P2,P2]);

Let us declare the displacement fields for each ANM order:

• 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).

• 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:

• 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 of all, the finite element displacement field $\left\{\widehat{\mathbf{u}}\right\}$ is computed using a linear system resolution (we used the MUMPS direct solver, which has proven is efficiency).

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

It is now necessary to assemble the nonlinear right-hand side vector $\left\{{\mathbf{F}}_{\mathbf{2}}^{\mathbf{nl}}\right\}$. Like it is done in the ANM algorithm (1), we define intermediate tensors (written in a vectorial form) $\left\{{\mathbf{S}}_1\right\}$ and $\left\{{\mathbf{S}}_2^{nl}\right\}$ 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[Snlxx,Snlyy,Snlzz,deuxSnlxy,deuxSnlxz,deuxSnlyz](Norder+1);

The first step is to compute $\left\{{\mathbf{S}}_1\right\}$:

• [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 $\left\{{\mathbf{S}}_2^{nl}\right\}$:

• [Snlxx,Snlyy,Snlzz,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 $\left\{{\mathbf{F}}_{\mathbf{2}}^{\mathbf{nl}}\right\}$:

• varf PbFnl2 ([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))’*
• ([Snlxx,Snlyy,Snlzz,deuxSnlxy,deuxSnlxz,deuxSnlyz]) )
• + on(lencasmid,u1=0.,v1=0.,w1=0.);
• Fnlu[]=PbFnl2(0,Vh);

Let us now enter the loop with respect to the order $p\ge 2$. All the FreeFEM++ implementation details are given in .

Let us now compute the norm $L_2$ of ${\mathbf{u}}_1$ and ${\mathbf{u}}_N$ to determine $a_{max}$ (17):

• 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 in evaluating the Taylor series of the displacement nodal vector, computing the residual error, and actualizing both ${\mathbf{u}}_0$ and ${\lambda }_0$.

Let us mention that the FreeFEM++ scripts for the FEM simulation of non linear solid mechanics using ANM will be available soon in the FreeFEM++ web site (https://freefem.org/).

6. Numerical Results

In order to validate the FreeFEM++ program, the first section will focus on a simulation example based on a simple geometry. This is a clamped beam subjected to a vertical conservative surface traction. This problem has already been studied by Zahrouni [37] in the case of shell finite elements and a concentrated force. For validation purpose, we will add some simulation results obtained with the commercial software Abaqus. In the following section, we will briefly show that the FreeFEM++ simulation tool is able to simulate more complex structures like film/substrate systems under axial compression.

6.1. Clamped Beam Subjected to a Conservative Vertical Surface Traction

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 a conservative vertical traction is applied to the left face (see Figure 1).

The beam has been meshed using tetraedral P2 Lagrangian finite elements (see Figure 2).

The beam is supposed to be an homogeneous isotropic elastic material with Young’s modulus of ${10}^5\mathrm{M}\mathrm{Pa}$, and 0 Poisson’s ratio. A FreeFEM++ script has been developed in order to simulate the non linear elastic behavior of the clamped beam. For building the finite element mesh of the beam, we choose 40 elements for its length, and 4 elements for its 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\mathrm{M}\mathrm{Pa}$. Let us recall that the loading parameter $\lambda$ appears in the weak form related to the applied force (right term of the second equation of (6)). Figure 3 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 a very good agreement. Large deflection up to $u_z=8.5$ have been obtained with only 7 ANM steps. At the end of each step, the normalized residual error is close to ${10}^{-6}$. Figure 3 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.

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

6.2. Film-Substrate Systems Under Uniaxial Compression

We are now going to consider a more complex nonlinear mechanical problem: the film-substrate system. Surface morphological instabilities of a soft material with a stiff thin surface layer have been extensively studied over the past years.

When depositing a stiff thin film on a polymeric soft substrate, the developed residual compressive stresses in the film during the cooling process due to the large mismatch of thermal expansion coefficient between the film and the substrate, are relieved by buckling with wrinkling patterns, which was pioneered by Bowden et al. [38] in 1998. Let us notice that wrinkles of a thin stiff film deposited on a soft substrate can be widely used in industry for lots of applications, like the micro/nano fabrication of flexible electronic devices with controlled morphological patterns [39] [40], to the mechanical property measurement of material characteristic [41].

The film-substrate system being studied consists of a stiff film deposited on a soft rectangle parallelepiped substrate (Figure 4).

Film substrate systems has already been studied with the help of ANM using a coupling between shell finite element (Büchter-Ramm [42]) and hexahedron finite elements by Fan Xu et al. ([25,26]). In this section we will present a numerical model with ANM using only tetrahedron finite element (because of 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 modelled film-substrate system (only one fourth of the structure will be modelled). Moreover, on the interfaces ${\Sigma }_d$, and, ${\Sigma }_l$ the vertical displacement is supposed zero. The length, width, and height of the substrate are respectively $1.5$ $\mathrm{m}$$\mathrm{m}$, $0.75$ $\mathrm{m}$$\mathrm{m}$, and, $0.1$ $\mathrm{m}$$\mathrm{m}$. The thickness of the film is $0.001$ $\mathrm{m}$$\mathrm{m}$.

The substrate and the film are supposed to be homogeneous isotropic and elastic materials. Young’s modulus of the substrate (resp. film), is $E_s=1.8$ (resp. $E_f=1.3\mathrm{e}5$), and Poisson ratio of the substrate (resp. film) is ${\nu }_s=0.48$ (resp. ${\nu }_f=0.3$). Let us notice that the Young’s modulus of the film is much larger than that of the substrate ($E_f/E_s=0.72\times {10}^5$), which will allow to assume linear elasticity for the substrate domain while nonlinear geometrical elasticity will be assumed for the film domain [43]. The number of elements for both 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 in 360000 tetraedra, and 1575639 degrees of freedom. We choose an ANM strategy with quite small steps ($\delta ={10}^{-6}$), and the ANM order is 15. One ANM step requires almost 40 minutes of CPU time when the FreeFEM++ program is parallized on 6 MPI processors while using the multifrontal MUMPS solver.

According to recent works [34] regarding acceleration convergence improvement of the ANM algorithm, MMPE extrapolation, adaptatives steps, and Newton Riks correction have been used. Let us give some details about these improvements. 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 to reduce Newton Riks’s corrections.

The residual before and after the application of the acceleration MMPE are presented in (Table 2). Seven starting points $\mathbf{u}\left(r{a}_{max}\right)$ with $r\in \left[0.7,1.3\right]$ 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 to improve the normalized residual error.

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

Many bifurcations can be seen on the diagram of Figure (Figure 5). At the first bifurcation for $\lambda =0.048$, uniformly distributed one-dimensional wrinkles appears. Next the curve goes back for $\lambda =0.09$ and five limit points are obtained that corresponds to the increase of a single wrinkle near the boundary and to a change of the number of wrinkles [34]. Small steps accumulate close to the limit points, which indicates that these singularities are in fact perturbed bifurcation points [44].

Figure 5. Bifurcation diagram: vertical displacement of the center of the film-substrate system (top surface) with $E_f/E_s=0.72\times {10}^5$.

Figure 5. Bifurcation diagram: vertical displacement of the center of the film-substrate system (top surface) with $E_f/E_s=0.72\times {10}^5$.

Figure 6. Bifurcation diagram: vertical displacement of the center of the film-substrate system (top surface) with $E_f/E_s=0.72\times {10}^5$: comparison ANM full algorithm and Newton-Raphson continuation algorithm wit Newton-Riks like corrections.

Figure 6. Bifurcation diagram: vertical displacement of the center of the film-substrate system (top surface) with $E_f/E_s=0.72\times {10}^5$: comparison ANM full algorithm and Newton-Raphson continuation algorithm wit Newton-Riks like corrections.

Figure 7. Deformed of the film top surface just before the first bifurcation (ANM step 20).

Figure 7. Deformed of the film top surface just before the first bifurcation (ANM step 20).

Figure 8. Deformed of the film top surface just before the second bifurcation (ANM step 34).

Figure 8. Deformed of the film top surface just before the second bifurcation (ANM step 34).

6.3. Spherical Film-Substrate Under Thermal Loading

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

Figure 9 is a plot of the mesh of one eighth of the whole spherical film/substrate system (three symmetry planes are used to obtain the whole spherical structure). We can see the tetrahedron geometry of the mesh and the presence of smaller element within the film part of the mesh. The spherical film-substrate geometry has been meshed with two elements in the film thickness, an average elements length $h=0.024mm$ and consists in 164528 tetrahedra (689358 degrees of freedom). We choose four different values of the parameter $C_s$ (0.2, 2.6, 4.9, 21.2), and, the mesh is the same for the four parameters.

Let us first make the numerical simulation for $C_s=0.2$. Figure 10 is the plot of the maximum normalized deflection versus the thermal loading parameter $\lambda$. We can observe two limit or turning points, the bifurcation curve being subcritical: the first one is close to step 10, and the second close to step 40. Then, close to step 100, another bifurcation happens and the curve is going backward. The ANM continuation does not allow to jump on another probable existing branch: future research will be done to solve this lack. Figure 11 shows the deformation (normalized deflection) for an half sphere at the steps 10, 34, 50, and 80. As in [3], we can observe isolated dimples, with a greater localization when deflection amplitude increased.

Second, let us now consider the numerical simulation for $C_s=2.6$. Figure 12 is the plot of the maximum normalized deflection versus the thermal loading parameter $\lambda$. We can observe first a limit or turning point, the bifurcation being weakly subcritical. We can observe small back-and-forth motions and the ANM continuation still goes on. Figure 13 shows the deformation (normalized deflection) for an half sphere, exhibiting buckyball deformation patterns.

Third, let us consider the case of $C_s=4.9$. Figure 14 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 is showing zigzag and the progression slows down. For $C_s=4.9$ Figure 15 we can observe a periodicity in the pattern, and a transition between buckyball and a symmetry breaking of buckyball pattern, around step 120 (the deflection amplitude of some dimples decreases).

Finally, let us consider $C_s=21.2$. Figure 16 is the plot of the maximum normalized deflection versus the thermal loading parameter $\lambda$. The first bifurcation is also supercritical as previously, and appears close to step 10. Then the bifurcation curve behavior is almost linear and reaches greater thermal loading values than in the previous case. Figure 17, we can observe a smaller periodicity in the patter and also a transition between buckyball and a symmetry breaking of buckyball pattern.

7. Conclusions

In this paper, we have presented a new implementation of the Asymptotic Numerical Method within the user-friendly finite element environment FreeFem++. ANM continuation method, based on a perturbation technique, has proven to be a very efficient technique for the numerical simulation of nonlinear instability problems.The ANM has been previously implemented in FORTRAN, Matlab, and C++ to solve various problems of fluid or solid mechanics.

On the one hand, 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. On the other hand, this new numerical tool written in FreeFem++ will be very helpful for new developments, like those concerning convergence acceleration [34] or multiscale finite elements [45,46].

A first computation using the new developed FreeFem++ code, consisting in a cantilever (Saint Venant-Kirchhoff material) submitted to a vertical surface force, where geometrical nonlinearities are taken into account, has shown a very good comparison with the commercial software Abaqus. The two other applications concern film/substrate systems. 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 applications have demonstrated the efficiency of the implementation of ANM in FreeFem++ for the following of equilibrium paths when solving solid mechanical problems.

In the present state, the achieved implementation concerns geometrically nonlinear elasticity, but the extension to other models of solid and fluid mechanics seems rather easy. Among the contributions of FreeFem++, let’s mention the treatment of multiphysics modelings or of large scale problems thanks to its interface with the PETSc library. In this paper, the MPI/MUMPS parallel version has been applied.

Finally, a supplementary article shows the details of the FreeFEM++ implementations.