Lionel Gendre
email: gendre@lmt.ens-cachan.fr

Olivier Allix
email: allix@lmt.ens-cachan.fr

Pierre Gosselet
email: gosselet@lmt.ens-cachan.fr

François Comte

Non-intrusive and exact global/local techniques for structural problems with local plasticity

Keywords: finite element, elastic-plastic, successive approximations, nonlinear

This paper introduces a computational strategy to solve structural problems featuring nonlinear phenomena that occur within a small area, while the rest of the structure retains a linear elastic behavior. Two finite element models are defined: a global linear model of the whole structure, and a local nonlinear "submodel" meant to replace the global model in the nonlinear area. An iterative coupling technique is then used to perform this replacement in an exact but non-intrusive way, which means the model data sets are never modified and the computations can be carried out with standard finite element software. Several ways of exchanging data between the models are discussed and their convergence properties are investigated on two examples.

Introduction

In the aircraft industry, it is a common task to perform a finite element (FE) analysis on a complex structure that mostly evolves in a linear elastic way, but exhibits confined plasticity (or other nonlinear phenomena) in a small critical region. In most FE software, such an analysis is usually carried out by dividing the loading history into several load increments, and by solving nonlinear equilibrium equations at each increment, using Newton's method or one of its variants. When the problem size is too large or the loading history is too complex, this approach can lead to unaffordable computational costs; in such cases, two categories of dedicated techniques can be used, taking advantage of the small extent of the nonlinear area.

The first class of techniques is called submodeling and consists in a global linear computation followed by a local nonlinear "zoom" centered on the critical zone and driven by the global displacements [START_REF] Kelley | Mesh requirements for the analysis of a stress concentration by the specified boundary displacement method[END_REF][START_REF] Cormier | Aggressive submodelling of stress concentrations[END_REF][START_REF] Voleti | Global-local analysis of large-scale composite structures using finite element methods[END_REF][START_REF] Srinivasan | Identifying global/local interface boundaries using an objective search method[END_REF] or stresses [START_REF] Jara-Almonte | The specified boundary stiffness and force (SBSF) method for finite element subregion analysis[END_REF]. Although widely available and numerically efficient, those methods suffer from a strong limitation: they usually consist in a single one-way data transfer that ignores the global influence of local plastic-ity. As a consequence, they cannot assess phenomena such as stress redistributions, and tend to introduce uncontrolled errors due to the inaccuracy of local boundary conditions. To overcome this problem, some authors have suggested adding a global correction step to obtain a full iterative procedure [START_REF] Mao | A refined global-local finite element analysis method[END_REF][START_REF] Whitcomb | Iterative global-local finite element analysis[END_REF] or using static condensation to improve global results [START_REF] Hirai | An efficient zooming method for finite element analysis[END_REF][START_REF] Hirai | An exact zooming method[END_REF]. However, most of those methods were designed as structural zooming or mesh refinement techniques for linear elasticity, and are not always suitable for nonlinear problems.

The second possibility is to use multiscale methods [START_REF] Mandel | Convergence of a balancing domain decomposition by constraints and energy minimization[END_REF][START_REF] Farhat | FETI-DP: a dual-primal unified FETI method. I. A faster alternative to the two-level FETI method[END_REF], some of which have been specifically designed to handle local nonlinearity: for example, the LaTIn micro/macro approach [START_REF] Ladevèze | A micro/macro approach for parallel computing of heterogeneous structures[END_REF][START_REF] Guidault | A two-scale approach with homogenization for the computation of cracked structures[END_REF] or nonlinear localization techniques [START_REF] Cresta | Nonlinear localization strategies for domain decomposition methods: application to post-buckling analyses[END_REF][START_REF] Pebrel | A nonlinear dual domain decomposition method: application to structural problems with damage[END_REF]. Those methods are robust and efficient, but they use complex formulations that are intrusive with respect to the traditional FE framework. This means they cannot be carried out using commercial FE software and industrial data sets, and it would be long and difficult to implement them into an industrial environment; therefore, they should be considered as a medium-term solution, rather than an immediately usable technique.

In this paper, we propose an intermediate way to analyze complex structures that contain small nonlinear areas, when full nonlinear computations are too expensive. This method is designed to be both robust and non-intrusive, which means it does not introduce any error based on approximate assumptions and is entirely based upon standard FE computations. Moreover, in addition to handling local nonlinearity, the local model can be used to introduce geometric details, mesh refinements or specific constitutive laws that are absent from the global model; it can also be analyzed using a separate piece of code, which may contain features that are not implemented in the global FE solver. In this sense, the approach also works as a flexible and exact structural reanalysis [START_REF] Akgün | Fast exact linear and nonlinear structural reanalysis and the Sherman-Morrison-Woodbury formulas[END_REF] and solver coupling technique.

Summary of the method

The analysis technique is based on three main ideas. First, the computations are performed on two models: a global linear elastic model representing the whole structure, and a local nonlinear "submodel" describing the nonlinear area only. The local model is usually defined after a first global linear analysis, by the means of an appropriate model error indicator that depends on the nonlinear constitutive law. This framework is similar to traditional submodeling techniques, with one major difference: should the local model contain many enhancements or changes with respect to the global model, this would not cause inaccuracies in the final result.

Second, the local model is used to replacerather than to enrich or to refine -a part of the global model. In other terms, the overlap of models is not used, except on the boundary of the nonlinear area; the goal of the strategy is to discard the inaccurate area of the global model and replace it with the local model. For that purpose, the global/local solution is defined as:

u = u local in
"#!$%&'%()*#) *+%)(#,-*.#!)#/

Fig. 1 Summary of the iterative exchange technique

This framework ensures that if the sequence of global/local solutions (as defined above) is convergent, then its limit is the "reference solution", ie. what would be obtained by actually performing the model substitution and running a full nonlinear analysis; in this sense the technique is exact, since the error can be reduced as much as needed. It is also non-intrusive because it only uses traditional FE analyses with widely available input and output (displacements, concentrated loads and reaction forces at nodes) and requires no modification of existing FE data sets.

From a practical point of view, this approach can have two advantages over a direct nonlinear analysis. First, it can be computationally more efficient, thanks to the local handling of nonlinearity (as explained in [START_REF] Cresta | Nonlinear localization strategies for domain decomposition methods: application to post-buckling analyses[END_REF][START_REF] Pebrel | A nonlinear dual domain decomposition method: application to structural problems with damage[END_REF]) combined with appropriate convergence acceleration techniques. Second, as Section 3.2 will show, it allows to introduce complex constitutive laws or geometric details in the local model, and perform an exact structural reanalysis without having to modify the global model. This can be helpful for testing various local modifications of a part, since no remeshing is required. Moreover, the local model could very well be analyzed using a separate piece of code which is independent from the global solver, and therefore contain all sorts of unusual refinements such as nonlocal constitutive laws, that would be difficult to implement in the global FE software. This paper focuses mostly on the efficient solving of nonlinear problems; structural reanalysis and solver cou-pling capabilities of the approach will be discussed in future works.

The rest of the paper is structured as follows. In Section 3, we derive the governing equations of the method and investigate some essential theoretical properties. In Section 4, two convergence acceleration techniques that greatly improve the method's efficiency are presented. Finally, Sections 5 and 6 illustrate the behavior and performances of the technique in a commercial FE software environment on a simple 2D test case and a more complex 3D model.

Basic theory

In this paper, we consider an elastic problem with local plasticity. However, this framework is not restrictive; other kinds of local nonlinearity could as well be considered, such as local buckling or contacts occuring inside the domain of interest.

The reference problem

Consider an elastic-plastic statics problem in which all plasticity is supposed to occur within a small region Ω I , called the domain of interest. In the remaining region Ω C = Ω \Ω I , called the complement domain, we assume the constitutive law is purely linear elastic. In addition, all displacements and strains are supposed to be small and no contact is considered; therefore the only source of nonlinearity is the plasticity in Ω I . Also, only one loading increment is considered; multiple load increments are handled by repeating the procedure in an incremental scheme.

This problem is assumed to be well-posed. It can be written as follows (see Figure 2): find a displacement field u ∈ U and a stress field σ ∈ S such that

div σ + f d = 0 in Ω (1) σ = F ε(u), X K : ε(u) in Ω I in Ω C (2)
u = u d on ∂ 1 Ω (3) σ • n = F d on ∂ 2 Ω (4)
where U and S are the respective spaces of regular displacements and stresses; u d , f d and F d respectively denote prescribed displacements, body forces and surface tractions; n is an outerpointing normal; F denotes the (incremental) elastic-plastic constitutive law with internal variables X, and K is Hooke's tensor.

% $ % % & % " % $ ! % & $ " % ' $ " $ " % & $ ! % ' $ !
Fig. 2 The reference elastic-plastic problem

Global-local formulation

As explained previously, we intend to solve this problem using two distinct models: a global model that represents the whole structure but is possibly inaccurate on Ω I , and a local model that represents Ω I alone.

For that purpose, a common approach could be to define the solution as the sum of the global solution and a local enrichment term [START_REF] Dhia | Multiscale mechanical problems: the Arlequin method[END_REF][START_REF] Düster | A combination of an h-and a p-version of the finite element method for elastic-plastic problems[END_REF][START_REF] Mao | A refined global-local finite element analysis method[END_REF]. Since this form may not be very convenient in the presence of nonlinearity, a different choice is made and the local solution is used to replace a part of the global solution, which means the solution is built as:

s(x) = u(x), σ (x) = s L (x) if x ∈ Ω I s G (x) if x ∈ Ω C (5)
where the superscripts G and L respectively denote quantities from the global model (defined on Ω ) or the local model (defined on Ω I alone). By introducing this decomposition into Equations [START_REF] Akgün | Fast exact linear and nonlinear structural reanalysis and the Sherman-Morrison-Woodbury formulas[END_REF][START_REF] Dhia | Multiscale mechanical problems: the Arlequin method[END_REF][START_REF] Cormier | Aggressive submodelling of stress concentrations[END_REF][START_REF] Cresta | Nonlinear localization strategies for domain decomposition methods: application to post-buckling analyses[END_REF], it is clear that the local solution must verify every equation written on inner points of Ω I (inner equilibrium and elastic-plastic constitutive law), as well as boundary conditions written on

∂ 1 Ω I = ∂ 1 Ω ∩ ∂ Ω I (prescribed displace- ments) and ∂ 2 Ω I = ∂ 2 Ω ∩ ∂ Ω I (
prescribed tractions). Defining U I , the space of regular displacement fields on Ω I that match the boundary conditions on ∂ 1 Ω I , and U 0 I , the corresponding linear space (ie. with zero-valued conditions on ∂ 1 Ω I ), this group of equations can be rewritten in the following weak form:

                   u L ∈ U I ∀v * ∈ U 0 I , Ω I Tr σ L ε v * dΩ = Ω I f d • v * dΩ + ∂ 2 Ω I F d • v * dΓ + Γ σ L n i • v * dΓ σ L = F ε(u L ), X (6)
where Γ is the surface between Ω I and Ω C , known as the interface; the subscript i in n i indicates that the normal should point out of Ω I . In a similar way, the global solution must verify the equations written on inner points of Ω C and the corresponding boundary conditions, which can be written as follows:

                   u G ∈ U C ∀v * ∈ U 0 C , Ω C Tr σ G ε v * dΩ = Ω C f d • v * dΩ + ∂ 2 Ω C F d • v * dΓ + Γ σ G c n c • v * dΓ σ G = K : ε(u G ) (7)
where U C and U 0 C are similar to U I and U 0 I , but concern Ω C . The subscript c in n c also indicates that the normal should point out of Ω C , and the notation σ G c is used because our method makes the global stress field become discontinous across Γ , as Section 3.4 will show; since [START_REF] Guidault | A two-scale approach with homogenization for the computation of cracked structures[END_REF] is a subproblem written on Ω C , we need to consider the global stress value on the complement domain's side. Note that the global solution is completely unconstrained inside Ω I .

Finally, the regularity of the constructed solution implies that the global and local displacements and tractions must match on Γ :

u L = u G on Γ (8) σ L n i + σ G c n c = 0 on Γ (9)
Defining U 0 Γ as the space of regular displacement fields on Γ that are equal to zero on Γ ∩ ∂ 1 Ω if such part of the boundary exists, the latter equation can be rewritten as:

∀v * ∈ U 0 Γ , Γ σ L n i + σ G c n c • v * dΓ = 0 (10)
Equations ( 6), ( 7), ( 8) and ( 10) define what we call the global-local formulation.

Notice that U 0 I and U 0 C have the same trace on Γ , which is U 0 Γ . Therefore, by summing the two integrals in ( 6) and ( 7) and choosing a test field that is continuous across Γ and zero-valued on ∂ 1 Ω , the Γ terms vanish thanks to [START_REF] Jara-Almonte | The specified boundary stiffness and force (SBSF) method for finite element subregion analysis[END_REF]; the result is the usual weak form of the reference problem defined by Equations (1-4), which proves that the (continuous) global-local formulation is equivalent to the reference problem.

However, this formulation may induce cheaper computations since ( 6) is nonlinear but concerns a small area, whereas (7) concerns a large domain but contains linear equations only: this framework allows to handle nonlinearity locally, within Ω I , instead of handling it globally. As shown in [START_REF] Cresta | Nonlinear localization strategies for domain decomposition methods: application to post-buckling analyses[END_REF][START_REF] Pebrel | A nonlinear dual domain decomposition method: application to structural problems with damage[END_REF], this approach can lead to a reduced number of global solver calls, and therefore save computational costs.

The two finite element models

To address this formulation, two standard finite element models are used (see Figure 3). The global model represents the whole structure Ω with a linear elastic constitutive law, whereas the local model represents Ω I alone with a nonlinear, elastic-plastic constitutive law. Those models are considered as given; determining the extent of the area that needs to be reanalyzed is not a trivial question (for linear elastic problems, a method is proposed in [START_REF] Srinivasan | Identifying global/local interface boundaries using an objective search method[END_REF]) and would require an appropriate model error estimator (such as those used in the HPDM and some related works [START_REF] Oden | Analysis and adaptive modeling of highly heterogeneous elastic structures[END_REF][START_REF] Zohdi | A method of substructuring large-scale computational micromechanical problems[END_REF][START_REF] Zohdi | A domain decomposition method for bodies with microstructure based upon material regularization[END_REF], again for linear elasticity), but this is beyond the scope of this paper. At present, for the sake of simplicity, the following assumptions are made:

-Ω C is a set of uncut elements of the global mesh (Γ doesn't cut any global element); -Ω I is the exact region covered by the local mesh, and Γ is part of its boundary; both meshes match on Γ : same nodes, edges, faces (in 3D), degrees of freedom and shape or basis functions. Therefore there is no need for a separate discretization of Γ .

It is important to note that there is no "model compatibility" constraint inside Ω I : the two representations of Ω I can be completely different in terms of geometries, discretizations and constitutive laws, given that the meshes match on Γ . This property makes it easier for the user to define the local mesh and allows the method to be used for structural reanalysis and solver coupling purposes, as explained previously; those are natural advantages of the substitution formulation introduced in (5), in particular when compared to additive enrichment approaches [START_REF] Düster | A combination of an h-and a p-version of the finite element method for elastic-plastic problems[END_REF][START_REF] Mao | A refined global-local finite element analysis method[END_REF]. In the future, we intend to make the approach even more flexible by extending it to non-matching meshes, which would allow to define the local mesh in complete independence from the global discretization of Γ .

The global/local formulation is then discretized using the two meshes. Admissibility in Ω I (6) is discretized on the local mesh's FE displacement approximation space, which results in the local FE model ; admissibility in Ω C (7) is discretized on the global FE space, which results in can actually be enforced exactly, and traction equilibrium ( 10) is discretized on the "interface space", ie. the restriction of FE approximation spaces to Γ 's degrees of freedom. Finally, the reference mesh is the mesh that would be obtained by actually removing Ω I 's elements from the global mesh, and replacing them with the local mesh; the reference FE problem is the discretization of Equations 1-4 on this reference mesh, and the reference (FE) solution is the solution to this problem. The three meshes are represented on Figure 3.

! ! ! " " ! # ! " " " #$/"5 #$/"5

The exchange algorithm

To solve the problem in a non-intrusive way, an iterative algorithm is used, which consists in exchanging data between the two models in order to couple the local model with the global representation of Ω C . A simple description of this algorithm is given in Section 2; from a more formal point of view, it can be explained as follows. Starting from a global elastic solution that verifies [START_REF] Guidault | A two-scale approach with homogenization for the computation of cracked structures[END_REF], the following three computations are repeated:

1. A local solution u L that verifies Ω Iadmissibility ( 6) and the displacement transmission condition ( 8) is computed; this is a nonlinear FE problem formulated on the local model. 2. At this point, [START_REF] Jara-Almonte | The specified boundary stiffness and force (SBSF) method for finite element subregion analysis[END_REF] should not be verified. The following linear form is defined:

r(v * ) = - Γ σ L n i + σ G c n c • v * dΓ (11)
By replacing v * with each FE basis function attached to degrees of freedom of Γ (which are the same for the global and local meshes), the residual mentioned in Section 2 is obtained. This quantity measures the lack of equilibrium of the interface in the sense of finite elements. If its magnitude is "small enough", iterations are stopped. 3. Otherwise, global correction is performed in two steps: first, a global corrective term ∆ u G is computed. For that purpose, the global model is loaded on Γ with the residual, and all other prescribed displacements and loads are set to zero, as shown in Figure 4 (the corrective term is therefore "zero-admissible" on Ω C ). The corresponding weak form is:

           ∆ u G ∈ U 0 C ∀v * ∈ U 0 C , Ω C Tr ∆ σ G ε v * dΩ = r(v * ) ∆ σ G = K : ε(∆ u G ) (12)
This is a linear elastic FE problem formulated on the global model; the global mesh and constitutive law are reused without any modification. Then the global solution is updated by simply letting

u G ← u G + ∆ u G (13)
and the process is repeated from Step 1.

This algorithm has similarities with the "Iterative Global-Local" technique proposed by Whitcomb in [START_REF] Whitcomb | Iterative global-local finite element analysis[END_REF], with a slightly different formulation.

Several important remarks can be made; first, the " % 56

$ ! $ " $ % & $ % ' $
Fig. 4 The global corrective problem residual defined in Step 2 is the following FE vector:

r Γ = - Γ (σ G c n C + σ L n I ) • N Γ dΓ ( 14
)
where N Γ is the vector of the FE basis functions of Γ . This quantity has some essential properties:

it defines a FE load vector concentrated on degrees of freedom of Γ ; it can be computed by simply taking the sum of the FE nodal reaction forces exerted on each of the two subdomains across Γ ; it is equal to zero if and only if the FE equilibrium of Γ (10) is verified.

Also, it can be noticed that after the local analysis (ie. at the beginning of Step 2), the globallocal solution defined by [START_REF] Düster | A combination of an h-and a p-version of the finite element method for elastic-plastic problems[END_REF] verifies Equations ( 6), ( 8) and ( 7) (because each of the successive corrective terms defined by ( 12) is "zero-admissible" in Ω C ). The only equation that can be violated is the interface's equilibrium [START_REF] Jara-Almonte | The specified boundary stiffness and force (SBSF) method for finite element subregion analysis[END_REF]. Therefore, the residual is equal to zero if and only if the global/local solution is equal to the reference solution, that verifies every equation of the problem.

By analogy with domain decomposition methods [START_REF] Tallec | Domain-decomposition methods in computational mechanics[END_REF], this technique can be referred to as primal -since global/local displacements are always continuous across Γ , whereas forces are not.

Some convergence properties

From the remarks above, two important properties of the method can be derived. First, a norm of the residual can be used as a simple convergence indicator, for example when compared to the same norm of the global load vector. The validity of this indicator is investigated in Section 5.3.

Second, it is clear that if this algorithm converges, then the limit of the residual must be zero (a nonzero limit would imply ∆ u G = 0 since the global stiffness operator is positive definite); therefore, if the global/local solution converges, then it converges to the reference solution. This ensures the robustness of the method, since the solution error can be reduced as much as needed.

It is more difficult to determine if the algorithm is actually convergent. Some hints can be obtained by writing the global correction equations ( 12) and ( 13) as matrix equations. Denoting by K G the global stiffness operator (which is linear, symmetric and positive-definite) and by u G the vector of global nodal displacements, the global correction can be written as:

u G ← u G + K G -1 r G ( 15
)
where r G is the global FE load vector obtained from the residual r Γ (with the corresponding values on interface degrees of freedom and zeros everywhere else). This equation shows that the algorithm can be interpreted as a modified Newton's method formulated on the equation r G = 0ie. on reaction force equilibrium [START_REF] Jara-Almonte | The specified boundary stiffness and force (SBSF) method for finite element subregion analysis[END_REF]. Therefore, it has the same behavior in terms of convergence. In particular, the algorithm should always be convergent when the material stress/strain curve is monotonic and concave, which is often the case in classical plasticity with positive hardening. This also means that if the local problem is significantly nonlinear, then convergence is likely to be slow, since K G is never updated. In the next section, two convergence accelerations techniques are proposed to overcome this limitation.

Convergence acceleration techniques

The essential idea behind convergence acceleration is that the "basic" algorithm described above is based on a modified Newton's method. As a result, its convergence rate can often be improved by using an updated global operator KG at each global correction, instead of always using the same global stiffness operator K G :

u G ← u G + KG -1 r G (16)
However, it can be difficult to access and modify stiffness matrix terms in commercial FE software; in addition, by doing so, the global matrix would have to be factored again. As a consequence, the solution to [START_REF] Oden | Analysis and adaptive modeling of highly heterogeneous elastic structures[END_REF] should not be computed directly, but rather deduced from the solution to [START_REF] Mao | A refined global-local finite element analysis method[END_REF]. This can be done conveniently with the Sherman-Morrison and Woodbury formulas [START_REF] Akgün | Fast exact linear and nonlinear structural reanalysis and the Sherman-Morrison-Woodbury formulas[END_REF], which are efficient when K G and KG differ only by a small rank corrective matrix. Two particular techniques are presented for that purpose.

A Quasi-Newton global correction

A first possible class of techniques is the Quasi-Newton family, which consists in updating the global operator at each iteration so that the updated operator is "secant" with respect to the previous two solutions and right hand side vectors. Those techniques are particularly suitable for our problem because:

update terms have very small ranks, which keeps computations simple; no additional information transfer from the local model is required; update terms are directly built from solutions and residuals.

Many Quasi-Newton update formulas have been proposed. The SR1 (Symmetric Rank One) formula is used in this study; in our framework, it is easier to implement than rank two methods such as BFGS or DFP, and gives better results than nonsymmetric methods such as the "good" or "bad" Broyden formulas, that do not preserve the global stiffness operator's symmetry. The classical SR1 update formula is:

K k = K k-1 - (∆ r k + K k-1 ∆ u k )(∆ r k + K k-1 ∆ u k ) T (∆ r k + K k-1 ∆ u k ) T ∆ u k ( 17
)
where subscript k indicates the iteration number

(k ≥ 1), ∆ u k = u k -u k-1 and ∆ r k = r k -r k-1 . Since K k-1 ∆ u k = r k-1
(no relaxation is used), this equation can be rewritten as:

K k = K k-1 - r k r k T r k T ∆ u k (18)
which shows that only the interface terms need to be updated, since the residual's only nonzero values are on Γ . Therefore, the formula can be applied directly to the global stiffness matrix, using K 0 = K G (if Ω I 's inner stiffness were to be modified, the formula would have been applied to the local model, then transferred to the global model using some kind of homogenization).

In order to solve ( 16), the corresponding inverse update equation is formed by applying the Sherman-Morrison formula, and right-multiplied by the residual r n to avoid manipulating matrices. The result is:

∀n ≥ 1, ∀k ≤ n, K -1 k r n = K -1 k-1 r n + K -1 k-1 r k r T k K -1 k-1 r n r T k ∆ u k -K -1 k-1 r k (19)
This formula is used as follows. At the first global correction (n = 1), the residual r 1 is first sent to the global solver like before (K -1 0 r 1 ), and then the updated solution is obtained by applying the formula once, with k = 1. More generally, at iteration n, first the global solver is called once, and then the formula is applied n successive times, with k ∈ [|1, n|]. This requires storing up to three vectors per iteration: the successive residuals r k and update terms ∆ u k = K -1 k-1 r k-1 , and also the K -1 k-1 r k terms.

The main advantage of this technique is that the global solver is still called once per iteration: no additional call is required, because the updates are deduced from previous results. Consequently, extra computational costs are negligible whereas convergence can be greatly accelerated. Moreover, the procedure only manipulates vectors defined on interface degrees of freedom, and is very easy to implement.

Tangent global correction

In order to obtain faster convergence than the Quasi-Newton approach, another possibility is to assess the tangent behavior of the structure, resulting in a global/local version of Newton's classical method. However, this raises several difficulties: first, displacements and residuals alone are no longer sufficient, and information about the local model's tangent stiffness must be obtained, which is not always easy. Second, since the meshes can be completely different inside Ω I , a homogenization technique must be used to transfer this information to the global model. Finally, by doing so, the corrections that must be applied to the global matrix are unlikely to have a small rank, and ultimately the additional costs arising when applying Woodbury's formula might exceed the benefit of accelerating convergence. Therefore, care should be taken to keep data transfers simple enough, to ensure the technique is computationally efficient.

To overcome these difficulties, we can take advantage of the way the models are coupled. All exchanges concern quantities defined on Γ alone, which means surface coupling is used. As a consequence, it can be proved that convergence rates only depend on the Schur complement of the updated global correction operator on Γ , and not of its detailed values inside Ω I . More specifically, if the local problem is assumed to be linear, then performing a global correction with an operator having the same Schur complement as the reference problem's operator will lead to convergence in just one iteration -which is one possible definition of a tangent Newton method.

This can be proved by writing the FE matrix equations in condensed form and retaining only interface degrees of freedom. The global FE problem is written as:

S G u G Γ = S G c + S G i u G Γ = b G c + b G i = b G (20)
where subscripts c and i respectively denote the contributions of Ω C and Ω I to the Schur complements S G and to the condensed right hand side vectors b G . Similarly, assuming the local problem is linear (but different from the global problem !) its condensed equation is:

S L u L Γ = b L + g L C→I (21)
where g L C→I denotes the reaction forces exerted by Ω C on the local model through the interface. Now let u G Γ be an initial global solution, with an arbitrary value: since the local analysis is performed with the prescribed displacement condition u L Γ = u G Γ , Equation ( 21) gives:

g L C→I = S L u G Γ -b L (22)
The reaction forces exerted by Ω I on Ω C within the global model are defined in a similar way:

g G I→C = S G c u G Γ -b G c ( 23
)
Taking the sum of those reaction forces, which should be balanced, gives the residual:

r Γ = -g L C→I + g G I→C = b L + b G c -S L + S G c u G Γ = b R -S R u G Γ (24)
In this expression, the Schur complement and condensed right hand side vector of the reference problem (obtained by substituting the local problem into the global problem) can be recognized. This proves two things:

-The limit of u G Γ , if it exists, is the reference solution u R Γ as shown in Section 3.5.

-Performing an updated global correction with

any global operator whose Schur complement equals S R = S L + S G c will make the algorithm converge in one iteration in this linear case:

u G Γ ← u G Γ + S R -1 r Γ = u G Γ + S R -1 b R -S R u G Γ = u G Γ + S R -1 b R -u G Γ = S R -1 b R = u R Γ ( 25
)
This is true for any initial estimate u G Γ , no matter how different are the two representations of Ω I (ie. S L and S G i ). As a result, in the general, nonlinear case, a global tangent operator can be constructed at the end of the local analysis, by performing static condensation on the local tangent problem (after removing all external loads) and then updating only the interface terms of K G . This is done by first computing the difference ∆ S = S L -S G i , and then applying a particular form of Woodbury's matrix identity [START_REF] Akgün | Fast exact linear and nonlinear structural reanalysis and the Sherman-Morrison-Woodbury formulas[END_REF] to perform a non-intrusive "tangent global correction":

K G + R T ∆ SR -1 = K G -1 -K G -1 R T ∆ S I + RK G -1 R T ∆ S -1 RK G -1 ( 26
)
where I is the identity matrix, R is the restriction operator that keeps the interface values of a vector u G and conversely, its transpose R T extends an interface vector to the whole global model by appending zeros. The implementation is more complex than the Quasi-Newton technique because additional global solver calls can no longer be avoided. By right-multiplying Equation ( 26) with a residual r G = R T r Γ , it appears that this technique requires to compute K G -1 R T first -which means, to prescribe unit concentrated forces on each interface degree of freedom and to collect the corresponding global displacements. This operation needs to be performed only once, and the whole set of load vectors can be processed simultaneously; then no further global solver calls are needed, as the arising residuals can be expressed on the basis of elementary global responses that have been computed previously.

Using this method, a global tangent operator can be obtained, which leads to faster convergence than the Quasi-Newton method. However, memory requirements can be very high on large models because many global load vectors need to be processed at the same time.

A simple 2D example

The different variants of the analysis strategy described above were tested on a simple 2D model with Abaqus/Standard, version 6.7-1, which was used for both global and local analyses. The two meshes and subdomains are shown on Figure 5; the structure was clamped on the bottom and lower side edges, and the rest of the boundary was traction-free. Pressure and centrifugal loads were applied monotonically in one increment. The global model was linear elastic; in the local model, a refined mesh and an elastic-plastic constitutive law (with linear isotropic hardening) were introduced, while the geometry, material density, loading... were similar to those of the global representation of Ω I . Parameters were chosen so that the elastic limit would be significantly exceeded at a couple of integration points inside Ω I while the inelastic zone would remain small.

Both models were designed as a typical submodeling data set; performing the strategy on them only required minimal modification of the input files, and the algorithm itself was implemented as a Python script, writing prescribed conditions on Γ to text files, submitting jobs and reading displacements and reaction forces from output databases. As a consequence of the non-intrusive framework, no other software or piece of code was used except for convergence acceleration, which required the ability to handle vectors and perform dot products -and also, for the tangent update, to handle full symmetric matrices and solve small linear systems. Those simple operations were performed within the script itself; everything else was actually performed through the Abaqus scripting interface.

Initial submodeling errors

At the beginning of the method, a global linear elastic analysis is performed. Then global interface displacements are extracted, and a local elasticplastic analysis is carried out with prescribed displacements on Γ . At this point, the global/local solution (obtained by substitution) is similar to what would be obtained with classical submodeling.

Figure 6 displays the global/local Von Mises equivalent stress around the interface. As one could predict, a discrepancy in colors can be observed between both sides of Γ , which means stresses are discontinuous across the interface and the solution is inaccurate. To measure this inaccuracy, the reference solution was computed, then Fig. 5 Global and local meshes of the 2D test case Fig. 6 The Von Mises stress discontinuity obtained with submodeling displacement and stress errors with respect to the reference were estimated as:

η u = ||u GL -u R || ||u R || (27)
η σ = ||σ GL -σ R || ||σ R || (28)
where the norm || || is the Euclidean norm of the vector collecting all available quantities (ie. displacements at all nodes and stresses at all integration points). An error of about 10% for displacements and 5% for stresses was obtained, which illustrates the need to enforce a better coupling between the two subdomains. Note that for other nonlinear phenomena such as local buckling, such an approach can lead to much greater errors (see [START_REF] Cresta | Nonlinear localization strategies for domain decomposition methods: application to post-buckling analyses[END_REF]).

Convergence rates

Now we investigate how the errors decrease along the iterations. Figure 7 represents the interface displacement error plotted against the iteration number; this is defined as

η u Γ = ||u GL Γ -u R Γ || ||u R Γ || (29)
ie. the displacement error restricted to the interface only; this requires a much lighter amount of output than the full solution error, and is meaningful because the exchanges are driven with interface displacements. Note that in the general case, this error is unknown. The three curves correspond to the three variants of the method: non-accelerated, Quasi-Newton (secant) and tangent.

As expected, the basic variant of the strategy seems to converge to the reference solution, but the process is quite slow: according to Figure 7, it

- . & / 0 1 2 , 3 4 .- .561
.560

.56/ .56&

.56. 789:)%+;+<9#"+"9989 =7+;+<9#"+"9989 <)>*">?+;+<9#"+"9989 @?"9)?$8>+>#:A"9 @>?"9B)C"+D$EF%)C":">?+"9989 Fig. 7 Compared convergence rates of the three variants, for the 2D case takes about 7 iterations before obtaining a solution error less than 0.1% whereas the nonlinear reference computation, which has a similar size, ran for 4 iterations before a similar accuracy was attained. The Quasi-Newton acceleration offers significant improvements over the non-accelerated method, since convergence is much faster whereas the computational cost of one iteration is almost the same: only 3 iterations are needed to bring the solution error down to 0.1% (the initial global and local analyses are counted as iteration zero). In comparison, the full tangent update performs even better at the beginning, and the 0.1% threshold is almost reached after 2 iterations; however, the convergence rate then deteriorates and the Quasi-Newton method becomes more accurate. A possible explanation for this phenomenon could be that the nonintrusive update technique, based on Woodbury's matrix identity, can sometimes lead to numerical instability [START_REF] Yip | A note on the stability of solving a rank-p modification of a linear system by the Sherman-Morrison-Woodbury formula[END_REF]. Since single precision storage had to be used for technical reasons, the quality of the "tangent" operator might not be sufficient to improve convergence rates any further. In addition, this technique can be costly since it requires to analyze many global load vectors at the same time (as many as interface degrees of freedom), which could easily exceed memory limits of today's computers on large problems.

It can be noticed that displacement errors do not actually converge to zero. Their limit was found to depend on the convergence threshold of the elastic-plastic constitutive law, and can be reduced by setting this threshold lower (at the expense of local computational costs). It should also be noted that the Quasi-Newton method actually used less global solver calls than the reference analysis to bring the solution error under 0.1% (3 iterations versus 4), which confirms that it can be more efficient to handle local nonlinearity locally that globally [START_REF] Cresta | Nonlinear localization strategies for domain decomposition methods: application to post-buckling analyses[END_REF][START_REF] Pebrel | A nonlinear dual domain decomposition method: application to structural problems with damage[END_REF].

Convergence indicators

The interface displacement error discussed above provides a good measurement of the average error; however, there is no way to compute it when the reference interface displacement is unknown. In this section, two possible convergence indicators are discussed.

The first quantity that can be used to assess convergence is the residual's relative norm, or relative residual. It is computed after each local analysis and is defined as /56/ 789:*%,;,<"=$>#*% 789:*%,;,?@*+A*@$8A 789:*%,;,B9#","9989 C7,;,<"=$>#*% C7,;,?@*+A*@$8A C7,;,B9#","9989 D@"9*@$8A,A#:E"9

η r = ||r Γ || ||b Γ || (30)
<"%*@$F","9989

Fig. 8 Evolution of the "true error" and convergence indicators, with and without acceleration, for the 2D case where ||b Γ || is an estimation of the condensed right hand side vector's norm (here, the norm of global reaction forces was used as a rough approximation). In Section 3.5, we proposed to use this quantity as a convergence indicator, since it is equal to zero if and only if the exact reference solution is reached.

On Figure 8, this indicator is compared with the "true error" (ie. the interface displacement error) for the non-accelerated and the Quasi-Newton variants of the method. The dotted curves show the evolution of the interface displacement errors. The thick plain curves represent the relative residuals. We can see that for the non-accelerated method, both quantities are almost equal: they have similar magnitudes and decrease at the same constant rate. For the Quasi-Newton method, the two decrease rates are slightly different, and a factor of almost twenty is found between the two quantities.

Judging from this curve, the relative residual seems to provide an upper bound for the solution error; however, the general rule is not that straightforward. Let us consider the linear case discussed previously; from Equation (24), it can be noticed that the residual is just the difference between the two sides of the condensed reference problem (ie. it is actually the "residual", in the sense of iterative methods). Therefore, the algorithm behaves like a basic iterative method trying to cancel the residual by passing it through a preconditioner, which can either be S G (for the non-accelerated version) or the "tangent" or "secant" Schur complements (for the accelerated versions). As a consequence:

-The relative residual can provide bounds for the solution error, but their quality depend on the condition number of the reference problem's total Schur complement S R :

1 κ(S R ) η r ≤ η u Γ ≤ κ(S R ) η r (31)
Therefore both quantities should be of similar orders of magnitude when S R is wellconditioned, which is the case here. Usually, Schur complements are better conditioned than full stiffness matrices, but in practice, these bounds may not always be sufficient to provide an accurate error measurement, as shown above.

-The successive residuals and displacement errors decrease such that

(u G k+1 -u R ) = G•(u G k - u R ) and, since the operators are symmetric, r G k+1 = G T • r G k
, where:

G = I -S G -1 S R (32)
for the non-accelerated variant (G is called the iteration matrix). As a result, the decrease rate of both quantities is bounded by the spectral radius of G, which directly depends on the relative "change in stiffness" introduced in the local model.

This shows that the residual's norm can provide information about how fast the solution error decreases, but it does not directly measure the solution error, and in particular it does not provide bounds for stresses or displacements, as shown in [START_REF] Whitcomb | Iterative global-local finite element analysis[END_REF]. To obtain more information, we introduce a second quantity called the stagnation indicator; it is computed after each global correction and is defined as

η ∆ u = ||∆ u Γ || ||u Γ || (33)
where ∆ u Γ is the current global corrective term, computed with any of the three variants of the method; this quantity converges to zero when interface displacements stabilize, hence its name. Its evolution is represented by the thin plain curves on Figure 8. It can be seen that for the nonaccelerated variant, this indicator decreases at the same rate as the "true error", but has lower values; this is a typical behavior of the modified Newton method in such cases (corrective terms tend to be "too small" because the global stiffness operator S G is stiffer than the reference operator S R ). For the Quasi-Newton variant, the stagnation indicator is almost equal to the "true error" during the first iterations (then the error stops decreasing, since the constitutive law's inaccuracy is no longer negligible). This is a consequence of the SR1 update technique's efficiency: previous corrective terms are reused in an optimal way, and the modified ∆ u Γ is "as close as possible" to the error.

Therefore, this indicator could be used in conjunction with the relative residual (and, if available, the constitutive law's local threshold) to form a simple stop criterion, like in classical Newton methods. More advanced estimators could also be used, such as the model error estimators used in the HPDM and some related works [START_REF] Oden | Analysis and adaptive modeling of highly heterogeneous elastic structures[END_REF][START_REF] Zohdi | A method of substructuring large-scale computational micromechanical problems[END_REF][START_REF] Zohdi | A domain decomposition method for bodies with microstructure based upon material regularization[END_REF].

3D results

In order to assess the strategy's performance and ease of use in industrial situations, addi-Fig. 9 Global and local meshes of the 3D test case tional tests were ran on a larger 3D test case provided by Snecma. This case consists in two Abaqus/Standard models represented on Figure 9 and is based on an aircraft engine's turbine blade, with an extremely simplified geometry for confidentiality purposes. The global model is subjected to a complex set of loads comprising centrifugal forces, a pressure field on the blade's surface, and a precomputed temperature field; normal displacements are blocked on the slanted surfaces at the bottom, and the constitutive law is linear elastic as usual. The local model is located at the bottom of the blade and contains the plastic area, with reasonable margins; it features an elastic-plastic constitutive law and its mesh is locally identical to the global mesh (no refinement was introduced).

The two models had been prepared as a typical submodeling data set and as explained previously, few modifications of the input files were needed to run the strategy. The global problem totalled about 500,000 DOFs, with each global solver call tak-

, - & . / 0 + 1 2 3 -, -451 -45+
-450

-45/ -45.

-45&

-45-6789 (%*:*;"<$=#(% 6789 (%*:*>?()@ (?$7@ 6789 (%*:*A8#"*"8878 B6*:*;"<$=# (% B6*:*>?()@(?$7@ B6*:*A8# "*"8878 C?"8(?$7@ *@#9 D"8

;"%(?$E"*"8878< Fig. 10 Evolution of the "true error" and convergence indicators, with and without acceleration, for the 3D case ing about 3 minutes for one single load; the local problem contained about 80,000 DOFs and the interface contained 6,400 DOFs. Figure 10 represents the evolution of the "true" (interface displacement) error and the two convergence indicators discussed previously, for the non-accelerated and the Quasi-Newton variants. It can be noticed that the initial solution error (ie. the error at "iteration zero", which corresponds to one-way submodeling) has a value of about 1%, which is already quite low; this means stress redistributions are not very critical on this test case and submodeling is sufficient to attain that level of accuracy. However, if more accurate results are needed, several iterations may be needed. Convergence rates are similar to the 2D results; it only takes 2 iterations of the Quasi-Newton variant to reach a 0,1% solution error, versus 7 iterations of the non-accelerated variant (note that here, the solution error does apparently converge to zero because the constitutive law's threshold was set to a lower value). Convergence indicators also show similar behaviors: for the non-accelerated method, they all decrease at the same constant rate and the stagnation indicator underestimates the "true error". For the Quasi-Newton method, the relative residual overestimates the solution error by a factor up to ten, and once again the stagnation indicator gives an accurate estimate of the "true error". Note that for this test case, the tangent vari-ant couldn't even be tested, because it would have required to analyze 6,400 global load vectors simultaneously (each with 500,000 DOFs), which wasn't affordable in reasonable times.

Conclusions

In this paper, an iterative technique was presented, which allows to analyze linear elastic structures containing small plastic areas and to handle nonlinearity in a separate local model while converging to the whole elastic-plastic solution. Its formulation is non-intrusive, which means it can be used in a standard FE analysis environment. The simplest variant of this technique is based on a modified Newton's method, and basically has the same convergence properties; since it can be quite slow, two accelerations techniques were proposed, respectively based on a Quasi-Newton and a full tangent Newton algorithms. The Quasi-Newton variant is simple to implement, and significantly improves the convergence rate. The full tangent variant seems even more promising, but requires to process many local and global load vectors at once, and therefore may be too expensive to be used on large models in its current form.

To improve those results, two possibilities are foreseen. The first one is to modify the full tangent variant to avoid performing a complete con-densation of both models on their interface, but rather a reduced condensation on a small basis that includes rigid body motions, similarly to the micro/macro LaTIn approach [START_REF] Ladevèze | A micro/macro approach for parallel computing of heterogeneous structures[END_REF][START_REF] Guidault | A two-scale approach with homogenization for the computation of cracked structures[END_REF]. This should ensure a proper transmission of the resultant reaction forces and moments between the nonlinear and linear subdomains, which should be enough to ensure an efficient coupling, according to Saint-Venant's principle. Hence, convergence could be accelerated without having to analyze many load cases at once, which should result in more acceptable memory requirements.

The second possibility is to replace the prescribed displacement condition, which is currently used to drive the local analysis, with a mixed (Robin) condition that uses both displacements and stresses. This should improve the local analysis' accuracy from the first iteration; 1D experiments proved that such conditions can lead to extremely fast convergence. Such a technique once again requires to estimate the Schur complements of each subdomain; combining it with the reduced condensation discussed above could lead to very high efficiency.

Finally, two possible extensions of the method are considered. Currently, the strategy processes each load increment one by one; following another idea of the LaTIn method [START_REF] Ladevèze | A micro/macro approach for parallel computing of heterogeneous structures[END_REF], it could be more convenient and efficient to handle the whole loading history at once, and therefore to solve whole evolution problems at each iteration. Also, extending the exchange technique to non-matching meshes would probably involve technical difficulties, but would certainly allow this technique to address a much more general class of locally nonlinear problems, and to be a more flexible structural reanalysis and model coupling tool.

Fig. 3

 3 Fig. 3 Global, local and reference FE models (darker elements are elastic-plastic)

  the nonlinear area u global everywhere else Finally, an iterative exchange technique is used to couple the local model with the accurate area of the global model, as represented on Figure 1.

	"
	! )	! (	"
	*+,-.+'/,0 0'-123.,4,56-	789'2'/+'94 +,-'0932	! '	! (
	"
			Starting from a global elastic solution, each itera-
			tion goes as follows:
			1. Local analysis. A full nonlinear analysis is
			carried out on the local model, with current
			global displacements prescribed as a boundary
			condition.
			2. Residual computation. The residual, a load
			vector measuring the unbalance in nodal re-
			action forces between the linear and nonlin-
			ear areas, is formed. If its magnitude is small
			enough, iterations are stopped.
			3. Global correction. Otherwise, the residual is
			applied on the global model as an additional
			inner load, which causes it to deform as if local
			details were present. The global model's lin-
			ear elastic constitutive law and boundary con-
			ditions are retained. Optional convergence ac-
			celeration techniques can be used at this step,
			then the global solution is updated and the pro-
			cess is repeated.

This work is supported by Snecma and is part of the MAIA-MM1 research program.