Kristofer Jovanov
email: k.jovanov@tudelft.nl

Roeland De Breuker
email: r.debreuker@tudelft.nl

Mostafa M Abdalla
email: m.abdalla@tudelft.nl

Christophe Blondeau
email: christophe.blondeau@onera.fr

A linear aerodynamics-based preconditioner for high-fidelity aeroelastic analysis and sensitivity analysis

scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. Introduction

The trend of employing high-fidelity models in multidisciplinary design optimization (MDO) has become more pronounced with increasing computational resources, the advent of code parallelism and improving solution algorithms. [START_REF] Maute | Coupled Analytical Sensitivity Analysis and Optimization of Three-Dimensional Nonlinear Aeroelastic Systems[END_REF][START_REF] Barcelos | A Schur-Newton-Krylov solver for steady state aeroelastic analysis and design sensitivity analysis[END_REF][START_REF] Kennedy | Parallel Solution Methods for Aerostructural Analysis and Design Optimization[END_REF][START_REF] Kenway | Scalable Parallel Approach for High-Fidelity Steady-State Aeroelastic Analysis and Adjoint Derivative Computations[END_REF][START_REF] Blondeau | Recent Achievements Towards Aero-Structure Gradient Computation using High-Fidelity CFD-CSM in the ONERA elsA Software[END_REF] State-of-the art solvers tailored for single-disciplinary problems are commonly coupled through well-defined interfaces to facilitate the solution of multidisciplinary analysis problems. If subject to gradient-based optimization, additional design sensitivity analysis is required to compute the gradients of the responses with respect to any design variables. In the context of aerostructural design, these responses typically vary from structural weight and stresses to aerodynamic lift and drag. These quantities, and their associated gradients, are often regarded as inputs for the "black box" optimizer. The bottleneck in highfidelity MDO is the high computational cost for acquiring these inputs. Much effort is therefore invested in improving the analysis and sensitivity analysis of high-fidelity fluid-structure interation (FSI).

Maute et al. [START_REF] Maute | Coupled Analytical Sensitivity Analysis and Optimization of Three-Dimensional Nonlinear Aeroelastic Systems[END_REF] solve the quasi-static aeroelastic problem using a nonlinear block Gauss-Seidel (NLBGS) method augmented by a stationary under-relaxation technique. This method, which is fundamentally a loosely coupled approach, couples the nonlinear Euler flow solver to a linear finite element model through a spring analogy mesh deformation strategy. The sensitivity problem is solved by applying the same method, modified to a linear problem. This straightforward approach offers a high level of software modularity.

As highlighted in a follow-up paper by Barcelos et al., [START_REF] Barcelos | A Schur-Newton-Krylov solver for steady state aeroelastic analysis and design sensitivity analysis[END_REF] the NLBGS method lacks robustness and efficiency. Moreover, the convergence behaviour is highly influenced by the relaxation parameter, whose optimal value is in general unknown. To address this, they propose the Schur-Newton-Krylov (SNK) solver. This approach is based on a Schur complement formulation applied to the interface variables of the aerodynamic mesh. The nonlinear quasi-static FSI problem is then solved by a Newton-Krylov method were the exact Jacobian is used. The same Krylov solver is applied for the subsequent linear sensitivity problem. It is shown that this method is more robust and efficient than the previous method described by Maute. Kenway et al. [START_REF] Kenway | Scalable Parallel Approach for High-Fidelity Steady-State Aeroelastic Analysis and Adjoint Derivative Computations[END_REF] solve the steady-state aeroelastic problem by a matrix-free coupled Newton-Krylov (CNK) method. The aeroelastic Jacobian is approximated by finite differences and the successive linear subproblems are solved by FGMRES. [START_REF] Saad | A Flexible Inner-Outer Preconditioned GMRES Algorithm[END_REF] This approach shows reduced computation times compared to the conventional NLBGS approach, which in this work is augmented by a dynamic Aitken relaxation strategy. The linear sensitivity problem, in contrast to the nonlinear static problem, is solved by the exact aeroelastic Jacobian. The coupled linear system is solved by GMRES, [START_REF] Saad | GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF] where the aerodynamic block is preconditioned by ILUfactorization and the structural block is preconditioned by an exact factorization. No preconditioning is applied to the off-diagonal coupling terms. They show that this method converges faster compared to the loosely coupled approach. However, the main shortcoming of an ILU-preconditioner for the aerodynamic block is that the exact flux Jacobian need to be explicitly stored. This can result in vast memory requirements for large-scale problems. [START_REF] Nielsen | Recent improvements in aerodynamic design optimization on unstructured meshes[END_REF] Despite the advancements described above, high-fidelity numerical methods still remain computationally expensive. In particular transonic flow simulations where the Euler or the RANS equations are employed are mainly regarded as computationally expensive. Consequently, these higher order aerodynamic methods are confined to the later stages of the design process where they are used for verification purposes. [START_REF] Jameson | Re-engineering the design process through computation[END_REF] The early design stages, in contrast to the later ones, often rely on panel methods that deliver accurate results at a low computational cost. Panel methods, however, are not equally efficient in modelling nonlinear flow phenomena, e.g. shocks and flow separation. This implies that feasible designs might be filtered out in the early design stages due to the seemingly inferior transonic flow modelling capabilities of panel methods. It is therefore imperative that high-fidelity models are adopted early in the design process. This is only viable if these models are computationally inexpensive.

In an attempt to rectify the aerodynamic loads (aeroloads) obtained by panel methods and still maintain a low computational cost, several correction strategies have been developed. [START_REF] Giesing | Correction Factor Techniques for Improving Aerodynamic Prediction Methods[END_REF][START_REF] Palacios | Assessment of strategies for correcting linear unsteady aerodynamics using CFD or experimental results[END_REF][START_REF] Jovanov | Accelerated convergence of static aeroelasticity using low-fidelity aerodynamics[END_REF][START_REF] Jovanov | Accelerated convergence of high-fidelity aeroelasticity using low-fidelity aerodynamics[END_REF] The basic principle is to correct the aerodynamic influence coefficient (AIC) matrix or the aeroloads based on data obtained by higher order computational fluid dynamics (CFD) simulations or wind tunnel tests. These correction methods have proven fruitful and are actively used in industry. However, correction methods for steady-state flows are mostly applied to the equilibrium equations. Correcting sensitivities for gradient-based optimization has not received equal amount of attention.

In this paper, we propose the application of a panel method as a preconditioner for solving the highfidelity static aeroelastic analysis and sensitivity analysis problem. Initially, the high-fidelity equations are expressed in Schur complement form applied to the structural variables. A quasi-Newton method is then used to solve the nonlinear static aeroelastic problem, where the Newton correction is approximately solved by a fixed-point strategy. The fixed-point problem is preconditioned by a low-fidelity Schur complement obtained by the inexpensive panel method. The same principle of solving the fixed-point problem is applied to the subsequent linear sensitivity problem. The main difference is that the linear system in the sensitivity analysis problem is converged to a low tolerance to ensure accurate gradients for the optimizer.

Our proposed method is similar to the work by Heil. 14 However, our approximated Schur complement requires information from the structural model and is therefore interdisciplinary. Furthermore, in this paper we investigate transonic flow problems of moderate size by employing the Euler equations for a 1.11 million cell model. The structural model is linear finite element model with 1812 degrees of freedom. The focus in this paper is purely on analysis. Its application in optimization is envisioned in a follow-up paper.

The remainder of this paper is organized as follows: the computational modules and routines necessary for the aeroelastic analysis and sensitivity analysis are described in Section II. Three different aeroelastic analysis methods, of which one is our proposed method, are described in Section III. The sensitivity analysis with our proposed approach is presented Section IV, followed by the comparison results in Section V. The paper is concluded by a brief discussion and reflection on possible improvements and suggestions for future work in Section VI.

II. Computational modules

Several computational modules must be established to enable aeroelastic analysis and sensitivity analysis, with each module performing individual subtasks and being capable of exchanging information in a systematic fashion. This section provides an overview of the modules and resources necessary for the aeroelastic analysis and the subsequent sensitivity analysis.

Organizing the work in well-defined modules might increase the overall computational effort as opposed to monolithic codes operating in single-language environments. However, this approach allows for a higher level of software modularity which is a coveted feature in industry. It enables engineers to choose among well-established solvers with minimal effort dedicated to adapting the existing code.

In the following subsections the (A) aerodynamic solvers, the (B) structural solver, the (C) coupling module and the (D) low-fidelity Schur complement will be described, respectively.

A. Aerodynamic Solvers

Two aerodynamic solvers are necessary in this work: one solving a set of low-fidelity (LoFi) equations and another solving a set of high-fidelity (HiFi) equations. The former falls into the category of panel methods where vortex ring elements are used to approximate the aerodynamic pressure. The flat two-dimensional panels are distributed on the camber surface of the wing and not on the actual surface (see Figure 1a). It should be mentioned that this particular panel distribution can generate a small offset in the results for thick wing configurations. The discrete linearized potential flow equations can be expressed in residual form as

R a (γ, x p ) = Aγ -b = 0 (1)
where γ is the unknown vector of vortex strengths, A is the dense AIC matrix and b is the Neumann boundary condition vector. Both A and b are functions of the panel coordinates x p . In order to avoid confusion with the HiFi solver later in the text, terms associated with the LoFi solver will be marked with a tilde notation. The AIC matrix is purely a function of geometry and can be assembled through the application of the Biot-Savart law. Once the linear system (1) is solved, the aeroloads can be obtained in a post-treatment step by the Kutta-Joukowski theorem

f a = ρ ∞ V ∞ × γ (2)
where ρ ∞ and V ∞ are the free-stream density and free-stream velocity vector, respectively. A Prandtl-Glauert coordinate transformation rule is applied to account for compressibility effects in high-subsonic flows. The computational routines for the LoFi solver are in accordance with, and can be further studied in, Katz and Plotkin. [START_REF] Katz | Low-Speed Aerodynamics[END_REF] In addition to aeroloads, the LoFi solver must prove capable of generating the LoFi Schur complement, S. The formulation of this term will be described in detail in Subsection D.

The main shortcoming of the LoFi solver is that its application is limited to attached flow conditions in low-subsonic to moderately high-subsonic flows. Panel methods, in general, lack the predictive capability of modelling nonlinearities in transonic flows. To address this shortcoming, a more advanced aerodynamic solver is required.

The HiFi aerodynamic solver is the multi-block structured flow solver elsA. [START_REF] Cambier | An Overview of the Multi-Purpose elsA Flow Solver[END_REF] In this work we solve the compressible Euler equations to account for recompression shocks. While the Euler equations might not be acknowledged as "high-fidelity" by many researchers, there are clear differences between the Euler equations and the linearized potential flow equations. It is based upon these relative differences that we have permitted to label the Euler equations in this work as high-fidelity. The discrete steady-state equations of the Euler flow problem can be expressed in residual form as

R a (w, x) = 0 ( 3
)
where w is the unknown vector of conservative variables, with five unknowns per cell, and x are the mesh coordinates. The aerodynamic residual, or the flux balance, is cancelled out when the solution converges to a steady-state equilibrium. elsA solves the nonlinear steady-state problem by an implicit backward-Euler (BWE) time integration scheme. In addition to solving the steady-state problem for the static aeroelastic analysis, a linearization of these equations are necessary for the subsequent sensitivity analysis. Consequently, four additional components are necessary from the flow solver: 1. ∂R a /∂w: This is the exact linearization of the discrete residual (3) with respect to the conservative variables, also known as the exact flux Jacobian. If the underlying spatial discretization scheme for the steady-state equations is of second order, then the flux Jacobian is a 13-banded sparse matrix formulated by an extended neighbours-of-neighbours stencil.

2. ∂ R a /∂w: This is an accurate approximation of the Jacobian matrix used for the implicitation of the convective fluxes in the BWE scheme. Usually a standard first order upwind approximation of the differential of the flux balance is preferred for better numerical conditioning. The approximate Jacobian matrix can then be expressed in the factorized form: ∂ R a /∂w = (L+D+U ). The resolution of system (4) relies on a method of relaxation which consists in calculating successive approximate solutions of the exact system using forward and backward sweeps through the computational domain. [START_REF] Peter | Large Stencil Viscous Flux Linearization for the Simulation of 3D Compressible Turbulent Flows with Backward-Euler Schemes[END_REF] 3. ∂R a /∂x: This is referred to as the geometric flux Jacobian. Similar to the exact flux Jacobian, the geometric flux Jacobian maintains a sparse structure. It describes how the flux balance changes with respect to aerodynamic grid perturbations.

4. ∂f a /∂w: This term can be obtained by linearizing the surface load pressure integration scheme. Only cells adjacent to the surface affect this term and as such the term can be decomposed to (∂f a /∂w b ) × (∂w b /∂w), where w b are the wall conservative variables that meet the physical boundary conditions (typically the wallslip condition for Euler equations). They are adequately extrapolated from the numerical scheme variables w in the adjacent cells.

All the components described above are only formulated as matrix-vector products in order to save memory. This feature is desirable for iterative solvers, such as GMRES, which do not require an explicit matrix formulation but can handle matrix-vector products to solve large sparse linear systems. Moreover, all the components in elsA are computed analytically by hand-differentiating the code.

To obtain the matrix-vector products, modifications to the existing elsA source code were necessary. The opt-module in elsA, dedicated to aerodynamic shape optimization, implements a fixed-point iterative scheme for the aerodynamic sensitivity problem. An extensive overview of the gradient computation capabilities is covered in Blondeau et al. [START_REF] Blondeau | Recent Achievements Towards Aero-Structure Gradient Computation using High-Fidelity CFD-CSM in the ONERA elsA Software[END_REF] The preconditioned fixed-point iterative scheme for the direct method gradient computation can be formulated as

V ∆t I + ∂ R a ∂w Implicit Operator dw dp (n+1) - dw dp (n) Correction = - ∂R a ∂x dx dp - ∂R a ∂w dw dp (n) Right -Hand Side (4)
where a pseudo-time term ∆t in the implicit operator is introduced to stabilize the convergence behaviour by increasing the diagonal dominance of ill-conditioned systems. V are the cell volumes and dx/dp is a geometric gradient term that can be determined by the coupling module in Subsection C. It is apparent that the unknown vector, dw/dp, converges when the correction, and hence the right-hand side (RHS), approaches zero.

The modification sequence in the elsA/Opt module can be narrowed down to three steps:

1. First, the computation of (∂R a /∂x) × (dx/dp) is dropped in the RHS of ( 4). This operation is only required once per gradient computation. Components 3 and 4, which are required for the aeroelastic sensitivity analysis, are obtained in separate routines.

2. The RHS at this point is zero. As it is a restart procedure, the term dw/dp is initialized to an external vector of our choice. The RHS construction routine is then called which results in the product of the exact flux Jacobian matrix ∂R a /∂w by the prescribed vector.

3. At this stage the RHS is equivalent to the flux Jacobian matrix-vector product. We then reset the RHS to a new input vector, v, and (4) now reduces to the following linear system:

V ∆t I + ∂ R a ∂w y = v ( 5
)
where y is the matrix-vector product of the inverse flux Jacobian approximation. Hence, by relaxing the linear system (5) with multiple sweeps the inverse of the flux Jacobian approximation is multiplied by v. Typically 4-8 cycles, resulting in 8-16 sweeps, are sufficient. An increased number of cycles improves the quality of the preconditioner. However, this comes at the expense of an increased computation time. As an additional control parameter on the diagonal dominance of the preconditioner, it is the users choice to adjust or even cancel the pseudo-time term in the implicit operator.

With the new capabilities at our disposal, the aerodynamic sensitivity problem (4) can alternatively be solved by an external preconditioned GMRES algorithm. The right-hand side is thus solved directly in its original form

∂R a ∂w dw dp = - ∂R a ∂x dx dp (6)
which is preconditioned by the approximate flux Jacobian. A simple test case is conducted to demonstrate the difference in convergence behaviour between the external GMRES solver and the standard fixed-point LU-SSOR solver. Our test case is identical to the gradient computation test case in Blondeau et al. [START_REF] Blondeau | Recent Achievements Towards Aero-Structure Gradient Computation using High-Fidelity CFD-CSM in the ONERA elsA Software[END_REF] The computational Euler model is a 1.11 million cell mesh divided into five blocks (see Figure 1b). An upwind Roe scheme with a MUSCL interpolation associated to a Van Albada limiting function is applied. The Mach number is 0.734 and the incidence angle is 2.8 • . The freestream flow condition is set to ρ ∞ = 1.225m 3 /kg and p ∞ = 101325N/m 2 . The design variable used to construct the right-hand side in Equation ( 6) controls the section camber at 50% span and affects linearly the region from root to 90% of wing span.

Initially, the steady-state problem is converged to machine precision before the gradient computation is initialized. We implement a preconditioned GMRES solver with a restart subspace of 20 Krylov vectors. This corresponds to 0.82 Gigabytes in our case which is a reasonable memory storage. The number of cycles for the preconditioner, in both methods, is fixed at 8. Figure 2 shows a clear advantage for the GMRES solver which has a convergence rate roughly six times faster compared to the LU-SSOR method. The computation times for the two methods follows the same trend as in Figure 2. The main additional computation time that emanates from the GMRES algorithm is the formulation of the Kylov basis by the Arnoldi iteration which is very small compared to the matrix-vector products. Hence, by implementing non-intrusive modifications to the current source code of elsA we have shown the benefits of implementing a matrix-free GMRES solver. Moreover, the LU-SSOR preconditioner does not require the explicit storage of the flux Jacobian as opposed to ILU-preconditioned GMRES solvers. [START_REF] Luo | A Fast, Matrix-free Implicit Method for Compressible Flows on Unstructured Grids[END_REF] This implies a significant advantage for large-scale CFD problems. [START_REF] Nielsen | Recent improvements in aerodynamic design optimization on unstructured meshes[END_REF] However, it should be mentioned that the factorization error of LU-SSOR is larger than ILU. An ILU-preconditioner is therefore anticipated to yield higher converge rates than LU-SSOR.

B. Structural Solver

The structural solver in this work is the finite element analysis program Nastran. [START_REF] Macneal | The Nastran theoretical manual[END_REF] The main task of Nastran is to assemble and output the structural stiffness matrix, K, to be used externally for the aeroelastic analysis and the sensitivity analysis. If the deformations are of moderate size the structure can be assumed to exhibit a linear-elastic behaviour. Hence, the resulting stiffness matrix need only be generated once for the initial finite element model and stored in memory. In the event of large structural deformations the assumption of linear elasticity no longer holds and the stiffness matrix need to be continuously updated. In this paper only linear elasticity is considered. The discrete structural equations can then be expressed as

R s (u) = Ku -f s = 0 (7)
where u are the structural degrees of freedom and f s are any externally applied loads. In addition to the stiffness matrix, pseudo-loads are required for the subsequent sensitivity analysis. The pseudo-loads can be derived by linearizing (7) with respect to structural design variables p:

dR s dp = K du dp + ∂K ∂p u - df s dp P seudo-Loads = 0. ( 8
)
This equation is derived, and its components are explained in detail, in the Nastran Design Sensitivity and Optimization User's Guide. [START_REF] Raymond | MSC.Nastran 2012, Design Sensitivity and Optimization Users Guide[END_REF] What is of interest to us, prior to the sensitivity analysis, is the explicit dependence of the structural residual with respect to structural design variables. The only term that satisfies this requirement in Equation ( 8) is the stiffness gradient

∂R s ∂p = - ∂K ∂p u * (9)
where u * represents the structural displacements at static aeroelastic equilibrium. The remaining terms in (8) are implicit functions of the residual and are only available when the sensitivity analysis has converged. The gradient term in Equation ( 9) is the right-hand side for the sensitivity analysis in Section IV and is of size n f dof × n p , where n f dof is the number of unconstrained structural degrees of freedom and n p is the number of structural design variables.

The stiffness matrix and the pseudo-loads can be readily accessed by adding Direct Matrix Abstraction Program (DMAP) [START_REF] Raymond | MSC.Nastran 2013, DMAP Programmer's Guide[END_REF] alters in the case control section of the bulk data files. The stiffness matrix requires a SOL 101 solution sequence whereas the pseudo-loads require a SOL 200 static analysis solution sequence.

The DMAP alters can then export the matrices during the execution of Nastran. The exported op4 files are converted to convenient sparse matrix triplet format using external fortran scripts to be used in the aeroelastic analysis and sensitivity analysis.

C. Coupling module

The coupling module enables the deformation of the aerodynamic mesh and the load interpolation on the structural nodes. As the aerodynamic models and the structural model in this work are non-conforming, there is a necessity for a robust and efficient mesh deformation algorithm. To this end, we apply an algebraic coupling routine based on radial basis functions. The basic principle is to formulate a coupling matrix H that interpolates the displacements of the aerodynamic grid based on the structural deformations:

u a = Hu ( 10
)
where u a are the aerodynamic grid displacements. The coupling matrix in this work is derived in accordance with the methodology in Beckert and Wendland [START_REF] Beckert | Multivariate interpolation for fluid-structure-interaction problems using radial basis functions[END_REF] for multivariate data interpolation. We use the globally supported thin-plate spline function for the HiFi aerodynamic model which results in a dense interpolation matrix. To circumvent the need for large amounts of memory, the coupling matrix is only formulated as a matrix-vector product. Hence, it is never stored explicitly. Moreover, the mesh deformation code is executed on multiple processors. As the radial basis function approach does not require connectivity information, parallelizing the code becomes a trivial task.

A simple test case was conducted to assess the parallel scalability. Mesh deformation is performed on the 1.11 Euler model described in Subsection A. These computations, including all others in this work, were performed on Intel Xeon E5-2660 v3 Processors with a 2.60GHz CPU frequency. The mesh deformation wall time was clocked for 2, 4, 8, 16 and 32 processors. The results in Figure 3a show very good scalability up to approximately 16 processors. Hence, the CPU time for one mesh deformation on a 1.11 million cell model can be reduced from 86 seconds on 2 processors to 18 seconds on 16 processors. Another method, in addition to parallelism, of improving the computation time is to limit the number of surface points in the formulation of the coupling matrix. There are numerous data point reduction algorithms that can be used to reduce the size of the coupling matrix without significantly affecting the quality of the deformation. [START_REF] Rendall | Efficient mesh motion using radial basis functions with data reduction algorithms[END_REF][START_REF] Rendall | Reduced surface point selection options for efficient mesh deformation using radial basis functions[END_REF] In this work, due to the structured mesh, we simply apply every 10th surface point to formulate the coupling matrix (see Figure 3b). This has been found to yield a favourable trade-off between mesh deformation accuracy and computation time.

In previous work [START_REF] Jovanov | Accelerated convergence of static aeroelasticity using low-fidelity aerodynamics[END_REF][START_REF] Jovanov | Accelerated convergence of high-fidelity aeroelasticity using low-fidelity aerodynamics[END_REF] we used the principle of equivalence of work which allows for the transposed coupling matrix to be used as a load transfer operator

f s = H T f a . (11)
This, however, is not always a favourable option. Based on our experience this operator can result in unrealistic oscillatory loads in regions where structural nodes are tightly clustered. We therefore resort to a simple, yet highly effective and robust, nearest-neighbour scheme for multivariate interpolation. [START_REF] Thvenaz | Interpolation Revisited[END_REF] In this work, the load transfer operator is denoted by T , such that

f s = T f a . (12)

D. Low-fidelity Schur complement

As emphasized in Subsection A, the LoFi Schur complement needs to be readily accessible by the panel method. Only its derivation is discussed in this subsection. Its purpose will be covered in Section IV. The LoFi Schur complement can be derived by formulating a static LoFi aeroelastic problem. Consequently, the aerodynamic and structural residual in Equations ( 1) and ( 7), respectively, can be concatenated and expressed as an aeroelastic residual

R ael = R a (γ, u) R s (u, γ) = 0. ( 13
)
If the aeroelastic residual is linearized with respect to the unknown variables the following linear system is obtained:

∂ Ra ∂γ ∂ Ra ∂u ∂Rs ∂γ ∂Rs ∂u ∆γ ∆u = - R a R s (14)
where the incremental vector of vortex strengths and structural displacements is also known as the Newton correction. This is in fact results in Newton's method for solving the static aeroelastic problem where the tangent operator is the aeroelastic Jacobian. The diagonal blocks are inherently single-disciplinary and are therefore determined individually by their respective computational solvers. The off-diagonal blocks include coupling matrices that require information from both computational solvers as well as the coupling module. The Schur complement of the upper-diagonal block in the tangent operator is:

S =   ∂R s ∂u - ∂R s ∂γ ∂ R a ∂γ -1 ∂ R a ∂u   (15)
The four blocks required to construct the LoFi Shur complement are studied in detail:

1. ∂R s /∂u is nothing but the reduced stiffness matrix K which is already provided as an output by the structural solver, as described in Subsection B.

2. ∂R s /∂γ describes how the structural residual change with respect to vortex strength perturbations. This matrix can be decomposed as:

(∂R s /∂ f a )×(∂ f a /∂γ).
The first term is the load transfer operator T that interpolates loads from the aerodynamic panels on the structural nodes. The second term can be easily obtained by differentiating the Kutta-Joukowski theorem (2).

3. ∂ R a /∂γ is the AIC matrix, A. This can easily be assembled by the LoFi solver, inverted and stored in memory.

4. ∂ R a /∂u can be decomposed as: (∂ R a /∂x p ) × (∂x p /∂u). The first term requires a linearization of the aerodynamic residual (1) with respect to the panel coordinates. Both the AIC matrix and the Neumann boundary condition vector are affected by the panel coordinates. This term, in particular the AIC derivatives, is cumbersome to obtain and can result in excessive computation times for aerodynamic models with fine panel grids. To overcome this shortcoming, a parallel approach with analytically obtained gradients is implemented. The second term in the decomposition is a displacement interpolation operator H which is readily available by the coupling module. For the panel method an infinite-plate spline [START_REF] Harder | Interpolation using surface splines[END_REF] method is used to conform the panels with the structural displacements.

III. Aeroelastic Analysis

The main purpose of the static aeroelastic analysis is to determine the state of equilibrium for the FSI problem. This part is crucial in design optimization as it allows for important quantities such as aerodynamic drag and lift, and structural stresses to be evaluated. All these responses can be represented by J = J (x, w, u; p) (16) where J can be any response. A response is a collective term for both the objective function and other possible constraint functions. The design variables can impose an explicit influence on the response, as well as an implicit influence through the solution of the underlying governing equations. To obtain these responses for a given set of design variables the static aeroelastic problem need to be solved. The static HiFi aeroelastic problem can be defined in a two-field formulation as

R ael = R a (w, u; p) R s (u, w; p) = 0. ( 17
)
In contrast to the static LoFi problem (13), we now want to satisfy a state of equilibrium between the Euler flow equations and the structural equations. The static aeroelastic problem is solved when the flow and structural subproblem converge within tolerances ǫ a1 and ǫ s1 , respectively:

R (n+1) a 2 ≤ ǫ a1 R (0) a 2 (18)
R (n+1) s 2 ≤ ǫ s1 R (0) s 2 . ( 19
)
For any iteration with index n, both criteria must be met simultaneously before the problem is considered converged. There are numerous possibilities for solving the nonlinear problem (17). [START_REF] Maute | Coupled Analytical Sensitivity Analysis and Optimization of Three-Dimensional Nonlinear Aeroelastic Systems[END_REF][START_REF] Barcelos | A Schur-Newton-Krylov solver for steady state aeroelastic analysis and design sensitivity analysis[END_REF][START_REF] Kennedy | Parallel Solution Methods for Aerostructural Analysis and Design Optimization[END_REF][START_REF] Kenway | Scalable Parallel Approach for High-Fidelity Steady-State Aeroelastic Analysis and Adjoint Derivative Computations[END_REF] In this work three solution methods are highlighted in the following subsections.

Nonlinear Block Gauss-Seidel (NLBGS) method

The first method is a conventional NLBGS method. It is initiated by approximately solving the steadystate flow problem (3) by performing n a flow subiterations for the initial shape. The aeroloads are then interpolated on the finite element model (12) and the linear structural equations (7) are solved. Finally, the aerodynamic mesh is updated (10) to align with the deformed structure whereby the routine is restarted and iterated until the aeroelastic problem converges. This method is enhanced by a dynamic relaxation parameter, θ, obtained by Aitken acceleration. [START_REF] Irons | A Version of the Aitken Accelerator for Computer Iteration[END_REF] The NLBGS method is described in Table 1.

Table 1: The NLBGS algorithm (0) Set initial values for u (0) , w (0) , x (0) and θ (0) . Iterate n = 0, ..., ∞

(1) Approximately solve the flow problem, R a (w (n+1) , x (n) ) ≈ 0, by n a flow subiterations.

(2) Interpolate the aeroloads on the structural model: f

(n+1) s = T f (n+1) a .
(3) Solve the structural subproblem: R s (u (n+1) ) = 0.

(4) Evaluate the increment ∆u (n+1) = u (n+1)u (n) and update the relaxation parameter by Aitken acceleration θ

(n+1) = θ (n) 1 - (∆u (n+1) -∆u (n) ) T ∆u (n+1) (∆u (n+1) -∆u (n) ) T (∆u (n+1) -∆u (n)
) .

(5) Relax the structural solution:

u (n+1) = 1 -θ (n+1) u (n) + θ (n+1) u (n+1) .
(6) Update the aerodynamic mesh by: x (n+1) = x (0) + Hu (n+1) .

Schur Quasi-Newton (SQN) method

The second method is a quasi-Newton algorithm. The term "quasi" implies that some form of approximation is applied to the Jacobian. The method can be derived by linearizing the aeroelastic residual (17) with respect to the unknown variables w and u, respectively ∂Ra ∂w ∂Ra ∂u ∂Rs ∂w ∂Rs ∂u

∆w ∆u = - R a R s (20)
where the variable increments represent the Newton correction which is used to update the variables for the subsequent iteration

w (n+1) u (n+1) = w (n) u (n) + ∆w ∆u . (21)
The design variables are constant throughout the aeroelastic analysis and are therefore not considered in the linearization. A Schur complement formulation applied to the structural variables allows for the elimination of the aerodynamic residual

∂R s ∂u - ∂R s ∂w ∂R a ∂w -1 ∂R a ∂u S ∆u = -R s ( 22
)
where S is the HiFi Schur complement. This method, which is an exact Newton method, exhibits a quadratic rate of convergence in the proximity of the solution. Unfortunately, for large-scale problems the HiFi Schur complement can be difficult and computationally expensive to determine. In addition, matrix inversion or a complete factorization of the flux Jacobian becomes highly impractical. To circumvent this, the HiFi Schur complement can be replaced by its computationally inexpensive LoFi counterpart (15):

S∆u = -R s . (23)
Consequently, the Schur complement becomes an approximation, hence the name: quasi-Newton. Since it is expressed in Schur complement form, we refer to this method as the Schur quasi-Newton (SQN) method. This method can become a powerful tool if the approximation is close to the exact Schur complement and it is computationally cheap to obtain.

A similar quasi-Newton approach is performed by Gerbeau and Vidrascu 28 for transient FSI problems in blood flows. They propose to neglect the expensive mesh motion and the nonlinear terms in the fluid model. In our work, however, the LoFi solver is a stand-alone module that is completely independent from the HiFi solver. The routines for the PQN solver are described in Table 2. Note that in this method the structural residual is only evaluated and not solved, in contrast to the NLBGS method.

Table 2: The SQN algorithm (0) Set initial values for u (0) , w (0) and x (0) . Iterate n = 0, ..., ∞

(1) Approximately solve the flow problem, R a (w (n+1) , x (n) ) ≈ 0, by n a flow subiterations.

(2) Interpolate the aeroloads on the structural model: f

(n+1) s = T f (n+1) a .
(3) Evaluate the structural residual: R s = Ku (n)f (n+1) s .

(4) Solve the approximated Newton correction: S∆u = -R s .

(5) Update the structural solution by the correction: u (n+1) = u (n) + ∆u. ( 6) Update the aerodynamic mesh by: x (n+1) = x (0) + Hu (n+1) .

Preconditioned Schur Quasi-Newton (PSQN) method

The third method, which is our proposed method, capitalizes on the elsA modifications described in Subsection A to solve (22) by reformulating it into a fixed-point problem

S ∆u (k+1) -∆u (k) = -S∆u (k) + R s . ( 24
)
The fixed-point formulation ( 24) yields a subiteration sequence with index k. The linear system is solved by a direct method whereas the Schur matrix-vector product, S∆u (k) , is approximated by GMRES. Based on this information two questions emerge: how much do we need to converge the correction and how much do we need to converge the Schur matrix-vector product? An accurate correction is anticipated to result in fewer Newton iterations. This might seem advantageous at first. However, to achieve the quadratic convergence rate of Newton's method, in addition to being close to the solution, the correction must be solved exactly. This implies that the flow problem must converge to a low tolerance by performing sufficient amount of subiterations, n a . This is not a viable option as the overall computation time needed to perform an exact Newton method would vastly surpass that of an approximated Newton method.

In our approach only one fixed-point iteration is performed for the Newton correction and the residual norm in GMRES is converged to ǫ g = 1e -1 . To distinguish from the exact Schur matrix-vector product, this approximated version can be denoted as S∆u (k) , such that Equation ( 24) can be formulated as

S ∆u (k+1) -∆u (k) = -S∆u (k) + R s . ( 25
)
It should be noted that S is based on an exact linearization of the Euler flow model whereas S is based on the panel method. Since the aeroelastic problem is solved in Schur complement form and S is used as a preconditioner to approximate the correction, this method is referred to as a preconditioned Schur quasi-Newton (PSQN). The details of the algorithm are explained in Table 3. This method resembles the SQN. The main difference is that the Newton correction is enhanced by the additional steps 5-7. Note that the structural stiffness matrix, the coupling operators and the elsA components described in Subsection A, are used to clarify the construction of S∆u.

Table 3: The PSQN algorithm (0) Set initial values for u (0) , w (0) and x (0) . Iterate n = 0, ..., ∞

(1) Approximately solve the flow problem, R a (w (n+1) , x (n) ) ≈ 0, by n a flow subiterations.

(2) Interpolate the aeroloads on the structural model: f ,n) . ( 7) Determine the fixed-point correction, ∆u (k+1,n) , by solving Equation (24).

(n+1) s = T f (n+1) a . (3) Evaluate the structural residual: R s = Ku (n) -f (n+1) s . (4) Solve S∆u (0,n) = -R s prior to initializing the fixed-point iteration. Iterate k = 0, ..., ∞ (5) Approximately solve ∂Ra ∂w ∆w (k+1,n) ≈ ∂Ra ∂x H∆u (k,n) by GMRES where ǫ g = 1e -1 . (6) Subsequently compute S∆u (k,n) = K∆u (k,n) -T ∂fa ∂w ∆w (k+1
Terminate for k = 0.

(8) Update the structural solution by the correction: u (n+1) = u (n) + ∆u (k+1,n) . ( 9) Update the aerodynamic mesh by: x (n+1) = x (0) + Hu (n+1) .

IV. Sensitivity Analysis

Sensitivity analysis allows for the evaluation of gradients that are used by the optimizer to determine the optimum. A response function, J , can be expressed as a function of one or several design variables. The main purpose of sensitivity analysis is to determine the total gradient of the response function with respect to the design variables, dJ /dp. In aeroelasticity, the response function can have an explicit dependence on the design variables, as well as an implicit dependence by the solution of the underlying governing equations and the mesh deformation problem. The total gradient is therefore obtained by linearizing ( 16)

dJ dp = ∂J ∂x dx dp + ∂J ∂w dw dp + ∂J ∂u du dp . (26)
The partial derivatives of the response function are typically straightforward to evaluate and can be provided directly by the respective solver or by external post-processing tools. The geometric gradient dx/dp is managed by the coupling module. To obtain the perturbed variables, dw/dp and du/dp, it is necessary to solve a linearization of the static aeroelastic problem at equilibrium. Hence, Equation ( 17) is linearized with respect to the design variables

The solution of this linear system is equivalent to the direct approach, indicating that the number of linear systems correspond to the number of design variables. This approach is therefore convenient when the optimization task is limited to a few design variables. To solve the linear system by GMRES a quality preconditioner is required. The task of the preconditioner is to improve the spectral properties of the system matrix, such that when applied to an iterative solver the number of iterations decrease. Let us assume that the system matrix in Equation ( 27) can be expressed as

A = A B 1 B 2 C (28)
where A, B 1 , B 2 and C are blocks that form a generalized saddle point system, A. A block LDUfactorization of this matrix would then yield

A B 1 B 2 C = I 0 B 2 A -1 I L A 0 0 S D I A -1 B 1 0 I U ( 29
)
where S is the Schur complement of A in the system matrix A. The diagonal matrix can then be our candidate preconditioner, P d = D. Indeed, it is highlighted in the review paper by Benzi et al., [START_REF] Benzi | Numerical solution of saddle point problems[END_REF] that the preconditioned system matrix M = P -1 d A has three distinct eigenvalues. This implies that a Krylov-based solver, e.g. GMRES, applied to the preconditioned linear system would converge in three iterations or less. Based on this premise, we adopt the block diagonal preconditioner to our linear sensitivity problem (27) in an attempt to replicate the favourable spectral properties. This yields in

P d = ∂Ra ∂w 0 0 S . (30)
The difficulty of applying this preconditioner is that the inverse of the large sparse multi-banded flux Jacobian, ∂R a /∂w, must be determined. To circumvent this impracticality the following approximated block diagonal preconditioner is constructed:

P d = ∂ Ra ∂w 0 0 S (31)
where ∂ R a /∂w is the aerodynamic preconditioner covered in Subsection A and S is the low-fidelity Schur complement covered in Subsection D.

Another important issue that emanates for the multidisciplinary linear system is associated with scaling. If the preconditioner is not properly scaled when applied to GMRES, the aerodynamic variables might become accurate while the structural variables are rendered inaccurate, or vice versa. [START_REF] Wathen | Preconditioning and convergence in the right norm[END_REF] We therefore redefine the linear problem ( 27) by a Schur complement formulation applied to the structural variables

S du dp = - ∂R s ∂p . (32)
Once the perturbed displacements du/dp are determined, the perturbed conservative variables dw/dp may be computed as a pure aerodynamic sensitivity problem by solving

∂R a ∂w dw dp = - ∂R a ∂u du dp . (33)
In fact, the perturbed conservative variables are implicitly computed during the formulation of the Schur matrix-vector product. Hence, this approach presents an opportunity to solve the coupled sensitivity equations absent the scaling problem. A similar approach to the PSQN solver is employed to solve (32). It is only modified to solve a linear problem, compared to the aeroelastic analysis that is nonlinear. This method is decribed in the following subsection.

Preconditioned Schur Complement (PSC) solver

We propose to use a preconditioned Schur complement solver by reformulating the linear problem (32) in fixed-point form

S du dp (k+1) - du dp (k)
= -S du dp

(k) + ∂R s ∂p (34)
which is identical to Equation ( 25) with the exception of the right-hand side. This is a favourable feature that can reduce the implementation costs. The routines in 4-7 in Table 3 for the static aeroelastic analysis can be reused for the sensitivity analysis by only modifying the right-hand side and performing multiple fixed-point iterations. To ensure that the total gradients for the optimizer are accurate, the convergence is monitored by the criterion:

dR (k+1) a 2 ≤ ǫ a2 dR (0) a 2 (35) dR (k+1) s 2 ≤ ǫ s2 dR (0) s 2 (36)
where dR s and dR a are the residual of Equations ( 32) and (33), respectively. Both criteria must be satisfied simultaneously before the solution is considered to be converged. This solution method, based on its similarities to the PSQN, is referred to as the preconditioned Schur complement (PSC) solver and is described in Table 4. The iteration index n is constant and is thus omitted for brevity. (2) Subsequently compute S du dp k+1) .

(k) = K du dp (k) -T ∂fa ∂w dw dp (
(3) Determine the fixed-point correction, du dp (k+1) , by solving Equation (34).

V. Results

Two test cases are conducted: one simulating subsonic flight at M ach = 0.60 and one simulating transonic flight at M ach = 0.85. The purpose is to assess the performance of our proposed method for a flight regime where the aeroloads are in agreement with the HiFi model (M ach = 0.60) and a flight regime where the aeroloads are presumed to have a higher offset (M ach = 0.85). The aerodynamic model and the solution parameters are identical to the test case in Subsection A. The impact of the number of implicit flow iterations, n a , on the aeroelastic solution will be investigated by performing two subcases: 1) n a = 100 and 2) n a = 300. A v-cycle multi-grid acceleration technique with two coarse grids and a fixed CFL number of 10 are used to speed up the convergence of the steady-state equations. The number of cycles for the LU-SSOR preconditioner is reduced to 4, whereas the number of cycles is maintained at 8 for the gradient computation in the Schur matrix-vector product. The flow problem is solved on five processors (one for each block) and the mesh deformation and the LoFi Schur complement formulation are split on 16 processors each.

The structural model is a classic wing box layout containing 1812 degrees of freedom. Member thicknesses and sections have been designed in a pre-processing optimization step (see Figure 4a). This model can be easily tuned in order to control flexibility and consequently aeroelastic effects. Small adjustments have been applied to the original model in Blondeau et al. [START_REF] Blondeau | Recent Achievements Towards Aero-Structure Gradient Computation using High-Fidelity CFD-CSM in the ONERA elsA Software[END_REF] RBE2 elements are added to the ribs with force nodes at z=0 (see Figure 4b). This allows the interpolated loads from the panel method to have a smooth and acceptable distribution. Three force nodes per rib are added for the load transfer. The first one is at 50% of the wing surface leading edge and the wing box leading edge. The second one is located in the centre of the wing box and the third one is at 50% of the wing surface trailing edge and the wing box trailing edge. The governing structural equations (7) in this work are solved by a direct method.

Comparison of aeroelastic analysis results

The static aeroelastic problem is solved individually by the three methods in Section III. It should be noted that a pure aeroelastic computation is performed with a clamped condition at the root. The equilibrium solution is therefore not a trimmed solution, which means that the aerodynamic lift does not equal the structural weight. The CPU time is monitored for the CFD analysis, the mesh deformation, the formulation of the LoFi Schur complement and the formulation of the approximated HiFi Schur matrix-vector product, respectively. The CPU time for the structural analysis is disregarded in this work as the structural equations are solved in less than 0.01 seconds. The convergence tolerances are set to ǫ a1 = ǫ s1 = 1e -6 .

The M ach = 0.60 test case yields a 7% tip displacement of the semispan. The best performing method for the n a = 100 subcase is the NLBGS (see Table 5). This can be explained by the moderate flexibility at this particular flight condition. No improvement is observed for the remaining methods and they all converge in 9 iterations. For the n a = 300 subcase, the SQN and the PSQN converge in fewer iterations compared to the NLBGS. The reason is that the aeroloads are evaluated at a higher precision, prompting the quadratic convergence of Newton's method to become more pronounced. This is illustrated in Figure 5b. However, the best performing method for the M ach = 0.60 test case in general is the NLBGS method. The improved robustness of the PSQN method is cancelled out by the increased CFD analysis time. The equilibrium tip displacement for the M ach = 0.85 test case is 13% of the semispan. The best performing method for n a = 100 is the SQN which shows a 36% CPU reduction to the equivalent subcase for NLBGS. Significant reduction in computation times is achieved by accelerating the convergence through the LoFi Schur complement, which is depicted in Figure 6a. The PSQN shows no advantage over the SQN in terms of number of iterations. When the flow subiteration count is increased to n a = 300, as for the M ach = 0.60 test case, an improved convergence rate is manifested for the PSQN. This can clearly be observed in Figure 6b. However, the improved convergence rate comes at the expense of an increased CFD analysis time. Hence, the overall top performer for the M ach = 0.85 case if the SQN method. The flow subiteration count could arguably be decreased further (n a < 100) to give the NLBGS the edge over the SQN. However, this is not advisable as it will result in many more coupling iterations, which are illustrated by the peaks in the figures. Moreover, the Mesh-to-CFD CPU ratio might increase up to where one mesh deformation is comparable to one CFD analysis. It is not desirable for a large portion of the CPU time of the aeroelastic problem to be devoted to mesh deformation.

Comparison of sensitivity analysis results

The design variable that is used to construct ∂R s /∂p for the test case is the thickness of the shell element on the upper skin closest to the root and the trailing edge. This term is constructed by multiplying the perturbed stiffness matrix with the equilibrium displacements (9) obtained by the previous analysis. The application for the analysis and sensitivity analysis in this work is structural sizing or aeroelastic tailoring through directional stiffness variation. This implies that the design variables have no explicit influence on the aerodynamic residuals, ∂R a /∂p = 0. For the sensitivity analysis we set ǫ a2 = ǫ s2 = 1e -8 . Two methods are compared in the sensitivity analysis. The first one is the PSC method as described in Section IV. The second one is identical with one minor difference: the LoFi Schur complement preconditioner is replaced by the structural stiffness matrix, such that Equation (34) is reduced to

K du dp (k+1) - du dp (k) = -S du dp (k) + ∂R s ∂p . (37)
The assumption for this modification is that no panel method is available to precondition the fixed-point linear system. This modification results in a looser coupling and is distinguished with a subscript as PSC K . The original version (34) is therefore identified for this test case with an added subscript as PSC S . Note that the SQN method from the aeroelastic analysis is not applicable in the sensitivity analysis. The direct application of the LoFi Schur complement would only result in LoFi gradients.

The LoFi Schur complement preconditioner indicates to have an impact on the convergence rate for both test cases, as shown in Table 6. Since we are solving a linear system, only one generation of the LoFi Schur complement is required. The CPU cost of this operation is negligible compared to the overall CPU cost. Moreover, the subcases of setting a fix number of flow iterations, n a , as in the aeroelastic analysis do not come into play in the sensitivity analysis. Here we are only concerned with gradient computation.

There is only a moderate improvement at M ach = 0.60 by using PSC S . A 17% CPU reduction is observed and the convergence behaviour in Figure 7a does not display much difference in the decaying residual between the two methods. For the more flexible M ach = 0.85 case, the application of PSC S results in more than 3 times faster convergence compared to PSC K . The addition of the LoFi Schur complement improves the stability and enables it to converge much faster, as illustrated in Figure (7b). This clearly indicates that S is more prominent choice of preconditioner. It should be noted that one fixed-point iteration for the M ach = 0.85 case is almost twice as costly as the M ach = 0.60 case. This might suggest that the M ach = 0.85 case experiences a higher degree of ill-conditioning.

The elegance of the PSC S is that the coupled equations ( 27) are practically transformed to an aerodynamic sensitivity problem. The mesh deformation operations, that would otherwise be necessary for each iteration, are reduced to a minimum. Only 9 and 11 mesh deformations are required for the M ach = 0.60 and M ach = 0.85 case, respectively, to converge the perturbed variables to a tolerance of 1e -8 . The bulk of the computation time is devoted to solving Equation (33) which is embedded in the Schur matrix-vector product, S du dp .

VI. Conclusions

We have described a panel code preconditioning method that facilitates the solution of high-fidelity aeroelastic analysis and sensitivity analysis problems. The panel code preconditioner is expressed in Schur complement form applied to the structural degrees of freedom. It is used to improve the spectral properties of the linear subproblems that appear in the aeroelastic analysis, when subject to Newton's method, and the subsequent sensitivity analysis. The linear subproblems are formulated in fixed-point iterative form and are solved approximately for the static aeroelastic problem and exactly for the sensitivity problem. The method applied for the static aeroelastic problem, which we refer to as the PSQN method, is not necessarily the optimal choice. It converges in less iterations compared to the NLBGS and the SQN method, but it also requires a higher precision on the aeroloads for each Newton iteration. It is the most robust method, in this regard, but not the most efficient in terms of CPU time.

The sensitivity analysis version of the PSQN, the PSC, demonstrates much faster convergence when preconditioned with the LoFi Schur complement. The convergence rate is up to three times faster compared to the case when the PSC is preconditioned with the structural stiffness matrix. Moreover, the reduction of the problem in Schur complement form effectively reduces the number of mesh deformations to a minimum. The number of mesh deformations correspond to the number of fixed-point iterations which varied from 9 to 11 for the two test cases conducted in this paper. These values can be further reduced if the solution is allowed to converge to a higher tolerance. In conclusion, the panel code preconditioner is successfully demonstrated to reduce the computation time of high-fidelity sensitivity analysis for fluid-structure interaction problems.

The PSC method described in this paper for the direct approach gradient computation is directly applicable to the adjoint approach. The full implementation of the method for the adjoint problem is envisioned in a follow-up paper. Furthermore, trim analysis will be performed to allow for response generation that take into account trimmed flight for gradient-based optimization.

  (a) Camber surface panels for the panel method. (b) Structured Euler mesh for elsA.

Figure 1 :

 1 Figure 1: Aerodynamic models of the ONERA M6 wing.

Figure 2 :

 2 Figure 2: Convergence plot of the direct method gradient computation.

  Parallel computation scalability. Simulations are performed on 2, 4, 8, 16 and 32 processors. (b) Surface point reduction on the ONERA M6 wing. Only black points are used for the coupling matrix.

Figure 3 :

 3 Figure 3: Coupling module improvements for accelerated mesh deformation.

( 0 )( 1 )

 01 Determine ∂Rs ∂p in a post-analysis step. Solve S du dp (0) = -∂Rs ∂p prior to initializing the fixed-point iteration. Iterate k = 0, ..., ∞ Approximately solve ∂Ra ∂w dw dp (k+1) ≈ ∂Ra ∂x H du dp (k) by GMRES where ǫ g = 1e -1 .

  (a) Original model. Shell thickness variation is displayed on the upper surface. (b) Modified model. Rigid RBE2 elements are added. Aeroloads are interpolated on the force (red) nodes.

Figure 4 :

 4 Figure 4: Finite element model.

  Residual convergence for na = 300.

Figure 5 :

 5 Figure 5: Convergence histories for the M ach = 0.60 test case.

  Residual convergence for na = 300.

Figure 6 :

 6 Figure 6: Convergence histories for the M ach = 0.85 test case.

  Residual convergence for M ach = 0.85.

Figure 7 :

 7 Figure 7: Convergence histories for the sensitivity analysis.

Table 4 :

 4 The PSC algorithm

Table 5 :

 5 Aeroelastic Analysis

	CPU time [s]

Table 6 :

 6 Sensitivity Analysis

	CPU time [s]

VII. Acknowledgements

The research leading to these results has received funding from the European Union Seventh Framework Programme FP7-PEOPLE-2012-ITN under grant agreement n 316394, Aerospace Multidisciplinarity-Enabling DEsign Optimization (AMEDEO) Marie Curie Initial Training Network. The authors would like to acknowledge Philippe Girodroux-Lavigne for his support on all actions that allowed us to use the flow solver elsA.


MD, PhD João Pedro Ferreira
email: j.ferreira@chru-nancy.fr

Msc Ana Cristina Oliveira

Msc Francisca A Saraiva

MD, Msc Francisco Vasques-Nóvoa

MD, PhD Adelino Leite-Moreira

Sodium-glucose co-transporter inhibitors in insulin-treated diabetes: a meta-analysis

Keywords: type 2 diabetes, insulin, sodium-glucose cotransporter 2 inhibitors, cardiovascular outcomes, meta-analysis

Background: Patients with insulin-treated type 2 diabetes (T2D) have a high risk of major adverse cardiovascular events. Sodium-glucose cotransporter inhibitors (SGLTi) improve outcomes without hypoglycemic risk.

Aims:

To study the effect of SGLTi in patients with T2D with and without background insulintreatment in outcome-driven RCTs.

Methods: Random effects models.

Results: A total of 54,374 patients with T2D were included in the analysis, of which 26,551 (48.8%) were treated with insulin. For 3P-MACE in patients without insulin treatment, the HR (95%CI) for the effect of SGLTi vs. placebo was 0.93 (0.81-1.05), with moderate heterogeneity (I2 =49.2%, Q statistic P =0.11). In insulin-treated patients, the HR (95%CI) was 0.88 (0.82-0.95), without evidence of heterogeneity (I2 =0.0%, Q statistic P =0.91). The pooled effect evidenced a 10% reduction of 3P-MACE with SGLTi (HR 0.90, 95%CI 0.85-0.96), without SGLTi-by-insulin interactionP=0.53. For the composite outcome of HF hospitalisation or cardiovascular death in patients without insulin treatment, the HR (95%CI) for the effect of SGLTi vs. placebo was 0.77 (0.61-0.92), with marked heterogeneity (I2 =66.8%, Q statistic P =0.02). In insulin-treated patients, the HR (95%CI) was 0.77 (0.68-0.86), without significant heterogeneity (I2 =31.7%, Q statistic P =0.25). The pooled effect evidenced a 23% reduction of HF hospitalisations or cardiovascular death with SGLTi (HR 0.77, 95%CI 0.68-0.85), without SGLTi -by-insulin interactionP=0.98.

Conclusion:

SGLTi reduce cardiovascular events regardless of insulin use. However, the treatment effect is more homogeneous among insulin-treated patients, supporting the use of SGLTi for the treatment of patients with T2D requiring insulin for glycemic control.

Introduction

Patients with insulin-treated diabetes have a high risk of major adverse cardiovascular events (MACE), frequently leading to symptomatic heart failure (HF), hospitalizations and death [START_REF] Herman | Insulin Therapy Increases Cardiovascular Risk in Type 2 Diabetes[END_REF] . The risk of patients with type 2 diabetes (T2D) treated with insulin is higher than the risk of patients with T2D treated with oral anti-diabetic drugs only (i.e., non-insulin-treated diabetes) [START_REF] Macdonald | Impact of diabetes on outcomes in patients with low and preserved ejection fraction heart failure: an analysis of the Candesartan in Heart failure: Assessment of Reduction in Mortality and morbidity (CHARM) programme[END_REF][START_REF] Laing | Mortality from heart disease in a cohort of 23,000 patients with insulin-treated diabetes[END_REF][START_REF] Rossello | Impact of insulin-treated diabetes on cardiovascular outcomes following high-risk myocardial infarction[END_REF] . In patients with T2D, insulin treatment is associated with long-standing diabetes, poor cardiometabolic profile (higher body mass index, abdominal obesity, heart rate, blood pressure, triglycerides and lower HDL cholesterol) and a higher prevalence of micro-and macrovascular complications [START_REF] Braffett | Association of Insulin Dose, Cardiometabolic Risk Factors, and Cardiovascular Disease in Type 1 Diabetes During 30 Years of Follow-up in the DCCT/EDIC Study[END_REF][START_REF] Holden | Glucoselowering with exogenous insulin monotherapy in type 2 diabetes: dose association with all-cause mortality, cardiovascular events and cancer[END_REF] .

Sodium-glucose cotransporter 1/2 or 2 inhibitors (SGLTi) have demonstrated to improve cardiovascular and kidney outcomes in patients with T2D and in patients with heart failure with reduced ejection fraction (HFrEF) with and without diabetes [START_REF] Zinman | Cardiovascular Outcomes, and Mortality in Type 2 Diabetes[END_REF][START_REF] Neal | Canagliflozin and Cardiovascular and Renal Events in Type 2 Diabetes[END_REF][START_REF] Wiviott | Dapagliflozin and Cardiovascular Outcomes in Type 2 Diabetes[END_REF][START_REF] Perkovic | Canagliflozin and Renal Outcomes in Type 2 Diabetes and Nephropathy[END_REF][START_REF] Zannad | SGLT2 inhibitors in patients with heart failure with reduced ejection fraction: a meta-analysis of the EMPEROR-Reduced and DAPA-HF trials[END_REF] . The reduction of hospitalisations due to HF is particularly striking with these agents [START_REF] Zelniker | SGLT2 inhibitors for primary and secondary prevention of cardiovascular and renal outcomes in type 2 diabetes: a systematic review and meta-analysis of cardiovascular outcome trials[END_REF][START_REF] Neuen | SGLT2 inhibitors for the prevention of kidney failure in patients with type 2 diabetes: a systematic review and meta-analysis[END_REF] . As patients with insulin-treated diabetes have a high-risk of cardiovascular complications, including HF, it is plausible that SGLTi provide major benefits to these patients by potentially reducing their risk of MACE, including HF hospitalisations and cardiovascular death. Furthermore, SGLTi may help controlling blood pressure, weight, and sodium retention, thus counteracting the deleterious effects associated with insulin-treatment without increasing the risk of hypoglycemic episodes [START_REF] Ferreira | Empagliflozin for Patients with Presumed Resistant Hypertension: a Post Hoc Analysis of the EMPA-REG OUTCOME Trial[END_REF][START_REF] Neal | Efficacy and safety of canagliflozin, an inhibitor of sodium-glucose cotransporter 2, when used in conjunction with insulin therapy in patients with type 2 diabetes[END_REF][START_REF] Rosenstock | Impact of empagliflozin added on to basal insulin in type 2 diabetes inadequately controlled on basal insulin: a 78-week randomized, double-blind, placebo-controlled trial[END_REF][START_REF] Rosenstock | Improved glucose control with weight loss, lower insulin doses, and no increased hypoglycemia with empagliflozin added to titrated multiple daily injections of insulin in obese inadequately controlled type 2 diabetes[END_REF][START_REF] Wilding | Dapagliflozin in patients with type 2 diabetes receiving high doses of insulin: efficacy and safety over 2 years[END_REF][START_REF] Wilding | Longterm efficacy of dapagliflozin in patients with type 2 diabetes mellitus receiving high doses of insulin: a randomized trial[END_REF][START_REF] Wilding | A study of dapagliflozin in patients with type 2 diabetes receiving high doses of insulin plus insulin sensitizers: applicability of a novel insulin-independent treatment[END_REF] .

Using updated data from large outcome-driven randomised controlled trials (RCTs) assessing the efficacy and safety of SGLTi versus placebo, we aim to study the effect of these agents in patients with T2D with and without background insulin-treatment.

Methods

Search strategy and selection criteria

This meta-analysis was conducted in accordance with the Preferred Reporting Items for Systematic Reviews and Meta-analyses (PRISMA) reporting guideline [START_REF] Liberati | The PRISMA statement for reporting systematic reviews and meta-analyses of studies that evaluate healthcare interventions: explanation and elaboration[END_REF] . We defined both MeSH terms (controlled language) and free text terms to express each component of PICO randomized controlled trial) AND (cardiovascular death OR cardiovascular events). We did not restrict the search by language, date or publication status. Additionally, we have included recently published trials that were not displayed on PubMed at the time [START_REF] Bhatt | Sotagliflozin in Patients with Diabetes and Recent Worsening Heart Failure[END_REF][START_REF] Bhatt | Sotagliflozin in Patients with Diabetes and Chronic Kidney Disease[END_REF] [START_REF] Bhatt | Sotagliflozin in Patients with Diabetes and Chronic Kidney Disease[END_REF] . Trial eligibility was confirmed by 2 independent reviewers (J.P.F., A.C.O.), data extraction was performed by A.C.O and J.P.F, and the study quality was performed by J.P.F., A.C.O., and F.A.S. We have also cross-checked the RCTs selected to the present analysis with other recently published meta-analyses and the concordance was perfect [START_REF] Zelniker | SGLT2 inhibitors for primary and secondary prevention of cardiovascular and renal outcomes in type 2 diabetes: a systematic review and meta-analysis of cardiovascular outcome trials[END_REF][START_REF] Neuen | SGLT2 inhibitors for the prevention of kidney failure in patients with type 2 diabetes: a systematic review and meta-analysis[END_REF][START_REF] Mcguire | Association of SGLT2 Inhibitors With Cardiovascular and Kidney Outcomes in Patients With Type 2 Diabetes: A Meta-analysis[END_REF] . The Evaluation of the Effects of Canagliflozin on Renal and Cardiovascular Outcomes in Participants With Diabetic Nephropathy (CREDENCE) [START_REF] Perkovic | Canagliflozin and Renal Outcomes in Type 2 Diabetes and Nephropathy[END_REF] , which was primarily a kidney outcomes trial, did not report treatment effects in insulin subgroups and was not included in the analysis. A detailed description of each individual trial is provided in the respective main reports as above referenced.

The analyses were performed to assess heterogeneity of treatment effects in outcomes in patients with and without background insulin treatment.

Outcomes

The primary outcome in EMPA-REG OUTCOME, CANVAS-Program, DECLARE-TIMI 58, and VERTIS-CV was the composite of time-to-first of myocardial infarction, stroke, or cardiovascular death (3P-MACE). The primary outcome in SOLOIST-WHF and SCORED was the total number of deaths from cardiovascular causes and hospitalizations and urgent visits for HF (first and recurrent). On top of the primary outcome, the EMPA-REG OUTCOME trial reported the treatment effect of empagliflozin in patients with and without insulin treatment for the outcome of cardiovascular death. The DECLARE-TIMI 58 and VERTIS-CV reported the treatment effect of dapagliflozin and ertugliflozin, respectively, in patients with and without insulin treatment for the composite outcome of time-to-first of cardiovascular death or HF hospitalisation. For optimizing the consistency of outcomes across trials, the EMPA-REG OUTCOME, CANVAS-Program, DECLARE-TIMI 58, and VERTIS-CV trials were meta-analyzed for the composite outcome of time-to-first 3P-MACE of myocardial infarction, stroke, or cardiovascular death. The EMPA-REG OUTCOME, DECLARE-TIMI 58, VERTIS-CV, SOLOIST-WHF, and SCORED trials were meta-analyzed for the outcomes of cardiovascular death (EMPA-REG OUTCOME) and composites of cardiovascular death or HF hospitalisations, as time-to-first event (DECLARE-TIMI 58 and VERTIS-CV) or first and recurrent (SOLOIST-WHF and SCORED). The studied trials and outcomes are displayed in Table 1.

Statistical analysis

A random-effect meta-analysis approach was used with heterogeneity assessed using the Cochran Q test statistic and Higgins and Thompson I2 [START_REF] Higgins | Quantifying heterogeneity in a meta-analysis[END_REF] . Heterogeneity was considered to be low, moderate, or high if I2 was <25%, 25% to 75%, or >75%, respectively. The hazard ratio (HR) and respective 95% confidence intervals (95%CI) were extracted from the insulin subgroup analyses reported in the publications of each individual study. All studies, except DECLARE-TIMI 58, reported the HR (95%CI) for patients with and without insulin treatment.

In DECLARE-TIMI 58 we have extracted the approximate HR (95%CI) from the forest plots.

The number and proportion of patients with events in insulin subgroups were inconsistently reported and we show these numbers in a table whenever available. The meta-analysis was directly implemented on the HR (exponential) scale. The assessment of whether background insulin treatment could modify the effect of SGLTi on outcomes was performed using interaction term for subgroups and reported with a P-value for interaction. Records

Results

The characteristics of the patients from the RCTs included in this analysis are depicted in the Table 2. A total of 54374 patients with T2D were included in the analysis, of which 26551 (48.8%) were treated with insulin. The mean age of all trial participants was 66 years, 19702 (36.2%) were women, and 43806 (80.6%) were Caucasian. The number of participants per trial ranged from 1222 in SOLIST-WHF to 17160 in DECLARE-TIMI 58. The median follow-up ranged from 0.8 years in SOLIST-WHF to 4.2 years in DECLARE-TIMI 58 (Table 2). The assessment of risk of bias by visual inspection of funnel plots and quantitative assessment suggested no indication of publication bias. The risk of bias in randomized trials as assessed by the RoB 2 tool was low.

The number of events per trial and insulin subgroup are presented in the Table 3.

The treatment effect of SGLTi vs. placebo in patients with and without insulin treatment for the composite outcome of myocardial infarction, stroke or cardiovascular death and the composite of heart failure hospitalisation or cardiovascular death (cardiovascular death alone in the EMPA-REG OUTCOME trial) are presented in the Table 3 and Figure 1.

For the composite outcome of myocardial infarction, stroke or cardiovascular death in patients without insulin treatment, the HR (95%CI) for the effect of SGLTi vs. placebo was 0.93 (0.81-1.05), with moderate heterogeneity (I2 =49.2%, Q statistic P =0.11). In insulintreated patients, the HR (95%CI) for the effect of SGLTi vs. placebo was 0.88 (0.82-0.95), without evidence of heterogeneity (I2 =0.0%, Q statistic P =0.91). The pooled effect evidenced a 10% reduction of 3P-MACE with SGLTi (HR 0.90, 95%CI 0.85-0.96), without significant interaction between subgroups of insulin (yes vs. no), interaction P =0.53. Figure 1,panel A.

For the composite outcome of HF hospitalisation or cardiovascular death (cardiovascular death alone in the EMPA-REG OUTCOME trial) in patients without insulin treatment, SGLTi reduced the event rate by 23% (HR 0.77, 95%CI 0.61-0.92) but the treatment effect was heterogeneous (I2 =66.8%, Q statistic P =0.02). In insulin-treated patients, SGLTi significantly reduced HF hospitalisation or cardiovascular death (HR 0.77, 95%CI 0.68-0.86) without significant heterogeneity (I2 =31.7%, Q statistic P =0.25). The pooled effect evidenced a 23% reduction of HF hospitalisations or cardiovascular death with SGLTi (HR 0.77, 95%CI 0.68-0.85), without significant interaction between subgroups of insulin (yes vs. no), interaction P =0.98. Figure 1,panel B.

Discussion

The present meta-analysis shows that SGLTi reduced the composite outcome of myocardial infarction, stroke or cardiovascular death and the composite outcome of HF hospitalisation or cardiovascular death regardless of insulin use at baseline. However, the treatment effect was more homogeneous and consistent among insulin-treated patients, highlighting the need for adding a SGLTi to the therapeutic armamentarium of patients with T2D requiring insulin for glycemic control.

In patients with T2D, insulin is used to lower blood glucose levels in patients that fail to achieve targets or that do not tolerate oral glucose-lowering therapies. However, the intensity of glucose lowering should be tailored to each individual patient as it may increase hypoglycemic episodes and mortality [START_REF] Turnbull | Intensive glucose control and macrovascular outcomes in type 2 diabetes[END_REF][START_REF] Gerstein | Effects of intensive glucose lowering in type 2 diabetes[END_REF] . Furthermore, insulin-treatment may lead to weight gain, fluid retention, higher blood pressure, and increase the risk of hypoglycemia without reducing cardiovascular complications [START_REF]Intensive blood-glucose control with sulphonylureas or insulin compared with conventional treatment and risk of complications in patients with type 2 diabetes (UKPDS 33). UK Prospective Diabetes Study (UKPDS) Group[END_REF][START_REF] Knatterud | Effects of hypoglycemic agents on vascular complications in patients with adult-onset diabetes. VII. Mortality and selected nonfatal events with insulin treatment[END_REF][START_REF] Lebovitz | Insulin: potential negative consequences of early routine use in patients with type 2 diabetes[END_REF] . Patients with T2D who require insulin, usually have a poor glycemic and metabolic control, and have a high risk of cardiovascular complications, including HF hospitalisations and death from cardiovascular causes [START_REF] Rossello | Impact of insulin-treated diabetes on cardiovascular outcomes following high-risk myocardial infarction[END_REF][START_REF] Shen | Insulin treatment and clinical outcomes in patients with diabetes and heart failure with preserved ejection fraction[END_REF] . For these reasons, there is an urgent need to reduce the risk profile of insulin-treated patients, while also improving their metabolic components, reduce weight, and decrease the burden of hospitalisations and adverse cardiovascular outcomes. Inhibitors of the sodium glucose co-transporter have many beneficial effects beyond glucose control, and, indeed, may help to reduce weight, control blood pressure, limit sodium retention, reduce hospitalisations (HF-related in particular) and cardiovascular death, counteracting the negative effects of insulin. As shown in the present study, SGLTi reduce cardiovascular complications regardless of insulin treatment. However, the treatment effect was more homogeneous among insulintreated patients leading to consistent reductions of cardiovascular complications, comprising HF hospitalisations and cardiovascular death across multiple populations with T2D, including patients with chronic kidney disease and worsening heart failure [START_REF] Zinman | Cardiovascular Outcomes, and Mortality in Type 2 Diabetes[END_REF][START_REF] Neal | Canagliflozin and Cardiovascular and Renal Events in Type 2 Diabetes[END_REF][START_REF] Wiviott | Dapagliflozin and Cardiovascular Outcomes in Type 2 Diabetes[END_REF][START_REF] Bhatt | Sotagliflozin in Patients with Diabetes and Recent Worsening Heart Failure[END_REF][START_REF] Bhatt | Sotagliflozin in Patients with Diabetes and Chronic Kidney Disease[END_REF][START_REF] Cannon | Cardiovascular Outcomes with Ertugliflozin in Type 2 Diabetes[END_REF] . This effect was seen regardless of the type of SGLTi i.e., the SGLT1/2 inhibitor sotagliflozin and the remainder SGLT2 inhibitors had similar effects to reduce cardiovascular events regardless of insulin treatment, but with a more homogeneous effect among insulin-treated patients.

A recently published meta-analysis showed that in placebo-controlled RCTs SGLT2i reduced cardiovascular events and HF hospitalisations, independently of the presence of atherosclerotic cardiovascular disease, kidney function, and glycated hemoglobin levels [START_REF] Mcguire | Association of SGLT2 Inhibitors With Cardiovascular and Kidney Outcomes in Patients With Type 2 Diabetes: A Meta-analysis[END_REF] . It is likely that insulin-treated patients also had higher burden of cardiovascular disease, worse kidney and higher glycated hemoglobin levels; however, beyond the higher cardiovascular risk (with potentially larger treatment benefit), insulin-treated patients may experience a more consistent benefit from SGLTi due to the effects and mechanisms associated with insulin administration, such as weight gain and sodium retention that may be counterbalanced by SGLTi. Supporting this hypothesis, smaller "mechanistic" studies testing the effects of SGLT2i among insulin-treated patients also showed a beneficial effect of these agents. For example, a subsample of patients enrolled in the CANVAS-Program and receiving insulin at baseline were randomized to receive placebo (n =690), canagliflozin 100 mg (n =692), or canagliflozin 300 mg (n =690). Compared with the placebo both canagliflozin doses reduced glycated hemoglobin, fasting plasma glucose, body weight and blood pressure [START_REF] Neal | Efficacy and safety of canagliflozin, an inhibitor of sodium-glucose cotransporter 2, when used in conjunction with insulin therapy in patients with type 2 diabetes[END_REF] . Similar findings were observed with dapagliflozin and empagliflozin in patients with inadequately controlled T2D on basal insulin [START_REF] Rosenstock | Impact of empagliflozin added on to basal insulin in type 2 diabetes inadequately controlled on basal insulin: a 78-week randomized, double-blind, placebo-controlled trial[END_REF][START_REF] Wilding | Dapagliflozin in patients with type 2 diabetes receiving high doses of insulin: efficacy and safety over 2 years[END_REF] and obese difficult-to-treat patients with T2D requiring high doses of insulin [START_REF] Rosenstock | Improved glucose control with weight loss, lower insulin doses, and no increased hypoglycemia with empagliflozin added to titrated multiple daily injections of insulin in obese inadequately controlled type 2 diabetes[END_REF] . It should be noted, that SGLT2i also reduced insulin requirements without increased risk of hypoglycemia [START_REF] Rosenstock | Improved glucose control with weight loss, lower insulin doses, and no increased hypoglycemia with empagliflozin added to titrated multiple daily injections of insulin in obese inadequately controlled type 2 diabetes[END_REF][START_REF] Wilding | Dapagliflozin in patients with type 2 diabetes receiving high doses of insulin: efficacy and safety over 2 years[END_REF][START_REF] Wilding | Longterm efficacy of dapagliflozin in patients with type 2 diabetes mellitus receiving high doses of insulin: a randomized trial[END_REF] .

Limitations

Several limitations are acknowledged. Differences in eligibility criteria, patient characteristics, trial duration, and outcome definition may have affected the treatment effects and increased heterogeneity in the meta-analysis. However, evidence of some heterogeneity was only found among patients without background insulin treatment. The trials here reported, randomised patients to either a SGLTi or placebo but insulin treatment was not part of the randomisation procedure; therefore, these subgroup analyses should be regarded with caution because the treatment effect may be influenced by insulin itself, but also by the patients` characteristics, their risk profile, comorbidities (including heart failure and chronic kidney disease) or the quality of healthcare. Notwithstanding, the homogeneity of effect across trials among insulin-treated patients suggests that SGLTi are a beneficial therapy for these patients. For purposes of the present analysis, we were limited to outcomes that had been reported by subgroups of insulin use and some outcome definitions varied between trials (e.g., time-to-first and total heart failure hospitalisations), however the treatment effects were similar across outcomes and this is not a major source of bias in our report. The number needed-to-treat to benefit is not reported because the data required for a reliable calculation of this estimate was not available in all trials.

Conclusions

The present meta-analysis shows that SGLTi reduce cardiovascular events regardless of insulin use. However, the treatment effect is more homogeneous among insulin-treated patients, supporting the use of SGLTi for the treatment of patients with T2D requiring insulin for glycemic control.

  expression: P) Population, patients with diabetes; I) Intervention, SGLT inhibitors; C) Comparison, placebo and O) Outcomes, MACEs. The search query included the following and it was last searched on 22 nd of November 2020: ("Diabetes Mellitus"[Mesh] OR diabetes) AND ("Sodium-Glucose Transporter 2 Inhibitors"[Mesh] OR sglt inhibitors) AND ("Placebos"[Mesh] OR placebo) AND ("Randomized Controlled Trial" [Publication Type] OR

  identified in each database were imported and managed through Microsoft Excel®. We assessed publication bias by visual inspection of funnel plots, with the ascertainment for potential asymmetry of published results by Egger's regression test and Duval and Tweedie's trim-and-fill method. The statistical analyses were performed using Stata® software version 16 (StataCorp. 2019. Stata Statistical Software: Release 16. College Station, TX: StataCorp LLC.) and the package meta. The risk of bias in randomised trials was assessed with the RoB 2 tool 27 . Two-sided P values <0.05 were considered significant.

F

  .A. Saraiva and F. Vasques-Nóvoa are supported by Universidade do Porto/FMUP and FSE-Fundo Social Europeu, NORTE 2020-Programa Operacional Regional do Norte, NORTE-08-5369-FSE-000024-Programas Doutorais.

Figure 1 .

 1 Figure 1. Treatment effects by insulin treatment subgroups

Acknowledgments

A.C.Oliveira is supported by FCT-Fundação para a Ciência e a Tecnologia, FSE-Fundo Social Europeu, Portugal 2020, NORTE 2020-Programa Operacional Regional do Norte, Centro 2020-Programa Operacional Regional do Centro, Alentejo 2020-Programa Operacional Regional do Alentejo (UI/BD/150647/2020).

Disclosures

J.P.F. is consultant for Boehringer Ingelheim. All the other authors have no relevant conflicts of interest to report.