Dimitri Serrano-Velarde

Eutelsat

Hector Fenech
email: hfenech@eutelsat.com

Georges Rodriguez-Guisantes

A NOVEL DIMENSIONING METHOD FOR HIGH THROUGHPUT SATELLITE DESIGN

This work describes a novel methodology for the dimensioning of a Ka-Band high throughput satellite for broadband communications. The method is based on the optimization of performance for a forward link, as a function of a set of input criteria and a given envelope of available power. This approach is based on a spacecraft architecture using a multi-beam coverage implementing frequency re-use. Among the input criteria, we use the percentage of covered service area with a certain type of earth stations and the service availability. The proposed methodology is adaptable to any kind of service area. A DVB-S2 air interface with an Adaptive Coding and Modulation (ACM) is used as a reference. The method, employing iterative advanced link budget calculations including carrier-to-interference at antenna level, provides the highest capacity given a batch of antenna and pragmatic feed design.

Introduction

High Throughput Satellites (HTS) orientate the satellite systems towards multi-beam coverages at high frequency bands in order to offer broadband access on a large area while maximizing the usage of frequency reuse. Indeed, the terrestrial broadband offer does not cover large territories and a satellite broadband alternative is considered as a complement rather than a direct competitor to fiber optics [START_REF] Beretta | Keynote Speech for the 16th Ka Band Confernce[END_REF]. There are numerous ways of designing satellite architectures and different kinds of trade-offs based on different parameters. Those parameters range from the ground segment with the level of complexity of a user terminal and a Gateway to the space segment with the whole definition of the spacecraft system and antennas. In order to provide a clarification of the choices, the logic presented hereafter aims at showing a novel methodology which helps achieving the best trade-off between a set of input criteria and opens an efficient dialog between system engineers and antenna engineers. This will provide a new level of cooperation integrating these sometimes opposite engineering approaches and will accelerate the designing phase of the satellite. Also, the resources can be quickly focused on a selected set of technical solutions that will be analyzed in detail. In the first part, we are going to show the issues to which the methodology is going to provide answers, then the steps for the reasoning, and finally present some technical solutions suggested for different scenarios.

Problem Description

When a satellite operator orders a new satellite, it invests effort to assess the different transmission scenarios and services by taking into account all the criteria. For instance, it is challenging to align the system requirements with the antenna requirements as one impact the other and vice versa. So it is important to achieve a hierarchy in those requirements to reach the best trade-off. A general approach has been defined by a satellite manufacturer, but in the frame of beam hopping and for a very industrial approach [START_REF] Anzalchi | Beam Hopping in Multi-Beam Broadband Satelllite Systems: System Simulation and Performance Comparison with Non-Hopped Systems, Advanced satellite multimedia systems conference (asma)[END_REF]. Also another approach has been developed in [START_REF] Blanco | Design Methodology For High Capacity SATCOM Systems[END_REF] but takes into account very specific traffic constraints. Our methodology is based on the case where the satellite uses a multi-beam coverage with frequency re-use and spatial separation. Starting with a set of high level inputs (available DC power, antenna reflector size, etc.), the engineer has to choose first a number of spots and the frequency plan i.e. frequency/polarization organization over a given Service Area (SA). It then achieves the best architectural trade-off through an iterative approach by increasing or decreasing the number of spots for instance. The chosen initial goals are the service area filling percentage (CP) and the availability (AV) of the communication system. Indeed, these criteria are very important for a satellite operator as the main aim is to guarantee the most robust communication system as long as possible and to be able to cover the largest area so as to reach bigger market opportunities or to strengthen already existent markets. Besides, this set of parameters, can be considered as "macro parameters" and can be provided by the commercial and marketing needs before the design and they offer the possibility to test the system in accordance with the objectives of a satellite operator. Finally, the optimization criterion is the achieved maximum capacity. In order to ease the procedure without loss of generality, the link budget calculations are focused on the Forward Link (from Gateway towards the end users) and more specifically on the Downlink (Satellite to terminals) as this link segment is the most constraining and systems dimensioning for the capacity.

The key steps of the method itself are the following: first of all it is necessary to define the input parameters (frequency plan, DC power, terminal, etc.) and a service area covered by a certain number of spots. Then it is necessary to generate the EIRP or T G radiation diagrams so that by using these diagrams, it is possible to calculate first the I C performances at antenna level and then an advanced link budget in order to obtain the raw performances in terms of

I N C  .
Then through an iterative approach, it is possible by increasing or decreasing the number of spots to find out which antenna configuration satisfies the AV and CP criteria in terms of total capacity. Ultimately, the chosen system will be the one which provides the highest capacity while maintaining the initial goals. This algorithm has been applied to several communication scenarios over different areas and a system has been achieved each time. Also by changing the AV and CP criteria, the method is flexible enough to provide a new solution. The link budget model has been build upon data provided mostly by two references [START_REF] Maral | [END_REF] and [START_REF] Morgan | Communications Satellite Handbooks[END_REF] so as to be adapted to specific needs later on. Finally, an improved model including mass and cost is also developed as an extension of this method.

Proposed Models

To begin with, the model is based on several "modules" and input data which are used for all the calculations.

 Antenna Model

In order to start the reasoning, it is necessary to generate the antenna system with software tools and convert the radiation diagrams into usable formats for the link budget, in other words EIRP density or T G matrices. The adopted reflector model is designed by taking into account the illumination taper and spillover losses but also by neglecting the scan aberration. It is based on the model proposed by Peter Balling [START_REF] Balling | Spacecraft Multi-Beam and Contoured-Beam Antennas[END_REF]. First, the aperture distribution from a single feed is calculated by:

n a r p p r g ) ) ( 1 ( * ) 1 ( ) ( 2     with:
 a, the radius of the aperture  p, the relative edge illumination  n is typically 1, but may be changed in the illumination exponent field

The element beams are approximated by a linear combination of Bessel functions:

)) , 1 ( * ) * , 1 ( * ( * * ) , ( 2 1 j j j kax n BF c ax k BF c a k F      with:  k, the propagation constant   2  j x , the distance from beam center, ) ) ( ) (( 2 2 j j j v v u u x      BF(n,x
), Bessel functions given by:

       1 ) 0 , ( ) ( !* * 2 ) , ( n BF x x J n x n BF n n n
 The coefficients 1 c and 2 c depend on the edge taper p, and are normalized so that 2 j F yields directivity. Balling uses an analytic approximation to determine directivity. We prefer perform a PO integration of an on-focus feed to determine the directivity, which is more accurate.

 Interference Model

Once these radiation files are created, it is necessary to define the cells for this multi-beam architecture. The cell definition is also necessary, as each user will be linked to a certain spot and the satellite has to be able to organize the traffic on every spot. Based on this cell division, it is possible to calculate the I C ratio as a performance indicator at antenna level. Indeed, for a multi-beam antenna system implementing frequency reuse and spatial separation, it is important to calculate the interferences generated by all the beams on each other, as the more spots one includes, the more interference is generated. In a general way, it is possible to define the calculation of the I C as follows.

Let's consider a certain spot k with a directivity of k C . Considering, that the architecture is a multibeam coverage with frequency reuse, there will be two types of interferers. Depending on the side of the transmission, there will be the interferers in co-polarization, which are at the same frequency and the same polarization but used in different spatial spots. Their contribution can be quantified as:

   N q Co q co x D I 0 )) ( ( with:
 q, the identifier for an interferring spot  N, the total number of interferers in Co-polarization  x, the point defined by a coordinate system

 Co q D ,

the directivity in Co-polarization for the interferer q

There are also the interferers related to the cross polarization, which are the spots at the same frequency but with a directly opposite polarization. Their contribution is as follows:  Link Budget and Design Model If the performances are acceptable for the coverage area, this data is injected into an advanced link budget model. This model takes into account all the different station parameters, link parameters and satellite parameters. Also, the different attenuations, both technical (such as the interferences) and natural (such as rain, clouds, ...), are taken into account for the final link calculation. Once the data has been calculated, it is possible to apply the logic itself. To begin with, one of the results of the link budget is the performance of the I N C  under clear sky conditions i.e. no rain attenuations, but with all the other attenuations. Corresponding to the proper service area, a performance criterion is applied, in order to identify the users able to receive the lowest type of modulation and able to maintain a satellite communication. By applying a set of rain attenuations based on the models recommended by the ITU it is possible to extrapolate via an iterative procedure for the link budget, the availabilities that can be achieved. Combining this data with the percentage of coverage allows obtaining a solid view of what antenna structure is the most efficient and interesting to go into detail. In order to underline this approach, the next section will focus on two given antenna scenarios and show how the method can be applied.

Results

There will be two scenarios that are going to be discussed in the following part. One will be based on four 2.8 meter antenna reflectors with a service area centered over France and a second one based on four 2 meter antenna reflectors with a service area centered over Central Africa. Both scenarios will use an orbital position of 0°E as an arbitrary position.

 Scenario 1

This scenario aims at providing a service over an elliptical polygon covering mainly France. The satellite is placed at an orbital position of 0°E. The antenna system is composed of four identical reflectors of 2.8m (state of the art) and the number of included spots is allowed to vary over the area. The frequencies used are 19.7GHz and 20.2GHz in Ka-Band with four colors. The air interface is the DVB-S2 standard. With Figure 2, it is possible to see the performance of each antenna system and highlight a clear trend: the more spots you include for a given service area, the more your power is divided per spot and the more interference is generated. Consequently, lesser coverage percentage is achieved for every availability. By applying the AV criteria of 99.96% and the CP criteria of 98%, the method highlights which antenna systems (with a fixed antenna aperture of 2.8m) is the most optimized in terms of capacity as shown in Figure 3. Based on this set of figures, it appears that the antenna system providing the highest capacity while achieving the criteria is the system with a beam spacing of 0.19°(72spots) which corresponds to an approximate spot size of 0.21°. Also, a trend is highlighted for the service area: the more spots are included the more capacity is available but beyond a given number of spots, the criteria are no more respected. So, if a satellite engineer wants to develop a Ka-Band multi-beam Satellite on this kind of region, more fine tuned analysis could be focused around a spot architecture of seventy-two spots corresponding to a beam spacing of 0.19°.

 Scenario 2

This scenario is based on an elliptical polygon covering Central Africa. The satellite is placed at an orbital position of 0°E. The antenna system is composed of four identical reflectors of 2m and the number of stacked spots varies over the area. The downlink frequencies used are 19.7GHz and 20.2GHz in Ka-Band. The air-interface is the DVB-S2 standard. No Eb threshold of -1dB and an operational margin of 3dB, it is possible to obtain Figure 5 showing the behavior of the different antenna systems with a varying number of spots in terms of coverage percentage versus given availabilities. With Figure 5 it is possible to see the performance of each antenna system and also the same trend as in Scenario 1. Compared to the previous scenario, the curves have more losses due to the high attenuation region close to the equator causing a clear "break" at an availability of around 99%. By applying the AV criteria of 98% and the CP criteria of 97%, the method will highlight which antenna system with what kind of beam spacing is the most optimized in terms of capacity as shown in Figure 6. Based on this set of figures, an antenna system with a beam spacing of 0.5°(137 spots) appears to offer the best trade off in terms of capacity. On a side note, due to the high attenuations and inter-spot interferences, it appears clearly that the criteria have to be reassessed compared to the first scenario.

Further Work

Further work will be spend on using real antenna data generated with more accurate models taking into account the complete antenna system including real clusters and all kind of losses. These calculations shall provide a solid view on the behavior of the methodology. Also a new module will be added to this model taking into account the mass and cost of the different antenna structures and add a new layer of design.

Conclusions

It has been demonstrated that the methodology gives very consistent and reliable results for the different scenarios. Indeed, by using this methodology it is possible to focus quickly and efficiently the design optimizations around only a few cases and find the best trade-off. As shown in the two scenarios, no matter which area is studied, a possible trade-off can be found according to the requirements set by the system designer. This shows also how adaptive the methodology is and that the reasoning itself is system defining and not area dependent.



  the identifier for an interferring spot  M, the total number of interferers in the orthogonal polarization  x, the point defined by a coordinate system  Cx p D , the directivity in the orthogonal polarization of interfering spot p Combined all together, the It is also important to include the Beam Pointing Error (BPE), which can influence heavily on the performances.

Figure 1

 1 Figure 1 shows an example of EIRP coverage and the given

Figure 1 :

 1 Figure 1: EIRP and C/N+I performance maps over France

Figure 2 :

 2 Figure 2: Coverage Percentage vs Availability, Scenario 1

Figure 3 :

 3 Figure 3: Capacity vs Number of Spots, Scenario 1

Figure 4

 4 Figure 4 shows an example of EIRP coverage and the given

Figure 4 :

 4 Figure 4: EIRP and C/N+I performance maps over Central Africa

Figure 5 :

 5 Figure 5: Coverage Percentage vs Availability, Scenario 2

Figure 6 :

 6 Figure 6: Capacity vs Number of Spots, Scenario 2


P Kerfriden

J C Passieux

S P A Bordas

Local/global model order reduction strategy for the simulation of quasi-brittle fracture

Keywords: Adaptive Model Order Reduction (MOR), Local/global Approach, Nonlinear Fracture Mechanics, Proper Orthogonal Decomposition (POD), Newton/Krylov Solver

This paper proposes a novel technique to reduce the computational burden associated with the simulation of localised failure. The proposed methodology affords the simulation of damage initiation and propagation whilst concentrating the computational effort where it is most needed, i.e. in the localisation zones. To do so, a local/global technique is devised where the global (slave) problem (far from the zones undergoing severe damage and cracking) is solved for in a reduced space computed by the classical Proper Orthogonal Decomposition, while the local (master) degrees of freedom (associated with the part of the structure where most of the damage is taking place) are fully resolved. Both domains are coupled through a local/global technique. This method circumvents the difficulties associated with model order reduction for the simulation of highly non-linear mechanical failure and offers an alternative or complementary approach to the development of multiscale fracture simulators.

Introduction

Simulating damage initiation and subsequent global structural failure is one of the most active topics in computational mechanics. Several mathematical models and numerical methods have been developed over the years to assess various limit states such as failure due to permanent deformations, cracks or decohesion/delamination, e.g. in composite materials. Yet, these models, be they damage based or relying on discrete cracks are generally computationally expensive, as they require a fine scale description of the structural and material properties. Therefore, today's engineers are not able to use these state-of-the-art models for routine design. For important recent advances in the treatment of material failure (e.g. discontinuous fracture [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF], advanced damage models [START_REF] Allix | Modeling and simulation of crack propagation in mixed-modes interlaminar fracture specimens[END_REF][START_REF] Jefferson | A model for cementitious composite materials based on micro-mechanical solutions and damage-contact theory[END_REF], damage plasticity models [START_REF] Ladevèze | An enhanced mesomodel for laminates based on micromechanics[END_REF][START_REF] Grassl | Damage-plastic model for concrete failure[END_REF] or their combination [START_REF] Rabczuk | A geometrically non-linear threedimensional cohesive crack method for reinforced concrete structures[END_REF], etc.) to become useful in practice, it is thus important to devise techniques which are able to significantly reduce the computational effort required without sacrificing accuracy.

Historically, reducing the computational time associated with solving nonlinear problems in solid mechanics has mainly been addressed by developing homogenisation techniques [START_REF] Fish | Computational plasticity for composite structures based on mathematical homogenization: Theory and practice[END_REF][START_REF] Feyel | FE 2 multiscale approach for modelling the elastoviscoplastic behaviour of long fibre SiC/Ti composite materials[END_REF]. In this case, the material properties associated with a material point in a coarse representation of the structure is obtained by averaging of the fine scale material behaviour over a "representative volume element" which suppresses the need to resolve the fine scale explicitly over the (possibly large) structure. These methods are widely adopted (see e.g. for a recent description [START_REF] Xiao | Two-scale asymptotic homogenisation-based finite element analysis of composite materials[END_REF]) and are rather effective in decreasing the computational costs associated with solving fine scale problems.

However, the application of homogenisation-based methods to failure simulations is not straightforward, as the assumptions required to prove the separation of scales are not satisfied. This may explain why such methods are still under-developed for fracture problems. The ability to coarsegrain fracture information would however lead to a further decrease of computational time by concentrating the computational effort where it is needed most. In the wake of a propagating crack, for example, the full resolution used to resolve fine scale effects close to the crack front is no longer necessary, and a coarse description may be sufficient. This issue is being intensively investigated by several teams, namely [START_REF] Belytschko | Multiscale aggregating discontinuities: A method for circumventing loss of material stability[END_REF] who developped a coarse graining method called Multiscale Aggregating Discontinuities, and [START_REF] Nguyen | Homogenization-based multiscale crack modelling: From micro-diffusive damage to macro-cracks[END_REF][START_REF] Cazès | A thermodynamic method for the construction of a cohesive law from a nonlocal damage model[END_REF][START_REF] Verhoosel | Computational homogenization for adhesive and cohesive failure in quasi-brittle solids[END_REF] who were able to derive an homogenisation scheme for lumping the fine scale deffects of Representative Volume Elements (RVEs) into cohesive cracks introduced at the coarse scale. The solutions proposed so far are however highly problem-dependent and usually do not inherit the solid mathematical foundations from homogenisation theory.

We propose here to follow an alternative route which relies on very recent advances in model order reduction strategies for highly nonlinear structural problems [START_REF] Barrault | An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations[END_REF][START_REF] Ryckelynck | A priori hyperreduction method: an adaptive approach[END_REF][START_REF] Niroomandi | Model order reduction for hyperelastic materials[END_REF][START_REF] Verdon | Reduced-order modelling for solving linear and non-linear equations[END_REF][START_REF] Kerfriden | Bridging proper orthogonal decomposition methods and augmented newton-krylov algorithms: an adaptive model order reduction for highly nonlinear mechanical problems[END_REF]. We believe that, though rarely studied in the context of failure simulation in nonlinear mechanics, this avenue provides a possibly interesting alternative (or complementary) approach to homogenisation-based solution procedures. This paper is in substance a significant improvement of the work proposed in [START_REF] Kerfriden | Bridging proper orthogonal decomposition methods and augmented newton-krylov algorithms: an adaptive model order reduction for highly nonlinear mechanical problems[END_REF], dedicated to the fast prediction of the onset of structural failure by adaptive model order reduction. We aim here at preserving the computational efficiency and accuracy obtained in [START_REF] Kerfriden | Bridging proper orthogonal decomposition methods and augmented newton-krylov algorithms: an adaptive model order reduction for highly nonlinear mechanical problems[END_REF] when the size of the damaged zone becomes significant compared to the overall structure.

The definition of model order reduction can vary depending on authors. We speak here of an approximation of the unknown field in a coarse space spanned by a basis of global vectors (as opposed to locally supported shape functions, in finite elements for example). Model order reduction (MOR) finds its roots in the systems engineering community, where the issue of obtaining a fast and reliable approximation of transfer functions has been extensively addressed over the years. The interested reader can refer to the review proposed in [START_REF] Antoulas | Approximation of large-scale dynamical systems: an overview[END_REF]. In the context of solid mechanics, most studies have focused on the reduction of dynamics problems. Starting from classical modal basis [START_REF] Rayleigh | The Theory of Sound[END_REF] in a Ritz-Galerkin framework, two main classes of model order reduction methods have been proposed. The first propose a truncation of the modal basis and an approximation of the highfrequency behaviour by static corrections (see for instance [START_REF] Dickens | A critique of mode acceleration and modal truncation augmentation methods for modal response analysis[END_REF]). The second use the idea of domain partitioning. Slave degrees of freedom are eliminated by Guyan reduction (or static condensation). The famous Craig-Bampton method [START_REF] Craig | Coupling of substructures for dynamic analysis[END_REF] is an evolution of this idea, using truncated modal bases to enhance static condensation. These techniques have been extended to the approximation of nonlinear dynamical behaviour (an interesting review is proposed in [START_REF] Kerschen | Nonlinear normal modes, part i: A useful framework for the structural dynamicist[END_REF]).

The literature concerning the application of model order reduction to the approximation of the (quasi-)static behaviour of nonlinear structures is a lot sparser. Different routes have been investigated to develop such reduced models. The first family of methods consists in taking advantage of the a priori knowledge of a set of representative solutions to fasten the solution process. These representative responses, called snapshots [START_REF] Sirovich | Turbulence and the dynamics of coherent structures. part I: coherent structures[END_REF], can be for instance a set of solution vectors obtained for given values of structural parameters. The Proper Orthogonal Decomposition (POD) [START_REF] Karhunen | Über lineare methoden in der wahrscheinlichkeitsrechnung[END_REF][START_REF] Loeve | Probability Theory. 3rd edition[END_REF] can classically be used to extract a relevant global basis from these snapshot vectors, which in turn defines a coarse space where the solution to the structural problem with a different set of parameters can be cheaply searched. The second family is not exactly based on reduction by projection. Spectral approaches referred to as Proper Generalized Decomposition (PGD) are a family of solvers capable of building the solution to a multidimensional reference problem as a sum of products of functions of separate variables. The dimensions can be space or time, but also structural parameters. Once this solution is obtained offline, the online computation simply consists in the evaluation of the separated functions, which are known explicitly. These novel methods are very appealing and were recently applied to linear, stochastic and/or time dependant problems [START_REF] Nouy | A generalized spectral decomposition technique to solve a class of linear stochastic partial differential equations[END_REF][START_REF] Chinesta | Recent advances and new challenges in the use of the proper generalized decomposition for solving multidimensional models[END_REF] (including linear predictions of nonlinear solvers [START_REF] Ladevèze | The latin multiscale computational method and the proper generalized decomposition[END_REF]), and to non-linear partial differential equations [START_REF] Nouy | Generalized spectral decomposition method for stochastic non linear problems[END_REF][START_REF] Chinesta | Recent advances and new challenges in the use of the proper generalized decomposition for solving multidimensional models[END_REF]. Their extension to irreversible constitutive laws such as the ones considered in damage mechanics does not seem to be straightforward, and, as far as we know, has never been adressed before (it is therefore a very interesting possible extension of the present piece of work).

We focus here on the first family of solvers. To the authors' knowledge, the application of model order reduction methods to failure simulations has virtually never been addressed. The main reason for this is that the solution fields in the zones of the structure where damage or failure initiate and propagate are highly parameter-dependent, and cannot be known in advance. Hence, precomputed reduced bases are unable to represent the structural behaviour beyond the onset of strain localisation, damage and fracture. In the particular case of linear elastic fracture mechanics, the authors of [START_REF] Galland | A global model reduction approach for 3d fatigue crack growth with confined plasticity[END_REF] have proposed a mesh morphing technique that allows to keep an autosimilar description of the singularity during propagation. In this particular case, the use of a snapshot POD is relevant and allow for interesting computational savings. Corrective POD algorithms have also been proposed to alleviate this in the case of plasticity or damage models [START_REF] Ryckelynck | A priori hyperreduction method: an adaptive approach[END_REF][START_REF] Kerfriden | Bridging proper orthogonal decomposition methods and augmented newton-krylov algorithms: an adaptive model order reduction for highly nonlinear mechanical problems[END_REF]. Both these techniques are based on global corrections of the initial reduced basis by Krylov algorithms. Though efficient for global nonlinearities, it was shown in [START_REF] Kerfriden | Bridging proper orthogonal decomposition methods and augmented newton-krylov algorithms: an adaptive model order reduction for highly nonlinear mechanical problems[END_REF] that in the case of damage assessment, the number of corrections increases when nonlinearities localise, i.e.: when failure occurs. The purpose of this paper is to circumvent this phenomenon.

The first stage of the proposed method is to only reduce the set of balance equations which exhibits a "smooth" nonlinearity. In the context of the simulation of damage, one part of the degrees of freedom, corresponding to highly damaged zones, will be considered as master degrees of freedom, which will be fully solved for, while the remaining will be approximated as a linear combination of reduced basis vectors, which may be precomputed and corrected on the fly. In this respect, the proposed method has strong links with Craig-Bampton methods, or model order reduction by substructure such as the one proposed in [START_REF] Rixen | A dual craig-bampton method for dynamic substructuring[END_REF][START_REF] Markovic | Partitioning based reduced order modelling approach for transient analyses of large structures[END_REF][START_REF] Ladevèze | The latin multiscale computational method and the proper generalized decomposition[END_REF]. A similar approach was also proposed in [START_REF] Ammar | Coupling finite elements and reduced approximation bases[END_REF]. The second stage of our developments consists in developing an efficient solver for the linearized system of equations obtained by a Newton solution algorithm for the reduced problem. The system is condensed on the master degrees of freedom, and solved by a Krylov algorithm. The global reduced basis is used to efficiently solve the consensed problem: an augmentation technique ensures that the iterations are performed in a space which is orthogonal to the one spanned by the reduced basis. As a result, the augmented Krylov iterations correspond to a two-level solver, the coarse level being the space spanned by the restriction of reduced basis vectors to the fully resolved set of degrees of freedom.

The paper is organized as follows. In the first section, we recall the basics of projection-based model order reduction applied to the solution of nonlinear quasi-static problems. In section 3, we formally describe the proposed local/global reduction algorithm. This technique is applied to the reduction of damage models in section 4. Finally, we show in the last section that this technique permits to drastically reduce the number of corrections of the reduced basis required to obtain a desired level of accuracy.

Model order reduction by projection: general concepts

Projection-based model order reduction can be introduced as a means to obtain, at low computational cost, a good approximation to a set of nonlinear equations of the form:

F (X i ) i∈ 0,N = 0 (1)
Where F is a vectorial, possibly nonlinear, function of a set of vectorial state variables X. The expression of this function can take various forms, depending on the application.

Problem statement in nonlinear continuum mechanics

We consider a structure occupying a continuous domain Ω with boundary ∂Ω. It is subjected to prescribed displacements U D on its boundary ∂Ω u and to prescribed tractions F D on the complementary boundary ∂Ω f = ∂Ω\∂Ω u , over time interval [0, T ]. Let u be the unknown displacement field, which belongs to the space U of kinematically admissible fields:

U = u ∈ H 1 (Ω) | u |∂Ωu = U D (2)
Let U 0 be the associated vector space.

U 0 = u ∈ H 1 (Ω) | u |∂Ωu = 0 (3)
Under the assumptions of quasi-static evolution of the structure and small perturbations, the weak form of the equilibrium and constitutive laws read, at any time t ∈ [0, T ]: ∀u ∈ U 0 , find u ∈ U such that:

Ω σ : (u ) dΩ = Ω f D .u dΩ + ∂Ω f F D .u dΓ σ = σ (u) |τ τ ≤t ( 4
)
where σ is the Cauchy stress tensor and (u) is the symmetric part of the displacement gradient.

The constitutive relation between σ and (u) is nonlinear and described using internal variables (plasticity, damage for instance). It is assumed local and rate-independent.

Space and time discretization

Let us perform a standard finite element approximation of the space of unknown displacement field U (and a similar approximation of the space of test functions U 0 ):

U h (Ω) = u(x) | u(x) = nn i=1 N i (x) u i (5)
where n n is the number of nodes, x is the position vector, (N i ) i∈ 1,nn are the finite element shape functions associated to each node of the finite element mesh, and (u i ) i∈ 1,nn are the nodal values of the displacement field.

The nonlinear solution strategy used to solve the problem over time is a classical time discretization scheme for quasi-static and rate-independent problems. This procedure consists in finding a set of consecutive solutions at times (t n ) n∈ 0,nt .

The introduction of the finite element approximation and time discretization into equation (4) at any time t n+1 of the analysis leads to the following nonlinear vectorial equation:

F Int ∆U, U |tm m∈ 0,n + F Ext = 0 (6)
where the vector of increment in the nodal unknowns ∆U ∈ R nu (n u is the number of scalar nodal unknowns) is defined by ∆U = U |tn+1 -U |tn , F Int ∈ R nu and F Ext ∈ R nu are respectively the internal forces resulting from the discretization of the internal virtual work (left-hand side of the first equation of system (4)) and the external forces resulting from the discretization of the external virtual work (right-hand side of the first equation of system (4)). In the following, the dependency of the internal forces with respect to the history of the unknown fields will not be written explicitly, unless necessary. The set of successive solution vectors to problem (6) at times (t n ) n∈ 0,nt , obtained by a tangent Newton algorithm, will be considered as the reference solution. The algorithms used in the following should provide a solution that is closer to the reference solution, but for a cheaper computational cost. Therefore, the relevance of (i) the initial space and time discretizations and (ii) the chosen Newton solution algorithm will not be discussed in this paper.

Reduced solution of the balance equations by projection

Principle

The solution vector is searched for in a space of small dimension (several orders smaller than the number of finite element degrees of freedom). Let us call C the matrix whose columns form a basis of this space:

C = C 1 C 2 ... C nc (7)
where n c is the dimension of the reduced space, and (C k ) k∈ 1,nc ∈ (R nu ) nc are the basis vectors. Applied to the reduction of problem [START_REF] Rabczuk | A geometrically non-linear threedimensional cohesive crack method for reinforced concrete structures[END_REF], the increment in the solution field is approximated by:

∆U = C α (8)
where we introduced the reduced state variables α ∈ R nc . Problem ( 6) is now overconstrained. One usually tackle this problem by prescribing a Galerkin orthogonality condition: the residual of equation ( 6) R = F Int (U |tn + C α) + F Ext is constrained to be orthogonal to any test vector δU = C δα belonging to the space spanned by the reduced basis vectors (C k ) k∈ 1,nc :

C T R = 0 (9)
The reduced form of problem ( 6) is:

C T F Ext + F Int (U |tn + C α) = 0 (10)
2.2.2 One particular choice of projection space: the proper orthogonal decomposition (POD)

Various choices of projection spaces are proposed in the literature. One of the most successful, as far as nonlinear simulations are concerned, is the snapshot-POD [START_REF] Sirovich | Turbulence and the dynamics of coherent structures. part I: coherent structures[END_REF]. This particular technique requires the knowledge of a representative family of solutions to the global problem. This set of vectors (S k ) k∈ 1,ns is called snapshot. The aim is to find an orthonormal basis (C k ) k∈ 1,nc , of dimension n c smaller than n s such that the distance between spaces Im(C) and Im(S) is minimum. This distance can be quantified by the following metric:

J(C 1 , . . . , C nc ) = ns j=1 S j - nc i=1 C i T S j C i 2 (11)
which must be minimized under the constraint of orthonormality of (C k ) k∈ 1,nc . Problem [START_REF] Nguyen | Homogenization-based multiscale crack modelling: From micro-diffusive damage to macro-cracks[END_REF] can be reformulated as an unconstrained formulation:

L(C 1 , . . . , C nc , λ 11 , . . . , λ ncnc ) = J(C 1 , . . . , C nc ) + nc i=1 nc j=i λ ij C i T C j -δ ij (12)
where (λ ij ) i,j∈ 1,nc are Lagrange multipliers and δ ij the Kronecker symbol. The n c (n c +1)/2 optimality conditions for maximizing the Lagrangian with respect to the lagrange multiplier naturally ensure the orthogonality of the basis:

∂L ∂λ ij = 0 → C i T C j = δ ij (13)
Whereas the n c optimality conditions for minimizing the Lagrangian with respect to the basis vectors reads:

∂L ∂C i = 0 → ns j=1 S j C i T S j = λ ii C i and λ ij = 0 for i = j (14)
Defining λ i = λ ii and S = S 1 S 2 ... S ns ∈ R nu × R ns , equations ( 14) can be written as the following n u × n u eigenvalue problem:

S S T C i = λ i C i ( 15
)
This problem is equivalent to computing the Singular Value Decomposition (SVD) of S, where the singular values s i are such that s 2 i = λ i ≥ 0. If n s n u (which is the case in general when applying the snapshot-POD), replacing problem [START_REF] Ryckelynck | A priori hyperreduction method: an adaptive approach[END_REF] by the eigenvalue problem on the n s × n s operator S T S is computationally cheaper, and reads:

C i = λ -1 2 i S V i with S T S V i = λ i V i (16)
The error associated with the truncation of the SVD at rank n c is quantifyed by the distance J(C 1 , ..., C nc ) which can be shown to be equal to the sum of the truncated eigenvalues:

J(C 1 , ..., C nc ) = ns i=nc+1 λ i (17)
This result provides a precious tool for choosing the most relevant basis. Indeed, only the basis vectors C i associated with the largest λ i (in practice, those such that λ i /λ max > ε, ε being a small parameter depending on the required accuracy) are selected as reduced basis vectors. This ensures that for an given number n c of basis vectors, the approximation of the snapshot space is maximised in the sense of (11).

Nonlinear solution algorithm

At any time step, problem [START_REF] Belytschko | Multiscale aggregating discontinuities: A method for circumventing loss of material stability[END_REF] can be solved by classical Newton algorithms (successive linearizations). Iteration i + 1 of the tangent Newton method consists in solving the following linear predictor:

K i T,R δα i+1 = -R i R ( 18
)
where the increment in the reduced state variables is δα i+1 = α i+1α i , and the residual and tangent operator obtained from the knowledge of the solution at iteration i read (see [START_REF] Kerfriden | Bridging proper orthogonal decomposition methods and augmented newton-krylov algorithms: an adaptive model order reduction for highly nonlinear mechanical problems[END_REF] for more details):

     R i R = C T F Int (U |tn + C α i ) + F Ext K i T,R = C T ∂F Int (U |tn + C α) ∂α α=α i (19)
3 Local/global model order reduction strategy

Projection-based model order reduction techniques introduce a global approximation of the displacement. As such, even adaptive versions [START_REF] Ryckelynck | On the a priori model reduction: overview and recent developments[END_REF][START_REF] Kerfriden | Bridging proper orthogonal decomposition methods and augmented newton-krylov algorithms: an adaptive model order reduction for highly nonlinear mechanical problems[END_REF] of these methods are not well suited to the analysis of localised nonlinearities. We propose in the following to use reduction techniques by projection for the "weakly nonlinear" equations of reference problem (6), namely equations for which reduced order modelling (ROM) is relevant, while the remaining equations will be solved directly.

Displacement approximation

The principle of the proposed local/global strategy is to split the unknown solution vector into two parts, only one of them being approximated as a linear combination of reduced basis vectors. We shall use superscript (r) for the approximated "slave" part of the solution vector, while superscript (f ) (which stands for "fully resolved") will be used to denote its complementary "master" part. More precisely, let ∆U be the unknown increment vector at time t n+1 , the numbering being reorganised as follows:

∆U = ∆U (r) ∆U (f ) = E (r) E (f ) ∆U (20)
In the above notations, ∆U (r) ∈ R nr and ∆U

(f ) ∈ R n f , with n r + n f = n u and n f much smaller than n r . E (r) ∈ {0, 1} nr × {0, 1} nu and E (f ) ∈ {0, 1} n f × {0, 1}
nu are two boolean extractors (rectangular matrices with one 1 per line, the other coefficients being null).

∆U (r)
is approximated as a linear combination of global vectors

(C k ) k∈ 1,nc ∈ (R nu ) nc : ∆U (r) = E (r) C α (21)
where α ∈ R nc is the vector of reduced degrees of freedom. In the initial numbering, the unknown vector can be expressed in the form:

∆U = ∆U (r) + ∆U (f ) where    ∆U (r) = E (r) T ∆U (r) ∆U (f ) = E (f ) T ∆U (f ) (22)
Introducing projector P (r) = E (r) T E (r) , the approximation of the displacement increment finally reads:

∆U = P (r) C α + E (f ) T ∆U (f ) (23)
For the sake of legibility, we define the new state vector X at time t n+1 :

X = α ∆U (f ) (24)
Using the above definition, the approximation of the displacement increment is now:

∆U = A X where A = P (r) C E (f ) T (25)

Locally reduced set of balance equations

At that point, we need to define how the overconstrained balance equations ( 6) will be solved for. Let us restrict our numerical investigations to the Galerkin procedure. The balance equation are therefore required to be orthogonal to any test vector:

δU = A α δU (f ) (26)
where δα ∈ R nc and δU (f ) ∈ R n f are arbitrary test state variables. The orthogonality requirement applied to the initial nonlinear set of equations ( 6) yields the following nonlinear system:

R R (X) = A T F Int ∆U(X) + U |tn + F Ext = 0 (27)

Newton solution scheme

At each time step of the time discretization scheme, the reduced problem ( 27) is solved by a Newton algorithm. The (i + 1) th Newton iteration consists in finding δX i+1 satisfying the following set of

n c + n f linearized equation: ∂R R (X) ∂X X=X i δX i+1 = -R i R (28)
where δX i+1 = X i+1 -X i and R i R is the residual of problem ( 27) computed at iteration i:

R i R = R R (X i ) ( 29
)
The expression of the tangent operator

K i T,R = ∂R R (X) ∂X X i
is obtained by differentiation of equation ( 27) with respect to the set of reduced state variables, in any arbitrary direction δX such that:

δX = δα δU (f ) (30)
This differentiation reads:

K i T,R δX = ∂R R (X) ∂X X i δX (31)
Using the tangent operator of the full system of equations ( 6), and assuming that ∂F Ext ∂∆U = 0, on can write:

K i T,R δX = A T ∂F Int (∆U) ∂∆U ∆U=A X i ∂∆U ∂X ∆U=A X i δX (32)
Using the notation

K i T = ∂F Int ∂∆U ∆U(X i )
, the previous equation yields the expression of the reduced tangent operator as a function of the tangent of the full system of equations:

K i T,R δX = A T K i T A δX (33)
Introducing the above expression of the tangent into equation ( 28), and expanding the expression of operator A, one finally obtains the following linearized system:

C T P (r) K i T P (r) C C T P (r) K i T E (f ) T E (f ) K i T P (r) C E (f ) K i T E (f ) T δX i+1 = - C T P (r) E (f ) F Int ∆U(X i ) + U |tn + F Ext (34)
Notice that problem [START_REF] Ammar | Coupling finite elements and reduced approximation bases[END_REF] is considerably reduced compared to a direct linearisation of (6). Indeed,

K i T,R ∈ R nc × R n f
where n c n u and, applied to the analysis of localised nonlinearities, we can reasonably expect that n f is at least one order of magnitude smaller than n u , which means that only a few of the balance equations are strongly nonlinear. This assumption will be validated in section 4.

This framework is very general. It does not rely on domain decomposition but on a splitting of the balance equations. Yet, this technique will be particularized in section 4, and we will show that using the ideas inspired from domain decomposition methods [START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF][START_REF] Kerfriden | A three-scale domain decomposition method for the 3D analysis of debonding in laminates[END_REF] and local/global methods [START_REF] Gendre | Non-intrusive and exact global/local techniques for structural problems with local plasticity[END_REF][START_REF] Allix | A relocalization technique for the multiscale computation of delamination in composite structures[END_REF] provides efficient splitting in the case of damage mechanics.

Two-level solution of the linearized systems

The successive solution to (28) could be obtained by a direct solver. This would not affect the results given in the following. However, we believe that, for this particular application, an iterative solver has several advantages:

• Iterative solvers, such as a conjugate gradient (see [START_REF] Saad | Iterative Methods for Sparse Linear Systems Iterative Methods for Sparse Linear Systems[END_REF] for instance), do not require any factorisation. Factorisations usually deteriorate the sparsity of algebraic systems, which results in an increase in the memory requirements. Newton-Krylov methods do not, in theory [START_REF] Knoll | Jacobian-free newton-krylov methods: a survey of approaches and applications[END_REF], even require the assembly of the stiffness operator.

• Direct solvers are less versatile than iterative ones, and are usually a limiting point to extend numerical solution schemes to parallel computing.

• Solving exactly the successive linearized systems ( 28) is not required. This idea has led to intensive work on Inexact Newton Methods [START_REF] Dembo | Inexact newton methods[END_REF], but this issue will not be addressed in this paper. The interested reader can refer to [START_REF] Kerfriden | A three-scale domain decomposition method for the 3D analysis of debonding in laminates[END_REF][START_REF] Kerfriden | Bridging proper orthogonal decomposition methods and augmented newton-krylov algorithms: an adaptive model order reduction for highly nonlinear mechanical problems[END_REF] for studies on similar topics in the scope of domain decomposition and model order reduction.

• More importantly, a basic direct solver such as Cholesky factorizations will not make use of any mechanically relevant information obtained during the previous iterations of the nonlinear solver. In the chosen iterative algorithm, the part of the displacement field that is "fully resolved" will be solved for efficiently using the reduced basis as a preconditioner, following the ideas proposed for instance in [START_REF] Risler | Iterative accelerating algorithms with krylov subspaces for the solution to large-scale nonlinear problems[END_REF][START_REF] Knoll | Jacobian-free newton-krylov methods: a survey of approaches and applications[END_REF][START_REF] Kerfriden | A three-scale domain decomposition method for the 3D analysis of debonding in laminates[END_REF].

We use here a Krylov algorithm to efficiently find a solution to linearized system of equation [START_REF] Chinesta | Recent advances and new challenges in the use of the proper generalized decomposition for solving multidimensional models[END_REF]. Three stages of preconditioning are performed: a condensation, a projection and a crude preconditioning of the resulting system.

Condensation

Let us condense the linearized problem [START_REF] Ammar | Coupling finite elements and reduced approximation bases[END_REF] on the master (i.e.: "fully resolved") degrees of fredom:

S (f ) P δU (f ) = R (f ) C (35)
where the superscript corresponding to the current and previous Newton iterations have been dropped for the sake of legibility. S (f ) P is the primal Schur complement and R (f )

C the condensed residual. They are defined by:

     S (f ) P = K (f f ) T,R -K (f r) T,R K (rr) T,R -1 K (rf ) T,R R (f ) C = K (f r) T,R K (rr) T,R -1 R (r) R -R (f ) R (36)
Where we have used the usual block notations:

K T,R = K (rr) T,R K (rf ) T,R K (f r) T,R K (f f ) T,R (37)
R R = R (r) R R (f ) R (38)
The state variables corresponding to the reduced part of the displacement vector increment have been eliminated from the system of equations and can be retrieved as follows:

δα = -K (rr) T,R -1 R (r) R + K (rf ) T,R δU (f ) (39)
K (rr) T,R -1
is an operator of very small size, which can be explicitly computed prior to applying the iterative Krylov solver. However, in the case where a very large number of reduced basis vectors is used to approximate the unknown field in the smooth region, computing the Schur-complement/vector products corresponding to the construction of the successive Krylov basis vectors will become computationally more efficient.

Iterative solution to the condensed problem

The classical projected (or augmented) conjugate gradient [START_REF] Dostál | Conjugate gradient method with preconditioning by projector[END_REF][START_REF] Farhat | Implicit parallel processing in structural mechanics[END_REF] is applied to the approximate solution of (35) (the linearized operator is assumed symmetric, positive and definite). The chosen augmentation space is the space spanned by the reduced basis vectors Im(C (f ) ) with the restriction

C (f ) = E (f ) C.
The starting point of augmented Krylov solvers is to decompose the unknown solution increment into two supplementary spaces:

δU (f ) = δU (f ) C + δU (f ) K where    δU (f ) C ∈ Im C (f ) δU (f ) K ∈ Im(C (f ) ) ⊥ S (f ) P = Ker C (f ) T S (f ) P (40)
⊥ S (f ) P designing the S (f ) P -orthogonality, which is ensured by introducing the classical projector:

P = I d -C (f ) (C (f ) T S (f ) P C (f ) ) -1 C (f ) T S (f ) P ( 41
)
and writing that δU

(f ) K = P δU (f )
K . This separation of the search space into two subspaces Im C (f ) and Im(P) in direct sum leads to the following uncoupled equations:

δU (f ) C = C (f ) (C (f ) T S (f ) P C (f ) ) -1 C (f ) T R (f ) C P T S (f ) P P δU (f ) K = P T R (f ) C (42)
The first line of equation ( 42) is a coarse initialization of the projected conjugate gradient. The second line is the linear prediction of the full problem projected on Im(C (f ) ) ⊥ S (f ) P . This system is symmetric and can be solved by a preconditioned conjugate gradient. Hence one solves iteratively:

M -1 P T S (f ) P P δU (f ) K = M -1 P T R (f ) C (43)
where

M -1
is the left preconditioner (symmetric, definite and positive). In our test cases, M is a diagonal matrix whose entries are the elements of the diagonal of S (f ) P . Algorithm 1 gives the classical implementation of an augmented preconditioned conjugate gradient.

Algorithm 1: Augmented preconditioned conjugate gradient for the solution to the condensed linearized reduced problem Compute

             S (f ) P = K (f f ) T,R -K (f r) T,R K (rr) T,R -1 K (rf ) T,R R (f ) C = R (f ) R -K (f r) T,R K (rr) T,R -1 R (r) R C (f ) = E (f ) C (C (f ) T S (f ) P C (f ) ) -1 δU (f ) C = C (f ) (C (f ) T S (f ) P C (f ) ) -1 C (f ) T R (f ) C ; R CG,0 = R (f ) C -S (f ) P δU (f ) C = P T R (f ) C ; Z 0 = P M -1 R CG,0 , W 0 = Z 0 and δU (f ) K,0 = 0; for j = 1, . . . , n do α j-1 = (R CG,j-1 , W j-1 )/(S (f ) P W j-1 , W j-1 ) δU (f ) K,j = δU (f ) K,j-1 + α j-1 W j-1 R CG,j = R CG,j-1 -α j-1 S (f ) P W j-1 Z j = P M -1 R CG,j β j = (S (f ) P W j-1 , Z j )/(W j-1 , S (f ) P W j-1 ) W j = Z j -β j W j-1 end

Interpretation and comments on the expected additional costs

The three preconditioning techniques that we have detailed in the previous sections define a cheap structural preconditioner. A schematic representation of this algorithm is given in figure [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF].

From an algebraic point of view, it can easily be proved that the condensation has the same effect on the conjugate gradient iterations as a left block preconditioner for system [START_REF] Ammar | Coupling finite elements and reduced approximation bases[END_REF]:

M Condens -1 = K (rr) T,R -1 0 0 0 (44)
From a mechanical point of view, condensing the problem on the master degrees of freedom ensures that any solution provided by the iterative solver exactly solves the reduced part of the balance equations.

The augmentation of the conjugate gradient can also be seen as a left preconditioner for the condensed system of equations (second line of system (42)), which is associated with an initialization of the iterative alogrithm. The initialization δU (f ) C (first line of system [START_REF] Dembo | Inexact newton methods[END_REF]) is the best solution of the condensed problem in the restriction of the pre-computed reduced basis to the master degrees of freedom (in the sense of a (S (f ) P ) -1 -norm). Taking into account the previous comments on the interpretation of the condensation step, this initialisation is exactly the solution to the linearized

δU = E (f ) T δU (f ) + P (r) C δα δU f C = δU (f ) C + δU (f ) K j ← j + 1 j = 0
F D Z j = P M -1 R CG,j
Coarse initialisation in the reduced space system of equation ( 18) on the master degrees of freedom. In other words, the initialization of the iterative algorithm provides the solution of the reduced linearized problem which would have been obtained without local/global enhancement. No significant additional costs compared to the basic projection-based reduced order modelling are involved so far. In the case where the solution which is looked for is correctly approximated in the reduced space (e.g.: "smooth" nonlinearity, or unloading of the structure), no conjugate gradient iteration is required.

δU f C ∈ Im E (f ) C ||R CG,j || ||R (f ) C || < R CG,j δU (f ) K,0 = 0 δU (f ) K,j = δU (f ) K,j-1 + α j W j
If conjugate gradient iterations are required, the augmentation acts as a global correction of the successive search directions (see figure [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF]). It ensures that any solution of the form [START_REF] Saad | Iterative Methods for Sparse Linear Systems Iterative Methods for Sparse Linear Systems[END_REF] still satisfies the reduced balance equations without local/global enhancement. Therefore, the complementary solution δU (f ) K can be interpreted as a local correction to the solution provided by the pre-computed reduced model. We will see in the examples that, as the conjugate gradient does not need to search for the part of the solution that is directly captured in the pre-computed reduced space, the convergence rate is very high.

The diagonal preconditioner M -1 classically accounts for local heterogeneities in the structure. A much more sophisticated preconditioner can of course be used (e.g.:incomplete Cholesky factorization). Yet, if the reduced basis is correctly pre-computed and/or updated, the computational effort required to obtain a good solution of the projected problem should be relatively low. Hence, the potentially expensive construction of an advanced preconditioner is not a priori justified. The example section will show that the number of conjugate gradient iterations required to compute a correction to the reduced model is indeed kept very low, even while only using this crude preconditioner. ) are the two extremities of the bar as defined in figure [START_REF] Allix | Modeling and simulation of crack propagation in mixed-modes interlaminar fracture specimens[END_REF]. The displacement field in is supposed constant in a cross-section. Therefore, the displacement of any point M ∈ Ω b can be expressed in the following form:

u |M ∈Ω b = (u(x).n kl ) n kl = u(x) n kl (45)
where each cross section is parametrised by a scalar x ∈ [0, P k P l ], defined by x = P k M .n kl . We use a constitutive law based on classical damage mechanics to describe the behaviour of the lattice structure. The lineic strain energy of bar b reads:

e d|x = 1 2 E(1 -d) S b 2 |x (46)
where |x = u ,x (x) is the strain measured in the direction of the bar, (S b ) b∈ 1,n b , is the section of the bar, and d is a damage variable which ranges from 0 (undamaged material) to 1 (completely damaged material). The local state equations, derived from the strain energy read:

     N |x = ∂e d ∂ = E(1 -d) S b |x Y |x = - ∂e d ∂d = 1 2 E S b 2 |x (47)
where we have introduced the axial force N = S b n T kl .σ.n kl , and the thermodynamic force Y . A non-local evolution law is defined to link the damage variables to the thermodynamic forces: We suppose that the lattice is not subjected to any volume force, that the material of each elementary bar is isotropic and homogeneous, and that the section of each bar is constant. As a consequence, the resulting displacement field evolves linearly along the direction of each bar.

d |x = max x∈[0, P k P l ] min 1, max τ ≤t α(Y |x,τ ) β (48)
P k P l M x n kl
In terms of finite element discretization, each bar will be considered as a linear element

u(x) = 1 - x P k P l u |P l + x P k P l u |P k (49)
Because of the assumptions made previously on the loading, geometry and material properties, the finite element solution yields the "exact" solution to the damageable lattice problem.

Lattice structure and loading

We solve the problem defined in figures [START_REF] Jefferson | A model for cementitious composite materials based on micro-mechanical solutions and damage-contact theory[END_REF]. The sections and Young's modulus are unitary. The material parameters are β = 0.5 and α = √ 2. Any bar parallel to axis (O, x), (O, y) or (O, z) has a unitary length.

The Dirichlet boundary conditions applied to the structure read, for any node (P i ) (with i ∈ 1, n p ) of the lattice such that (P i ) ∈ (O, x, y):

u |Pi .z = 0 (50)
As represented in figure (3), homogenous Neuman boundary conditions, in the -z direction, are applied to any point P i that belongs to plane P defined by: Increasing the value of the applied forces leads to the initiation and propagation of a damaged zone as shown in figure (4). This phenomenon is unstable from a structural point of view, and the damage state represented in this figure is past the limit point at which the force applied to the structure reaches its maximum value. An arc-length algorithm with local control is used to solve the damage propagation beyond this limit point. The basic idea of this algorithm is to look for the amplitude of the force under the constraint that the maximum increment in the damage variables over a time step is fixed (see [START_REF] Riks | The application of newton's methods to the problem of the application of newton's methods to the problem of elastic stability[END_REF][START_REF] Allix | Modeling and simulation of crack propagation in mixed-modes interlaminar fracture specimens[END_REF] or [START_REF] Kerfriden | Bridging proper orthogonal decomposition methods and augmented newton-krylov algorithms: an adaptive model order reduction for highly nonlinear mechanical problems[END_REF] for additional details). The reference solution is obtained by solving the full problem (6) in 30 of these damage increments of the Newton/arc-length algorithm. The final state is represented in figure [START_REF] Ladevèze | An enhanced mesomodel for laminates based on micromechanics[END_REF].

P = {M = (x, y, z) | x ∈ [7,
Our purpose is to check the adequacy of the proposed local/global reduction technique when trying to obtain at cheap costs the load/deflection curve (pre and post failure phases) of the structure under the loading conditions described previously. The results provided in the following sections give the value of the norm of the external forces vector (also called loading factor) obtained by the arc-length method. The displacement used to plot the load/deflection curves is the average vertical displacement of the points at which forces are applied. We need at that point to particularize the general method developed in section 3. The splitting between reduced degrees of freedom ∆U (r)

and the complementary ones ∆U

(f )
is made by considering that the part of the structure undergoing strong damage variations will not be correctly solved for when projected on a pre-computed global basis (see [START_REF] Gendre | Non-intrusive and exact global/local techniques for structural problems with local plasticity[END_REF][START_REF] Kerfriden | A three-scale domain decomposition method for the 3D analysis of debonding in laminates[END_REF] for similar assumptions in a domain decomposition-based local/global approaches). This statement is purely based on mechanical intuition and will be validated in the following examples.

The following procedure is adopted. The splitting will be fixed over a time step, to avoid the stagnation or divergence of the Newton process. At the end of a time step, the element undergoing the maximum damage increment is spotted, and a sphere of radius ρ s , centred at the isobarycenter of the element is created. Every degree of freedom belonging to a node located inside this sphere is set as a fully resolved degree of freedom.

This procedure is repeated on the remaining elements (those which do not have all their degrees of freedom in the previous list).

The algorithm is stopped if one of the two following statements is satisfied.

• The maximum damage increment in the remaining elements is lower than k Dam times the global maximum. In our examples, k Dam is set to 0.5. This ensures that the fully resolved degrees of freedom are actually connected to elements undergoing a significant increase in their damage state compared to the remaining of the structure.

• The number of fully resolved degrees of freedom exceeds k LocGlo times the total number of unknowns. k LocGlo is set to 0.1 in our examples. This statement will be satisfied if the former one is not, which means that the structure undergoes a virtually homogeneous increase in its damage state. In that case, the local/global procedure might be inefficient.

To illustrate this procedure, two different reduced domains are represented in figure [START_REF] Grassl | Damage-plastic model for concrete failure[END_REF]. The dark elements are connected to at least one fully resolved node. The first one is obtained at the end of the third time step and is used to solve the fourth time step of the analysis, while the second one is used to solve the thirtieth and last time step. The first figure is obtained during the initiation phase. Notice that the maximum damage increments are found in the zone where the the load is applied, and in the "pillar" that is closer to this zone, which undergoes a significant bending loading. In the later case, corresponding to unstable propagation phase, the damage increment localises in the zone where the load is applied.

This empirical choice for the local/global splitting is only justified by its efficiency. However we only considered one type of model for brittle fracture. We expect the optimal splitting to depend on the chosen model and on the stage of the structural failure. A more suitable criterion should be investigated in future work, based on algebra to circumvent this dependency. For instance, one could consider the decrease rate of the singular values of the correlation matrix restricted to an assumed smoothed domain and loosely maximise it, with a greedy-like algorithm, by successive small modifications of this domain.

Modification of the reduced basis

At the end of a time step, the reduced basis is enriched in order to take into account the information provided by the local iterative solver. The procedure used here is simply to add the successive solutions to the initial snapshot, and sort them by a singular value computation. This step can be optimized, following our previous work [START_REF] Kerfriden | Bridging proper orthogonal decomposition methods and augmented newton-krylov algorithms: an adaptive model order reduction for highly nonlinear mechanical problems[END_REF]. The issue here is that the reduced basis is only locally enhanced, which might results in ill-posed reduced problems if not correctly addressed. This question will not be further discussed in this paper, which focuses on the benefits provided by the local/global approach.

A posteriori reduced basis

Reference solution to construct the a posteriori reduced basis

The first test consists in using as a snapshot the consecutive reference solutions obtained by a Newton/arc-length solution algorithm performed on the initial system of equations, without reduction. In other terms, we try to reproduce the snapshot by reduced order modelling. It terms of realistic applications, this example is not of particular interest as one would need to know the solution in advance in order to obtain the followinf results. However, studying this academic case is interesting to understand the behaviour of the proposed strategy.

The snapshot is composed of n s = 30 vectors. This information is compressed by using a singular value decomposition as detailed in section 2.2.2. We plot, in figure [START_REF] Rabczuk | A geometrically non-linear threedimensional cohesive crack method for reinforced concrete structures[END_REF], the normalised error in the snapshot representation

ν SVD =   ns j=1 S j - nc i=1 C i T S j C i 2   1 2   ns j=1 S j 2   1 2 = ns i=nc+1 λ i 1 2 ns i=1 λ i 1 2 (52)
as a function of the order of truncation of the SVD, where the eigenvalues (λ I ) i∈ 1,ns of the correlation matrix are now sorted in decreasing order. The decrease rate of this criterion is relatively low, which suggests that the information from the snapshot is poorly suited to compression. For reference, in [START_REF] Ryckelynck | A priori hyperreduction method: an adaptive approach[END_REF], the author estimates that the snapshot space is correctly approximated by the space spanned by the first n c singular vectors if ν SVD ≤ 10 -8 . In our case, this would lead to the selection of the 30 singular vectors to define a relevant reduced model. The singular vector basis is here truncated at order n c = 3 in order to yield a significant (thus observable) error in the solution In solid lines dark line, the set of 30 reference solutions are defined as the snapshot. In grey lines, a second snapshot operator is defined by concatenating the 30 reference solutions and the 30 solutions of the nearby problem. Both analyses are repeated when excluding the process zone (defined a priori by the set of points {M = (x, y, z) | x ∈ [5, 13] and y ∈ [6, 13]})) from the snapshot definition (i.e.: the snapshot now is the restriction of the solutions to the part of the domain which does not undergo significant damage). The resulting error curves are plotted in dashed lines. when using the basic POD method. This error can then be compared to the one obtained when using the proposed local/global scheme.

In figure [START_REF] Fish | Computational plasticity for composite structures based on mathematical homogenization: Theory and practice[END_REF], the reference load/deflection curved is compared to the one obtained when using the basic POD on the first hand, and when using the local/global reduction technique on the other hand. The peak load obtained with the basic POD is slightly underestimated, by less than 2%, but the overall behaviour of the structure is correctly captured. The local/global algorithm permits to further increase the accuracy of the approximated solution. The error in the maximum load factor is now less than 0.5 %.

We now analyse the effect of preconditioning the local iterations of the conjugate gradient by the projection technique detailed in section 3.4.2. The convergence curves that are plotted in figure (8, left) are obtained at the third Newton iteration of the fiftieth time step of the analysis (i.e.: during the damage propagation phase). The norm of the condensed system (35) of equation, normalized by the norm of the condensed right-hand term is observed as a function of the number of iterations of the conjugate gradient solver. One can see that the number of iteration to convergence is significantly decreased when using the augmented Krylov solver.

An other interesting observation is that the convergence curve decreases in a monotonic manner, which means that the error in the global residual keeps decreasing as a function of the cost of the local conjugate gradient. This is not, of course, a proof, and this conclusion is probably problem dependent, but the same behaviour as been observed in all our numerical experiments. A probable reason for this is that the augmentation deflates the spectrum of the condensed operator, removing its highest eigenvalues and forcing the Krylov algorithm to its super-convergence phase (one can refer to [START_REF] Gosselet | Méthodes de décomposition de domaine et méthodes d'accélération pour les problèmes multichamps en mécanique non-linéaire[END_REF] for a proper analysis of the effect of eigenvectors-based augmentation on the convergence of Krylov solvers). On the left-hand side, the reference solution is used as a snapshot. On the right-hand side, the snapshot used to define the reduced model is the solution of the nearby problem.

Nearby solution to construct the a posteriori reduced basis

In the second test case, we use as a snapshot the consecutive solutions obtained by solving a nearby problem, which differs from the reference one only by the position of the applied forces: homogenous Neumann boundary conditions are applied to any point P i (with i ∈ 1, n P ) such that P i ∈ P defined by:

P = {M = (x, y, z) | x ∈ [9, 11] and y ∈ [9, 11] and z = 11} ( 53
)
This snapshot is for instance a particular output of a parametric study meant to obtain the response of the system subjected to various load cases. One would expect to be able to reuse the information generated during this simulation to fasten the solution process of the next one. The difference between the snapshot solution and the reference solution that is looked for is illustrated in figure [START_REF] Xiao | Two-scale asymptotic homogenisation-based finite element analysis of composite materials[END_REF]. It should be noted that far away from the part of the structure where Neumann boundary conditions are applied, the solutions are qualitatively similar. Yet, they are very different within this process zone, which qualitatively justifies the local/global splitting proposed in section (4.3.1). The same effect appears in figure (6, which shows that the combination of the information from the two load cases can be significantly compressed when removing the process zone from the Only the solutions to the eighth load increment are represented in both cases. On the left-hand side, the snapshot is the reference solution which is looked for. On the right-hand side, the snapshot is the solution to a slightly different problem ("nearby problem"), where the position of the applied load has been modified. Figure [START_REF] Belytschko | Multiscale aggregating discontinuities: A method for circumventing loss of material stability[END_REF] compares the load/deflection curves obtained in the reference case, and when using successively the basic POD strategy and the local/global reduction technique. Obviously, the POD strategy looses its relevance very quickly, as the damage localises in the part of the structure where the damage would localise during the snapshot simulation. The resulting model is too stiff, and the maximum strength of the structure is overestimated. When using the local/global approach, with the parameters described previsouly for the fine analysis of the damaged areas, a much more accurate load/deflection curve is obtained. Yet the peak load is still overestimated by 5%. This remaining inaccuracy is due to the fact that the infomation from the snapshot is not correct in the reduced zone. Global corrections of the reduced basis are still necessary. This issue is discussed in the next section.

In figure (8, right) we analyse the effect of preconditioning the local iterations of the conjugate gradient by the projection technique. As in the previous test case, the convergence curves are observed at the third Newton Newton iteration of the fiftieth time step of the analysis. Again, the number of iteration to convergence is significantly decreased when using the augmentation. We would expect this precontionning technique to be less efficient in this particular case. Yet it is not the case, which is probably due to the fact that the global update described in section 4.3.2 is mechanically relevant.

Adaptive reduced basis

We now try to reduce the inaccuracy observed in figure [START_REF] Belytschko | Multiscale aggregating discontinuities: A method for circumventing loss of material stability[END_REF] by using global corrections of the reduced model. The adaptation algorithm used here was introduced in [START_REF] Kerfriden | Bridging proper orthogonal decomposition methods and augmented newton-krylov algorithms: an adaptive model order reduction for highly nonlinear mechanical problems[END_REF]. We just recall the basics of this particular strategy, and explain how it can be used to enhance the proposed local/global scheme.

"On-the-fly" global corrections of the reduced basis

The principle of these global corrections, is to update the reduced basis used to define reduced problem [START_REF] Belytschko | Multiscale aggregating discontinuities: A method for circumventing loss of material stability[END_REF]. It is, for now, not linked to the global/local technique described previously in this paper. The updates are done "on-the-fly", at any iteration i + 1 of the Newton algorithm used to solve [START_REF] Belytschko | Multiscale aggregating discontinuities: A method for circumventing loss of material stability[END_REF] at a given time step t n+1 of the analysis. If the reduced problem is sufficiently converged, and if the residual of the norm of the initial problem evaluated at iteration i is estimated too high, then the following linearized system is considered:

Ki T δ U = -R i (54)

Ki

T is an approximation of the tangent K i T , while R i is the residual of the initial system of equations ( 6) evaluated in U |tn + C i α i . Operator C i designates the reduced basis considered at iteration i of the Newton process.

The solution to this problem is searched for in two supplementary spaces by a projected Krylov algorithm:

δ

U = δ U C + δ U K where δ U C ∈ Im(C i ) δ U K ∈ Im(C i ) ⊥ Ki T = Ker(C i T Ki T ) ( 55
)
This problem is solved coarsely by a projected Krylov solver (10 -1 is the typical value of the stopping criterion), othogonally to the current reduced basis. The resulting solution, which is Ki T -orthogonal to the C i , is added to the reduced basis:

C i+1 = C i δ U K δ U K (56)
Newton prediction i + 1 can now be solved after updating residual R i R and reduced stiffness K i T,R with respect to this new basis vector.

Various features are added to this algorithm to efficiently sort the successively added basis vectors, but the main idea described above is sufficient to understand the following results.

Coupling with the local/global approach

Once the global norm of the residual of locally reduced problem ( 27) is sufficiently converged, and if the global residual is juged too high, a correction vector to the global reduced basis is calculated, following the algorithm given in the previous subsection. The residual and tangent operator of problem [START_REF] Nouy | A generalized spectral decomposition technique to solve a class of linear stochastic partial differential equations[END_REF] are updated and condensed tangent problem [START_REF] Ryckelynck | On the a priori model reduction: overview and recent developments[END_REF] is solved with the projected conjugate gradient algorithm. This algorithm converges in a very few iterations as the new basis vector, added to the augmentation space, is very close to the solution which is looked for. The slight mismatch is only due to:

• the approximation of the tangent used to provide the global correction vector.

• the coarse solution of (54).

• the approximation of the displacement in the reduced zone of the local/global approach.

Results

The test case described in section 4.4.2 is solved using successively the local/global technique on its own, the global correction algorithm on its own, and the local/global reduction technique coupled with global corrections. The target norm of the relative global residual norm is 10 -1 . Corrections are allowed only if the relative norm of the residual of locally reduced problem ( 27) is below 10 -3 . It is clear from the load-displacement curve in Figure 11 that the global corrections significantly improve the results of the local/global algorithm. [START_REF] Nguyen | Homogenization-based multiscale crack modelling: From micro-diffusive damage to macro-cracks[END_REF]. An illustration of the POD basis and its successive corrections is given in figure [START_REF] Verhoosel | Computational homogenization for adhesive and cohesive failure in quasi-brittle solids[END_REF].

More importantly, the number of relatively expensive global corrections required to achieve the target in terms of global residual norm is very low when using the local/global approach, as opposed to the one observed when using only global corrections. This effect is shown in figure [START_REF] Cazès | A thermodynamic method for the construction of a cohesive law from a nonlocal damage model[END_REF]. On can also notice that using only the global correction methods lead to an increase in the computational costs in the damage propagation phase of the analysis, which has been reported in [START_REF] Kerfriden | Bridging proper orthogonal decomposition methods and augmented newton-krylov algorithms: an adaptive model order reduction for highly nonlinear mechanical problems[END_REF]. With the local/global reduction approach, this effect disappears. The same trends are observed in [START_REF] Kerfriden | A three-scale domain decomposition method for the 3D analysis of debonding in laminates[END_REF], where a local/global algorithm was developed in a domain decomposition framework to avoid unnecessary computations far away from the process zones. This suggests that the result obtained here may be extended, and that the number of global corrections required to achieve a given level of accuracy does not depend on the local nonlinearity (i.e: independent on the damage state).

Conclusions and discussion

The work described in this paper focuses on decreasing the computational effort required in solving large scale damage and fracture problems, where small scale phenomena must be taken into account to accurately represent the global structural behaviour. Although most investigations on this issue rely on the extension of homogenization frameworks, we followed an alternative route, choosing projection-based model order reduction as a starting point. The proposed local/global reduction technique was demonstrated on a damaging 3D frame structure, but we expect that the algorithm will carry forward to the more general case of the initiation and propagation of cracks in structures.

The main conclusions of the paper are:

• The proposed local/global model order reduction approach can significantly improve the relevancy of using a global reduced basis to approximate the displacement field in the case of localised failure;

• The number of global corrections of the reduced model required to obtain a given level of accuracy is drastically reduced when the balance equations which exhibit the highest level of nonlinearity are excluded from the reduction. In most of the applications for which this methodology is intended, the additional computational costs are kept minimal for two reasons:

• The damaged zone where no reduction is performed is small compared to the size of the structure in realistic engineering problems;

• The local fine scale solution is efficiently obtained by Krylov corrections of a good initialisation provided by the restriction of the reduced model to the damaged zone.

Yet, some issues need to be addressed before global/local reduction algorithms can be efficiently and safely used to tackle engineering problem:

• Cost of the construction of the reduced systems. A system reduction is needed to ensure the computational efficiency of the method. It is now well known that projection-based model order reduction only reduces the size of the algebraic systems to solve, and that the construction of these systems is still expensive. Indeed, it involves the evaluation of the internal forces to construct the tangent stiffness and the residual of the reduced problem, which requires global updates of the constitutive behaviour, and integrations over the whole domain. System approximations permit to solve this issue. They basically consist in obtaining a cheap, but reliable approximation of the internal forces. In [START_REF] Kerfriden | Bridging proper orthogonal decomposition methods and augmented newton-krylov algorithms: an adaptive model order reduction for highly nonlinear mechanical problems[END_REF], we used the so-called hyperreduction [START_REF] Ryckelynck | A priori hyperreduction method: an adaptive approach[END_REF] (also introduced as the missing-point approximation by [START_REF] Astrid | Missing point estimation in models described by proper orthogonal decompo-sition[END_REF] in a different context). We found that adaptive reduced basis methods need to be carefully tailored when used within this framework. Indeed, one cannot rely on the fact that the residual of the full-scale problem is known at any iteration of the nonlinear solver anymore: the number of global evaluation of the internal forces must be kept minimal. In addition, we observed that the stability of the particular system approximation that we used is highly dependent on the choice of its parameters (reduced integration domain), especially when the damage mechanism localises. We foresee that excluding the most damaged zones from the reduced domain can significantly improve the stability of this scheme. Some numerical investigations are required to confirm this intuitive speculation.

• Quality control. Each step in the algorithm is a source of error: (i) choice of the set of degrees of freedom for which the solution is searched in a reduced space (ii) choice of the reduced basis (iii) choice of the accuracy for the local solution (iv) choice of the system reduction (v) choice of the time stepping parameters. These choices should be made within an encompassing framework to ensure that the error associated with each of them can be controlled in the same manner. Therefore, one needs to find a relevant error criterion and build a control process, capable of ensuring that the above choices introduce the same order of error in the solution so that the global error be not dominated by any of them in particular.

In terms of further applications, this method should be validated using more advanced models for brittle fracture, in particular nonlocal continuum damage or softening plasticity models, for which the distinction between local and global effects is not trivial. The extension to ductile failure and/or geometric nonlinearities should also be possible, the global corrections being expected to efficiently handle the "smooth" part of the nonlinearities, while the local/global scheme should treat separately the localisation of defects.

We believe that this local/global reduction framework is an important milestone to obtain an optimal reduction method applicable to the parallel simulation of failure in solids. Indeed, as alluded to in the introduction of this paper, several authors have proposed to perform a systematic model order reduction of the expensive local problems (so-called "problems by substructure") classically involved in domain decomposition methods [START_REF] Rixen | A dual craig-bampton method for dynamic substructuring[END_REF][START_REF] Markovic | Partitioning based reduced order modelling approach for transient analyses of large structures[END_REF][START_REF] Ladevèze | The latin multiscale computational method and the proper generalized decomposition[END_REF]. The application of such ideas to failure assessment has never been tackled so far, as it first requires an efficient reduction scheme for the local problems, such as the one proposed in this paper, in order to guarantee reasonable scalabilities and load balancing. of the Royal Academy of Engineering and of the Leverhulme Trust for Bordas' Senior Research Fellowship "Towards the next generation surgical simulators" (2009-2010) as well as the support of EPSRC under grant EP/G042705/1 Increased Reliability for Industrially Relevant Automatic Crack Growth Simulation with the eXtended Finite Element Method. The support of the British Council Exchange Programme is also gratefully acknowledged by the authors.

  Condensation of the linearized reduced problem on the master degrees of freedom Computation of the residual of the projected problem ("non-corrected" search direction) Coarse scale correction of the search direction by orthogonality with respect to the reduced space Orthogonalization and Minimisation Computation of the displacement of the slave degrees of freedom Computation of the displacement of the master degrees of freedom Convergence check YES NO A priori reduced basis F D New problem

Figure 1 :

 1 Figure 1: Graphical representation of the two-level solution algorithm.

4

  Application to the adaptive reduction of damage problems 4.1 Damageable lattice model We focus on a very simple damageable lattice structure, made of bars under traction or compression (figure (3)). Elementary bar indiced b occupies domain Ω b such that Ω = b∈ 1,n b Ω b , where n b is the number of bars. The direction of bar number b is denoted by n kl = P k P l / P k P l where P k and P l (k < l and (k, l) ∈ 1, n P 2

Figure 2 :

 2 Figure 2: Definition of the lattice problem.

Figure 3 :

 3 Figure 3: Reference solution obtained at the eighth load increment of the full-scale simulation (maximum value of prescribed forces before unstable damage propagation). Darker bars have a higher damage state.

9 ]

 9 and y ∈ [8, 10] and z =

Figure 4 :

 4 Figure 4: Final damage obained at the end of the full-scale simulation, performed in 30 load increments. The bars that are completely damaged are not represented.

4. 3 Figure 5 :

 35 Figure 5: Splitting obtained at two different stages of the simulation performed by the local/global reduction method (left: load increment 4, right: load increment 30). Dark bars are connected to at least one fully resolved node. As time increases, the degrees of freedom for which no reduction is performed localises in the region of failure.

Figure 6 :

 6 Figure6: Measure of the error in the snapshot representation as a function of the error of truncation of the singular value decomposition. In solid lines dark line, the set of 30 reference solutions are defined as the snapshot. In grey lines, a second snapshot operator is defined by concatenating the 30 reference solutions and the 30 solutions of the nearby problem. Both analyses are repeated when excluding the process zone (defined a priori by the set of points {M = (x, y, z) | x ∈[START_REF] Grassl | Damage-plastic model for concrete failure[END_REF][START_REF] Verhoosel | Computational homogenization for adhesive and cohesive failure in quasi-brittle solids[END_REF] and y ∈ [6, 13]})) from the snapshot definition (i.e.: the snapshot now is the restriction of the solutions to the part of the domain which does not undergo significant damage). The resulting error curves are plotted in dashed lines.

Figure 7 :

 7 Figure 7: Load/deflection curves obtained one using successively direct numerical simulation, simple proper orthogonal decomposition and local/global reduction technique. The snapshot used here is the solution of the direct numerical simulation.

Figure 8 :

 8 Figure 8: Convergence curves of the augmented conjugate gradient used to solve the local problem, at the third Newton iteration of load increment 15. On the left-hand side, the reference solution is used as a snapshot. On the right-hand side, the snapshot used to define the reduced model is the solution of the nearby problem.

Figure 9 :

 9 Figure 9: Two set of reference solutions obtained by direct numerical simulations and used as snapshot vectors in our examples.Only the solutions to the eighth load increment are represented in both cases. On the left-hand side, the snapshot is the reference solution which is looked for. On the right-hand side, the snapshot is the solution to a slightly different problem ("nearby problem"), where the position of the applied load has been modified.

Figure 10 :

 10 Figure 10: Load/deflection curves obtained one using successively direct numerical simulation, simple proper orthogonal decomposition and local/global reduction technique. The snapshot used here is the solution of the direct numerical simulation to the nearby problem.

  SVD analysis.

Figure 12 :

 12 Figure 12: Number of global corrections performed when using the adaptive model order reduction on the one hand, and when using adapative local/global model order reduction on the other hand.

15 A 6

 156 posteriori POD basis (decreasing order of the associated singular value) Global corrections of the reduced basis (ordered as a function of time)

Figure 13 :

 13 Figure 13: A posteriori POD basis (3 vectors) and successive global corrections (4 vectors computed respectively at time steps 1, 3, 10 and 18). Two different views of each basis vector are represented for better qualitative understanding.

Acknowledgements

Pierre Kerfriden would like to thank Pierre Gosselet (LMT/ENS Cachan, France) for his precious advices. Stéphane Bordas and Pierre Kerfriden would like to acknowledge the financial support