P Feiler

L Wrage

J Hansson

System Architecture Virtual Integration: A Case Study

Keywords: virtual integration, embedded software system, safety-critical, validation & verification

Aerospace industry is experiencing exponential growth in the size and complexity of onboard software. It is also seeing a significant increase in errors and rework of that software. All of those factors contribute to greater cost; the current development process is reaching the limit of affordability of building safe aircraft. An international consortium of aerospace companies with government participation has initiated the System Architecture Virtual Integration (SAVI) program, whose goal is to achieve an affordable solution through a paradigm shift of-integrate then build. A key concept of this paradigm shift is an architecturecentric approach to analysis of virtually integrated system models with respect to multiple operational quality attributes such as performance, safety, and reliability. By doing so early and throughout the life cycle at different levels of fidelity, system-level faults are discovered earlier in the life cycle-reducing risk, cost, and development time. The first phase of this program demonstrated the feasibility of this new development process through a proof of concept demonstration and a return on investment analysis, which are the topics of this paper.

Introduction

Aerospace industry is experiencing exponential growth in size, complexity, errors, rework and cost of their onboard software. The current development process is reaching the limit of affordability of building safe aircraft. The size in terms of source lines of code (SLOC) has doubled every four years. The cost 27M SLOC of software is estimated to exceed $10B (see Figure 1). The Aerospace Vehicle Systems Institute (AVSI) has launched an international, industry-wide multi-phase program called System Architecture Virtual Integration (SAVI), whose members are acknowledged in Section 7. The objective of SAVI is to reduce cost/cycle-time and risk by using early and repeated "virtual" integration and analysis. SAVI builds on model-based engineering as practiced by the system and software engineering communities, where models of different aspects of a system are developed and analyzed (see Figure 2). However, industrial practice has shown that such independently developed models tend to result in multiple versions of the "truth." SAVI improves current practice by ensuring that analytical models are consistent with each other and the evolving architecture throughout the lifecycle. To achieve these goals, the SAVI paradigm necessitates

• an architecture-centric multi-aspect reference model approach as the single source of "truth",

• a component-based framework in support of model-based and proof-based engineering,

• a model bus for consistent model interchange between repositories and tools, and

• an architecture-centric acquisition and development process throughout the system life cycle that is supported by industrial standards and tool infrastructure.

Key Concepts of SAVI

SAVI is driven by three key concepts: a single source architecture reference model, a model repository and model bus, and a model-based development process that encompasses both system engineering and embedded software system engineering.

Single Source Architecture Reference Model

A key concept of virtual integration is the use of an annotated architecture model with well-defined semantics as the single source for architecture analysis, as illustrated in Figure 3. The annotations provide relevant analysis-specific information (e.g., fault rates or security properties) that allow autogeneration of an analytical models.

Any changes to the architecture throughout the life cycle are reflected in this model and automatically propagate to all dimensions of analysis. For example, substitution of a faster processor to accommodate a high workload is reflected in schedulability analysis, impacts end-to-end response time, and requires re-evaluation power consumption against capacity as well as its impact on the mass of the system.

Model Repository and Model Bus

A second key concept is the use of a model repository and model bus, illustrated in Figure 4. The model repository contains the annotated architecture reference model, as well as detailed models that are refinements of architecture components. For example, details of physical system components can be modeled with Modelica, control system components with Simulink, and discrete application behavior with UML statecharts. This architecture reference model may reside in a single model repository or it may be distributed across different model repositories, e.g., those of different suppliers and the system integrator. The model bus is a data interchange mechanism a mechanism that operates with a standardized model representation. It supports interchange between model repositories and translation into -as much as possible standardized -representations acceptable to different analysis and generation tools. For example, it supports the interchange of annotated architecture models in a standardized XML format. Similarly, transformation specifications provide the translation from an annotated architecture model to specific analytical model formats, such as those for timing models, fault trees, or security models.

Development Process

The SAVI development process must support collaboration between

• system engineers-whose primary focus is the architecture of the physical system

• embedded software system engineers-whose focus is the interaction between the physical system architecture, the computer platform architecture, and the task and communication architecture of the embedded application

• other engineers-whose focus is the architecture and detailed design of physical components, computer hardware components, and software components. Thus, the system is expressed in a combination of modeling notations whose models are mapped into a single underlying reference model in the model repository with minimal redundancy of information that must be maintained manually. Figure 5 illustrates the relationships between such models. To support the development process, the repository also includes requirements, test data, and analysis results. Standardized interchange formats such as the Requirements Interchange Format (RIF) and AP233 for system engineering artifacts are emerging. Therefore, SAVI will define the data structures needed in the model repository for information storage and analysis, and data transformations needed for data interchange and to leverage ongoing efforts in standard organizations.

Standardized interchange representations for these analytical models facilitate tools chains (integration of multiple tools) to enable new and more effective integration checks and analyses by leveraging existing tools and minimizing transformations into tool-specific and proprietary representations. Examples of such emerging interchange formats include ISO/IEC 15909 for Petri nets to interface with Petri net analysis tools [START_REF] Iso | Software and System Engineering -Highlevel Petri nets -Part 2: Transfer Format[END_REF], and FIACRE for statebased behavior specification to interface with different model checking tools [START_REF] Berthomieu | Fiacre: an Intermediate Language for Model Verification in the TOPCASED Environment[END_REF]. SAVI virtual integration activities augment traditional design reviews by a model-based approach. Subsystem requirements are recorded in an initial system model during request for proposals and made available to suppliers. Suppliers respond to the proposals by including subsystem models with interface specifications that include non-functional properties. This allows the system integrator

• to validate supplier model compatibility, interface consistency, and initial resource allocations during proposal evaluation,

• to assure compatibility of functionality with interface specifications including non-functional properties during preliminary design integration, and

• to verify system-level non-functional properties such as performance or safety during critical design integration. Early and continuous virtual model integration based on standardized representations insures that • errors are detected as early as possible with minimal leakage to later phases,

• models with well-defined semantics facilitate autoanalysis and generation to identify and eliminate inconsistencies,

• automated compatibility analyses at the architecture level scale easily,

• industrial investment in tools is leveraged through well-defined interchange formats.

Proof-of-Concept Phase

To establish cost-effective management and limit risks of the SAVI program, the Proof-Of-Concept (POC) project has been executed as the first of multiple phases with the following goals:

• Document the main differences between a conventional acquisition process and the projected SAVI acquisition process and identify potential benefits of the SAVI acquisition process.

• Evaluate the feasibility and scalability of the multiaspect model repository and model bus concepts central to the SAVI project.

• Assess the cost, risk, and benefits of the SAVI approach through a return on investment (ROI) study and development of a SAVI development roadmap.

Proof of Concept Demonstration Requirements

The SAVI POC team established a prioritized set of requirements, which are summarized in Table 1.

Note that emphasis was placed on validation early in the development process. Architecture models were developed in AADL because it specifically supports modeling safetycritical, software-reliant systems and is extensible to support analysis of operational quality attributes along multiple dimensions [START_REF]Architecture Analysis & Design Language (AADL): SAE International Standards document AS5506A[END_REF].

Figure 6 shows a drawing of the aircraft system provided to the POC team. It shows major physical subsystems, some providing aircraft capability, such as navigation or landing gear, and others providing physical resources to the subsystems, such as the electrical power, hydraulics, and fuel. We have elaborated the flight guidance system (FGS) into a Tier 2 model representing the distributed computer platform (physical view) and the embedded application subsystems (logical view) of the IMA subsystem. This elaboration is not a separate model, but a refinement of the FGS system component using the AADL extends mechanism.

Because of this refinement, we can now specify a Tier 1 variant and a Tier 2 variant of the aircraft model and instantiate both for analysis from a single source.

Figure 8: IMA Computer Platform

Figure 8 shows a portion of the physical view-that is, devices to represent sensors and actuators to the physical system, buses to represent networks such as ARINC429, and systems to represent processing units and communication units. The symmetry reflects the dual redundant nature of the IMA platform.

Figure 9 shows a portion of the logical view as a collection of embedded application subsystems. We have used AADL port groups and connections to model interaction between subsystems. Port groups represent a collection of individual port connections, which suppliers later elaborate through port group types.

Figure 9: IMA Embedded Application Subsystems

We have also included two end-to-end flows in order to analyze the stick-to-surface response time when operating in direct mode (with a maximum allowed latency of 150 milliseconds) and with flight guidance and autopilot involved (with maximum allowed latency of 25 milliseconds). IMA computer resources (MIPS for processors, memory size for RAM and ROM, and bandwidth for networks), in addition to weight and electrical power are specified as properties of processor and bus components. Similarly, computer resource budgets are assigned to the application subsystems and endto-end latency requirements to the flows. In our demonstration, seven of the IMA subsystems are contracted out to suppliers. The suppliers first refine the subsystem specification with interface details such as

• detailed communication properties,
• properties of communicated data, and

• mapping into protocols such as ARINC429.

Proof of Concept Scenarios

The requirements for the POC demonstration include several development use scenarios:

• Aircraft system modeling

• Modeling of the IMA system as embedded software system

• Subcontracting support between airframer and suppliers

• Subsystem development by supplier

• Virtual integration testing by airframer 4.1 Tier 1 Aircraft System Modeling Scenario

The first scenario illustrates analysis early in the life cycle based on the Tier 1 model, whose results are revalidated throughout the life cycle as the fidelity of the model increases.

Analysis Reporting

Since AADL is a strongly typed language, valid models guarantee a certain level of system architecture consistency without specialized analysis. For example, the landing gear system type has been specified to require access to the power system and the hydraulic system. The AADL semantic checker ensures that the correct bus is connected to the bus access feature and that all features that require a connection are connected. Physical characteristics relevant to system engineers, such as mass and electrical power, are represented by AADL properties. Information regarding the mass of a system is typically kept in a spreadsheet that must be manually updated and analyzed from time to time. Instead, we associate mass information with the AADL model and drive the analysis from the model. By doing this, we can analyze the mass of the aircraft for the Tier 1 model and then revisit the analysis with more details about the mass at the Tier 2 level. The analysis examines the net weight, gross weight, and weight limit of the physical system components and connections. In the case of the Tier 1 model, we can run an analysis that adds up the gross values of the Tier 1 elements and compares the total against the limit for the aircraft. We have implemented the analysis as a model bus transformation that translates the relevant data into the CSV (comma separated values) file format for import into a spreadsheet. The electrical power information is recorded in the AADL model as PowerCapacity, PowerSupply, and PowerBudget properties. The analysis, in this case implemented as an OSATE plug-in, compares the power supplied to the power system (an instance of bus type PowerSystem) from both the engine and the auxiliary power unit against the power system capacity as well as the power budgets of components drawing power from the power system. This analysis can be revisited when the Tier 2 model of the IMA is available to look at the power distribution of the IMA power subsystem to the computer hardware components and compare the demand against the power budget assigned in the Tier 1 model.

Use Scenario

In our use scenario, the analysis reports that 1. The power supplied by the engine and the auxiliary unit exceeds the capacity of the power system.

2. The power budgets also exceed the capacity but are less than the power supply. As a remedy, we choose a higher capacity variant of the power system (from the component specification library, an AADL package) and rerun the electrical power analysis. We must also rerun the mass analysis from the same model to ensure that the change in the power system has not exceeded any weight limits. The AADL model of the system drives various highlevel quantitative system analyses. From the same model, we can perform analyses for hydraulic pressure, fuel flow, and airflow once the component specifications are annotated with the relevant properties. The models can be made more realistic by adding by refining the AADL component specifications or by associating detailed physical models using a specialized notation, e.g., Modelica, with a component.

Tier 2 Embedded Software System Modeling Scenario

The Tier 2 model refines the IMA part of the system into a networked computing platform and an interacting set of application subsystems and blackbox subsystems, which will get subcontracted to suppliers. However, before doing so, the airframer will analyze the Tier 2 model to revalidate the mass and electrical power results, by taking into account the more detailed architecture specification, and to validate properties specific to the elaborated IMA subsystem (i.e., computer resource usage and endto-end flow response time). The analysis shows that the IMA power subsystem draws less power from the main power system and that the power consumption by the computer platform is at 60% of the locally available power. Therefore, we could consider reversing the earlier main power system upgrade.

Analysis Reporting

The computer resource analysis comes in two variants:

1. Budget totals against capacity totals 2. Budgets of deployed application components against the target resource once deployment decisions have been recorded This computer resource analysis was demonstrated with an OSATE plug-in and can also be supported through a spreadsheet interface similar to the mass analysis.

Use Scenario

In our use scenario, the analysis reports that the MIPS budget totals exceed the total capacity of all processors. It also indicates that only a subset of the application components has a MIPS budget. The memory budget totals reflect 70% of the components with budgets and are well within memory capacity. In our use scenario, we reduce the budgets and expect suppliers to meet them with their Tier 3 models. Such a decision can be justified based on historical data, if available. At this point in the life cycle, or at a later stage, the system architect may make a first attempt at an allocation of major application subsystems to hardware. A variant of the computer resource analysis will consider the deployment in its results. Sampling jitter and changes in latency due to implementation decisions regarding the runtime architecture can affect the stability of control systems [START_REF] Cervin | Control Loop Timing Analysis Using TrueTime and Jitterbug[END_REF]. The end-to-end latency analysis [START_REF] Feiler | Impact of Runtime Architectures on Control System Stability[END_REF] at the Tier 2 level takes into account

• processing latency in the stick and surface (represented by a latency property on the flow specification in the respective AADL device type)

• communication latency associated with the connections involved in the flow

• processing and sampling latency of IMA subsystems involved in one of the two end-to-end flows The latency analysis [START_REF] Feiler | Flow Latency Analysis with the Architecture Analysis & Design Language (AADL)[END_REF] calculates the minimum worst-case latency for the two flows (a lower bound that can only increase as the model is refined) and reports that in direct mode the response time requirement is met (121ms versus a 150ms requirement), while the IMA-based response time is almost twice the original requirement (46 ms versus the required 25 ms). The main contributor to the increased latency value is the sampling latency of the partition. We could reduce the latency by not sampling (moving to a data-driven architecture), double the partition execution rate (doubling our processor resource requirements), or renegotiate the response time requirement as an inherent property of the chosen runtime architecture. In our scenario, we pursue the latter option.

Airframer-Supplier Subcontracting Scenario

In support of subcontracting, we have organized the AADL model into a number of separate AADL packages that are version controlled through a model repository. For the demonstration, the AVSI organization hosted this model repository on a Subversion server with POC demo team members playing the airframer and the supplier roles located at two sites in the U. S. and two sites in Europe.

Repository Organization

The repository has different access-controllable public and internal areas for the airframer and the subcontractors, as shown in Figure 11. The standardized AADL XMI representation [START_REF]Architecture Analysis & Design Language (AADL) Annex Volume 1: Annex C: AADL Meta model & XML Interchange Format Annex[END_REF] enables inter-repository model interchange.

Use Scenario

In the scenario, as part of a request for proposals (RFP), the airframer makes available the AADL specification of the desired subsystems with a possibly partial interface specification (system types) including resource budget properties as well as expected latency requirement on flow specifications through the subsystem, and an interface control document in the form of port group type specifications.

We have done so for seven subsystems to be contracted out. The suppliers respond, in this scenario, with a completed subsystem specification including details about the exchanged data and its mapping into the ARINC429 protocol. The airframer verifies overall consistency by virtually integrating the suppliers' AADL subsystem specifications as a variant of the Tier 2 model and performing functional integrity checks. Figure 12 shows reported inconsistencies that are traceable to the model. The suppliers refine their subsystem AADL models to model their architecture and reflect implementation decisions. This example included three suppliers expanding their subsystems. For the air data computer, we have included UML diagrams, Ada code, a test harness and automatic build scripts, as shown in Figure 13. The supplier allocates threads of the application task model to computer hardware and performs scheduling analysis. The scheduling protocol property of the processor determines which scheduling analysis algorithm is used. Analysis results indicate that the subsystem is schedulable on its internal hardware with 45% utilization. Scheduling analysis uses the worst-case execution time of threads. This figure may be an estimate early in the development that is scaled for processors of different speed. Once the source code exists and has been benchmarked on different processors the benchmark figures replace the estimates in the model resulting in higher fidelity results.

Virtual System Integration Testing Scenario

At various stages of the development, each supplier delivers an AADL model of their subsystem architecture to the airframer. These updated models identify properties of the subsystems pertinent to integration, but do not necessarily include detailed design descriptions. The airframer virtually integrates them into an aircraft model refined down to Tier 3, while the supplier independently validates that the model properties of the architectural components reflect the detailed design.

Analysis Reporting

The models include properties which can be used for traceability to requirements. The airframer queries the model repository to determine which subcontracted subsystems are involved in satisfying certain requirements. This allows the airframer to focus on virtually integrating and analyzing the system with respect to requirements of greatest concern.

Use Scenario

The airframer revisits the mass and power analysis to include the Tier 3 details in the results. Similarly, computer resource analysis aggregates periods and execution times of the subsystem task models, compares them against the assigned budgets and rolls them up for a comparison against capacity. The analysis results show that many subsystems stay within the reduced budgets from the Tier 2 analysis. As part of the scenario, the airframer revisits the end-to-end latency analysis, now taking into account any latency or latency jitter contributed by the subsystem task models that exceed the expected latency as recorded in the flow specification of the subsystems in the Tier 2 model. When task models are first delivered by the supplier, the airframer reruns the end-to-end latency analysis, which now takes the task model into account in its latency calculation. The airframer discovers that the minimum worst-case end-to-end latency for the IMA mode has increased considerably to 185ms on a synchronous hardware platform and to 196ms for processors operating on independent clocks. Examination of the detailed analysis report reveals a low-rate thread in one of the subsystems contributing a sampling latency of 100ms. Without virtual integration, such a problem would only be discovered during system integration test. The airframer allocates the application tasks from the different supplier task models to the various processors and performs scheduling analysis. For a given deployment configuration of a three-processor system, the analysis reports that all deadlines are met with processor utilizations of 54%, 55%, and 75%. The airframer can validate the analysis results by applying a different scheduling analysis tool. In both cases, a model transformation is performed to generate a timing model in the representation acceptable to the analysis tools from the AADL model. A resource allocation tool [START_REF] Deniz | On Resource Allocation in Architectural Models[END_REF] provides an option to explore alternative deployment configurations, showing that the system would be schedulable with 97% utilization on two processors, and suggesting a three processor allocation that better balances the task load, and supporting a quick what-if analysis of a four processor system to reduce the average processor utilization to a target of 50%. Finally, the airframer performs network bandwidth analysis on the aircraft model with Tier 3 detail and a specific deployment configuration. This analysis identifies all application task connections that are routed over a particular network and determines the data volume from the size of the data communicated through ports and their transfer rate. This data volume is then compared against the bandwidth budget assigned at the Tier 2 level and against the capacity of each the network.

Return on Investment Analysis

In this section we summarize the result of an ROI analysis [9] that focuses on cost avoidance due to early detection of design defects through virtual integration. Rework cost is dominantly driven by the cost of managing defects injected in the requirements and design phases but which are detected late in the system development, most often in the integration and test phases, where the cost to repair the defect is one to two orders of magnitude higher. Virtual integration reduces and prevents this down-stream flow. Current development processes allow 70% of faults to be introduced early in the life cycle, while 80% of them are not caught until integration test or later with a repair cost of 16x or higher. Figure 14 shows percentages for fault introduction, discovery, and cost factor of repair [START_REF] Rti | The Economic Impacts of Inadequate Infrastructure for Software Testing[END_REF][START_REF] Galin | Software Quality Assurance: From Theory to Implementation[END_REF][START_REF] Boehm | Software Engineering Economics[END_REF]. If we can use the SAVI approach of architecture-centric virtual integration and analysis to discover a reasonable percentage of system-level faults earlier in the process, we can expect cost savings larger than the additional investment in modeling and analysis. We estimate the size and complexity of an avionics software system based on data from previously built aircraft. We use two scenarios that reflect the current aircraft software (27 and 30 MSLOC), and a synthetic software system (60 MSLOC), which would not be affordable using current methods. Cost was estimated using the COCOMO II model [START_REF]COCOMO II[END_REF]. Based on these "as-is" costs, we compute cost reduction for the "to-be" scenario, i.e., after introduction of the SAVI methodology. To estimate avoided rework cost, we calculate an estimate for the cost to manage defects in the various phases. This is based on published empirical data from case studies and on observations providing data on the number of defects that would be introduced in a system of a defined size and complexity, and the nominal cost of removing defects. Table 2 shows the resulting cost to remove a defect of a given type relative to the total cost of defect removal. For example, requirements defects account for 79% of rework cost and 62% of rework cost occur during integration. Two scenarios of rework cost were considered: 30% and 50% of the total development cost. Similarly, we considered two scenarios for defect removal efficiency due to SAVI: 33% and 66% as applied to requirement defects. This is based on using Hayes' fault taxonomy [42] and estimating removal efficiency for the per fault category. Finally, ROI and NPV were calculated in light of an estimated cost of $86M to implement SAVI. In summary, the ROI analysis showed the following: The nominal cost reduction for a 27MSLOC system is $2,391M (out of an estimated $9,176M) occurring at 70% reuse, 50% rework, and with a removal efficiency of 66%. The arithmetic and logarithmic ROIs are 26.8 and 3.33, with an NPV of $1,076M. Studying the most conservative of the scenarios, the smallest cost avoidance is $717M, occurring at 70% reuse, 30% rework cost, with a removal efficiency of 33% for a 27MSLOC system. The arithmetic and logarithmic ROIs are 7.3 and 2.12, with an NPV of $263M.

Conclusion

In this paper, we have focused on reporting the experiences derived from AVSI's SAVI project, utilizing an architectural reference model to achieve virtual integration. Virtual integration demonstrated early and repeated quantitative analysis at various levels of fidelity, validation of architecture consistency across subcontracted subsystem interfaces and independent protocol mappings, and discovery of intricate operational system-level faults due to design problems in the runtime architecture. The resulting early discovery of system-level faults reduces risk, lowers system life cycle costs, and improves quality. The SAVI POC demonstrated how multi-tier modeling and analysis across levels, coverage of system engineering and embedded software system analysis, propagation of changes across multiple analysis dimensions, maintenance of multiple model representations in a model repository, autogeneration of analytical models via model bus, interfacing of multiple tools to perform the same analysis, as well as distributed team development via repository to support airframe manufacturer and supplier interaction are possible through the SAVI concepts. In addition, ROI analysis has shown that SAVI can lead to significant cost savings by avoiding defects in early development phases (requirements and design), by reducing corresponding costly rework in the later phases.

Figure 1 :Figure 2 :

 12 Figure 1: Estimated Onboard SLOC Growth and associated Cost

Figure 4 :Figure 3 :

 43 Figure 4: Model Repository and Model Bus

Figure 5 :

 5 Figure 5: Collaborative Engineering

Figure 6 :

 6 Figure 6: Aircraft System DrawingFigure7shows a portion of the corresponding Tier 1 AADL model. In this model, we have represented the physical subsystems as AADL systems that can later be refined and the physical resources as AADL buses. Each aircraft subsystem is represented by a separate AADL system type, and each type of resource by a bus type. Bus access connections represent the physical connection between subsystems and their resources. The model elements include properties about physical characteristics (e.g., mass) of the subsystem. In addition, each bus type has a resource capacity property, and the bus access features (connection points) have resource supply or budget properties, e.g., the engine contributing electrical power and the cockpit drawing electrical power.

Figure 7 :

 7 Figure 7: AADL Model of Aircraft System

Figure 10 :

 10 Figure 10: Spreadsheet Based Mass Analysis

Figure 11 :

 11 Figure 11: Distributed Model Repository

Figure 12 :

 12 Figure 12: Functional Integrity Checking 4.4 Tier 3 Supplier Subsystem Development Scenario

Figure 13 :

 13 Figure 13: Supplier Subsystem Model & Analysis

Figure 14 :

 14 Figure 14: Benefits of Early Fault Discovery

Table 1 :

 1 Prioritized POC Requirements

	# Requirement	Category
	1 Establish Model Bus infrastructure	Process
	2 Establish Model Repository	Process
	Infrastructure
	3 Inform return on investment (ROI)	Process
	estimates through POC performance &
	results
	4 Analyses be conducted across the system Analysis
	5 Two or more analyses must be	Analysis
	conducted
	6 Analyses be conducted at multiple levels	Analysis
	of abstraction
	7 Analyses must validate system model	Analysis
	consistency at multiple levels of
	abstraction
	8 Analyses must be conducted at the	Analysis
	highest system level abstraction
	9 Model infrastructure must contain	Model
	multiple model representations
	10 Model infrastructure must contain	Model
	multiple communicating components

Table 2 :

 2 Relative Defect Removal Cost

		Phase in which defect is removed
	Defect	Require-
	Type	ments DesignCode Test Integration Sum
	Require-
	ments	0.03% 0.21% 1.87% 28.11%	48.73% 78.97%
	Design	0.04% 0.37% 5.79%	9.75% 15.96%
	Code	0.19% 1.28%	3.10% 4.57%
	Test	0.17%	0.26% 0.43%
	Integra-
	tion		0.09% 0.09%
	Sum	0.03% 0.26% 2.44% 35.36%	61.92% 100.00%
	Increasing defect detection and removal efficiency would positively influence rework cost. The cost avoided is computed as

Several areas were not fully explored and will be addressed by SAVI phase 2:

• support for multiple architecture modeling notations (SAE AADL & OMG SysML), and integration with mechanical system modeling (e.g., Modelica) and control system modeling (Simulink),

• validation focus on a system-level non-functional property such as safety, reliability, or security,

• end-to-end validation of systems from requirements to models and system implementation,

• scalability of the model repository and commercial tool support.

Acknowledgements

System Architecture Virtual Integration (SAVI) is an industry initiative by a number of aerospace companies and government organizations to improve the engineering practice for software-reliant aircraft systems under the umbrella of the Aerospace Vehicle Systems Institute (AVSI). The proof-ofconcept phase of this initiative was carried out by representatives of the following member companies and organizations: AVSI, Boeing, Airbus, Lockheed Martin, Rockwell Collins, BAE Systems, GE Aviation, U.S. Army, and Federal Aviation Administration (FAA). The Software Engineering Institute (SEI) contributions were funded by SAVI members through AVSI. In particular the authors would like to acknowledge the contributions by the other members of the proof-of-concept (POC) demonstration project, which is the focus of this case study report, to the definition of a to-be process and the development of a return on investment (ROI) model. Keith Appleby (BAE Systems), John Glenski (Rockwell Collins), Jean-Jacques Toumazet (Airbus), Joe Shultz (GE Research), and Dave Redman (AVSI) actively participated in creating the aircraft model for the demonstration and the Microsoft Excel spreadsheet version of the mass analysis. We would also like to thank the other members of the full SAVI team who participated in shaping the requirements for the POC demonstration and use cases.


Mainak Bhattacharyya

Amelie Fau

Udo Nackenhorst

David Néron

Pierre Ladevèze

A LATIN-based model reduction approach for the simulation of cycling damage

Keywords: LATIN method, Proper Generalised Decomposition, Non-linear solid mechanics, Damage, Reduced order model

The objective of this article is to introduce a new method including model order reduction for the life prediction of structures subjected to cycling damage. Contrary to classical incremental schemes for damage computation, a non-incremental technique, the LATIN method, is used herein as a solution framework. This approach allows to introduce a PGD model reduction technique which leads to a drastic reduction of the computational cost. The proposed framework is exemplified for structures subjected to cyclic loading, where damage is considered to be isotropic and micro-defect closure effects are taken into account. A difficulty herein for the use of the LATIN method comes from the state laws which can not be transformed into linear relations through an internal variable transformation. A specific treatment of this issue is introduced in this work.

Introduction

In engineering problems, the failure of structures is often governed by the creation and growth of micro-voids and micro-cracks. These micro-level discontinuities can coalesce and grow into macro-cracks which may propagate through structures resulting in failure. They are quantified in continuum scale as a homogenised quantity described by an internal variable. Damage, in its mechanical sense, is the creation and growth of micro-voids or micro-cracks. The representation of damage as a scalar variable was first introduced by Kachanov [START_REF] Kachanov | Introduction to continuum damage mechanics, mechanics of elastic stability[END_REF]. Many works have been followed since then, notably by Lemaitre and Desmorat [START_REF] Lemaitre | Coupled elasto-plasticity and damage constitutive equations[END_REF][START_REF] Lemaitre | Engineering damage mechanics: ductile, creep, fatigue and brittle failures[END_REF]. For ductile materials, isotropic damage can be physically interpreted as a volume density of micro-voids, and also as reduction in the stiffness of the material [START_REF] Lemaitre | Engineering damage mechanics: ductile, creep, fatigue and brittle failures[END_REF]. Mostly, the isotropic damage variable is defined between 0 and 1 [START_REF] Kachanov | Introduction to continuum damage mechanics, mechanics of elastic stability[END_REF][START_REF] Lemaitre | Coupled elasto-plasticity and damage constitutive equations[END_REF][START_REF] Lemaitre | Engineering damage mechanics: ductile, creep, fatigue and brittle failures[END_REF][START_REF] Murakami | Continuum damage mechanics: a continuum mechanics approach to the analysis of damage and fracture, solid mechanics and its applications[END_REF], where 0 represents a virgin material and 1 represents complete failure. It is, however, also possible to introduce a damage variable B Mainak Bhattacharyya mainak.bhattacharyya@ibnm.uni-hannover.de 1 IBNM, Leibniz Universität Hannover, Appelstraße 9a, 30167 Hannover, Germany 2 LMT, ENS Paris-Saclay, CNRS, Université Paris-Saclay, 61 avenue du Président Wilson, 94235 Cachan, France between 0 (undamaged) and ∞ (completely damaged) [START_REF] Bellenger | Plastic and viscoplastic damage models with numerical treatment for metal forming processes[END_REF]. The phenomenon of low cycle fatigue (LCF) is governed by macro-plasticity (or visco-plasticity), where the dissipative energy for damage exists at the macro-scale [33]. For high cycle fatigue (HCF), the structure is macroscopically elastic, however plasticity and damage exist at the microscale. In the presence of cyclic loading, damage during tension is much higher than during compression because of the closure of the micro-voids and the material recovers its stiffness, partially or completely [START_REF] Lemaitre | Engineering damage mechanics: ductile, creep, fatigue and brittle failures[END_REF]. For solving this type of problems numerically, most of the developments in the field of continuum damage mechanics are based on a classical time incremental framework [START_REF] Lemaitre | Engineering damage mechanics: ductile, creep, fatigue and brittle failures[END_REF]33]. Numerical challenges rely on the accuracy of the required iterative algorithm as well as on its stability when damage increases [START_REF] De Souza Neto | Computational methods for plasticity[END_REF]o ro nt h e efficient computation of the response under cyclic loading [START_REF] Cojocaru | A simple numerical method of cycle jumps for cyclically loaded structures[END_REF][START_REF] Guennouni | Sur une méthode de calcul de structures soumises à des chargements cycliques : l'homogénéisation en temps[END_REF][START_REF] Maitournam | Détermination de la réponse asymptotique d'une structure anélastique sous chargement thermodynamique cyclique[END_REF][START_REF] Stolz | Optimal control approach in nonlinear mechanics[END_REF].

For the last decades, model order reduction has been appealing for many mechanical problem to tackle large size problems, real-time computations or parametric studies [START_REF] Chatterjee | An introduction to the proper orthogonal decomposition[END_REF][START_REF] Galvis | Spectral multiscale finite element for nonlinear flows in highly heterogeneous media: a reduced basis approach[END_REF]. For that purpose, the idea of model reduction is to seek the solution of a given problem in a reduced-order basis, whose dimension is much smaller than the size of the original high-dimensional model. This idea takes advantage of the redundancy of information that usually exists when describing the solution. Model reduction techniques based on the Proper Orthogonal Decomposition (POD) [START_REF] Atwell | Proper orthogonal decomposition for reduced basis feedback controllers for parabolic equations[END_REF][START_REF] Carlberg | Efficient nonlinear model reduction via a least-squares petrov-galerkin projection and compressive tensor approximations[END_REF][START_REF] Chatterjee | An introduction to the proper orthogonal decomposition[END_REF][START_REF] Kunish | Pod-based feedback control of the burgers equation by solving the evolutionary hjb equation[END_REF], involve a learning phase which consists in solving the full-order problem at some particular time instants and/or parameter values arbitrarily chosen. These solutions, called snapshots, are used to build reduced-order basis, by means of truncated POD, and a reduced-order model is then generated (e.g. by Galerkin projection) to solve the problem for the entire time and/or parameter domain. Since the 80's, the Reduced-Basis (RB) approach has been developed and consists of a greedy algorithm which selects the most relevant calculations to be performed on the parametric space in order to enrich the reduced-order basis [START_REF] Grepl | Efficient reducedbasis treatment of nonaffine and nonlinear partial differential equations[END_REF][START_REF] Maday | The reduced-basis element method: application to a thermal fin problem[END_REF][START_REF] Nguyen | Certified real-time solution of parametrized partial differential equations, handbook of materials modeling[END_REF][START_REF] Prud'homme | Reliable real-time solution of parametrized partial differential equations: reduced-basis output bound methods[END_REF][START_REF] Rozza | Reduced-basis methods for elliptic equations in sub-domains with a posteriori error bounds and adaptivity[END_REF][START_REF] Rozza | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations[END_REF][START_REF] Veroy | Certified real-time solution of the parametrized steady incompressible navier-stokes equations: rigorous reduced-basis a posteriori error bounds[END_REF]. Another appealing family of model reductions which has received a growing interest during the last decade is based on the Proper Generalised Decomposition (PGD). Roughly, PGD also consists in seeking the solution of a problem in a relevant reduced-order basis but this one is generated automatically and on-the-fly by a greedy algorithm, simultaneously with the successive approximations of the solution. The interested reader is referred to [START_REF] Ladevèze | Nonlinear computational structural mechanics, mechanical engineering series[END_REF][START_REF] Ladevèze | On reduced models in nonlinear solid mechanics[END_REF] for structural non-linear problems and the works of Chinesta and co-authors to solve parameterised problems [START_REF] Aguado | Real-time monitoring of thermal processes by reduced-order modeling[END_REF][START_REF] Ammar | On the space-time separated representation of integral linear viscoelastic models[END_REF][START_REF] González | Computational vademecums for the real-time simulation of haptic collision between nonlinear solids[END_REF]]. An extensive review of the literature on this method can be found in [START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF].

Computation of continuum damage problems may lead to solution bifurcation or to strain localisation. Because such behaviours are highly sensible to any modification of the model, computational strategies using model order reduction for damage are challenging and may be hazardous. It is generally recommended to control with circumspection, the accuracy of the solution provided by computations based on POD and to prefer some adaptive schemes as A Priori Hyper Reduction Method [START_REF] Ryckelynck | A robust adaptive model reduction method for damage simulations[END_REF] or POD coupled with Newton-Krylov algorithms [21]. PGD also offers such a flexibility as used recently by [START_REF] Metoui | The proper generalized decomposition for the simulation of delamination using cohesive zone model[END_REF], which introduced the coupling of PGD with a cohesive zone for delamination, or [14] for multi-scale computations for simulating problems with a rate-dependent damage model. These different authors have shown that model order reduction are promising to forecast damage evolution with limited computational cost.

In the current article, a different approach is proposed, based on the LATIN method [START_REF] Ladevèze | Nonlinear computational structural mechanics, mechanical engineering series[END_REF] which looks for an approximation of the solution on the entire time-space domain at every iteration. Therefore, this method is said to be nonincremental in time although discretisation schemes in space and in time are used as usual. As it handles an approximation of the whole time-space domain at every iteration, this approach is very suitable to include model order reduction techniques even for non-linear computations. A drastic decrease of the computational cost compared to a classical approach would allow to tackle sophisticated cases such as cyclic loading for example.

The LATIN-PGD has been developed for solving plasticity and visco-plasticity problems, even with cycling loading. From the first works [START_REF] Cognard | A large time increment approach for cyclic viscoplasticity[END_REF][START_REF] Ladevèze | The large time increment method for the analyse of structures with nonlinear constitutive relation described by internal variables[END_REF][START_REF] Ladevèze | Nonlinear computational structural mechanics, mechanical engineering series[END_REF], performances and robustness have improved thanks to [START_REF] Ladevèze | On reduced models in nonlinear solid mechanics[END_REF][START_REF] Ladevèze | On a multiscale computational strategy with time and space homogenization for structural mechanics[END_REF][START_REF] Ladevèze | On multiscale computational mechanics with time-space homogenization[END_REF]38,[START_REF] Relun | A model reduction technique based on the PGD for elastic-viscoplastic computational analysis[END_REF], leading to a mature approach for classical (visco-)plastic problems with parameters. The aim of this paper is to extend the approach for (visco-)plastic problems with unilateral damage. A major difficulty comes from the state laws which can not be transformed into linear relations through an internal variable transformation in order to employ the usual LATIN-PGD framework. An extension of the algorithm is introduced for the treatment of this non-linearity. The powerful numerical framework offered by LATIN-PGD seems to be robust and efficient to also predict damage including intrinsic model order reduction.

The article is structured as follows. In Sect. 2, the general problem is introduced in time-space domain with corresponding admissibility conditions. Constitutive formulations for the state and internal variables are presented including damage evolutions. In Sect. 3 the innovative LATIN algorithm is presented, incorporating damage. In Sect. 4, the PGD technique is described. The stress and strain are mathematically separated in terms of a part that depends on the plastic strain and the other part that depends on damage. The former is represented in a separable form of space and time which is estimated using a hybrid technique. Finally, two numerical examples are detailed in Sect. 5 to analyse the performances of the new algorithm.

The reference problem

The reference problem is a quasi-static isothermal evolution of a structure defined over time-space domain [0, T ] × Ω, assuming small perturbation. The structure is subjected to prescribed body forces f d , to traction forces F d over a part ∂ 2 Ω of the boundary. The structure is also subjected to prescribed displacements u d over the complementary part ∂ 1 Ω (Fig. 1). The state of the structure is defined by the set of field variables s = ε p , ε e , Ẋ, Ḋ, σ , Z, Y , where ε p is the inelastic part of the total strain ε corresponding to the displacement field u, that satisfies the strain partition relation ε = ε e + ε p , with ε e being the elastic strain. X represents the internal variables related to hardening. D is the isotropic damage variable associated with damaged state of the material. -σ is the Cauchy stress, Z is the set of variables conjugate to X, i.e. the thermodynamic forces associated to X. Y represents the energy release rate associated with damage and is the conjugate variable of D.

f d F d ∂ 2 Ω Ω ∂ 1 Ω u d
Fig. 1 Reference problem in domain Ω

Constitutive relations

The formulation of constitutive relations involves the description of the state equations and the formulation of the evolution laws.

State equations

The equations of state are obtained from the Gibbs free energy function ψ which can be decoupled into an elastic damage part ψ e and a plastic hardening part ψ p , ρψ = ρψ e (σ , D) + ρψ p (X) ,

where ρ is the mass density. The difference in the elastic state equation during tension and compression is due to the fact that during compression the material regains some stiffness, as some of the micro-defects are closed, thereby increasing the effective area. This can be represented by the effective modulus of elasticity during tension Ẽ+ , which is given by

Ẽ+ = E (1 -D)
, where E is the true modulus of elasticity, and by the effective elastic modulus during compression Ẽgiven by Ẽ-= E (1 -hD). The closure parameter h has values between 0 (complete stiffness recovery) and 1 (no stiffness recovery). This phenomenon is described by an elastic free energy function given as

ρψ e = 1 + ν 2E σ + ij σ + ij 1 -D + σ - ij σ - ij 1 -hD - ν 2E σ kk 2 1 -D + -σ kk 2 1 -hD , ( 2
)
where ν is the Poisson ratio. Hardening in classical sense is classified into isotropic and kinematic hardening. Isotropic hardening is the phenomenon where there is an increase in the radius of the yield surface within the Π plane and the centre of the yield surface remains constant. Kinematic hardening is where the yield surface translates in the Π plane without any change in radius. For metals, pure isotropic hardening can not describe Bauschinger effect but pure kinematic hardening can. For practical cases, however a combination of the two (mixed hardening) is used [START_REF] Lemaitre | Mechanics of solid materials[END_REF]. For mixed hardening case the plastic free energy function can be written as

ρψ p (X) = ρψ p (α, r ) = 1 2 α ij Cα ij + g (r ) , (3)
where α is internal variable corresponding to kinematic hardening, r is the internal variable for isotropic hardening. β and R are conjugate to α and r respectively. C is a material parameter for kinematic hardening and g (r ) is a function describing the isotropic hardening. These definitions give the state equations as

ε e ij = ρ ∂ψ e ∂σ ij = 1 + ν E σ + ij 1 -D + σ - ij 1 -hD - ν E σ kk 1 -D + -σ kk 1 -hD δ ij , (4a)
β ij = ρ ∂ψ p ∂α ij = Cα ij , (4b)
R = ρ ∂ψ p ∂r = g ′ (r ) , (4c)
Y = ρ ∂ψ e ∂ D = 1 + ν 2E σ + ij σ + ij (1 -D) 2 + h σ - ij σ - ij (1 -hD) 2 - ν 2E σ kk 2 (1 -D) 2 + h -σ kk 2 (1 -hD) 2 . ( 4d
)
It has to be noticed that the linear relationship ε e = C -1 σ between stress and elastic strain field, with C being the Hooke tensor, is modified in Eq. (4a) to include unilateral damage. An exponential isotropic rule is generally considered as,

g ′ (r ) = R ∞ 1 -exp (-γ r ) , (5)
with R ∞ and γ being material parameters describing isotropic hardening.

Evolution laws

The starting point of the formulation of the evolution equations is the second law of thermodynamics, written as Clausius-Duhem inequality,

σ ij ε p ij -R ṙ -β ij αij + Y Ḋ 0. ( 6
)
The evolution equations for the internal variables are obtained from a dissipation potential φ, which is a convex function of the associated variables to ensure that Eq. ( 6) is satisfied. The dissipation potential can be decoupled as

φ = φ p + φ d , (7)
where φ p represents the dissipation due to plasticity and φ d represents the dissipation due to damage. For visco-plasticity φ p can be defined according to Norton's law as

φ p = k n + 1 f p n+1 + , (8)
where k and n are material dependent viscous coefficient and exponent respectively and f p is the yield function. Consid-

ering τ ij = σ D ij 1-D -β ij , with σ D ij
being the deviatoric part of the stress tensor, the yield function for J 2 plasticity in a Marquis-Chaboche model is defined as,

f p = 3 2 τ ij τ ij + a 2C β ij β ij -R -σ y . (9)
The second term is specific to Marquis-Chaboche model and indicates material softening, where a is a material parameter and σ y is the yield stress. The evolution equations can then be written as

ε p ij = ∂φ p ∂σ ij = k f p n + ⎡ ⎣ 3 2 τ ij 3 2 τ ij τ ij ⎤ ⎦ 1 1 -D , ( 10a
) αij =- ∂φ p ∂β ij =-k f p n + ⎡ ⎣ - 3 2 τ ij 3 2 τ ij τ ij + a C β ij ⎤ ⎦ , (10b) ṙ
=- ∂φ p ∂ R =-k f p n + (-1) . (10c)
Similar to plasticity, the damage potential φ d can be defined as

φ d = k d n d + 1 f d n d +1 + , (11)
with k d and n d being the damage coefficient and damage exponent respectively. f d is the damage yield function given by,

f d = Y -Y 0 , (12)
where Y 0 describes the damage threshold. The evolution equation for damage is then given by

Ḋ = ∂φ d ∂Y = k D f d n d + . (13)
A similar damage evolution equation had been proposed in [START_REF] Bellenger | Plastic and viscoplastic damage models with numerical treatment for metal forming processes[END_REF], where an unbounded damage variable defined between 0 and ∞ was used. The damage evolution law (13)i sa l s o thermodynamically consistent for bounded damage variable defined between 0 and 1. Due to the presence of macroplasticity, it can be considered that damage occurs only during plastic deformation and there is no evolution of damage in the elastic regime, which defines the threshold Y 0 as

Y 0 = σ 2 y 2E . ( 14
)

Admissibility conditions

For the admissibility conditions of the reference problem, the following spaces and the corresponding vector spaces (denoted with subscript 0) are defined:

-The space U of the kinematic admissibility fields u such that

u |t=0 = u 0 and u = u d on ∂ 1 Ω. (15)
-The space S of the statically admissible fields σ such that,

[0,T ]×Ω σ : ε u * dΩ dt = [0,T ]×Ω f d • u * dΩ dt + [0,T ]×∂ 2 Ω F d • u * dS dt , ∀ u * ∈ U 0 . (16)
-The space E of the kinematically admissible fields ε such that ∃u ∈ U , ε =∇ sym u, which in weak form can be written as,

[0,T ]×Ω σ * : ε dΩ dt = [0,T ]×∂ 1 Ω σ * n • u d dS dt , ∀ σ * ∈ S 0 . (17)

Partial normal formulation

To use efficiently the LATIN method as a solver, the state equations should be described by linear operators. Hence it is necessary to introduce at this point a "normal" formulation of the constitutive relations [START_REF] Cognard | A large time increment approach for cyclic viscoplasticity[END_REF][START_REF] Ladevèze | Nonlinear computational structural mechanics, mechanical engineering series[END_REF]. For that purpose, the state law (4c) is transformed into linear relations by a change of variable

R = R ∞ r (18)
where, R and r are the new isotropic variables, with

r = 2 γ 1/2 1 -exp - γ r 2 , R = R ∞ R R ∞ γ 1/2 2 2 - R R ∞ γ 1/2 2 . ( 19
)
Using these two new variables, the yield function becomes,

f p = 3 2 τ ij τ ij + a 2C β ij β ij -R ∞ R R ∞ γ 1/2 2 2 - R R ∞ γ 1/2 2 -σ y , (20)
and the evolution Eq. (10c) can be written as,

ṙ =-k f p n + R R ∞ γ 2 -γ 1/2 . ( 21
)
However, it has to be noted that the elasticity law (4a) cannot be transformed into a linear relation, due to the damage term, and it will require a specific treatment.

The LATIN framework

Solving the reference problem has three difficulties; first being satisfying the global equilibrium of the structure, the second being taking into account the non-linear elastic law, and finally to solve the non-linear evolution equations. The LATIN method is a non-incremental solver in time that tackles these sets of equations iteratively on the whole time-space domain. The algorithm is initialised by the solution of the reference problem considering the loading is elastic and then plastic and damage corrections are added iteratively. For that purpose, two manifolds are introduced, based on a subdivision of the set of equations that must be solved:

-the first is the space A which belongs to the manifold of the admissibility conditions [START_REF] Galvis | Spectral multiscale finite element for nonlinear flows in highly heterogeneous media: a reduced basis approach[END_REF][START_REF] González | Computational vademecums for the real-time simulation of haptic collision between nonlinear solids[END_REF][START_REF] Grepl | Efficient reducedbasis treatment of nonaffine and nonlinear partial differential equations[END_REF], the linear state laws (4b, 18), and also the non-linear state law for damage (4d), that although being non-linear, can be postprocessed from the stress tensor and the damage variable at the end of each iteration; -the second is the space Γ which belongs to the manifold of the evolution equations (10a, 10b, 13, 21) and the elastic state law (4a) which was not linearisable due to damage.

The exact solution is then given by

s ex ∈ A ∩ Γ . (22)
To find s ex , an iterative algorithm is used, which consists, at each iteration, of a non-linear local stage and a linear global stage, seeking alternatively for an approximation of the solution field s over Γ and A. The relation between the two spaces are governed by linear operators called search directions. The iterative algorithm can be represented as follows:

s 0 ∈ A -→ ŝ1/2 ∈ Γ ••• -→ ŝi+1/2 ∈ Γ -→ s i+1 ∈ A ••• -→ s ex . ( 23
)
An overview of the extended method including unilateral damage is schematised in Fig. 2.

Initialisation

The beginning of the algorithm starts with an elastic initialisation. The loading is considered to be elastic and a solution set s 0 ∈ A is calculated. Due to the elastic assumption, all the internal variables X, Z, D and ε p are set equal to zero. All the boundary conditions are taken into account in the initialisation stage. At each subsequent iteration, correction to the elastic solution is computed.

Local stage

The objective of the local stage is to solve the evolution equations for internal variables that are local in space and non-linear. The elastic state law, being non-linear is also treated in this stage. The solution field s is separated into two parts s p ={ε p , ε e , Ẋ, σ , Z} and s d ={ Ḋ, Y }. Both solution fields are solved using linear operators having different properties. Knowing a solution set s i ∈ A, the objective is to find ŝi+1/2 ∈ Γ such that the local search directions (directions of ascent) are satisfied

⎡ ⎢ ⎢ ⎣ ε p i+1/2 - ε p i -Xi+1/2 -Ẋi εe i+1/2 -ε e i ⎤ ⎥ ⎥ ⎦ + B + ⎡ ⎢ ⎢ ⎣ σ i+1/2 -σ i Ẑi+1/2 -Z i σ i+1/2 -σ i ⎤ ⎥ ⎥ ⎦ = 0, ( 24a
)
Di+1/2 -Ḋi + b + Ŷi+1/2 -Y i = 0. ( 24b
)
Here, B + and b + are the directions of ascent for solution sets s p and s d respectively. Following [START_REF] Ladevèze | Nonlinear computational structural mechanics, mechanical engineering series[END_REF], the choice of these search directions is given by

B + -1 = 0 , b + -1 = 0. ( 25
)
The solution of the search direction equations ( 24) along with the evolution Eqs. ( 10) and ( 13) and the non-linear elastic law (4a) gives easily the complete set of solution ŝi+1/2 . Knowing the solution at the local stage ŝi+1/2 , the solution set s i+1 will be obtained in the global stage.

Initialisation elastic computation

Local stage: elastic state law evolution equations

Global stage: mechanical equilibrium admissibility conditions other state laws

PGD: Update ∆ε p i+1 = m j=1 ∆ λj(t)ε p j (x) ∆σ i+1 = m j=1 ∆ λj (t)Cε p j (x) PGD: Add a pair ∆ε p i+1 = λm+1(t)ε p m+1 (x) ∆σ i+1 = λm+1(t)Cε p m+1 (x) σ 0 ,ε e 0 i =0,m =0 { εp i+1/2 , εe i+1/2 , Xi+1/2 , σi+1/2 , Ẑi+1/2 }, { Di+1/2 , Ŷi+1/2 } ∆σi+1 =∆σ i+1 +∆σi+1
∆σi+1 estimated from the local stage information

{∆ λj } m j=1 = f εp i ,σ i , σi+1/2 , εp i+1/2 ζi <ζ tol ? λ m+1 (t), εp m+1 (x) orthonormalise εp m+1 (x), correct {λj (t)} m+1 j=1 λm+1 satisfying? rejection { εp i+1 =ε p i +∆ε p i+1 ,ε e i+1 , Ẋi+1 ,σ i+1 = σ i +∆σ i+1 ,Z i+1 }, { Ḋi+1 ,Y i+1 } ζi <ζ end ? yes no yes i = i +1 no no m = m +1 yes { εp ,ε e , Ẋ, σ, Z}, { Ḋ, Y }

Global stage

In the global stage, apart from the state laws and admissibility conditions, the solution set s i+1 ∈ A must also satisfy the descent search directions, i.e.

⎡ ⎢ ⎣ ε p i+1 - ε p i+1/2 -Ẋi+1 -Xi+1/2 ε e i+1 -ε e i+1/2 ⎤ ⎥ ⎦ -B - ⎡ ⎣ σ i+1 -σ i+1/2 Z i+1 -Ẑi+1/2 σ i+1 -σ i+1/2 ⎤ ⎦ = 0, (26a)
Ḋi+1 -Di+1/2 -b -Y i -Ŷi+1/2 = 0, (26b)
where

B -= H -0 0 C -1 .
The choice of the operator H - is such that it belongs to the tangent space associated with the solution set ŝi+1/2 in the manifold Γ . It is defined as the second derivative of the visco-plastic dissipation potential φ p . For mixed hardening, this is of the form

H -= ⎡ ⎢ ⎣ H σ H σβ H σ R H βσ H β H β R H Rσ H Rβ H R ⎤ ⎥ ⎦ . ( 27
)
For the sake of simplicity, only an approximation of the tangent is computed and the off-diagonal terms are considered to be zero, which results to

H -= ⎡ ⎣ H σ 00 0 H β 0 00 H R ⎤ ⎦ = H σ 0 0 H Z , (28)
where,

H Z = H β 0 0 H R . (29)
At each space-time point, the operator H -obtained from the viscoplastic potential is a symmetric matrix, positive but semi-definite so its inverse is not always defined. Hence it is necessary to regularise it by

H -= H -+ cM -1 , ( 30
)
where ζ is the regularisation coefficient taken herein to be 0.15, and

M = ⎡ ⎣ C 00 0 C 0 00R ∞ ⎤ ⎦ . (31)
The search direction operator for damage b -, is taken to be zero, that is the damage variable does not change in the linear stage.

The new difficulty, compared to former works with the LATIN method, is that the elastic state law (4a) is non-linear due to the presence of damage, leading to solve a non-linear problem at the global stage. This point is particularly tricky as it prevents the introduction of a model reduction strategy at this stage. The idea proposed herein is to transform this nonlinear problem into separate linear equations by decomposing stress and total strain into two parts depending on plastic deformation and damage respectively.

The quantities of interest at this point σ i+1 , ε e i+1 and ε p i+1 are represented in a corrective form at iteration i + 1as

∆σ i+1 = σ i+1 -σ i ,∆ ε e i+1 = ε e i+1 -ε e i and ∆ ε p i+1 = ε p i+1 - ε p i . (32)
The stress correction in the global stage at iteration i + 1 is separated into a part ∆σ ′ i+1 that depends on the plastic deformation and another part ∆ σ i+1 that is obtained from the non-linear state law,

∆σ i+1 = ∆σ ′ i+1 + ∆ σ i+1 . ( 33
)
Similarly the total strain correction ∆ε i+1 = ∆ε e i+1 +∆ε p i+1

at iteration i + 1 is also separated into a part ∆ε ′ i+1 that depends on plasticity and a part ∆ε i+1 that depends on damage,

∆ε i+1 = ∆ε ′ i+1 + ∆ εi+1 . (34)
Now, the corrective terms given by Eq. ( 32)helptorewrite the search direction equation C -1 as

∆σ i+1 = C∆ε e i+1 -∆R i+1 , (35)
where ∆R i+1 represents a residual stress term at iteration i + 1 and is given by

∆R i+1 = (σ i -σ i+1/2 ) -C ε e i -εe i+1/2 . ( 36
)
From the previous separations along with the additive strain decomposition relation, it can be established that

∆σ ′ i+1 + ∆ σ i+1 = C ∆ε ′ i+1 -∆ε p i+1 + C ∆ε i+1 -∆ε R i+1 , (37)
where ∆ε R i+1 can be interpreted as a residual strain obtained from non-linear state law at iteration i + 1 and is given by

∆ε R i+1 = C -1 ∆R i+1 . ( 38
)
On the other hand, if only the plastic part is considered, the search direction defined by Eq. (26a) is written as

∆ε p i+1 -H σ ∆σ i+1 + Δi+1 = 0, (39)
with

Δi+1 = H σ ( σ i+1/2 -σ i ) - ε p i+1/2 - ε p i . ( 40
)
Due to the separation of the stress tensor, only the part depending on plasticity is included in Eq. [START_REF] Nguyen | Certified real-time solution of parametrized partial differential equations, handbook of materials modeling[END_REF], which can thereby be re-written as

∆ ε p i+1 -H σ ∆σ ′ i+1 + Δi+1 = 0. ( 41
)
The objective of the global stage is then to calculate (∆ε p i+1 ,∆σ i+1 ,∆ε e i+1 ) by solving the weak form of the equilibrium equation.

PGD formulation of the global stage

The Proper Generalised Decomposition (PGD) is based on the idea that any quantity of interest which is dependent on several independent variables can be approximated as a sum of products of one-variable functions [START_REF]PGD in linear and nonlinear computational solid mechanics[END_REF][START_REF] Ladevèze | Nonlinear computational structural mechanics, mechanical engineering series[END_REF]. This decomposition includes thus an error due to the assumption of a separable form. For instance a field v(x, t) dependent on space and time variables is approximated as

v(x, t) ≈ n i=1 v x i (x)v t i (t). (42)
The usage of PGD in the LATIN framework lies in the fact that the quantities of interest are defined by Eq. ( 32)atevery successive iteration on the whole time-space domain and are estimated from a linear problem, allowing easily to estimate them as separable forms. Although, this type of representation is general and can be applied to any given field [START_REF] Ladevèze | Nonlinear computational structural mechanics, mechanical engineering series[END_REF], it is only efficient to solve a problem if the solution is separable. If is not the case, for example in the case of moving loads, some adaptations can be proposed to make this type of approximation efficient [START_REF] Allier | Proper generalized decomposition computational methods on a benchmark problem: introducing a new strategy based on constitutive relation error minimization[END_REF].

The boundary conditions have been taken into account in the elastic initialisation, so here the solution set searched in terms of corrections has to be kinematically admissible to zero, i.e. belonging to A 0 . The static admissibility condition ( 16) can be written as

[0,T ]×Ω ∆σ ′ i+1 : ε u * dΩ dt = 0 , ∀ u * ∈ U 0 , (43)
with

∆σ ′ i+1 = C ∆ε ′ i+1 -∆ε p i+1 . (44)
Also the part that is related to damage can be written in a similar weak form as

[0,T ]×Ω ∆ σ i+1 : ε u * dΩ dt = 0 , ∀ u * ∈ U 0 , (45)
with

∆ σ i+1 = C ∆ εi+1 -∆ε R i+1 . (46)

Separable representation of the quantities of interest

The quantities of interest represented in corrective forms are currently ∆σ ′ i+1 , ∆ σ i+1 , ∆ ε p i+1 , ∆ε e i+1 . The plastic strain rate is written as

∆ε p i+1 = λ p (t) ε p x ⇒ ∆ ε p i+1 = λp (t) ε p x . ( 47
)
The total strain that depends on plastic deformation and the corresponding displacement field can also be separated as

∆u i+1 = λ u (t) ū x ⇒ ∆ε ′ i+1 = λ u (t) ε x . (48)
For the following development, the indices are dropped for simplicity and used if necessary. Introducing the PGD approximation, the test field in Eq. ( 43) becomes u * = λ u * ū + λ u ū * with ū * ∈ U 0 and λ u * does not have any condition. The space-time problem then can be separated into a space problem and a time problem. The time problem consists in finding λ u such that ∀λ u * , [0,T ]

λ u λ u * dt Ω Cε : ε ū dΩ = [0,T ] λ p λ u * dt Ω Cε p : ε ū dΩ. ( 49
)
This gives that λ u and λ p are proportional and the simplest choice is to assume λ u = λ p = λ. With this choice of the time function, the space problem simply consists in seeking ū ∈ U 0 such that

Ω Cε : ε ū * dΩ = Ω Cε p : ε ū * dΩ, ∀ū * ∈ U 0 . (50)
This allows to define an operator E [START_REF] Relun | A model reduction technique based on the PGD for elastic-viscoplastic computational analysis[END_REF] such that

ε = Eε p . (51)
Equation ( 44) can then be written as

∆σ ′ = λC (E -I) ε p = λC ε p , (52)
with C = C (E -I) and I being the identity matrix. Finally including PGD in the global stage, at LATIN iteration i + 1, the description of the quantities of interest depending on the plastic strain becomes

∆ ε p i+1 = λ (t) ε p x , (53a)
∆σ ′ i+1 = λ (t) Cε p x . (53b)
The quantities that depend on damage, namely ∆ σ i+1 and ∆ εi+1 , are obtained from the weak form of Eq. ( 45) along with Eq. ( 46)giving

[0,T ]×Ω C ∆ε i+1 -∆ε R i+1 : ε u * dΩ dt = 0 , ∀ u * ∈ U 0 , (54)
where ∆ε R i+1 is a known quantity. ∆ε i+1 can be calculated using the same operator E as

∆ εi+1 = E∆ε R , (55)
and the stress tensor depending on damage is calculated as

∆ σ i+1 = C∆ε R . (56)
Finally the total stress and elastic strain tensors in corrective terms are obtained as

∆σ i+1 = ∆σ ′ i+1 + ∆ σ i+1 , (57a
)

∆ε e i+1 = ∆ε ′ i+1 + ∆ εi+1 -∆ε p i+1 . (57b)

Hybrid method to construct the PGD reduced-order basis

For calculating the approximation of the corrective terms in a separable form, it is possible to enrich the reduced-order basis by calculating a new product of a function of space and a function of time, or to re-use the space functions generated at previous iterations and only update the time functions.

As updating time functions is less expensive, this strategy is considered at the beginning of the global stage of every LATIN iteration, then enriching the reduced-order basis is done only if necessary as shown on Fig. 2.

Update of the time functions

Considering that m pairs have been generated at the end of the global stage of LATIN iteration i, at LATIN iteration i +1 the objective is to re-use the spatial basis. This phase is thus equivalent as a POD computation on the current PGD basis. The quantities of interest which depend on plastic deformation can be written as

∆ε p i+1 = m j=1 ∆ λ j (t) ε p j x , ∆σ ′ i+1 = m j=1 ∆λ j (t) Cε p j x . (58)
These updates of the time functions are calculated from the minimisation of a mechanical residual which is defined by the norm of the search direction operator, i.e.

{∆λ i } m j=1 m j=1 ∆ λ j ε p j -H σ m j=1 ∆λ j C ε p j + ¯ i+1 H -1 σ . ( 59
)
The minimisation problem gives a multi-variable differential equation, which is solved using discontinuous Galerkin method of order zero. More information on discontinuous Galerkin method for solving minimisation problems in LATIN method can be found in [START_REF] Ladevèze | The latin multiscale computational method and the proper generalized decomposition[END_REF][START_REF] Relun | Multiscale elasticviscoplastic computational analysis[END_REF]. Subsequently the total stress and elastic strain are calculated using Eq. (57). After the calculation of all the quantities of interest, if the correction provided by the update of the time functions is not satisfactory, a new space-time pair will be added. For a minimisation problem both in space and in time, a criterion for adding a space-time pair was suggested in [START_REF] Relun | A model reduction technique based on the PGD for elastic-viscoplastic computational analysis[END_REF]. Currently, the criterion to enrich the reduced-order basis is estimated by the saturation of the error indicator. Similar criterion was proposed in [START_REF] Heyberger | Multiparametric analysis within the proper generalized decomposition framework[END_REF], a modified version is given as,

ζ = ξ i -ξ i+1 ξ i + ξ i+1 . ( 60
)
If this indicator is larger than a particular pre-defined value ζ tol , the approximation for this LATIN iteration is considered satisfactory. Then, the internal variables are calculated before pursuing with the subsequent LATIN iteration. At the contrary, if the indicator is lower than ζ tol , a space-time PGD pair is added to enrich the reduced-order basis.

Addition of space-time modes

To enrich the reduced-basis, a hybrid strategy is used, that combines a Galerkin formulation to approximate the space fields and a minimisation technique to estimate the corresponding time functions [START_REF] Ladevèze | On reduced models in nonlinear solid mechanics[END_REF]. The search direction equation during the stage of enrichment becomes

∆ ε p i+1 -H σ ∆σ ′ i+1 + Δi+1 = 0, (61)
with

Δi+1 = H σ σ i+1/2 -σ up i+1 - ε p i+1/2 - ε p , up i+1 , (62)
where ε p , up i+1 and σ up i+1 are the quantities obtained from the update stage. The quantities of interest are written in a separable form as

∆ε p i+1 = ∆ λm+1 (t) ε p m+1 x , ∆σ ′ i+1 = ∆λ m+1 (t) Cε p m+1 x . (63)
A fixed point algorithm is used to calculate the spacetime pair. All the superscripts and subscripts that have been used previously are dropped for simplicity for the following development to calculate the space function by a Galerkin technique. The strain partition relation combined with the state equation is written as

∆ ε′ = ∆ ε p + C -1 ∆ σ ′ , (64)
where ∆ε p is obtained from Eq. ( 61) and ∆σ ′ = λ σ with σ = λCε p . The kinematic admissibility condition in rate form is written as

[0,T ]×Ω ∆ ε′ : σ * dΩ dt = 0, ∀ σ * ∈ S 0 , ( 65
) which leads, introducing • = [0,T ]
• dt and σ * = λ σ ,to

Ω H σ λ 2 σ + λλ C -1 σ -Δλ : σ * dΩ = 0, ∀ σ * ∈ S 0 . ( 66
)
By defining ε such that

ε = W -1 σ -δ, (67)
with

W -1 = H σ λ 2 + λλ C -1 and δ = Δλ , it is estab- lished that Ω ε : σ * dΩ = 0, ∀ σ * ∈ S 0 . ( 68
)
The static admissibility to zero of σ reads

Ω σ : ε( ū * ) dΩ = 0, ∀ū * ∈ U 0 , (69)
which is rewritten using Eq. ( 67) and by introducing displacement ū ∈ U 0 such that ε( ū) = ε:

Ω Wε( ū) : ε( ū * ) dΩ =- Ω W δ : ε( ū * ) dΩ = 0, ∀ū * ∈ U 0 . ( 70
)
This problem is solved classically to obtain ū and thereby calculate the associated strain ε. The space function ε p m+1 is thereby calculated as,

ε p m+1 = 1 λ λ ε -λ λ C -1 W( ε + δ) . ( 71
)
Then the time function λ m+1 is calculated using a minimisation technique similar to the update stage

λ m+1 = arg min λ m+1 λm+1 ε p m+1 -H σ λ m+1 Cε p m+1 + ¯ i+1 H -1 σ . ( 72
)
The spatial basis is orthonormalised using Gram-Schmidt algorithm [START_REF] Ruhe | Numerical aspects of Gram-Schmidt orthogonalization of vectors[END_REF]. The new space function is orthonormalised with respect to the previously existing spatial bases, and in the numerical process all the former time functions are updated. The new time function is also modified and the corresponding space-time pair may be rejected if the corresponding modified time function has an insignificant norm. Finally, the global stage is concluded by estimating the hardening and damage variables.

Internal variables at the global stage

As the internal variables are local in space, they are calculated without any PGD approximation, but by simply solving first order ordinary differential equations in time locally at each Gauss point (GP). The hardening variables represented by X and Z are obtained by the state laws and search direction operators. The "normal" formulation allows to represent the hardening state laws as

Z i+1 = X i+1 , (73)
where is a linear operator describing the "normal" formulation. The search direction equation for hardening variables (26a) combined with the state equation ( 73) can be written as

-( Ẋi+1 -Xi+1/2 ) = H Z (Z i+1 -Ẑi+1/2 ) = H Z ( X i+1 -Ẑi+1/2 ). ( 74
)
For damage, the search direction operator b -is taken to be zero, leading to

Ḋi+1 = Di+1/2 , ( 75
)
and finally, the energy release rate Y i+1 is calculated through Eq. (4d).

Relaxation of the solution field and convergence criterion

To ensure convergence of the algorithm, the global step is modified by adding a relaxation step. From the solution set si+1 obtained at the end of the iterative PGD process, the approximation provided by the global stage is defined as

µs i+1 + (1 -µ) s i .
The relaxation parameter µ is chosen to be 0.8. The convergence of the iterative algorithm is determined by a relative LATIN indicator. This indicator is basically the distance between the local solution and the global solution, given by

ξ = ŝ p i+1/2 -s p i+1 ŝ p i+1/2 + s p i+1 , (76)
with

s p 2 = [0,T ]×Ω σ : H σ σ + Z : H Z Z + ε p : H -1 σ ε p + ε e : Cε e + Ẋ : H -1 Z Ẋ dΩ dt. ( 77
)
It has to be noted that the relaxation parameter µ, the regularisation coefficient ζ , and the error indicator tolerance ξ mostly affect the rate of convergence and not the convergence itself. The values of these parameters that have been used here are from the numerical experiments published in previous works [START_REF] Ladevèze | Nonlinear computational structural mechanics-new approaches and non-incremental methods of calculation, mechanical engineering series[END_REF][START_REF] Relun | A model reduction technique based on the PGD for elastic-viscoplastic computational analysis[END_REF].

Numerical examples

The innovative algorithm has been tested on academic examples in one-dimensional and two-dimensional cases.

Bar under traction

The first test problem considered is a bar in traction represented in Fig. 3. The geometry of the structure is defined by the length L = 1000 mm and area of cross section A = 100 mm 2 . The bar is constrained at x = 0. At x = L, a sinusoidal prescribed displacement loading of amplitude 1.210 -3 L is applied with a time period ∆T = 10 s, for 20 The discretisation in space is done with classical finite element scheme and 90 linear bar elements are used to discretise the structure such that each part has the same number of elements. For the time discretisation, the time step chosen is 0.1 s, so for the given loading there are 100 time elements per cycle and in total 2000 time elements. In the computational model, the criterion ζ tol to enrich the reduced-order basis is taken to be 0.1. The algorithm is stopped if the saturation parameter ζ is lower than 10 -4 . The search direction operators have been defined in [START_REF] Ladevèze | Nonlinear computational structural mechanics, mechanical engineering series[END_REF] and the formulation is extended here to incorporate damage: The convergence of the algorithm is measured by the quantity ξ defined by Eq. (76) and depicted in Fig. 4. The LATIN indicator ξ stagnates at the end which reflects that neither updating the reduced-order basis, nor its enrichment will improve the approximated solution. A finer time step size could decrease the final value of ξ . After the convergence of the algorithm, the quantity of interest ε p is obtained as the sum of products of three space-time modes, where the spatial basis is orthonormal (Fig. 5). Due to the rejection of the insignificant space-time modes after orthonormalisation a maximum of three modes is obtained. The time functions generated have increasing amplitude with respect to time.

H σ = kn f p n-1 1 (1 -D) 2 , (78a)
H β = kn f p n-1 -sign σ 1 -D -β + a C β 2 + k f p n a C , (78b)
H R = kn f p n-1 γ 1 - Rγ 1/2 2R ∞ 2 + k f p n γ 2R ∞ . (78c)
The quantities of interest (internal variables in rate form and the corresponding associated variables) for the three sections of the bar are depicted in Fig. 6. Due to the particular loading, the stress is constant all over the structure for a particular time, and with respect to time the stress amplitude decreases. The plastic strain rate and kinematic hardening variables have increasing amplitude with respect to time for mat. 1. For mat. 2 and mat. 3, the amplitudes decrease. The isotropic hardening variables show similar behaviour, i.e. for mat. 1 there is a monotonic increase in the amplitude of ṙ , and for mat. 2 and mat. 3, the amplitude of ṙ decreases. The associated variable R increases for all the three materials.

The quantities that indicate loss of stiffness in the material under the direct influence of unilateral condition are depicted in Fig. 7. The damage variable D does not practically increase during the compressive part of the loading and all the evolution takes place during the tensile part. The energy release rate Y during compression is much less than during tension. The elastic strain ε e also shows a higher value during tension than in compression. Finally after 20 cycles the values of D are 0.22, 0.18 and 0.15 for mat. 1, mat. 2 and mat. 3 respectively.

For general engineering problems both plastic deformation and damage are often highly localised phenomena. The presence of stress raisers (e.g. notches, holes, defects and such others) in most engineering components, concentrates the effect of plasticity and damage to a very limited region, while the rest of the structure remains mostly undamaged and elastic. To depict the localised nature of plasticity and damage, classical two-dimensional mechanical problems are considered in Sects. 5.2 and 5.3.

"L" shaped structure

The test problem considered here is an "L" shaped structure subjected to a concentrated load represented in Fig. 8. The geometry of the structure is defined by the length L = 120 mm, and width W = 20 mm. The structure is filleted at the inside with fillet radius r = 5 mm to avoid stress singularity. The thickness of the structure is taken to be 1 mm. A prescribed sinusoidal displacement U d (t) of amplitude 1.5 mm is applied with a time period ∆T = 10 s, for 5 cycles.

The material properties used are the same as in the onedimensional example of Sect. 5.1, and are given in Table 1, with σ y = 85 MPa. The discretisation in space is done with classical finite element scheme. 277 linear two-dimensional quadrilateral isoparametric plane stress elements with 4 Gauss points per element are used to discretise the structure, which generates 337 nodes. For the time discretisation, the time step chosen is 0.2 s. In the computational model, the criterion to enrich the reduced-order basis ζ tol is taken to be 10 -2 . The search direction operators H σ , H β and H R have to be defined for a two-dimensional problem. Theoretically, any symmetric positive definite matrices, used as the search direction operators, suffice to guarantee the convergence of the algorithm, however, the fastest route to convergence is to use the search direction operators associated to the tangent . The search direction operator has been specified in [START_REF] Relun | A model reduction technique based on the PGD for elastic-viscoplastic computational analysis[END_REF] for perfect plasticity problems. This idea is extended to include damage and hardening variables,

H σ = kn f p n-1 1 (1 -D) 2 3/2τ J 2 ⊗ 3/2τ J 2
+ k f p n 1 (1 -D) 2 3/2 J 2 I -3/2 τ ⊗τ J 2 J 2 2 , (79a)
H β = kn f p n-1 - 3/2τ J 2 + a C Iβ ⊗ - 3/2τ J 2 + a C Iβ + k f p n 3/2 J 2 I -3/2 τ ⊗τ J 2 J 2 2 a C I, (79b)
H R = kn f p n-1 γ 1 - Rγ 1/2 2R ∞ 2 + k f p n γ 2R ∞ , (79c)
with J 2 = 3 2 τ • τ . Unlike the one-dimensional problem introduced in Sect. 5.1, the calculation of the search direction operators is considerably expensive, so they are not calcu-lated at every iteration.The Reference Point Method, a hyper reduction technique, has been proposed in [START_REF] Capaldo | The reference point method, a hyper reduction technique: application to PGD-based nonlinear model reduction[END_REF], and will be introduced in the future. In the current stage, search direction operators are only updated if the difference between two successive global solutions becomes lower than a certain tolerance. This can be measured by the criterion ζ defined by Eq. ( 60) with the tolerance being 10 -2 . The decrease in the LATIN indicator ξ is shown in Fig. 9. A total of 35 modes are generated to obtain ξ = 9 × 10 -4 . An indicator value of 5 × 10 -2 is obtained at the end of iteration 1. After 21 iterations, the LATIN indicator becomes 9 × 10 -3 with 8 PGD modes. The indicator seems to saturate at ξ = 8 × 10 -4 , and the total number of modes added saturates to 40. This saturation of LATIN indicator at a relatively high value is due to the fact that a coarse temporal mesh has been used. For a finer temporal mesh this LATIN indicator can be lowered further, as shown in the next numerical example (Sect. 5.3).

The quantity of interest to describe the contribution of the plastic deformation is the accumulated plastic strain given by (80)

The distribution of the accumulated plastic strain at time t = T is depicted in Fig. 10. It is distributed on both sides of the vertical web, with the inner side being predominant and the maximum at the filleted corner of the structure.

Although the main part of the structure has not been permanently deformed, some plastic deformation can be observed on the inner side of the vertical web, and to a lesser degree, on the outer side of the web. The distribution of damage variable D at the end of the loading (t = T ) is represented in Fig. 11. The distribution is similar to the one of the accumulated plastic strain, i.e. maximum at the filleted corner. The spread of damage in this region of interest after each load cycle is also represented in Fig. 11. The evolution of damage with respect to time for The distribution of the residual von Mises stress at the end of the loading is shown in Fig. 13, which similar to the distribution of damage and plastic strain, is concentrated at the filleted corner.

Plate with a hole

The numerical example considered here is a rectangular plate of length L = 40 mm and width W = 60 mm, with a central circular hole of diameter φ = 20 mm. The plate is subjected to uniformly distributed sinusoidal displacements of amplitude 0.012 mm and time period ∆T = 10 s on both ends, as shown in Fig. 14, for 10 cycles. Due to symmetry of the structure only a quarter of the plate, with symmetric boundary conditions on the interior of the plate, is considered for the analysis. The thickness of the plate is considered to be 1 mm.

The material properties used are the same as in Table 1 with σ y = 85 MPa. The discretisation in space is done with classical finite element scheme. 129 linear two-dimensional quadrilateral isoparametric plane stress elements with 4 Gauss points per element are used to discretise the struc- The decrease in the LATIN indicator ξ for the virgin structure is shown in Fig. 16. A total of 34 modes are generated to obtain ξ = 9 × 10 -5 . Unlike Sect. 5.2, the LATIN indicator obtained is relatively low. This expected result is due to the usage of finer temporal mesh in this example and it confirms the fact that the LATIN indicator is dependent on temporal discretisation [START_REF] Ladevèze | Nonlinear computational structural mechanics-new approaches and non-incremental methods of calculation, mechanical engineering series[END_REF] (Fig. 17).

For the pre-damaged structure, compared to the virgin structure, 38 PGD modes are generated to obtain ξ = 9 × 10 -5 . Also a total of 83 iterations are needed, compared to 75 iterations needed for the virgin structure.

The distribution of damage variable D, for both the cases, at the end of the loading (t = T ) are represented in Fig. 18. The spread of damage for both cases after certain load cycles in the regions of interest are also depicted in Fig. 18.Itisquite evident from Fig. 18, that there is a shift of the maximum damaged point in the pre-damaged case, compared to the virgin case. The weakest Gauss point for the virgin structure GP 1 (Fig. 15) is shifted to GP 2 (Fig. 15) for the pre-damaged structure. The evolution of damage for GP 1 and GP 2 for the virgin and pre-damaged cases respectively is shown in Fig. 19, with a maximum of 0.05 for the virgin case and 0.07 for the pre-damaged case.

Conclusion

In this work, the LATIN-PGD has been developed for solving damage-visco-plasticity problems under cyclic loading. A major difficulty comes from the state laws which can not be transformed into linear relations through an internal variable transformation in order to use the classical LATIN-PGD framework. An extension of the algorithm has been intro-duced for the treatment of this non-linearity and seems to be robust and efficient.

Further works are in progress for the damage detection in the case of fatigue loading involving large number of cycles. For that purpose a new computational framework will be introduced, allowing a multi-scale description of the fields along the time and leading to a strong decrease of the computational cost.

Fig. 2

 2 Fig. 2 Scheme of the innovative LATIN algorithm including unilateral damage (i: iteration of LATIN method, m: number of PGD pairs involved in the separable form)

Fig. 3 A

 3 Fig. 3 A bar in traction

Fig. 4 Fig. 5

 45 Fig. 4 Evolution of the LATIN indicator with respect to the number of PGD pairs or LATIN iterations for the bar problem

1 Fig. 6

 16 Fig. 6 Evolution of quantities of interest in mat. 1 (blue line), mat. 2 (red line) and mat. 3 (green line) for the bar under cyclic loading. (Color figure

Fig. 7 Fig. 8

 78 Fig. 7 Evolution of quantities describing damage in mat. 1 (blue line), mat. 2 (red line) and mat. 3 (green line) for the bar under cyclic loading. (Color figure online)

Fig. 9

 9 Fig.[START_REF]PGD in linear and nonlinear computational solid mechanics[END_REF] Evolution of the LATIN indicator with respect to the number of PGD pairs or LATIN iterations for the "L" shaped structure

  ,τ : ε p x,τ 1/2 dτ, ∀t ∈ [0, T ] .

Fig. 11 Fig. 12

 1112 Fig. 11 Damage distribution in the "L" shaped structure at t = T and the spread of damage in the region of interest after certain load cycles

Fig. 13 2 Fig. 15 FE

 13215 Fig. 13 Distribution of residual von Mises stress in the "L" shaped structure at t = T

Fig. 14 AFig. 16 Fig. 17 Fig. 18 Fig. 19

 1416171819 Fig.14A plate with a hole subjected to distributed loads and the symmetric part considered for analysis

Table 1

 1 Material properties

	E	134,000 MPa
	ν	0.3
	R ∞	30 MPa
	γ	2MPa
	C	5500 MPa
	a	250 MPa
	k D	2.778 MPa -n d s -1
	n D	2
	K	1220 MPa s 1/n
	n	2.5
	kK -n
	h	0.2
	σ y	mat. 1	mat. 2	mat. 3
		80 MPa	82.5 MPa	85 MPa
	cycles. The structure is composed of three different elasto-
	viscoplastic materials denoted by mat. 1, mat. 2 and mat. 3.

The materials are distinguished by the yield stress σ y .T h e material considered is a Cr-Mo steel at 580 • C and its properties are given in Table 1

[START_REF] Lemaitre | Engineering damage mechanics: ductile, creep, fatigue and brittle failures[END_REF]

.

Acknowledgements

The authors are grateful to the German Research Foundation (DFG) for funding the research through International Research Training Group 1627.