# On the use of an high order perturbation method for numerical time integration in structural dynamics.

Bertille Claude, Gregory Girault, Bruno Leblé, Jean-Marc Cadou

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

## To cite this version:

Bertille Claude, Gregory Girault, Bruno Leblé, Jean-Marc Cadou. On the use of an high order perturbation method for numerical time integration in structural dynamics. Comptes Rendus. Mécanique, 2023, 351, pp.227-245. 10.5802/crmeca.195. hal-04409835

HAL Id: hal-04409835

https://hal.science/hal-04409835v1

Submitted on 22 Jan 2024

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

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

<figure><img src="image_2.png" /></figure>
# Comptes Rendus
## Mécanique

Bertille Claude, Grégory Girault, Bruno Leblé and Jean-Marc Cadou

**On the use of an high order perturbation method for numerical time integration in structural dynamics.**

Volume 351 (2023), p. 227-245

Published online: 19 June 2023

https://doi.org/10.5802/crmeca.195

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

This article is licensed under the
CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.
http://creativecommons.org/licenses/by/4.0/

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

Les Comptes Rendus. Mécanique sont membres du
Centre Mersenne pour l'édition scientifique ouverte
www.centre-mersenne.org
e-ISSN : 1873-7234
numerical damping features, allowing for larger time steps, and the use of adaptive time stepping schemes.

More recently, a serie of studies [6–9] were focused on implicit time integration for wave propagations in structural dynamics. Let cite the proposal of the ρ∞-Bathe method. This latter is based on a composite strategy for time integration containing as special cases the Bathe and Newmark methods.

This topic also emerges in non-smooth dynamics featured by a lack of regularity, such as contact or multibody problems. For this kind of problems, the previous time-stepping schemes are used [10]. Let note that alternative approaches based on event-driven schemes may be implemented [11, 12].

Previous approaches produce discrete time solutions. An alternative improvement is to compute a continuous solution in time. Thus, perturbation methods can be used and time variable becomes the perturbation parameter.

In this context, the Asymptotic Numerical Method is a semi-analytical method where the unknowns are expanded in power series with respect to the time variable. By grouping terms of the same power of the time variable, a set of linear systems are obtained with same matrices.

This strategy was employed to solve ordinary differential equations encountered in various fields of physics. Various problems have been solved, such as neutron kinetics equations [13, 14], trajectory of a conservative non-linear pendulum [15] or the famous Duffing equation [16].

The Asymptotic Numerical Method was also employed to solve linear and non-linear structural problems in different ways. In [17], the ANM was used as an implicit time solver to study dynamic buckling of a cylindrical shell. Discrete solution was computed with a generalized $\alpha$-method as time solver. An homotopy technique was also used to keep valid for several time steps the Jacobian matrix in order to save computational time. The buckling instant was perfectly detected. Unfortunatly, this strategy failed to compute the entire response once the buckling appeared, the Jacobian matrix becoming singular.

Probably first attempts to use ANM as an explicit time solver [18, 19], it was proposed to use this semi-analytical method to compute time-dependent displacement solution. Velocity and acceleration were easily computed by differentiation. Nevertheless, these studies were restricted to linear elastic structure excited by an harmonic external sollication.

Thus, it is proposed to apply the methodology presented in [16, 19] to compute transient responses for structural dynamics. Problems of the elastic wave propagation in an elastic bar and the dynamic buckling of a cylindrical shell are adressed combining an ANM based time solver with finite element spatial discretizations.

Solutions of transient dynamic problems are not analytical in times, then it may seem surprising to use time asymptotic expansions to solve this kind of problems. Nevertheless, as underlying by Deeb et al. in [20], there are two possibilities when using time perturbation method for solving non linear transient problems. The first one is that the series is convergent and then a continuation technique based on the evaluation of the validity range of the asymptotic expansions is sufficient to get the whole non linear solution. Padé approximants [21], for example, can then be used to increase the validity range and consequently decrease the step numbers of the continuation technique. The second possibility is that this validity range becomes very small or worstly null. In this case, a resummation procedure must to be applied to the series expansion and then lead to a time-analytical representation of the solution. Several resummation methods exist in the literature : for example Borel-Padé-Laplace (BPL) [22], Inverse factorial series [23] or generalized hypergeometric Meijer G-function [24]. Let us note that the BPL scheme seems to be today the only numerical method that can be used to solve transient linear or non linear problems. Aim of this study is the computation of semi-analytical continuous time solutions even for non linear problem without using any resummation technique.
Equations of motion are derived and the ANM algorithm is recalled in Section 1. Advantages and drawbacks of the method are exposed and influences of numerical parameters are evaluated on numerical precision and computation efficiency in Section 2. Finally, some recommandations are proposed in Section 3 to enhance the methodology.

# 1. Methods

## 1.1. *Formulation of the problem*

Let $t$ be the time variable and $\mathbf{u}$ the displacement unknown. Dot symbol stands for time derivative $\frac{d\mathbf{u}}{dt} = \dot{\mathbf{u}}$. For an elastic structure with structural damping and submitted to some external sollicitation, equations of motion $\mathcal{R}(\mathbf{u}, \dot{\mathbf{u}}, \ddot{\mathbf{u}}) = \mathbf{0}$ read:

$$ \mathcal{M}(\ddot{\mathbf{u}}) + \mathcal{L}(\mathbf{u}) + \mathcal{D}(\dot{\mathbf{u}}) + \mathcal{Q}(\mathbf{u}, \mathbf{u}) - \mathcal{F}(t) = \mathbf{0} \quad (1) $$

where $\mathcal{L}$ and $\mathcal{D}$ are linear operators standing for elastic and damping properties and $\mathcal{M}$ is a constant mass operator. The non-linear operator $\mathcal{Q}$ stands for the geometrically non-linear effects. Operator $\mathcal{F}$ contains the forcing of the structure. Vectors $\mathbf{u}$, $\dot{\mathbf{u}}$ and $\ddot{\mathbf{u}}$ stand, respectively, for displacement, velocity and acceleration fields of the structure.

For spatial discretization, a finite element method is applied to equation (1). Finite dimensional operators $\mathcal{M}$, $\mathcal{L}$, $\mathcal{N}$ and $\mathcal{F}$ are the discretized counterpart of previous operators in (1) and vector $\mathbf{q}$ is the vector of nodal unknowns associated with the variational formulation. Thus, the following set of ordinary differential equations $\mathcal{R}(\mathbf{q}, \dot{\mathbf{q}}, \ddot{\mathbf{q}}) = \mathbf{0}$ is obtained:

$$ \mathcal{M}\ddot{\mathbf{q}} + \mathcal{L}(\mathbf{q}) + \mathcal{D}(\dot{\mathbf{q}}) + \mathcal{Q}(\mathbf{q}, \mathbf{q}) - \mathbf{F}(t) = \mathbf{0} \quad (2) $$

To solve this problem, it is proposed to apply the ANM on the variable $t = t_0 + \hat{t}$, by considering $\hat{t}$ as the perturbation parameter and $t_0$ stands for the initial time.

## 1.2. *A perturbation-continuation based solver*

The use of the ANM to solve non-linear quasi-static problem in structural mechanics is well established. A high-order perturbation method is associated with a continuation technique in order to determine a semi-analytical solution [25, 26]. This procedure was applied in fluid mechanics to determine bifurcation points [27, 28] and to determine global and local elastic instabilities [29–31]. It is worth noticing that solution obtained with the ANM is independent from the discretization technique. It was shown that meshless method [32, 33], finite difference method [34] and finite element method [29, 35] can be used.

Perturbation in time is introduced as $t = t_0 + \hat{t}$. Vectors $\mathbf{q}$ and $\mathbf{F}$ are sought as power series of time $\hat{t}$ of order $N$ and it is supposed that $\mathbf{q}_0$ and $\mathbf{F}_0$ at instant $t_0$ are known:

$$ \left\{ \begin{array}{l} \mathbf{q}(\hat{t}) = \sum_{i=0}^{N} \hat{t}^i \mathbf{q}_i \\ \mathbf{F}(\hat{t}) = \sum_{i=0}^{N} \hat{t}^i \mathbf{F}_i \end{array} \right. \quad (3) $$

Expressions (3) are introduced in equation (2). After identifying same power of $\hat{t}$, a set of linear systems are obtained. A key point in applying ANM is the quadratic recast of eq. (2), specially the operator $\mathcal{Q}$. Full details can be found in [36, 37]. Expansion of external force $\mathbf{F}$ requires a specific treatment with introduction of a regularization parameter (see Section 1.3). This procedure was defined in previous works for calculation of stress in plastic buckling [38, 39].
Terms $\mathbf{q}_0$ and $\mathbf{q}_1$ supposed to be knowns, identification of equal power of $\hat{t}$ allows to define the 2nd order term:

$$2\mathbb{M}\mathbf{q}_2 + \mathbb{L}(\mathbf{q}_0) + \mathbb{D}(\mathbf{q}_1) + \mathbb{Q}(\mathbf{q}_0, \mathbf{q}_1) - \mathbf{F}_0 = \mathbf{0} \quad (4)$$

Defining the operator $\mathbb{L}_{\mathbf{q}_0}(\mathbf{q}_k) = \mathbb{L}(\mathbf{q}_k) + \mathbb{Q}(\mathbf{q}_0, \mathbf{q}_k) + \mathbb{Q}(\mathbf{q}_k, \mathbf{q}_0)$, terms of higher order ($k = 1, \dots, N-2$) can be computed:

$$(k+2)(k+1)\mathbb{M}\mathbf{q}_{k+2} + \mathbb{L}_{\mathbf{q}_0}(\mathbf{q}_k) + (k+1)\mathbb{D}\mathbf{q}_{k+1} + \sum_{r=1}^{k-1} \mathbb{Q}(\mathbf{q}_r, \mathbf{q}_{k-r}) - \mathbf{F}_k = \mathbf{0} \quad (5)$$

Solving problem (2) consist in solving eq. (4)-(5). To progress in time, a continuation technique is applied. The validity of one step is evaluated according to the following criteria adapted from works on buckling analysis with the ANM [26] where $\delta$ is a user-defined numerical parameter:

$$\hat{t}_{\max} = \left( \delta \frac{\|\mathbf{q}_1\|}{\|\mathbf{q}_N\|} \right)^{\frac{1}{N-1}} \quad (6)$$

A new continuation step is initialized from solution obtained at the actual step according to the following procedure:

$$\left\{ \begin{array}{l} t_0^{\text{new}} = t_0 + \hat{t}_{\max} \\ \mathbf{q}_0^{\text{new}} = \sum_{i=0}^{N} \hat{t}_{\max}^i \mathbf{q}_i \\ \mathbf{q}_1^{\text{new}} = \sum_{i=1}^{N} i \hat{t}_{\max}^{i-1} \mathbf{q}_i \\ \mathbf{F}_0^{\text{new}} = \sum_{i=0}^{N} \hat{t}_{\max}^i \mathbf{F}_i \end{array} \right. \quad (7)$$

Because solution $\mathbf{q}$ is expanded into power series, computation of velocity by differentiation is done very easily. Finally, let note that this process is initialized, for the very first step, by knowledge of the initial conditions $\mathbf{q}_{\text{ini}}$ and $\dot{\mathbf{q}}_{\text{ini}}$:

$$\left\{ \begin{array}{l} \mathbf{q}_0 = \mathbf{q}_{\text{ini}} \\ \mathbf{q}_1 = \dot{\mathbf{q}}_{\text{ini}} \end{array} \right. \quad (8)$$

### 1.3. Recast of the external forcing F

In transient dynamics, structures are submitted to external forcing $\mathbf{F}(t)$. Models of forcing are generally combination of Dirac, Heaviside, geometric and hypergeometric functions. To set these kind of functions in a quadratic framework, specific recast is necessary.

In this study, a ramp function (see Figure 6b) is considered as forcing on cylindrical shell and can be defined by:

$$st = f(t) + \eta \frac{f_c^2}{f_c^2 - (f(t))^2} f(t) \quad (9)$$

where $s = f_c/t_c$, $t_c$ and $f_c$ stand, respectively, for the ramp slope, the duration and the maximum intensity of forcing. The symbol $\eta$ stands for a regularization parameter.

The quadratic recast is inspired from previous works on yield stress calculation for elastic/perfectly plastic material [38, 39]. After definition of auxiliary variable $g(t) = f(t)^2$ and parameters $\xi = f_c^2$, $\zeta = s\xi$ and $\gamma = \eta\xi$, expression (9) is written:

$$(t_0 + \hat{t})\zeta - s(t_0 + \hat{t})g - (\zeta + \gamma)f + fg = 0 \quad (10)$$

This quadratic recast allows to define power series expansion of $f(t)$ and $g(t)$ and to identify terms of identical order in $\hat{t}$. Thus, the following sets of linear systems are obtained:

$$\text{Order 0:} \left\{ \begin{array}{l} g_0 - f_0^2 = 0 \\ \zeta t_0 - s t_0 g_0 - (\zeta + \gamma) f_0 + f_0 g_0 = 0 \end{array} \right. \quad (11)$$
$$ \text{Order 1:} \left\{ \begin{array}{l} g_1 - 2f_0f_1 = 0 \\ \zeta - st_0g_1 - sg_0 - (\zeta + \gamma)f_1 + f_0g_1 + f_1g_0 = 0 \end{array} \right. \qquad (12) $$

$$ \text{Order } k \in [2, N] : \left\{ \begin{array}{l} g_k - 2f_0f_k - \sum_{r=1}^{k-1} f_r f_{(k-r)} = 0 \\ -st_0g_k - sg_{(k-1)} - (\zeta + \gamma)f_k + f_0g_k + f_kg_0 + \sum_{r=1}^{k-1} f_r g_{(k-r)} = 0 \end{array} \right. \qquad (13) $$

These equations are invertible and require knowledge of $f_0$, from which $g_0$ is easily deduced.

As for all power series expansion, a validity range must be defined according to the previous definition (6). Two validity durations $\hat{t}_{\text{max}}^{q}$ and $\hat{t}_{\text{max}}^{f_k,g_k}$ can be defined from respective series and the lowest value must be considered. From a practical point of view, the forcing signal has smoother temporal variations than the structural response. Thus, the continuation technique is performed on the lowest value which is, in general, $\hat{t}_{\text{max}}^{q}$.

## 2. Numerical studies

In this section, two examples are considered. The first case-study is the well-known transient response of an elastic bar excited by an end force. Exact solution can be found in [4, 40]. The second case-study is related to transient dynamics of a cylindrical roof submitted to a point dead force. This problem can be considered as a reference one which was deeply examined [17, 41]. Numerical results are obtained with a spatial finite element discretization and the ANM based time integration procedure, previously introduced.

### 2.1. Elastic wave propagation in a rod

The model is considered as a rod of length $L = 1$ m which is clamped at its left hand-side on $x = 0$ m while its extremity on $x = L$ m is submitted to an Heaviside force $f(t) = F_0$ (for $t \ge 0$). Force intensity and cross-section of the rod are set to $F_0 = 1$ N and $A = 0.01$ m$^2$, respectively. A theoretical elastic material is defined according to Young modulus $E = 100$ Pa and mass density $\rho = 1$ kg.m$^{-3}$. Thus, longitudinal wave celerity $c_0 = \sqrt{E/\rho} = 10$ m.s$^{-1}$ is defined. Duration of simulation $T_{\text{max}} = 0.8$ s corresponds to the time required for both incident and reflected waves to travel 4 round trips at speed $c_0$. Four tracers are defined along the rod which abscissa are $x_1 = 0.25$ m, $x_2 = 0.5$ m, $x_3 = 0.75$ m and $x_4 = 1.00$ m (rod tip). The finite element model is composed of $N_{\text{FE}} = 20$ elements which are based on linear Lagrange shape function. For this study, the lumped mass matrix is adopted.

For comparison purpose, following author's write convention [4], the ANM based solver is compared to the Central Difference scheme (CD scheme, $\gamma = 1/2, \beta = 0$), the Purely Explicit scheme (PE scheme, $\gamma = \beta = 0$) and the “Average constant acceleration” scheme from the Implicit Newmark schemes family ($\gamma = 1/2, \beta = 1/4$). Constant time-step $\Delta t$ is also adopted and calculated as a percentage of the critical time-step $\Delta t_{\text{cr}}$. This latter corresponds to the time required for a longitudinal wave to travel through an element of length $L/N_{\text{FE}}$ at the constant wave propagation speed $c_0$. In this study, $\Delta t_{\text{cr}} = 5 \times 10^{-3}$ s.

In Figure 1, transient responses at locations $x_1$ and $x_4$ and obtained with the AMN solver and the CD-scheme are compared with the exact solution. Phenomena of wave displacement, reflection and superposition are observed. Because the mass matrix is lumped, the CD scheme allows to recover the exact solution for $\Delta t = \Delta t_{\text{cr}}$. In this figure, CD solution is computed with a ratio $\Delta t/\Delta t_{\text{cr}} = 0.99$ and exhibits the wave propagation. Nor attenuation neither periodicity error can be observed. Only sharp ripples of small amplitude can be noticed due to a time-step
<figure><img src="image_5.png" /></figure>

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

**Figure 1.** Computed displacement *u*(*t*) with ANM solver. Comparisons with the CD solver based solution and theoretical one [4, 42]. Study of the transient rod response.

different from the critical one. The ANM solution is obtained with a truncation order $N = 10$ and a continuation parameter $\delta = 10^{-8}$. Results show numerical solutions which preserve the wave propagation nature. But an attenuation of wave amplitudes, a periodicity error and oscillations are also noticed.

In order to better understand the ANM solver, comparisons are made with schemes relied on a Pure Explicit method and an implicit Newmark $\beta$-method.

In Figure 2, solutions measured in $x_1 = 0.25$ m and obtained with the ANM solver are compared with the ones computed with the PE solver. Five different values of the time-step $\Delta t$ varying from 20% to 0.1% of the critical time-step $\Delta t_{cr}$ are used to compute PE solutions. The lower time-step is, the lesser oscillatory the solutions are. Interestingly, the ANM solution tends to the PE one obtained with the lowest value of $\Delta t$. The ANM solution also presents the smallest amplitude of oscillations. In this simulation, values of the ANM parameters are identical: $\delta = 10^{-8}$ and $N = 10$ leading to 348 continuation steps and a mean value $< \hat{t}_{max} > = 2.3044 \times 10^{-3}$ s for the range of validity. This value can be compared to the value $\Delta t = 5.000 \times 10^{-6}$ s used for the finest PE simulation. The difference is of three orders magnitude.
<figure><img src="image_7.png" /><figcaption>(a) Displacement at location $x_1 = 0.25$ m.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b) Enlarged vision on the displacement at location $x_1 = 0.25$ m.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(c) Enlarged vision on the displacement at location $x_1 = 0.25$ m.</figcaption></figure>

<caption>**Figure 3.** Computed displacement $u(t)$ with ANM solver. Comparisons with solutions obtained with the mean acceleration scheme from implicit Newmark $\beta$-method ($\beta = 1/4$, $\gamma = 1/2$) with different values of time-step $\Delta t$. Study of the transient rod response.</caption>
<figure><img src="image_10.png" /><figcaption>(a) Displacement at location $x_1 = 0.25$ m.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(b) Enlarged vision on the displacement at location $x_1 = 0.25$ m.</figcaption></figure>

**Figure 4.** Computed displacement $u(t)$ with ANM solver. Truncation order set to $N = 10$. Influence of the ANM user-defined parameter $\delta$ (eq. (6)). Study of the transient rod response.

<table><caption>**Table 1.** Influence of the parameter $\delta$ on the ANM precision. The truncation order is set to $N = 10$. Study of the transient rod response.</caption><tr><td>δ parameter (see Eq. (6))</td><td>10<sup>-03</sup></td><td>10<sup>-04</sup></td><td>10<sup>-05</sup></td><td>10<sup>-06</sup></td><td>10<sup>-08</sup></td><td>10<sup>-10</sup></td></tr><tr><td>number of ANM steps</td><td>86</td><td>125</td><td>162</td><td>209</td><td>348</td><td>579</td></tr><tr><td>mean value &lt; t̂<sub>max</sub> &gt; (×10<sup>-3</sup> s)</td><td>9.3272</td><td>6.4303</td><td>4.9521</td><td>3.8389</td><td>2.3044</td><td>1.3829</td></tr></table>

In previous simulations, the user-defined ANM parameter $\delta$ is set to $10^{-8}$. In Figure 4, transient responses measured at location $x_1$ for the values $\delta = 10^{-10}, 10^{-8}, 10^{-6}, 10^{-5}, 10^{-4}$ and $10^{-3}$ are shown. For value of $\delta \ge 10^{-4}$, differences are noticeable. The wave propagation nature is still present but numerical oscillations are attenuated. For smaller values, no difference is observed. Thus, minimizing the $\delta$ value ensures a very high level of precision for the high-order prediction. The price to pay is decrease of the range of validity $\hat{t}_{max}$ and increase of ANM continuation steps. In Table 1, influence of this parameter $\delta$ is highlighted. The lower the value of $\delta$ is, the larger the number of ANM steps and the mean value $\langle \hat{t}_{max} \rangle$ are. For this example, the calculation with
$\delta = 10^{-5}$ required 162 continuation steps with an average value $\langle \hat{t}_{\text{max}} \rangle = 4.9521 \times 10^{-3}$ instead of 579 ANM steps and a mean value $\langle \hat{t}_{\text{max}} \rangle = 1.3829 \times 10^{-3}$ s with $\delta = 10^{-10}$. This behaviour is in agreement with past observations for the non-linear buckling analysis with ANM [26, 43, 44].

<figure><img src="image_12.png" /><figcaption>(a) Displacement at location $x_1 = 0.25$ m.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(b) Enlarged vision on the displacement at location $x_1 = 0.25$ m.</figcaption></figure>

**Figure 5.** Computed displacement $u(t)$ with ANM solver. Parameter $\delta$ set to $\delta = 10^{-5}$. Influence of the truncation order $N$ (for clarity purpose, for $N > 5$, only continuation points are plotted). Study of the transient rod response.
<table><caption>**Table 2.** Influence of the truncation order *N* on the ANM precision with value *δ* = 10⁻⁵. Study of the transient rod response.</caption><thead><tr><th>Truncation order <i>N</i></th><th>5</th><th>10</th><th>15</th><th>20</th><th>30</th></tr></thead><tbody><tr><td>number of ANM steps</td><td>916</td><td>162</td><td>83</td><td>56</td><td>34</td></tr><tr><td>mean value &lt; t&#x0302;<sub>max</sub> &gt; (×10<sup>-3</sup> s)</td><td>0.8742</td><td>4.9521</td><td>9.699</td><td>14.51</td><td>23.762</td></tr></tbody></table>

In the same way, influence of the truncation order *N* is studied. In Figure 5, transient responses measured at location *x*₁ are computed with the values *N* = 5, 10, 15, 20 and 30. In this figure, for convenience, only continuation points are plotted for the responses computed with *N* > 5. The wave propagation is well described and there is no observable difference between numerical predictions. Only very minor gaps can be detected: attenuation of the wave amplitude can be noticed with the smallest values of *N*. Data are summarized in Table 2. It is noticed that the higher truncation order is, the smaller number of continuation steps are and the greater mean value of *<t&#x0302;<sub>max</sub>>* is.

In the light of these previous analysis, it can be postulated that the explicit ANM solver allow to compute solution possessing the same properties of the Newmark solution. This later is unconditionally stable, without numerical damping and introducing periodicity error (as analysed in [4] for the Newmark integration methods). Some improvement can be reached with a pertinent choice of ANM parameter values. Increasing the truncation order and decreasing value of the *δ* parameter allow to decrease the number of continuation steps. Consequently, it is observed an increase of the mean value of *<t&#x0302;<sub>max</sub>>* which is greater than classical time-step value, outperforming the stability criterion associated with explicit time-integration schemes.

## 2.2. Non-linear elastic cylindrical shell submitted to a dead force

The second example concerns the Snap-through of a cylindrical shell whose geometry is given in Figure 6a. The material properties are *E*=2.1 10¹¹ N/m², *ν*= 0.25 and *ρ*= 1000 kg/m³. This example is the same that in Refs. [41] and [17].

The time evolution of the load at shell center is shown in Figure 6b. Due to symmetries, only a quarter of the shell is discretized with 16 shell elements based on geometrically exact element [45]. It is a quadrilateral eight nodes element with six degrees of freedom per node. Details relative to this finite element can be found in Ref. [45, 46]. This element is well-adapted to the perturbation method as shown in Ref. [39] and is the same element which was used in Ref. [17, 41] for this numerical example.

Time evolution of the vertical displacement of point M (see Figure 6a) previously computed in Ref. [17, 41] is plotted in Figure 7. In Ref. [41], Generalized Energy-Momentum Method was used to solve this non-linear dynamics shell problem with a time step equal to 1 ms. In [17], authors associated the Generalized Energy-Momentum Method to homotopy and perturbation methods with a time-step value set to 0.5 ms. The two computed time evolutions of the displacement were exactly the same as reported in Figure 7.
<figure><img src="image_14.png" /><figcaption>**Figure 6.** Geometry and load definitions of the cylindrical shell. Only 1/4 of the roof is discretized, with ρ=5m, L=2.5m, h=0.1m and θ = 30°. The nodal force is applied at the center of the cylindrical shell.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 7.** Time evolution of the displacement (Point M, Figure 6a) according to solutions given in Baguet (2011) [17] and Kuhl (1999) [41].</figcaption></figure>
Results obtained with the time perturbation method are now considered. For the following computations with ANM, the truncation order N and the prescribed tolerance $\delta$ (Eq. (6)) are equal to 20 and $10^{-6}$ respectively. These values have been chosen because on the one hand they lead to a good compromise between the step length (i.e. $\hat{t}_{\max}$) and the computation time of the right hand side and on the other hand they ensure a good accuracy of the solution (see Ref. [47]). The first difficulty with this approach is to correctly represent time evolution of the load which presents a non-smooth evolution at instant 0.2 s. Value of the regularization parameter $\eta$, introduced in Eq. (9), has to be well chosen to give accurate results. So, in Figure 8a, time evolution of the modified load (Eq. (9)) are plotted for four values of $\eta$ between $10^{-8}$ and $10^{-11}$. According to these results, it seems that a value of this regularization parameter lower than $10^{-10}$ should be chosen to correctly represent load evolution and its singularity at instant 0.2 s.

<figure><img src="image_16.png" /><figcaption>**Figure 8.** Influence of the regularization parameter $\eta$ (Eq. (9)) on the applied load and on the response curve. The solution obtained with the perturbation method is carried out with the following parameters : $N=20$ and $\delta = 10^{-6}$.</figcaption></figure>

This is confirmed by results shown in Figure 8b where time evolutions of the vertical displacement at the centre of the shell are plotted for seven values of the regularization parameter $\eta$. In this figure, it is observed that values of parameter $\eta$ greater than $10^{-8}$ do not permit to compute a solution with the correct dynamic buckling of the shell. Additionally, it seems that only values lower than $10^{-10}$ allow to compute the right instant for which buckling is triggered.

So in the sequel of the numerical study, value of this regularization parameter is set to $\eta = 10^{-15}$. This extra-small value has a great influence on the validity range (i.e. $\hat{t}_{\max}^{f_k,g_k}$) of power series expansion of the load expression. Nevertheless, as explained above, step length of the continuation method is governed by two ranges of validity $\hat{t}_{\max}$ : one for the unknowns $q$, noted $\hat{q}_{\max}$, and an other one for the load, noted $\hat{t}_{\max}^{f_k,g_k}$. As this latter, for all the numerical tests carried out in this study, is greater than the one for the unknowns (even in the case where $\eta = 10^{-15}$), this small value of the regularization parameter does not affect numerical simulations.

In Figure 8b, it is noticed that the time for which buckling appears obtained with the proposed method is lower than the ones obtained in Ref. [17, 41]. This difference is also illustrated in Figure 9, where solution obtained with Abaqus explicit software [48] is compared with solutions computed in Ref. [41] and by the proposed method. Computation with Abaqus is carried out by using 16 shell elements S4R (a 4 nodes doubly curved shell with reduced integration, hourglass control and finite membrane strains). In Figure 9, one can remark that the three numerical methods give three different instants for the panel buckling.
<figure><img src="image_17.png" /><figcaption>**Figure 9.** Comparison of time evolutions of the displacement (Point M, Figure 6a) to solutions given in Ref. [41], with Abaqus Explicit software [48] and with the proposed method with $\eta = 10^{-15}$, N=20 and $\delta = 10^{-6}$. For all the computations, 16 finite elements are considered.</figcaption></figure>

To explain this difference, several computations have been performed with Abaqus for a different number of shell elements (varying from 25 to 506). Results are exposed in Figure 10. From this plot, it is observed that solutions computed with Abaqus converge towards the one obtained with the perturbation method when number of finite elements is greater than 300. So, the ANM based solution seems to be the correct one. Let recall that this solution is obtained with only 16 shell elements which are also used in Ref. [41]. Consequently, the difference between solutions obtained with perturbation method and those of Ref. [41] is only due to the time integration method. One of advantages of a perturbation method is that the computed prediction is analytical, in time in this study. This property can explain the fact that, even with a small number of elements, the asymptotic expansions follow the correct solution emerging from the model. On the contrary, classical time integration schemes iteratively proceed with small time-steps introducing errors at each iteration. This can explain the difference on the buckling trigger instant between Abaqus predictions and ANM based ones.

Moreover, as shown in Figure 11, the solution obtained with the perturbation method is not sensitive to the number of finite elements considered. Indeed, the buckling trigger instant is always the same even if the number of elements is five times greater that the initial mesh.

It is known that classical time-integration methods also suffer from mesh dependency: CFL condition for explicit scheme (see [4, 42]) and spectral conditioning number associated to linear systems after spatial and time discretization for implicit scheme see, for example, analysis [9] for hyperbolic problem and [49] for parabolic problems).

To illustrate this drawback, value of the Abaqus time increment obtained with two different meshes is added on Figure 12a. This time increment is approximately equal to $4.5 \cdot 10^{-5}$ s and $2.8 \cdot 10^{-5}$ s for 100 and 306 elements, respectively. On this figure, time evolution of the parameter $\hat{t}_{max}$ (Eq. 6), is also plotted for three distinct meshing, with an increasing number of element ranging from 16 to 100. Let recall this parameter $\hat{t}_{max}$ evaluates the range of validity of asymptotic expansions. It is observed that $\hat{t}_{max}$ value is approximately the same for the three meshing. It coincides with the time-step value $4.5 \cdot 10^{-5}$ s used in the Abaqus simulations with 100 elements.
<figure><img src="image_18.png" /><figcaption>**Figure 10.** Comparison of the evolution of the displacement (Point M, Figure 6a) versus time according to the solutions given in Ref. [41], with Abaqus Explicit software [48] with several number of elements and with the proposed method with $\eta = 10^{-15}$, N=20 and $\delta = 10^{-6}$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 11.** Comparison of the evolution of the displacement (Point M, Figure 6a) versus time according to the solutions given with the Perturbation method with several number of elements with $\eta = 10^{-15}$, N=20 and $\delta = 10^{-6}$.</figcaption></figure>

Even if the value of parameter $\hat{t}_{max}$ evolves a lot during computations, it seems do not depend on the mesh size as the time-step does with classical time integration schemes.

Finally, it is proposed to evaluate influence of the truncation order *N* of asymptotic expansions to get the whole solution up to time *t* = *T*<sub>max</sub> = 0.3 s. So, in Figure 12b, the number of continuation steps is plotted with respect to the parameter *N*. It is also reported the total number of time
<figure><img src="image_20.png" /><figcaption>(a) Evolution of the parameter $\hat{t}_{\max}$ (Eq. (6)), $N=20$.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(b) Evolution of the number of steps.</figcaption></figure>

<caption>**Figure 12.** Evolutions of parameter $\hat{t}_{\max}$ (Eq. (6)) and number of ANM continuation steps with respect to the truncation order $N$. Continuation is performed until the end of simulation set to $T_{\max} = 0.3$ s. ANM parameter values are: $\eta = 10^{-15}$ and $\delta = 10^{-6}$.</caption>

iterations required with Abaqus to perform the same simulation, with meshing involving 306 and 506 finite elements. From these results, it is noticed that an increase of the truncation order leads to a decrease of the number of continuation steps. This is mainly due to the fact the validity range of the time, parameter $\hat{t}_{\max}$, depends on the truncation order. Greater is this latter, greater is the parameter $\hat{t}_{\max}$. This has been highlighted in the past when using asymptotic expansions in computational solid mechanics [26] or fluid dynamics [27].

In this dynamic buckling analysis, the perturbation method, based on power series of time, is compared with classical time integration methods available in Abaqus software. Computations with the perturbation method are carried out with 16 shell finite elements. To get same coherent results with Abaqus, it is required more than 300 finite elements (see Figure 10). This implies very small time-step values and large number of iterations. In that case, truncation order $N$ has to be greater than 10 to compute the solution with approximately the same number of continuation steps. As with a classical explicit time integration scheme, only a single matrix triangulation is needed for all the continuation step of the perturbation method. Nevertheless, each continuation step requires a number of backward/forward substitutions equals to $N-1$. The higher the truncation order is, the higher the computational time is. This can lead to a huge increase of the computational time especially for models involving a high number of degrees of freedom. One remedy, already used in computational mechanics, relies on the use of Padé approximants [21, 50]. In these past studies, a reduction of about approximately half the number of continuation step was observed. The Padé approximants have been tested in this study but results were not convincing. The reason seems mainly due to numerical instability issues. These instabilities appear in the computation of the Euclidian norm of vectors $q_i$ which dramatically increases with the truncation order $N$ and leads to numerical difficulties.

### 3. Conclusion

In this study, it was proposed a specific time-integration method in order to solve transient problems encountered in structural dynamics.

The new method relies on a high-order perturbation method in which the time variable is chosen as the perturbation parameter. Based on the use of a regularization parameter, it was shown how to recast equations in order to deal with arbitrary time variable external forcing.
Two classical problems were addressed. The first problem was the wave propagation in an elastic rod submitted to a constant end force. The second problem was the dynamic buckling of a cylindrical panel.

For both problems, results obtained with the ANM were analysed in the light of available results. For the first problem, these later were obtained with classical time integration schemes (Pure explicit method, Newmark $\beta$-method and Central Difference method). For the second problem, comparisons were carried out with results available in the literature and computed with the commercial Abaqus software.

For both problems, it was shown that explicit ANM predictions were very coherent. For the first problem, the wave propagation nature of the solution is preserved. For same meshing, the ANM solution is close to the implicit Newmark based solution computed with the smallest time-step value.

For the second problem, influence of the regularization parameter in the external forcing was established and it was demonstrated that a very small value was necessary. ANM solution also described the non-linear dynamics of the panel. Trigger instant of the buckling and post-bifurcated displacement were correctly obtained. Surprisingly, this study allowed for questioning some previous well-established results. More importantly, it was also demonstrated that the ANM solution is very similar to the Abaqus one which was based on a very refined meshing.

Finally, it was demonstrated that the ANM can be used as a pure explicit integration scheme without suffering from a stability criterion of CFL kind. Under the hypothesis of a correct choice of the continuation tolerance $\delta$ and regularization parameter $\eta$, it has seemed that the ANM allowed to intrinsically compute the right solution. To get solutions with the same level of precision, the alternative methods have been used with very fine meshing and very small time-steps.

However, in order the ANM based solver to become competitive, further developments are still necessary. As explained, the ANM scheme required high value of truncation order. But, the higher truncation order is, the higher computational time is, becoming harmful for models involving great number of unknowns. The use of a convergence accelerator could be a traditional solution. Unfortunately first attempt with Padé approximants failed. So, alternative solutions have to be sought. In order to enhance the range of validity of series and to decrease number of continuation steps, it could be envisaged use of the generalized factorial series or the Borel resummation technique [16, 23, 51]. Alternatively, use of special functions, as generalized hypergeometric Meijer G-function [24], could also be tested.

## Conflicts of interest

The authors have no conflict of interest to declare.

## References

[1] N. M. Newmark, "A method of computation for structural dynamics", *J. Eng. Mech. Div.* **85** (1959), no. 3, p. 67-94.

[2] T. J. R. Hughes, "Analysis of transient algorithms with particular reference to stability behavior", in *Computational methods for transient analysis*, Mechanics and Mathematical Methods. First Series: Computational Methods in Mechanics, vol. 1, North-Holland, 1983, p. 67-155.

[3] J. Chung, G. M. Hulbert, "A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized-$\alpha$ method", *J. Appl. Mech.* **60** (1993), no. 2, p. 371-375.

[4] M. Géradin, D. J. Rixen, *Mechanical vibrations: Theory and Application to Structural Dynamics*, John Wiley & Sons, 2014.

[5] D. Kuhl, M. A. Crisfield, "Energy-conserving and decaying algorithms in non-linear structural dynamics", *Int. J. Numer. Methods Eng.* **45** (1999), no. 5, p. 569-599.
[6] N. Gunwoo, S. Ham, K.-J. Bathe, "Performance of an implicit time integration scheme in the analysis of wave propagations", *Computers & Structures* **123** (2013), p. 93-105.

[7] N. Gunwoo, K.-J. Bathe, "The Bathe time integration method with controllable spectral radius: The $\rho_{\infty}$-Bathe method", *Computers & Structures* **212** (2019), p. 299-310.

[8] S.-B. Kwon, K.-J. Bathe, N. Gunwoo, "An analysis of implicit time integration schemes for wave propagations", *Computers & Structures* **230** (2020), article no. 106188.

[9] M. M. Malakiyeh, S. Shojaee, S. Hamzehei-Javaran, K.-J. Bathe, "New insights into the $\beta1/\beta2$-Bathe time integration scheme when L-stable", *Computers & Structures* **245** (2021), article no. 106433.

[10] D. Doyen, A. Ern, S. Piperno, "Time-integration schemes for the finite element dynamic Signorini problem", *SIAM J. Sci. Comput.* **33** (2011), no. 1, p. 223-249.

[11] V. Acary, "Projected event-capturing time-stepping schemes for nonsmooth mechanical systems with unilateral contact and Coulomb's friction", *Comput. Methods Appl. Mech. Eng.* **256** (2013), p. 224-250.

[12] O. Brüls, V. Acary, A. Cardona, "Simultaneous enforcement of constraints at position and velocity levels in the nonsmooth generalized-$\alpha$ scheme", *Comput. Methods Appl. Mech. Eng.* **281** (2014), p. 131-161.

[13] A. E. Aboanber, Y. M. Hamada, "PWS: an efficient code system for solving space-independent nuclear reactor dynamics", *Ann. Nucl. Energy* **29** (2002), no. 18, p. 2159-2172.

[14] A. E. Aboanber, Y. M. Hamada, "Power series solution (PWS) of nuclear reactor dynamics with Newtonian temperature feedback", *Ann. Nucl. Energy* **30** (2003), no. 10, p. 1111-1122.

[15] L. Guillot, B. Cochelin, C. Vergez, "A Taylor series-based continuation method for solutions of dynamical systems", *Dyn. Syst.* **98** (2019), no. 4, p. 2827-2845.

[16] C. Tayeh, G. Girault, Y. Guevel, J.-M. Cadou, "Numerical time perturbation and resummation methods for nonlinear ODE", *Nonlinear Dyn.* **103** (2021), no. 1, p. 617-642.

[17] S. Baguet, B. Cochelin, "Méthode Asymptotique Numérique adaptative pour la dynamique transitoire non-linéaire", in *10e colloque national en calcul des structures, Giens, France*, 2011, hal-00592804.

[18] N. Berrahma-Chekroun, M. Fafard, J. J. Gervais, "Resolution of the transient dynamic problem with arbitrary loading using the asymptotic method", *J. Sound Vib.* **243** (2001), no. 3, p. 475-501.

[19] M. Fafard, K. Henchi, G. Gendron, S. Ammar, "Application of an asymptotic method to transient dynamic problems", *J. Sound Vib.* **208** (1997), no. 1, p. 73-99.

[20] A. Deeb, A. Hamdouni, D. Razafindralandy, "Performance of Borel-Padé-Laplace integrator for the solution of stiff and non-stiff problems", *Appl. Math. Comput.* **426** (2022), article no. 127118.

[21] A. Elhage-Hussein, M. Potier-Ferry, N. Damil, "A numerical continuation method based on Padé approximants", *Int. J. Solids Struct.* **37** (2000), p. 6981-7001.

[22] D. Razafindralandy, A. Hamdouni, C. Allery, "Numerical divergent series resummation in fluid flow simulation", *Eur. J. Comput. Mech.* **17** (2008), no. 4, p. 431-451.

[23] D. Razafindralandy, A. Hamdouni, A. Deeb, "Considering inverse factorial series as time integration method", *AIP Conference Proceedings* **1798** (2017), article no. 020129.

[24] H. Mera, T. G. Pedersen, B. K. Nikolić, "Fast summation of divergent series and resurgent transseries from Meijer-G approximants", *Phys. Rev. D* **97** (2018), no. 10, article no. 105027.

[25] L. Azrar, B. Cochelin, N. Damil, M. Potier-Ferry, "An asymptotic-numerical method to compute the post-buckling behaviour of elastic plates and shells", *Int. J. Numer. Methods Eng.* **36** (1993), no. 8, p. 1251-1277.

[26] B. Cochelin, "A path-following technique via an asymptotic-numerical method", *Computers & Structures* **53** (1994), no. 5, p. 1181-1192.

[27] J.-M. Cadou, M. Potier-Ferry, B. Cochelin, "A numerical method for the computation of bifurcation points in fluid mechanics", *Eur. J. Mech. B Fluids* **25** (2006), p. 234-254.

[28] Y. Guevel, H. Boutyour, J.-M. Cadou, "Automatic detection and branch switching methods for steady bifurcation in fluid mechanics", *J. Comput. Phys.* **230** (2011), no. 9, p. 3614-3629.

[29] H. Boutyour, H. Zahrouni, M. Potier-Ferry, M. Boudi, "Bifurcation points and bifurcated branches by an asymptotic numerical method and Padé approximants", *Int. J. Numer. Methods Eng.* **60** (2004), p. 1987-2012.

[30] P. Vannucci, B. Cochelin, N. Damil, M. Potier-Ferry, "An Asymptotic-Numerical Method to compute bifurcating branches", *Int. J. Numer. Methods Eng.* **41** (1998), p. 1365-1389.

[31] F. Xu, M. Potier-Ferry, "A multi-scale modeling framework for instabilities of film/substrate systems", *J. Mech. Phys. Solids* **86** (2016), p. 150-172.

[32] A. Tri, H. Zahrouni, M. Potier-Ferry, "High order continuation algorithm and meshless procedures to solve nonlinear Poisson problems", *Eng. Anal. Bound. Elem.* **36** (2012), no. 11, p. 1705-1714.

[33] H. Tian, M. Potier-Ferry, F. Abed-Meraim, "A numerical method based on Taylor series for bifurcation analyses within Föppl-von Karman plate theory", *Mech. Res. Commun.* **93** (2018), p. 154-158.

[34] S. Mordane, "Calcul du problème de la houle non-linéaire et instationnaire par une méthode asymptotique numérique", PhD Thesis, Université HASSAN II, Faculté des sciences BEN M'SIK Casablanca, Morocco, 1995.
[35] B. Claude, L. Duigou, G. Girault, Y. Guevel, J.-M. Cadou, "Numerical comparison of eigenvalue algorithms for vibroacoustic problems", *Mech. Res. Commun.* **91** (2018), p. 39-45.

[36] M. Potier-Ferry, N. Damil, B. Braikat, J. Descamps, J.-M. Cadou, H. L. Cao, A. Elhage-Hussein, "Traitement des fortes non-linéarités par la méthode asymptotique numérique", *C. R. Acad. Sci. Paris* **324** (1997), no. 3, p. 171-177.

[37] L. Guillot, B. Cochelin, C. Vergez, "A generic and efficient Taylor series-based continuation method using a quadratic recast of smooth nonlinear systems", *Int. J. Numer. Methods Eng.* **119** (2019), no. 4, p. 261-280.

[38] H. Zahrouni, M. Potier-Ferry, H. Elasmar, N. Damil, "Asymptotic numerical method for nonlinear constitutive laws", *Revue Européenne des Éléments finis* **7** (1998), no. 7, p. 841-869.

[39] H. Abichou, H. Zahrouni, M. Potier-Ferry, "Asymptotic numerical method for problems coupling several nonlinearities", *Comput. Methods Appl. Mech. Eng.* **191** (2002), no. 51-52, p. 5795-5810.

[40] K. F. Graff, *Wave motion in elastic solids*, Courier Corporation, 2012.

[41] D. Kuhl, E. Ramm, "Generalized energy-momentum method for non-linear adaptive shell dynamics", *Comput. Methods Appl. Mech. Eng.* **178** (1999), no. 3-4, p. 343-366.

[42] E. Grosu, I. Harari, "Stability of semidiscrete formulations for elastodynamics at small time steps", *Finite Elem. Anal. Des.* **43** (2007), no. 6-7, p. 533-542.

[43] S. Baguet, B. Cochelin, "On the behaviour of the ANM continuation in the presence of bifurcations", *Commun. Numer. Methods Eng.* **19** (2003), no. 6, p. 45-471.

[44] B. Cochelin, M. Medale, "Power series analysis as a major breakthrough to improve the efficiency of asymptotic numerical method in the vicinity of bifurcations", *J. Comput. Phys.* **236** (2013), p. 594-607.

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

[46] H. Zahrouni, B. Cochelin, M. Potier-Ferry, "Computing finite rotations of shells by an asymptotic-numerical method", *Comput. Methods Appl. Mech. Eng.* **175** (1999), no. 1-2, p. 71-85.

[47] J.-M. Cadou, B. Cochelin, N. Damil, M. Potier-Ferry, "Asymptotic Numerical method for stationary Navier-Stokes equations and with Petrov-Galerkin formulation", *Int. J. Numer. Methods Eng.* **50** (2001), p. 825-845.

[48] K. Hibbitt, Sorenson, ABAQUS, "Theory and users manuals. version 5.8", Pawtucket, RI 02860-4847, USA.

[49] L. Zhu, Q. Du, "Mesh dependent stability and condition number estimates for finite element approximations of parabolic problems", *Math. Comput.* **83** (2014), no. 285, p. 37-64.

[50] G. A. Baker, P. Graves-Morris, *Padé approximants*, second ed., Encyclopedia of Mathematics and Its Applications, vol. 59, Cambridge University Press, 1996.

[51] D. Razafindralandy, A. Hamdouni, "Time integration algorithm based on divergent series resummation for ordinary and partial differential equations", *J. Comput. Phys.* **236** (2013), p. 56-73.