Data-driven model order reduction with surrogate elements for transient simulations

The Finite Element Method (FEM) is a widely used numerical simulation method in almost all scientific and engineering disciplines. Despite its capacity to address the most complex problems in any physical domain, the method has the disadvantage of requiring significant computational resources, particularly for systems with a high number of degrees of freedom (dof).

Achieving the goal of reducing computational costs involves several approaches. Recent research focuses on advancements in numerical solvers (Sato et al., 2023; Chen et al., 2021). Another promising area is the replacement of complete models with surrogates as described by Kudela and Matousek (2022) or Franke and Wagner (2024). While this significantly reduces computational cost, these surrogate models are typically limited to specific configurations. Another important area is called model order reduction (ROM). Methods such as proper orthogonal decomposition (POD), the Krylov subspace method or the modal reduction method are used (Bathe, 2007; Wriggers, 2008; Ingrid and Rottner, 1999; Schilders et al., 2008). All these methods have in common that the full system is approximated by a system with less dofs and after solving it is re-expanded to the full model. Furthermore, there are purely non-intrusive MOR techniques for finite element models. These include projection-based approaches (Mahdiabadi et al., 2021; Le Guennec et al., 2018), machine learning-based methods (Kneifl et al., 2024) or hybrid approaches combining both (Czech et al., 2022; Fresca et al., 2022).

Additionally, the method of dynamic substructuring with superelements, as described in Cammarata et al. (2019), Allen et al. (2020) can be employed. Here, dofs that are not relevant to other parts of the model are removed from the model by static condensation, see Wilson (1976). Similarly, the Guyan reduction addresses the same issue, as detailed by Guyan (1965). These substructures can now be solved once in an offline process and then reused on numerous occasions. Such a pre-solved part is called a superelement. This method is applicable to circular symmetric geometries where the connection to the other parts of the model or the boundary condition is known and unaltered. Consequently, it is not possible to adapt a superelement to other models, which limits their flexibility. This is intended to be circumvented by the proposed method of using machine learning as described in the following sections. It is important to note that substructuring does not introduce approximations beyond those inherent to the FE formulation itself, in contrast to the majority of machine learning approaches.

Furthermore, the multiscale FEM can also be employed to reduce computational costs, as demonstrated by Chung et al. (2023). In multiscale finite element simulations, a problem is initially solved based on a coarse discretization. Each component or element of the coarse solution is simulated again using a fine discretization. In this context, machine learning methods are described in Nguyen et al. (2023), Koeppe et al. (2020). In many instances, the calculation of the coarse discretization is not performed with regular finite elements, moreover, a method known as the unit displacement method is used (Pourazarm et al., 2011; Koeppe et al., 2020). A given part of the model in the macro scale is defined by a number of boundary nodes. For each of these nodes, an incremental element nodal force is evaluated for small unit deformations. This approach deviates from conventional FEM because shared shape functions across multiple elements, which are required to form the macroelement, are not available. In essence, the computation of the macroelement is analogous to the calculation of Craig–Bamptons interface modes, as elucidated by Craig and Bampton (1968). In recent research work, the fine discretization is modelled by artificial neural networks (ANN), as summarized in the review work of Bishara et al. (2023), or deep learning as demonstrated in Koeppe et al. (2020), Deng et al. (2024). However, at first the stiffness properties of the macroelements must be characterized. Huang et al. (2023) present a method for training an ANN to predict stiffness matrices for a coarse mesh and corresponding shape functions for substructure modeling. The deep NN uses the Young’s modulus at each node of the coarse mesh as input for both the stiffness matrices and shape functions. As a result, these models are adaptable to various materials within the offline-trained range of Young’s modulus. The stiffness matrix and shape functions are subsequently mapped to the specified material property. Beyond the context of substructuring, but also constructing stiffness matrices of quadrilateral finite elements, Jung et al. (2020) use deep learning and a large amount of strain data to corresponding material properties, displacements and geometries. Capuano and Rimoli (2018) employ machine learning techniques to establish a relationship between the element state and its forces in multiscale simulations. This approach eliminates the need for element integration. The stiffness matrix is generated by applying unit loads to the nodes of the coarse mesh, identical to standard multiscale simulations. However, the stiffness matrix is then calculated by machine learning.

The identification of stiffness matrices falls within the domain of inverse problems. In accordance with the classification proposed by Beck and Woodbury (1998), inverse problems can be divided into two categories, measurement problems, which pertain to parameter identification (Chamekh et al., 2009; Römer et al., 2024) and design problems, as exemplified by Fachinotti et al. (2020). Design problems are formulated for cases where the underlying data is explicitly and accurately defined. In the following we focus on design problems, according to our use case. Solving inverse problems of partial differential equations (PDE) with (physically informed) NN is a topic that has been extensively discussed in the literature (Zhang et al., 2019; Berg and Nyström, 2021; Badia et al., 2024; Raissi et al., 2019). With regard to the calculation of stiffness matrices, it is possible for an ANN to work at different level. Oishi and Yagawa (2017) propose a replacement of the standard element integration routine. Still working with Gauss-Legendre quadrature, the weights and integration points are adjusted to an optimum for different cases. These changes are controlled by deep learning networks. From a more comprehensive perspective, the data-driven computational mechanics paradigm approach (Kirchdoerfer and Ortiz, 2016; Nguyen et al., 2020), also searches a stiffness matrix. All elements within a simulation are fit to existing (experimental) stress strain curves while satisfying compatibility and equilibrium due to the Lagrange multiplier method. Hence, the measured material data can be utilized directly in computational models. The data-driven computational mechanics paradigm also applies to multiscale simulations, as proposed by Karapiperis et al. (2021). The identification of stiffness matrices using ANNs as an inverse problem is described in Meethal et al. (2023). In this context, the stiffness matrix of a physical model includes unknown components. It is therefore divided into a known part, an unknown part and two corresponding unknown cross-coupling components. The unknown portion of the stiffness matrix is determined using an ANN. Optimization of the unknown entries is achieved using the known force vector and displacement vector. However, it is important to note that the resulting predicted stiffness matrix is not compatible with other models and cannot be directly combined with other stiffness matrices.

In this work we propose a method, where the ideas behind the substructuring and the multiscale method are combined. Macroelements or superelements are created for repeated parts of a model, which can be used in a more flexible manner than substructuring or static condensation in different models. For this purpose, stiffness matrices of repeated parts of a large model are generated in a reduced manner. This implies that only a subset of the original number of dofs is utilized in the surrogate elements stiffness matrix. A significant focus is placed on positive semi-definite matrices, which align with conventional finite element. For this, we introduce artificial degrees of freedom (adofs). To reduce the training effort of the surrogate element matrices, the identification is conducted in silico on models with a low number of dofs. Subsequently, the surrogate element is used to model the repeating parts of large models as a forward problem. From a connectivity perspective, the stiffness matrix of the surrogate element can be treated in the same manner as a standard finite element, allowing it to be added to the dofs of adjacent elements. Consequently, no explicit coupling terms are required to connect two regions. The surrogate element behaves similar to macroelements known from multiscale FEM, wherein information on the region described by the element is only provided at the boundary. Nevertheless, the structural properties of the entire element are considered with high accuracy. As a benefit, calculation time is considerably reduced.

The proposed method is related to the work of Huang et al. (2023) and the work of Capuano and Rimoli (2018). Regarding the work of Huang et al. (2023), the ANN uses the material parameter as input, whereas our approach does not require input, as only the bias is trained. Hence, our architecture of the ANN is less complex, which benefits the training and prediction costs. A limitation of our approach is that a separate model must be trained for each material. This restriction applies unless the ANN is provided with material-specific inputs, which parametrize the selected material model of the training data. But in scenarios where the material is predetermined, this issue is not a significant concern. In addition, our approach is applied to transient simulations. Regarding Capuano and Rimoli (2018), the main differences compared to the proposed method here, lies in the application to solid elements and the direct training of the surrogate elements on simulation data. Hence, in our approach no additional load cases for unit forces or displacements methods are necessary. Furthermore, in both works all elements within the simulations are replaced by surrogate elements, which are named coarse elements and smart elements in their publications. Novel contributions of the proposed method are the combination of the initial discretization with surrogate elements within one simulation and the introduction of artificial degrees of freedom.

The article is organized as follows: The initial section presents the fundamental principles of FEM and surrogate modeling. Subsequently, in Section 2.2, a comprehensive analysis of the ANN responsible for generating the matrices is presented. Section 2.3 describes in detail the concept of adofs. This is followed by Section 2.4 and 2.5, which present approaches for the generation of the surrogate elements stiffness- and mass matrices. Section 2.6 and 2.7 describe the data-generation routines and the training process, respectively. Later, in Section 3, at first results for the low dimensional training examples are presented. Thereafter, the application of the surrogate elements to large examples, which utilize several surrogate elements for modeling the problems, are discussed. There, we also focus on the application of the surrogate elements to complex dynamic structures. Finally, a conclusion is drawn in Section 4.

In order to facilitate the surrogate modeling of mechanical structures, we propose the introduction of surrogate elements, which are described by stiffness matrices and mass matrices that depend on the adjacent elements. The following section will provide a detailed description of the surrogate elements. Furthermore, the architecture and training of the ANN generating the surrogate elements are discussed.

The discretized linear form of the partial differential equation (PDE) describing the motion of a mechanical system is given by

see Hughes (1987). Thereby, $u(t)∈RN×1$ and its derivatives contain the nodal displacements, velocities and accelerations. The variable N denotes the total number of dofs. Moreover, $f(t)∈RN×1$ is the vector of the nodal forces. Further, $K∈RN×N$⁠, $C∈RN×N$ and $M∈RN×N$ denote the stiffness, damping and mass matrix, respectively. If the problem is of a static nature, C, M, $u̇(t)$ and $u¨(t)$ are neglected. Also the time dependency of u and f is neglected, when loads and boundary conditions are not time dependent. Thus, in the following static cases, they appear without (t). The quadratic and symmetric matrices are created by adding element matrices, for example $Ke∈RNe×Ne$ with N_e as number of dofs per element, taking into account node connectivity. Element matrices are built with numerical element integration of shape functions. If modeling of the stiffness matrix is correct, all K^e are positive semidefinite (psd). This is due to the fact that these matrices represent the strain energy stored in the structure, which is modelled by the elements. Given that an energy can only be positive or zero, it follows that elements which are not psd cannot model physical systems in a satisfactory manner. This implies that the eigenvalues of the matrices are either zero or positive. Since all K^e are psd, the global stiffness matrix K, which is built by summing all K^e according to their corresponding dofs, is also psd. Before the kinematic boundary conditions are applied, Equation (1) cannot be solved for u, because K is singular and hence, K⁻¹ cannot be computed. After elimination of supported dofs, the mechanical problem is kinematically determined. That means K is now positive definite and can be inverted.

It is not possible to interpret a surrogate element in the same way as a standard finite element. Moreover, it can be employed to replace several standard finite elements. The construction of a surrogate element does not necessitate the use of shape functions, material models, transformation to local coordinate systems or numerical integration. Consequently, a surrogate element is defined as a region within the matrices that model a specific portion of a domain Ω. The modeling is not of a pure physical manner, instead it is based on a data-driven representation.

To describe the fundamentals of this approach, it is necessary to introduce different regions. In essence, a specific part of the physical domain is substituted with a surrogate element. This region is called the Replaced Region and can be found in Figure 1 as Ω_RR.

The region which persists is called Preserved Region, Ω_PR. This region is fully modelled in the discretization process with standard finite elements. Hence, results of that region can be fully interpreted after solving the system of equations. In contrast to that, Ω_RR is not fully discretized. Only the interface between the two regions, Γ_IR, is discretized. Furthermore, we state that Γ_IR ⊂ Ω_RR and Γ_IR ⊂ Ω_PR and therefore, Γ_IR = Ω_RR ∩Ω_PR and Ω = Ω_RR ∪ Ω_PR. The selection of Ω_RR should be guided by two primary criteria. First, this region should not contain areas of interest, as its results will not be directly calculated. Second, Ω_RR should represent a structural component that recurs multiple times within the model to maximize the reduction in system size.

A schematic assembly of the stiffness matrix K from the sub-areas is shown in Figure 2. There, it is assumed that the non-zero coefficients in the global matrix are sorted diagonally dominant. The stiffness matrix K can be subdivided in blocks corresponding to Ω_PR, Ω_RR and Γ_IR. The dofs of the elements, which are not affected of the model reduction process, are part of K_PR. Elements, which are not modelled in the reduced stiffness matrix, are part of K_RR. Due to the compatibility, there are dofs, which are part of K_PR as well as part of K_RR, see K_IR. The entries of K_IR created from elements in Ω_PR are preserved. In contrast, if they originate from Ω_RR, they are deleted.

The reduced stiffness matrix, denoted as $K̃$⁠, can be built from K by keeping K_PR and replacing K_RR with K_SE, the stiffness matrix of the Surrogate Element. To achieve this, one has to add K_SE to the corresponding parts of the interface regions, which are indexed with the subscript i, K_IR,i, see Figure 3.

The dimension of K_SE has to be at least equal to the sum of K_IR,i dimensions, which is N_IR. Optionally, the dimension of K_SE can be increased by inserting adofs from the matrix K_AD. Once K_SE is generated to represent a specific K_RR, the region Ω_RR can not be modified thereafter. Consequently, always the same number of finite elements is replaced by the surrogate element. The reduced mass matrix $M̃$ is built in an identical manner, whereas $C̃$ is obtained by the Rayleigh damping approach. A more detailed description on the reduced matrices is given in Sections 2.2, 2.4 and 2.5.

With regard to Equation (1), the system of equations of the reduced model is described by

where $K̃∈Rn×n$⁠, $C̃∈Rn×n$ and $M̃∈Rn×n$ as well as $ũ(t)∈Rn×1$⁠, respectively the derivations and $f̃(t)∈Rn×1$ describe the reduced matrices and vectors with n as number of dofs of the reduced approximation.

Standard surrogate models reproduce complete systems and therefore, their application is limited to the very same systems. In contrast, the proposed surrogate elements require only the interface elements to fit to the trained system if a different model is considered. Next to Γ_IR, an arbitrary region can be modelled. Hence, the offline trained surrogate elements are independent of the boundary conditions applied to an arbitrary model. In order to achieve a MOR character, it is necessary that the dimension of K_SE is smaller than the dimension of K_RR, so that n < N.

The stiffness and mass matrices are independently modelled using two separate dense layers, which are designed without interconnections. This separation is necessitated by the significantly different scales and distinct properties of the two matrices. An overview of the ANN is shown in Figure 4.

First we will focus on the stiffness matrix. The method will be introduced to be theoretically capable of handling geometrically nonlinear elastic problems for future work, although it will be applied only to linear elastic problems in this article.

As an input, the displacement vector from the dofs in K_IR of the previous iteration, schematically denoted as u_IR(−1), is utilized. This vector is processed by a dense layer and an activation function, see Goodfellow et al. (2016) 

with $kSE∈Rnb×NSE$ where n_b denotes the batch size and N_SE is the number of the surrogate elements dofs. The weight-matrix of the dense layer for the stiffness matrix is labelled W_K, whereas the bias-vector is denoted b_K. Furthermore, the hyperbolic tangent activation function is defined by tanh = (e^x − e^−x)/(e^x + e^−x). This function is selected because it produces both negative and positive values, which is essential for accurately modeling the stiffness properties.

From this point onward, the method is applied specifically to linear elastic modeling. In this context, no previous iterations are required as the system is solved once using the undeformed geometry, under the assumption of small deformations. Consequently, the input data is set to u_IR(−1) = 0. In general, this represents the displacements from a previous iteration in nonlinear problems, denoted as (−1). Hence, even for linear transient simulation u_IR(−1) is zero. The input of the ANN is multiplied element-wise with the weights W_K. As a result, the calculation of k_SE is independent of the weights W_K and depends solely on the bias b_K. In summary, this reduces the problem to a straightforward optimization task, where a single dense layer is sufficient.

In a similar way, the surrogate element mass matrix is constructed. Unlike the stiffness matrix, the mass matrix only allows positive entries. Hence, Equation (3) is adapted to utilize a sigmoid activation function σ = 1/(1 + e^−x),

where W_M is the mass matrix dense layer weight matrix and b_M the bias vector.

A more detailed description is required for the FEM block of the architecture, which is shown on the right in Figure 4. In this instance, the ANN is informed, as the system of equations comprising the discretized model is solved. Initially, the stiffness matrix $KSE∈RNSE×NSE$ of the surrogate element has to be built. This is explained in detail in Section 2.4. Next, $K̃$ is assembled from K_SE and K_PR. Similarly, if the analysis is not of a static nature, the mass matrix is assembled. After that, the reduced force vector $f̃$ is built from the original force vector f by adjusting dimensions. If required for transient simulations, the reduced damping matrix $C̃$ is calculated following the Rayleigh damping approach, see Wagner (2022),

where α_d denotes the mass proportional damping factor and β_d denotes the stiffness proportional damping factor.

In the static case, where $K̃ũ=f̃$ has to be solved for $ũ$⁠, the Cholesky decomposition, is used, see Wagner (2022). We choose this specific method for solving the system, since every physical system modelled correctly with FEM has to be psd, which is ensured by successfully solving using the Cholesky decomposition. Due to this, solving with a multifrontal method is also possible. Furthermore, another widely known solving method, the conjugate gradient method, also requires the use of psd matrices. While iterative methods like the preconditioned conjugate gradient and generalized minimal residual methods do not require positive-definite matrices, solution techniques relying on positive-definite matrices remain central to finite element simulations. This highlights the importance of our focus on psd surrogate element matrices. For transient problems, Equation (2) is frequently integrated over time using an unconditionally stable implicit Newmark Beta scheme; see Hughes (1987).

The surrogate element is trained on simple models, which contain only one surrogate element, see Section 2.6 and 2.7. Nevertheless, in the inference phase it should be possible to place the surrogate element several times in the same model in different orientations and locations. Since training of each orientation is not practical, the surrogate element is transformed from the training orientation to the desired orientation. For this, we use Euler angles ψ, θ, ϕ to rotate around the z, x′ and z″ axis with the rotation matrices R:

Equation (6) is performed on K_{SE, dof}, which denotes the parts of K_SE separated by the assignment of dofs in K_SE to their nodes. For a better readability, in the following the oriented matrices, K_{SE, dof, oriented}, are named K_SE. Equation (6) applies identically to M_SE.

The last step in the model is only necessary for the training process. From the vector of displacements, $ũ$⁠, only those dofs that originate from Ω_PR are extracted for the loss calculation. The forward propagation of the ANN is described in detail in Algorithm 1.

As depicted in Figure 3, it is possible to incorporate adofs as K_AD. Consequently, the number of instances of K_SE is increased, thereby enhancing the ability to adequately model Ω_RR. The placement of the adofs can be arbitrary; we chose rows at the right and consequently columns at the bottom of K_SE. Noteworthy, the dimension of the surrogates stiffness matrix is increased by the number of adofs N_AD, such that $KSE∈RNIR+NAD×NIR+NAD$⁠. As previously stated, K_SE is then added to an enlarged version of K_PR to construct $K̃$⁠. Due to this enlargement, there are no underlying entries of K_PR where K_AD is added. The initialization of K_AD with small random values can give rise to two potential issues:

1. The system of equations cannot be solved, despite sufficient boundary conditions.

2. The positive-definiteness of K_SE may be compromised when eliminating dofs associated with the boundary condition. On the beam example, this will happen when the boundary condition is only set on one side of the surrogate element. The introduction of adofs with a zero initialization will result in the formation of two non-connected regions, which renders the solution of the problem impossible.

Therefore, initialization of K_AD is done with

where $I∈RNAD×NAD$ is the identity matrix and f_i denotes an additional hyperparameter for the training process. Consequently, the unit of f_i is identical to the unit of K_AD, which depends on the unit system of the training dataset. It is recommended to chose values of the magnitude of the diagonal elements from K_PR for f_i.

Before training, adofs are not coupled with FE dofs. The coupling is created during the training process, whereby off-diagonal entries are assigned with values mostly ≠ 0 obtained from k_SE. These values are then assigned with a stiffness with respect to FE dofs. The adofs are not identifiable with a specific location within Ω, and therefore cannot be included in any subsequent interpretation of the results. This approach is analogous to the hierarchic FEM, where dofs assigned to higher order polynomials are also incapable of forming interpretable location and result pairs, as no additional nodes are created and the dofs are purely mathematical in nature, see Zienkiewicz et al. (1983).

The stiffness matrix of a surrogate element should generally possess the same properties as that of a standard finite element. Hence, after adding the surrogate element K_SE, $K̃$ still has to be psd. Furthermore, following the elimination of the boundary dofs, $K̃$ has to be positive-definite. Ensuring this property, it is possible to solve the system of equations with the Cholesky decomposition as mentioned in Section 2.2. Therefore, K_SE has to be symmetric and psd. In addition, to construct physical element stiffness matrices, the requirements can be enriched by positive non-zero main diagonal values, negative off-diagonal values and zero row sums. This implies that the sum of each row of the matrix must be zero.

The values in k_SE are out of the interval $[−1,1]≔{x∈R∣−1≤x≤1}$⁠. Hence, they have a different magnitude compared to the values in K_PR. To adjust that, k_SE is scaled by

with f_k as an dimensionless additional hyperparameter. Consequently, a minimum or maximum value, defined by $efk$⁠, is nonlinearly scaled by the dense layer. It is challenging to maintain accurate records of the units in question. Consequently, in this approach, the unit of k_{SE, scaled} is defined as a stiffness in the chosen unit system.

At first, we do not apply the method of artificial degrees of freedom. In order to ensure symmetry, only a lower triangular matrix of the surrogate element is formulated in the ANN. For this, the vector k_{SE, scaled} is resorted to build the lower triangular matrix $kSE, tri.∈RNSE×NSE$⁠. Furthermore, this approach resulted in reduced computational costs during the training and prediction phase of the ANN, in comparison to modeling a quadratic matrix. The number of output parameter of the dense layer, N_kSE, is reduced, such that $kSE∈RNkSE×1$ with

To ensure positive values on the main diagonal and symmetry, $KSEmm∈RNSE×NSE$ is calculated,

As a result, non-diagonal entries may still retain negative values. The unit of k_SE has to be expressed as the square root of a stiffness in the chosen unit system to result in a meaningful matrix K_SEmm. For psd, all eigenvalues λ of the unsupported $K̃$ have to be λ ≥ 0. The fulfilment of this requirement necessitates a specific approach, since K_SE is optimized by the ANN with regard to the displacements $ũ$⁠, rather than the positive eigenvalues.

In general, there are numerous methods for the optimization of eigenvalues as shown by Lewis (2003). The implementation of these methods into our workflow is challenging in certain instances. Hence, we seek a more straight-forward solution. In the field of comparing measurements, it is also necessary to deal with psd matrices. In this context, it is beneficial to search for the nearest psd matrix, as was done by Higham (1988). Unfortunately, this approach is not applicable in our problem, as the physical models are no longer adequately described following the search for a substitution matrix.

In order to generate psd matrices, the final requirement is assessed, which primarily should ensure force equilibrium. For this, the sum of the row-wise entries in an element stiffness matrix must be zero. We force the surrogate element to inherit that property. Hence, the resulting matrix corresponds to a Laplacian matrix, which is psd. For the row-wise zero sums, the final three columns of K_SE, and due to symmetry the final three rows too, are not modelled directly by k_SE. Moreover, $KSE∈RNIR×NIR$ is initialized as a zero matrix and K_SEmm is added to the left top corner of K_SE. This approach allows for the reduction of the number of ANN output parameter

The last columns are then populated with the negative column- and component-wise sums of K_SE. More precisely, in the three-dimensional representation, all dofs corresponding to the x, y and z-direction are added separately. Subsequently, in analogy, the last row is filled with the negative row-wise sum. Due to the summation process, the values of the last three rows and columns exhibit a different scale, especially entries in the lower right are most different. This distinctive configuration is prone to generating local anomalies in the results at these nodes.

By combining the method presented in Section 2.4 with adofs, this effect can be mitigated. It is possible to store the row- and column-wise sums on dofs that are not part of the interpretable output displacement vector. Therewith, the anormal values in the last three columns and rows can be transferred to dofs, which are not part of the physical displacements. Hence, they will not be evaluated. As an acceptable downside, the number of ANN output parameter increase once more, see Equation (9). Consequently, the surrogate elements matrix is now of the dimension $KSE∈RNIR+3×NIR+3$⁠, since N_AD = 3. Such a surrogate elements stiffness matrix is shown in Figure 5.

The procedure outlined in Section 2.2 and the steps described here are summarized in Algorithm 1, which illustrates the complete forward propagation of the network. Such a forward propagation is also part of the inference phase, when a simulation with surrogate elements is performed. The main steps in the linear case are the following:

1. Propagation of the network (export bias for linear models: u(−1) = 0)

2. Scaling and transformation of the networks output

3. Computations of surrogate elements matrices

4. Computation of global matrices

5. Static or dynamic solving

• 1: u_IR−1 ⇐0 #Input for linear modeling

• 2: K_SE ⇐0

• 3: $K̃⇐adapt_dimensions(KPR)$

• 4: $f̃⇐adapt_dimensions(f)$

• 5: k ⇐DENSE(u_IR−1)

• 6: k ⇐tanh(k)

• 7: $k⇐sign(k)efkk$ #Scaling

• 8: a, b = indices_lower_triangular_matrix(N_SE − 3)

• 9: K_SE(a, b) ⇐k

• 10: K_SE(b, a) ⇐k

• 11: $KSE⇐KSEKSET$

• 12: K_SE ⇐R(ψ, θ, φ)K_SE, dofR(ψ,θ,φ)^T #Correct orientation of the SE

• 13: K_{SE, idx_right_column} ⇐sum_comp_row(K_SE)

• 14: K_{SE, idx_bottom_row} ⇐sum_comp_col(K_SE)

• 15: h, p = indices_interface_regions(Γ_IR)

• 16: for i do #Add SE to conventional elements

• 17:  $K̃(h(i),p(i))⇐KSE$

• 18: end for

• 19: $L⇐chol(K̃)$

• 20: $α⇐L−1f̃$

• 21: $ũ⇐(LTα)−1$

• 22: $ũ⇐ũ(≠additional dofs)$ #Remove adofs

The mass matrix of finite element simulations describes the distribution of the mass of the domain Ω to the dofs. As a surrogate element should not change properties of the standard finite element matrices, a surrogate element mass matrix must also possess two essential properties. These are symmetry and strictly positive entries, due to the non-physicality of negative masses. The positive entries of the surrogate element mass matrix are ensured by applying a sigmoid activation function after the dense layer, see Equation (4). This vector is then multiplied by the materials density ρ. Hence, the surrogate elements are not dependent on the density in the training dataset. The problems under investigation are of a three-dimensional structural nature, so that each dof is assigned to a directional component, x, y or z. Unlike stiffness matrices, where x, y and z components are not independent due to the material law, this is the case for mass matrices. Taking advantage of this circumstance, the dimension of m_SE does not have to be $NSE2$⁠, it can be reduced. The values from m_SE are sorted in the following manner:

where m_SE,i is the i-th component of m_SE, with $mSE∈RNkSE/3×1$⁠. A symmetric M_SE is ensured due to the sorting process. The proposed approach is only applicable for N_SE being a multiple of three, which is suitable for the directional components in structural mechanics.

Assembling $K̃$ and $M̃$ is identical and therefore only explained for the mass matrices. First, conventional elements originating from Ω_PR are built and added to $M̃$⁠. Subsequently, M_SE is added to the entries of M_IR in $M̃$ and the optional adofs. A graphical representation of this procedure is provided in Figure 3.

Prior to the training of the surrogate element, training, validation and test data have to be generated in an in silico manner. For this, the full model is solved for several random distributed load cases. Subsequently, the displacement vectors u(t) can be stored as dataset alongside the force vectors f(t). It is important to note that at this stage, a consistent system of units must be selected. Once the training process has commenced, it cannot be altered.

With ANNs it is best to scale the input data to a prescribed range or standard deviation, see Goodfellow et al. (2016). In the context of our linear elastic use case, this is not necessary, as the input of the ANN, see Equations (3) and (4), is a constant zero vector, u(−1) = 0. Nevertheless, we strive for datasets, which lead to a similar size of error in the training process. This is not a trivial matter, since the calculation of the error is not done based on scalable input data of the ANN. In the data generation routine the exciting forces are drawn from an uniform distribution within the range [−1, 1]. The resulting vector u(t) is scaled to a prescribed absolute mean value of the displacements for each calculation in the dataset, so that every dataset has the same mean. Consequently, the force vector is also scaled by the same factor. We also investigated scaling of the maximum displacement to the same value, which results in inferior datasets.

Furthermore, it is important to generate balanced datasets. Thereby, the forces should act equally in all directions across all simulations in the datasets. To achieve this property, each component of a load is mirrored on every plane in the global coordinate system, thereby directing it towards every octant in the 3D domain. Hence, for each load vector, additional seven load vectors are generated.

To illustrate, the data generation process is elucidated in depth with regard to the case study presented in Section 3.1. In general, five types of load cases are under investigation for the beam example, see Figure 6. Noteworthy, no supports or loads are applied in Ω_RR, marked with yellow dashed lines in Figure 6, since they will be replaced with the surrogate element in the training process:

1. Type a: Random support and random force nodes on both halves

2. Type b: Support on both end faces, force on every node in a random layer

3. Type c: Support on both end faces, force on every node to apply a moment in a random layer

4. Type d: Support on one end face, other end face is supported only in the axis of the beam, force on every node in a random layer

5. Type e: Support on only one end face, force on every node in a random layer

In terms of dataset representativeness, we select load cases that represent the range of practical relevance. Additionally, space-filling considerations are not critical for the application of linear modeling, as the model input is always a zero vector. Hence, modeling is done by the bias, which is independent from the input. Due to linearity, each type of load case and its force directions need not be represented at different load levels, as the ANN-derived stiffness matrices remain constant. Consequently, the dataset size can be relatively small.

The main objective function of the generation of K_SE is to minimize the error between the displacements of the full and reduced model. In respect thereof, only the displacements of the dofs associated with K_PR are relevant. Optional dofs from K_AD are not evaluated. In the case of static models the error is calculated using the mean squared error (MSE) function:

where $nΩPR$ is the number of dofs originating from Ω_PR. Furthermore, a pure static dataset is used for training, see Section 2.6. The ANN is trained with the Adam optimizer, described by Kingma and Ba (2014).

For training of transient models two different approaches are available. At first, training can be conducted based on the stiffness matrices generated by the aforementioned routine. Hence, only the ANN of the mass matrices has to be optimized in addition. For this, a transient dataset is used. In the forward propagation of the transient training procedure, Equation (2) is solved by the Newmark beta time integration scheme for all timesteps n_t in the training dataset. Hence, Equation (13) has to be adapted for the error calculation, since in this procedure $ũPR(t)$ and u_PR(t) can be iterated by t = [1, n_t],

The second approach applies when there is no trained ANN for the stiffness properties. Therefore, both the stiffness and mass matrices must be considered simultaneously. Once more, a transient dataset is required. However, it should be noted that due to the more complex solving process of dynamic simulations, the approach where training is split into a static and dynamic part is less expensive.

In this section, examples of the proposed surrogate elements are presented. At first, only the stiffness matrices are evaluated by static simulations. The second part of the study involves the investigation of transient simulations involving the mass and damping matrix.

The first example to examine is that of a beam structure, in which the central part is replaced by a surrogate element, see Figure 1. The beam has a geometric shape of 100 mm ×10 mm ×10 mm. The linear elastic material model is parametrized based on the properties of structural steel, young’s modulus is E = 2.1 × 10⁵ MPa and the Poisson’s ratio is ν = 0.3. The full model comprises 540 dofs, where 24 are removed for the boundary condition. All elements are linear, regularly shaped hexahedrons. Element integration is conducted with two integration points per coordinate. Initially, there are 80 elements. For this example, 108 dofs are removed from the system to be later replaced by the surrogate element. In total this are 20 elements, which are marked with yellow dashed lines in Figure 6.

The data generated for the training process consists of 2,048 samples, distributed equally between type (a)–(e). As the model is small, generation of all data points takes 10 s on a single cpu. A total of 70% of the data is allocated for training, 20% for validation and 10% for testing, with the categories assigned randomly. The training dataset is used for training a model, whereas the validation dataset is used to check for overfitting and hyperparameter optimization. The test dataset is utilized to assess the performance of the fully trained model. The mean value for the displacements of all datasets is chosen to be u_mean = 0.1 mm. The overall training time is 122 min on a Quadro P2200 graphic card. Furthermore, all neural network related operations are implemented by PyTorch (Paszke et al., 2019). The parameters and the additional hyperparameter f_k are shown in Table 1.

Particularly, f_k, see Equation (8), is crucial for determining the appropriate stiffness, as low values result in an excessively flexible stiffness matrix. The learning rate α is not constant over the epochs. The initial learning rate is gradually decreased linearly to enhance the stability of training and improve the ANNs performance. As illustrated in Figure 7, training for 100 epochs is sufficient, since the loss already has converged at this point.

Reference solutions u_i and investigated solutions $ũi$ are compared using a local relative error measure to ensure comparability across different load cases

with i as node number. However, for training a MSE measure is used, since a relative reference is not necessary in the optimization process and MSE is the standard with ANNs in our field. For a global assessment of the results, the mean of the absolute values calculated by Equation (15), ϵ_m.a.r., is used.

A comprehensive evaluation of the test dataset for each load type reveals low ϵ_m.a.r. as shown in Table 2. Due to the random nature of load case (a) and the sufficient accuracy, it is evident that the surrogate elements are capable of generalizing to different types of loads and supports.

In addition to that, two test load cases, one of type (d) and one of type (b), see Figure 8, are investigated in detail. In the first example, all dofs of the left side of the beam are fixed supported, whereas on the right side, only the dofs in the beams axis are fixed clamped. For the second example, all nodes on the left and right side are fixed clamped. In both cases, the structure is excited on every force node by a force orthogonal to the beams axis with f_y = 55.6 N in the first and f_y = 1667 N in the second example.

The proposed method yields promising results, as illustrated in Figure 9. Displacements of both test examples cover the reference solution with a very low relative error. In the case of the fixed-sliding example, the reference solution and the surrogate element solution diverge by ϵ_m.a.r. = 0.51%, while in the case of the fixed-fixed example, the error is even lower, ϵ_m.a.r. = 0.48%.

In summary, the proposed surrogate element is capable of accurately reproducing the results obtained from a FE simulation.

In this section one more crucial aspect should be addressed. In the case of the beam example, which is fixed supported on both ends (load case d), it is basically possible to replace the surrogate element with four long finite elements. Hence, the model is fully modelled with standard finite elements and has the same number of dofs, and therefore, can be solved by identical computational costs as the model with surrogate elements. But, these four elements show a bad aspect ratio of five. In Figure 10 both versions can be seen. For the surrogate element, a maximum relative error of up to 3.5% can be calculated, whereas for the pure FE solution a higher relative error of 15% is present. It is noteworthy, that for the surrogate element solution, the error in the most deformed region is below 2%, with a majority of errors being below 1%. From a global perspective, ϵ_m.a.r. = 7.3% for the pure FE solution and ϵ_m.a.r. = 0.47% for the surrogate element solution. In summary, although from a computational point of view, the surrogate element is not advantageous in this comparison, but it clearly outperforms the coarse mesh in terms of precision.

The surrogate element demonstrates accurate results for models that are, from a geometric view, identical to the training data. The fundamental concept underlying the surrogate elements is the utilisation of such an artificial element in a multitude of distinct models. So it is to determine if the position of the surrogate element can be altered. Furthermore, if the stiffness behaviour of the region that is replaced with the surrogate element is accurately modelled, it should be possible to use the surrogate element on multiple occasions within the same model.

We use a longer beam with a length of l = 225 mm, where the surrogate element is applied five times, see Figure 11. The beam is again fixed at both end faces, a force of f_z = 222.2 N is applied on each force node.

The results obtained with the surrogate and FEM reference demonstrate a high degree of accuracy. Nevertheless, ϵ_m.a.r. = 3.7% is slightly higher than in other examples. Further studies with a training dataset, wherein the surrogate element is trained on examples with more than one surrogate element in a simulation, yielded more favourable outcomes in this instance. However, optimising in this manner would result in surrogate elements that are biased towards the specific training data. Conversely, our objective is to facilitate the flexible application of the surrogate elements. This example allows for the assessment of the benefits of such a surrogate element. While K_SE is used several times, it only has to be propagated once. Due to the reduction from N_{dof, full} = 1,242 to N_dof,SE = 702 dofs, computational time with our in-house FE solver for the long beam example drops from t_ref = 0.038 s to t_SE = 0.024 s. To ensure comparability, computational times are averaged over ten simulations, with both types of simulations executed on a single CPU processor. Thus, a speed-up factor of 1.6 can be achieved in this example. Overall, assessing the training time of 122 min, a dynamic simulation of $>5.2s$ with a timestep of 1 × 10⁻⁵ s amortizes the training costs.

In another, example the adaptability of the surrogate element should be discussed in greater detail. The structure under investigation initially comprises 604 elements and 1,332 nodes, see Figure 12. A force of $222.2,−222.2,0N$ is applied to each top node of the vertical structures. All dofs at the bottom are fixed. Hence, $K∈R3861×3861$ can be stated.

The surrogate element allows for the omission of 360 elements from the full model, which represents 1,944 dofs. This results in a reduced stiffness matrix with dimensions $K̃∈R1917×1917$⁠. Our in-house FE code requires 1.26 s to solve for u, whereas the same solver requires 0.12 s to solve for $ũ$⁠. Hence, a speed-up factor of ten can be achieved. The preserved dofs match accurately to the reference solution. The relative error is greatest at the nodes where the boundary conditions are located. This is due to the fact that small displacements lead to high relative error values, although the absolute error is minimal. Of greater importance are the regions with high displacements, for example at the top of the structure. In this region, the absolute relative error does not exceed ϵ_rel = 2%. When all nodes are considered the error is ϵ_m.a.r. = 2.0%.

One more aspect, which should be examined, is the influence of adjacent elements. The surrogate element is only applicable when Γ_IR and consequently K_IR are identical to those used in the training process. In the following example, the elements at these positions are altered to deviate from the ideal hexahedrons. The circumference of the beam is increased for the adjacent node layers, see Figure 13. The beam is once again fixed supported on both end faces. A force of $0,555.6,555.6N$ acts on each purple marked node.

Once more, the greatest deviations regarding ϵ_rel can be observed in regions with low displacements. Nevertheless, the solution obtained with the surrogate element shows once again high accuracy, as the global error is ϵ_m.a.r. = 0.80%. Because of that, one can attest the surrogate element a high application flexibility. To summarize, surrogate elements can be employed independently of an arbitrary mesh, provided that the interface regions are identical to those used for training.

In this section, the surrogate element is employed multiple times within one model in transient simulations. Hence, besides K_SE, M_SE is also generated by an ANN. First, a beam structure analogous to that presented in the previous section is trained. For this, the data generation routine from Section 2.6 is utilized. Hence, the surrogate element is trained on a single beam with different load cases and boundaries. The length of the partition, which is modelled by the surrogate element, is now eight elements. Consequently, it is longer compared to the one presented with the static simulations where the surrogate elements length is five conventional elements.

Following training, the surrogate element is applied to a scaled schematic model of a multistorey building, which is excited according the FF earthquakes NS acceleration curves, see Krishnamoorthy and Anita (2016). The full model comprises 26.7 × 10³ dofs and 4,068 elements. The surrogate model is built with 99 surrogate elements and 900 conventional elements. Consequently, there are only 8,019 conventional dofs and 99 ⋅ 3 = 297 adofs.

The results presented in Figure 14 show high accordance to the full solution.

Similarly, the evaluation of displacements of the marked evaluation nodes from Figure 14 over time is accurate, see Figure 15. In the majority of the timesteps, the relative error is below 2%. The evaluation nodes in the vicinity of the excitation source exhibit the highest error, with an average of approximately 4%. Due to a slight shift in the oscillating nature of the displacements, on some timesteps the relative error is higher. At t = 0 s and t = 0.02 s, where the displacements are close to zero, the relative error measure provides high discrepancies. However, this is not crucial, as the absolute error is actually quite low.

In this example, the damping parameters are set to α = 1 and β = 1 × 10^−4.  These values are not included in the transient training dataset, as there the damping is α = 2 and β = 1 × 10⁻⁷.  This shows, that the damping specifications are not dependent on the training. From a computational perspective, solving for the full model requires 11.9 × 10³ s, whereas solving the reduced model only takes 370 s, resulting in a speed-up factor of 32.

In a more complex example, a schematic gear is evaluated. But instead of training the surrogate element on the complete gear, it is trained only on a quarter disc. The boundaries and forces are applied according Figure 16. The gear is modelled with hexahedrons, based on 972 dofs in the full approximation and 720 + 3 dofs in the reduced approximation. The material is identical to the beam examples. Two types of simulations were conducted, one in which forces were applied to each face, and another in which all faces were excited simultaneously. In order to create a balanced dataset, all three directions of the forces are alternated, as explained in Section 2.6.

The model is trained in accordance with the split approach outlined in Section 2.7. Hence, K_SE is based on a static dataset, whereas the transient dataset is the foundation of M_SE. Two distinct damping sets are considered, α = [1,2] and β = [10^–5, 10^–6]. After evaluation of the dataset it is evident that the proposed method offers a significant advantage, as solving a full simulation for training takes 1.7 s. In contrast, solving the same model with the surrogate elements only takes 0.9 s. The use of the surrogate model resulted in a reduction in the time required by 47% for the example presented in Figure 16.

After training, the surrogate element derived from the quarter disc is applied to a gear. In contrast to the training model, the gear is now a full disc and on the outer boundary of the gear, teeth are present. Consequently, the excitation is now applied to the teeth, which represents a distinct load case compared to the training model. The gear is fixed at all nodes of the inner circle. Furthermore, the damping values, which were not included in the training dataset, are selected as α = [1] and β = [2 ⋅ 10^–6]. The full model is discretized by 5,328 dofs, while the reduced model contains only 4,320 + 4 ⋅ 3 dofs. Despite the differences between the models during training, the reduced simulation demonstrates high accuracy with the full model. The meshes of the cogwheel at t = 0.006 s, are shown in Figure 17.

The maximum absolute relative deviation of 2.3% is observed at the nodes of the excited teeth. This phenomenon can be attributed to the tilting of this teeth because of some erroneous flexibility of the base cylinder. Apart from this, the deviations are below an absolute relative error of 0.5%.

Furthermore, the deformation of the single nodes over time is consistent with the reference, as illustrated in Figure 18. For this, the nodes marked in green in Figure 17 are evaluated.

A relative mean error of −2% can be identified. Hence, the structure is somewhat over-flexible when utilising surrogate elements. In the event of a significant change of the exciting signal, the error increases in an oscillating manner. Although, the errors are oscillating, the nature of the displacements are met with high accuracy. The discrepancies between the surrogate and reference curves are not transient in nature, rather, they are the result of a minimal deficiency in the stiffness of the surrogate elements. This issue requires further investigation and resolution, potentially through the use of a more suitable dataset and additional training.

The utilisation of surrogate elements enables a 37% reduction of computational time for the full model from 490 s to 310 s of the surrogate model. By increasing the area of the surrogate element, it is possible to achieve a further reduction in computational costs. Once more, all computations were conducted on the same machine with the same in-house FE-code. In a use case, where the inner part of this model remains constant and only the regions near the teeth undergo iterative design changes, the training effort of 36.9 × 10³ s for the K_SE and 7.5 × 10³ s for M_SE becomes more negligible. Furthermore, when transient simulations are conducted over a longer time span, the amortisation of the training process is accelerated. To illustrate, if the presented cogwheel is simulated for approximately at least t_sim = 2.5 s, training and solving of the surrogate model is faster than solving the full model.

This paper addresses the training and application of surrogate elements, which are capable of replacing several finite elements in the sense of a load and boundary independent substructuring for static and transient simulations. Furthermore, the presented surrogate elements can be employed for machine learning-based modeling of macroelements in the field of multiscale simulations. The replacement of standard finite elements with surrogate elements results in a reduction in the dimensions of the discretised system of equations. Consequently, the computational costs for solving them is reduced up to 32 times for one of the presented examples. The stiffness matrices of a surrogate element generated by ANNs exhibit identical properties to those of standard finite elements. As a result, they can be seamlessly integrated with standard finite elements. Therefore, a strong focus is on the creation of artificial stiffness matrices, which are psd. Besides the stiffness matrices of the surrogate elements, also mass matrices for transient simulations are trained. Furthermore, the introduction of artificial degrees of freedom, which cannot be mapped to a specific geometric representation, is proposed. A data-driven in silico approach is employed for the training of the ANN. It is demonstrated that training surrogate elements on models with only a limited number of dofs is sufficient for the desired outcome. Experimental results demonstrate the high accuracy of computations using surrogate elements on complex structures for static and dynamic simulations, even when trained on simple models. Consequently, any structure can be considered, as long as the interface regions connecting the surrogate elements and the standard finite element regions are identical. The main contribution of this work is the presentation of a method that is highly flexible in its applicability. It fills the gap in surrogate modeling on the system level and on the element integration level, as it is with constitutive laws. The surrogate elements can be interpreted as substructures that are not dependent on boundary conditions and load cases in the offline generation process.

In future work, our aim is to incorporate the ability to model different material properties of linear materials. Also, it should be evaluated how the proposed approach scales for even larger problems, for example in simulations with very fine meshes where the surrogate elements replace much bigger regions. The ANNs input could consist of two vectors, one representing displacements and the other capturing material properties, such as the elasticity modulus corresponding to the chosen material model. Regarding accuracy of the method, it could be further enhanced by implementing more precise training procedures and datasets. Furthermore, the influence of artificial degrees of freedom on the accuracy of the surrogate modeling should be assessed in detail. It is also our intention to apply surrogate elements to geometrically nonlinear problems to fully exploit the potential of the presented method with non-constant input vectors. Furthermore, this approach could be extended to enable the replacement of specific regions of FE models that deal with contact.

This open access publication was funded by Ostbayerische Technische Hochschule Regensburg.