Industrial Digital Twins based on the non-linear LATIN-PGD

Digital Twins, which tend to intervene over the entire life cycle of products from early design phase to predictive maintenance through optimization processes, are increasingly emerging as an essential component in the future of industries. To reduce the computational time reduced-order modeling (ROM) methods can be useful. However, the spread of ROM methods at an industrial level is currently hampered by the difficulty of introducing them into commercial finite element software, due to the strong intrusiveness of the associated algorithms, preventing from getting robust and reliable tools all integrated in a certified product. This work tries to circumvent this issue by introducing a weakly-invasive reformulation of the LATIN-PGD method which is intended to be directly embedded into Simcenter SamcefTM\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^{\hbox {TM}}$$\end{document} finite element software. The originality of this approach lies in the remarkably general way of doing, allowing PGD method to deal with not only a particular application but with all facilities already included in such softwares—any non-linearities, any element types, any boundary conditions...—and thus providing a new high-performance all-inclusive non-linear solver.


Introduction
The digitization of the manufacturing industry is a topic of global interest known as "Industry 4.0". Companies are prone to generate a digital representation of their entire value chain in order to be more efficient with a more rational use of resources and energy. One talks about "Digital Twins" [1,2] which include two main components: ultra-high fidelity simulations and data collected during the life of the structure [3]. Our approach is more particularly focused on the first component-namely performing ultra-high fidelity simulations very quickly-and hence we do not intend to deal with data in this paper. Despite the current strong interest in data, one of the requirements for the Digital Twin edifice is still based on the simulation part. Reduced-order modeling (ROM) methods [4] are a key technology for the simulation component requiring a very fast and robust 3D simulation [5]. They have the advantage of drastically reducing CPU time without requiring any compromise on precision where high-fidelity solutions can always be obtained. The richness of the physical phenom-ena remain purely intact, only the way we are representing them-as separate variables format-differs and allows for interesting time savings. Additionally, ROM methods can lead to a better understanding of certain physical phenomena which were previously difficult or even impossible to simulate due to the curse of dimensionality. To be used on a large scale and in a reliable way, ROM methods must be integrated as close as possible into the tools that engineers use on a daily basis, which require the joint use with industrial software. Hence, one requirement is to develop non-intrusive ROM methodologies [6]. Some recent papers are interested in this issue [7][8][9][10][11][12]. Most of the time, it results from a coupling between industrial software used as a black-box and in-house software that drives the ROM methods for parametrized problems on specific applications. Another point lies in the management of time-dependent non-linear problems with ROM methods. In this context, the LATIN-PGD method [13,14] is a worthwhile one, which has been the subject of many works highlighting its assets: ostensible reduction of computation time [15,16], construction of virtual charts [7], ability to address all types of non-linearities [17,18] among others. Hence, the LATIN-PGD method is attractive to Siemens Digital Industries Software, as a promising tool in the Digital Twins landscape. As a software company, Siemens Digital Industries Software wants to integrate this method into the finite element industrial software Simcenter Samcef TM part of the Simcenter 3D TM package. The idea is to get the easiest possible implementation, even if it results in slightly intrusive parts at some places when necessary. This is one of the benefits of having access to the core routines of the software. Three main elements drive the approach: the robustness, the performance and the ease of use by external customers. This work proposes a strategy to use the Proper Generalized Decomposition (PGD) method [19] directly within an industrial software capable of handling any time-dependent non-linear problem in solid mechanics. We develop a weakly-invasive version of the LATIN-PGD method, which differs from the one most commonly used in recent years based on the internal variables description of materials [13,15,20]. Instead, it can be seen as an extension of the first functional version [21][22][23][24]. Numerous numerical treatments have been made to ensure a maximum generality of the presented approach designed to deal with any non-linearity (material, contact, large deformation,...). In addition to performance aspect, robustness remains an essential requirement to be able to deal with any industrial case. While until now the LATIN-PGD method has been used to solve specific applications (see appendix in [25] for related topics), the presented approach shows the possibilities offered by the method with an industrial software where no particular application is targeted. The novelty is the methodology put in place to use the PGD natively and as seamlessly as possible in an industrial software, thus providing an end-to-end process in one single software environment. The LATIN-PGD method treats the problem simultaneously in space and time. Regarding parameters of the Digital Twin, they are managed externally via a classical non-intrusive snapshots approach leveraging that the LATIN-PGD algorithm can be initialized from any previous computation [15,26]. In general, the (extra-)parameters cannot be known in advance for all possible particular applications as they vary from one to another. Hence, by robustness concerns, only a space-time PGD is managed by the Simcenter Samcef TM software itself because space and time are always present in non-linear computations. Then, if required, a space-time-parameters PGD can be obtained a posteriori following classical non-intrusive technics [8,10] based on snapshots. Simcenter Samcef TM is organized as different modules to be able to simulate different physics: non-linear mechanical problems, high-speed rotating machines, articulated mechanisms or thermal simulations. As these modules are all part of the same family, it is possible to switch from one analysis to another, combine analyses or conduct cosimulations. In this work, we focus on the module dedicated to non-linear mechanical solver as an example. However, because of its generality, our approach could be used in other modules (thermal one for example) or other finite element software. The only requirement to set up the LATIN-PGD method is to have a non-linear resolution scheme already in place (Newton-Raphson), which notably solves local non-linearities. This is the strength of the proposed version to reuse without any modification the full management of non-linear equations in the local stage. This paper focuses on the implementation-and the difficulties associated with-of the LATIN-PGD method in an industrial finite element software. Fundamentals as convergence and robustness studies will be presented in companion paper [27]. At the present time, what has been done fits into the framework of moderate displacements, outside instability zones. In what follows, after a brief reminder of the type of mechanical problems that are considered, the LATIN-PGD strategy and some details on its implementation in Simcenter Samcef TM software are provided. Finally, results on two industrial test-cases are given to highlight the assets.

Weakly-invasive LATIN-PGD within Simcenter Samcef TM software
A wide variety of problems in many fields of physics-mechanics, thermics, rotor dynamics-can be handled with the Simcenter Samcef TM software. In this study, we focus on the mechanical solver capable of solving any non-linear evolutionary problem including potentially materials non-linearities, contacts or large deformations. One of the goals is to combine ROM methods (by PGD in this work) with all possible non-linearities that already exist in the software.

Types of mechanical problems involved
Following a generic manner, any physical system subjected to a whole set of external sollicitations (displacements, efforts, pressures, temperature...) over a time interval I = [0, T ] is considered under quasi-static assumption. Using a classical finite element formulation, the problem to be solved can be systematically reduced to the following semi-discretized problem: where Q defines the space of kinematically admissible fields subject to sufficient regularity, q groups all the displacements unknown, F d gathers the external loadings, F i represents the generalized force associated with internal cohesion efforts and T symbolizes a characteristic operator describing all the material non-linearities. This operator T is built from the processing of local elements within the considered industrial software. It will act as a black-box in our approach as explained later.
Remark Under the notation T , internal variables can be hidden if some visco-plastic behavioral laws are considered. Let us note in particular that there is an irreversibility in time within the various non-linearities characterized by the operator T i.e. the expression at a given time t depends on the loading history according to a functional formalism.
In case of contacts, additional contact forces F c have to be taken into account, which requires the introduction of two additional local variables (v, λ) with v the trace of the displacements on the parts in contact and λ the Lagrange multipliers. Hence, the reference problem becomes: where in the finite element formulation, one can define an operator B such as v = Bq and the operator R represents the non-linear contact equations. This problem is classically solved incrementally in industrial software including Simcenter Samcef TM software. For quasi-static sollicitations, a time implicit scheme is commonly used which consists in solving for each time step a non-linear problem that can be treated by a Newton or a modified Newton-Raphson algorithm [28][29][30].

LATIN-PGD methodology
The novelty of our approach is to natively set up a space-time PGD method directly in an industrial finite element software (Simcenter Samcef TM here) in order to best integrate the reduced-order modeling tools into a unified workflow that can be easily handled by customers. One of the requirements is to integrate a new non-linear solver in the software which has to be non-incremental. We choose the non-linear LATIN solver [13,31] designed to work well with the PGD method. This is a non-incremental iterative solver producing at each iteration a better approximation of the solution over the whole space and time. While standard industrial non-linear solvers generally used are incremental by nature, we propose in this work a gateway to effortlessly implement a non-incremental solver directly within incremental one. Schematically, this non-incremental iterative algorithm can be represented as shown in Fig. 1 where two sub-spaces A d and are defined. Sub-space A d groups the linear equations of the problem-equilibrium and kinematic compatibility equations-while subspace groups the local non-linear equations-behavioral laws, contact equations... The solution of the problem (1-2) then results from the fulfillment of this iterative scheme with two search directions ϒ + and ϒ − representing the parameters of the method. More specifically, Eq. (1) can be identified as A d while Eq. (2) defines sub-space. The unknowns in this iterative scheme are the generalized displacements q and forces F i which differ from classical variables most recently used in the LATIN framework based on internal variables. This choice is guided by the finite element software. Usually starting from a linear elastic initialization q 0 (t), the algorithm leads us through a succession of local and global stages by means of the search directions to the terms of a convergent series q k (t) k∈N . We suppose to be at a given state of the algorithm where q n (t) is considered to be known.

Local stage: s n (t) −→ s n (t)
Given s n (t), a new solution s n (t) is sought in sub-space by solving the non-linear equations. One has: Constitutive relations We use a vertical search direction because this is what is used naturally in Newton-Raphson's algorithm. Indeed, at each iteration, one provides the displacement and one recovers the associated internal and external forces to be able to solve the equilibrium afterwards. Kinematics being imposed, this step is local in space and operates at each gauss point of each element. By making this choice for the ascent search direction, we can then reuse exactly the same routines as classically used. The only difference lies in the fact that these routines are no longer used once but several times in order to process all the time steps at each LATIN-PGD iteration.

Global stage: s n (t) −→ s n+1 (t)
Given s n (t) and an operator H that defines ϒ − as a parameter, a new solution s n+1 (t) is sought in A d sub-space by solving the equilibrium equations. One has: Kinematic conditions and equilibrium equation A d q n+1 (t) ∈ Q ∀t ∈ I (10) The choice of search direction operator H is the major parameter of the method. The optimal choice in terms of convergence rate would be to consider the tangent operator at each time steps, while a less costly choice could be to consider only the elastic Hooke operator.
In similar fashion, when contacts are involved, one considers s n = q n , F i,n , v n , λ n which leads to the following LATIN scheme: Local stage: s n (t) −→ s n (t)

Kinematic conditions and equilibrium equation
Both search directions H and h are the parameters of the method. One can interpret h as an interface rigidity [18]. The computation of the R operator comes directly from the hosting software as a black-box. No particular modification is required for the local stage. In this form, the LATIN-PGD algorithm can be seen as the analogue of the Newton-Raphson algorithm to which convergent loop and time steps loop have been swapped. This is a crucial aspect in terms of ease of implementation as the method could be seamlessly integrated into Simcenter Samcef TM software without breaking anything and without having to redevelop specific tools at this level. Thus, it is henceforth possible to choose the non-linear solver that will be used for the resolution: either Newton-Raphson or LATIN methods. In the latter case, as only linear equations on the whole space and time are involved in the global stage, one can use PGD technics [4,32].

Focus on the PGD
The PGD method [31,33] consists in approximating the correction calculated in the linear global stage by a finite sum of functions under separate variables: This low-rank approximation is the reduced-order basis of the problem (1-2), which is obtained using a greedy algorithm. More specifically, a first step consists in carrying out an update phase of the temporal modes already constituted, the spatial ones being kept orthogonal to each other according to the Hooke operator. Then, a second step consists in enhancing the reduced-order basis by adding some new modes if the update phase were proved to be not sufficiently effective. Thus, as the LATIN-PGD iterations progress, one comes to improve the solution until convergence. A criterion based on the energy norm averaged over the time interval enables to determine the state of convergence of the algorithm at each iteration and thus set a stop criterion when a certain threshold is reached. Some actions are currently underway to bring out these modes directly within the postprocessing graphical interface provided by Simcenter 3D TM platform in order to visualized them as easily as possible. The interest is twofold insofar as this data format under separated variables also allows to significantly compress the information on the outputs and therefore increases even more the execution speed of such post-processing interface.

Integration of the LATIN-PGD method into Simcenter Samcef TM software
From a software developer's point of view, only a few parsimonious modifications have been necessary to introduced the LATIN-PGD method in pre-existing Simcenter Samcef TM software notably through the use of a dedicated LATIN-PGD fortran module gathering all the new functionalities in an unique place. Everything that makes up the signature and the richness of Simcenter Samcef TM software-the way of treating behavioral laws, the way of managing contacts and so on-is fully recycled. In a sense, one of the strength of our approach lies in the privileged place granted to this local stage operating at the level of each element, each gauss point of each time steps, thus guaranteeing robustness and generality. A delicate point that has required particular attention concerns data storage since complete spatial-temporal approximations have to be stored at each LATIN-PGD iteration which is quite uncommon within incremental solvers. Nevertheless, the judicious use of pre-existing memory spaces, as well as the use of a few dynamic allocations ensure that the methodology can be carried out smoothly, taking advantage in particular of the fact that the solution can be written in a separable variables format. The result is a controlled memory increase. Seen in a macroscopic way, all this boils down to an interface management between different mathematical spaces grouping together the equations to solve, therefore to establish an efficient communication between the various parts of the code. Another interesting thing is the notion of having a multi-fidelity solver [34]: typically in design offices, one often only needs to know an idea of the solution, a trend, especially in the upstream phases of the design. This allows a better understanding of the overall functioning of a new system before establishing afterwards which studies should then be carried out in greater depth by calculating high-fidelity solutions, always at low cost.

Management of parametrized aspects
In an industrial environment, it is not uncommon to perform these non-linear spacetime computations several times depending on some sets of parameters. The notion of parameter should be taken in a broad sense. By parameter, we mean for instance material data, boundary conditions, number of cycles, geometry variations (morphing),... In what follows, these parameters are managed in a completely external way taking into account an interesting property of the LATIN-PGD method: the algorithm can be initialized with any previous kinematically admissible field. Hence, we do not explicitly refer to the parameters in the notations although our methodology is designed to deal with models involving parameters as illustrated in the numerical results. There is nothing new here, we use the classical techniques as stated in [15,26]. More precisely, if one notes μ the set of p parameters μ = (μ 1 , ..., μ p ), then each LATIN-PGD computation with Simcenter Samcef TM software provides a space-time PGD approximation q (μ) for one value of the parameters such as: Next, if desired, one can then compress the information from these different reducedorder basis to reconstitute a posteriori a space-time-parameters PGD [8,10]:

Some industrial illustrations
In the following, two industrial test-cases are presented. They both highlight, on the one hand, the possibilities of the LATIN-PGD method incorporating within Simcenter Samcef TM finite element software to build efficiently reduced-order basis and, on the other hand, the gain in computation time compared to classical approach.

Combustion engine piston
The first test-case deals with a combustion engine piston submitted to pressure forces on the upper part. Due to the symmetry of the problem, both in terms of loadings and geometry, only a quarter of the piston is modeled and the boundary conditions are set accordingly. In addition, since the piston can only move vertically in its chamber, zero normal boundary conditions are fixed on the outer cylindrical part. Finally, one leads to the resolution of an evolutionary mechanical problem involving 1 million of degrees of freedom, subjected to cyclic pressure loading in order to simulate the back and forth motion of the piston. The material is non-linear Norton law whose equivalent cumulative plastic strain rateṗ follows the formula: where J 2 (σ) represents the equivalent Von Mises stress, k enables to take into account an isotropic hardening and (K, n, m) are three parameters depending on the temperature T of the combustion chamber which is assumed to be uniform for a sake of simplicity. We propose in the following to carry out a parametric study in order to determine the whole space-time solution according to different values of the temperature. In particular, we take out the curve showing the evolution of the cumulative plastic strain during cyclic  Here, a rather naive path have been pursued consisting in starting from the closest neighbor in the parametric space which proves to be very effective when the gradient evolution remains relatively low (typically the range from 515 to 525 • C) and requires a little more time otherwise. In [26] one can find strategies to follow for a more judicious selection of the path. However, we note that we always outperform the classical resolution except for the first calculation which required 50 min-the time to build a first reduced-order basis through PGD from an elastic-solution. Of course, one could probably gain even more by implementing a first phase of estimation or extrapolation of the reduced-order basis as an initialization of the algorithm for a new calculation. The better the estimation, the shorter the computation time, but the fact remains that a reliable high-fidelity solution will always be obtained by proceeding again a few iterations with the LATIN-PGD algorithm if necessary.

Turbine blade-disk assembly
The second test-case concerns a turbine blade-disk assembly involving material nonlinearities-plasticity and creep-and contact non-linearities at the junction between the disk and the blade. The assembly is subjected to centrifugal forces as well as pressure ones on the surface of the blade. Additionally, a temperature gradient is applied according to the radial coordinate. Due to the rotating machine loads and geometry, cyclic symmetry conditions are also prescribed on the lateral parts. This geometry consists of 1 million of degrees of freedom, including 40,000 Lagrange multipliers, resulting from different types of volumic elements-tetrahedrons or hexahedrons-and specific elements at the interfaces between the different parts of the mesh to write kinematic conditions. Because In what follows, a two parameters study is undertaken where the parameters are the Young modulus E and the thermal expansion coefficient α of the blade. Values of these parameters are assumed to be insufficiently known, mainly related to the manufacturing process. Due to loading, a small variation in the value of these coefficients can affect the maximum Von Mises stress in the contact area where creep occurs over long periods of time, which can damage the structure faster than expected. Although other parameters may be of influence, we focus on these two parameters for a sake of simplicity. Figure 6 presents the stress map obtained for a given point (E = 220 MPa and α = 14 · 10 −6 K −1 ) in the parametric space. One gives the results of the LATIN-PGD method included in Simcenter Samcef TM software and the remaining error with respect to the reference. The latter is taken as the classical solution from the commercial software (using Newton-Raphson method). One can note a very good concordance of the solutions as the relative error remains less than one percent in the area of interest for 10 PGD modes. Figure 7 shows the response surface of the maximum Von Mises stress in the contact area according to the two parameters. This chart highlights a non-linear variation according to the Young modulus E and an almost linear one on the second parameter α. In order to achieve this, 70 computations have been carried out by scanning the parametric space sequentially in a uniform manner. Thanks to both ROM methods and LATIN-PGD algorithm allowing to restart the computation from a close previous converged point in the parametric space, only 3 h have been necessary to build the virtual chart while more than 21 h would have been necessary using classical Simcenter Samcef TM software based solely on Newton-Raphson algorithm. This leads to 7 times faster computation. This test-case highlights one of the advantages having the LATIN-PGD method incorporated natively within an industrial software insofar as one can then exploit all the capabilities of the hosting software at the local stage. As an example, 3 different behavior laws (depending on the temperature) have been considered here, as well as different types of elements, in addition to the contact behavior.

Conclusions
In this work, the LATIN-PGD method have been directly implemented into Simcenter Samcef TM finite element software, which allows to build and then enrich reduced-order model as desired at the industrial level. Here, a very general way of looking is adopted with PGD by making the best use of all the local behavior pre-existing in industrial softwares. The PGD method becomes in fact integrated as close as possible to the solver, enabling to combine performance, robustness and a wide range of applications that can be carried out. Moreover, a first time saving on parametric studies has been demonstrated but a significantly additional gain could still be made by exploiting hyper-reduction techniques [35][36][37] and in particular the Reference Point Method [14,38] well suited for LATIN-PGD method. For each quantity of interest, the latter consists in storing the information in an even more condensed manner at a few cleverly chosen points, an analytical formula being subsequently used to reconstruct the field over the whole space and time. These methods are currently being investigated and should enable to go further in the reduction of computation time insofar as up to two thirds of the computation time can be spent in the local stage reflecting the fact that the current algorithm complexity is not independent of the number of degrees of freedom. Furthermore, one also keep in mind that the gain is directly correlated to the complexity of the problem: thus, considering a dozen of parameters and several hundreds or thousands of loading cycles, the performance would be even more better. Examples on material and contact non-linearities have been presented above, but this very general methodology can deal with any other type of non-linearity already present in the industrial Simcenter Samcef TM software. This is the strength of this approach to be able to treat any type of finite element, any geometry, any boundary condition, any loading or even any non-linearities all using reduced-order modeling methods in order to circumvent prohibitive computation times. However, at the present time, computations are limited to moderate displacements, outside instability zones. One of our next task will be to overcome this limitation.