# A multiple-solution method for non-linear structural mechanics

Pierre-Alain Boucard, Pierre Ladevèze

► To cite this version:

Pierre-Alain Boucard, Pierre Ladevèze. A multiple-solution method for non-linear structural mechanics. NMCM'98 - 7th International Conference on Numerical Methods in Continuum Mechanics., Oct 1998, Stará Lesná, Slovakia. hal-01715170

HAL Id: hal-01715170

https://hal.science/hal-01715170

Submitted on 12 Oct 2020

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# A multiple-solution method for non-linear structural mechanics

PIERRE-ALAIN BOUCARD* - PIERRE LADEVEZE*

The numerical solution process for complex, time-dependent non-linear problems (material stamping, cyclic viscoplasticity, crashes, etc.) requires, even with the most powerful computers, a computation time which turns out to be prohibitive. Moreover, with respect to those structural analyses involving multiple solution sequences (parametric studies, probabilistic analyses, flawed structures, etc.), each series of data necessitates performing a full calculation, even if the problems to be solved do resemble one another. The goal of the work presented herein is to develop a strategy that is well-suited to multiple-solution problems. Such a strategy is to be based on the LATIN method and, more specifically, on its capacity to reuse the solution to a given problem in order to solve similar problems. The goal then is both to assess the major obstacles encountered during this computational procedure and to define and apply a strategy that enables minimizing total computing costs with respect to multiple-solution problems.

Keywords: structural mechanics, non-linear computations, multiple solution

## 1. INTRODUCTION

Whenever the analysis of a structure involves multiple solutions (parametric studies, probabilistic analyses, flawed structures, etc.), each set of data utilized necessitates a separate, full-scale calculation. As such, a good number of similar types of problems wind up getting solved. The choice of an efficient and well-adapted computational method is thus of vital importance. The LATIN method [5], [6] is non-incremental in nature and, consequently, would appear to represent a more promising approach, inasmuch as the

set of principles inherent in this method tend to be more applicable than are the majority of conventional incremental algorithms (stochastic finite element approaches, perturbation-related methods [1], [7], and Monte Carlo simulations [8]). The objective of the work presented herein is to demonstrate the potential of the LATIN method with respect to multiple-solution problems. For the special case of buckling, a strategy adapted to multiple solutions has been developed which relies on the LATIN method and, more specifically, on its capacity to reuse the solution to a given

<footnote>* LMT Cachan
(E.N.S. de Cachan / C.N.R.S. / Université Paris 6)
61 avenue du Président Wilson
94235 CACHAN CEDEX (FRANCE)
e-mail: Pierre-Alain.Boucard@lmt.ens-cachan.fr</footnote>
problem in order to solve similar problems.
A full presentation of the LATIN method
can be found in [5], [6]. The solution
process for buckling problems, along with
a set of pertinent examples, is provided in
[2], [3].

## 2. REVIEW OF THE LATIN METHOD

### 2.1 Notation

The structure being studied occupies a
domain Ω of the ordinary three-dimensional
Euclidean space, whose corresponding
vectorial space is denoted by E³. The
vectors of E³, displacement for instance,
are denoted by U. The current point of Ω is
noted M. The linear operators of E³, strain
and stress, are denoted by ε and σ,
respectively. These quantities, in a given
base, are represented by a matrix (3,3); this
notation is also used for the internal
variables. For more complex operators,
such as Hooke's elasticity tensor, they are
written as: K. To define the work, one uses
the trace operator denoted by "Tr", thus:

$$ \mathrm{Tr}[\sigma \varepsilon] = \sum_{i=1,2,3} [\sigma]_{ij} [\varepsilon]_{ij} $$

where $[\sigma]_{ij}$ and $[\varepsilon]_{ij}$ are the
components of the operators $\sigma$ and $\varepsilon$,
respectively.

### 2.2 Reference problem

For the whole time interval, denoted by
[0, T], the structure being studied is
subjected to surface loading F<sub>d</sub> on a part
∂₂Ω of the boundary, ∂Ω, of Ω, and to
volume loading f<sub>d</sub>. On the complementary
part, ∂₁Ω, to ∂Ω, the displacement U<sub>d</sub> is
assigned. The method is carried out by
means of small perturbations in order to
simplify the presentation of the concept.

This method has been conducted in [5] for
large transformations and has been
developed in [2], [3] for the analysis of
both buckling and post-buckling.

Let us consider small displacement
problems and quasi-static loading. U<sup>[0,T]</sup>,
S<sup>[0,T]</sup> denote the spaces where
displacement U and stress σ respectively,
defined on Ωx[0, T], are searched. The
subspaces of displacement kinematically
admissible for homogeneous equations are
such that:

$$ U_0 = \{ \underline{U}(\underline{M}) ; \underline{U}(\underline{M}) |_{\partial_1 \Omega} = 0 \} $$

$$ U_0^{[0,T]} = \{ \underline{U}(\underline{M},t) \in \underline{U}^{[0,T]} ; \forall t \in [0,T] \ \underline{U}(\underline{M},t) \in \underline{U}_0 \} $$

The material state is completely
determined by the values of strain, the
inelastic strain εP and additional internal
variables X. In general, X is a column of
R<sup>m</sup>. The conjugate quantity of X is denoted
by Y; it is also a column of R<sup>m</sup>.

Then, the structure's evolution is
described by the following problem:

Find U(M,t)) and σ(M,t)
M∈Ω, t∈[0,T] such that:

• kinematic constraints:

$$ \underline{U}(\underline{M},t) \in \underline{U}^{[0,T]}, \quad \underline{U}|_{\partial_1 \Omega} = \underline{U}_d $$

• equilibrium equation:

$$ \sigma \in S^{[0,T]}, \\ \forall \underline{U}^* \in \underline{U}_0^{[0,T]} \\ \int_{\Omega} d\Omega \mathrm{Tr}[\sigma \varepsilon(\underline{U}^*)] = \int_{\Omega} d\Omega \underline{f}_d \cdot \underline{U}^* \\ + \int_{\partial 2\Omega} dS \underline{F}_d \cdot \underline{U}^* $$

• state laws

$$ \dot{\varepsilon} = \dot{\varepsilon}^P + \dot{\varepsilon}^e, \quad \sigma = \mathbf{K}\varepsilon^e $$

$$ \underline{Y} = \underline{G}(\underline{X}) $$
• evolution laws

$$ \frac{d}{dt} \begin{bmatrix} \varepsilon \mathbf{P} \\ -\mathbf{X} \end{bmatrix} = \mathbf{B} \begin{bmatrix} \boldsymbol{\sigma} \\ \mathbf{Y} \end{bmatrix} \\ \varepsilon \mathbf{P}_{(t=0)} = 0 \quad \mathbf{Y}_{(t=0)} = 0 $$

K, G are given material characteristics. B is a given positive operator which depends on the material.

## 2.3 Principles of the method

Let us introduce $S = (\dot{\varepsilon}\mathbf{P}, \dot{\mathbf{X}}; \sigma, \mathbf{Y})$ for which each component is defined on $\Omega \times [0, T]$. The corresponding space is denoted by $S^{[0,T]}$. S is what we are looking for; it is defined by a set of equations. The principle of the method consists of constructing a sequence of elements $S_n$, defined over $\Omega \times [0,T]$. The entire loading is then handled in a single time interval. The method is based on the following three principles:

**P1:** the problem's equations are divided into two groups which serve to define two sets: $\Gamma$ and $A_d$, where:

• $A_d$: linear and ultimately global equations,

• $\Gamma$: local in space variable and ultimately non-linear equations

**P2:** The method then proceeds by means of a two-stage iterative sequence:

One iteration is composed of both a local stage, during which a solution belonging to $\Gamma$ is built, and a global stage during which a solution belonging to $A_d$ is sought. Thus, the solution to the problem being posed lies at the intersection of these two solution sets. In order to resolve these two stages, two search directions $\mathbf{E}^+$ and $\mathbf{E}^-$, parameters of the method, are to be

used.

**P3:** The problem's unknowns are represented by a product sum of both time and space functions.

To derive a formulation of the normal material (linear state laws, $\mathbf{Y} = \Lambda \mathbf{X}$), the only non-linear relations are the evolution laws; moreover, these laws are local in space variable.

In the case of normal material models, the splitting used is the following:

$$ (i) \quad \begin{cases} \text{kinematic constraints} \\ \text{equilibrium equations} (\forall t \in [0, T]) \\ \boldsymbol{\sigma} = \mathbf{K} (\varepsilon - \varepsilon_p), \mathbf{Y} = \Lambda \mathbf{X} \end{cases} $$

$$ (ii) \quad \begin{cases} \frac{d}{dt} \begin{bmatrix} \varepsilon \mathbf{P} \\ -\mathbf{X} \end{bmatrix} = \mathbf{B} \begin{bmatrix} \boldsymbol{\sigma} \\ \mathbf{Y} \end{bmatrix} \\ \varepsilon \mathbf{P}_{(t=0)} = 0 \text{ and } \mathbf{X}_{(t=0)} = 0 \end{cases} $$

The first group is associated with the free energy and the second one with the mechanical dissipation. Let $A_d$ be the linear subspace of S satisfying the first group of equations, i.e. the admissibility conditions and the state equations. The S satisfying the second group of equations, i.e. the evolution laws is a non-linear subspace denoted by $\Gamma$. So, the problem is to find the intersection of $A_d$ and $\Gamma$.

The second step of the LATIN method is to use a two-stage iterative scheme (see Figure 1).

The iterative procedure starts with an initial $S_0$ belonging to $A_d$, for example the elastic response of the structure over $[0, T]$. Then, successively, elements belonging to $\Gamma$ and to $A_d$ are built up until the practical convergence.
<figure><img src="image_2.png" /><figcaption>Figure 1: Iteration scheme</figcaption></figure>

Each iteration, i.e. the building of a new admissible element $S_{n+1}$ from a given admissible element $S_n$, has two stages:

* local stage:

find $\hat{S}_{n+1/2} \in \Gamma$ such that

$$ \hat{S}_{n+1/2} = (\Gamma) \cap (S_n + E^+) $$

* linear global stage:

find $S_{n+1} \in A_d$ such that

$$ S_{n+1} = (A_d) \cap (\hat{S}_{n+1/2} + E^-) $$

Where $E^+$ and $E^-$ are parameters. Practically, the search directions $E^+$ and $E^-$ are defined by two linear operators $H^+$ and $H^-$ such that:

$$ E^+ = \{ \Delta s / \frac{d}{dt} \begin{bmatrix} \Delta \epsilon^p \\ -\Delta \underline{X} \end{bmatrix} + H^+ \begin{bmatrix} \Delta \sigma \\ \Delta \underline{Y} \end{bmatrix} = 0 ; \\ \Delta s \in S^{[0,T]} \} $$

$$ E^- = \{ \Delta s / \frac{d}{dt} \begin{bmatrix} \Delta \epsilon^p \\ -\Delta \underline{X} \end{bmatrix} - H^- \begin{bmatrix} \Delta \sigma \\ \Delta \underline{Y} \end{bmatrix} = 0 ; \\ \Delta s \in S^{[0,T]} \} $$

$H^+$ and $H^-$ depend on the time and space variables.

The local stage then leads to a local non-linear, evolution problem in space variables. The resolution of this stage turns out to be rather straightforward. In contrast, the global stage requires the

solution to a global linear problem in space variable that's been parameterized by the time interval over which specific solution techniques are necessary in order to limit this particular stage's cost.

## 2.4 Processing of the global stage - Representation of the unknowns

The global stage consists of searching for a new element $S_{n+1}$ that's better than $S_n$; similarly, the difference between these two elements of $A_d$ could be chosen as an unknown during this stage. Through working on the corrections of the unknowns, two new quantities naturally get introduced:

$$ \Delta \hat{S} = \hat{S} - S_n = (\Delta \hat{\epsilon}^p, \Delta \hat{\underline{X}} ; \Delta \hat{\sigma}, \Delta \hat{\underline{Y}}) $$

and

$$ \Delta S_n = S_{n+1} - S_n = (\Delta \dot{\epsilon}_n^p, \Delta \dot{\underline{X}}_n, \Delta \sigma_n, \\ \Delta \underline{Y}_n) $$

It should be pointed out at this juncture that the element $S_{n+1}$ must belong to $A_d$ in order to remain within the framework of the LATIN method.

The direct solution to global linear problems which have been parameterized by time within the procedure of the global
stage can lead to sizable computing times. As a means of lowering these computing times, representations of both time and space unknowns are introduced in order to solve the global stage at a low computing cost. The separation of the functional dependencies in both time and space has yielded satisfactory results for quasi-static loadings [3], [5], [6]. The corrections are then sought by means of superimposing solutions of the simple-loading type. Such solutions are recognized as being good approximations of non-linear, quasi-static problem solutions.

The representation for each variable is as follows:

$$ \dot{\epsilon}_{n+1}^{P}(\underline{M}, t) - \dot{\epsilon}_{n}^{P}(\underline{M}, t) = \sum_{i=1}^{m} a_{i}(t) A_{i}(\underline{M}) $$

$$ \dot{\underline{X}}_{n+1}(\underline{M}, t) - \dot{\underline{X}}_{n}(\underline{M}, t) = \sum_{i=1}^{m} c_{i}(t) \underline{C}_{i}(\underline{M}) $$

$$ \sigma_{n+1}(\underline{M}, t) - \sigma_{n}(\underline{M}, t) = \sum_{i=1}^{m} b_{i}(t) \underline{B}_{i}(\underline{M}) $$

$$ \underline{Y}_{n+1}(\underline{M}, t) - \underline{Y}_{n}(\underline{M}, t) = \sum_{i=1}^{m} d_{i}(t) \underline{D}_{i}(\underline{M}) $$

In practice, m is limited to either 1 or 2. The time functions are scalar functions which have been entirely determined by the algorithm employed, with the finite element-type of representations being used for the discretization of the space functions. This representation is valuable in analyzing the results of non-linear simulations, especially if the solution in both time and space can be represented with a small number of function couples. In this respect, it becomes essential to seek the best function couples in order to represent the corrections in the aim of reducing the number of corrections. It should be

highlighted herein that for buckling problems, the approach adopted in [4], although different, is oriented along the same lines: just a small number of basic functions is indeed sufficient, provided the functions have been selected carefully.

## 2.5 Preliminary stage

This stage is to be performed prior to the global stage and consists of seeking to improve the solution without having to recalculate a new space function. Involved herein therefore is the efficient utilization of the existing basis of space functions in order to represent the corrections. The preliminary stage gets inserted into the algorithm in the following manner: in practice, two types of iterations can be encountered (see Figure 2):

- the preliminary stage doesn't provide for sufficient quality improvement, and the global stage must be therefore be solved so as to obtain new space functions, in which case the global stage entails costly computations. Furthermore, it becomes necessary to orthogonalize those space functions serving as the basis herein.

- the preliminary stage does enable, by means of determining the corrections to the existing time functions, not having to generate a new space function: this iteration then winds up being rather inexpensive since it leads to solving a system of differential equations whose size is the same as the number of existing space functions. The coefficients of this system of equations are simply the integrals over the half-domain of known quantities.
<figure><img src="image_3.png" /><figcaption>Figure 2: Insertion of the preliminary stage</figcaption></figure>

### 3. MULTIPLE-SOLUTION METHOD

The LATIN method makes it possible to obtain a representation of the problem's solution in the form of a product sum of both time and space functions. As such, it could be considered that the method serves to automatically build an optimal basis for representing the solution. The notion herein then is to use this representative basis for determining the solution to a problem that's similar to the one for which this basis was built in the first place. In this

context, a similar problem would be one whereby an initial consideration would suggest that both the geometry and the limit conditions are identical (this requirement ensures respecting the "zero-admissibility" condition of the space functions), yet whereby the material characteristics can vary throughout the structure.

This procedure allows utilizing the fact that the algorithm can be initialized by any solution belonging to the set $\mathbf{A}$d (though usually an elastic solution). The notion herein therefore is to initialize the process
associated with the similar problem (the "flawed" structure) by the results of the computation carried out on the "flawless" structure. In this manner, a basis of space functions with a strong mechanical content is available immediately as of the first iteration, thereby making it possible to represent the solution.

In this instance, the preliminary stage plays a vital role: it enables both verifying that the basis of the space functions is well-adapted to the target problem and searching for new time functions which allow deriving the solution to the "flawed" problem. In the best-case scenario, when the basis is sufficient, no new space function will be generated and the problem's solution will thus be obtained at low cost. The alternative would be the automatic computation of new space functions in order to enhance the initial basis: the number of such functions however must remain less than the number in the initial basis in order to ensure the method's efficiency. Should the solutions to both the "flawed" and "flawless" problems be close enough, the solution to the latter problem can then be derived at a significantly lower cost than that of a full-scale computation.

## 4. EXAMPLES

The selected scope of application for testing this method is the buckling of both beams and plane elastic arcs. Additional details on the formulation used herein can be found in [2], [3]. The first example pertains to a straight beam embedded at one extremity and subjected to a displacement

being applied at the other extremity. The structural flaws introduced consist of variations of between -60% and +100% in the Young's modulus on different elements. The influence of a particular flaw on the value of the critical buckling load ($F_c/F_{c0}$ ratio) can thereby be examined. The results have been presented in the two figures: Figure 3 and 4.

A similar computation was carried out on a half-embedded bar; the results are presented in Figure 5 and 6.

For both of the cases analyzed, figures 4 and 6 display the number of space functions added during the computations at the level of the initial basis (6 time-space function couples). This number provides an indicator of the total computing cost inasmuch as this phase constitutes the algorithm's most costly stage. It can be observed that for the majority of cases processed, at most one space function gets added. The basis of initial functions therefore enables conducting many computations "with flaws" at a much lower cost than that of a full-scale computation: within the cases presented herein, the computing time necessary for obtaining the solution on the "flawed" bar is between 10% and 20% of that associated with a full-scale computation.

## 5. CONCLUSION

These initial numerical examples have exhibited the algorithm's very satisfactory behavior when applied to the case of multiple solutions in the analysis of buckling. The basis of space functions generated by an initial solution on the
<figure><img src="image_4.png" /><figcaption>Figure 3: Influence of a defect</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Number of space functions generated</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: Influence of a defect on the critical load</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Number of space functions generated</figcaption></figure>
perfect structure has been very well-adapted to the set of computations conducted on "flawed" structures, provided that these computations don't exert excessive perturbations on the response. Moreover, the preliminary stage makes it possible to fully utilize the resources contained within this basis, thereby ensuring greatly reduced computing costs.

This initial study has provided the set of elements which enable, for both a given load and a given distribution of flaws, determining the probability of the

structure's failure by means of buckling [1], [7]. Lastly, the approach employed is quite general in nature and should be applicable to a good number of other non-linear problems which require multiple solutions.

An abbreviated version of this paper was presented at the VII conference Numerical Methods in Continuum Mechanics, High Tatras, Stara Lesna, Slovakia, October 6-9, 1998 and published in its Proceedings.

## 6. REFERENCES

[1] BENAROYA H., REHAK M., "Finite element methods in probabilistic structural analysis: A selective review", Applied Mechanics Reviews, Vol. 41, No. 5, pp. 201-213, 1988.

[2] BOUCARD P.A., LADEVÈZE P., POSS M., ROUGÉE P., "A non-incremental approach for large displacement problems", Computers and Structures Vol. 64, 1-4, pp. 449-508, 1997.

[3] BOUCARD P.A., LADEVÈZE P., "A Non-Incremental Approach for Post-Critical Analysis", Numerical Methods in Engineering 96, J.-A. Désidéri et al. Eds., J. Wiley, pp. 164-170, 1996.

[4] COCHELIN B., AZRAR L., DAMIL N., POTIER-FERRY M., "An asymptotic-numerical method to compute the post-buckling beahviour of elastic plates and shells", International Journal for Numerical Methods in Engineering, Vol. 36, 1251-1277, (1993).

[5] LADEVÈZE P., "Non Linear Mechanics of Structures - New approaches and non-incremental methods of calculation", Springer, 1998, also available in French (Hermès, 1996).

[6] LADEVÈZE P., New advances in the LArge Time INcrement method, New Advances in Computational Structural Mechanics, Ladevèze P., Zienkiewicz O.C. Eds., pp. 3-21, 1991.

[7] MACIAS O.F. - LEMAIRE M., "Eléments Finis stochastiques et Fiabilité Application en mécanique de la rupture", Revue française de génie civil, Vol. 1 - No. 2, 1997.

[8] SHINOZUKA M., ASTILL J., "Random eigenvalue problems in structural mechanics", AIAA Journal, Vol. 10, No. 4, pp. 456-462, 1972.