Nonintrusive parametric NVH study of a vehicle body structure

Abstract A reduced order model technique is presented to perform the parametric Noise, Vibration and Harshness (NVH) study of a vehicle body-in-white (BIW) structure characterized by material and shape design variables. The ultimate goal is to develop a methodology which allows to efficiently explore the variation in the design space of the BIW static and dynamic global stiffnesses, such that the NVH performance can be evaluated already in the preliminary phase of the development process. The proposed technique is based on the proper generalized decomposition (PGD) method. The obtained PGD solution presents an explicit dependency on the introduced design variables, which allows to obtain solutions in 0.1 milliseconds and therefore opens the door to fast optimization studies and real-time visualizations of the results in a pre-defined range of parameters. The method is nonintrusive, such that an interaction with commercial software is possible. A parametrized finite element (FE) model of the BIW is built by means of the ANSA CAE preprocessor software, which allows to account for material and geometric parameters. A comparison between the parametric NVH solutions and the full-order FE simulations is performed using the MSC-Nastran software, to validate the accuracy of the proposed method. In addition, an optimization study is presented to find the optimal materials and shape properties with respect to the NVH performance. Finally, in order to support the designers in the decision-making process, a graphical interface app is developed which allows to visualize in real-time how changes in the design variables affect pre-defined quantities of interest.


Introduction
The noise, vibration and harshness (NVH) performance of a vehicle has an ever increasingly stronger impact on the customer perception of ride comfort and brand quality, and it has become one of the most prioritized attributes when purchasing a new car. To be competitive in the global market, car manufacturers need to improve the NVH properties of their products without deteriorating other targets, often conflicting, such as crashworthiness, light-weight, safety, ecological impact and styling.
The automotive NVH optimization represents a multidisciplinary field of research in continuous evolution. To give some examples, it needs experimental techniques for the measurements of noise and vibrations (Panza 2015), methods based on psychoacoustics to evaluate the human perception of the discomfort (Griffin 2007), and strategies to correlate numerical and experimental results (Abdullah et al. 2017).
Several solutions have been proposed to improve the NVH experienced by the occupants of the cabin. Most of the times, the attention is focused on the source of noise and vibrations (Qatu 2012). In the past, this was mainly represented by the engine and the powertrain. Nowadays, with much quieter, the focus is mainly shifted to other sources of noise, such as road excitation and tyre performance (B€ acker, Gallrein, and Roller 2016;Uhlar, Heyder, and K€ onig 2021). Nevertheless, many times the most effective way of improving the NVH performance is to act on the material and geometric properties of the global body-in-white (BiW) structure. In fact, the NVH response is particularly sensitive to changes in the design parameters, as it strictly depends on the static and dynamic global stiffness of the vehicle body structure.
The dynamic properties, which provide important information about the vibrational behavior of the BiW, are usually obtained in terms of natural frequencies and mode shapes by performing the standard FE modal analysis. The standard FE static analysis, instead, allows to evaluate the global static stiffness and extract indicators of the ride comfort.
A common approach employed by the automotive industry to find the optimal BiW design, with respect to specific targets, is to perform optimization studies based on the design of experiments (DoE). This consists of generating a random set of possible design configurations and evaluating their performance by using full-order simulations. The obtained results are then used to build a statistical surface response (also called surrogate model), which approximates the behavior of the structure, and it can be used to perform efficient optimization studies (Kiani and Yildiz 2016). Due to the large number of configurations to be tested in real applications and the high cost of each static and dynamic simulation involved, standard DoE studies are computationally too expensive to be performed in the early phase of the design process. For this reason, the NVH analysis is traditionally performed at a later stage, risking to encounter late undesired issues. To speed up the design cycle and to avoid the unsustainable waste of time and resources in prototyping non-optimal products, the automotive industry is highly interested in new software (De Cuyper et al. 2007) and advanced simulation-based methodologies able to predict the vehicle performance already during the preliminary phase of the development process (Jans et al. 2006).
To this end, significant efforts have been placed in developing new techniques that substitute the detailed FE model, formed by shell elements, with an equivalent simplified FE model of the BiW characterized by beams, joints and panels (Qin et al. 2017;Mundo et al. 2009; Van der Auweraer et al. 2005). This approach leads to a significant reduction of the degrees of freedom, such that the calculation time of each simulation is much less and fast optimization studies can be performed. On the other hand, these methods require a preprocess step where full-order simulations are performed to calibrate the properties of the equivalent simplified model. Moreover, the transition from the optimized simplified model to the optimized detailed FE model is not straightforward and it generally introduces errors.
A valuable alternative which allows to work with detailed FE models, whilst reducing the computational effort, is to employ reduced order modeling (ROM) techniques. The basic idea behind ROMs is that the essential behavior of complex systems can be accurately described by simplified low-order models. Typical reduced techniques available in the literature (e.g., Krylov-based methods (Freund 2003), POD (Feeny and Kappagantu 1998), reduced basis (Rozza, Huynh, and Patera 2007)) are projection-based methods which aims at computing a reduced set of basis, able to accurately approximate the solution. These methods are also referred as a posteriori ROM techniques, as they typically need to first solve the full-order problem for a set of representative design configurations, based on which the set of reduced basis is built. The obtained reduced model can be then used to evaluate the response of new configurations at a significantly lower computational cost.
Several works have been recently proposed which combine ROM techniques with surrogate models. Such choice is particularly advantageous when optimization studies are performed in the context of crash simulations (Rocas et al. 2021(Rocas et al. , 2022Le Guennec et al. 2018), where designers have to deal with complex parametric problem involving material and geometric design variables, high displacements and material non-linearities. As already mentioned, the noise and vibration properties of a BiW structure are typically evaluated by performing standard linear static and dynamic FE analyses. A ROM technique which proved to work very well in the context of linear problems is the proper generalized decomposition (PGD) method (Ammar et al. 2006;Chinesta, Ladeveze, and Cueto 2011;Chinesta, Keunings, and Leygue 2014).
Unlike standard a posteriori ROM approaches, that require special attention on how to select appropriate snapshots to obtain accurate reduced models, the PGD is an a priori technique that computes the reduced basis without relying on previously computed full-order solutions associated with arbitrary samples of the parametric space. This is possible thanks to the main assumption of the method, that is to treat the parameters as extra-coordinates and approximate the solution of the resulting high-dimensional problem as a sum of "rank-one" terms. Each of these "rank-one" terms is given by the product of basis functions depending explicitly on the coordinates of the problem (spatial and parametric coordinates). This compact approximated expression of the solution, also known as computational vademecum, is particularized, in the online phase, for any set of the design variables at a negligible computational cost and with very low computational resources. As a consequence, it can be uploaded on light computational devices (such as tablets or smartphones), such that the visualization of the results, optimization studies or inverse analysis can be performed in real-time. In the last years, the method has been successfully tested in the most diverse areas of application, such as flow problems (Dumon, Allery, and Ammar 2011;Leblond and Allery 2014;Ib añez et al. 2017;D ıez, Zlotnik, and Huerta 2017;Garc ıa-Blanco et al. 2017;Giacomini et al. 2021), thermal problems (Ghnatios et al. 2012;Aguado et al. 2015;Huerta, Nadal, and Chinesta 2018), solid mechanics (de Almeida 2013; Reis et al. 2020), fracture mechanics (Giner et al. 2013;Garikapati et al. 2020), geophysical problems (Zlotnik et al. 2015;Signorini, Zlotnik, and D ıez 2017), elastic metamaterials, coupled magneto-mechanical problems (Sibileau et al. 2018;Barroso et al. 2020) and dynamic problems (Quesada et al. 2014;Gonz alez, Cueto, and Chinesta 2014;Germoso et al. 2016;Malik et al. 2018;Quaranta et al. 2019).
Having in mind the big potential that the PGD method could exhibit in the context of a NVH study, the authors recently developed an extended version of the standard PGD to perform the static (Cavaliere et al. 2021) and dynamic (Cavaliere et al. 2022) stiffness analysis of a 3D structure characterized by material and geometric design variables. In the static case, the PGD was coupled with the inertia relief method (PGD-IR), which is an approach widely used by the industry to solve the static analysis of unconstrained structures, such as BiW structures. In the dynamic case, the PGD was coupled with the inverse power method (PGD-IPM) to develop a parametric eigensolver able to extract the lowest natural frequencies and corresponding mode shapes. In addition, the authors tried to overcome a typical drawback of the standard PGD, that is its challenging application to geometrically parametrized models. This is due to the fact that one essential requirement of the PGD method is that the input quantities of the problem (such as the stiffness and mass matrix) are to be expressed with an explicit and separated dependence on the parameters. This is clearly not an easy task, especially when complex shape parametrisations, like the one which are typical in the automotive context, are considered. Moreover, even if an explicit expression could be found, the standard PGD would require an intrusive modification of the FE formulation of the problem, as described in several works (Leygue and Verron 2010;Bognet et al. 2012;Heuz e, Leygue, and Racineux 2016;Courard et al. 2016;Chamoin and Thai 2019;Sevilla, Zlotnik, and Huerta 2020). Clearly, this would make the method unemployable in an industrial context, where commercial software with inaccessible source codes are used and complex shape parametrisation are unavoidable. to overcome these limitations, the authors developed (Cavaliere et al. 2021) a hybrid algebraic approach, which requires a preprocess step where the input data are only sampled (assembled), without solving any full-order problem, and then expressed in the separated format. This approach preserves the nonintrusivity of the methodology, as it requires an interaction with the commercial software but does not need any modification of the FE formulation.
The goal of this work is to validate the feasibility of the proposed PGD-IR and PGD-IPM method in an industrial context, and not only with respect to a simple academic 3D structure as in the previous works. To achieve that, it was essential to construct a parametrized FE model of the structure by using morphing and optimization tools, such that typical BiW topology parametrisations can be easily handled, keeping the method nonintrusive. Moreover, the potential of the proposed method in the post-process phase is shown. First, a multi-objective optimization study is presented to find a set of optimal design configurations with respect to the NVH performance. Finally, the most attractive feature of the method is presented, which is the development of graphical interface app that allows to visualize in real-time how changes in the design variables affect the global response of the BiW. It is important to underline that this represents a powerful instruments which empowers the designers, who can discard suboptimal solutions already from the preliminary phase of the development process, alleviating the work load during the subsequent phases of the development process.
The remainder of the paper is structured as follows: Sec. 2 reviews the standard FE formulation for the static and dynamic global stiffness analysis of BiW structure, that is the reference full-order problem. In Sec. 3, the problems are redefined in the corresponding parametric framework and the basics of the encapsulated PGD toolbox are introduced. The core of the paper is in Sec. 4, where the NVH study of a BiW structure with material and geometric parameters is finally performed. In particular, three phases are described: the preprocess which concerns the parametrisation of the model and the preparation of the PGD input data; the offline computation, which solves the parametric static and dynamic problem by means of the PGD-IR and PGD-IPM algorithms and discusses the results; finally, the post-process is presented which consists of an optimization study and the development of the graphical interface app for the real-time visualization of the results. Section 5 gives the final conclusions.

Standard FE formulation of the NVH analysis
In real applications, the noise and vibration properties of a BiW structure are usually evaluated by performing standard static and dynamic FE analyses. Differently from other disciplines of the automotive development process, such as crashworthiness, geometric and material non-linearities are not considered during the NVH study.
The dynamic study consists of extracting the lowest natural frequencies and corresponding vibrational modes by means of the modal analysis. This allows to identify and optimize the first torsional and bending modes, which are good indicators of the BiW vibrational behavior. The ride comfort properties, instead, are related to the static torsional and bending global stiffness of the BiW structure, which is usually considered in its free-free (unconstrained) configuration to reproduce realistic conditions. As it is well known, if no boundary conditions are imposed, the structure undergoes rigid body motions and the static analysis cannot be performed due to the singularity of the stiffness matrix. A method that is widely used in industry to solve this issue is the inertia relief (IR) technique (Wijker 2004), which uses the rigid body modes to equilibrate the system and make it solvable.
From a mathematical point of view, the IR analysis can be seen as a limit case of the dynamic problem for the natural frequency tending to zero, thus the vibrational modes correspond to the rigid body modes. As a consequence, both the modal and the IR analyses are derived from the elastodynamic problem formulation, which in the discretized form and in absence of damping reads M € UðtÞ þ KUðtÞ ¼ FðtÞ: (1) Here, K and M represent the stiffness and mass matrices, while UðtÞ, € UðtÞ and FðtÞ are, respectively, the time-dependent displacement, acceleration and nodal vector of applied forces.
In the following, Eq. (1) is particularized to derive the standard formulation of the modal analysis and the IR method, which represent the basis of the parametric NVH study proposed in this work. It is worthy to mention that, for the sake of simplicity, this work focuses on the evaluation of the torsional static and dynamic properties. Nevertheless, the extension to the bending stiffness or other quantities of interests (QoIs) is straightforward.

Dynamic case: the modal analysis
The basic equation of the standard modal analysis is a reduced form of the elastodynamic problem in Eq. (1) under free vibration conditions, namely when no external loads are applied M € UðtÞ þ KUðtÞ ¼ 0: By following the standard procedure, the solution UðtÞ is assumed to have a harmonic dependence on time through the expression where x is the angular velocity, w is the phase angle and / is the deformed shape of the structure. The acceleration € UðtÞ is obtained by differentiation, and after carrying out the appropriate substitutions, Eq. (2) transforms into Equation (4) represents a generalized eigenvalue problem, which provides n dof eigenmodes / i , for i ¼ 1, 2, :::, n dof , being n dof the number of degrees of freedom of the structure. These eigenmodes are called natural mode shapes, as they represent the deformation of the structure when vibrating in its ith mode. Each eigenmode is associated with the eigenfrequency x i through the Rayleigh quotient while the natural frequency f (measured in Hertz) depends of the eigenfrequency through the expression If K and M are symmetric, M is positive definite and K is at least semi-positive definite, some mathematical properties hold which ensure the uniqueness of the eigenmodes. First, modes associated with different frequencies are orthogonal with respect to both M and K. That is / > i M/ j ¼ 0 and / > i K/ j ¼ 0 for i and j such that x i 6 ¼ x j : Moreover, they are normalized with respect to the mass M such that / > i M/ i ¼ 1, and consequently / > i K/ i ¼ x 2 i , for i ¼ 1, 2, :::, n dof : To solve the eigenvalue problem of Eq. (4), a numerical eigensolver needs to be employed. Several methods are available in the literature (Bai et al. 2000;Quarteroni, Sacco, and Saleri 2010;Golub and Van Loan 2013). Usually, the choice depends on the mathematical structure and the number of eigenpairs of interest. For instance, commercial software MSC-Nastran (Nastran 2004), which is used as a reference in this work, provides seven different methods for the real eigenvalue analysis. Among them, the Lanczos method (Lanczos 1950) is usually recommended for this case of study.
In this work, the choice is driven by two main reasons. First, the goal is not to compute the whole set of n dof eigenmodes, but a few fundamental modes corresponding to the lowest eigenfrequencies, which provide the essential information needed to assess the structural dynamic response. Moreover, the objective is the extension of the modal analysis to the parametric framework, which benefits from using the algorithmically simplest eigensolver available. Among others, the well-known inverse power method (IPM) fulfills all these requirements.
The IPM is an iterative algorithm, allows to find an approximation of the eigenvector / 1 associated to the smallest eigenfrequency x 1 by iteratively solving the system for ¼ 1, 2, 3::: starting with an initial guess / 0 1 : At every iteration , the quantity / 1 is normalized with respect to its M-norm given by ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi , such that unicity is preserved. After a certain number of iterations, the vector is expected to converge to the eigenmode / 1 : The corresponding eigenvalue x 2 1 can be easily computed according to the Rayleigh quotient, as indicated in Eq. (5). Subsequent eigenpairs are usually computed by employing a deflation technique. In standard algorithms, this is usually done by removing the already computed eigenvectors from the original matrix, while keeping the others unchanged.
Here, for the sake of easing the generalization to the parametric case, an equivalent strategy that uses Lagrange multipliers, introduced in (Cavaliere et al. 2022), is adopted. Instead of reducing the original system (7), an extended system of equations, of dimension ðn dof þ 1Þ Â ðn dof þ 1Þ, is solved where k is the Lagrange multiplier. By construction, the new sought eigenvector / 2 is automatically orthogonal to the previous one. When the first n eigenfrequencies, x 1 , x 2 , :::, x n , are already computed and the corresponding eigenvectors are collected in the n dof Â n matrix U n ¼ / 1 , / 2 , :::, / n ½ , the next eigenmode / nþ1 is obtained by solving the ðn dof þ nÞ Â ðn dof þ nÞ system of equations where k is the n Â 1 vector of Lagrange multipliers. An important aspect to consider in this work is that, since the BiW structure is assumed to be in the unconstrained configuration (with no loads and no constrains), its stiffness matrix K is singular. This means that the first six (in 3D) eigenmodes coincide with the rigid body modes of the structure and the corresponding first six eigenfrequencies are null. As a consequence, the strategy adopted is to first compute the rigid body modes and collect them in the n dof Â 6 matrix U 6 : Then, the following eigenvectors can be simply computed by solving Eq. (9). The rigid body modes can be constructed geometrically or they can be computed analytically, as explained in Sec. 2.3.

Static limit case: the inertia relief method
The setup of a standard FE analysis to evaluate the global static torsional stiffness of a BiW structure is shown in Fig. 1 (left picture). The test consists of loading the BiW model with a couple of parallel and opposite forces applied at the front and rear shock towers of the car frame, such that the resulting torsional moment is equal to 1 Nm. The QoI of this problem is the equivalent torsional stiffness (ETS), which is calculated as a function of the front and back twisting rotations of the car body when the constant torque is applied, namely Here the two twisting angles a AB and a CD represented in Fig. 1 (right picture) are given by the following relative vertical displacements where u z ðPÞ denotes the displacement in the vertical z direction at point P and L PQ denotes the distance between the points P and Q. The BiW is assumed to be in its free-free configuration, such that realistic results can be reproduced. As is well known, this causes a rigid body motion that prevents the application of the standard static analysis. The inertia relief (IR) method is usually employed to overcome this issue. The approach is based on the idea that, due to the mass of the system, the rigid body acceleration generates an inertial load that counteracts the external forces and deforms elastically the body. If the inertial load is such that the system is statically equilibrated (resultant forces and moments are zero), then the static analysis can be performed by applying any set of isostatic constraints.
The basic equation of the IR method is derived from the elastodynamic problem (1), assuming that the applied load F and the resulting displacement field U are constant in time, such that € U represents the steady state rigid body acceleration To solve the system of Eq. (12), the inertial load, M € U, that balances out the external forces, F, has to be computed. To this end, the rigid body acceleration is first expressed as a linear combination of the rigid body modes where U represents the rigid body matrix (each column is a rigid body mode) and the vector a contains the corresponding acceleration coefficients. Assuming that U is known, the acceleration vector a can be finally calculated by pre-multiplying Eq. (12) by U > and imposing that the system is equilibrated, that is Equation (14) leads to the final expression of the unknown acceleration vector a which provides the equilibrated forces With the system of equations given in Eq. (12) being equilibrated, the solution U can be found by imposing any set of isostatic constraints and solving Once the displacement field is known, the sought ETS can be easily calculated according to Eqs. (10) and (11).

Rigid body modes computation
By definition, a rigid body mode is a displacement vector that does not cause any static force. As a consequence, a simple way to analytically compute the six (in 3D) rigid body modes U ¼ ½/ 1 , / 2 , :::, / 6 is to impose their definition given by the system KU ¼ 0: First, n r reference degrees of freedom (n r ¼6 in 3D) are chosen such that, if they were constrained, no rigid body motion would occur. Then, the total number of n dof degrees of freedom is then partitioned into a reference set, r, and the remaining set, l, such that the matrix U can be computed by solving the partitioned version of KU ¼ 0, namely As K ll is symmetric and positive definite, U l set can be expressed in terms of U r : where U r is usually assumed to be the identity matrix of dimension n r Â n r , such that each column of the matrix U r represents a unit translation or rotation in the direction of the corresponding reference degrees of freedom. As mentioned above, the computation of the rigid body modes is the first common step of both the IR and IPM algorithms. Figure 2 gives a global overview of the static and dynamic analysis described above, which serves as starting point to extend the same algorithms to the parametric framework.

Reduced-order solver for the parametric NVH analysis
Let us consider a set of n p parameters, denoted by l ¼ ½l 1 , l 2 , :::, l n p > 2 M & R n p , describing the material properties (e.g., elastic modulus, density, etc.) and the geometric characterization of the BiW shape. Each parameter l j belongs to a predefined interval M j , such that the multidimensional parametric domain M is defined as the Cartesian product of the sectional intervals, From a conceptual point of view, the extension of the IR and IPM algorithms from the nonparametric to the parametric framework is as simple as rewriting all the quantities outlined in Fig. 2 with their parametric dependency. Standard numerical methods (e.g., finite elements, finite volumes, finite differences) would generally require to solve the problem for every possible combination of the parameters, which is usually not affordable for real industrial applications.
The alternative considered here is to treat the parameters as extra-coordinates. This means that the parametric input data, KðlÞ and MðlÞ, and the generalized solutions of the IR and IPM algorithms, UðlÞ and / i ðlÞ, are defined in the multidimensional domain D ¼ X Â M, where X and M denote the spatial and parametric domains, respectively. If standard mesh-based methods are adopted, the computational cost increases exponentially with the discretized parametric meshes, leading to the curse of dimensionality phenomenon. If, for example, n p ¼ 5 parameters are considered and each parametric space is discretized with m ¼ 10 nodes, the total number of degrees of freedom would be equal to n dof Â m n p ¼ n dof Â 10 5 : ROM techniques are usually employed to circumvent this issue by simplifying the numerical complexity of the problem. In this work, the underlined high-dimensional NVH problem is solved by means of the PGD method, which is introduced next.

PGD fundamentals
The PGD technique is based on the idea that the solution of a problem depending on a set of parameters can be approximated as a finite sum of terms. Each of these terms is given by the product of functions depending separately on the coordinates of the problem (spatial and parametric coordinates).
Let us consider the parametric solution UðlÞ of the IR problem. According to the PGD assumptions, it is approximated by 2 Þ ::: u i n p ðl n p Þ, where the vectors U i represent the spatial dimension, whereas the functions u i j ðl j Þ represent the dependency of the solution on each parameter. If the spatial and parametric functions are Figure 2. Overview of the NVH study in the non-parametric framework. The static analysis is performed by means of the IR method. The IPM algorithm is shown for the computation of the lowest (next) non-zero eigenvector / i , assuming that the matrix U of already computed modes is available. Conceptually, if all the input and output are expressed with their parametric dependency, the parametric version of the algorithm is identical. normalized, then an amplitude b i U is introduced, which indicates the relevance of the ith term of the sum on the final solution.
The solver scheme based on the PGD method consists of two main steps, as outlined in Fig. 3. A greedy algorithm is used to enrich the final solution by computing sequentially every term. The number of terms needed to reach a good accuracy is unknown a priori, and the enrichment automatically stops when a user-defined accuracy is reached.
To find the unknown spatial and parametric functions of each new term, an iterative scheme needs to be employed due to the non-linearity introduced by the products of unknown spatial and parametric functions in Eq. (19). The alternating direction scheme is normally employed in PGD problems, due to its simplicity and robustness. It consists in sequentially solving the problem for each unknown function, assuming that all the others are known, until a stationary solution is reached. It is worth emphasizing that, despite the introduction of the non-linearity, the computational cost of the problem increases linearly with the number of introduced parameters and not exponentially as it would be with standard numerical methods.
The output U PGD ðlÞ of the PGD solver, given in Eq. (19), is often defined as a generalized explicit solution, in the sense that it contains the solution for any combination of the parameters in the predefined design space. The most important advantage of the method is that, once the parametric output is available, the variation of the solution with respect to the parametric changes can be evaluated in real-time by simply performing a linear combination of terms, thus at a negligible computational cost.

PGD operations: the encapsulated PGD toolbox
The parametric IR and IPM algorithms require several parametric algebraic operations, such as products, sums, divisions, square roots. In this work, the encapsulated PGD Toolbox (D ıez et al. 2020) is utilized to perform such operations. The toolbox is a collection of PGD-based routines able to solve basic algebraic operations between parametric objects (i.e., scalar, vectors, matrices). The routines are implemented in an open-source Matlab environment and can be downloaded at https://git.lacan.upc.edu/zlotnik/algebraicPGDtools. The most important advantage of the toolbox from the point of view of a user, is that it works as a black box. This means that the user just needs to define the parametric input data and call the algorithm that performs the desired operation, which automatically gives back the parametric solution. This is an important feature for its use in an industrial context. The essential requirement of the proposed toolbox, as for any other PGD-based methods, is to define the input quantities in the separated format described in Eq. (19). However, being able to find a separated analytical representation of the input data is not always possible, especially when geometric parameters are considered in the problem. Furthermore, if an analytical expression existed, it would require to access the source code of the FE software and modify the problem formulation, as described in several works (Leygue and Verron 2010;Bognet et al. 2012; Heuz e, Leygue, and Racineux 2016;Courard et al. 2016;Chamoin and Thai 2019;Sevilla, Zlotnik, and Huerta 2020). Such an "intrusive" approach is not useful in practice in the industrial context, where commercial software with inaccessible source codes is typically employed.
For this reason, the authors developed a nonintrusive algebraic approach described in Cavaliere et al. (2021), which is able to find a separated expression of the input data for any material or geometric parametrisation. The main idea is to first sample the input data (without solving the problem) for every possible combination of the parametric values of interest and then express them in the separated format. In the presented static and dynamic problems, the stiffness and mass matrices are the required input data.
If each parametric dimension l j 2 M j for j ¼ 1, 2, :::, n p is discretized using m j nodal values, then the sampling of the parametric matrices consists of evaluating K i and M i in the set of m tot points used to discretise the parametric domain M :¼ M 1 Â M 2 Â Á Á Á Â M n p , where m tot ¼ Q n p j¼1 m j : This might be perceived as a computationally expensive step. However, it should be noted that the computational cost of assembling the FE matrices is small compared to the cost of a full-order simulation. As a result, this technique preserves efficiency, especially because the assembly process for different parametric values can be done in parallel.
Once the matrices K i and M i are sampled for the i ¼ 1, 2, :::, m tot combinations, the corresponding parametric functions k i j ðl j Þ and m i j ðl j Þ are defined such that Since m tot is usually a high number, a data compression based on the PGD rationale (D ıez et al. 2018) is employed. The compression routine, that is also included in the toolbox, is able to accurately approximate the original quantities with a much smaller number of PGD terms, such that the computational cost of the next operations can be substantially reduced.

Parametric version of the NVH algorithm
Based on the description of the encapsulated PGD toolbox given above, this section introduces the proposed parametric IR and IPM algorithms. In summary, the first step consists in sampling the input stiffness and mass matrices and expresses them in the PGD separated format. Then, the algebraic parametric operations involved in both algorithms are performed by simply calling the appropriate encapsulated PGD routine. The output of each operation is automatically defined in the separated PGD format, such that it can be directly used as input of the next operation. Once all the steps are performed, the final parametric solutions U PGD ðlÞ and / PGD i ðlÞ can be stored and used in the post-process phase. Figure 4 shows a global flowchart describing the sequence of assembled equations which need to be solved. Each algebraic operation is performed by employing the corresponding PGD function from the encapsulated PGD toolbox. For each step, the involved parametric operations are specified such that the reader can understand the complexity and the potential of the toolbox.
It is worth noting that certain operations (e.g., sum, products, concatenation) are as simple as adding, multiplying or concatenating the PGD terms defining each quantity. In contrast, operations like the square root, division or compression required the development of new algorithms based on the PGD framework (Cavaliere et al. 2022;D ıez et al. 2018). From an implementation point of view, an overloading of the arithmetic operators allows to simply use the standard Matlab symbols to call the algebraic operations contained in the encapsulated library, making the method highly user-friendly.
It is important to mention that the linear solver, the division, the square root and the compression PGD operations need two kind of tolerances to be set up. One tolerance for the iterative solver scheme, and another to stop the enrichment of the PGD terms. The compression, which is performed every time the number of PGD terms undergoes a substantial increment (e.g., after products), usually needs stricter tolerances, such that the accuracy is not compromised. From Fig. 4, it can be observed that the PGD-IR algorithm has a much simpler structure than the PGD-IPM. In fact, the IPM solver needs several iterations to be performed until the solution convergences to the eigenmode. Moreover, this has to be repeated for every new eigenmode of interest. Despite the approximations introduced by each PGD operation, the example in the next section will show a successful application to a realistic industrial case study for both the static and dynamic analyses. In particular, the study considers a simplified BiW structure characterized by four parameters, namely the thicknesses and cross sections of the C-pillars and the rear long member.

Numerical application: parametric NVH study of a BiW structure
The proposed method is finally applied to perform the NVH study of a simplified BiW structure. Figure  5 shows the geometry and the FE mesh of the model, which is formed by 3, 819 nodes, each one supporting six degrees of freedom (three translations and three rotations). Isoparametric triangular and quadrilateral elements based on the Mindlin-Reissner shell formulation (CTRIA3 and CQUAD4 respectively in MSC-Nastran) are used to discretise the model. All the car components are characterized by isotropic linear elastic materials (MAT1 in MSC-Nastran) with properties described in Table 1.
In this example, four parameters are introduced as extra-coordinates of the problem, which are the thicknesses and the cross sections of the C-pillars and the rear long members shown in the right picture of Fig. 5. Note that, although from a physical point of view the thickness is as a geometric parameter, in the shell element formulation it is treated as a material property. The cross section clearly represents a geometric design variable.
The proposed methodology consists of the following three stages: 1. Preprocess: A FE model of the BiW is constructed and parametrized such that the input data of the multidimensional problem can be sampled and expressed in a parametric format; 2. Offline computation: The static and dynamic global stiffness analysis are performed in an offline stage by means of one computation which uses the encapsulated PGD toolbox to obtain the parametric results; 3. Post-process: The parametric solution can be used for several purposes, such as efficient optimization studies and real-time evaluation of the results.  These three stages are described in detail in the following sections.

Preprocess: Parametrization of the BiW model
The preprocess starts with the preparation of a parametrized FE model of the BiW structure. In this work, this is done by means of the ANSA CAE preprocessor software, which contains a powerful optimization task tool able to organize the set-up of an optimization study. The first step consists of defining the design variables and their range of variation. In particular, the design variables of the problem are denoted by l ¼ ½l 1 , l 2 , l 3 , l 4 , where l 1 and l 2 represent the cross sections of the C-pillar and the rear long member, while l 3 and l 4 are the thicknesses of the same components. The first two geometric parameters are defined by means of a morphing tool available in ANSA, which allows to manage the shape changes. More precisely, the position of the nodes changes without changing the element connectivity. Each BiW component affected by the change is selected and subdivided into an action area and a transition area as depicted in Fig. 6. The action area is actively affected by the change in the cross section, while the transition part is used to smooth the deformation. A reference cross section with variable height and width of the action area can be chosen to guide the geometry variation along the whole component. In this example, the shape of the section is preserved, which means that the width and height vary together in a user-defined range. The reference undeformed cross section is associated to a current value of the design variable equal to 0. The cross section of the C-pillar is assumed to vary its height and width in a user-defined range of M 1 ¼ ½À20, 20 mm. Analogously, the rear long member cross section is defined such that the correspondent design variable l 2 varies in the range M 2 ¼ ½À10, 10 mm.
The material variables are parametrized by selecting the property of interest in the optimization tool and defining the ranges of variations. In this case, the variables l 3 and l 4 , representing the thicknesses of the two parametric BiW components, are defined in the ranges M 3 ¼ ½1:4, 1:8 mm and M 4 ¼ ½0:5, 1:3 mm, respectively. Table 2 summarizes the design variables definition. Figure 6. Definition of the geometric design variables in Beta-Ansa. The cross section of the components changes in the action area, the modification is smoothly absorbed through a transition area (left). The reference cross section changes its height and width while keeping the shape (right). Each of the four parametric spaces is discretized by means of nine equidistant nodal values, which means that the total number of parametric combinations is given by m tot ¼ 9 4 ¼ 6, 561 different configurations.
According to the description given in Sec. 3.2, the next step requires a sampling of the input data (mass and stiffness matrices) for each combination. To this end, a list with all the 6,561 parametric combinations can be uploaded into the ANSA optimization tool. Then, a design of experiments study is set up to automatically generate the input files in the format of the desired commercial software, which in this work is MSC-Nastran.
By using a special Nastran language (DMAP), all the input files generated by the design of experiment study are run such that the mass and stiffness matrices are just assembled and stored, without solving any static or dynamic problem. The stored files are then uploaded into Matlab and expressed in the PGD format introduced in Eq. (20). To finalize the preprocess phase, a data compression is performed to reduce the number of PGD terms. In this example, after performing the PGD compression imposing an accuracy of 10 À8 , the number of terms needed to approximate the stiffness K PGD ðlÞ and mass M PGD ðlÞ matrices reduces respectively to n K ¼ 66 and n M ¼ 20, instead of the initial 6,561 terms.

Offline computation: nonintrusive parametric NVH solver
This section shows the results of the proposed nonintrusive parametric NVH solver. The first common step of the PGD-IR and PGD-IPM algorithms is represented by the computation of the parametric rigid body matrix. As explained in the previous sections, the rigid body modes can be computed as the kernel of the stiffness matrix K PGD ðlÞ: First, a set, r, of reference degrees of freedom that represents a set of isostatic constraints is defined.
Once the set of total n dof is partitioned into the r and the left l-set of degrees of freedom, the encapsulated PGD toolbox is employed to solve the parametric version of Eqs. (17) and (18) to obtain the PGD rigid body matrix U PGD ðlÞ: Finally, the parametric static and dynamic solutions can be computed.

Static analysis
Once the rigid body matrix U PGD ðlÞ has been computed, the remaining steps of the PGD-IR algorithm shown in Fig. 4 can be performed by employing the encapsulated PGD routines. The final results are the parametric displacement vector which contains all the solution for all the combinations of the four parameters. In this example, N ¼ 75 terms were necessary to obtain the approximate PGD solution. In particular, the solution was considered sufficiently accurate when the amplitude b N of the last calculated PGD term was four orders of magnitudes smaller than the amplitude b 1 associated to the first PGD term.
To better illustrate the structure of a PGD solution, Fig. 7 shows the first three terms of the sum in Eq. (21). The amplitude factors b i , which are marked with a box for each term, are a measure of how much the ith term contributes to the final PGD solution. The spatial terms U i show the deformation induced by each PGD term. In the final solution, each spatial term is scaled by the parametric functions, which take a different value depending on the chosen set of parameters ½l 1 , l 2 , l 3 , l 4 : Once the parametric solution U PGD ðlÞ is available, it can be particularized in real-time for any combination of the parameters and the ETS can be easily calculated according to Eq. (10). To validate the proposed method, a comparison between the PGD solution and the standard FE solution is performed. For this comparison, the full-order problem was solved for all the 6,561 combinations by running the linear solver of MSC-Nastran in combination with the IR option to circumvent the singularity of the stiffness matrix. Note that, for this example it is feasible to compute the FE solution at every parametric point. However, if the number of parameters increases, computing all the FE solutions becomes unfeasible, whereas the proposed PGD approach is still a viable option.
The accuracy of the PGD with respect to the full-order computations is measured as the relative error between the PGD and Nastran ETS solutions in the L 2 ðM 1 Â M 2 Â M 3 Â M 4 Þ norm, that is and in terms of the maximum error, calculated as the L 1 ðM 1 Â M 2 Â M 3 Â M 4 Þ norm: The results obtained with the PGD method are in perfect agreement with MSC-Nastran, with a maximum error equal to 5 Â 10 À3 and an error measured in the L 2 -norm equal to 8:83 Â 10 À4 : To better understand how changes in the defined design variables affect the static response of the car, Fig. 8 shows the variation of the ETS with respect to the four parameters in terms of isosurfaces. In particular, each plot shows the variation in the parametric space defined by the first three parameters ðl 1 , l 2 , l 3 Þ when the forth parameter l 4 is fixed. The plots show a substantial range of variation in the ETS, which is around 6,000 Nm/degree. Moreover, it can be observed that the fourth parameter l 4 , corresponding to the thickness of the rear long members, has the biggest influence on the ETS variation. When l 4 increases, the isosurfaces tend to be parallel to the ðl 1 , l 2 Þ space, meaning that the variation due to changes in the cross sections gets reduced with the thicknesses increment.

Dynamic analysis
The mode shapes associated to the three smallest non-zero frequencies are computed using the proposed PGD-IPM eigensolver to identify the first torsional mode of the BiW structure under analysis.
The parametric input data of the problem is represented by the stiffness and mass matrices already sampled and expressed in the PGD format during the preprocess phase. Due to the unconstrained configuration, the matrix of rigid body modes U PGD ðlÞ is essential for the computation of the subsequent non-rigid eigenmodes. As shown in Fig. 4, a parametric guess vector is chosen and the IPM system of equations is iteratively solved by employing the encapsulated PGD toolbox. At every iteration, the obtained parametric eigenvector is normalized. The iteration stops when convergence is reached, that is when a quantity E / is smaller than a user-defined tolerance.
Here E / is defined as where b i / new and b j / old represent, respectively, the amplitudes of the PGD terms describing two eigenmodes obtained by two consecutive iterations.
Once convergence is reached, the sought ith eigenvector is obtained, in separated form, as The corresponding eigenfrequency can be calculated according to the Rayleigh quotient For the BiW structure considered here, the first and third non-rigid eigenvectors represent two different kind of bending modes, independently on the combination of the parameters. Similarly, the second non-rigid eigenvector always represents a torsional mode, which is the one of interest. Of course, when more complex models are analyzed, the order of the modes can easily change with the parameters, so their identification would represent an important task.
Once the sought smallest torsional eigenvalue x PGD t ðlÞ is available, the corresponding torsional natural frequency f PGD t ðlÞ can be computed as : As it was done in the static case, a comparison between the PGD solution and the standard FE solution of the modal analysis is performed. In this case, a full-order real eigenvalue analysis was performed in MSC-Nastran for all the 6,561 parametric combinations. Once again, the results are in perfect agreement. The relative errors are calculated as in Eqs. (22) and (23) by substituting ETS with the frequency, leading to a L 2 -norm error equal to 1:86 Â 10 À4 and a maximum relative error of 1:70 Â 10 À3 : Finally, the variation of the torsional frequency is also depicted in terms of isosurfaces in Fig.  9. Each plot refers to a fixed value of l 4 and shows the variation in the 3D Cartesian space defined by ðl 1 , l 2 , l 3 Þ: Also in this case, the thickness of the rear long member (described by l 4 ) proves to have the biggest influence on the QoI. The vertical character of the isosurfaces suggests that, differently from the static case, the thickness of the C-pillar (corresponding to l 3 ) represents the less influencing parameter. The plot also shows how small variations in the material and geometric parameters of the two components can lead to variations of the torsional frequency in the range of 2-3 Hz, which can substantially change the perception of vibration for the occupants of the vehicle.
It is worth emphasizing that the described algorithm is solved by means of just one offline computation. The resulting eigenpairs represent the compact version of all the possible solutions for every value of the parameters.

Post-process: optimization study and real-time visualization
One of the most interesting features of the PGD method is that, once the offline process is finished, obtaining the solution for a given value of the parameters takes 0.1 milliseconds. This opens the possibility to perform efficient optimization studies and visualize the results in real-time.
To show the potential of the proposed methodology, a multi-objective optimization analysis is presented next. The goal is to find the optimal combinations of the parameters such that the ETS Figure 9. Isosurfaces showing the variation of the smallest torsional natural frequency f 2 with respect to the parameters l 1 , l 2 and l 3 . Each plot refers to one specific value of parameter l 4 . and torsional frequency are maximized whilst the mass of the two car components is minimized. Three objective functions are defined as g 1 ðlÞ ¼ Mðl 1 , l 2 , l 3 , l 4 Þ, g 2 ðlÞ ¼ ETSðl 1 , l 2 , l 3 , l 4 Þ, g 3 ðlÞ ¼ f t ðl 1 , l 2 , l 3 , l 4 Þ, where Mðl 1 , l 2 , l 3 , l 4 Þ represents the total mass of the C-pillars and the rear long members, depending on their variable geometries and thicknesses. Clearly, this quantity is strictly related to the production cost. ETSðl 1 , l 2 , l 3 , l 4 Þ is the parametric output obtained by means of the proposed PGD-IR algorithm and f t ðl 1 , l 2 , l 3 , l 4 Þ is parametric torsional frequency calculated by the PGD-IPM eigensolver. The explicit dependency of the three functions on the parameters permits to easily compute the Pareto front of the multiple objective functions by means of a genetic algorithm. In this work, this is done by using the global optimization toolbox released by Matlab. The obtained Pareto points are shown in Fig. 10. According to the definition, they represent a tradeoff between the objective functions, meaning that each point is considered optimal if no objective function can be improved without compromising at least one other objective. The cloud of sampled points in the plot represents the mass and ETS coordinates corresponding to each of the 6,561 parametric combinations considered initially. It is clear that the optimization study allows to drastically reduce the number of configurations which would be considered by the designers in the final decision-making process.
It is worth noting that, in this example, the Pareto front was computed by assigning the same weight to both objective functions. Nevertheless, it is straightforward to obtain other fronts if the user wants to put more emphasis on one of the objective functions.
The ultimate goal of this work is the development of a computational tool which can be accessed by the designers to visualize the effect of the parametric changes on the global response of the BiW in real-time. To this end, a standalone desktop app has been built by using the Matlab App Designer software. Figure 11 shows the developed graphical user interface. An interactive version of the Pareto front depicted in Fig. 10 is shown in the app. The sliders in the red box allow to adjust the set of parameters to one of the 6,561 combinations given by the discretized parametric space. The corresponding point is simultaneously visualized in the Pareto Front plot, whilst the numerical values of the parameters, together with the mass and QoI coordinates, appear in the tables on the bottom right.
Since the optimal solutions obtained through the multi-objective optimization do not necessarily coincide with one of the 6,561 combinations represented by the cloud of points, an extra box is added which allows to scroll through a list of all the Pareto points (red box in Fig. 11). Also in this case, the point information is reflected in the Pareto front plot and the numerical values are updated in the corresponding tables. If the user is interested in visualizing all the Pareto points, it can be done by simply clicking on the "Display all optimal points" button. Finally, the static and dynamic deformation of the BiW corresponding to the selected parametric choice is updated in real-time in the left top plot.
To summarize the developed PGD method for the parametric NVH analysis, Fig. 12 offers an overview of all the steps described in this section.

Conclusions
The goal of this work was to propose a new ROM method to speed up the design process of a car structure. More precisely, a parametric global static and dynamic stiffness analysis was performed by means of a PGD-based methodology to optimize the NVH performance of a BiW structure characterized by material and geometric design variables. Two material and two shape design variables were studied, which correspond to the thicknesses and cross sections of the Cpillar and rear long member components. A coupling of the PGD with the IR method, which allows the static analysis of the BiW in its unconstrained configuration, was used to compute the parametric static equivalent torsional stiffness (ETS) of the BiW. The parametric lowest torsional frequency under free vibration condition was obtained, instead, by employing a coupled PGD-IPM (inverse power method) eigensolver.
The proposed technique makes use of the encapsulated PGD toolbox developed by D ıez et al. (2020), which enables to perform the algebraic operations between multidimensional data contained in the PGD-IR and PGD-IPM algorithms, as already presented by the authors in recent works (Cavaliere et al. 2021(Cavaliere et al. , 2022. The main contribution of this publication was to unify the PGD-IR and PGD-IPM in a global algorithm (also referred as PGD-NVH) and extend its application to a realistic case of study. To this end, a parametrized detailed FE model formed by shell elements was built in the ANSA CAE preprocessor commercial software, which allowed to deal with complex shape changes and prepare the input data of the problem in the separated PGD format.
Being able to deal with a detailed FE parametrized model represents an ideal solution for industry. In this way, designers do not have to select representative design configurations or develop equivalent simplified FE models, as it is usually needed by other methods. Moreover, the same parametric model could be used in other disciplines of the development processes, easing multidisciplinary design strategies. Once the parametric FE model was ready, the proposed parametric NVH solver was executed during an offline computation in an in-house Matlab environment, acting as a black-box, such that a nonintrusive interaction with commercial FE packages is possible. Finally, a post-process of the obtained parametric results was presented.
The accuracy of the method was measured by comparing the two quantities of interest (ETS and torsional frequency) with the corresponding full-order results, which resulted into a maximum relative error in the order of 10 À3 : It is important to mention that the PGD results were obtained by performing only one offline computation for the static and dynamic problems and then particularizing the results for any parametric combination in real-time. On the contrary, a total of 13,122 full-order simulations (6,561 for the static and 6,561 for the dynamic case) were needed to sample the results by means of standard methods. Another important observation is that computational time is not particularly significant in the PGD context, since the main goal is to provide a method which is able to explore an arbitrary large parametric space with only one offline computation. In other words, the PGD result describes the structural behavior for an infinite number of configurations in a predefined design space, which is difficult to compare to a standard design optimization study where a set of configurations need to be chosen. Nevertheless, in a previous work (Cavaliere et al. 2021), it was proven that PGD can save the 20% of computational time if compared to an equivalent standard static analysis. A future outlook is to apply the proposed method to a real automotive application. In that case, a realistic comparison of the proposed PGD-NVH efficiency with respect to standard NVH studies can be done.
Exploiting the explicit dependency of the QoIs on the parameters, a multi-objective optimization study was performed by using a genetic algorithm. The study allowed to identify a set of optimal Pareto points which drastically reduced the combination of design variables to take into account in the final decision-making process. Finally, a graphical interface app was developed by using the Matlab App Designer software, providing an interactive visualization of the results, such that the designers can check in real-time the effects of variables on the global static and dynamic behavior of the BiW structure. The developed app is just an example of the potential of this method. In fact, the information contained in the app could be modified and adapted to the needs of the specific problem, representing the kind of support that the industry urgently needs to optimize the development process.