# Effective reduced-order approximation for fast and robust MCRE-based parametric identification of nonlinear history-dependent material laws

Mainak Bhattacharyya, Ludovic Chamoin

► To cite this version:

Mainak Bhattacharyya, Ludovic Chamoin. Effective reduced-order approximation for fast and robust MCRE-based parametric identification of nonlinear history-dependent material laws. Computational Mechanics, 2026, 10.1007/s00466-026-02754-1. hal-05519744

HAL Id: hal-05519744

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

Submitted on 19 Feb 2026

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Effective reduced-order approximation for fast and robust MCRE-based parametric identification of nonlinear history-dependent material laws

Mainak Bhattacharyya<fnref n="1" /> · Ludovic Chamoin<fnref n="1" />

Received: 16 November 2025 / Revised: 15 January 2026 / Accepted: 22 January 2026

© The Author(s) 2026

## Abstract

This article essentially addresses the numerical frugality of model updating procedures using reduced order modelling. Nonlinear material behaviour is tackled in the article with the focus being on elasto(visco)-plasticity and elasto(visco)-plastic-damage. A proper generalised decomposition based formulation is introduced that separates the governing equations into sundered space and time problems, thereby providing low fidelity approximations. A proper orthogonal decomposition based model reduction method is also introduced to tackle variable Hooke's tensor for softening material behaviour. The proposed methodologies are exemplified through several numerical tests in order to assess and validate performance.

**Keywords** Parametric identification · Modified Constitutive Relation Error · Reduced order methods · Nonlinear constitutive laws · Generalised standard materials

## 1 Introduction

In the context of structural health monitoring, it is essential to detect defects early on and track their growth, ensuring the integrity of structures. Such a philosophy has applications in several industrial sectors, such as transport or energy (to name a few). This necessitates the combination of numerical simulations and physical measurements, which gave rise to dynamic data driven application systems (DDDAS) (see [1, 2]), involving methods enabling real-time, dynamic information exchange between the physical system and its corresponding numerical counterpart or virtual twin (see [3] for details on twins). DDDAS attempts not only to continuously predict the evolution of the relevant physical phenomena, but also dynamically update the computational

model in real-time. Such an approach assists in system adjustment as well as system diagnosis.

The coupling between measured data and constitutive modelling is generally performed by using the data to identify the parameters of the constitutive model. The identification procedure based on inverse approach exploits the affluent measurement data, which can be obtained through either digital image correlation (DIC) (see [4]), vibration measurement e.g. by means of accelerometers, or optic fibre sensing (see [5]). Out of pure machine learning strategies, the parametric identification strategy essentially involves model updating based on the adjustment of the model parameters that reduces in some manner the difference between the model prediction and the measurement data. Several physics-based inverse analysis methods for model identification from full-field measurements are detailed in [6], such as the virtual field method (VFM) in [7], reciprocity gap method (RGM) in [8], finite element model updating method (FEMU) in [9], equilibrium gap method (EGM) in [10], constitutive relation error (CRE) method in [11], or modified constitutive relation error (MCRE) method in [12].

In the case of the MCRE method, the full knowledge of the measurement data and boundary conditions is not required (see [13]). Also, MCRE is able to handle spurious data due to measurement noise quite effectively (see [12–14]). The principle of MCRE consists in dividing the

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

<footnote>Mainak Bhattacharyya and Ludovic Chamoin contributed equally to this work.</footnote>

<footnote>Ludovic Chamoin
ludovic.chamoin@ens-paris-saclay.fr</footnote>

<footnote>Mainak Bhattacharyya
mainak.bhattacharyya@ens-paris-saclay.fr</footnote>

<footnote n="1">¹ LMPS - Laboratoire de Mécanique Paris-Saclay, Université Paris-Saclay, CentraleSupélec, ENS Paris-Saclay, CNRS, 4 Avenue des Sciences, Gif-sur-Yvette 91190, France</footnote>
relations which characterise the problem into a reliable set (mechanical equilibrium) and a less reliable set (measurement data, constitutive relation, and/or boundary conditions). The reliable set is exactly satisfied, and the less reliable set is approximately satisfied through a functional minimisation. In case of CRE, only the constitutive relation is considered to be less reliable. A variant of this, referred to as optimal control (OC) method (see [15]), uses only measurement data as the less reliable set, while the reliable set includes both mechanical equilibrium and constitutive relations. MCRE is a hybrid approach, relying on both physics-based and data-based information, providing a good compromise between purely data-based approaches which require high data quantity, and model-based approaches which require low data quantity. Nevertheless, all three approaches provide estimates of the mechanical states, by means of the computations of admissible fields. [14, 16, 17] depict the usage of MCRE in the case of nonlinear parametric identifications leads to a minimisation at the mechanical model level translating to a non-standard finite element formulation, with a global (space-time) nonlinear optimisation problem. In this article, the MCRE procedure is employed for the inverse analysis.

In the nonlinear range, the resolution of the MCRE procedure incorporates an iterative procedure similar to the LATIN method (see [18]). In this article, the intent is to propose a proper generalised decomposition (PGD) based reduced order model (ROM) in space and time to be integrated in the LATIN type methodology for the MCRE process. The advantage of PGD is numerical cost effectiveness as ROMs provide low dimensional approximations for high fidelity models (see [19]); this is thus a good option to tend to real-time inverse applications with MCRE. Also, the LATIN method is inherently conducive to PGD implementation in space and time, as it is a non-incremental but iterative procedure (see [18, 20]). [21, 22] have used parametric PGDs for model updating of linear problems. A LATIN-PGD approach was used in [23] in space and time for nonlinear model updating using the aforementioned OC method. In the current article, space-time LATIN-PGD is used to solve the MCRE problem. However, unlike [23] where PGD was based on a minimisation formulation, in the current article a Galerkin formulation is used because of its computational frugality and formulation simplicity (see [24] for details on PGD formulations).

Also, as damageable materials are also studied in this article, a proper orthogonal decomposition (POD) based ROM strategy is used to approximate the variable Hooke tensor. In case of POD, certain full-order problems are computed in order to extract relevant information which can then be used to approximate the solution more efficiently (see [25]). The usage of POD in the field of mechanics deals with the creation and truncation of POD basis and coefficients in a least square sense from the snapshots of the full-order solution (see [26]).

<figure><img src="image_3.png" /><figcaption>Fig. 1 Representative schematic of the identification problem</figcaption></figure>

The article starts with the problem definition (section 2), describing the governing equations. Section 3 thereafter describes the basic concept of the MCRE method, in a classical way. Section 4 provides the PGD based formulation for the MCRE procedure, which is then extended in section 5 to damageable behaviour by coupling with a POD based method. Finally the article is concluded through Section 7.

## 2 Problem definition

To elucidate the basic theology of the MCRE procedure, we consider the domain $\Omega$ as shown in fig. 1. It consists of a boundary $\partial\Omega$, divided into three separate boundaries $\partial_u\Omega$, $\partial_f\Omega$, $\partial_{\emptyset}\Omega$, such that $\partial\Omega = \partial_u\Omega \cup \partial_f\Omega \cup \partial_{\emptyset}\Omega$. Dirichlet boundary conditions are applied on $\partial_u\Omega$, Neumann boundary conditions are applied on $\partial_f\Omega$, while $\partial_{\emptyset}\Omega$ contains no boundary information. Also, there exists a sub-domain $\Omega_m$, such that $\Omega_m \subseteq \Omega$, which represents the region where the kinematic data is measured.

For the direct mechanical problem, consider a time-dependent scenario where the time variable $t$ is defined such that $t \in [0, T]$, with 0 and $T$ the lower and upper temporal bounds. A prescribed displacement field $u_d$ is applied on $\partial_u\Omega$, and a prescribed traction field $f_d^s$ is applied on $\partial_f\Omega$. A body force field $f_d^v$ can also be considered to be applied in $\Omega$.

The kinematic admissibility defines the space $U_{ad}$ of displacement field $u$ satisfying the Dirichlet boundary conditions, that is

$$ \underline{u}(t)|_{\partial_u \Omega} = \underline{u}_d(t), \quad \forall t \in [0, T]. \qquad (1) $$

The static admissibility defines the space $S_{ad}$ of stress field $\underline{\sigma}$ (Cauchy stress) satisfying the equilibrium equation

$$ \begin{aligned} & \iint_{\Omega \times [0, T]} \underline{\sigma} : \underline{\varepsilon}(\underline{v}) \, d\Omega dt \\ &= \iint_{\Omega \times [0, T]} \underline{f}_d^v \cdot \underline{v} \, d\Omega dt + \iint_{\partial_f \Omega \times [0, T]} \underline{f}_d^s \cdot \underline{v} \, dS dt, \quad \forall \underline{v} \in U_{ad,0}, \end{aligned} \qquad (2) $$
where $\underline{\varepsilon}(\underline{v})$ is the infinitesimal strain tensor associated with the test or virtual displacement field $\underline{v}$, such that $\underline{\varepsilon}(\underline{v}) = \nabla^s \underline{v}$. The vector space $U_{ad,0}$ is the space associated with the test field vanishing on $\partial\Omega \setminus \partial f\Omega$.

The constitutive behaviour, dictating the behaviour of the material is given by

$$ \underline{\sigma} = \frac{\partial \psi}{\partial \underline{\varepsilon}^e} \left( \text{or } \underline{\varepsilon}^e = \frac{\partial \psi^*}{\partial \underline{\sigma}} \right), \quad \underline{Z} = \frac{\partial \psi}{\partial \underline{X}} \left( \text{or } \underline{X} = \frac{\partial \psi^*}{\partial \underline{Z}} \right) \quad (3a) $$

$$ \underline{\sigma} = \frac{\partial \phi}{\partial \dot{\underline{\varepsilon}}^p} \left( \text{or } \dot{\underline{\varepsilon}}^p = \frac{\partial \phi^*}{\partial \underline{\sigma}} \right), \quad \underline{Z} = -\frac{\partial \phi}{\partial \dot{\underline{X}}} \left( \text{or } \dot{\underline{X}} = -\frac{\partial \phi^*}{\partial \underline{Z}} \right). \quad (3b) $$

which refers to the formulation of generalized standard materials. Here, eq. (3a) refers to state equations obtained from the free energy function $\psi$, or the dual free energy function $\psi^*$. Equation (3b) describes evolution laws obtained from the dissipation potential $\psi$, or the dual dissipation potential $\psi^*$. The free energy function and the dissipation potential are convex scalar functions, and are related to their dual counterparts through the Legendre-Fenchel transform. The strain tensor $\underline{\varepsilon}$ is additively decomposed into an elastic part $\underline{\varepsilon}^e$ and an inelastic part $\underline{\varepsilon}^p$, i.e. $\underline{\varepsilon} = \underline{\varepsilon}^e + \underline{\varepsilon}^p$. The additional (inter nal) variables that describe the material state are represented by the set of all primal variables $\underline{X}$ and their associated dual variables $\underline{Z}$. The free energy function and the dissipation potential are parametrised by the set of parameters $\underline{\theta}$.

### 3 Concept of MCRE for history-dependent materials

The MCRE aims to identify the parameters $\underline{\theta}$ of constitutive relations that fit the best to the experimental data (considered here to be the locally measured displacement field $\tilde{\underline{u}}$). This means that

$$ \underline{u}(t)|_{\Omega_m} = \tilde{\underline{u}}(t), \quad \forall t \in [0, T]. \quad (4) $$

Here, the Dirichlet and Neumann boundary conditions are considered to be known, however the methodology can be easily extended for unknown boundary conditions.

The general idea of the MCRE-type minimisation is to split the equations of section 2 following the reliability of information. The constitutive relation is considered as the unreliable equation and the associated constitutive relation error is defined as

$$ \xi^2 = \int_0^T \int_{\Omega} \eta_{\psi} d\Omega dt + \int_0^T \int_0^T \int_{\Omega} \eta_{\phi} d\Omega d\tau dt, \quad (5) $$

where the residual errors are defined as

$$ \eta_{\psi} = \psi(\underline{\varepsilon}^e, \underline{X}) + \psi^*(\underline{\sigma}, \underline{Z}) - \underline{\sigma} : \underline{\varepsilon}^e - \underline{Z} \cdot \underline{X}, \quad (6a) $$

$$ \eta_{\phi} = \phi(\dot{\underline{\varepsilon}}^p, -\dot{\underline{X}}) + \phi^*(\underline{\sigma}, \underline{Z}) - \underline{\sigma} : \dot{\underline{\varepsilon}}^p + \underline{Z} \cdot \dot{\underline{X}}. \quad (6b) $$

The CRE concept essentially involves minimising the CRE $\xi$ over admissibility spaces to find the parameters $\underline{\theta}$ of the constitutive model. In the first idea of the method, observations were directly included in the definition of the kinematically admissible space, leading to weak robustness to measurement noise. To improve this robustness, the MCRE procedure does not impose the observations in the admissible space. This defines the MCRE functional as

$$ \zeta^2 = \xi^2 + \frac{\gamma}{2} \int_0^T \int_{\Omega_m} \|\underline{u} - \tilde{\underline{u}}\|^2 d\Omega_m dt, \quad (7) $$

where $\gamma$ is a scaling factor.

The inverse problem thus consists of finding the optimal parameters such that

$$ \underline{\theta}_{opt} = \min_{\underline{\theta}} \left( \min_{\underline{x} \in \mathcal{A}_{ad}} \zeta^2(\underline{x}, \underline{\theta}) \right) \quad (8) $$

where $\mathcal{A}_{ad}$ is the admissibility space such that $\mathcal{A}_{ad} = U_{ad} \times S_{ad}$, while $\underline{x}$ represents the whole set of state variables, i.e.

$$ \underline{x} = \left\{ \underline{\varepsilon}^e \left( = \underline{\varepsilon} - \underline{\varepsilon}^p \right), \underline{\varepsilon}^p, \underline{\sigma}, \underline{X}, \underline{Z} \right\} $$

The minimisation is performed sequentially and iteratively. For a given iteration $i$:

* The first step (also referred to as the basic problem) involves computing the optimal admissible state $\underline{x}$ for a given parameter set $\underline{\theta}_{i-1}$, i.e.

$$ \underline{x}_i = \min_{\underline{x} \in \mathcal{A}_{ad}} \zeta^2(\underline{x}, \underline{\theta}_{i-1}). \quad (9) $$

* The second step (also referred to as the identification problem) involves computing the updated material parameter by minimising the MCRE for the previously calculated admissible set, i.e.

$$ \underline{\theta}_i = \min_{\underline{\theta}} \zeta^2(\underline{x}_i, \underline{\theta}). \quad (10) $$

#### 3.1 The basic problem

The basic problem is the most computationally expensive part, and the resolution is performed with a strategy similar to the large time increment (LATIN) method. This is a non-incremental (i.e. global in time), and iterative method, well-suited for nonlinear problems. The method generally starts
with an approximation of the quantities of interest for the whole spatial-temporal domain. Thereafter, at each iteration, the quantities of interest are improved for the whole spatial and temporal domain till a convergence is reached. In general LATIN method is based on three basic principles (see [18]).

* The first principle is the separation of difficulties. The governing equations are separated into two groups. A group of equations which are local in space, possibly non-linear; and a group of equations which are linear but possibly global in space.

* The second principle is a two-step iterative algorithm, where at each iteration the solution is constructed alternatively for the first group of equations and then for the second group of equations.

* The last principle is to solve the global problem defined on the entire space-time domain using proper generalised decomposition (PGD) based reduced order approximation such that the numerical cost is drastically lowered.

In the context of MCRE, separation of difficulties is achieved by splitting $\zeta^2$ into two parts:

$$ \zeta_{\psi}^2 = \int_0^T \int_{\Omega} \eta_{\psi} d\Omega dt + \frac{\gamma}{2} \int_0^T \int_{\Omega_m} \|\underline{u} - \tilde{\underline{u}}\|^2 d\Omega_m dt, \quad (11a) $$

$$ \zeta_{\phi}^2 = \int_0^T \int_0^t \int_{\Omega} \eta_{\phi} d\Omega d\tau dt. \quad (11b) $$

This separation enables minimising each term alternatively. The iterative minimisation is performed through a two-step methodology.

### 3.1.1 Global step

This step consists in the minimisation of eq. (11a), representing the state equations residual and measurement discrepancy. This minimisation is linear and global in space, under admissibility constraints. The kinematic admissibility is enforced in the discretised space as essential boundary condition. The static admissibility is enforced through a Lagrangian as

$$ \mathcal{L} = \zeta_{\psi}^2 - \int_0^T \left[ \int_{\Omega} \underline{\sigma} : \underline{\varepsilon}(\underline{\lambda}) d\Omega - \int_{\Omega} \underline{f}_d^v \cdot \underline{\lambda} d\Omega - \int_{\partial_f \Omega} \underline{f}_d^s \cdot \underline{\lambda} dS \right] dt, \quad (12) $$

where $\underline{\varepsilon}(\underline{\lambda})$ is the strain associated to the Lagrange multiplier field $\underline{\lambda}$.

The residual on state equations $\eta_{\psi}$ defined in eq. (6a) requires an inconvenient dual Legendre-Fenchel transforma-

tion (see [14]), which is necessary for the expression of $\zeta_{\psi}^2$. Therefore, a displacement field $\underline{v}$ (defined up to a rigid body motion) is introduced by duality such that

$$ \underline{\sigma}_{\underline{v}} = \frac{\partial \psi}{\partial \underline{\varepsilon}_{\underline{v}}^{e}}, \quad \underline{\sigma}_{\underline{v}} \in \mathcal{S}_{ad}. \quad (13) $$

This therefore defines the variables associated to $\underline{v}$ as $\underline{\kappa}_v = \{\underline{\varepsilon}_{\underline{v}}^e (= \underline{\varepsilon}_{\underline{v}} - \underline{\varepsilon}_{\underline{v}}^p), \underline{\varepsilon}_{\underline{v}}^p, \underline{\sigma}_{\underline{v}}, \underline{X}_{\underline{v}}, \underline{Z}_{\underline{v}}\}$, and variables associated to $\underline{u}$ as $\underline{\kappa}_u = \{\underline{\varepsilon}_{\underline{u}}^e (= \underline{\varepsilon}_{\underline{u}} - \underline{\varepsilon}_{\underline{u}}^p), \underline{\varepsilon}_{\underline{u}}^p, \underline{\sigma}_{\underline{u}}, \underline{X}_{\underline{u}}, \underline{Z}_{\underline{u}}\}$. Equation (12) can be rewritten as

$$ \begin{aligned} \mathcal{L} = & \int_0^T \int_{\Omega} \left[ \psi(\underline{\varepsilon}_{\underline{u}} - \underline{\varepsilon}_{\underline{u}}^p, \underline{X}_{\underline{u}}) - \psi(\underline{\varepsilon}_{\underline{v}} - \underline{\varepsilon}_{\underline{v}}^p, \underline{X}_{\underline{v}}) \right. \\ & \left. - \underline{\sigma}_{\underline{v}} : \left( (\underline{\varepsilon}_{\underline{u}} - \underline{\varepsilon}_{\underline{u}}^p) - (\underline{\varepsilon}_{\underline{v}} - \underline{\varepsilon}_{\underline{v}}^p) \right) - \underline{Z}_{\underline{v}} \cdot (\underline{X}_{\underline{u}} - \underline{X}_{\underline{v}}) \right] d\Omega dt \\ & + \frac{\gamma}{2} \int_0^T \int_{\Omega_m} \|\underline{u} - \tilde{\underline{u}}\|^2 d\Omega_m dt \\ & - \int_0^T \left[ \int_{\Omega} \underline{\sigma}_{\underline{v}} : \underline{\varepsilon}(\underline{\lambda}) d\Omega - \int_{\Omega} \underline{f}_d^v \cdot \underline{\lambda} d\Omega - \int_{\partial_f \Omega} \underline{f}_d^s \cdot \underline{\lambda} dS \right] dt. \end{aligned} \quad (14) $$

Equation (14) is a compromise between the residual on state equations and the discrepancy with measurements. Its minimisation is a linear (because elasticity is assumed to be linear here) and global in space problem. During the global step, internal variables $\underline{X}_{\underline{u}}$, $\underline{X}_{\underline{v}}$, $\underline{Z}_{\underline{u}}$, $\underline{Z}_{\underline{v}}$, $\underline{\varepsilon}_{\underline{u}}^e$, $\underline{\varepsilon}_{\underline{v}}^p$ are frozen to the value obtained at the last local step.

Stationarity of eq. (14), i.e. search of the saddle point, results into the following set of equations

$$ \begin{aligned} & \int_0^T \int_{\Omega} \mathbf{c}(\underline{\varepsilon}_{\underline{v}} - \underline{\varepsilon}_{\underline{v}}^p) : \underline{\varepsilon}(\delta\underline{\lambda}) d\Omega dt - \int_0^T \int_{\Omega} \underline{f}_d^v \cdot \delta\underline{\lambda} d\Omega dt \\ & - \int_0^T \int_{\partial_f \Omega} \underline{f}_d^s \cdot \delta\underline{\lambda} dS dt = 0, \end{aligned} \quad (15a) $$

$$ \int_0^T \int_{\Omega} \mathbf{c} \left[ (\underline{\varepsilon}_{\underline{v}} - \underline{\varepsilon}_{\underline{v}}^p) - (\underline{\varepsilon}_{\underline{u}} - \underline{\varepsilon}_{\underline{u}}^p) - \underline{\varepsilon}_{\underline{\lambda}} \right] : \underline{\varepsilon}(\delta\underline{v}) d\Omega dt = 0, \quad (15b) $$

$$ \begin{aligned} & \int_0^T \int_{\Omega} \mathbf{c} \left[ (\underline{\varepsilon}_{\underline{u}} - \underline{\varepsilon}_{\underline{u}}^p) - (\underline{\varepsilon}_{\underline{v}} - \underline{\varepsilon}_{\underline{v}}^p) \right] : \underline{\varepsilon}(\delta\underline{u}) d\Omega dt \\ & + \gamma \int_0^T \int_{\Omega_m} (\underline{u} - \tilde{\underline{u}}) \cdot \delta\underline{u} d\Omega_m dt = 0. \end{aligned} \quad (15c) $$

Here, the free energy function $\psi$ is considered to be

$$ \begin{aligned} & \psi(\underline{\varepsilon} - \underline{\varepsilon}_{\underline{v}}^p, \underline{X}_{\underline{v}}) = \psi_{\underline{\varepsilon}}(\underline{\varepsilon} - \underline{\varepsilon}_{\underline{v}}^p) + \psi_{\underline{X}}(\underline{X}_{\underline{v}}), \\ & \text{with } \psi_{\underline{\varepsilon}} = \frac{1}{2}(\underline{\varepsilon} - \underline{\varepsilon}_{\underline{v}}^p) : \mathbf{C}(\underline{\varepsilon} - \underline{\varepsilon}_{\underline{v}}^p), \end{aligned} \quad (16) $$

where $\mathbf{C}$ is the Hooke tensor.
Consider a particular LATIN iteration j; the following quantities are represented in terms of corrections, i.e.

$$
\begin{align}
\Delta \underline{u}_j &= \underline{u}_j - \underline{u}_{j-1}, \quad \Delta \underline{v}_j = \underline{v}_j - \underline{v}_{j-1}, \quad \Delta \underline{\lambda}_j = \underline{\lambda}_j - \underline{\lambda}_{j-1}, \nonumber \\
\Delta \underline{\varepsilon}_{\underline{u} j} &= \underline{\varepsilon}_{\underline{u} j} - \underline{\varepsilon}_{\underline{u} j-1}, \quad \Delta \underline{\varepsilon}_{\underline{v} j} = \underline{\varepsilon}_{\underline{v} j} - \underline{\varepsilon}_{\underline{v} j-1}, \nonumber \\
\Delta \underline{\varepsilon}_{\underline{\lambda} j} &= \underline{\varepsilon}_{\underline{\lambda} j} - \underline{\varepsilon}_{\underline{\lambda} j-1}, \nonumber \\
\Delta \underline{\varepsilon}_{\underline{u} j}^p &= \underline{\varepsilon}_{\underline{u} j}^p - \underline{\varepsilon}_{\underline{u} j-1}^p, \quad \Delta \underline{\varepsilon}_{\underline{v} j}^p = \underline{\varepsilon}_{\underline{v} j}^p - \underline{\varepsilon}_{\underline{v} j-1}^p. \tag{17}
\end{align}
$$

Using these corrective terms, eq. (15b) can be rewritten as

$$
\begin{align}
& \int_0^T \int_{\Omega} \mathbf{C} \left[ \begin{aligned}[t]
    & \left( \underline{\varepsilon}_{\underline{v} j-1} - \underline{\varepsilon}_{\underline{v} j-1}^p \right) \\
    & - \left( \underline{\varepsilon}_{\underline{u} j-1} - \underline{\varepsilon}_{\underline{u} j-1}^p \right) - \underline{\varepsilon}_{\underline{\lambda} j-1} \end{aligned} \right] : \underline{\varepsilon} \left( \delta \underline{v} \right) d\Omega dt \nonumber \\
& + \int_0^T \int_{\Omega} \mathbf{C} \left[ \begin{aligned}[t]
    & \left( \Delta \underline{\varepsilon}_{\underline{v} j} - \Delta \underline{\varepsilon}_{\underline{v} j}^p \right) - \left( \Delta \underline{\varepsilon}_{\underline{u} j} - \Delta \underline{\varepsilon}_{\underline{u} j}^p \right) \\
    & - \Delta \underline{\varepsilon}_{\underline{\lambda} j} \end{aligned} \right] : \underline{\varepsilon} \left( \delta \underline{v} \right) d\Omega dt = 0. \tag{18}
\end{align}
$$

Now, considering at LATIN iteration $j-1$, eq. (15b) becomes

$$
\int_0^T \int_{\Omega} \mathbf{C} \left[ \begin{aligned} & \left( \underline{\varepsilon}_{\underline{v} j-1} - \underline{\varepsilon}_{\underline{v} j-1}^p \right) \\ & - \left( \underline{\varepsilon}_{\underline{u} j-1} - \underline{\varepsilon}_{\underline{u} j-1}^p \right) - \underline{\varepsilon}_{\underline{\lambda} j-1} \end{aligned} \right] : \underline{\varepsilon} \left( \delta \underline{v} \right) d\Omega dt = 0. \quad (19)
$$

Using eq. (19), eq. (18) yields

$$ (\Delta \underline{\varepsilon}_{\underline{v} j} - \Delta \underline{\varepsilon}_{\underline{v} j}^p) - (\Delta \underline{\varepsilon}_{\underline{u} j} - \Delta \underline{\varepsilon}_{\underline{u} j}^p) - \Delta \underline{\varepsilon}_{\underline{\lambda} j} = 0, \quad (20) $$

up to a rigid body motion taken here as 0. Using eq. (20), eq. (15) can be rewritten for LATIN iteration j as

$$
\begin{align}
& \int_0^T \int_{\Omega} \mathbf{C} \left[ \begin{aligned}[t]
    & \left( \underline{\varepsilon}_{\underline{v} j-1} - \underline{\varepsilon}_{\underline{v} j-1}^p \right) + \left( \Delta \underline{\varepsilon}_{\underline{u} j} - \Delta \underline{\varepsilon}_{\underline{u} j}^p \right) + \Delta \underline{\varepsilon}_{\underline{\lambda} j} \end{aligned} \right] \nonumber \\
& : \underline{\varepsilon} \left( \delta \underline{\lambda} \right) d\Omega dt \nonumber \\
& - \int_0^T \int_{\Omega} \int_{\partial_f \Omega} \int_{\underline{f}_d}^s \delta \underline{\lambda} d\Omega dt - \int_0^T \int_{\partial_f \Omega} \int_{\underline{f}_d}^s \delta \underline{\lambda} dS dt = 0, \tag{21a}
\end{align}
$$

$$
\begin{align}
& \int_0^T \int_{\Omega} \mathbf{C} \left[ \begin{aligned}[t]
    & \left( \underline{\varepsilon}_{\underline{u} j-1} - \underline{\varepsilon}_{\underline{u} j-1}^p \right) - \left( \underline{\varepsilon}_{\underline{v} j-1} - \underline{\varepsilon}_{\underline{v} j-1}^p \right) - \Delta \underline{\varepsilon}_{\underline{\lambda} j} \end{aligned} \right] \nonumber \\
& : \underline{\varepsilon} \left( \delta \underline{u} \right) d\Omega dt \nonumber \\
& + \gamma \int_0^T \int_{\Omega_m} \left( \underline{u}_{j-1} + \Delta \underline{u}_j - \tilde{\underline{u}} \right) \cdot \delta \underline{u} d\Omega_m dt = 0. \tag{21b}
\end{align}
$$

We now introduce a discretisation through finite elements such that

$$ \underline{u} = \mathbb{N}u, \quad \varepsilon = \mathbb{B}u, \quad \mathbb{K} = \left[ \int_{\Omega} \mathbb{B}^T \mathbf{C} \mathbb{B} d\Omega \right], \quad \mathbb{L} $$

$$ = \left[ \int_{\Omega} \mathbb{B}^T \mathbf{C} d\Omega \right], \qquad (22) $$

where $\mathbb{N}$ is the shape function matrix, $\mathbb{B}$ involves derivatives of the shape functions, $\mathbb{K}$ is the stiffness matrix, and $\mathbb{L}$ is a similar matrix operator that projects through integration the strains at Gauss points to the nodes. The time discretisation is such that $\{t_0 = 0, \dots, t_f = T\}$. The discretised system of equations for a particular time step $t_l$ becomes

$$ \left[ \begin{array}{cc} \mathbb{K}_{aa} & \mathbb{K}_{aa} \\ (\gamma \Pi^T \Pi)_{aa} & -\mathbb{K}_{aa} \end{array} \right] \left[ \begin{array}{c} \Delta \mathbf{u}_{a,l,j} \\ \Delta \lambda_{a,l,j} \end{array} \right] = \left[ \begin{array}{c} \mathbb{P} \\ \mathbb{Q} \end{array} \right], \qquad (23) $$

with

$$ \mathbb{P} = \mathbb{L}_{ag} \Delta \varepsilon_{ul,j}^p + \mathbb{F}_{a,l} - \mathbb{L}_{ag} \left( \varepsilon_{vl,j-1} - \varepsilon_{vl,j-1}^p \right) \quad (24a) $$

$$
\begin{align}
\mathbb{Q} ={}& \gamma \left( (\Pi^T)_{ao} \tilde{\mathbf{u}}_l - (\Pi^T \Pi)_{ao} \mathbf{u}_{l,j-1} \right) \nonumber \\
& - \mathbb{L}_{ag} \left( \varepsilon_{ul,j-1} - \varepsilon_{ul,j-1}^p - \varepsilon_{vl,j-1} + \varepsilon_{vl,j-1}^p \right). \tag{24b}
\end{align}
$$

$\mathbf{u}$ is the discretised form of $\underline{u}$, i.e the vector of nodal unknowns. $\varepsilon$ is the discretised form of $\underline{\varepsilon}$ (similarly, $\varepsilon^p$ is that of $\underline{\varepsilon}^p$), described at the Gauss points. $\mathbb{F}$ corresponds to the nodal forces, and $\lambda$ is the discretised version of $\underline{\lambda}$, described at the nodes.

For the initialisation of the LATIN iterations, i.e. $j = 0$, it is usual to consider an elastic problem. The discretised set of equations for the elastic initialisation reads

$$ \left[ \begin{array}{cc} \mathbb{K}_{aa} & \mathbb{K}_{aa} \\ (\gamma \Pi^T \Pi)_{aa} & -\mathbb{K}_{aa} \end{array} \right] \left[ \begin{array}{c} \mathbf{u}_{a,0} \\ \lambda_{a,0} \end{array} \right] = \left[ \begin{array}{c} \mathbb{F}_a - \mathbb{K}_{ad} \mathbf{U}_d \\ \gamma \left( (\Pi^T)_{ao} \tilde{\mathbf{u}} - (\Pi^T \Pi)_{ad} \mathbf{U}_d \right) \end{array} \right]. \quad (25) $$

Here the subscript $d$ denotes the nodes where the displacement is prescribed as $\mathbf{U}_d$, such that $\mathbf{u}_j = \{\mathbf{U}_d, \mathbf{u}_{a,j}\}^T$, where the subscript $a$ represents the unprescribed nodes. For the Lagrange multiplier field, the delineation is $\lambda_j = \{0, \lambda_{a,j}\}^T$. The corrective terms are demarcated as $\Delta \mathbf{u}_j = \{0, \Delta \mathbf{u}_{a,j}\}^T$ and $\Delta \lambda_j = \{0, \Delta \lambda_{a,j}\}^T$. Similarly, the force vector is delimited as $\mathbb{F} = \{\mathbb{F}_d, \mathbb{F}_a\}^T$, where $\mathbb{F}_a$ can contain both free edge and applied force. The operator $\Pi$ essentially projects the nodal unknowns to the measurement grid (see [12] for details), and the subscript $g$ represents all the Gauss points.

Once $\mathbf{u}_0$ and $\lambda_0$ are obtained, $\mathbf{v}_0$ is obtained as

$$ \mathbf{v}_0 = \mathbf{u}_0 + \lambda_0, \qquad (26) $$

where $\mathbf{v}$ is the discretised form of $\underline{v}$.

Similarly, the resolution of $\Delta \mathbf{u}_j$, $\Delta \lambda_j$ from eq. (23) results in the computation of $\Delta \varepsilon_{u j}$ and $\Delta \varepsilon_{\lambda j}$, and thereby $\Delta \varepsilon_{v j}$ is
calculated through eq. (20). Δv<sub>a,j</sub> is thereby computed as

$$
[\mathbb{K}_{aa}] \{ \Delta v_{a,j} \} = \{ \mathbb{K}_{aa} \Delta u_{a,j} + \mathbb{K}_{aa} \Delta \lambda_{a,j} \\
+ \mathbb{L}_{ag} (\Delta \varepsilon_{v,j}^p - \Delta \varepsilon_{u,j}^p) \}, \quad (27)
$$

defining Δv<sub>j</sub> = {0, Δv<sub>a,j</sub>}<sup>T</sup>. u<sub>j</sub>, v<sub>j</sub>, λ<sub>j</sub>, ε<sub>u,j</sub>, ε<sub>v,j</sub>, ε<sub>λ,j</sub> are
thereby calculated from eq. (17).

3.1.2 Local step

In this step, the minimisation of ζ<sub>φ</sub><sup>2</sup> (eq. (11b)) corresponds
to the integration of evolution laws at each Gauss point. The
integration of the evolution laws is performed with an Euler
scheme with the total strains ε<sub>≡u</sub> and ε<sub>≡v</sub> frozen to the values
obtained from the last global step. The initial conditions on
internal variables are enforced for the first time step. In the
rate-independent case, the dissipation potential is the indica-
tor of a convex domain and is not differentiable, whereas in
the rate-dependent case, this potential is differentiable.

The solution strategy is similar to an elastic predictor-
(visco)plastic corrector type methodology, with the total
strains (ε_u, ε_v) obtained from the last global stage. The initial
condition has to be imposed in this stage.

Consider, for instance, a particular time step t<sub>l</sub> and a par-
ticular LATIN iteration j. For any particular Gauss point,
with w = {u, v}, the elastic prediction involves initialising
the variables as

$$
\epsilon_{wl,j}^{p} = \epsilon_{wl-1,j}^{p}, \quad \chi_{wl,j} = \chi_{wl-1,j}, \quad (28a)
$$

$$
Z_{wl,j} = \frac{\partial \psi_X}{\partial \underline{X}} |_{X_{wl,j}}, \sigma_{wl,j} = \frac{\partial \psi_{\varepsilon}}{\partial (\underline{\varepsilon} - \underline{\varepsilon}^p)} |_{\varepsilon_{wl,j-1}, \varepsilon_{wl,j}^p}. \quad (28b)
$$

The next step is to check the yield criterion. If the yield
function f<sub>w,y</sub> ≥ 0, then the variables are updated by solving
the following set of equations,

$$
\sigma_{wl,j} = \frac{\partial \psi_{\epsilon}}{\partial \left(\underline{\epsilon} - \underline{\epsilon}^p\right)} |_{\epsilon_{wl,j-1}, \epsilon_{wl,j}^p} Z_{wl,j} = \frac{\partial \psi_X}{\partial \underline{X}} |_{X_{wl,j}}, \quad (29a)
$$

$$
\begin{align}
\epsilon_{wl,j}^p &= \epsilon_{wl-1,j}^p + \Delta t_l \frac{\partial \phi^*}{\partial \underline{\sigma}} |_{\sigma_{wl,j}, Z_{wl,j}} \chi_{wl,j} \nonumber \\
&= \chi_{wl-1,j} - \Delta t_l \frac{\partial \phi^*}{\partial \underline{Z}} |_{\sigma_{wl,j}, Z_{wl,j}}, \tag{29b}
\end{align}
$$

with the time step size Δt<sub>l</sub> = t<sub>l</sub> − t<sub>l−1</sub>. Here the Euler implicit
method is used to integrate the first order ODEs, however
other methods can also be used. The nonlinear equations (eq.
(29)) are to be solved simultaneously using any nonlinear

solver. Here X, Z, σ are the discretised versions of X, Z, σ,
defined at the Gauss points.

It has to be mentioned that for the initialisation step
(j = 0), the internal variables ε<sup>p</sup><sub>u0</sub>, ε<sup>p</sup><sub>v0</sub>, X<sub>u0</sub>, X<sub>v0</sub>, Z<sub>u0</sub>, Z<sub>v0</sub>
are considered to be 0, and the stress variables are obtained
as

$$
\sigma_{u0} = \frac{\partial \psi_{\epsilon}}{\partial (\underline{\epsilon} - \underline{\epsilon}^p)} |_{\epsilon_{u0}, \epsilon_{u0}^p}, \quad \sigma_{v0} = \frac{\partial \psi_{\epsilon}}{\partial (\underline{\epsilon} - \underline{\epsilon}^p)} |_{\epsilon_{v0}, \epsilon_{v0}^p}. \qquad (30)
$$

The stopping criterion (for the iterative process) is essen-
tially defined through a LATIN error indicator ℓ obtained
as

$$
\ell = \left[ \frac{\int_0^T \int_{\Omega} \left( \left( \underline{\epsilon}_{\equiv v j} - \underline{\epsilon}_{\equiv v j-1} \right) : \left( \underline{\epsilon}_{\equiv v j} - \underline{\epsilon}_{\equiv v j-1} \right) \right) d\Omega dt}{\frac{1}{2} \int_0^T \int_{\Omega} \left( \underline{\epsilon}_{\equiv v j} : \underline{\epsilon}_{\equiv v j} + \underline{\epsilon}_{\equiv v j-1} : \underline{\epsilon}_{\equiv v j-1} \right) d\Omega dt} + \frac{\int_0^T \int_{\Omega} \left( \underline{\epsilon}_{\equiv u j} - \underline{\epsilon}_{\equiv u j-1} \right) : \left( \underline{\epsilon}_{\equiv u j} - \underline{\epsilon}_{\equiv u j-1} \right) d\Omega dt}{\frac{1}{2} \int_0^T \int_{\Omega} \left( \underline{\epsilon}_{\equiv u j} : \underline{\epsilon}_{\equiv u j} + \underline{\epsilon}_{\equiv u j-1} : \underline{\epsilon}_{\equiv u j-1} \right) d\Omega dt} \right]^{1/2}, \quad (31)
$$

and if this is less than a certain pre-determined value, the
LATIN iteration is considered to be converged.

3.2 The identification problem

In the proposed approach, the identification problem of the
inverse problem is the minimisation of the MCRE functional
(ζ²) with respect to parameters. In order to solve this minimi-
sation problem, the various optimisation methods available
in the literature can be used. A convenient approach is to use
the optimal admissible set (as given in section 3.1), thereby
it is natural to solve the minimisation problem eq. (10) using
the gradient descent approach. The update of the parameters
is given as

$$
\underline{\theta}_i = \underline{\theta}_{i-1} - \hbar \frac{d\zeta^2(x_i, \underline{\theta})}{d\underline{\theta}}. \quad (32)
$$

Here, ħ is the step-size of the gradient descent. The gradient
is computed through the adjoint-state approach and eq. (32)
can be rewritten as

$$
\underline{\theta}_i = \underline{\theta}_{i-1} - \hbar \frac{\partial \zeta^2(x_i, \underline{\theta})}{\partial \underline{\theta}}. \quad (33)
$$

The succinct procedure of the total inverse problem is given
in algorithm 1.
**Algorithm 1 MCRE method for inverse analysis**

**Initialise:** Parameter set $\theta_0$, iteration $i = 0$

**while** $\zeta^2 \ge \epsilon_{tol}$ **do**

Update $i = i + 1$

Use $\theta_{j-1}$ for basic problem

**Elastic initialisation:**

Initialise $j = 0$ and $\ell$

Solve for $u_0, u_0, \lambda_0$

Compute $\epsilon_{u0}, \epsilon_{v0}, \epsilon_{\lambda0}, \sigma_{u0}, \sigma_{v0}$

Set $\epsilon_{u0}^p, \epsilon_{v0}^p, X_{u0}, X_{v0}, Z_{u0}, Z_{v0} = 0$

**while** $\ell \ge \ell_{tol}$ **do**

**Local step:**

Update $j = j + 1$

Use $\epsilon_{uj-1}, \epsilon_{vj-1}$

Solve for $\epsilon_{uj}^p, \epsilon_{vj}^p, X_{uj}, X_{vj}, Z_{uj}, Z_{vj}, \sigma_{uj}, \sigma_{vj}$ through eqs.

(28) and (29)

**Global step:**

Represent $\Delta u_j, \Delta v_j, \Delta \lambda_j, \Delta \epsilon_{uj}, \Delta \epsilon_{vj}, \Delta \epsilon_{\lambda j}, \Delta \epsilon_{uj}^p, \Delta \epsilon_{vj}^p$ using eq. (17)

Solve for $\Delta u_j, \Delta v_j, \Delta \lambda_j$ using eqs. (25) and (27)

Compute $\Delta \epsilon_{uj}, \Delta \epsilon_{vj}, \Delta \epsilon_{\lambda j}$

Compute $u_j, v_j, \lambda_j, \epsilon_{uj}, \epsilon_{vj}, \epsilon_{\lambda j}$ using eq. (17)

Calculate $\ell$ from eq. (31)

**end while**

Output $x_i = \{ u_j, v_j, \epsilon_{uj}, \epsilon_{vj}, \epsilon_{uj}^p, \sigma_{uj}, \sigma_{vj}, \epsilon_{vj}^p, X_{uj}, X_{vj}, Z_{uj}, Z_{vj} \}$

Compute $\zeta^2(x_i, \theta_{i-1}), \frac{\partial \zeta^2(x_i, \theta)}{\partial \theta}|_{\theta_{i-1}}$

Update $\theta_j$ from eq. (33)

**end while**

## 4 PGD for the basic problem

The main purpose of reduced order modelling (ROM) is to provide a frugal resolution technique. The idea of model reduction is to seek the solution in a reduced-order basis, whose dimension is much smaller than the size of the original high-dimensional model. As mentioned before, the basic problem is the most expensive part of the inverse problem, thereby the idea herein is to propose a reduced order methodology that reduces its cost.

Proper generalised decomposition (PGD) 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. PGD is a quite obvious choice among reduced order approaches as it is extremely conducive to LATIN-type approach and does not require any separate training phase for generating the reduced order basis.

In the context of the inverse problem described in section 3, the idea is to separate the space-time problem into segregated space and time problems. This is essentially done by representing the quantities of interest into separated variable forms (as functions of space and time).

Consider first the global stage of the basic problem, especially the quantities of interest $\Delta \lambda_j$ and $\Delta u_j$, for a particular LATIN iteration $j$. Usually there are two options, either add

new reduced order bases (i.e add new space-time pairs), or update the existing reduced order bases (i.e update only the time functions).

### 4.1 Addition of space-time modes

Considering that the variables $\lambda_{j-1}$ and $u_{j-1}$ are approximated as

$$ \begin{gathered} \underline{\lambda}_{j-1} = \underline{\lambda}_0 + \sum_{l=1}^{\omega-1} \overline{\Lambda}^l(\Omega) \alpha_{\lambda}^l(t), \\ \underline{u}_{j-1} = \underline{u}_0 + \sum_{j=1}^{\zeta-1} \overline{U}^j(\Omega) \alpha_{u}^j(t), \end{gathered} \quad (34) $$

the variables $\Delta \lambda_j$ and $\Delta u_j$ are represented as

$$ \Delta \underline{\lambda}_j (\Omega, t) = \overline{\Lambda}^{\omega} (\Omega) \alpha_{\lambda}^{\omega}(t), \quad \Delta \underline{u}_j (\Omega, t) = \overline{U}^{\zeta} (\Omega) \alpha_{u}^{\zeta}(t), \quad (35) $$

where $\omega$ and $\zeta$ are the numbers of modes. ($\overline{\Lambda}, \overline{U}$) represent space functions, while $(\alpha_{\lambda}, \alpha_{u})$ represent time functions. It has to be mentioned that according to eq. (35), only one new space-time is being added per quantity of interest for a given LATIN iteration $j$ (i.e. $\omega$th mode for $\Delta \lambda_j$ and $\zeta$th mode for $\Delta u_j$). The subscript $j$ is dropped for $\Delta \lambda_j$ and $\Delta u_j$, and the subscripts $\omega$, $\zeta$ are dropped for $\overline{\Lambda}, \overline{U}, \alpha_{\lambda}, \alpha_{u}$, and $\alpha_{\lambda}$ for the subsequent developments, for the sake of clarity.

Considering now eq. (21), $\delta \underline{\lambda}$ and $\delta \underline{u}$ are represented through separated variable form as

$$ \delta \underline{\lambda} = \overline{\Lambda} \delta \alpha_{\lambda} + \alpha_{\lambda} \delta \overline{\Lambda}, \quad \delta \underline{u} = \overline{U} \delta \alpha_{u} + \alpha_{u} \delta \overline{U}. \quad (36) $$

Also, the corresponding strains are written as

$$ \Delta \underline{\epsilon}_{\equiv \lambda} (\Omega, t) = \overline{\underline{\epsilon}}_{\equiv \lambda} (\Omega) \alpha_{\lambda}(t), \quad \Delta \underline{\epsilon}_{\equiv u} (\Omega, t) = \overline{\underline{\epsilon}}_{\equiv u} (\Omega) \alpha_{u}(t), \quad (37) $$

where $\overline{\underline{\epsilon}}_{\equiv \lambda}, \underline{\epsilon}_{\equiv u}$ are space functions.

Using eq. (36), eq. (21) can be re-written as a space problem of the form

$$ \begin{aligned} & \int_{\Omega} \mathbf{C} \left[ \overline{\underline{\epsilon}}_{\equiv u} \beta + \overline{\underline{\epsilon}}_{\equiv \lambda} \kappa \right] : \underline{\epsilon} (\delta \overline{\Lambda}) d\Omega \\ & + \int_{\Omega} \mathbf{C} \int_{0}^{T} \left( (\underline{\epsilon}_{v}^{p} j - \underline{\epsilon}_{v}^{p} j - \Delta \underline{\epsilon}_{u}^{p} j) \alpha_{\lambda} \right) dt : \underline{\epsilon} (\delta \overline{\Lambda}) d\Omega \\ & - \int_{\Omega} \int_{0}^{T} (\underline{f}_{d}^{v} \alpha_{\lambda}) dt \cdot \delta \overline{\Lambda} d\Omega \\ & - \int_{\partial f \Omega} \int_{0}^{T} (\underline{f}_{d}^{s} \alpha_{\lambda}) dt \cdot \delta \overline{\Lambda} dS = 0, \end{aligned} \quad (38a) \\ \int_{\Omega_m} (\overline{U} \eta) \cdot \delta \overline{U} d\Omega_m - \int_{\Omega} \mathbf{C} (\underline{\epsilon}_{\equiv \lambda}) : \underline{\epsilon} (\delta \overline{U}) d\Omega \end{aligned} $$
$$
\begin{equation}
\begin{aligned}
& + \int_{\Omega} \mathbf{c} \int_{0}^{T} \left( \left( \begin{aligned}[t]
    & \underline{\epsilon}_{u, j-1} - \underline{\epsilon}_{u, j-1}^{p} \\
    & - \underline{\epsilon}_{v, j-1} + \underline{\epsilon}_{v, j-1}^{p} \end{aligned} \right) \alpha_{u} \right) dt : \underline{\epsilon} (\delta \bar{U}) d\Omega \\
& + \gamma \int_{\Omega_m} \int_{0}^{T} \left( \left( \underline{u}_{j-1} - \tilde{u} \right) \alpha_{u} \right) dt \cdot \delta \bar{U} d\Omega_m = 0,
\end{aligned}
\tag{38b}
\end{equation}
$$

with

$$
\begin{align}
\beta &= \int_{0}^{T} (\alpha_{u} \alpha_{\lambda}) dt, \quad \kappa = \int_{0}^{T} (\alpha_{\lambda} \alpha_{\lambda}) dt, \nonumber \\
\eta &= \gamma \int_{0}^{T} (\alpha_{u} \alpha_{u}) dt, \tag{39}
\end{align}
$$

and a time problem,

$$
\begin{equation}
\begin{aligned}
& \int_0^T [\alpha_u \Upsilon + \alpha_\lambda \Gamma] \delta \alpha_\lambda dt \\
& \quad + \int_0^T \int_\Omega \left( \mathbf{C} \left( \underline{\epsilon}_{\equiv u, j-1} - \underline{\epsilon}_{\equiv v, j-1}^p - \Delta \underline{\epsilon}_{\equiv u, j}^p \right) : \underline{\bar{\epsilon}}_{\equiv \lambda} \right) d\Omega \delta \alpha_\lambda dt \\
& \quad - \int_0^T \int_\Omega \left( \underline{f}_d^v \cdot \underline{\bar{\Lambda}} \right) d\Omega \delta \alpha_\lambda dt \\
& \quad - \int_0^T \int_{\partial_f \Omega} \left( \underline{f}_d^s \cdot \underline{\bar{\Lambda}} \right) dS \delta \alpha_\lambda dt = 0,
\end{aligned}
\tag{40a}
\end{equation}
$$

$$
\begin{equation}
\begin{aligned}
& \int_0^T (\alpha_u \Psi) \delta \alpha_u dt - \int_0^T (\alpha_\lambda \Phi) \delta \alpha_u dt \\
& \quad + \int_0^T \int_\Omega \left( \mathbf{C} \left( \underline{\epsilon}_{\equiv u, j-1} - \underline{\epsilon}_{\equiv u, j-1}^p \right. \right. \\
& \qquad \left. \left. - \underline{\epsilon}_{\equiv v, j-1} + \underline{\epsilon}_{\equiv v, j-1}^p \right) : \underline{\bar{\epsilon}}_{\equiv u} \right) d\Omega \delta \alpha_u dt \\
& \quad + \gamma \int_0^T \int_{\Omega_m} \left( \left( \underline{u}_{j-1} - \underline{\tilde{u}} \right) \cdot \underline{\bar{U}} \right) d\Omega_m \delta \alpha_u dt = 0,
\end{aligned}
\tag{40b}
\end{equation}
$$

with

$$
\begin{align}
\Upsilon &= \int_{\Omega} \left( \mathbf{C} \underline{\underline{\epsilon}}_{\equiv u} : \underline{\underline{\epsilon}}_{\equiv \lambda} \right) d\Omega, \quad \Gamma = \int_{\Omega} \left( \mathbf{C} \underline{\underline{\epsilon}}_{\equiv \lambda} : \underline{\underline{\epsilon}}_{\equiv \lambda} \right) d\Omega, \nonumber \\
\Psi &= \gamma \int_{\Omega_m} \left( \underline{\underline{U}} \cdot \underline{\underline{U}} \right) d\Omega_m, \quad \Phi = \int_{\Omega} \left( \mathbf{C} \underline{\underline{\epsilon}}_{\equiv \lambda} : \underline{\underline{\epsilon}}_{\equiv u} \right) d\Omega. \tag{41}
\end{align}
$$

Using the finite element discretisation, the space problem
reads

$$
\begin{bmatrix} \beta \mathbb{K}_{aa} & \kappa \mathbb{K}_{aa} \\ \eta (\Pi^T \Pi)_{aa} & -\beta \mathbb{K}_{aa} \end{bmatrix} \begin{bmatrix} \bar{U}_a \\ \bar{\Lambda}_a \end{bmatrix} = \begin{bmatrix} M_a \\ N_a \end{bmatrix}, \quad (42)
$$

with

$$
M = \left[ \left( \mathbb{L} \Delta \epsilon_{ul,j}^p + F_l - \mathbb{L} \left( \epsilon_{vl,j-1} - \epsilon_{vl,j-1}^p \right) \right) \alpha_{\lambda l} \right], \quad (43a)
$$

$$
N = \int (\gamma (\Pi^T \tilde{u}_l - \Pi^T \Pi u_{l,j-1}))
$$

$$
-\mathbb{L} \left( \epsilon_{ul,j-1} - \epsilon_{ul,j-1}^p - \epsilon_{vl,j-1} + \epsilon_{vl,j-1}^p \right) \alpha_{ul} \right], \quad (43b)
$$

where $l = \{0, \dots, f\}$, and $\lfloor \rfloor$ represents integrations with respect to time. Equation (39) can also be written in discretised form, similar to eq. (43). These integrations can be solved by any numerical integration method (for instance the trapezoidal rule). Also, the subscript $a$ follows the same definition introduced in section 3.1. $\bar{U}, \bar{\Lambda}$ are the discretised forms of $\underline{U}, \underline{\Lambda}$ at the nodes, and $\alpha_u, \alpha_\lambda$ are the discretised forms of $\alpha_u, \alpha_\lambda$ at the time steps.

The time problem can similarly be written for a particular
time step $t_l$ as

$$
\begin{bmatrix} \Upsilon & \Gamma \\ \Psi & -\Phi \end{bmatrix} \begin{bmatrix} \alpha_{ul} \\ \alpha_{\lambda l} \end{bmatrix} = \begin{bmatrix} G_l \\ H_l \end{bmatrix}, \quad (44)
$$

with

$$
\begin{equation}
\begin{split}
\Upsilon &= \bar{\Lambda}^T \mathbb{K} \bar{U}, \quad \Gamma = \bar{\Lambda}^T \mathbb{K} \bar{\Lambda}, \quad \Psi = \gamma \bar{U}^T \Pi^T \Pi \bar{U}, \quad \Phi = \bar{U}^T \mathbb{K} \bar{\Lambda}, \\
\end{split}
\tag{45}
\end{equation}
$$

and

$$
G_l = \bar{\Lambda}^T \left( \mathbb{L} \Delta \epsilon_{ul,j}^p + F_l - \mathbb{L} \left( \epsilon_{vl,j-1} - \epsilon_{vl,j-1}^p \right) \right), \quad (46a)
$$

$$
\begin{equation}
\begin{split}
H_l = \bar{U}^T & (\gamma (\Pi^T \tilde{u}_l - \Pi^T \Pi u_{l,j-1}) \\
& - \mathbb{L} (\epsilon_{ul,j-1} - \epsilon_{ul,j-1}^p - \epsilon_{vl,j-1} + \epsilon_{vl,j-1}^p)).
\end{split}
\tag{46b}
\end{equation}
$$

The space problem (eq. (42)) and the time problem (eq. (44))
are solved through a fixed point iteration. The iteration is
initialised by arbitrary time functions and thereafter eqs. (42)
and (44) are solved iteratively till convergence is achieved.

The new space functions are orthonormalised with respect
to the previously existing spatial bases through the Gram-
Schmidt method, in order to prevent over expansion of the
number of basis functions with low contribution for some of
them, that may lead to a useless increase of the computation
cost. Also, in the numerical process all former time func-
tions are updated. The new time functions are thus modified
and the newly added space-time pair may be rejected if the
corresponding modified time function has an insignificant
norm. Interested readers may consult [27] for details on this
process.

4.2 Update of time functions

Other than expanding the reduced order bases, the second
option is to update the existing bases. In that case, considering
the variables $\lambda_{j-1}$ and $u_{j-1}$ described according to eq. (34),
the corrective terms defined in eq. (35) are approximated as

$$
\begin{align}
\Delta \underline{\lambda}_{j-1} &= \sum_{l=1}^{\varpi-1} \underline{\Lambda}^l (\Omega) \Delta \alpha_{\lambda}^l (t), \\
\Delta \underline{u}_{j-1} &= \sum_{j=1}^{\zeta-1} \underline{U}^j (\Omega) \Delta \alpha_{u}^j (t). \tag{47}
\end{align}
$$

Here, {$\Delta\alpha_{\lambda}^{l}$}$_{l=1}^{\varpi-1}$, {$\Delta\alpha_{u}^{l}$}$_{j=1}^{\zeta-1}$ are corrections to the time functions, such that

$$
\alpha_{\lambda}^{l \text{ up}} = \alpha_{\lambda}^{l} + \Delta\alpha_{\lambda}^{l}, \text{ with } l = \{1, \dots, \varpi - 1\}, \quad (48a)
$$

$$
\alpha_{u}^{j \text{ up}} = \alpha_{u}^{j} + \Delta\alpha_{u}^{j}, \text{ with } j = \{1, \dots, \zeta - 1\}. \quad (48\text{b})
$$

The updated functions are represented by {$\alpha_{\lambda, up}^{l}$}$_{l=1}^{\varpi-1}$,
{$\alpha_{u, up}^{j}$}$_{j=1}^{\zeta-1}$, and the total quantities of interest become

$$
\underline{\lambda}_j = \underline{\lambda}_0 + \sum_{l=1}^{\varpi-1} \underline{\Lambda}^l(\Omega) \alpha_{\lambda}^l(t), \quad \underline{u}_j = \underline{u}_0 + \sum_{j=1}^{\zeta-1} \underline{U}^j(\Omega) \alpha_{u}^j(t), \tag{49}
$$

where the time functions are updated functions as per eq. (48).

The problem to be solved is the time problem with
known space functions with the quantities of interest being
{Δα̇<sub>λ</sub>}<sub>l=1</sub><sup>ω-1</sup>, {Δα̇<sub>u</sub>}<sub>j=1</sub><sup>ζ-1</sup>. The time problem of eq. (40) can
be rewritten as

$$
\begin{align}
& \int_0^T \left[ \delta \Delta \alpha_{\lambda}^1 \dots \delta \Delta \alpha_{\lambda}^{\varpi-1} \right] \nonumber \\
& \left( \begin{bmatrix} \gamma_{1,1} & \dots & \gamma_{1,\zeta-1} \\ \vdots & \ddots & \vdots \\ \gamma_{\varpi-1,1} & \dots & \gamma_{\varpi-1,\zeta-1} \end{bmatrix} \begin{bmatrix} \Delta \alpha_u^1 \\ \vdots \\ \Delta \alpha_u^{\zeta-1} \end{bmatrix} \right. \nonumber \\
& \left. + \begin{bmatrix} \Gamma_{1,1} & \dots & \Gamma_{1,\varpi-1} \\ \vdots & \ddots & \vdots \\ \Gamma_{\varpi-1,1} & \dots & \Gamma_{\varpi-1,\varpi-1} \end{bmatrix} \begin{bmatrix} \Delta \alpha_{\lambda}^1 \\ \vdots \\ \Delta \alpha_{\lambda}^{\varpi-1} \end{bmatrix} \right. \nonumber \\
& \left. + \int_{\Omega} \left( \mathbf{C} \left( \varepsilon_{\equiv v}^{p} j-1 - \varepsilon_{\equiv v}^{p} j-1 - \Delta \varepsilon_{\equiv u}^{p} j \right) : \begin{bmatrix} \underline{\varepsilon}_{\equiv \lambda}^{1} \\ \vdots \\ \underline{\varepsilon}_{\equiv \lambda}^{\varpi-1} \end{bmatrix} \right) d\Omega \right. \nonumber \\
& \left. - \int_{\Omega} \left( \underline{f}_d^v \cdot \begin{bmatrix} \underline{\Lambda}^1 \\ \vdots \\ \underline{\Lambda}^{\varpi-1} \end{bmatrix} \right) d\Omega - \int_{\partial_f \Omega} \left( \underline{f}_d^s \cdot \begin{bmatrix} \underline{\Lambda}^1 \\ \vdots \\ \underline{\Lambda}^{\varpi-1} \end{bmatrix} \right) dS \right) dt = 0, \tag{50a}
\end{align}
$$

$$
\int_0^T \left[ \delta \Delta \alpha_u^1 \dots \delta \Delta \alpha_u^{\zeta-1} \right] \left( \begin{bmatrix} \Psi_{1,1} & \dots & \Psi_{1,\zeta-1} \\ \vdots & \ddots & \vdots \\ \Psi_{\zeta-1,1} & \dots & \Psi_{\zeta-1,\zeta-1} \end{bmatrix} \begin{bmatrix} \Delta \alpha_u^1 \\ \vdots \\ \Delta \alpha_u^{\zeta-1} \end{bmatrix} \right)
$$

$$
- \begin{bmatrix} \Phi_{1,1} & \dots & \Phi_{1,\varpi-1} \\ \vdots & \ddots & \vdots \\ \Phi_{\zeta-1,1} & \dots & \Phi_{\zeta-1,\varpi-1} \end{bmatrix} \begin{bmatrix} \Delta \alpha_{\lambda}^{1} \\ \vdots \\ \Delta \alpha_{\lambda}^{\varpi-1} \end{bmatrix} \\
+ \int_{\Omega} \left( \mathbf{C} \left( \varepsilon_{\equiv u}^{p} j-1 - \varepsilon_{\equiv u}^{p} j-1 - \varepsilon_{\equiv v}^{p} j-1 + \varepsilon_{\equiv v}^{p} j-1 \right) : \begin{bmatrix} \underline{\varepsilon}_{\equiv \lambda}^{1} \\ \vdots \\ \underline{\varepsilon}_{\equiv \lambda}^{\varpi-1} \end{bmatrix} \right) d\Omega \\
+ \gamma \int_{\Omega_m} \left( \underline{u}_{j-1} - \tilde{u} \right) \cdot \begin{bmatrix} \underline{U}_1 \\ \vdots \\ \underline{U}_{\zeta-1} \end{bmatrix} d\Omega_m \right) dt = 0, \quad (50\text{b})
$$

with

$$
\begin{align}
\Upsilon_{l,j} &= \int_{\Omega} \left( \mathbf{C}_{\equiv u}^{\varepsilon^j} : \underline{\varepsilon}_{\equiv \lambda}^{l} \right) d\Omega, \quad \Gamma_{l,l'} = \int_{\Omega} \left( \mathbf{C}_{\equiv \lambda}^{\varepsilon^{l'}} : \underline{\varepsilon}_{\equiv \lambda}^{l} \right) d\Omega, \nonumber \\
\Psi_{j,j'} &= \gamma \int_{\Omega_m} \left( \underline{U}^{j'} \cdot \underline{U}^j \right) d\Omega_m, \nonumber \\
\Phi_{j,l} &= \int_{\Omega} \left( \mathbf{C}_{\equiv \lambda}^{\varepsilon^j} : \underline{\varepsilon}_{\equiv u}^{j} \right) d\Omega, \tag{51}
\end{align}
$$

where the indices are *l*, *l'* = {1, ..., ϖ − 1} and *j*, *j'* =
{1, ..., ς − 1}.

The corresponding discretised equation for a time step t_l
can be written as

$$
\begin{bmatrix} [\Upsilon] & [\Gamma] \\ [\Psi] - [\Phi] \end{bmatrix} \begin{bmatrix} \{\Delta\alpha_u\}_l \\ \{\Delta\alpha_\lambda\}_l \end{bmatrix} = \begin{bmatrix} \{G\}_l \\ \{H\}_l \end{bmatrix}, \quad (52)
$$

with

$$
\begin{align*}
\Upsilon_{l,j} &= \bar{\Lambda}_l^T \mathbb{K} \bar{U}_j, \\
\Gamma_{l,l'} &= \bar{\Lambda}_l^T \mathbb{K} \bar{\Lambda}_{l'} , \quad \Psi_{j,j'} = \gamma \bar{U}_j^T \Pi^T \Pi \bar{U}_{j'} , \quad \Phi_{j,l} = \bar{U}_j^T \mathbb{K} \bar{\Lambda}_{jl} \tag{53}
\end{align*}
$$

and

$$
G_{l,l} = \bar{\Lambda}_l^T \left( \mathbb{L} \Delta \varepsilon_{ul,j}^p + \mathbb{F}_l - \mathbb{L} \left( \varepsilon_{vl,j-1} - \varepsilon_{vl,j-1}^p \right) \right), \quad (54a)
$$

$$
H_{j,l} = \bar{U}_j^T \left( \gamma \left( \Pi^T \tilde{u}_l - \Pi^T \Pi u_{l,j-1} \right) - \mathbb{L} \left( \varepsilon_{ul,j-1} - \varepsilon_{ul,j-1}^p - \varepsilon_{vl,j-1} + \varepsilon_{vl,j-1}^p \right) \right). \quad (54\text{b})
$$

This step is essentially inexpensive as the space functions
are known and only a single time problem has to be solved for
a given LATIN iteration *j*. However, just updating the time
functions is usually not enough for the reduced order approx-
imation. Therefore, sometimes it is necessary to expand the
reduced order basis according to section 4.1. For a particular
iteration *j*, after the time functions are updated, a criterion
*ρ* is calculated as

$$
\varphi_j = \frac{\|e_j - e_{j-1}\|}{\frac{1}{2} (\|e_j\| + \|e_{j-1}\|)} \text{, with}
$$
$$e_j = \left[ \frac{\int_0^T \int_{\Omega} \left( (\underline{\sigma}_v - \underline{\sigma}_u)_j : (\underline{\sigma}_{\bar{v}j} - \underline{\sigma}_{\bar{u}j}) \right) d\Omega dt}{\frac{1}{2} \int_0^T \int_{\Omega} (\underline{\sigma}_{\bar{v}j} : \underline{\sigma}_{\bar{v}j} + \underline{\sigma}_{\bar{u}j} : \underline{\sigma}_{\bar{u}j}) d\Omega dt} \right]^{1/2} . \quad (55)$$

If $\wp_j$ is less than a predetermined value, new space-time modes are added for the next iteration $j + 1$.

For a given iteration $j$, the quantities are written as

$$\underline{\lambda}_j = \underline{\lambda}_0 + \sum_{j'=1}^{j} \Delta\underline{\lambda}_{j'}, \underline{u}_j = \underline{u}_0 + \sum_{j'=1}^{j} \Delta\underline{u}_{j'}. \quad (56)$$

Denoting by $\omega$ and $\zeta$ the total numbers of modes generated till this iteration, then eq. (56) can be rewritten as

$$\underline{\lambda}_j = \underline{\lambda}_0 + \sum_{i=1}^{\omega} \bar{\underline{\Lambda}}^i \alpha_{\lambda}^i, \underline{u}_j = \underline{u}_0 + \sum_{j=1}^{\zeta} \bar{\underline{U}}^j \alpha_{u}^j. \quad (57)$$

Defining the cumulative difference operator as $\bar{\Delta}\underline{\lambda} = \sum_{j'=1}^{j} \Delta\underline{\lambda}_{j'}$, eqs. (56) and (57) indicate that

$$\bar{\Delta}\underline{\lambda} = \sum_{i=1}^{\omega} \bar{\underline{\Lambda}}^i \alpha_{\lambda}^i, \bar{\Delta}\underline{u} = \sum_{j=1}^{\zeta} \bar{\underline{U}}^j \alpha_{u}^j. \quad (58)$$

The corresponding strains are written as

$$\underline{\varepsilon}_{\underline{\lambda}j} = \underline{\varepsilon}_{\underline{\lambda}0} + \bar{\Delta}\underline{\varepsilon}_{\underline{\lambda}}, \underline{\varepsilon}_{\underline{u}j} = \underline{\varepsilon}_{\underline{u}0} + \bar{\Delta}\underline{\varepsilon}_{\underline{u}}, \quad (59)$$

with

$$\bar{\Delta}\underline{\varepsilon}_{\underline{\lambda}} = \underline{\varepsilon}(\bar{\Delta}\underline{\lambda}), \bar{\Delta}\underline{\varepsilon}_{\underline{u}} = \underline{\varepsilon}(\bar{\Delta}\underline{u}). \quad (60)$$

The strain corresponding to $\underline{v}$ is obtained as

$$\underline{\varepsilon}_{\underline{v}j} = \underline{\varepsilon}_{\underline{v}0} + \bar{\Delta}\underline{\varepsilon}_{\underline{v}}, \quad (61)$$

with

$$\bar{\Delta}\underline{\varepsilon}_{\underline{v}} = \bar{\Delta}\underline{\varepsilon}_{\underline{v}}^p - \bar{\Delta}\underline{\varepsilon}_{\underline{u}}^p + \bar{\Delta}\underline{\varepsilon}_{\underline{u}} + \bar{\Delta}\underline{\varepsilon}_{\underline{\lambda}}. \quad (62)$$

The corresponding displacement $\underline{v}$ is calculated by solving the discretised equation

$$[K_{aa}] \{\bar{\Delta}\underline{v}_a\} = \{K_{aa}\bar{\Delta}\underline{u}_a + K_{aa}\bar{\Delta}\underline{\lambda}_a + L_{ag}(\bar{\Delta}\underline{\varepsilon}_v^p - \bar{\Delta}\underline{\varepsilon}_u^p)\}, \quad (63)$$

and thereafter using

$$\underline{v}_j = \underline{v}_0 + \bar{\Delta}\underline{v}. \quad (64)$$

For the local stage, no ROM is used and the process is as discussed in section 3.1.2, so is the stopping criterion of the LATIN iteration. It has also to be mentioned that for the initialisation part of the LATIN iteration, no PGD is used, and the methodology is as described in section 3.1. The PGD procedure has been described in algorithm 2.

<caption>**Algorithm 2** The basic problem procedure using PGD</caption>

**Elastic initialisation:**
- Intialise $j = 0$ and $\ell$
- Solve for $u_0, u_0, \lambda_0$
- Compute $\varepsilon_{u0}, \dot{\varepsilon}_{v0}, \varepsilon_{\lambda0}, \sigma_{u0}, \sigma_{v0}$
- Set $\varepsilon_{u0}^p, \varepsilon_{v0}^p, X_{u0}, X_{v0}, Z_{u0}, Z_{v0} = 0$
- Set $c_r = 0$
- Intialise $\omega, \zeta = 0$
**while** $\ell \ge \ell_{tol}$ **do**
**Local step:**
- Update $j = j + 1$
- Use $\varepsilon_{u j-1}, \varepsilon_{v j-1}$
- Solve for $\varepsilon_{u j}^p, \varepsilon_{v j}^p, X_{u j}, X_{v j}, Z_{u j}, Z_{v j}, \sigma_{u j}, \sigma_{v j}$ through eqs. (28) and (29)
**Global step:**
- Represent $\Delta u_j, \Delta v_j, \Delta \lambda_j, \Delta \varepsilon_{u j}, \Delta \varepsilon_{v j}, \Delta \varepsilon_{\lambda j}, \Delta \varepsilon_{u j}^p, \Delta \varepsilon_{v j}^p$ using eq. (17)
**if** $c_r = 0$ **then**
**Expand bases:**
- Update $\omega = \omega + 1, \zeta = \zeta + 1$
- Intialise $\alpha_{\lambda, i'=0}, \alpha_{u, i'=0}$
**for** $i' = 1$ to $5$ **do**
- Solve for $\bar{\Lambda}_{i'}, \bar{U}_{i'}$ using $\alpha_{\lambda, i'-1}, \alpha_{u, i'-1}$ through eq. (42)
- Normalise $\bar{\Lambda}_{i'}, \bar{U}_{i'}$
- Solve for $\alpha_{\lambda, i'}, \alpha_{u, i'}$ using $\bar{\Lambda}_{i'}, \bar{U}_{i'}$ through eq. (44)
**end for**
- Output $\bar{\Lambda}^{\omega}, \bar{U}^{\zeta}, \alpha_{\lambda}^{\omega}, \alpha_{u}^{\zeta}$
**Orthonormalisation:**
- Orthonormalise $\bar{\Lambda}^{\omega}$ with respect to $\{\bar{\Lambda}^i\}_{i=1}^{\omega-1}, \bar{U}^\zeta$ with respect to $\{\bar{U}^j\}_{j=1}^{\zeta-1}$
- Update $\{\alpha_{\lambda}^i\}_{i=1}^{\omega}, \{\alpha_{u}^i\}_{j=1}^{\zeta}$ based on the projections
**if** $\|\alpha_{\lambda}^{\dot{\epsilon}}\|/\|\alpha_{\lambda}^i\| < \ell_{tol}$ **then**
- $e = \{u, \lambda\}, \bar{E} = \{\bar{\Lambda}, \bar{U}\}, \dot{e} = \{\omega, \zeta\}$
- Reject $\alpha_{\lambda}^{\dot{\epsilon}}, \bar{E}^{\dot{\epsilon}}$
- Update $\dot{e} = \dot{e} - 1$
**end if**
**else**
**Update bases:**
- Solve for $\{\Delta\alpha_{\lambda}^i\}_{i=1}^{\omega}, \{\Delta\alpha_{u}^i\}_{j=1}^{\zeta}$ using $\{\bar{\Lambda}^i\}_{i=1}^{\omega}, \{\bar{U}^j\}_{j=1}^{\zeta}$ through eq. (52)
- Update $\{\alpha_{\lambda}^i\}_{i=1}^{\omega}, \{\alpha_{u}^i\}_{j=1}^{\zeta}$ using eq. (48)
**end if**
- Compute $\bar{\Delta}\lambda, \bar{\Delta}u$ through eq. (58) and $\bar{\Delta}\varepsilon_{\lambda}, \bar{\Delta}\varepsilon_{u}$ through eq. (60)
- Compute $\bar{\Delta}\varepsilon_{v}$ using eq. (62) and $\bar{\Delta}v$ using eq. (63)
Compute $u_j, v_j, \lambda_j, \varepsilon_{uj}, \varepsilon_{vj}, \varepsilon_{\lambda j}$ using eqs. (57), (59), (61) and (64)

Calculate $\ell$ from eq. (31) an $\wp$ from eq. (55)

if $\wp > \wp_{tol}$ **then**

Update $c_r = 1$

else

Update $c_r = 0$

end if

end while

## 5 Extension to softening behaviour

The previously described methodology is valid when the Hooke tensor $C$ is constant. In case of damage, the elasticity (Young) modulus $E$ is affected by the damage variable $D$ as $\tilde{E} = E(1-D)$ (where $\tilde{E}$ is the effective modulus), therefore the Hooke tensor is also affected. As the damage variable is a function of space and time, so is the Hooke tensor, therefore an alternative strategy for the problem formulation and resolution is necessary.

For the free energy function, consider the linear elastic case of eq. (16), which in the presence of damage becomes

$$ \psi(\underline{\varepsilon} - \underline{\varepsilon}^p, D, \underline{X}) = \psi_{\varepsilon}(\underline{\varepsilon} - \underline{\varepsilon}^p, D) + \psi_X(\underline{X}), \quad \text{with} \\ \psi_{\varepsilon} = \frac{1}{2}(\underline{\varepsilon} - \underline{\varepsilon}^p) : C(1-D)(\underline{\varepsilon} - \underline{\varepsilon}^p). \qquad (65) $$

The associated or dual variable of $D$ is obtained as

$$ \check{Y} = \frac{\partial \psi_{\varepsilon}}{\partial D} = -\frac{1}{2} (\underline{\varepsilon} - \underline{\varepsilon}^p) : C (\underline{\varepsilon} - \underline{\varepsilon}^p). \qquad (66) $$

The variable $-\check{Y}$ is half of the variation in elastic energy generated by a damage variation at constant stress (see [28] for details), and is referred as strain energy release rate. For the sake convenience, a new variable $Y$ is introduced such that $Y = -\check{Y}$. This variable $Y$, i.e. the strain energy release rate, is considered to be the new dual (or thermodynamic force) of $D$, (see [29, 30] for details). It has also to be mentioned that $Y$ can also be expressed in terms of stress, damage, triaxiality function and triaxility ratio, as given in [30]. The corresponding evolution law is given as

$$ \dot{D} = -\frac{\partial \phi^*}{\partial \check{Y}}, \quad \text{or} \quad \dot{D} = \frac{\partial \phi^*}{\partial Y}. \qquad (67) $$

The variable set $\kappa$ is then expanded to $\kappa = \{\underline{\varepsilon}^e (\underline{\varepsilon} - \underline{\varepsilon}^p), \underline{\varepsilon}^p, \underline{\sigma}, \underline{X}, \underline{Z}, D, Y\}$, and the residual errors defined in eq. (6) becomes

$$ \eta_{\psi} = \psi(\underline{\varepsilon}^e, \underline{X}, D) + \psi^*(\underline{\sigma}, \underline{Z}, Y) \\ - \underline{\sigma} : \underline{\varepsilon}^e - \underline{Z} \cdot \underline{X} - YD, \qquad (68a) $$

$$ \eta_{\phi} = \phi(\underline{\dot{\varepsilon}}^p, -\underline{\dot{X}}, \dot{D}) + \phi^*(\underline{\sigma}, \underline{Z}, Y) \\ - \underline{\sigma} : \underline{\dot{\varepsilon}}^p + \underline{Z} \cdot \dot{X} - Y \dot{D}. \qquad (68b) $$

## 5.1 Classical formulation

To establish the methodology for the resolution of the basic problem, the attempt is first to develop the MCRE formulation without any model reduction. Following the descriptions of eq. section 3.1, the variables associated to $\underline{v}$ are $\kappa_v = \{\underline{\varepsilon}^e(\underline{\varepsilon} - \underline{\varepsilon}^p), \underline{\varepsilon}^p, \underline{\sigma}, \underline{X}, \underline{Z}, D_v, Y_v\}$, and variables associated to $\underline{u}$ are $\kappa_u = \{\underline{\varepsilon}^e(\underline{\varepsilon} - \underline{\varepsilon}^p), \underline{\varepsilon}^p, \underline{\sigma}, \underline{X}, \underline{Z}, D_u, Y_u\}$.

Equation (14) is rewritten as

$$ \mathcal{L} = \int_0^T \int_{\Omega} \left[ \psi(\underline{\varepsilon} - \underline{\varepsilon}^p, D_u, \underline{X}) - \psi(\underline{\varepsilon} - \underline{\varepsilon}^p, D_v, \underline{X}) \right. \\ \left. - \underline{\sigma}_v : \left( (\underline{\varepsilon} - \underline{\varepsilon}^p) - (\underline{\varepsilon}_v - \underline{\varepsilon}^p) \right) \right. \\ \left. - Y_v (D_u - D_v) - \underline{Z}_v \cdot (\underline{X}_u - \underline{X}_v) \right] d\Omega dt \\ + \frac{\gamma}{2} \int_0^T \int_{\Omega_m} \|\underline{u} - \tilde{\underline{u}}\|^2 d\Omega_m dt \\ - \int_0^T \left[ \int_{\Omega} \underline{\sigma}_v : \underline{\varepsilon} (\underline{\lambda}) d\Omega - \int_{\Omega} \underline{f}_d^v \cdot \underline{\lambda} d\Omega \right. \\ \left. - \int_{\partial_f \Omega} \underline{f}_d^s \cdot \underline{\lambda} dS \right] dt. \qquad (69) $$

Similar to section 3.1.1, internal variables $\underline{X}_u, \underline{X}_v, \underline{Z}_v, \underline{Z}_u, \underline{\varepsilon}^p, \underline{\varepsilon}^v, D_u, D_v, Y_u, Y_v$ at the global step are frozen to the value obtained at the last local step.

From eq. (69), eq. (15) can be rewritten as

$$ \int_0^T \int_{\Omega} \mathbf{C}_v (\underline{\varepsilon} - \underline{\varepsilon}^p) : \underline{\varepsilon} (\delta\underline{\lambda}) \, d\Omega dt \\ - \int_0^T \int_{\Omega} \underline{f}_d^v \cdot \delta\underline{\lambda} \, d\Omega dt - \int_0^T \int_{\partial_f \Omega} \underline{f}_d^s \cdot \delta\underline{\lambda} \, dS dt = 0, \qquad (70a) $$

$$ \int_0^T \int_{\Omega} \mathbf{C}_v \left[ (\underline{\varepsilon}_v - \underline{\varepsilon}^p) - (\underline{\varepsilon}_u - \underline{\varepsilon}^p) - \underline{\varepsilon}_\lambda \right] : \underline{\varepsilon} (\delta \underline{v}) \, d\Omega dt = 0, \qquad (70b) $$

$$ \int_0^T \int_{\Omega} \left[ \mathbf{C}_u (\underline{\varepsilon}_u - \underline{\varepsilon}^p) - \mathbf{C}_v (\underline{\varepsilon}_v - \underline{\varepsilon}^p) \right] : \underline{\varepsilon} (\delta \underline{u}) \, d\Omega dt \\ + \gamma \int_0^T \int_{\Omega_m} (\underline{u} - \tilde{\underline{u}}) \cdot \delta \underline{u} \, d\Omega_m dt = 0, \qquad (70c) $$
$$
\begin{align}
\varepsilon_{wl,j}^p &= \varepsilon_{wl-1,j}^p + \Delta t_l \frac{\partial \phi_p^*}{\partial \underline{\sigma}} |_{\sigma_{wl,j}, Z_{wl,j}, D_{wl,j}}, \nonumber \\
X_{wl,j} &= X_{wl-1,j} - \Delta t_l \frac{\partial \phi_p^*}{\partial \underline{Z}} |_{\sigma_{wl,j}, Z_{wl,j}, D_{wl,j}}, \nonumber \\
D_{wl,j} &= D_{wl-1,j} + \Delta t_l \frac{\partial \phi_d^*}{\partial Y} |_{Y_{wl,j}}, \tag{80b}
\end{align}
$$

where it is assumed that $\phi^*(\underline{\sigma}, \underline{Z}, D, Y) = \phi_p^*(\underline{\sigma}, \underline{Z}, D) + \phi_d^*(Y)$, according to [31]; however other variations are also possible. Here $Y$ is the discretised form of $Y$.

For the initialisation step ($j = 0$), the internal variables $Y_{u0}, Y_{v0}$ are computed using eq. (66), while other internal variables $\varepsilon_{u0}^p, \varepsilon_{v0}^p, X_{u0}, X_{v0}, Z_{u0}, Z_{v0}, D_{u0}, D_{v0}$ are considered to be 0, and the stress variables are obtained from eq. (30).

The stopping criterion for the LATIN iterations is similar to section 3.1.2.

## 5.2 ROM-based formulation

Similar to the formulations introduced in section 4, the intention herein is to propose a PGD based formulation with separated space and time problems for the global step. Consider the space and time problems given in section 4.1. Equation (38), can be rewritten for damageable behaviour as

$$
\begin{align}
& \int_{\Omega} \left[ \left( \int_{0}^{T} \mathbf{C}_v \, \mathbf{u}_\alpha \, \mathbf{\alpha}_\alpha \, dt \right) \underline{\underline{\varepsilon}}_{\underline{\underline{\alpha}}_u} + \left( \int_{0}^{T} \mathbf{C}_v \, \mathbf{u}_\alpha \, \mathbf{\alpha}_\alpha \, dt \right) \underline{\underline{\varepsilon}}_{\underline{\underline{\alpha}}_\lambda} \right] \nonumber \\
& \qquad : \underline{\underline{\varepsilon}}(\delta \underline{\underline{\Lambda}}) \, d\Omega \nonumber \\
& \qquad + \int_{\Omega} \left( \int_{0}^{T} \mathbf{C}_v \, \left( \underline{\underline{\varepsilon}}_{v \, j-1} - \underline{\underline{\varepsilon}}_{v \, j-1}^p - \Delta \underline{\underline{\varepsilon}}_{v \, j}^p \right) \mathbf{\alpha}_\lambda \, dt \right) : \underline{\underline{\varepsilon}}(\delta \underline{\underline{\Lambda}}) \, d\Omega \nonumber \\
& \qquad - \int_{\Omega} \left( \int_{0}^{T} \mathbf{f}_d^v \, \mathbf{\alpha}_\lambda \, dt \right) \cdot \delta \underline{\underline{\Lambda}} \, d\Omega \nonumber \\
& \qquad - \int_{\partial_f \Omega} \left( \int_{0}^{T} \mathbf{f}_d^s \, \mathbf{\alpha}_\lambda \, dt \right) \cdot \delta \underline{\underline{\Lambda}} \, dS = 0, \tag{81a} \\
& \int_{\Omega} \left( \int_{0}^{T} (\mathbf{C}_u \, \mathbf{u}_\alpha - \mathbf{C}_v \, \mathbf{u}_\alpha) \, \mathbf{\alpha}_u \, \mathbf{\alpha}_u \, dt \right) \underline{\underline{\varepsilon}}_{\underline{\underline{\alpha}}_u} : \underline{\underline{\varepsilon}}(\delta \underline{\underline{\Lambda}}) \, d\Omega \nonumber \\
& \qquad + \gamma \int_{\Omega_m} \left( \int_{0}^{T} \mathbf{\alpha}_u \, \mathbf{\alpha}_u \, dt \right) \underline{\underline{\overline{U}}} \cdot \delta \underline{\underline{\overline{U}}} \, d\Omega_m \nonumber \\
& \qquad - \int_{\Omega} \left( \int_{0}^{T} \mathbf{C}_v \, \mathbf{\alpha}_u \, \mathbf{\alpha}_\lambda \, dt \right) \underline{\underline{\varepsilon}}_{\underline{\underline{\alpha}}_\lambda} : \underline{\underline{\varepsilon}}(\delta \underline{\underline{\Lambda}}) \, d\Omega \nonumber \\
& \qquad + \int_{\Omega} \left( \int_{0}^{T} \mathbf{C}_u \, \left( \underline{\underline{\varepsilon}}_{u \, j-1} - \underline{\underline{\varepsilon}}_{u \, j-1}^p - \Delta \underline{\underline{\varepsilon}}_{u \, j}^p \right) \mathbf{\alpha}_u \, dt \right) : \underline{\underline{\varepsilon}}(\delta \underline{\underline{\Lambda}}) \, d\Omega \nonumber
\end{align}
$$

$$
\begin{align}
& - \int_{\Omega} \left( \int_{0}^{T} \mathbf{C}_v \, \left( \underline{\underline{\varepsilon}}_{v \, j-1} - \underline{\underline{\varepsilon}}_{v \, j-1}^p - \Delta \underline{\underline{\varepsilon}}_{v \, j}^p \right) \mathbf{\alpha}_u \, dt \right) \nonumber \\
& \qquad : \underline{\underline{\varepsilon}}(\delta \underline{\underline{\Lambda}}) \, d\Omega \nonumber \\
& + \gamma \int_{\Omega_m} \left( \int_{0}^{T} \left( \underline{\underline{u}}_{j-1} - \underline{\underline{\tilde{u}}} \right) \mathbf{\alpha}_u \, dt \right) \cdot \delta \underline{\underline{\Lambda}} \, d\Omega_m = 0. \tag{81b}
\end{align}
$$

The corresponding discretised system becomes

$$
\begin{align}
& \left[ \begin{matrix} \mathbb{J}_{vaa,j} & \mathbb{O}_{vaa,j} \\ (\gamma \Pi^T \Pi)_{aa} \varrho + \mathbb{W}_{uaa,j} - \mathbb{W}_{vaa,j} - \mathbb{J}_{vaa,j} \end{matrix} \right] \left[ \begin{matrix} \underline{\underline{\overline{U}}}_a \\ \underline{\underline{\Lambda}}_a \end{matrix} \right] \nonumber \\
& = \left[ \begin{matrix} \mathbb{M}_a \\ \mathbb{N}_a \end{matrix} \right], \tag{82}
\end{align}
$$

with

$$
\begin{align}
\mathbb{J}_{v,j} &= \llbracket \mathbb{K}_{v,l,j} \alpha_{u,l} \alpha_{\lambda,l} \rrbracket, \quad \mathbb{O}_{v,j} = \llbracket \mathbb{K}_{v,l,j} \alpha_{\lambda,l} \alpha_{\lambda,l} \rrbracket, \tag{83a} \\
\mathbb{W}_{u,j} &= \llbracket \mathbb{K}_{u,l,j} \alpha_{u,l} \alpha_{\lambda,l} \rrbracket, \tag{83b}
\end{align}
$$

$$
\mathbb{W}_{v,j} = \llbracket \mathbb{K}_{v,l,j} \alpha_{u,l} \alpha_{u,l} \rrbracket, \quad \varrho = \llbracket \alpha_{u,l} \alpha_{u,l} \rrbracket, \tag{83c}
$$

$$
\begin{align}
\mathbf{M} &= \llbracket \mathbb{L}_{v,l,j} \Delta \varepsilon_{ul,j}^p + \mathbf{F}_l \nonumber \\
&\quad - \mathbb{L}_{v,l,j} (\varepsilon_{vl,j-1} - \varepsilon_{vl,j-1}^p) \alpha_{\lambda,l} \rrbracket, \nonumber \\
\mathbf{N} &= \llbracket \gamma (\Pi^T \tilde{u}_l - \Pi^T \Pi \mathbf{u}_{l,j-1}) \nonumber \\
&\quad - \mathbb{L}_{ul,j} (\varepsilon_{ul,j-1} - \varepsilon_{ul,j-1}^p - \Delta \varepsilon_{ul,j}^p) \nonumber \\
&\quad + \mathbb{L}_{vl,j} (\varepsilon_{vl,j-1} - \varepsilon_{vl,j-1}^p - \Delta \varepsilon_{ul,j}^p) \alpha_{u,l} \rrbracket. \tag{83d}
\end{align}
$$

The time problem (eq. (40)) can also be rewritten as

$$
\begin{align}
& \int_0^T (\alpha_u \Upsilon + \alpha_\lambda \Gamma) \delta \alpha_\lambda dt \nonumber \\
& + \int_0^T \left( \int_\Omega \mathbf{C}_v \, \left( \underline{\underline{\varepsilon}}_{v \, j-1} - \underline{\underline{\varepsilon}}_{v \, j-1}^p - \Delta \underline{\underline{\varepsilon}}_{v \, j}^p \right) : \underline{\underline{\varepsilon}}_u \, d\Omega \right) \delta \alpha_\lambda dt \nonumber \\
& - \int_0^T \left( \int_\Omega \mathbb{f}_d^v \cdot \underline{\underline{\Lambda}} \, d\Omega \right) \delta \alpha_\lambda dt - \int_0^T \left( \int_{\partial_f \Omega} \mathbb{f}_d^s \cdot \underline{\underline{\Lambda}} \, dS \right) \delta \alpha_\lambda dt = 0, \tag{84a}
\end{align}
$$

$$
\begin{align}
& \int_0^T \alpha_u \Psi \delta \alpha_u dt - \int_0^T \alpha_\lambda \Phi \delta \alpha_u dt \nonumber \\
& + \int_0^T \left( \int_\Omega \mathbf{C}_{u,j} \left( \underline{\underline{\varepsilon}}_{u,j-1} - \underline{\underline{\varepsilon}}_{u,j-1}^p - \Delta \underline{\underline{\varepsilon}}_{u,j}^p \right) : \underline{\underline{\varepsilon}}_u \, d\Omega \right) \delta \alpha_u dt \nonumber \\
& - \int_0^T \left( \int_\Omega \mathbf{C}_{v,j} \left( \underline{\underline{\varepsilon}}_{v,j-1} - \underline{\underline{\varepsilon}}_{v,j-1}^p - \Delta \underline{\underline{\varepsilon}}_{v,j}^p \right) : \underline{\underline{\varepsilon}}_u \, d\Omega \right) \delta \alpha_u dt \nonumber \\
& + \gamma \int_0^T \left( \int_{\Omega_m} \left( \underline{\underline{u}}_{j-1} - \underline{\underline{\tilde{u}}} \right) \cdot \underline{\underline{\overline{U}}} \, d\Omega_m \right) \delta \alpha_u dt = 0, \tag{84b}
\end{align}
$$

with

$$
\begin{align*}
\Upsilon (t) &= \int_{\Omega} \mathbf{C}_v \, \mathbf{u}_{\underline{\underline{\alpha}}_u} : \underline{\underline{\varepsilon}}_{\underline{\underline{\alpha}}_\lambda} \, d\Omega, \quad \Gamma (t) = \int_{\Omega} \mathbf{C}_v \, \mathbf{u}_{\underline{\underline{\alpha}}_\lambda} : \underline{\underline{\varepsilon}}_{\underline{\underline{\alpha}}_\lambda} \, d\Omega, \\
\Psi (t) &= \int_{\Omega} (\mathbf{C}_u \, \mathbf{u}_\alpha - \mathbf{C}_v \, \mathbf{u}_\alpha) : \underline{\underline{\varepsilon}}_{\underline{\underline{\alpha}}_u} : \underline{\underline{\varepsilon}}(\delta \underline{\underline{\Lambda}}) \, d\Omega + \gamma \int_{\Omega_m} \underline{\underline{\overline{U}}} \cdot \underline{\underline{\overline{U}}} \, d\Omega_m,
\end{align*}
$$
$$ \Phi(t) = \int_{\Omega} \mathbf{C}_v \underline{\underline{\varepsilon}}_{\lambda} : \underline{\underline{\varepsilon}}_{\mu} d\Omega. \quad (85) $$

The discretised version (eq. (44)) then becomes

$$ \begin{bmatrix} \gamma_l & \Gamma_l \\ \Psi_l - \Phi_l \end{bmatrix} \begin{bmatrix} \alpha_{ul} \\ \alpha_{\lambda l} \end{bmatrix} = \begin{bmatrix} G_l \\ H_l \end{bmatrix}, \quad (86) $$

with

$$ \begin{aligned} \gamma_l &= \bar{\Lambda}^T \mathbb{K}_{vl,j} \bar{U}, \quad \Gamma_l = \bar{\Lambda}^T \mathbb{K}_{vl,j} \bar{\Lambda}, \\ \Psi_l &= \bar{U}^T (\mathbb{K}_{ul,j} - \mathbb{K}_{vl,j} + \gamma \Pi^T \Pi) \bar{U}, \quad \Phi_l = \bar{U}^T \mathbb{K}_{vl,j} \bar{\Lambda}. \end{aligned} \quad (87) $$

and

$$ G_l = \bar{\Lambda}^T (\mathbb{L}_{vl,j} \Delta \varepsilon_{ul,j}^p + F_l - \mathbb{L}_{vl,j} (\varepsilon_{vl,j-1} - \varepsilon_{vl,j-1}^p)), \quad (88a) $$

$$ \begin{aligned} H_l = \bar{U}^T & \left( \gamma \left( \Pi^T \tilde{u}_l - \Pi^T \Pi u_{l,j-1} \right) \right. \\ & \left. - \mathbb{L}_{ul,j} \left( \varepsilon_{ul,j-1} - \varepsilon_{ul,j}^p - \Delta \varepsilon_{ul,j}^p \right) \right. \\ & \left. + \mathbb{L}_{vl,j} \left( \varepsilon_{vl,j-1} - \varepsilon_{vl,j}^p - \Delta \varepsilon_{ul,j}^p \right) \right). \end{aligned} \quad (88b) $$

The time problem for updating the reduced order bases (eq. (52)) becomes

$$ \begin{bmatrix} [\Upsilon]_l & [\Gamma]_l \\ [\Psi]_l - [\Phi]_l \end{bmatrix} \begin{bmatrix} \{\Delta\alpha_u\}_l \\ \{\Delta\alpha_\lambda\}_l \end{bmatrix} = \begin{bmatrix} \{G\}_l \\ \{H\}_l \end{bmatrix}, \quad (89) $$

with

$$ \begin{aligned} \gamma_{l,j,l} &= \bar{\Lambda}^l T \mathbb{K}_{vl,j} \bar{U}^j, \quad \Gamma_{l,i',l} = \bar{\Lambda}^l T \mathbb{K}_{vl,j} \bar{\Lambda}^{i'}, \\ \Psi_{j,j',l} &= \bar{U}^j T (\mathbb{K}_{ul,j} - \mathbb{K}_{vl,j} + \gamma \Pi^T \Pi) \bar{U}^j, \\ \Phi_{j,l} &= \bar{U}^j T \mathbb{K}_{vl,j} \bar{\Lambda}^l, \end{aligned} \quad (90) $$

and

$$ G_{l,l} = \bar{\Lambda}^l T (\mathbb{L}_{vl,j} \Delta \varepsilon_{ul,j}^p + F_l - \mathbb{L}_{vl,j} (\varepsilon_{vl,j-1} - \varepsilon_{vl,j-1}^p)), \quad (91a) $$

$$ \begin{aligned} H_{j,l} = \bar{U}^j T & \left( \gamma \left( \Pi^T \tilde{u}_l - \Pi^T \Pi u_{l,j-1} \right) \right. \\ & \left. - \mathbb{L}_{ul,j} \left( \varepsilon_{ul,j-1} - \varepsilon_{ul,j}^p - \Delta \varepsilon_{ul,j}^p \right) \right. \\ & \left. + \mathbb{L}_{vl,j} \left( \varepsilon_{vl,j-1} - \varepsilon_{vl,j}^p - \Delta \varepsilon_{ul,j}^p \right) \right). \end{aligned} \quad (91b) $$

The expansion of the reduced order bases, i.e. eqs. (82) and (86), can be resolved as per section 4.1. Similarly, eq. (89) can be solved as per section 4.2 to update the reduced order bases.

However, the calculation of $J_{v,j}, O_{v,j}, W_{u,j}, W_{v,j}$ (eq. (83)) becomes expensive not only because of the integration, but also because $\mathbb{K}_{u,l,j}$ and $\mathbb{K}_{v,l,j}$ have to be computed at each $t_l$, where $l = \{0, \dots, f\}$. Also, for the time problem (eq. (86)), the matrix elements $\gamma_l, \Gamma_l, \Psi_l, \Phi_l$ are to be calculated at each time step $t_l$, which increases the computational cost further.

Consequently, it is necessary to reduce the fidelity of these matrices. To do so, after the local step as given in eqs. (79) and (80), singular value decompositions are performed on the snapshots $l-D_u j$ and $l-D_v j$, from which $t$ and $\mu$ modes are selected such that

$$ I - D_{u,j} = \sum_{i'=1}^{t} \bar{D}_{u}^{i'} \alpha_{d_{u}}^{i'}, \quad (92a) $$

$$ I - D_v j = \sum_{j'=1}^{\mu} \bar{D}_v^{j'} \alpha_{d_v}^{j'}, \quad (92b) $$

where $\bar{D}_u$ and $\bar{D}_v$ are the space functions, and $\alpha_{d_u}, \alpha_{d_v}$ are the time functions. It has to be mentioned that the number of modes to be selected is based on the magnitude of the singular values. The operators defined in eq. (77) can be rewritten in terms of spatial integrals for each mode as

$$ \begin{aligned} \bar{K}_{v, j', j} &= \left[ \int_{\Omega} \mathbb{B}^T \mathbf{C} (1 - \bar{D}_{v}^{j'}) \mathbb{B} d\Omega \right], & \text{assembled} \\ \bar{K}_{u, i', j} &= \left[ \int_{\Omega} \mathbb{B}^T \mathbf{C} (1 - \bar{D}_{u}^{i'}) \mathbb{B} d\Omega \right], & \text{assembled} \\ \bar{L}_{v, j', j} &= \left[ \int_{\Omega} \mathbb{B}^T \mathbf{C} (1 - \bar{D}_{v}^{j'}) d\Omega \right], & \text{assembled} \\ \bar{L}_{u, i', j} &= \left[ \int_{\Omega} \mathbb{B}^T \mathbf{C} (1 - \bar{D}_{u}^{i'}) d\Omega \right], & \text{assembled} \end{aligned} \quad (93) $$

which define eq. (83) as

$$ \begin{aligned} J_{v,j} &= \sum_{j'=1}^{\mu} \bar{K}_{v,j',j} \lfloor \alpha_{d_v}^{j'} \alpha_{u,l} \alpha_{\lambda,l} \rfloor, \quad O_{v,j} \\ &= \sum_{j'=1}^{\mu} \bar{K}_{v,j',j} \lfloor \alpha_{d_v}^{j'} \alpha_{\lambda,l} \alpha_{\lambda,l} \rfloor, \quad (94a) \end{aligned} $$

$$ \begin{aligned} W_{u,j} &= \sum_{i'=1}^{\ell} \bar{K}_{u,i',j} \lfloor \alpha_{d_v}^{i'} \alpha_{u,l} \alpha_{u,l} \rfloor, \quad W_{v,j} \\ &= \sum_{j'=1}^{\mu} \bar{K}_{v,j',j} \lfloor \alpha_{d_v}^{j'} \alpha_{u,l} \alpha_{u,l} \rfloor, \quad (94b) \end{aligned} $$

$$ M = \sum_{j'=1}^{\mu} \bar{L}_{v,j',j} \lfloor (\Delta \varepsilon_{ul,j}^p - \varepsilon_{vl,j-1} + \varepsilon_{vl,j-1}^p) \alpha_{d_v}^{j'} \rfloor $$
<figure><img src="image_4.png" /><figcaption>Fig. 2 Rectangular plate with circular hole under traction</figcaption></figure>

$$
\begin{align}
& \alpha_{\lambda l} \rfloor + \lfloor \underline{F}_l \alpha_{\lambda l} \rfloor, \tag{94c} \\
& N = \lfloor \gamma (\Pi^T \tilde{u}_l - \Pi^T \Pi u_{l,j-1}) \alpha_{ul} \rfloor \nonumber \\
& \quad - \sum_{i'=1}^{l} \overline{\Pi}_{ul}^{i',j} \lfloor (\varepsilon_{ul,j-1} - \varepsilon_{ul,j-1}^p - \Delta\varepsilon_{ul,j}^p) \alpha_{d_{ul}l}^{i'} \alpha_{ul} \rfloor \nonumber \\
& \quad + \sum_{j'=1}^{\mu} \overline{\Pi}_{vj}^{j',j} \lfloor (\varepsilon_{vj,j-1} - \varepsilon_{vj,j-1}^p - \Delta\varepsilon_{ul,j}^p) \alpha_{d_{vj}l}^{j'} \alpha_{ul} \rfloor. \tag{94d}
\end{align}
$$

As $t$, $\mu \ll f$, the computation cost is reduced drastically.
Now, the quantities defined in eqs. (87) and (88) become

$$
\begin{align}
\gamma_l &= \sum_{j'=1}^{\mu} \left( \bar{\Lambda}^T \bar{\mathbb{K}}_{v,j',j} \bar{U} \right) \alpha_{d_v l}^{j'}, \nonumber \\
\Gamma_l &= \sum_{j'=1}^{\mu} \left( \bar{\Lambda}^T \bar{\mathbb{K}}_{v,j',j} \bar{\Lambda} \right) \alpha_{d_v l}^{j'}, \nonumber \\
\Psi_l &= \sum_{i'=1}^{l} \left( \bar{U}^T \bar{\mathbb{K}}_{ul}^{i',j} \bar{U} \right) \alpha_{d_v l}^{i'} \nonumber \\
&\quad - \sum_{j'=1}^{\mu} \left( \bar{U}^T \bar{\mathbb{K}}_{vj,j'} \bar{U} \right) \alpha_{d_v l}^{j'} + \left( \gamma \bar{U}^T \Pi^T \Pi \bar{U} \right), \nonumber \\
\Phi_l &= \sum_{j'=1}^{\mu} \left( \bar{U}^T \bar{\mathbb{K}}_{vj,j'} \bar{\Lambda} \right) \alpha_{d_v l}^{j'}, \tag{95}
\end{align}
$$

and

$$
G_l = \sum_{j'=1}^{\mu} \bar{\Lambda}^T \bar{\Pi}_{v,j',j} \left( \Delta\epsilon_{ul,j}^p - \epsilon_{vl,j-1} + \epsilon_{vj,j-1}^p \right) \alpha_{d_v l}^{j'}
$$

$$
+ \bar{\Lambda}^T \underline{F}_l \qquad (96a)
$$

$$
\begin{equation}
\begin{split}
H_l ={}& \gamma \bar{U}^T \left( \Pi^T \tilde{u}_l - \Pi^T \Pi u_{l,j-1} \right) \\
& - \sum_{i'=1}^{l} \bar{U}^T \bar{\Pi}_{ul}^{i',j} \left( \epsilon_{ul,j-1} - \epsilon_{ul,j-1}^p - \Delta\epsilon_{ul,j}^p \right) \alpha_{d_{ul}l}^{j'} \\
& + \sum_{j'=1}^{\mu} \bar{U}^T \bar{\Pi}_{vj}^{j',j} \left( \epsilon_{vj,j-1} - \epsilon_{vj,j-1}^p - \Delta\epsilon_{ul,j}^p \right) \alpha_{d_{vj}l}^{j'}.
\end{split}
\tag{96b}
\end{equation}
$$

It has to be mentioned that the solution cost of eqs. (95) and (96) is more than that of eqs. (87) and (88) (i.e. without SVD). However, as the solution of eq. (94) is much less than eq. (83), the overall process is cost effective compared to the resolution without SVD. To finish the development, eqs. (90) and (91) can be reformulated very easily similar to eqs. (95) and (96).

6 Numerical results

For the numerical verification of the methodologies described
before, a classical two-dimensional plate with a hole problem
is considered (see fig. 2). The geometry of the plate is defined
by the width $W = 100$ mm, the length $L = 200$ mm, and the
hole diameter $d = 40$ mm. A uniformly distributed, linearly
increasing traction $f_d = 150$ $(t/T)$ N/mm is applied on top
and bottom surfaces, where $t \in [0, T]$ with $T = 200$ s.

Due to symmetry, one-fourth of the plate is consid-
ered for the study with symmetric boundaries as shown
in fig. 2. Now, two types of material behaviour (MB) are
considered for the following studies; one is the classical
elasto-viscoplasticity with isotropic hardening (MB₁), and
the other is the elasto-viscoplasticity coupled with isotropic
hardening and isotropic damage (MB₂).

• MB₁- For elasto-viscoplasticity coupled with isotropic hardening, the free energy function is given as

$$
\psi(\underline{\epsilon} - \underline{\epsilon}^p, r) = \psi_{\epsilon}(\underline{\epsilon} - \underline{\epsilon}^p) + \psi_{r}(r), \quad (97)
$$

with

$$
\psi_{\epsilon} = \frac{1}{2} (\underline{\epsilon} - \underline{\epsilon}^p) : \mathbf{C} (\underline{\epsilon} - \underline{\epsilon}^p), \psi_r = \frac{1}{2} r R_{\infty} r. \quad (98)
$$

These provide the state equations

$$
\underline{\sigma} = \mathbf{C} (\underline{\epsilon} - \underline{\epsilon}^p), \quad R = R_{\infty} r. \qquad (99)
$$

Here $R_\infty$ is the hardening parameter that relates the drag
strain (primal variable) $r$ to the drag stress (dual variable)
<figure><img src="image_5.png" /><figcaption>**Fig. 3** Error curves for the MB₁ solution</figcaption></figure>

R. The Hooke tensor C(E, v) is described by the Young modulus E and Poisson's ratio v. The pseudo dissipation potential in this case is given as per Norton's law as

$$ \phi^* = \frac{1}{n+1} \left\langle \frac{f_y}{K} \right\rangle^{n+1}, \qquad (100) $$

where the yield function given by von Mises criterion can be written as

$$ f_y = |\underline{\sigma}^D|_{eq} - R - \sigma_y, \qquad (101) $$

with $|\circ|_{eq} = \sqrt{3/2\circ:\circ}$, $\sigma_y$ being the yield stress defining the elastic limit, and $\underline{\sigma}^D$ being the deviatoric part of the stress tensor. K and n are the viscous coefficient and exponent, respectively. The evolution laws then become

$$ \underline{\dot{\epsilon}}^p = K^{-n} \langle f_y \rangle^n \left( \frac{3}{2} \frac{\underline{\sigma}^D}{|\underline{\sigma}^D|_{eq}} \right), \dot{r} = K^{-n} \langle f_y \rangle^n. \qquad (102) $$

It has to mentioned here that the cumulative plastic strain p, defined as

$$ \dot{p} = \sqrt{\frac{2}{3} \underline{\dot{\epsilon}}^p : \underline{\dot{\epsilon}}^p}, \qquad (103) $$

is basically the same as the drag strain, i.e

$$ \dot{p} = \dot{r}. \qquad (104) $$

* MB₂- For elasto-viscoplasticity coupled with isotropic hardening and isotropic damage, the free energy function

is given as

$$ \psi(\underline{\epsilon} - \underline{\epsilon}^p, D, r) = \psi_{\epsilon}(\underline{\epsilon} - \underline{\epsilon}^p, D) + \psi_{r}(r), \qquad (105) $$

with

$$ \psi_{\epsilon} = \frac{1}{2} (\underline{\epsilon} - \underline{\epsilon}^p) : C(1-D) (\underline{\epsilon} - \underline{\epsilon}^p), \quad \psi_{r} = \frac{1}{2} r R_{\infty} r. \qquad (106) $$

These provide the state equations

$$ \begin{aligned} \underline{\sigma} &= C(1-D)(\underline{\epsilon} - \underline{\epsilon}^p), \\ Y &= \frac{1}{2}(\underline{\epsilon} - \underline{\epsilon}^p) : C(\underline{\epsilon} - \underline{\epsilon}^p), \\ R &= R_{\infty}r. \end{aligned} \qquad (107) $$

The dual potential is separated as (according to [32])

$$ \psi^*(\underline{\sigma}, R, D, Y) = \psi_p^*(\underline{\sigma}, R, D) + \phi_D^*(Y), \qquad (108) $$

with

$$ \begin{aligned} \psi_p^*(\underline{\sigma}, R, D) &= \frac{1}{n+1} \left\langle \frac{f_y}{K} \right\rangle^{n+1}, \\ \phi_D^*(Y) &= \frac{1}{s+1} \left\langle \frac{Y - Y_{th}}{S} \right\rangle^{s+1}, \end{aligned} \qquad (109) $$

where $Y_{th}$ is the threshold value beyond which damage evolution occurs, and s, S are the material parameters dictating the damage law. The yield function for this case is given by
<figure><img src="image_6.png" /><figcaption>Fig. 4 Initial space-time PGD modes describing the Lagrange multiplier $\lambda$ for $MB_1$</figcaption></figure>

The evolution laws are given by

$$f_y = \left| \frac{\sigma^D}{1-D} \right|_{eq} - R - \sigma_y. \qquad (110)$$
$$\dot{\underline{\sigma}}^p = K^{-n} \langle f_y \rangle^n \left( \frac{3}{2} \frac{\underline{\sigma}^D}{\left| \frac{\underline{\sigma}^D}{1-D} \right|_{eq}} \right) \frac{1}{1-D},$$
<figure><img src="image_7.png" /><figcaption>Fig. 5 Initial space-time PGD modes describing the displacement $\underline{u}$ for MB$_{1}$</figcaption></figure>

$$ \dot{r} = K^{-n} \langle f_y \rangle^n, \dot{D} = S^{-s} \langle Y - Y_{th} \rangle^s, \quad (111) $$

and in the presence of damage

$$ \dot{p} = \frac{\dot{r}}{1 - D}. \quad (112) $$

252 four-node quadrilateral elements, with four Gauss points per element, are used to discretise the structure to be studied. This generates a total of 288 nodes, with two degrees of freedom per node, and 1008 Gauss points. For the time discretisation, 101 time steps are used to discretise the time frame [0, 200 s].
<figure><img src="image_8.png" /><figcaption>Fig. 6 Displacement fields ($\mathbf{u}$, $\mathbf{v}$) at the final time step for MB₁</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 7 Fields of drag strains ($r_u$, $r_v$) at the final time step, and maximum drag strain evolutions for MB₁</figcaption></figure>

At first, direct problems are solved for both material behaviours in order to generate the synthetic data $\tilde{u}$. The reference material parameters used for the direct problem solutions are $E^{ref} = 200$ GPa, $\nu^{ref} = 0.3$, $\sigma_y^{ref} = 400$ MPa, $R_\infty^{ref} = 10$ GPa, $K^{ref} = 400$ MPas¹/ⁿ, $n^{ref} = 7$, $s^{ref} = 2.5$, $S^{ref} = 1.5$ MPas¹/ₛ, $Y_{th}^{ref} = 0.45$ MPa.

<figure><img src="image_10.png" /><figcaption>Fig. 8 Fields of von Misses stresses ($\sigma_u^{vm}$, $\sigma_v^{vm}$) at the final time step, and maximum von Misses stress evolutions for MB₁</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 9 Fields of residual errors ($\eta_\psi$, $\eta_\phi$) at the final time step for MB₁</figcaption></figure>

Thereafter, for the inverse analyses, the intend is to perform two separate studies. The first one is to solve the basic problem using eq. (4) for MB₁ and using eq. (5) for MB₂. The second study is to perform complete inverse simulation using the total MCRE process including the ROMs.

## 6.1 Solving the basic problem

The discretisation for this part of the study is the same as that of the direct problem. The scaling factor $\gamma$ is considered to be 10⁷, based on Morozov's principle (see [17] for details). Reference basic problem solutions (without any reduced order modelling) are performed using the reference material parameters, first for MB₁ through the process described in eq. section 3.1, and then for MB₂ through the procedure described in section 5.1. Thereafter ROM based formulations
<figure><img src="image_12.png" /><figcaption>**Fig. 10** Error curves for the MB₂ solution</figcaption></figure>

are used to solve the basic problem, according to section 4 for MB₁, and according to section 5.2 for MB₂.

For MB₁ ROM-based solution, the error curves are shown in fig. 3. They basically represent the LATIN error indicator (defined in eq. (31)) with respect to the LATIN iterations. Section 6.1 also depicts the evolution of the error indicator with respect to the number of PGD modes generated. A total of 9 PGD modes were generated to approximate both the displacement field and the Lagrange multiplier. As per the CPU time is concerned, the ROM-based solution time was observed to be 60% of the solution time of the reference problem.

The initial four PGD modes used to approximate $\lambda$ are shown in fig. 4, and the same for $\underline{u}$ is depicted in fig. 5.

The final displacement field for MB₁ is depicted in fig. 6, and the displacement fields $\underline{u}$ and $\underline{v}$ have errors of $3.5 \times 10^{-15}\%$ and $2.9 \times 10^{-13}\%$ with respect to the reference solutions. The error $er_\circ$ of any quantity $\circ$ is given as

$$er_\circ = \left[ \frac{\int_0^T \int_\Omega ((\circ - \circ_{ref}) : (\circ - \circ_{ref})) d\Omega dt}{\int_0^T \int_\Omega (\circ_{ref} : \circ_{ref}) d\Omega dt} \right]^{1/2} \quad (113)$$

The drag strains $r_u, r_v$ at the final time step, for MB₁, are shown in fig. 7. This figure also shows the evolution of the maximum drag strains with respect to time. The errors for $r_u, r_v$ are obtained to be $3.1 \times 10^{-11}\%$, $5.5 \times 10^{-11}\%$, and their corresponding duals $R_u, R_v$ provide the same error values.

The von Misses stresses $\sigma_u^{vm}, \sigma_v^{vm}$ at the final time step for MB₁, along with the time evolutions of their maximum

values, are depicted in fig. 8. The stresses $\sigma_{\underline{u}}, \sigma_{\underline{v}}$ have errors of $3.2 \times 10^{-13}\%$ and $3.7 \times 10^{-13}\%$.

The errors corresponding to the plastic strains $\varepsilon_{\underline{u}}^p, \varepsilon_{\underline{v}}^p$ for MB₁ are obtained to be $3.1 \times 10^{-11}\%$, $5.3 \times 10^{-11}\%$, and $3.4 \times 10^{-13}\%$, $4.9 \times 10^{-13}\%$ for the total strains $\varepsilon_{\underline{u}}, \varepsilon_{\underline{v}}$.

The residual errors, defined in eq. (6) is plotted at the end of loading in fig. 9, and as expected the magnitudes are quiet low (as the basic problem is solved using reference parameters), and maximum values are concentrated at the plastified zone.

Similarly, the same exercise is conducted for MB₂. The error curves, for both classical solution and ROM-based solution are shown in fig. 10. 15 PGD modes were generated to approximate both $\underline{u}$ and $\underline{\lambda}$. The CPU time for the ROM-based solution was found to be 70% of the classical resolution.

The separated PGD modes (first four) for both $\underline{u}$ and $\underline{\lambda}$ are depicted in figs. 11 and 12.

The displacement fields $\underline{u}, \underline{v}$ for MB₂ have depicted errors of $6.3 \times 10^{-13}\%$ and $1.2 \times 10^{-11}\%$, and the final field is shown in fig. 13.

The von Misses stresses $\sigma_u^{vm}, \sigma_v^{vm}$ at the final time step for MB₂, along with the time evolutions of their maximum values are depicted in fig. 14. The stresses $\sigma_{\underline{u}}, \sigma_{\underline{v}}$ have errors of $1.5 \times 10^{-11}\%$, $5.5 \times 10^{-10}\%$.

The cumulative plastic strains $p_u, p_v$ at the final time step, for MB₂, are shown in fig. 15. This figure also shows the evolution of the maximum drag strains with respect to time. The errors corresponding to the plastic strains $\varepsilon_{\underline{u}}^p, \varepsilon_{\underline{v}}^p$ for MB₂ are obtained to be $2.9 \times 10^{-11}\%$, $4.9 \times 10^{-11}\%$, and $1.6 \times 10^{-11}\%$, $6.9 \times 10^{-10}\%$ for the total strains $\varepsilon_{\underline{u}}, \varepsilon_{\underline{v}}$.

The errors for $r_u, r_v$ are obtained to be $2.8 \times 10^{-11}\%$, $4.8 \times 10^{-11}\%$ for MB₂, and their corresponding duals $R_u, R_v$ provide the same error values. $r_u, r_v$ at the final time step for
<figure><img src="image_13.png" /><figcaption>Fig. 11 Initial space-time PGD modes describing the displacement $\underline{u}$ for MB$_2$</figcaption></figure>

MB$_2$, along with the time evolutions of their maximum values, are depicted in fig. 16.

The damage variables $D_u, D_v$ at the final time step for MB$_2$, along with the time evolutions of their maximum values, are depicted in fig. 17. Errors corresponding to $D_u, D_v$ are $1.3 \times 10^{-11}\%$, $1.7 \times 10^{-11}\%$, and errors for $Y_u, Y_v$ are $2.6 \times 10^{-11}\%$, $8.3 \times 10^{-10}\%$. It has to be mentioned that the

8 SVD modes were generated (for both $D_u, D_v$) from the snapshots of the damage variables (as given in section 5.2).

The residual errors, defined in eq. (68) is plotted at the end of loading in fig. 18, and similar to fig. 9 the magnitudes are quiet low, and maximum values are concentrated at the plastified (and damaged) zone.
<figure><img src="image_14.png" /><figcaption>Fig. 12 Initial space-time PGD modes describing the Lagrange multiplier $\lambda$ for MB$_2$</figcaption></figure>

From these analyses, it is seen that that admissible states provide excellent accuracy for the ROM-based solutions of the basic problem, compared to the classical solutions. Also the ROM-based solutions were found to be more frugal compared to the classical solutions.

## 6.2 Identification of material parameters

In this part, the first exercise is to solve the total inverse problem, with the basic problem being solved through section 4 for MB$_1$, and through section 5 for MB$_2$. The elastic parameters ($E, \nu$) for both cases are considered to be known parameters.
<figure><img src="image_15.png" /><figcaption>Fig. 13 Displacement fields ($u, v$) at the final time step for MB2</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 14 Fields of von Misses stresses ($\sigma_{u}^{vm}, \sigma_{v}^{vm}$) at the final time step, and maximum von Misses stress evolutions for MB2</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 15 Fields of cumulative plastic strains ($p_u, p_v$) at the final time step, and maximum cumulative plastic strain evolutions for MB2</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 16 Fields of drag strains ($r_u, r_v$) at the final time step, and maximum drag strain evolutions for MB2</figcaption></figure>

Figure 19 depicts the evolution of parameters for MB1 at different gradient descent iterations, along with the corresponding MCRE functional values. It is clear from fig. 19 that the parameters converge to their reference values with successive iterations of the gradient descent algorithm, and the MCRE functionals decrease to their minimal values. Each of the material parameters shown in fig. 19 are identified individually, with the rest considered to be known and equal to reference values.

A similar exercise is performed for MB2, as depicted in fig. 20. Similar to the previous case, each of the material parameters shown in fig. 20 are identified individually, with the rest considered to be known and equal to reference values.

Now, for the next set of numerical tests, percentage white noises are added as a spurious input in the measurement data as

$$ \tilde{u} = \tilde{u}_0 + \delta \tilde{u}, \qquad (114) $$
<figure><img src="image_19.png" /><figcaption>Fig. 17 Fields of damage variables ($D_u, D_v$) at the final time step, and maximum damage evolutions for MB$_2$</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 18 Fields of residual errors ($\eta_\psi, \eta_\phi$) at the final time step for MB$_2$</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 19 Convergence curves of material parameters with respect to the MCRE for MB$_1$</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 20 Convergence curves of material parameters with respect to the MCRE for MB$_2$</figcaption></figure>

where $\tilde{u}_0$ is the measurement data from the direct problem and $\delta\tilde{u}$ is a sample of Gaussian white noise. Different levels of Gaussian white noises (0.5% and 1%) are added to the non-spurious synthetic data from the direct problem resolutions. For each noise level, 10 different random samples (i.e. $\delta\tilde{u}$) are used, and then the ROM based MCRE inverse analysis is performed to determine the parameters using the spurious data. The corresponding means and standard deviations at each noise level, of the errors with respect to the reference parameters are calculated (see figs. 21 and 22).

It is clear from fig. 21, that the accuracies and the precisions decrease with increase in noise level, also from fig. 21, it can be seen that $\sigma_y$ is most precisely and accurately identified, and $R_\infty$ is most susceptible to noise.

Similar to MB$_1$, it is clear from fig. 22, that the accuracies and the precisions decrease with increase in noise level. Also
<figure><img src="image_23.png" /><figcaption>**Fig. 21** Accuracy and precision of identified material parameters with respect to the reference parameters for MB₁</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Fig. 22** Accuracy and precision of identified material parameters with respect to the reference parameters for MB₂</figcaption></figure>

from fig. 22, it can be seen that $\sigma_y$ and $Y_{th}$ are the best identified parameters, and $R_\infty$ is the worst identified parameter.

These exercises were conducted considering the complete time domain, to establish the proof of concept of the ROM based basic problems in the overall framework of the inverse analysis.

## 7 Conclusion

This article essentially depicted the accuracy and frugality of PGD reduced order models in the solution of MCRE-based inverse approaches for nonlinear model updating. The methodology has been tested on viscoplastic and coupled damage-viscoplastic material identifications, and the veracity along with numerical expense has been observed to be appreciable. For damageable behaviours, a POD-based model reduction method has also been coupled, to treat the variable Hooke tensor efficiently. The overall attempt to introduce ROM in MCRE has been quite successful. A further work will consist in implementing the proposed strategy in a sequential model updating method, using a modified dual Kalman filter introduced in [21, 33] and which leans on MCRE. In addition, an application of the proposed strategy with real data, as well as an extension to the dynamics/vibration regime, are forecast in the framework of an ongoing project dealing with real-time SHM [1].

**Funding** Open access funding provided by CentraleSupélec. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Grant agreement No. 101002857).

**Data Availability** Data, materials, and code that support the findings of this study are available from the authors upon reasonable request.

## Declarations

**Competing interests** The authors have no competing interests to declare that are relevant to the content of this article.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
## References

1. Chamoin L (2021) Merging advanced sensing techniques and simulation tools for future structural health monitoring technologies. Project Repository J 10:124–127

2. Blasch EP, Darema F, Bernstein D (2022). In: Blasch EP, Darema F, Ravela S, Aved AJ (eds) Introduction to the Dynamic Data Driven Applications Systems (DDDAS) Paradigm. Springer International Publishing, Cham, pp 1–32

3. Chinesta F, Cueto E, Abisset-Chavanne E, Duval JL, Khaldi FE (2020) Virtual, Digital and Hybrid Twins: A New Paradigm in Data-Based Engineering and Engineered Data. Arch Comput Methods Eng 27(1):105–134

4. Hild F, Roux S (2006) Digital Image Correlation: from Displacement Measurement to Identification of Elastic Properties - a Review. Strain 42(2):69–80

5. Chamoin L, Farahbakhsh S, Poncelet M (2022) An educational review on distributed optic fiber sensing based on Rayleigh backscattering for damage tracking and structural health monitoring. Meas Sci Technol 33(12):124008

6. Avril S, Bonnet M, Bretelle AS, Grédiac M, Hild F, Ienny P et al (2008) Overview of Identification Methods of Mechanical Parameters Based on Full-field Measurements. Exp Mech 48(4):381–402

7. Grédiac M, Toussaint E, Pierron F (2002) L'identification des propriétés mécaniques de matériaux avec la méthode des champs virtuels, une alternative au recalage par éléments finis. Comptes Rendus Mécanique 330(2):107–112

8. Bui HD, Constantinescu A, Maigre H (2005). In: Mróz Z, Stavroulakis GE (eds) The Reciprocity Gap Functional for Identifying Defects and Cracks. Springer Vienna, Vienna, pp 17–54

9. Ereiz S, Duvnjak I, Fernando J-AJ (2022) Review of finite element model updating methods for structural applications. Structures 41:684–723

10. Amiot F, Périé JN, Roux S (2013) 12. Equilibrium Gap Method. John Wiley & Sons, Ltd, pp 331–362

11. Ladevèze P, Chamoin L (2016). In: Chamoin L, Díez P (eds) The Constitutive Relation Error Method: A General Verification Tool. Springer International Publishing, Cham, pp 59–94

12. Ben Azzouna M, Feissel P, Villon P (2015) Robust identification of elastic properties using the Modified Constitutive Relation Error. Comput Methods Appl Mech Eng 295:196–218

13. Huang S, Feissel P, Villon P (2016) Modified constitutive relation error: An identification framework dealing with the reliability of information. Comput Methods Appl Mech Eng 311:1–17

14. Benady A, Baranger E, Chamoin L (2024) Unsupervised learning of history-dependent constitutive material laws with thermodynamically-consistent neural networks in the modified Constitutive Relation Error framework. Comput Methods Appl Mech Eng 425:116967

15. Bhattacharyya M, Feissel P (2023) A Bayesian inference approach for parametric identification through optimal control method. Int J Numer Meth Eng 124(14):3145–3165

16. Marchand B, Chamoin L, Rey C (2019) Parameter identification and model updating in the context of nonlinear mechanical behaviors using a unified formulation of the modified Constitutive Relation Error concept. Comput Methods Appl Mech Eng 345:1094–1113

17. Nguyen HN, Chamoin L, Ha MC (2022) mCRE-based parameter identification from full-field measurements: Consistent framework, integrated version, and extension to nonlinear material behaviors. Comput Methods Appl Mech Eng 400:115461

18. Ladevèze P (1999) Nonlinear Computational Structural Mechanics: New Approaches and Non-Incremental Methods of Calculation. Mechanical Engineering Series, Springer, New York

19. Chinesta F, Ammar A, Leygue A, Keunings R (2011) An overview of the proper generalized decomposition with applications in computational rheology. J Non-newton Fluid Mech 166(11), 578–592

20. Ladevèze P, Passieux JC, Néron D (2010) The LATIN multiscale computational method and the Proper Generalized Decomposition. Comput Methods Appl Mech Eng 199(21–22):1287–1296

21. Marchand B, Chamoin L, Rey C (2016) Real-time updating of structural mechanics models using Kalman filtering, modified constitutive relation error, and proper generalized decomposition. Int J Numer Meth Eng 107(9):786–810

22. Chamoin L, Allier PE, Marchand B (2016) Synergies between the constitutive relation error concept and PGD model reduction for simplified V&V procedures. Adva Model Simulation Eng Sci 3(18):1–26

23. Bhattacharyya M, Feissel P (2025) A LATIN-PGD reduced order approximation dedicated to the solution of an optimal control based identification strategy for non-linear constitutive parameters. Int J Solids Struct 309:113189

24. Nouy A (2010) A priori model reduction through Proper Generalized Decomposition for solving time-dependent partial differential equations. Comput Methods Appl Mech Eng 199(23):1603–1626

25. Chatterjee A (2000) An introduction to the proper orthogonal decomposition. Curr Sci 78(7):808–817

26. Pinnau R (2008) In: Schilders WHA, van der Vorst HA, Rommes J, (eds). Model Reduction via Proper Orthogonal Decomposition. Berlin, Heidelberg: Springer, Berlin Heidelberg; pp. 95–109

27. Bhattacharyya M (2018) A model reduction approach in space and time for fatigue damage simulation [Thesis]. Universität Paris Saclay (COmUE)

28. Lemaitre J, Chaboche JL (1990) Mechanics of Solid Materials. Cambridge University Press, Cambridge

29. Lemaitre J (1996) A Course on Damage Mechanics, 2nd edn. Springer-Verlag, Berlin, Heidelberg

30. Lemaitre J, Desmorat R (2005) Engineering Damage Mechanics: Ductile, Creep. Springer-Verlag, Berlin, Heidelberg, Fatigue and Brittle Failures

31. Ladevè P, Moës N, Douchin B (1999) Constitutive relation error estimators for (visco)plastic finite element analysis with softening. Comput Methods Appl Mech Eng 176(1):247–264

32. Bhattacharyya M, Fau A, Nackenhorst U, Néron D, Ladevèze P (2018) A LATIN-based model reduction approach for the simulation of cycling damage. Comput Mech 62(4):725–743

33. Diaz M, Charbonnel PÉ, Chamoin L (2023) A new Kalman filter approach for structural parameter tracking: Application to the monitoring of damaging structures tested on shaking-tables. Mech Syst Signal Process 182:109529

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.