Enhancing CFD predictions in shape design problems by model and parameter space reduction

In this work we present an advanced computational pipeline for the approximation and prediction of the lift coefficient of a parametrized airfoil profile. The non-intrusive reduced order method is based on dynamic mode decomposition (DMD) and it is coupled with dynamic active subspaces (DyAS) to enhance the future state prediction of the target function and reduce the parameter space dimensionality. The pipeline is based on high-fidelity simulations carried out by the application of finite volume method for turbulent flows, and automatic mesh morphing through radial basis functions interpolation technique. The proposed pipeline is able to save 1/3 of the overall computational resources thanks to the application of DMD. Moreover exploiting DyAS and performing the regression on a lower dimensional space results in the reduction of the relative error in the approximation of the time-varying lift coefficient by a factor 2 with respect to using only the DMD.

Reduced order modeling (ROM) is nowadays a quite popular and consolidated technique, applied to several fields of engineering and computational science thanks to the remarkable computational gain for the solution of parametric partial differential equations (PDEs).The ROM goal is to reduce the dimension of the studied system letting unaltered some important properties of the original problem, resulting in a more efficient computation.Such methods are frequently applied when many solutions for different parameters are required, for example in the context of parametric optimal control problems, uncertainty quantification, and shape optimization.
For parametric reduced order models, the most common approach is to sample the solution manifold by creating a solutions database corresponding to different parameters, using a high-dimensional discretization, then combine the latter to identify the intrinsic lower dimension of the problem.For parametric reduced order models see [21,40,41], while for a more applications oriented overview we suggest [49,42,43].
For parametric time-dependent problems, a proper orthogonal decomposition approach can be applied to reduce the dimensionality of the system, as in [17,23].In this work we propose a novel data-driven approach for parametric dynamical systems, combining dynamic mode decomposition (DMD) with active subspaces (AS) property.These two relatively new methodologies provide a simplification of the dynamical system, and an analysis of the input parameter space of a given target function, respectively.Exploiting AS property we are able to obtain an estimation of the importance of the parameters of such function, as well as a reduction in the number of parameters.Moreover the methods are equation-free, being based only on input/output couples and do not make assumptions on the underlying governing equations.
We define a generic scalar output s(µ, t) ∈ R that depends both on time t and on the parameters of the model µ ∈ D ⊂ R k , with k denoting the dimension of the parameter space.We denote the state of the parametric system at time t with s t (µ) ∈ R. The solution manifold in time is approximated using the DMD in order to obtain an approximation of the linear map A defined as: ( It is easy to note that using (1) we have the possibility to forecast a generic future state of the parametric system.
To numerically compute the linear operator A, we need to sample the parameter space D, and for each time store the quantity of interest for each parametric configuration.Formally, considering a set of parameter samples with dimension N s, the discrete vector referring to the system state at time t results: Collecting several time states s i (µ) for i = 1, . . ., m, we compute the operator A with a best-fit approach such that s t+1 ≈ As t .Once computed the future prevision, we are able to exploit the relation between the input parameters µ i and the related outputs s future (µ i ) to approximate the output for any new parameter.In this work we use a Gaussian Process Regression (GPR) [56,20], but any regression or interpolation method can be used.We underline that the chosen regression model has to be fitted for any forecasted time we want to analyse.
Since the high dimensionality in the parameter space can cause a decrease in the accuracy of the solution approximation, we couple the regression with the AS property in order to perform a sensitivity analysis of function s t (µ).AS indeed is able to provide an approximation g of a scalar function f , where the input parameters of g are a linear combination of the original parameters of f .The coefficients of such combination give information about the importance of the original parameters.In this work, we use this information to reduce the dimension of the parameter space -in which we build the regression -by not considering the parameters whose AS coefficients are smaller than a certain threshold, that is they are almost zero.
The developed methodology is tested on an aeronautics application given by the flow past an airfoil profile.As output of interest we considered the lift coefficient and the parameters vector µ describes geometrical transformations according to the morphing technique proposed in [22].The fluid dynamics problem is described using the incompressible Navier-Stokes equations with turbulence modeling.These are discretized using a finite volume approximation.The deformed meshes corresponding to different input parameters are automatically obtained exploiting a Radial Basis Function (RBF) mesh morphing technique.
This work is structured as follows: in section 2 we present the general parametric problem over which we apply the proposed numerical pipeline, providing some information about the geometrical deformation.In section 3 and section 4 we present the DMD and AS methods, respectively, while in section 5 we show the numerical setting of the problem and the results obtained.Finally in section 6 we propose some final remarks and highlight possible future developments.

The parametric problem
Let be given the unsteady incompressible Navier-Stokes equations described in an Eulerian framework on a parametrized space-time domain holds.Here, Γ = Γ in ∪ Γ 0 ∪ Γ out is the boundary of Ω(µ) and it is composed by three different parts Γ in , Γ out and Γ 0 (µ) that indicate, respectively, inlet boundary, outlet boundary, and physical walls.The term f (x) depicts the stationary non-homogeneous boundary condition, whereas k(x) denotes the initial condition for the velocity at t = 0. Shape changes are applied to the domain Ω, and in particular to its boundary Γ 0 (µ) corresponding to the airfoil wall.Such shape modifications are associated to numerical parameters contained in the vector µ ∈ R k which, in the numerical examples shown in this work has dimension k = 10.As said, the only portion of the domain boundary subject to shape parametrization is the physical wall of the airfoil Γ 0 (µ), which in the undeformed configuration corresponds to the 4-digits, NACA 4412 wing profile [1,25].To alter such geometry, we adopt the shape parametrization and morphing technique proposed in [22], where k shape functions are added to the airfoil profiles.Let y u , and y l be the upper and lower ordinates of a NACA profile, respectively.We express the deformation of such coordinates as where the bar denotes the reference undeformed state, which is the NACA 4412 profile.
The parameters µ ∈ D ⊂ R 10 are the weights coefficients, c i and d i , associated with the shape functions s i .The range of each parameter will be specified in section 5.The explicit formulation of the shape functions can be found in [22], we report them in Figure 1.After the reference profile is deformed, we also apply the same morphing to the mesh coordinates by using a radial basis functions (RBF) interpolation method [7,38,36].With this approach the movement s of all the points which do not belong to the moving boundaries is approximated by an interpolatory radial basis function: where x bi are the coordinates of points for which we know the boundary displacements, for this particular case the points located on the wing surface.N b is the number of control points on the boundary, ξ is a given basis function, q(x) is a polynomial.The coefficients β i and the polynomial q(x) are obtained by the imposition of interpolation conditions where d bi is the displacement value at the boundary points and by the additional requirement: In the present case, we select basis functions for which it is possible to use linear polynomials q(x).For more informations concerning the selection of the order of polynomials see [3].Finally the values of the coefficients β i and the coefficients δ i of the linear polynomials q can be obtained by solving the linear problem: where M b,b ∈ R N b ×N b is a matrix containing the evaluation of the basis functions ξ bibj = ξ( x bi − x bj ), and P b ∈ R N b ×(d+1) is a matrix where d is the spatial dimension.Each row of this matrix, that contains the coordinates of the boundary points, is given by row i (P b ) = 1 x bi .Once the system of ( 9) is solved one can obtain the displacement of all the internal points using the RBF interpolation: where x ini are the coordinates of the internal grid points.The computation of the displacement of the grid points entails the resolution of a dense system of equations that has dimension N b + d + 1. Usually, the number of boundary points N b is much smaller with respect to the number of grid points N h .

Dynamical systems approximation by dynamic mode decomposition
Dynamic mode decomposition (DMD) is an emerging reduced order method proposed by Schmid in [44] for the analysis of dynamical systems.Approximating the linear infinite-dimensional Koopman operator [28], DMD decomposes the original system into few main features, the so called DMD modes, that evolve linearly in time, even if the original system has nonlinear behaviour.This means that, other than individuating recurrent patterns in the evolution of the system, DMD provides a real-time midcast/forecast of the output of interest.An important advantage of such method is the complete data-driven nature: the algorithm relies only on the system output, without the necessity of any information regarding the model or equations used.Dynamic mode decomposition has been successfully employed in naval hull shape optimization pipelines [13], for online real-time acquisitions in a wind tunnel experiment [59], and in meteorology [4], among others.We also mention the higher order DMD extension [31,32].
In the following paragraph, we provide just an algorithmic overview of the method.For an exhaustive explanation of DMD, its applicability, and possible extensions, we suggest [29,6].
We define the linear operator A such that where x k+1 ∈ R N and x k ∈ R N are the vectors containing the system outputs at two sequential instants.Thus, the operator A : R N → R N expresses the dynamics of the system.In order to construct it using only data, we need to collect m equispaced in time outputs x i for i = 1, . . ., m -from now on called snapshots -then arrange them in two matrices: X = x 1 . . .x m−1 and Y = x 2 . . .x m .Since the corresponding columns in X and Y are sequential snapshots, we are able to use (11) to represent the relationship between X and Y, such that Y = AX.Minimizing the error Y − AX we obtain the linear operator, which however has very large dimension, especially when the studied system requires a fine discretization.To reduce the dimensionality, a POD approach is adopted.The matrix X is decomposed using the singular value decomposition as: where the matrix U contains the orthogonal left-singular vectors.We can then project the operator onto the space spanned by the left-singular vectors to get the reduced operator Ã.It is possible to note that the reduced operator does not require the construction of the high-dimensional one: where the † refers to the Moore-Penrose pseudo-inverse.We can now reconstruct the eigenvectors and eigenvalues of the matrix A thanks to the eigendecomposition of Ã as ÃW = WΛ.In particular the elements in Λ correspond to the nonzero eigenvalues of A, while the real eigenvectors, the so called exact modes [53], can be computed as Φ = YVΣ −1 W. Thus, being A = ΦΛΦ † , we can approximate the evolution of the system x k+1 = ΦΛΦ † x k .Moreover, it is easy to demonstrate that the approximation of a generic future snapshots can be computed as: In this work we compute the DMD modes of the matrix composed by the value of the time-varying lift coefficient for a set of given geometrical parameters.Then we can predict the future state of the coefficient and, using a regression method, approximate the target function at untried new parameters.All the DMD computation have been carried out by the Python package PyDMD [15].
Active subspaces have also been proven as a useful tool to enhance model order reduction techniques such as proper orthogonal decomposition (POD) with interpolation for structural and fluid dynamics problems [16], and POD-Galerkin methods for a parametric study of carotid artery stenosis [47].
Here we briefly introduce the active subspaces property for functions not depending on time, for the details and estimates regarding the method we refer to [8].For the actual computations to find AS we used the Python Active subspaces Utility Library [12].
Let µ ∈ R k the parameters of our problem, f be a parametric scalar function of interest f (µ) : R k → R, and ρ : R k → R + a probability density function representing uncertainty in the input parameters.Active subspaces are a property of the pair (f, ρ).They are defined as the leading eigenspaces of the second moment matrix of the target function's gradient and constitutes a global sensitivity index more general than coordinate-aligned derivative-based ones [58].
The second moment matrix of the gradients C, also called uncentered covariance matrix of the gradients of f with respect to the input parameters, is defined as where E[•] is the expected value, and C is symmetric thus it admits a real eigenvalue decomposition that reads: where W indicates the orthogonal matrix containing the eigenvectors of C as columns, and Λ is a diagonal matrix composed by the non-negative eigenvalues arranged in descending order.We can decompose the two matrices as follows where M < k has to be properly selected by identifying a spectral gap.In particular, we define the active subspace of dimension M as the principal eigenspace corresponding to the eigenvalues prior to the gap.Then we can map the full parameters to the reduced ones through W 1 .We define the active variable as , and the inactive variable as η = W T 2 µ ∈ R k−M .In practice the matrix C is constructed with a Monte Carlo procedure.
AS stipulates that the directional derivatives in directions belonging to the kernel of W T 1 are significantly smaller that those belonging to the range of the same matrix.Moreover this assumptions are made in expectation rather then in absolute sense [57].
Since in this way we are considering a linear combinations of the input parameters, we can associate the eigenvectors elements to the weights of such combinations, thus providing a sensitivity of each parameter.We underline that if a weight is almost zero, that means f does not vary along that direction on average.
We can use the active variable to build a ridge function g [33] to approximate the function of interest, that is In this work we want to study the behaviour of a target function f (µ, t) : R k × R + → R that depends on the parameters µ and on time t as well.This results in extending the active subspaces property to dynamical systems, that means having to deal with time-dependent uncentered covariance matrix C(t), and corresponding eigenvectors w i (t).Efforts in this direction has been done in [9] for a lithium ion battery model, in [34] for long term model of HIV infection dynamics, and more recently an application of dynamic mode decomposition and sparse identification to approximate one-dimensional active subspaces in [2].In these works they refer to dynamic active subspaces (DyAS) as the time evolution of the active subspaces of a time-dependent quantity of interest.
DyAS are useful to assess the importance of each input parameter at given times and to study how the weights associated to the inputs evolve.In the following we are going to compute the AS for a set of equispaced times t i .If some of the parameters are almost zero in the entire time window we can safely ignore them in the construction of the Gaussian process regression.

Computational pipeline
In the present section we will discuss the numerical experiments carried out to test the DyAS analysis and present the results obtained.As reported in section 2, each high fidelity simulation is based on a parametric fluid dynamic model governed by the Reynolds Averaged Navier-Stokes (RANS) equations.Thus, a number of flow simulations have been carried out selecting different samples in the parametric space to test the performance -in terms of lift coefficient -of different airfoil shapes.The simulations made use of both the RANS solver provided in the OpenFOAM [54] finite volumes library, and of the DMD acceleration methodology described in section 3. Once the lift coefficients output were available for all the samples tested in the input parameters space, the DyAS analysis was applied to assess possible parameter redundancy.The elimination of the redundant parameters detected in the DyAS analysis allowed for the generation of a surface response model based on a lower dimensional space, which has been finally tested against the original RANS model accelerated through DMD, and against the surface response model based on the original input parameter space.The following sections will further detail each part of the computational pipeline just outlined.

Parametric shape deformation
The fluid dynamics problem is resolved using the finite volume method.The wing is immersed in a rectangular domain according to Figure 2. The reference mesh counts 46500 hexahedral cells and is constructed using the blockMesh utility of the OpenFOAM library. Figure 2 depicts a detail of the grid in proximity of the wing.The meshes in the deformed configuration have been obtained starting from the reference configuration using a radial basis function smoothing algorithm similar to the one implemented in [5].A single deformation corresponds to a sample µ in the parameter space D := [0, 0.03] 10 ⊂ R 10 .Therefore all the deformed meshes share the same number of cells and the same mesh topology.In particular Wendland [55] second order kernel functions with radius r RBF = 0.1 m have been used.The control points of the RBF procedure have been placed on each mesh boundary point located onto the wing surface.Since the outer boundary points are fixed we decided to neglect them from the RBF computation using a smoothing function defined in such a way that the RBF contribution reduces to zero after a certain distance from a focal point [26].
Particularly, the focal point has been placed in the geometric center of the airfoil chord segment and the distance from the focal point after which the RBF contribution is neglected is set to r out = 7 m.In Figure 3 we depict the envelope of all the tested configurations, and the flow velocity streamlines for a particular sample in the parameter space.A uniform and constant velocity equal to u in = 1 m/s is set at the inlet boundary, while the constant value of the kinematic viscosity is set to ν = 2e−5 m 2 /s.This configuration, considering a chord length D = 1 m, corresponds to Reynolds number Re = 50000.As well known, a flow characterized by Reynolds number of such magnitude requires turbulence modeling to be numerically simulated with reasonable computational effort.In the present work, turbulence has been modeled using a RANS approach with a Spalart-Allmaras turbulence model [45].The pressure velocity coupling is resolved in a segregated manner making use of the PIMPLE algorithm which merges the PISO [24] and the SIMPLE [39] algorithm.The time step used to advance the simulation in time is set constant and equal to ∆t = 1e − 3 s.The convective terms have been discretized using a second-order upwinding scheme, while the diffusion terms are discretized using a linear approximation scheme with non-orthogonal correction.The time discretization is resolved using a second order backward differentiation formula.The simulation is advanced in time until the flow has reached stationary behavior.For the present problem, setting a total simulation time T s = 30 s is sufficient to reach a solution which is reasonably close to the steady state one.

Parameter space reduction
The present section will discuss the application of DyAS to the problem of the two dimensional turbulent flow simulation past airfoil sections with parameterized shape.Such a fluid dynamic problem is relevant in several engineering fields, as it is encountered in a number of industrial applications, ranging from aircraft and automotive design, to turbo machinery and propeller modeling.
A few plots describing the DyAS results for the lift coefficient output are presented in Figure 4, 5, 6, and 7.The plots in the figures are aimed at representing the evolution of the active subspace effectiveness and composition over the time dependent flow simulations.More specifically, the left diagram in each figure plots the lift coefficient at each sample point tested, as a function of the first active variable obtained through a linear combination of the sample point coordinates in the parameter space, that is f (µ, t) against W T 1 µ.Presenting the components of the first eigenvector of the uncentered covariance matrix, the right plot in each figure indicates the weights used in such linear combination to obtain the first active variable.In summary, the right diagram in each Figure suggests the impact of each of the original parameters on the first active variable, while the left diagram is an indicator of how well a one dimensional active subspace is able to represent the input to output relationship.Following the evolution of these two indicators it is possible, at each time instant, to assess how effective the one dimensional parameter dimension reduction is, and what is the sensitivity of the reduced lift coefficient output to variations of the original parameters.The plots in Figure 4, 5, 6, and 7 show the results of the DyAS at the fixed time instants t = 6, 10, 14, 18, respectively.A first look at the right plots for each time steps, suggests that the contribution of the parameters corresponding to the bump shape functions s 1 , and s 5 , for both the top and the bottom part of the airfoil profile are almost negligible.This means the lift coefficient is almost insensitive to variations of these 4 parameters.Alternatively, it can be said that the output function is on average almost flat along directions corresponding to the axes corresponding to parameters c 1 , c 5 , d 1 , and d 5 .
Figure 4 and 5 present the characterization of the one dimensional active subspace at time t = 6 s and t = 10 s, respectively.We can clearly see that the  lift coefficient is perfectly approximated along the identified direction, and such direction (the eigenvector elements) is almost the same at t = 6 s and t = 10 s.This should not completely surprise as both time instants are included in an initial acceleration phase during which the air coming from the inflow boundary is reaching the airfoil.Given the domain arrangement described in Figure 2, the flow velocity around the impulsively started airfoil leading edge is expected to reach the inflow value at time t = 10 s.For such reason, we will focus the description on the plots for t = 10 s, although the considerations can be immediately reproduced for previous time steps.The left plot in Figure 5 suggests that at this meaningful instant, the first active subspace represents the input to output relationship with remarkably good accuracy.In fact, only a single output value corresponds to each active variable value.In other words, when plotted against the first variable, the output appears like a curve -a line in the present case.A look at the right diagram suggests that the shape parameters having the most impact on the lift generated by the airfoil are c 3 , c 4 , d 3 and d 4 ,  which are the ones associated to shape functions with peaks located around the middle of the airfoil chord.The positive values of the eigenvector components associated to c 3 , c 4 , d 3 and d 4 , along with the positive slope of the curve in the left plot in Figure 5 suggest that, at this particular time instant, higher values of lift can be obtained by increasing the airfoil thickness in the mid-chord region.Similar considerations can be drawn from Figure 6, which refers the the DyAS analysis carried out at t = 14 s.Here, the points in the left diagram do not completely cluster on top of a single valued curve as was the case for the previous time step considered.Compared to what has been observed at t = 10 s, the data clearly indicate that at t = 14 s an input to output relationship obtained using only a one dimensional active subspace will lead to less accurate lift coefficient predictions.Yet, the points in the plot are still all located within a rather narrow band surrounding a regression line having positive slope.Thus, all the considerations on the lift coefficient sensitivity with respect to variations of the shape parameters that can be inferred from the right plot, will still hold at least by a qualitative standpoint.Here, the eigenvector components suggest that the most influential parameters on the lift coefficient are c 3 , d 3 and d 4 , while c 2 and d 2 affect the output in lesser but not negligible fashion.Compared to the previous case the importance of coefficient c 4 on the output is significantly reduced.We recall that c 4 is associated with increased y coordinates of the airfoil suction side past the mid-chord region.Thus, we might infer that in the acceleration phase higher lift values are obtained not only increasing the front thickness, but also lowering the camber line in the region past mid-chord.
Figure 7 shows the results of the DyAS analysis at t = 18 s, when the flow approaches the final regime solution.Following the trend observed for t = 14 s, the left plot in the figure indicates that a one dimensional is not completely able to represent the input to output relationship in a satisfactory fashion.With respect to the previous plots, the output values are here located in an even wider band around a regression line with positive slope.Again, on one hand this increasingly blurred picture suggests that higher dimensional active subspaces are required to reproduce the steady state solution with sufficient accuracy; on the other hand, the diagram still suggests a quite definite trend in the output, which can be exploited for qualitative considerations.Quite interestingly, at the present time step the eigenvector component corresponding to the c 4 coefficient has negative sign.Given the positive slope of the input to output relationship in the left plot of Figure 7, this implies that increases in the airfoil ordinates on the top side in the region past the mid-chord result in lift loss.Thus, this seems to suggest that an airfoil with a higher camber line curvature, combined with a thicker leading edge region might result in increased lift.This should not surprise, as a similar kind of airfoil would result in a higher downwash due to the increased camber line curvature, yet being able to avoid stall by means of a thicker and rounder leading edge.Thus, the DyAS analysis at different time steps shows that as the impulsively started airfoil moves from an acceleration phase to a steady state regime solution, the shape modifications leading to increased lift transit from a purely symmetric increase of the thickness in the mid-chord region, to a non-symmetric modification of the camber line united with a symmetric leading edge thickness increase, respectively.Such behavior is indicated by the sign of c 4 coefficient in the eigenvector characterizing the one dimensional active subspace, which is likely detecting that at steady state, regime solution, airfoils with higher camber line curvature and thicker leading edges produced higher downwash.
We underline that the eigenvector components of all the time instants presented corresponding to the coefficients c 1 , c 5 , d 1 , and d 5 are almost zero.This means that on average the lift coefficient is almost flat along these directions.We are going to exploit this fact by freezing these parameters and constructing a GPR on a reduced parameter space.

GPR approximation and prediction of the lift coefficient
The previous analysis pointed out the presence of several input parameters with minimal average influence on the target function.Making use of such consideration we construct a response surface which only depends on the remaining parameters.Both for the full parameter space and the reduced one, we use a Gaussian process regression with a RBF kernel implemented in the open source The relative error is computed on 100 test samples, using the high-fidelity lift coefficient to train the regression for t ≤ 20 s, while for t > 20 s the DMD forecasted states are used for the training.
Python package GPy [19].We then compare the performance of the two regression strategies by computing the relative error over a test data set composed by 100 samples.The error is computed as the Euclidean norm of the difference between the exact and the approximated solution over the norm of the exact solution.The training set is composed by the same 70 samples, in 10 dimensions for the GPR over the original parameter spaces, and in 6 dimensions for the reduced one.Up to t = 20 s the training is done using the high-fidelity simulations.To speed up the convergence to the regime state (t = 30 s) we applied the DMD to get the future-state prediction of the lift.In Figure 8 we compare the two GPR performance at each of the time steps analyzed in the simulations.Until 12 s, the regressions behave in a very similar fashion, while from 15 s the accuracy gain obtained by distributing the 70 samples in a lower dimensional space becomes significant.The error gap between the 6 and 10 dimensional response surface in fact, consistently increases from 1% at 15 s to more than 4% at steady state.

Conclusions and perspectives
We presented a computational pipeline to improve the approximation of the time-varying lift coefficient of a parametrized NACA airfoil.The pipeline comprises automatic mesh deformation through RBF interpolation, high-fidelity simulation with finite volume method of turbulent flow past the airfoil, global sensitivity analysis exploiting AS, and future state prediction via DMD reduced order method.This resulted in more accurate Gaussian process regression of the lift coefficient even if in a reduced parameter space.
After the creation of the high-fidelity solutions database the application of AS highlighted a possible reduction of the parameter space due to negligible contributions of 4 different parameters.We exploit this reduction to construct a GPR over a smaller parameter space, thus improving its performance.Since the training of the regression model is done over 6 dimension instead of 10, given the same high-fidelity database dimension, the GPR is able to better approximate the solution manifold.This results in better lift coefficient predictions for new untried parameters.We also applied DMD to have future-state prediction of the target function up to 30 seconds and proved that the effective gain of the new GPR is preserved also for any time after the 20 seconds simulated with FV.In particular from 13 seconds the actual gain is significant, at 15 seconds we have an increased performance of 1% in the relative error.Evolving in the future the error drop increases up to more than 4% at regime.This computational pipeline can be seen as a parametric dynamic mode decomposition for some extent.Moreover, the sensitivity analysis has a negligible computational cost with respect to the creation of the offline high-fidelity database.
Future developments can be the study of adaptive sampling strategies exploiting a generic n-dimensional active subspace, and the coupling of different model order reduction methods.It would be interesting to use this non-intrusive setting as a preprocessing tool to reduce the number of simulations required to build a reduced basis space which is later used in an intrusive manner [46].We think this new computational pipeline can be of much interest in the context of shape optimization and dynamical systems.

Figure 1 :
Figure 1: Airfoil shape functions with respect to the profile abscissa.The leading edge corresponds to x = 0.

Figure 2 :
Figure 2: Sketch of the computational domain used to solve the fluid dynamics problem in its reference configuration.The left picture reports a schematic view on the domain with the main geometrical dimensions.The right plot reports a zoom on the mesh in the proximity of the wing.

Figure 3 :
Figure 3: The left picture reports in light blue the envelope of all the tested configurations used during the training stage.The right picture depicts the flow velocity streamlines for one particular sample inside the training set µ = [0.0071;0.0229; 0.0015; 0.0015; 0.0087; 0.0107; 0.0033; 0.0130; 0.0247; 0.0280].

Figure 4 :
Figure 4: On the left the sufficiency summary plot for the lift coefficient at time t = 6.0 seconds.On the right the first eigenvector components at the corresponding parameters.

Figure 5 :
Figure 5: On the left the sufficiency summary plot for the lift coefficient at time t = 10.0 seconds.On the right the first eigenvector components at the corresponding parameters.

Figure 6 :
Figure 6: On the left the sufficiency summary plot for the lift coefficient at time t = 14.0 seconds.On the right the first eigenvector components at the corresponding parameters.

Figure 7 :
Figure 7: On the left the sufficiency summary plot for the lift coefficient at time t = 18.0 seconds.On the right the first eigenvector components at the corresponding parameters.

Figure 8 :
Figure 8: The relative error of the approximated outputs at different times.The relative error is computed on 100 test samples, using the high-fidelity lift coefficient to train the regression for t ≤ 20 s, while for t > 20 s the DMD forecasted states are used for the training.