A data-driven reduced-order surrogate model for entire elastoplastic simulations applied to representative volume elements

This contribution discusses surrogate models that emulate the solution field(s) in the entire simulation domain. The surrogate uses the most characteristic modes of the solution field(s), in combination with neural networks to emulate the coefficients of each mode. This type of surrogate is well known to rapidly emulate flow simulations, but rather new for simulations of elastoplastic solids. The surrogate avoids the iterative process of constructing and solving the linearized governing equations of rate-independent elastoplasticity, as necessary for direct numerical simulations or (hyper-)reduced-order-models. Instead, the new plastic variables are computed only once per increment, resulting in substantial time savings. The surrogate uses a recurrent neural network to treat the path dependency of rate-independent elastoplasticity within the neural network itself. Because only a few of these surrogates have been developed for elastoplastic simulations, their potential and limitations are not yet well studied. The aim of this contribution is to shed more light on their numerical capabilities in the context of elastoplasticity. Although more widely applicable, the investigation focuses on a representative volume element, because these surrogates have the ability to both emulate the macroscale stress-deformation relation (which drives the multiscale simulation), as well as to recover all microstructural quantities within each representative volume element.

domain for different scenarios, and (3) the range of variations for which the framework was trained.Predecessors to these studies 24,25 , e.g.Nagel et al. 32 , also did not reveal these points.
On the other hand, the recent study of Im et al. 26 clearly mentions the limitation to small deformations and rotations.Furthermore, plasticity-induced localization only appears at one or a few location(s) that do not change during the course of the simulations.Moreover, the applied boundary conditions seem to be limited to a single (history-dependent) 'variable' (e.g. a horizontal excitation or a vertical force).The numerical limitations of the investigation 26 are furthermore characterized by the fact that only eight modes suffice to obtain highly accurate results.
Aim.The aim of the current contribution is to numerically reveal the capabilities and limitations of RNN-POD frameworks, by considering: 1. the entire, current displacement field as the output of interest, 2. a relatively simple RNN (with a single central layer of gated recurrent units (GRUs), sandwiched between two layers of conventional neurons), 3. large deformations with plasticity-induced localizations at several locations in the domain, as well as large variations of the imposed deformations, and 4. the application to RVE computations with periodic boundary conditions in order to serve multiscale analyses.
However, embedding the RNN-POD framework in a multiscale simulation is not considered yet.Also, the accuracy of the consistent tangent moduli is not quantatively considered yet-only the accuracy of the macroscale stress-deformation response, and the accuracy of the microstructural field of the plastic multiplier.
Outline.The remainder of this paper is organized as follows.The next section concisely discusses how to efficiently train the surrugate for an RVE, by removing the rigid body rotations from the macroscale deformation gradient tensor (which is imposed on the boundaries of the RVE).The direct numerical simulations are concisely discussed in section "Direct numerical simulations".Section "POD-based model-order-reduction" describes a conventional POD-based reduced-order-model (ROM) followed by the description of the neural network architecture in section "Recurrent neural network".Section "Data collection" discusses the generation of training data.The results are divided over two sections: the results of the RNN are presented in section "RNN results" and the mechanical results, which are of course of primary interest, are presented in section "Mechanical results".A conclusion is presented in section "Conclusion".

Removing rigid body rotations from the RVE and its surrogate
Surrogate models based on neural networks require a substantial amount of training data.The rather slow, conventional POD-based ROM (which in turn also needs training data, but substantially less than the RNN-POD surrogate), must generate this large amount of data.Therefore, it is convenient if the amount of training data can be restricted as much as possible.
In the case of RVE computations, which are governed by an applied macroscopic deformation, it is more efficient to train the emulator for 'pure' macroscopic deformations (meaning that rigid body rotations are extracted from the macroscopic deformation) than for the actual macroscopic deformations.Extracting the rigid body rotations from the macroscale deformation reduces the number of inputs of the RNN-POD emulator from four to three in 2D (and from nine to six in 3D).The current section shows that the 'pure' macroscale stressdeformation relation (i.e.without rigid body rotations) can be rapidly transformed into the actual macroscale stress-deformation relation (with rigid body rotations).

Macroscopic stress-deformation relation.
First-order multiscale (FE 2 ) simulations comprise interwoven simulations at two length scales [33][34][35][36] .The macroscale simulation, with its modelling domain at the size of the entire application or product, uses an implicit constitutive relation (as well as its sensitivity), which is obtained from microstructural RVE simulations (with its modelling domain at the size of the microstructure).The macroscale deformation gradient tensor at time t, F M (t) (at a quadrature point of the macroscale discre- tization), is imposed on an RVE.This results in a displacement field within the RVE, in turn yielding a stress field inside the RVE.The macroscale (homogenized) 1 st Piola Kirchoff (PK) stress tensor, P M (F M (t), z µ (t)) (which depends on the macroscale deformation gradient tensor and on a yet unspecified number of plastic variables), can be computed using the field of the microstructural 1 st (PK) stress tensors that occurs inside the RVE, where F µ and − → X µ denote the (microstructural) deformation gradient tensor in the RVE and the reference location vector in the RVE, respectively.V 0 denotes the reference volume of the RVE (i.e.area in 2D).
As mentioned before, the implicit relation P M (F M (t), z M (t)) drives the simulation at the macroscale.The sensitivity of the macroscale stress with respect to the macroscale deformation gradient tensor, i.e. ∂P M ∂F M , is therefore also required at the macroscale.
The issue of interest here is that macroscale deformation gradient tensor F M embeds rigid body rotations and we want to extract them from the deformation (and consequently, also from the macroscale stress).To (1) www.nature.com/scientificreports/this purpose, we multiplicatively decompose the macroscale deformation gradient tensor at time t in a (proper orthogonal) rotation tensor, R M (t) and a (symmetric) 'pure' deformation gradient tensor, U M (t) , according to: where ( • ) denotes an inner product.Consequently, if a pure macroscopic deformation, U M (t) , is imposed to the RVE (instead of F M (t) ), the resulting macroscale stress tensor, PM (U M (t), z M (t)) (which is still obtained according to the volume average in Eq. ( 1)), is related to the macroscale stress tensor that embeds rigid body rotations, i.e.P M (F M (t), z M (t)) , as follows: The sensitivities of the two stress tensors with respect to the macroscale deformation gradient tensor are related as follows: where superscript T denotes the transpose, superscript T 12 denotes a conjugate transpose where the first and second base vectors change sequence, and ( : ) denotes the double inner product.
This indicates that if the RVE is exposed to macroscale deformations without rigid body rotations, the resulting macroscale stress tensor (as well as its sensitivity) can straightforwardly be related to the relevant macroscale stress tensor (as well as its sensitivity) that is affected by rigid body rotations.

Direct numerical simulations
In the current section, we discuss the RVE simulations.The RVE simulations are the direct numerical simulations (DNS) which the POD-based ROM, and particularly the RNN-POD emulator must mimic.
The DNS are 2D RVE computations with stiff, elastic particles embedded in an elastoplastic matrix (see Fig. 1).In relation to the previous section, this entails that our simulations are microstructural simulations (superscript µ in the previous section) and the macroscale (superscript M in the previous section) is completely neglected.Thus, unless specified otherwise, all quantities discussed in the current section are microstructural quantities, but superscript µ is omitted in order to ease the notation.
Constitutive model.We consider plane strain conditions and employ isoparametric bilinear quadrilateral finite elements with four quadrature points.An F-bar method is utilized to alleviate locking due to the incompressibility of the plastic deformation.Within the F-bar framework, the volume change of the deformation at a quadrature point is replaced with the volume change at the center of the element.Thus, the deformation gradient tensor that is considered for the constitutive model in each quadrature point, F , is written as: (2) where F denotes the deformation gradient tensor in the quadrature point, J = det(F) the volume change in the quadrature point, and J c the volume change in the center of the element.The deformation gradient tensor that is used in the actual constitutive model, F , is multiplicatively decom- posed into an elastic (subscript e) and a plastic (subscript p) deformation gradient tensor: The 1 st Piola-Kirchhoff (PK) stress tensor is defined according to the following derivative: with the following strain energy density: where E and ν denote Young's modulus and Poisson's ratio, respectively.Furthermore: I e = tr(F T e • F e ) and J e = det(F e ).
In order to assess if plastic yielding occurs, we first relate the 1 st PK stress tensor to the Mandel stress, M , according to: then we compute the deviatoric Mandel stress according to: (where I denotes the identity tensor) and subsequently, we use the deviatoric Mandel stress in the following von Mises-type yield function: where material parameters M 0 , h and m denote the initial yield stress, the hardening modulus and an exponential hardening parameter, respectively.Furthermore, denotes the plastic multiplier, with an initial value of zero.
The evolution of the plastic multiplier is linked to the evolution of the plastic deformation according to the following associated flow rule: The Karush-Kuhn-Tucker conditions close the constitutive model: The same constitutive description is used for the elastic particles and the elastoplastic matrix.However, in order to guarantee that the particles only deform elastically, their initial yield stress ( M 0 ) is set to a substantially large value that will never be reached during the simulations.Furthermore, Young's modulus of the elastic particles is set to twenty times that of the elastoplastic matrix.The actual parameter values are presented in section 6.

Periodic RVE simulations.
Because we are interested in RVE simulations, the boundary conditions are dictated by the 'pure' macroscale deformation, U M (t) .The fact that the pure macroscale deformation depends on pseudo-time t, not a single U M is provided, but the entire sequence for all time increments.
The current boundary conditions of the RVE simulations are governed by the current, 'pure' deformation gradient tensor, U M .We impose the macroscale pure deformation via periodic boundary conditions.This entails that the displacement of a node on the bottom boundary (node a) is linked to the displacement of a node on the top boundary (node b) according to: where u and − → X denote a nodal displacement and a nodal reference location, respectively.The same constraints are present between nodes on the left boundary and nodes on the right boundary.To enable the application of the displacement constraints between pairs of nodes on opposing boundaries, we employ a periodic mesh (see Fig. 1).
(5) F = J In the RVE simulations, the left bottom node obtains displacement boundary conditions of zero.This entails that the displacements of the other corner nodes are known and can thus also be prescribed using displacement boundary conditions.However, the displacements of the remaining nodes on the boundaries, are not known in advance and hence, the aforementioned constraint must be included as an actual constraint (in fact, as two scalar constraints in a 2D setting as here).We have chosen to enforce the constraints using Lagrange multipliers.
After linearization, the DNS implies that the following system of linear equations must be constructed and solved for each iteration, for each time increment: where column u collects the displacement components of all nodes at an intermediate solution, column z the plastic variables in all quadrature points ( and the components of F p ), column g the Lagrange multipliers, column c the constraints due to the periodic boundary conditions, column f ext the components of the reaction forces, column f int the components of the internal forces and matrix K int the derivatives of the internal force components with respect to the displacement components.du and dg denote the corrections to u and g , which are to be computed.

POD-based model-order-reduction
Now that the DNS are discussed, we continue the current section with the reduced-order model (ROM).The ROM is not the end goal; it only serves to reduce the number of output variables that will be emulated by the final framework, i.e. the RNN-POD surrogate.
Projection-based model-order-reduction (MOR) projects the original n u kinematic variables, u , onto the n m variables of the ROM ( α ) using n m modes according to: where φ i of length n u denotes the i th mode and scalar α i denotes its associated modal coefficient that is to be computed by the ROM (whereas the final RNN-POD surrogate will rapidly emulate these modal coefficients).
and α collect all the modes and the modal coefficients in an n u × n m snapshot matrix and a column of length n m , respectively.
In the method of proper-orthogonal-decomposition (POD), modes are the left-singular (orthonormal) vectors corresponding to the largest singular values of an n u × n s matrix that stores all relevant n s training solu- tions (as columns) computed by the DNS.In other words, the ROM needs training data predicted by the DNS, whereas the final RNN-POD surrogate needs training data predicted by the ROM.However, the ROM requires much less training data than the RNN-POD surrogate.
In RVE simulations, the modes are not directly extracted from the training solutions (computed by the DNS).Rather, each solution field is first decomposed into a homogeneous (i.e.affine) field, ū , and a fluctuating field, ũ: Because the macroscopic deformation is applied to the boundaries of the RVE, homogeneous field ū is known.The modes are thus only used to interpolate the fluctuating part of the displacement field, yielding the following expression: where ω of length 3 (in 2D) is completely dictated by the (known) components of U M (three components in 2D simulations) and (of size n u × 3 in 2D) homogeneously interpolates the displacement components.
According to Eqs. ( 17) and ( 18), the homogeneous deformations are first subtracted from the training solutions, and SVD is applied to the resulting (i.e.fluctuating) displacement field.
The decomposition generates modes that are themselves periodic and hence, periodicity does not need to be actively enforced in the ROM (nor in the RNN-POD surrogate).Consequently, the system of linear equations that needs to be solved in terms of update dα in the online POD reduced-order-model reads: where dω = 0 during Newton's method.
A relevant note must be made with regards to the use of projection-based MOR for elastoplastic simulations.The reason is that POD-based MOR requires only a few modes (five to twenty) for (non-linear) elastic simulations to obtain a good accuracy, but a substantially larger number of modes for elastoplastic simulations.In fact, studies as those of Zhang et al. 30 and Daniel et al. 31 adaptively change the modes during the course of an elastoplastic ROM simulation in order to reduce the number of modes (and hence, the number of degrees of freedom; the n m coefficients in α).It is furthermore worth noting that the discussed ROM strategy only reduces the number of kinematic variables relative to the DNS (from n u variables in the DNS to n m variables in the ROM).This entails that solving the systems of linear equations is faster (cf.Eq. ( 15) and Eq. ( 19)).However, the computational efforts for constructing the systems of linear equations remains the same.In order to accelerate the construction of the systems, hyperreduction-type strategies may be utilized [37][38][39] , but our ROM does not exploit them.

Recurrent neural network
In this section, we discuss the recurrent neural network (RNN) that rapidly emulates the model coefficients, α , for each increment.This circumvents the iterative process required for the conventional ROM of the previous section (see Eq. ( 19)).However, the conventional ROM must be performed many times beforehand, in order to generate the training data for the RNN.
Network architecture.The overall architecture of the employed RNN is graphically presented in Fig. 2 and is the same as that of Wu et al. 9 .It consists of a gated recurrent unit (GRU), sandwiched between two feed-forward neural networks (FNNs).The input, i.e. the three independent components of U M t at time t, is fed to FNN I , which transforms it into u ′ .The output of FNN I , u ′ , together with the hidden states of the previous increment, h t−1 , is transformed by the GRU into output v ′ and the next hidden states, h t .In turn, FNN O takes the GRU's output, v ′ and transforms it into the modal coefficients of the current time increment, i.e. α t .
Thus, the RNN outputs are evaluated from the problem inputs ( U M t ) and from the hidden states of the previous increment ( h t−1 ).These hidden states are updated by the RNN for the evaluation of the next increment ( h t ).The user can choose the number of hidden states that the RNN emulates during its training.These auxiliary quantities, i.e. the hidden states, play the role of internal state variables in constitutive models although they lack a physical meaning.Instead, they correspond to a reduced version of the internal variables of the physical problem.
Relevant to realize is that the two sets of output of the GRU, i.e. the output of interest to the current increment ( v ′ ) and the updated hidden states necessary for the next increment ( h t ), are one and the same, i.e. v ′ = h t .This entails that the GRU must emulate the same quantity in such a way that (i) quantity v ′ is accurate for the current increment, and (ii) the same quantity ( v ′ = h t ) is accurate for the next increment.It may be clear that this limits the capabilities of the GRU to some extent.For instance, the RNN of Liu et al. 15 and Bhattacharya et al. 16 uses two independent quantities: one for the output of interest to the current increment (the equivalent to v ′ ) and one for the hidden states of interest to the next increment (the equivalent to h t , v ′ � = h t ).Nevertheless, the GRU as exploited in the current contribution was demonstrated to be sufficiently accurate to emulate macroscale stress-deformation relations by Wu et al. 9 .
The architecture of both FNNs is rather straightforward, because they both consist of a single layer of neurons.This is depicted in Fig. 3 for FNN I .All neurons use the leaky ReLU activation function.This entails that the i th output of FNN I , u ′ i , reads: where the three components of U M t are denoted by (U M t ) 1 , (U M t ) 2 and (U M t ) 3 .W 1i , W 2i and W 3i denote weights that are to be identified (i.e.'trained' , 'learned').Bias b i must also be identified.The leaky ReLU activation func- tion reads: The relation between the entire output of FNN I and its entire input may then be written as: with: The leaky ReLU activation function in Eq. ( 22) must be applied to each component of the column independently.
Figure 2. Schematic of the RNN.The three components of U M t (at time increment t) are the input variables and the n m coefficients of α t (at time increment t) are the output variables.FNN I and FNN O are each a single layer feed-forward neural network.The GRU's input is the output of FNN I (denoted by u ′ ), as well as the hidden states the GRU has predicted for the previous time increment, h t−1 (i.e. at time t − 1 ).The GRU's primary output, v ′ , serves as the input for FFN O , as well as updated hidden states, h t , that the GRU will use as input for the next time increment, i.e. at t + 1. FNN I 's output, u ′ , is the input for the GRU.The GRU, graphically presented in Fig. 4, transforms the input together with the previous hidden states, h t−1 , into output v ′ , as well as into the next hidden states, h t .We stress again that the GRU's output and its updated hidden states are the one and same ( v ′ = h t ).
Fig. 4 illustrates that the GRU has a so-called 'reset gate' (blue) and an 'update gate' (green).The relations between the two outputs of the update gate and the input and the previous hidden states read: where matrices W and column b A denote matrices with weights and a column with biases, respectively.All their components must be trained/identified/learned.Furthermore, σ (•) denotes the sigmoid activation function:  On the other hand, the relation between the input, the previous hidden states and the output of the reset gate reads: where tanh(•) denotes the hyperbolic tangent activation function: which means for Eq. ( 27) that it must be applied to each component of the column independently.
Finally, the GRU's output of current interest, which equals the updated hidden states, reads: The new hidden states, h t , are stored for the next time that the RNN is called upon (since they play a similar role as internal variables in elastoplastic simulations).On the other hand, the GRU's output serves as the input of FNN O , which only contains a single layer of neurons (employing the leaky ReLU activation function).Consequently, FNN O must contain the same number of neurons as the number of modal coefficients in α , i.e. n m .(Note that n O , denoting the number of hidden states as well as the number of output variables of the GRU, is independent of n m .) Learning phase.Now the RNN is formulated, its parameters (i.e. the weights and biases) must be identified (i.e.'learned' or 'trained' in NN terminology).The number of parameters is substantial, because FNN I and FNN O each comes with a weight matrix and a column of biases (i.e.4n I for FNN I and (n O + 1)n m for FNN O ).Furthermore, the GRU comes with six weight matrices and three columns with biases (i.e.(3n The total number of parameters is thus 4n I + n m + (n m + 3n O + 3n I + 3)n O (where once again, n I denotes the number of neurons in FNN I , n O the number of hidden states, which equals the number of neurons in FNN O , and n m the number of coefficients, i.e. the length of α).
For the ease of the notation, we store all parameters of the RNN in column p .First, the conventional reduced- order-model is performed numerous times and all coefficients predicted for each time increment, for each training simulation, are stored.Both input variables ( U M ) and output variables ( α ) are then normalised in view of training the RNN.We will denote the normalized output variables as αi t , where superscript i indicates the training simulation number and subscript t refers to the time increment number.The RNN must replicate these coefficients.To this purpose, one could minimize the following least square error: where n s denotes the number of training simulations performed with the conventional ROM, n incr denotes the number of increments of each training simulation, and α i t denotes the RNN's prediction of the coefficients.Furthermore, � • � denotes the L 2 -norm.
We use the ADAM optimizer 40 to identify (i.e.'train' or 'learn') the weights and biases (stored in p ). ADAM combines stochastic gradient descent with root mean square propagation.The stochastic gradient descent part of the optimizer requires the gradient of L with respect to parameters p .Because the RNN's prediction depends directly on parameters p , but also indirectly, via the dependency on the previous hidden states ( h i t−1 ), this dependency must be incorporated in the gradient.
The stochastic part of the gradient descent algorithm implies that it does not consider the entire least squares error every iteration of the optimization.Instead, only a part of the n s training simulations are considered for each iteration of the optimization, which is called a batch.Thus, if ten batches are considered, it takes ten iterations to consider all n s training simulations in the optimization (which is called an epoch in NN terminology).Because of the dependency on the hidden states, the data of each training simulation (i.e. a sequence of n incr columns α ) must be present in one batch in its entirety.In fact, the correct sequence of the incremental data must be incorporated in the training.
Thus, the (stochastic) gradient determines in which direction the RNN's parameter values are updated.On the other hand, the step size ('learning rate' in NN terminology), which governs how much the RNN's parameter values are updated in the direction of the (stochastic) gradient, are determined by the first moment and second moment estimates of the root mean square part of ADAM.( 26) coefficients of all increments are stored.Furthermore, the sequence of the modal coefficients for all increments must be accounted for in the training.Finally, we summarize the hyperparameters for the ADAM optimizer.Each mini-batch considers approximately 1% of the training data.In the ADAM optimizer, we set the initial learning rate to 0.001 (i.e.'γ ' in ADAM terminology) and the exponential decays for the first order and second order estimates to 0.9 and 0.999, respectively ('β 1 ' and ' β 2 ' in ADAM terminology, respectively).A cut-off value of 10 −8 is used to prevent divisions by zero ('ǫ ' in ADAM terminology).The initial values for the weights and biases are randomly generated and all hidden states are initially set to -1 (as suggested by Wu et al. 9 ).

RNN results
We now investigate the accuracy of the RNN for (i) different numbers of layers in FFN O , (ii) different numbers of neurons in FFN I , and (iii) different numbers of hidden states (where the number of hidden states equals the number of neurons in each layer of FFN O ).We investigate 64 possible configurations.The number of layers in FFN O is varied between 1, 2, 3 and 4. The number of neurons in FFN I is varied between 100, 250, 500 and 1,000.The number of hidden states (and hence, the number of neurons in each layer of FFN O ) is varied 400, 800, 1,600 and 2,500.
We investigate the accuracy of the 64 RNNs after they are trained for 60,000 epochs.After 60,000 epochs, the RNN with the smallest loss function is then assumed to be the best possible RNN from the available options.We train this final RNN for a total of 450,000 epochs after which we have observed that its loss function does not reduce further.
The loss functions of the 64 possible RNNs after 60,000 epochs are presented as bar charts in Fig. 6.The RNN with the smallest loss function after 60,000 epochs is the one with 1 layer in FFN O , 100 neurons in FFN I and 1,600 hidden states (and hence 1,600 neurons in FFN O ).

Mechanical results
Although the previous section has demonstrated the capabilities of the RNN, this does not automatically entail that the overall accuracy of the emulator is the same.In this section, the mechanical results of the RNN-POD surrogate are investigated and compared to those of the DNS and the conventional ROM.The comparison focuses on the macroscopic stress-deformation responses, but for the cyclic loading path, we also compare the microstructural field of the plastic multiplier at the end of the simulation/emulation ( in "POD-based model-order-reduction"section).
The comparison is divided in four parts: we first focus on the accuracy for a cyclic loading verification simulation/emulation, then on the accuracy for a random loading verification simulation/emultion, subsequently on the accuracy for three simulations/emulations for non-proportional load paths, and finally we compare the computational times.Cyclic loading.Fig. 10 shows the macroscale 1st PK stress components for a cyclic verification simulation/ emulation (that was not used to train the ROM, nor the RNN).We can see that the MOR results match those of the DNS fairly accurately (albeit not perfectly), indicating that the number of 100 modes is sufficiently large.The results of the RNN-POD surrogate also match those of the DNS and those of the MOR fairly accurately.
However, the problem is that the stress-deformation results emulated by the RNN-POD surrogate clearly fluctuate around the results of the ROM.This lack of smoothness, i.e. the fact that RNN-POD surrogate does not predict ∂P M ∂F M accurately, is expected to substantially affect the convergence behaviour of a true multiscale simulation.
We nevertheless also study the accuracy of the RNN-POD surrogate for a microstructural field.To this purpose, the final fields of the plastic multiplier ( ) predicted by the different approaches for the same cyclic load case are presented in Fig. 11.The diagrams clearly show that the difference of the field emulated by the surrogate stems from the inaccuracy of the modal representation of the displacement field, and not from the inaccuracy of the RNN.Random loading.We continue with results for random loading scenarios.Fig. 12 presents the macroscale, in-plane 1st PK stress components for a random loading verification simulation.Comparing the different predic- tions with the naked eye demonstrates that the surrogate is fairly accurate for the stringent scenarios of random loading.
However, this changes if we introduce an error for the macroscopic stress and compare it to that of the cyclic loading verification simulation.To this purpose, we introduce the following average relative error for the inplane stress components: where the jjth stress component for the ith increment computed by the DNS is denoted by Pi jj and the same component resulting from the conventional ROM or the RNN-POD surrogate is denoted by Pi jj .
(31) e =  The logarithm of this error is presented in Fig. 13 as the first two bars on the left for the cyclic loading verification simulation (blue: ROM, red: RNN-POD surrogate) and in the subsequent two bars for the random loading verification simulation.
Several observations can be made.First, the conventional ROM for the random case is significantly less accurate than for the cyclic case.In fact, the conventional ROM for the random case is substantially less accurate than the RNN-POD emulator for the cyclic case.The additional increase in error of the RNN-POD emulator is of the same order of magnitude compared to the conventional ROM.
The results discussed so far have focused on cyclic loading and random loading; the scenarios for which the surrogate was trained.At this moment, we can already conclude that the RNN-POD surrogate cannot be used in true multiscale simulations.Although the surrogate is sufficiently accurate for cyclic loading, the resulting  smoothness of the resulting macroscopic stress-deformation response is insufficient to drive a multiscale simulation.On top of that, the macroscpic stress-deformation response for random loading is simply inaccurate.
Nevertheless, it is interesting to investigate scenarios for which the RNN is not trained.The reason is that we have assumed that training for random loading and cylic loading will make the RNN capable to emulate any scenario (as previously demonstrated for an RNN that only captures the macroscale stress-deformation relation 9 ).
Non-proportional loading.We now focus on the three non-proportional load paths in Fig. 14.The homogenized stress components predicted for these load paths are presented in Figs. 15 for the black load path, Fig. 16 for the yellow load path, and Fig. 17 for the green load path.We first concentrate on Fig. 15 which presents the 1st PK stress components for the black load path in Fig. 14.It is important to realize that the black load path was discretized using 1,000 increments.Thus, the size of an increment of the black load path, which looks somewhat similar to the cyclic loading simulations used for training, is substantially larger than that of the cyclic loading training simulations, but substantially shorter than that of the random loading training simulations.
The homogenized stress components in Fig. 15 emulated by the RNN-POD framework are clearly substantially worse than those of the investigated cyclic loading verification simulation (Fig. 10).This is also clear from Fig. 13, in which the 5 th (conventional ROM) and 6 th bar (RNN-POD surrogate) present the logarithm of the error of the stress (Eq.( 31)).The error is however substantially smaller than for the random loading verification simulation.
One could argue that the emulated stress-deformation curves (Fig. 15) are so poor because the incremental stepsize used to discretize the associated load path (black in Fig. 14) is different from the incremental stepsizes of the cyclic loading training simulations and of the random loading training simulations.For this reason, the yellow and green load paths in Fig. 14 are made such that when discretized, their incremental stepsize is exactly the same as that of the random loading training simulations (i.e.0.02 in the U M 11 − U M 12 − U M 22 space).However, Figs. 16 and 17 (together with four bars on the right in Fig. 13) demonstrate that the homogenized stress components of the RNN-POD framework for the yellow and green load paths are even worse than those predicted for the black load path.Thus, restricting the stepsize to that of the random loading training simulations (computed with the conventional ROM) has no positive effect on the accuracy of the RNN-POD emulator.
This leads to the conclusion that training the surrogate for random loading and cyclic loading scenarios (two scenarios on opposite sides of the spectrum of possible loading scenarios) does not mean that other loading scenarios will be adequately emulated by the surrogate.As mentioned before, this contrasts previous observations 9 .
Acceleration.So far we have only discussed accuracy and have ignored speed up.Table 1 summarizes the computational times of the different types of simulations/emulations.Of course, these numbers highly depend on the implementations, hardware, spatial dimensions, type of elements employed in the DNS, number of elements employed in the DNS and more.Consequently, substantial higher speed ups may be realized than reported in Table 1.
Though the data preparation and training the RNN required a total of two weeks of computational time, Table 1 shows that the surrogate is approximately 100 times faster than the DNS and 22 times faster than the conventional MOR for random loading.For cyclic loading on the other hand, the surrogate is only 13 times faster than the DNS, whilst the conventional ROM is hardly faster than the DNS.The difference in time savings of the surrogate for cyclic and random loading are because the DNS and conventional MOR require more iterations for random loading than for cyclic loading.The reasons are that (1) the loading paths for cyclic loading are substantially shorter than those for random loading whilst the same number of increments is employed, and (2) the previous plastic state is assumed at the start of each increment of a cyclic loading simulation (i.e.DNS and conventional ROM) in order to increase the speed of the cyclic loading simulations.

Conclusion
An RNN-POD surrogate model is investigated that emulates the displacement field of each increment of a rateindependent elastoplastic simulation.A recurrent neural network (RNN) in the surrogate is responsible for the emulation.In order to reduce the number of output variables of the RNN, proper-orthogonal-decomposition (POD) is exploited to efficiently reconstruct each displacement field as the sum (of a relatively small number) of orthogonal modes, each multiplied by a single coefficient.Consequently, the RNN only has to emulate the coefficient of each mode (instead of the displacement vector of each finite element node).
We have applied the framework to an elastoplastic RVE.Whereas many neural networks 1,8,9,14 aim to rapidly emulate the homogenized stress-deformation of an RVE (which drives multiscale simulations), they lose all microstructural information.The RNN-POD surrogate has the potential to both drive multiscale simulations, as well as to recover any microstructural field of interest.
Because the surrogate emulates the displacement field in the entire RVE domain, the non-linear governing equations of the RVE simulations do not need to be constructed, nor solved.Instead, only the new plastic variables must be computed once per increment (together with any other quantity of interest).The time savings are thus substantial and we report a maximum acceleration factor of 100 relative to the direct numerical simulations.
Instead of a conventional feed-forward neural network (FNN), an RNN is used in the surrogate.The reason is that RNNs predict both the primary output of interest, as well as hidden states, which affect the RNN's next prediction.Thus, the role of hidden states in RNNs imitates the role of plastic variables in elastoplastic finite element simulations.
Although we report acceleration factors of up to 100, the accuracy of the surrogate leaves much to be desired.The reconstruction of the displacement field using the POD modes definitely reduces the accuracy (in terms of both macroscale stress-deformation relation, as well as the recovery of a microstructural field), but the macroscale stress-deformation relations predicted by the conventional ROM seem to be just as smooth as those of the DNS.Hence, we believe that the inaccuracy introduced by the POD modes would not affect the convergence of a multiscale simulation.
On the other hand, the RNN introduces both a substantial inaccuracy in the macroscale stress-deformation curves, as well as a substantial reduction of the smoothness of the macroscale stress-deformation curves.Only the macroscopic stress-deformation curves for the cyclic load case are accurate (but lack smoothness).We want to emphasize that this is in contrast to the results of Wu et al. 9 , in which the same RNN was used to only emulate the macroscale stress-deformation curves.Apparently, an RNN that mimics an elastoplastic constitutive response may not be sufficient to emulate the modal coefficients of an elastoplastic ROM.
In conclusion, we can state that the RNN-POD surrogate, as presented here, cannot emulate the microstructural fields sufficiently accurate to truly drive multiscale simulations.Potential avenues to accomplish this aim are: • A wider variety of load paths and of incremental step sizes of the macroscale deformation may be incorpo- rated in the training data of the RNN.• Another type of RNN may increase the accuracy.For instance, the RNN-POD framework of Im et al. 26 uses two long-short-term-memory (LSTM) units in parallel (whereas we have only used a single GRU) Another possible avenue are the RNNs of Liu et al. 15 and Bhattacharya et al. 16 .In contrast to our GRU and the LSTM of Liu et al. 26 in which the quantity of current interest and the updated hidden states are one and the same, the output of current interest and the updated hidden states are truly separate in the RNNs of Liu et al. 15 and Bhattacharya et al. 16 , potentially increasing the accuracy.• In order to drive multiscale simulations, the sensitivity of the macroscale stress-deformation relation must be sufficiently smooth.The presented results have revealed that the RNN emulation (and not the POD construction of the displacement field in the ROM) is responsible for the lack of smoothness.A potential avenue may therefore be to not only train the RNN to emulate the displacement field itself, but also to train the RNN's sensitivity according to the sensitivity of the macroscale stress-deformation relation.

Figure 3 .
Figure 3. Schematic of FNN I , that consists of a single layer.The three components of U M t (at time increment t) are the input variables and the n I variables in u ′ are the output variables.

Figure 4 .
Figure 4. Schematic of the gated recurrent unit (GRU).Input u ′ is the output of FNN I and the GRU's output, v ′ , is the input of FNN O .Output v ′ equals the hidden states for the next calculation of the RNN, h t .The reset gate is presented in blue (with output R ) and the update gate in green (with outputs U and 1− U).

Figure 6 .
Figure 6.The loss function for the verification data after training for 60,000 epochs for three neural network parameters: (i) the number of hidden layers in the FFN O (top-left: 1 layer, top-right: 2 layers, bottom-left: 3 layers, bottom-right: 4 layers), (ii) the number of neurons in the single layer of FFN I ('N') and (iii) the number of hidden states (the different colors/patterns).The size of the bar corresponds to the value of the loss function (i.e. a large bar corresponds to a large value of the loss function).Please note that some loss functions are so small that they are difficult to observe.

Figure 7 .
Figure 7.The evolution of the loss function for both training data (blue) and verification data (red) for the first 50,000 epochs (a) and the last 75,000 epochs for which it was trained (b).The loss function for the verification data (red) clearly does not reduce during the last 75,000 epochs (b).

Figure 8 .
Figure 8. Cyclic loading verification simulations: Some RNN predictions (crosses) and the actual values (lines).The colors distinguish the four verification simulations.

Figure 9 .
Figure 9. Random loading verification simulations: Some RNN predictions (crosses) and the actual values (lines).The colors distinguish two verification simulations.

Figure 10 .
Figure 10.Macroscale, in-plane 1 st PK stress components as functions of the deformation for a cyclic loading verification simulation.Black: DNS results, dotted blue: conventional ROM results, dashed red: RNN-POD surrogate results.

,Figure 11 .
Figure 11.The plastic multiplier ( ) for a cyclic loading verification simulation.Top-left: the DNS results, top-right: the difference between the ROM results and the DNS results, bottom-left: the difference between the surrogate results and the ROM results, bottom-right; the difference between the surrogate results and the DNS results.

Figure 12 .
Figure 12.Components of the macroscale 1st Piola-Kirchhoff stress as functions of the number of increments for a random loading verification simulation.

Figure 13 .
Figure 13.The logarithm of the error of Eq. (31) for the cyclic loading verification simulation (left two bars), the random loading verification simulation (3th and 4 th bar), the black non-proportional load path in Fig. 14 (5th and 6 th bar), the yellow non-proportional load path in Fig. 14 (7th and 8 th bar) and the green non- proportional load path in Fig. 14 (9th and 10th bar).The blue bars correspond to the errors of the conventional ROM and the red ones to those of the RNN-POD surrogate.

Figure 14 .
Figure 14.Three non-proportional load paths.The black path is discretized using 1,000 increment, whereas the other two are discretized such that each stepsize is the same as in the random loading simulations used for training.

Figure 15 .
Figure 15.Components of the macroscale 1st Piola-Kirchhoff stress as functions of the deformation for the non-proportional cyclic load path in black in Fig. 14.

Figure 16 .
Figure 16.Components of the macroscale 1 st Piola-Kirchhoff stress as functions of the deformation for the non-proportional cyclic load path in yellow in Fig. 14.

Figure 17 .
Figure 17.Components of the macroscale 1st Piola-Kirchhoff stress as functions of the deformation for the non-proportional cyclic load path in green in Fig. 14.

Figure 1 .
The meshed RVE with the elastoplastic matrix in green and the stiff, elastic particles in red.

Table 1 .
Average computational times for data preparation and the online verification simulations/emulations.