Integration of PGD-virtual charts into an engineering design process

This article deals with the efficient construction of approximations of fields and quantities of interest used in geometric optimisation of complex shapes that can be encountered in engineering structures. The strategy, which is developed herein, is based on the construction of virtual charts that allow, once computed offline, to optimise the structure for a negligible online CPU cost. These virtual charts can be used as a powerful numerical decision support tool during the design of industrial structures. They are built using the proper generalized decomposition (PGD) that offers a very convenient framework to solve parametrised problems. In this paper, particular attention has been paid to the integration of the procedure into a genuine engineering design process. In particular, a dedicated methodology is proposed to interface the PGD approach with commercial software.


Introduction
Due the rapid and constant increase in computing power, particularly with the development of High Performance Computing, it is possible to run simulations of complex structures, which was unimaginable a few decades ago. It is, then, possible for some structures to run a parametric simulation (for all the values of parameters) and to choose the right set of parameters for the structure to design afterwards.
Nevertheless the tendency nowadays is to introduce more and more physics into modelling, i.e. dealing with very complex material laws, geometric non-linearities due, especially, to the development of composite materials in many engineering fields (aeronautics, automotive...). Simulations of these models can take a substantial amount of time (weeks, months). Thus considering each new structure, i.e. a new set of parameters, as a new problem during the design stage leads to highly expensive simulations. Consequently, structures tend to be oversized in order to reduce the design time (the design stage is stopped before finding the optimal set of parameters) which is an issue, for instance, in aeronautics where weight reduction is an engineering challenge.
To fix ideas, the optimisation problem would be the minimisation of the mass M of an aeronautical structure with respect to a set of design parameters α = (α 1 ,...,α n ) ∈ A = A 1 ×···×A n ⊂ R n . Meanwhile, the Von Mises stress σ VM must remain under a given threshold σ 0 in order to guaranty the safety of the structure. The constrained optimisation problem can be written as follows: Find α * = α * 1 ,...,α * n ∈ A = A 1 ×···×A n such that which can be rewritten by restraining the space A to A σ 0 ,the space where the parameters associated to structures, whose Von Mises stress do not exceed σ 0 , belong. Hence, the optimisation problem becomes: Find α * = α * 1 ,...,α * n ∈ A σ 0 such that α * = arg min However, the difficulty lies in the definition of the search space A σ 0 because it implies evaluating the stress field σ for all α ∈ A, which can be out-of-reach of classic methods in terms of computation time.
The aim of this work is not to derive any new approach on the optimisation process itself because numerous algorithms are already available but to propose a strategy to build efficiently the fields needed for evaluation (stress, strain, displacement, quantities of interest, etc.) It is worth noticing that a non negligible part of structural design is repetitive since many structures are similar in shape. So it could be helpful to store the results obtained for previous structures to use them for new ones. That is why the idea of handbooks is still of topical interest. Handbooks were, for centuries, powerful decision support tools in many fields such as architecture, engineering (Vademecum des Mechaniker [7] by Christoph Bernoulli) and science in general.
A novel idea is to update the notion of handbooks to current engineering requirements. These new kind of handbooks are called virtual charts [17]o rcomputational vademecum [11]. Varying the values of a set of geometric parameters, it is possible to construct a whole family of structures from a reference one. Results of numerical simulations are stored once and for all in virtual charts for families of structures. Hence the end-user has solely to seek the sets of parameters that meet the specifications among all the sets considered in a certain range. The construction of a virtual chart may be time consuming but, since it is done once and for all, it is worth loosing time at this stage with regards to the one won during the practical use, above all for repetitive tasks.
However, as it was previously mentioned, the construction of virtual charts is still out of reach of standard Finite Elements approaches using "brute force" due to the prohibitive computation time.
The path, which is followed in this work, is to build the charts by mean of a reduced-order modelling technique. The community of model reduction is very active and many strategies are, nowadays, available. One of the most popular approaches relies on the use of the proper orthogonal decomposition (POD), which is based on some preliminary computations, called snapshots, of the high-fidelity problem for given values of the parameters (see e.g. [6,15,19,20,31]). The reduced basis method (RB, see e.g. [27,30]) adds an automatic selection of these snapshots by a greedy algorithm based on some error indicators. Finally, the (PGD, see e.g. [4,16]or [10] for a review of the method) follows a different path as it builds progressively an approximated separated representation of the solution, without assuming any snapshot or basis.
The PGD was considered, in our former works, in the context of the LATIN method [16], in particular for the analysis of elastic-viscoplastic problems [29], multiscale problems [12], multimodel problems [22], multiphysics problems [24] or parametrised studies [23]. The PGD has also been widely developed by Chinesta and co-authors, who proposed very efficient implementations of the method (see [9] for an overview). To focus only on some very recent works, one can cite the examples of real-time simulations in surgery [2,14], real-time monitoring of thermal process [1] or the simulation of viscoelastic models [5].
Concerning the case of geometric variations for shape optimization, the literature in the field of model reduction is poorer as it involves some technical difficulties that will be discussed along this paper. Some POD based approaches deal with this issue such as [21,30], where geometric parameters are taken into account in the RB framework and [28], where the authors handled it through a combination of POD (called Principal Component Analysis here) and Diffusion approximation. The PGD approach is, obviously, a natural framework to deal with geometric variations by handling geometric parameters in the separated variables decomposition and then solve the problem for any set. The first demonstration of this approach has been done in [3] on rather academic geometries of thermal problems and reused in [32] very recently.
The aim of the present work is to address more complex problems encountered in industry such as axisymmetric geometries and shapes defined by non-straight lines such as splines and to take care of the implementation in the engineering process as well as the coupling with engineering software. In our knowledge, the coupling between PGD and commercial software is quite new in the literature due to the high degree of intrusiveness of the PGD method and marrying the new algorithms with the techniques and tools that are used in engineering design offices is a clearly mandatory. This work is a first demonstrator of what can be envisaged in the future to introduce PGD in industry.
The article is organised as follows: in Sect. 2, we present the engineering case ; in Sect. 3, the notion of virtual chart is defined and the Proper Generalized Decomposition is introduced ; Sect. 4 deals with geometric parameters and how the PGD has to be modified to take them into account ; Sect. 5 treats the integration of the PGD into an engineering design process using commercial software ; finally, Sect. 6 exemplified the approach on the demonstrator proposed by AIRBUS Defence & Space. In Appendix, some discussions on the PGD can be found for the novice reader.

Presentation of the engineering structure
The problem that will be studied in detail in Sect. 6 is an engineering case. During the design process, different parameters can be played with such as loads, material parameters (Young's modulus, Poisson's ratio...) and geometric parameters. For our study, we focus on geometric parameters. As shown in Fig. 1, the considered structure is an axisymmetric geometry in isotropic, linear elastic material, whose boundary is defined by splines and parametrised by the positions of control points 1 and 2. In other words, there is a total of four geometric parameters (two per point). The set of parameters is denoted as α = (α 1 ,α 2 ,α 3 ,α 4 ) = (r 1 , z 1 , r 2 , z 2 ), which belongs to design space A. One of the quantities of interest sought for design is the size of the plastic zones. For the sake of demonstration, the computations performed in this article are linear elastic and the plastic zones are estimated by defining a threshold and considering the zones whose strain is superior to this threshold (see Sect. 6.2 for further details). An example of virtual chart representing this quantity of interest will be constructed and exposed in Sect. 6.2.
In the current approach, the problem is solved, after discretisation in space, through a FE procedure governed by an equation of the form: where The brute force approach consists in solving (3) for each set of parameters and checking whether the structure meets the mechanical requirements or not. If not, a new simulation is run and so on. The issue is that numerous simulations may be needed before finding the right set of parameters increasing the computation time and conception costs. That is why reduced order modelling and, in particular, PGD was considered. Difficulties will be encountered due to axisymmetric modelling and the complex geometry considered (splines). A methodology to take care of them will presented.

Virtual chart
Coming back to the continuous framework, a PGD-virtual chart of the displacement u, depending on a set of parameters α ∈ A, can be seen as a separated variables approximation of u: where λ k (α) = m j=1 λ j k α j and n is the number of "modes" used in the approximation. λ j k are scalar functions and k are vectors. All is remaining for the final user is to particularise this solution for a given set of parameters.
For parameters such as material parameters (Young's moduli, Poisson's ratios, diffusion coefficients...), the PGD method furnishes a separated variables representation thanks to a greedy algorithm [18]. For geometric parameters, the method cannot be applied directly as is and must be slightly modified as it was exemplified in [3] (see Subsect. 4).
The construction of the virtual chart, using a greedy algorithm, is classic and recalled in Appendix 1 for the reader who is not familiar with the method. The important point is that, roughly, the greedy algorithm is based upon a variational formulation of the form: Find u ∈ I ⊗ V such that where ⊗ stands for the tensor product. The whole PGD algorithm is exposed in Algorithm 1. An example of application for material parameters is given in Appendix 2.

Geometric parameters
Geometric parameters are of different nature regarding to other parameters as, for instance, material parameters (Young's moduli, Poisson's ratios...) due to the fact that their variations affect directly the integration domain Ω (α).L e t us come back to the problem shown in (5). This time, we will assume that the geometry depends on a set of parameters α = (α 1 ,...,α m ) ∈ A = A 1 ×···×A m . Equation (5) becomes:

Algorithm 1: PGD greedy algorithm
..× A m Level 1: Greedy algorithm; 1 while R (u n ) > threshold (with R (u n ) being the residual) do 2 Level 2: Fixed-point algorithm; 3 Initialisation of the (m + 1)-uplet λ 1 ,...,λ m , ; 4 for k = 1 to k max do 5 for j = 1 to m do 6 solve the parametric problem associated to λ j 7 solve the spatial problem associated to 8 Update of the functions This complicates the calculations of the different integrals since they cannot be computed independently from one another any more. The dependency of the integration domain on the geometric parameters has to be, in some way, skirted. The solution, introduced in [3], is to use a change of variables, i.e. to define a reference structure Ω 0 and, from the current structure Ω (α), where the problem is initially defined, to come back to that reference structure Ω 0 thanks to a geometric transformation T α (see Fig. 2). Thus, The methodology is recalled, hereafter, for the sake of consistency and to highlight the technical difficulties of dealing with axisymmetric modelling. One can write (6) applying the change of variables: In this formulation, the dependency on the geometric parameters has been transferred from the integration domain to the integrand through the Jacobian. Now a PGD approximation of the field u can be sought using the PGD procedure already exposed in Sect. 3.

Proper generalized decomposition for geometric parameters
Let us come back to the problem defined in Sect. 2.O n e applies the principle of stationary potential energy and integrates it on the geometric parameter space, which reads after discretisation in space: where N is the number of degrees of freedom of the discretisation in space. From this formulation, the PGD-approximation of displacement vector U is calculated: At iteration n, the following test function vector is used for the computation of the new (m + 1)-uplet: withŨ * ∈ R N and λ i * ∈ I i for i = 1..m. Using Algorithm 1 defined previously, one solves alternatively the linear systems: Let us remark that the computation of the various integrations can be greatly facilitated if the stiffness matrix K is written in a separated variables form, which is the aim of the next subsection.

Stiffness matrix computation in the axisymmetric case
Let us start writing the expression of the stiffness matrix K e (α) of an element for the current structure. Applying a change of variables, it will be defined on the reference structure: where -B: matrix of the partial derivatives of the shape functions r : radius -C: Hooke's matrix (here independent on α because we are not dealing with material parameters) -B (α, X 0 ) = B (T α (X 0 )) andr (α, X 0 ) = r (T α (X 0 )) It can be immediately noticed that axisymmetric modelling adds an extra difficulty with respect to a classic two dimensional modelling sincer , which depends on the geometric parameters, appears in the integrand. To obtain a separated variables form of K e (α), one must first writeB in the separated way. Let us write the matrixB: where p corresponds to the number of shape functions used and: It has to be recalled that the problem is defined on the reference structure Ω 0 . Therefore N I,r and N I,z are unknown on the contrary to N I,r 0 and N I,z 0 . Obviously, it is possible to express N I,r and N I,z in function of N I,r 0 and N I,z 0 : Hence, Introducing the adjugate matrix of J T α : (18) becomes, defining Adj = adj J T α T : ∀ I ∈{1 .. p}, Jacobian matrix J T α of the geometric transformation depends on both space variables and geometric parameters. Therefore obtaining a separated variable representation is not trivial. It can be considered to obtain a separated form of the Jacobian through a higher-order singular value decomposition (HOSVD [13]) but due to the large number of HOSVD needed (four for each elements), it would highly increase the computation time. The proposed solution is to avoid the dependency on space variables of the Jacobian matrix. To do so, the simplest way is to use solely affine geometric transformations with respect to space variables. During the computation of the Jacobian matrix, the space variables will disappear naturally. For that purpose, the structure is divided into sub-domains in which is defined an affine geometric transformation. This solution was also proposed in [3] where geometries are decomposed in triangular sub-domains.
In Sect. 6, a possible subdivision is detailed for the engineering case. In particular, it will be shown that triangles can be limited in some cases of complex geometries such as splines and hence are not a panacea. An affine geometric transformation is written in the section as: with X 0 = (r 0 , z 0 ). Hence the Jacobian and Adj matrices become: The Jacobian does not depend on the space variables any more but it is still not obvious to give a separated variables representation of matrixB. For this purpose, let us define the following matrices: and coefficients: So thatB can be conveniently written: Due to the definition of B 5 , one can remark: which leads to the following expression of K e : This expression of K e satisfies only partially the desired objective, i.e.obtaining a separated variables formulation. Partially becauser and 1 /r contain both space variables and geometric parameters. Since the choice of an affine geometric transformation was made,r is already separated and is written asr = a r r 0 + b r z 0 + c r .
The case of 1 /r is slightly trickier since 1 /r has no immediate separated form. For this reason, a linear approximation of 1 /r is computed in a neighbourhood of the barycentre (r bar , z bar ) of the sub-domain to which the element belongs: In Sect. 6, an evaluation of the error for this approximation will be given in the case of the industrial demonstrator. Finally the elementary stiffness matrix K e can be written in a fully separated formulation: The different matrices that appear in the decomposition of K e (α) are to be computed for each element and the coefficients depending on the geometric parameters for each sub-domain. The number of terms belonging to decomposition of the stiffness matrix is fixed with respect to the discretisation of the parameter space A.

Integration into the engineering design process
As it was previously said, particular attention has been paid to inserting the methodology developed into a genuine engineering design process. Due to engineering requirements and, in particular, the coupling with other parts of the system that affect the loadings applied on the structure presented Fig. 1, the actual design process is performed using the commercial software SAMCEF.
It is important to note that a single optimisation of the structure with respect to geometric parameters should have certainly been done efficiently using the embedded SAM-CEF optimisation procedures. However, this work is the first step of a more complex design process in which some other parameters will be taken into account. Further works are in progress to take also into account various types of loadings that can be applied on the structure and, then, the total number of situations considered justify the offline building of a virtual chart of solutions.
The proposed PGD approach must, then, be included in the design workflow with a limited number of modifications to take advantage of the existing facilities. The main problem of the approach detailed in Subsect. 4.2 is the high degree of intrusiveness of the PGD. Indeed, the linear system defined in (14) is not a classic finite elements problem since the stiffness matrix and the load vector are averaged on the parameter space A. To limit the intrusiveness in the design workflow, different approaches have been envisaged, in collaboration with SAMTECH, the SIEMENS subsidiary which develops SAMCEF software. The first was a complete embedding of the PGD approach in SAMCEF, calculating the averaged stiffness matrices and averaged load vectors (see Eq. (14)): directly in the commercial code. This approach has been given up because it would have need the development of a new structure of fields in SAMCEF which has been considered a too consequent job. Finally, the technique which has been chosen and is demonstrated herein is to "encapsulate" SAMCEF into an in-house MATLAB code that deals with the particular operators associated to PGD. To sum up, on the one hand, SAMCEFs utilities are used to defined the geometry, the partition, the loads, the mesh and the materials, as well as to solve the space finite elements systems mandatory to build the space modes i of the PGD. On the other hand, the MATLAB in-house code is used to build the specific PGD operators Eqs. (35) and (36), as well as the problems corresponding to parametric modes λ i . The communication between the two codes has needed the development by SAMTECH of specific interface functions which takes advantage of the possibility offered by SAMCEF to stop it at a defined step and to restart it from this step afterwards. The whole procedure is detailed in Algorithm 2 and a scheme of the communication between SAMCEF and MATLAB is given in Fig. 3. The overall process of coupling is not optimised and this work only constitutes a demonstrator of the approach that will be extended in the future.

12
Transfer stiffness matrices and load vectors to SAMCEF

Geometry partition
In order to apply what has been exposed in Sect. 3, a division of the geometry of the engineering structure has to be chosen. The division has to verify the sufficient condition defined before, viz. that the geometric transformation associated to each sub-domain must be affine with respect to the space variables. An example of a possible division of the geometry is represented in Fig. 4. The level of error due to the approximation of 1 r Eq. (33) goes from 1.28 % (sub-domain 1) up to 6.32 % (sub-domain 2). The sub-domains are mainly triangles where the geometric transformations are quite immediate to compute. To construct an affine geometric transformation in the neighbourhood of the fillet is trickier. In fact, the geometry must be accurately fitted, so triangles are to be avoided. One can think about using a high number of triangles to fit the fillet but here we face a problem of computation time because the more sub-domains there are, the longer it takes to compute the stiffness matrix (see Subsect. 4.2).
The problem is that six-nodes triangles cannot be considered either, even though they will accurately fit the curved geometry because the geometric transformation associating a six-nodes triangle to another one is quadratic and not affine.
That is why it was decided to use modified three-nodes triangles with two straight edges and a curved one as shown in the close-up of Fig. 4. To do so, the fillet in the reference structure Ω 0 is modelled as an arc of a circle. Let us suppose that for each geometric configuration, the positions of points A and B (used here as control points), together with the tangents at these points, are known and point C is "free". First let us seek the inverse geometric transformation T −1 α (from the current structure Ω (α) to the reference structure Ω 0 ) associated to the displacements of points A and B which is also affine.  The affine inverse geometric transformation can be written as follows: It remains to determine the six coefficients λ r , γ r , δ r , λ z , γ z and δ z . Four equations are given by (37) for the positions of A and B are known for both current and reference structures.
As it was previously mentioned, the fillet boundary is modelled by an arc of a circle further referred as (Ŵ 0 ) in the reference structure. Thus, Substituting (37)in(39), it leads to: Let us derive it with respect to r : (λ r r + γ r z + δ r − r centre ) λ r + γ r dz dr Which gives the last two missing equations: The six coefficients λ r , γ r , δ r , λ z , γ z and δ z of the inverse geometric transformation being determined, the geo-

Results
Let us come back to the engineering case presented in Sect. 2. Figure 5 shows the variations of geometry, we are dealing with. The convergence curve of the simulation can be seen in Fig. 6. It is possible to improve the convergence rate of the method following strategies presented in [26], what will be done in future developments. Anyway a PGD approximation of the displacement was obtained with 17 modes with an error on the residual of 5.6 % and from it, it was also deduced a PGD approximation of the strain ε. Comparisons between the PGD solution and FE simulations show a good accuracy of the method developed (Fig. 7).
The question now is: how to explore efficiently the virtual charts? Indeed, this is a real difficulty because one must still have in mind that the final user need a decision support tool  as easy to use as a handbook. He must be able to look over all the possible geometries and solutions in an easy and practical way. For that purpose, we exploit the plug-in PXDMF for PARAVIEW developed at Ecole Centrale de Nantes which allows to plot separated variables data [8]. Given the modes associated to each considered field (displacement, strain...), PARAVIEW reconstructs the full fields in real-time by doing summations and multiplications. It allows the visualisation of the evolution of the fields varying each geometric blueparameter.
The criteria used for the design is the size of the plastic zones. For the sake of demonstration in this paper, the calculations performed were linear elastic and the plastic zones are roughly evaluated by defining a threshold and considering the elements whose strain is superior to this value. However the proposed strategy could be extended to truly non-linear computation using the framework proposed in [23]. Figure  8 shows an example of the exploration of the virtual chart using the PXDMF plugin. It is important to note that, the virtual chart being computed, this exploration can be done in real time by the user.

Conclusions
In this work, geometric parameters were introduced into the Proper Generalized Decomposition following [3] and extended to complex geometries defined by splines and axisymmetric modelling. The approach allows to efficiently build a virtual chart of the solution that can be used, afterwards, in an optimisation process. Particular attention has been paid to integrating the methodology developed into a engineering design process. Due to the high degree of intrusiveness of the PGD, SAMCEF was interfaced with an in-house MATLAB code. The specific PGD operators are now computed by the MATLAB code and the resolution is done by SAMCEF after receiving these operators. We recall that a single optimisation of the structure with respect to geometric parameters should have certainly been done efficiently using the embedded SAMCEF optimisation procedures. However, this work is the first step of a more complex design process in which some other parameters, such as loadings, will be taken into account, which justifies the offline building of a virtual chart of solutions. This target demonstrator will be composed of parts, whose behaviours can be highly non-linear (large strains, large displacements, non-linear materials) leading to many configurations of loadings that will be taken into account.
at one side and subjected to an upward load of 1000 N at the other side. The structure is partitioned into three different sub-domains and the Young's modulus of each domain, being considered as a parameter, can vary from 50 to 1000 GPa. The Poisson's ratio is homogeneous and equal to 0.3. Fig. 9 The test structure is partitioned into three sub-domains, the Young's modulus of each domain being considered as a parameter A PGD approximation of the generalised displacement field u is sought: which gives after space discretisation: The modes, i.e. the functions λ 1 i , λ 2 i , λ 3 i andŨ i are computed thanks to Algorithm 1. The four first modes associated to the Young's moduli E 1 , E 2 and E 3 are plotted in Fig. 10 and the four first spatial modes are represented in Fig. 11.    (c) (d) Fig. 11 First spatial modes. a First spatial modeŨ 1 , b second spatial modeŨ 2 , c third spatial modeŨ 3 , d fourth spatial modeŨ 4